Ответ 1

Теория

Сначала обратите внимание на разницу, использующую арифметику с одной точностью.

• [Уравнение 8] Минимальное значение f может быть больше. Поскольку числа с двойной точностью являются супер-множеством чисел с одной точностью, ближайший single до кратного 2/pi может быть только дальше, чем ~ 2.98e-19, поэтому число начальных нулей в фиксированной арифметике представление f должно быть не более 61 ведущих нулей (но, вероятно, будет меньше). Обозначим эту величину fdigits.

• [Уравнение до 9] Следовательно, вместо 121 бит y должен быть точным до fdigits + 24 (отличные от нуля значащие бит в одной точности) + 7 (дополнительные защитные биты) = fdigits + 31 и не более 92.

• [Уравнение 9] "Поэтому вместе с шириной показателя x 2/pi должно содержать 127 (максимальный показатель single) + 31 + fdigits или 158 + fdigits и не более 219 бит.

• [ПОДРАЗДЕЛ 2.5] Размер A определяется числом нулей в x до двоичной точки (и не изменяется на ход до single), а размер C определяется уравнением до 9.

  • Для больших x (x >= 2 ^ 24), x выглядит следующим образом: [24 бит, M нулей]. Умножая его на A, размер которого является первым битом M 2/pi, приведет к целому числу (нули x будут просто смещать все в целые числа).
  • Выбор C для начала с M+d бит 2/pi приведет к тому, что продукт x*C будет иметь размер не более d-24. В двойной точности d выбрано равным 174 (а вместо 24 - 53), так что произведение будет иметь размер не более 121. В single достаточно выбрать d, чтобы d-24 <= 92, точнее, d-24 <= fdigits+31. То есть d можно выбрать как fdigits +55 или не более 116.
  • В результате B должен иметь размер не более 116 бит.

Поэтому мы имеем две проблемы:

• Вычисление fdigits. Это связано с чтением ссылки 6 из связанной бумаги и ее пониманием. Не может быть так просто.:) Насколько я вижу, это единственное место, где используется nearpi.c.
• Вычислить B, соответствующие биты 2/pi. Поскольку M ограничен снизу на 127, мы можем просто вычислить первые 127 + 116 бит 2/pi в автономном режиме и сохранить их в массиве. См. Wikipedia.
• Вычисление y=x*B. Это включает в себя умножение x на 116-битное число. Здесь используется Раздел 3. Размер блоков выбирается равным 24, потому что 2 * 24 + 2 (умножение двух 24-битовых чисел и добавление 3 таких чисел) меньше точности double, 53 (и потому, что 24 делит 96). Мы можем использовать блоки размером 11 бит для арифметики single по тем же причинам.

Примечание. Трюк с B применим только к числам, показатели которых положительны (x >= 2 ^ 24).

Подводя итог - во-первых, вам нужно решить проблему с точностью double. Ваш код Matlab также не работает в точности double (попробуйте удалить single и вычислить sin(2^53), потому что ваш twooverpi имеет только 53 значащих бита, а не 175 (и в любом случае вы не можете напрямую умножать такие точные числа в Matlab). Во-вторых, схема должна быть адаптирована для работы с single, и снова ключевая проблема представляет собой 2/pi достаточно точно и поддерживает умножение высокоточных чисел. Наконец, когда все работает, вы можете попытаться найти лучшее fdigits, чтобы уменьшить количество бит, которое вы должны хранить и умножать.

Надеюсь, что я не совсем ушел - комментарии и противоречия приветствуются.

Пример

В качестве примера давайте вычислим sin(x) где x = single(2^24-1), у которого нет нулей после значимых бит (M= 0). Это упрощает поиск B, поскольку B состоит из первых 116 бит 2/pi. Поскольку x имеет точность 24 бит и B 116 бит, продукт

y = x * B

будет иметь 92 бит точности, если требуется.

В разделе 3 в связанном документе описывается, как выполнять этот продукт с достаточной точностью; тот же алгоритм можно использовать с блоками размера 11 для вычисления y в нашем случае. Будучи тяжелым, я надеюсь, что я извинился за то, что не делал этого явно, вместо этого полагаясь на Matlab символическую математическую панель инструментов. Этот набор инструментов предоставляет нам функцию vpa, которая позволяет нам определять точность числа в десятичных разрядах. Итак,

vpa('2/pi', ceil(116*log10(2)))

будет производить приближение 2/pi не менее 116 бит точности. Поскольку vpa принимает только целые числа для аргумента точности, мы обычно не можем точно определить двоичную точность числа, поэтому мы используем следующий лучший.

Следующий код вычисляет sin(x) в соответствии с документом в точности single:

x = single(2^24-1);
y = x *  vpa('2/pi', ceil(116*log10(2)));    % Precision = 103.075
k = round(y);
f = single(y - k);
r = f * single(pi) / 2;
switch mod(k, 4)
    case 0 
        s = sin(r);
    case 1
        s = cos(r);
    case 2
        s = -sin(r);
    case 3
        s = -cos(r);
end
sin(x) - s                                   % Expected value: exactly zero.

(Точность y получается с помощью Mathematica, который оказался намного лучшим числовым инструментом, чем Matlab:))

В libm

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

См. файл s_sin.c для обертки (таблица 2 из связанной бумаги отображается как оператор switch в конце файла) и e_rem_pio2.c для кода сокращения аргументов (особый интерес представляет массив, содержащий первые 396 шестнадцатеричных цифр 2/pi, начиная с строки 69).