Синус косинус модульная расширенная точность арифметики

Question

Я видел во многих реализациях синуса / косинуса так называемую расширенную модульную арифметику точности. Но для чего это нужно? Например, в реализации cephes, после сокращения до диапазона [0,pi/4], они делают эту модульную арифметику точности, чтобы улучшить точность.

Ниже код:

z = ((x - y * DP1) - y * DP2) - y * DP3;

где DP1, DP2 и DP3 - некоторые жестко закодированные коэффициенты. Как найти этот коэффициент математически? Я понимаю цель "арифметики модульного расширения" для большого числа, но вот каково ее точное назначение?

Answer 1

В контексте сокращения аргументов для тригонометрических функций вы смотрите на сокращение аргументов Коди-Уэйта, метод, представленный в книге: Уильям Дж. Коди и Уильям Уэйт, Руководство по программному обеспечению для элементарных функций, Прентис-Холл, 1980. Цель состоит в том, чтобы добиться для аргументов до определенной величины точного сокращенного аргумента, несмотря на вычитающее аннулирование в промежуточных вычислениях. Для этой цели соответствующая константа представлена ​​с большей, чем собственная точность, с использованием суммы нескольких чисел убывающей величины (здесь: DP1, DP2, DP3), так что все промежуточные продукты, кроме наименее значимого, могут быть вычислены без ошибки округления.

Рассмотрим в качестве примера вычисление sin (113) в IEEE-754. binary32 (одинарная точность). Типичное сокращение аргумента будет концептуально вычислять i=rintf(x/(π/2)); reduced_x = x-i*(π/2), binary32 ближайший к π/2 номер 0x1.921fb6p+0, Мы вычисляем i=72, продукт округляется до 0x1.c463acp+6, что близко к аргументу x=0x1.c40000p+6, Во время вычитания некоторые ведущие биты отменяются, и мы получаем reduced_x = -0x1.8eb000p-4, Обратите внимание на завершающие нули, введенные перенормировкой. Эти нулевые биты не несут полезной информации. Применяя точное приближение к приведенному аргументу, sin(x) = -0x1.8e0eeap-4тогда как истинный результат -0x1.8e0e9d39...p-4, Мы получаем с большой относительной ошибкой и большой ошибкой ulp.

Мы можем исправить это, используя двухэтапное сокращение аргумента Коди-Уэйта. Например, мы могли бы использовать pio2_hi = 0x1.921f00p+0, а также pio2_lo = 0x1.6a8886p-17, Обратите внимание на восемь завершающих нулевых битов в представлении одинарной точностиpio2_hiчто позволяет нам умножать на любое 8-битное целое i и до сих пор есть продукт i * pio2_hi представляется точно как число с одинарной точностью. Когда мы вычисляем ((x - i * pio2_hi) - i * pio2_lo), мы получаем reduced_x = -0x1.8eafb4p-4, и поэтому sin(x) = -0x1.8e0e9ep-4довольно точный результат.

Лучший способ разделить константу на сумму будет зависеть от величины i нам необходимо обработать максимальное количество битов, подлежащих вычитающему аннулированию для заданного диапазона аргументов (основываясь на том, как близкие целые кратные числа π/2 могут получить целые числа), и соображения производительности. Типичные случаи использования в реальной жизни включают двух-четырехступенчатые схемы сокращения Коди-Уэйта. Доступность слитного множественного сложения (FMA) позволяет использовать составляющие константы с меньшим количеством завершающих нулевых битов. Смотрите эту статью: Сильви Болдо, Марк Даумас и Рен-Канг Ли, "Формально проверенное сокращение аргументов с помощью слитного умножения-сложения". IEEE Transactions on Computers, 58: 1139–1145, 2009. Для работающего примера, использующего fmaf() Вы можете посмотреть код в одном из моих предыдущих ответов.