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

Question

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

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

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

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

Answer 1

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

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

M.Пейн и Р. Ханек.Уменьшение радиана для тригонометрических функций.Информационный бюллетень SIGNUM, 18: 19–24, 1983

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

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

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

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

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

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

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