В различных контекстах, например, для сокращения аргументов для математических функций, нужно вычислить (a - K) / (a + K), где a - аргумент положительной переменной, а K - постоянная. Во многих случаях K является степенью двух, что является прецедентом, относящимся к моей работе. Я ищу эффективные способы более точно вычислить этот коэффициент, чем это можно сделать при прямом разделении. Можно предположить аппаратную поддержку плавного многократного добавления (FMA), так как эта операция обеспечивается всеми основными архитектурами CPU и GPU в настоящее время и доступна в C/С++ с помощью функций fma() и fmaf().
Для удобства исследования я экспериментирую с арифметикой float. Поскольку я планирую перенести подход на арифметику double, то никакие операции, использующие более высокую, чем внутренняя точность как аргумента, так и результата, могут быть использованы. Мое лучшее решение до сих пор:
/* Compute q = (a - K) / (a + K) with improved accuracy. Variant 1 */
m = a - K;
p = a + K;
r = 1.0f / p;
q = m * r;
t = fmaf (q, -2.0f*K, m);
e = fmaf (q, -m, t);
q = fmaf (r, e, q);
Для аргументов a в интервале [K/2, 4.23*K] приведенный выше код вычисляет фактор, почти правильно округленный для всех входов (максимальная ошибка чрезвычайно близка к 0,5 ulps) при условии, что K является степенью 2, а там не является переполнением или недостаточным промежуточным результатом. Для K, не являющегося степенью двух, этот код еще точнее, чем наивный алгоритм, основанный на делении. Что касается производительности, этот код может быть быстрее, чем наивный подход на платформах, где обратная точка с плавающей запятой может быть вычислена быстрее, чем деление с плавающей запятой.
Я делаю следующее наблюдение, когда K= 2 n: когда верхняя граница рабочего интервала возрастает до 8*K, 16*K,... максимальная ошибка увеличивается постепенно и начинается медленно аппроксимировать максимальную ошибку наивного вычисления снизу. К сожалению, то же самое, как представляется, не верно для нижней границы интервала. Если нижняя граница падает до 0.25*K, максимальная ошибка улучшенного метода выше равна максимальной ошибке наивного метода.
Существует ли метод вычисления q = (a - K)/(a + K), который может достичь меньшей максимальной ошибки (измеренной в ulp по сравнению с математическим результатом) по сравнению с наивным методом и приведенную выше кодовую последовательность на более широком интервале, в частности, для интервалов, нижняя граница которых меньше 0.5*K? Эффективность важна, но может быть допущено несколько операций, чем в приведенном выше коде.
В одном из ответов ниже было указано, что я мог бы повысить точность, возвращая фактор как неоценимую сумму двух операндов, то есть в виде пары хвоста q:qlo, то есть аналогичной хорошо известной двойной - float и double- double. В моем коде выше это означало бы замену последней строки на qlo = r * e.
Этот подход, безусловно, полезен, и я уже рассматривал его использование для логарифма с расширенной точностью для использования в pow(). Но это принципиально не помогает с желательным расширением интервала, на котором расширенное вычисление обеспечивает более точные коэффициенты. В конкретном случае, на который я смотрю, я хотел бы использовать K=2 (для одиночной точности) или K=4 (для двойной точности), чтобы сохранить интервал первичной аппроксимации узким, а интервал для a примерно равен [0, 28]. Практическая проблема, с которой я столкнулся, заключается в том, что для аргументов < 0,25 * К точность улучшенного деления не намного лучше, чем с наивным методом.
