Вычисление x по модулю y, где y не может быть представлено как число с плавающей запятой

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

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

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

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


floating-point math language-agnostic modulo
person R.. GitHub STOP HELPING ICE    schedule 25.06.2011    source источник
comment
Итак, как таблица решает вашу проблему? поскольку 2Pi (трансцедентное число) не представимо каким-либо компактным образом, как вы ожидаете, что 2^n mod 2Pi будет представимо? И с чего бы это давало более правильный результат?   -  person zvrba    schedule 25.06.2011
comment
Это можно представить с ошибкой не более 1ulp. С другой стороны, x mod y, где y имеет ошибку не более 1ulp, имеет ошибку до x/y * 1ulp. В частности, если x на MANT_BITS порядков больше, чем y, результат по модулю с использованием аппроксимации для y не имеет никакого значения, т. е. полностью ошибочен.   -  person R.. GitHub STOP HELPING ICE    schedule 25.06.2011
comment
@R ..: Если X на много порядков больше, чем Y, то в скольких реальных ситуациях X mod Y все равно будет иметь смысл? Если X является верхним компонентом суммы Кахана, я вижу, что его точное значение имеет смысл, но в противном случае идея точного вычисления такого модуля выглядит как утверждение, что дыра, измеренная как 25/32 (с точностью до 1/ 64) не следует описывать как имеющий примерно такой же размер, как 2-миллиметровый штифт, а вместо этого следует описывать его как имеющий зазор 0,0061515748. Если бы отверстие было точно 25/32, измерение было бы правильным, но...   -  person supercat    schedule 04.06.2014
comment
... отсутствие информации, указывающей, что размер отверстия на самом деле точно кратен 1/64, любая предполагаемая точность вычисленной разницы в размерах не имеет смысла.   -  person supercat    schedule 04.06.2014

Ответы (1)


arrow_upward
10
arrow_downward

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

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

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

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

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

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

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

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

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

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

person njuffa    schedule 26.06.2011