@Ромуальд_7

Да, конечно, расскажу. Делалось это с использованием MatLab, естественно
Как строится расчёт синуса для точности float (если не CORDIC)? Через интервальные аппроксимации параболами (см. рисунок). Например, точность float (ошибка <1.0E-08) обеспечивается всего лишь ~150 интервалами (если разбиваем от 0 до pi/2), при этом, точности 1.0E-12 достигнуть не удаётся и при числе интервалов в 5000. Более того - точность не растёт при увеличении числа интервалов свыше 3000, а падает. Очевидно, что точности double таким образом достичь не удастся.
Для сходящихся функций (синус, косинус) прекрасно работает ряд Тейлора. Точно так же делим участок [0, pi/2] на N интервалов и раскладываем функцию синуса в ряд Тейлора в точке центра каждого интервала, после чего записываем коэффициенты получившегося полинома (с заранее заданной степенью i[SUB]max[/SUB]) в таблицу:

N	k0 (*x0)	ki (*xi)	kimax (*ximax)
0	2.361868486512384E-02	...	1.651448768486486E-07
1	...	...	...

Мне удалось добиться ошибки в <3.0E-16 (по данным MatLab) на 125 интервалах с полиномом 6 степени.
Отдельно стоит сказать о том, как я нахожу в таблице нужные коэффициенты полинома. Поскольку участок равномерно разбит на N интервалов, можно использовать простую линейную функцию вида y = k*x, где y - номер интервала, а x - аргумент функции sin(x). Отмечу, однако, что мой "бенчмарк" C++ сказал мне, что "MatLab не прав", и выдал максимальную ошибку в 1.2E-14, что, в целом, меня устроило.
С тангенсом я не стал разбираться и бахнул tan(x) = sin(x) / cos(x). Дорого, но я использую эту функцию довольно редко.
Функции же арксинуса и арккосинуса расходятся, поэтому Тейлор тут мдееееее.. Да тут, как оказалось, что угодно будет мдееееее. Оказывается, арк-функции вычисляют через арктангенс, который сходится по формулам, использующим корень. Но разложить функцию арктангенса через ряд Тейлора, как я понял, невозможно, поэтому пользуются другими "аппроксиматорами". Я решил использовать аппроксимацию Паде: степень полинома числителя - 3, знаменателя - 5, число интервалов - 159. Эта функцию определена для любого x, кроме inf и NaN (нужно для inf - бахайте туда x=квинтиллион). Порядок точности получился тем же, что и для синуса с косинусом. Оговорки есть только с интервальным разбиением - здесь не одно равномерное разбиение, а три (x: 0-10, 10-20, 20-40) и пять определённых вручную "широких" интервалов, неравномерно уходящих в "бесконечность".
Ссылок нет, поскольку когда я это всё делал, большей теорией, чем та, что я сейчас написал, я не обладал. Неделя экспериментов, и программная реализация тригонометрии для double готова. В шейдеры, кстати, я загружал таблицы через SSBO; хотел захардкодить, но компилятор оПтИмИзИрУеТ константы, если они кажется ему слишком большими. Ох уж эти 10 Кб информации. Спасибо, умный компилятор.

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

Вернуться к обсуждению:
Непредсказуемые изменения точности расчётов внутри kernel OpenCL

Миниатюры

 

1