Вычисление x mod y, где y не представляется в виде с плавающей запятой

В качестве канонического примера рассмотрим проблему уменьшения аргумента для тригонометрических функций, как при вычислении x mod 2π в качестве первого шага вычисления sin (x). Эта проблема затруднена тем, что вы не можете просто использовать fmod, потому что y (2π в примере) не представляется.

Я придумал простое решение, которое работает для произвольных значений y, а не только 2π, и мне любопытно, как он сравнивает (по производительности) с типичными алгоритмами сокращения аргументов.

Основная идея состоит в том, чтобы хранить таблицу, содержащую значение 2 n mod y для каждого значения n в диапазоне log2 (y), до максимально возможного показателя с плавающей запятой, а затем используя линейность модульной арифметикой, суммируйте значения в этой таблице над битами, которые установлены в значении x. Он равен N ветвям и не более N дополнений, где N - количество бит мантиссы в вашем типе с плавающей точкой. Результат не обязательно меньше y, но он ограничен N * y, и процедура может быть применена снова, чтобы дать результат, ограниченный log2 (N) * y или fmod, можно просто использовать в этой точке с минимальным ошибка.

Можно ли улучшить это? И типичные алгоритмы сокращения тригонометрических аргументов работают для любого y или только для 2π?

Ответы

Ответ 1

Современные реализации тригонометрических функций в математических библиотеках корректно работают во всем домене ввода. Они делают это, представляя некоторую константу, связанную с π, такую ​​как 2/π, с достаточной точностью для используемого формата с плавающей запятой.

Например, для сокращения триггерной функции в двойной точности IEEE необходимо представить константу примерно до 1150 бит, которая будет использоваться в том, что по сути является вычислением с фиксированной точкой. Этот подход был впервые предложен авторами следующей работы:

М. Пейн и Р. Ханек. Уменьшение радиан для тригонометрических функций. SIGNUM Newsletter, 18: 19-24, 1983

Оригинальная идея с тех пор была изменена и уточнена другими авторами; возможны варианты с плавающей точкой и с целым числом. Библиотека FDLIBM предоставляет здесь полностью обработанный пример:

http://www.netlib.org/fdlibm/e_rem_pio2.c

В следующей статье автора FDLIBM описан подход, используемый в этом коде

http://www.validlab.com/arg.pdf К.С. Ng. Уменьшение аргумента для огромных аргументов: полезно для последнего бит

Обратите внимание, что нет необходимости переносить промежуточное вычисление на 1150 бит. Так как при сокращении ведущие биты отменяют вычисление, нужно только включать гораздо меньшую группу бит внутри полной константы. Из-за необходимости арифметики с несколькими точками это по-прежнему довольно дорогостоящая операция.

Для уменьшения аргумента тригонометрической функции в более жестко ограниченных диапазонах возможны другие, более экономичные схемы, особенно если аппаратное обеспечение поддерживает FMA (плавное многократное добавление).

Методы, используемые для сокращения тригонометрических аргументов, кажутся доступными для сокращений с помощью произвольных высокоточных констант.