Эффективно вычислять (a - K) / (a + K) с улучшенной точностью
В различных контекстах, например, для сокращения аргумента для математических функций, необходимо вычислить(a - K) / (a + K)
, гдеa
является аргументом положительной переменной иK
постоянная Во многих случаях,K
является степенью двойки, которая является вариантом использования, относящимся к моей работе. Я ищу эффективные способы вычисления этого коэффициента более точно, чем это можно сделать с помощью простого деления. Можно предположить аппаратную поддержку слияния с многократным добавлением (FMA), поскольку в настоящее время эта операция предоставляется всеми основными архитектурами ЦП и ГП и доступна в C / 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 ульпс), при условии, чтоK
является степенью 2, и в промежуточных результатах нет переполнения или недостатка. ЗаK
этот код не является степенью двойки, но он все же более точен, чем простой алгоритм, основанный на делении. С точки зрения производительности этот код может бытьБыстрее чем наивный подход на платформах, где обратные числа с плавающей точкой могут быть вычислены быстрее, чем деление с плавающей точкой.
Я делаю следующее наблюдение, когдаK
= 2n: Когда верхняя граница рабочего интервала увеличивается до8*K
, 16*K
, ... максимальная ошибка постепенно увеличивается и начинает медленно приближаться к максимальной ошибке наивного вычисления снизу. К сожалению, то же самое не кажется верным для нижней границы интервала. Если нижняя граница падает до0.25*K
максимальная погрешность улучшенного метода выше равна максимальной погрешности наивного метода.
Существует ли метод для вычисления q = (a - K) / (a + K), который может обеспечить меньшую максимальную ошибку (измеряется вULP по сравнению с математическим результатом) по сравнению как с наивным методом, так и с вышеуказанной кодовой последовательностью, в более широком интервале,в частности, для интервалов, нижняя граница которых меньше0.5*K
? Эффективность важна, но скорее всего допустимо несколько больше операций, чем используется в приведенном выше коде.
В одном ответе ниже было указано, что я мог бы повысить точность, возвращая частное в виде неоцененной суммы двух операндов, то есть в виде пары голова-хвостq:qlo
то есть похож на хорошо известныйfloat
и дваждыdouble
форматы. В моем коде выше это означало бы изменение последней строки наqlo = r * e
.
Этот подход, безусловно, полезен, и я уже обдумывал его использование для логарифма с повышенной точностью для использования вpow()
, Но это принципиально не помогает с желаемым расширением интервала, в котором расширенные вычисления обеспечивают более точные коэффициенты. В конкретном случае я смотрю, я хотел бы использоватьK=2
(для одинарной точности) илиK=4
(для двойной точности), чтобы сохранить интервал первичного приближения узким, а интервал дляa
примерно [0,28]. Практическая проблема, с которой я сталкиваюсь, состоит в том, что для аргументов <0,25 * K точность улучшенного деления не намного лучше, чем при использовании наивного метода.