Точные и быстрые вычисления для чисел с плавающей точкой на примере функции синуса. Часть 3: fixed-point

[Refridgerator]

Кто сказал, что значения tsx должны быть равны точно коэффициентам ряда Тейлора. Напрашивается достаточно очевидная идея проварьировать значения коэффициентов, чтобы улучшить точность аппроксимации и убрать постоянную составляющую ошибки. Грамотно сделать такую вариацию достаточно сложно.

А если воспользоваться математическим пакетом — то очень даже просто, потому что там обычно присутствует инструментал для аппроксимации. В WM, в частности, есть функция FindFit, которой можно аппроксимировать что угодно чем угодно (в том числе и рациональным полиномом, который для многих функций даёт большую точность).

[Tiriet]

всего 10 используеых коэффициентов, у каждого можно поменять только три-четыре младших бита максимум, полный перебор потребует 4^10- миллион вариантов. Примерно 2е суток на машине автора перебирать пришлось бы. Если не трогать младшие коэффициенты (например, до пятого включительно), то весь перебор ужмется до трех минут.

[Refridgerator]

Я перебирал 30 коэффициентов и Вольфрам справлялся за несколько секунд, не минут. Методом наименьших квадратов так и ещё быстрее.

[Tiriet]

Мне кажется, вы спутали МНК и градиентный спуск (МНК — это не про поиск максимумов и минимумов функций многих переменных), это раз, и второе- у нас целевая функция- это средняя ошибка округления на большом множестве чисел- она сама по себе вычисляется долго (ТС указал на х86- 0,174с на вычисление одного значения- быстро перебирать варианты не получится), при этом коэффициенты- де-факто целочисленные, и потому к ним градиентный спуск не применим (градиента для целых чисел нету), мы уже имеем решение, близкое к идеальному (где градиентный спуск неэффективен), это решение очень близкое к идеальному- отличающесеся в одном-двух младших битах, поэтому мы даже не сможем вычислить вектор оптимального приращения, и решение возле текущей аппроксимации осциллирует- и градиентный и покоординатный спуски на таких функциях ломаются. Вот и остается чистый перебор множетства значений медленно вычисляемой функции.

[dorex]

Спасибо за комментарий! Всё абсолютно правильно. Более того добавлю, что метод наименьших квадратов использовать нельзя. Нужно минимизировать максимальное значение отклонения. А градиент от этого посчитать ещё более невозможно. В следующих статьях я покажу, как можно аналитически аппроксимировать последний коэффициент. Последние коэффициенты, очевидно можно варьировать больше, чем последние биты. Там есть красивый математических подход.

[Refridgerator]

Мне кажется, вы спутали МНК и градиентный спуск

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

коэффициенты- де-факто целочисленные, и потому к ним градиентный спуск не применим (градиента для целых чисел нету)

Так коэффициенты можно округлять уже после, а ещё лучше — рационализировать. Понятно, что таким образом (равно как и МНК) наилучшего результата не получить — но в качестве первичной аппроксимации сойти вполне может, и результат таки может оказаться лучше Тейлора.

[byman]

«Всё-таки long double достаточно экзотический вид данных и его поддержка в современных процессорах не является приоритетной».
А если в процессоре нет целочисленного умножения 64*64? :) Можно конечно через 32*32,
но это в 10 раз медленнее чем double*double.
Правильно ли я понимаю, что если в сложении-вычитании double можно производить вычисления большей точности за счет double.double, то с умножением такой фокус не проходит? Умножили и хвост безнадежно потерян. Кстати, в процессоре можно легко сделать вычисления double+double c сохранением обоих частей результата. Аналогично можно было бы сделать и с умножением. т.е. по аналогии с целыми числами int*int = long long.
Был бы какой-то толк от подобных инструкций?

[dorex]

1) А если в процессоре нет аппаратной поддержки FP арифметики? Почему вы говорите, что

но это в 10 раз медленнее чем double*double.

Если у Вас есть какой-то конкретный процессор возьмите его и сделайте эксперимент.
2) a) Есть такая суперштука как FMA. С помощью этих комманд можно проводить точное умножение.
b) Если мантисса заполнена на половину, то результат умножения будет точным. Поэтому использование битовых масок решит проблему. Правда понадобится много комманд.

[Refridgerator]

если в сложении-вычитании double можно производить вычисления большей точности за счет double.double, то с умножением такой фокус не проходит?

Проходит, есть алгоритмы для умножения и double на double, и double-double на double-double. Тут в архиве qd-2.3.22.tar.gz есть и описание в pdf, и реализация на си и фортране.

[byman]

Спасибо. Именно об этом и был мой вопрос. Попробую посмотреть, может чего и пойму :)
А то жалко иметь хороший аппаратный double, а расчеты производить на счетах, потому что точнее выходит :)

[byman]

1. я как раз про свой конкретный случай :) в моем процессоре можно 2 double умножения за такт, но ни одного int64*int64. Эксперимент с sin_e7 примером я провел. Реализовал команду умножения int64*int64 через последовательность умножений int32*int32 и сложений.
2. насколько я помню, FMA избегает ненужного округления перед сложением. Результат будет немного точнее. Но у меня нет и FMA :) Только умножить и сложить.
3. Да, я использовал такой подход в синусе. Тот еще геморрой :)

[dorex]

Если Вам интересна практическая реализация то почему бы Вам не адаптировать s_sin? Вам нужна большая точность? Если так то самым простым образом можно решить эту проблему увеличением числа членов в разложении синуса и косинуса:

...
  s = x + (dx + x * xx * (sn3 + xx * sn5)); // здесь
  c = x * dx + xx * (cs2 + xx * (cs4 + xx * cs6)); // и здесь
  SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
...

[dorex]

P.S. По моим прикидкам следующего члена будет достаточно (cs7 и cs8) соответственно. Вопрос, разумеется остаётся к качеству самой таблички. Я, если честно, не проверял.
P.P.S. Проверьте существующие значения cs и sn. Они могут отличатся от ряда Тейлора. Если так, то рассчитайте их сами точно.

[byman]

Спасибо. C s_sin я уже активно работаю :) У меня имеется своя реализация sin-cos (некоторые подробности реализации я изложил в части 2). Интерес к s_sin у меня как к генератору эталонного результата. Также в s_sin мне понравились некоторые приемы работы с double. Моя реализация очень мала по обьему кода и по размеру таблицы (80 байт) и имеет очень хорошую скорость. Длина полинома у меня на два члена длиннее чем в s-sin. Цель — улучшить точность, незначительно потеряв в скорости. Вполне возможно, что я реализую на ассемблере и s_sin, как альтернативный вариант в библиотеке.

[KaminskyIlya]

Друзья, не подскажете какой лучше алгоритм выбрать для вычисления синуса с точностью до 200-300 бит (размер мантиссы). Входной аргумент: 0 <=x <1 в радианах. Может есть хорошие библиотеки по работе с вещественными числами настраиваемой точности (типа, хочешь 200 бит мантиссы - на тебе)? За совет буду благодарен. Или тут можно отделаться арифметикой с фиксированной запятой и работать с рациональным представлением чисел?

Если использовать ряд Тейлора, то не будет ли оптимальнее (для заявленной точности) вычислять sin(x/2) и cos(x/2)? Эти ряды, по идее, будут сходиться быстрее...