Для простой и эффективной реализации быстрых математических функций с разумной точностью полиномиальные минимаксные аппроксимации часто являются методом выбора. Минимаксные аппроксимации обычно генерируются с использованием варианта алгоритма Ремеза. Различные широко доступные инструменты, такие как Maple и Mathematica, имеют встроенные функции для этого. Сгенерированные коэффициенты обычно вычисляются с использованием высокоточной арифметики. Хорошо известно, что простое округление этих коэффициентов до точности машины приводит к субоптимальной точности в результирующей реализации.
Вместо этого выполняется поиск тесно связанных наборов коэффициентов, которые точно представляются в виде машинных чисел для создания аппроксимации, оптимизированной машиной. Двумя соответствующими документами являются:
Николас Бризебарре, Жан-Мишель Мюллер и Арно Тиссеранд, "Вычислительные машинные эффективные полиномиальные аппроксимации", ACM Transactions on Mathematical Software, Vol. 32, № 2, июнь 2006 г., стр. 236-256.
Николя Бризебарре и Сильвен Чевляр, "Эффективные полиномиальные L∞-аппроксимации", 18-й симпозиум IEEE по компьютерной арифметике (ARITH-18), Монпелье (Франция), июнь 2007 г., стр. 169-176.
Реализация LLL-алгоритма из последней статьи доступна как команда fpminimax() инструмента Sollya. Насколько я понимаю, все алгоритмы, предложенные для генерации оптимизированных машиной аппроксимаций, основаны на эвристике, и поэтому, как правило, неизвестно, какая точность может быть достигнута с помощью оптимального приближения. Мне непонятно, влияет ли доступность FMA (fused multiply-add) для оценки приближения на ответ на этот вопрос. Мне кажется наивно, что он должен.
В настоящее время я рассматриваю простую полиномиальную аппроксимацию для арктангенса на [-1,1], которая оценивается в арифметике с одной точностью IEEE-754, используя схему Хорнера и FMA. См. Функцию atan_poly() в коде C99 ниже. Из-за отсутствия доступа к машине Linux на данный момент я не использовал Sollya для генерации этих коэффициентов, но использовал свою собственную эвристику, которую можно было бы свободно описать как смесь самого крутого приличного и симулированного отжига (чтобы избежать застревания на локальных минимумах), Максимальная ошибка моего машинного оптимизационного многочлена очень близка к 1 ulp, но в идеале я хотел бы, чтобы максимальная ошибка ulp была ниже 1 ulp.
Я знаю, что я мог бы изменить свои вычисления, чтобы повысить точность, например, используя ведущий коэффициент, представленный более чем с точностью до одной точности, но я хотел бы сохранить код в точности как есть (то есть, как просто насколько это возможно), регулируя только коэффициенты, чтобы обеспечить максимально точный результат.
"Проверенный" оптимальный набор коэффициентов был бы идеальным, ссылки на соответствующую литературу приветствуются. Я сделал литературный поиск, но не смог найти какую-либо статью, которая значительно улучшает состояние искусства за пределами Sollya fpminimax(), и ни один из них не изучает роль FMA (если таковая имеется) в этой проблеме.
// max ulp err = 1.03143
float atan_poly (float a)
{
float r, s;
s = a * a;
r = 0x1.7ed1ccp-9f;
r = fmaf (r, s, -0x1.0c2c08p-6f);
r = fmaf (r, s, 0x1.61fdd0p-5f);
r = fmaf (r, s, -0x1.3556b2p-4f);
r = fmaf (r, s, 0x1.b4e128p-4f);
r = fmaf (r, s, -0x1.230ad2p-3f);
r = fmaf (r, s, 0x1.9978ecp-3f);
r = fmaf (r, s, -0x1.5554dcp-2f);
r = r * s;
r = fmaf (r, a, a);
return r;
}
// max ulp err = 1.52637
float my_atanf (float a)
{
float r, t;
t = fabsf (a);
r = t;
if (t > 1.0f) {
r = 1.0f / r;
}
r = atan_poly (r);
if (t > 1.0f) {
r = fmaf (0x1.ddcb02p-1f, 0x1.aee9d6p+0f, -r); // pi/2 - r
}
r = copysignf (r, a);
return r;
}
