Программирование игр, создание игрового движка, OpenGL, DirectX, физика, форум
GameDev.ru / Программирование / Форум / Вычислить Sin и Cos (3 стр)

Вычислить Sin и Cos (3 стр)

Страницы: 1 2 3 4 58 Следующая »
FordPerfectПостоялецwww9 апр. 201819:20#30
Мысли вслух.

Вычисление аппроксимаций функций обычно балансирует между несколькими факторами. Основной баланс обычно - "скорость/точность", но есть и другие моменты, например: затраты памяти (таблицы), дополнительная логистика (код по инициализации этих таблиц). Плюс "скорость" - тоже не столь однозначна. Например функции для вычисления одного синуса и пачки синусов (на SIMD) могут выглядеть весьма по-разному.

Для баланса "скорость/точность" - есть граница: баланс можно выкрутить до упора - "точнее некуда". В терминах IEEE-754 это называется correctly-rounded result - ближайшее (в выбранном направлении округления) представимое значение выбранного типа.
Задача correctly-rounded cos/sin решена для float (float32) и double (float64). Например, математическая библиотека GCC предоставляет эти функции.
Исходный код:
float (несколько файлов). Вообще-то не уверен, correctly-rounded ли здесь.
double

По-хорошему их надо бы сопровождать доказательствами, но не нашёл. Возможно, IBM подарили только код.

Кстати аппаратная функция x87, FSIN/FCOS, не является correctly-rounded (и вообще, местами сильно неточная), и медленнее ручной sinf():
https://randomascii.wordpress.com/2014/10/09/intel-underestimates… -quintillion/

Одна из причин - то, что период cos/sin - иррациональный. Это затрудняет range reduction (см. напр. http://www.imada.sdu.dk/~kornerup/papers/RR2.pdf ).
Собственно, примерно по этой причине IEEE-754 рекомендует реализовывать sinPi, cosPi. У которых период целый.

Если не требуется correctly-rounded, то вопрос - какой точности достаточно?
Это может быть несколько ulp (units in the last place). Вот пример библиотеки, которая считает так:
http://www.netlib.org/fdlibm/

Старая годная дока:
http://www.research.scea.com/research/pdfs/RGREENfastermath_GDC02.pdf
Там озвучивается следующая мысль. Если считаете аппроксимацию полиномом - имеет смысл убедиться, что коэффициенты точно представимы в выбранном типе; иначе если например первый коэффициент неточный, то это означает, что остальные посчитаны исходя из неверных предположений. На этом легко потерять несколько ulp.

Как считать.
Вычисление обычно распадается на 2 этапа: range reduction и polynomial approximation.
По-поводу range reduction была ссылка выше. Ну и вообще она бывает нетривиальной. Поэтому зачастую имеет смысл брать функции, где она тривиальна - sinPi, cosPi, pow2. Если хочется чуть более точной, но не совсем правильной range reduction, то там есть свои трюки (разбиение периода на несколько констант, например).

Приближение полиномом.
По Тейлору (и Маклорену). Основная мысль - за исключением случаев, когда точно знаете, что делаете,- использование Тейлора в полиномиальном приближении - сомнительная идея. Тейлор хорош тем, что его легко считать, и легко оценивать точность, поэтому он удобен для референсных значений, с арифметикой произвольной точности.
Сами функции обычно лучше считать по-другому.
Как именно?
Начать можно тривиально с полиномиальной интерполяции. Посчитать (достаточно точно, хоть бы и Тейлором) значение функции в нескольких точках и построить полином, проходящий через эти точки. У него есть свои проблемы (феномен Рунге), и точность неидеальная, но на практике вполне может быть достаточно.
Если хочется большего - есть другие методы.
Я временами использовал минимизацию L2-нормы отклонения от требуемой функции. Её преимущество в том, что для многих функций её можно посчитать аналитически.
Но приличные люди используют минимизацию L. Для этого используется алгоритм Ремеза.
Если тянуть алгоритм Ремеза влом, можно использовать разложение по базису многочленов Чебышёва.

Для подобных вещей существует специальная тулза (доказывающая оценки ошибок, умеющая correctly-rounded и вообще учитывающая плавучку):
http://sollya.gforge.inria.fr

Если требования к точности ниже - можно идти дальше.
По скорости таблицу обогнать сложно (если не float sin(float x) {return 0.0f;}).
Но это для одного значения. Для нескольких одновременно - таблица может хуже SIMD-ифицироваться.
Размер таблицы может быть небольшим.

Вот здесь:
High-Speed Function Approximation Using a Minimax Quadratic Interpolator
говорят, что таблицы на 64 набора квадратичных коэффициентов достаточно для 1 ulp, а кубических - для correctly-rounded.

Точности 10-3 можно достичь нечётным полиномом 5-й степени (3 умножения, 4 сложения).

Мои темы:
https://gamedev.ru/code/forum/?id=204567
https://gamedev.ru/flame/forum/?id=220141

Правка: 9 апр. 2018 19:24

FordPerfectПостоялецwww9 апр. 201820:13#31
FordPerfectПостоялецwww9 апр. 201822:02#32
eDmk
Проблема ряда Тейлора - что он очень неравномерно распределяет погрешность ( O((x-x0)n) ) - очень малая погрешность вблизи точки, и резко улетает с удалением, собственно о чём Suslik говорил.
FordPerfectПостоялецwww9 апр. 201822:33#33
Пример:
http://rextester.com/MBUN65239
Как можно заметить, на [-pi/2; pi/2] у Тейлора ошибка в 100 раз меньше (~27), чем на [-pi; pi].
FordPerfectПостоялецwww9 апр. 201822:57#34
Suslik
В #18 ещё вычисление через степень забавное.

eDmk
> https://albertveli.wordpress.com/2015/01/10/zx-sine/
А где там таблицы?

Так-то Тейлором можно хоть sin(20000) вычислять. Благо радиус сходимости это позволяет. Если точности хватает.

Lion007Постоялецwww9 апр. 201823:36#35
каирады, нисколько не пытаясь прервать увлекательной дискуссии о методах вычисления синусов и косинусов...
а про банальный сопроцессор забыли? а он ведь - как для родных, одной инструкцией и синус и косинус считал! :) http://www.felixcloutier.com/x86/FSINCOS.html
кто его знает, что там - внутри... но, между прочим, хотя это и не впоне очевидно - сопроцессор считает лучше, чем повсеместно вытеснившее его sse. медленнее, но точнее - у него числа десятибайтовые, а не восьмибайтовые даблы... факт известный, но, как мне кажется, незаслуженно забытый.
eDmkУчастникwww10 апр. 20180:04#36
FordPerfect
Переделал на Delphi.
Можно в коллекцию или это не ваше произведение?
+ Показать

FordPerfect
>А где там таблицы?

Разве не таблицы?

+ Показать

Правка: 10 апр. 2018 0:26

eDmkУчастникwww10 апр. 20180:05#37
>Lion007
Тут заговор против студентов.

Правка: 10 апр. 2018 1:01

FordPerfectПостоялецwww10 апр. 20180:32#38
eDmk
> Может подскажите что не так в коде?
В коде на Python используются целые числа произвольной точности.

Но код всё-равно получился неожиданно забавным.
Потому что внимание.
Чего я хотел: посчитать Тейлора в fixed point с 14 знаками после точки (и дофига - до).
Если бы я считал отдельно числитель и знаменатель каждого члена ряда Тейлора - всё бы так и работало.
Ошибка (абсолютная) на каждом делении была бы порядка 10-14, и на всём суммировании была бы абсолютная ошибка порядка 100000*10-14=10-9.
Ок.
Но я считаю тупо накапливая множитель.
И соответственно уже на первом умножении получаю floor(200003*1014/6) вместо 200003*1014/6.
В смысле, ошибка порядка 10-14, которая дальше вносится во все остальные члены, некоторые из которых больше 10100, т. е. эта ошибка изменяет значение ряда больше, чем на единицу.

А оно - всё равно работает.
Вопрос - какого фига?!

У меня есть некоторые гипотезы, но пока не очевидно.

Вот добавил более влобную версию (и считаю sin(2000), чтоб во время влезть):
http://rextester.com/UDQV85800

Правка: 10 апр. 2018 0:41

eDmkУчастникwww10 апр. 20180:53#39
FordPerfect
Данная версия достаточно точная. Ошибка расхождения наблюдаются в 11 разряде и то не везде.
По большей части совпадает с процессором. Но она очень медленная.

Sin_Teylor(30) = 0,49999999999999 - Тоже немного напрягает. Должно быть 0.5

Ваш вариант на Delphi:

+ Показать

Правка: 10 апр. 2018 0:55

FordPerfectПостоялецwww10 апр. 20181:55#40
> http://www.research.scea.com/research/pdfs/RGREENfastermath_GDC02.pdf
Robin Green - Faster Math Functions.
Ссылка не рабочая. Есть по адресу http://www.ue.eti.pg.gda.pl/~wrona/lab_dsp/cw04/fun_aprox.pdf .
Реально вкусная дока - годное введение в тему для широкого круга читателей.

eDmk
> Ваш вариант на Delphi
Оно достаточно отличается от моего. И не особо осмысленно в таком виде.
Поинт моего кода был - прямо посчитать сумму ряда Тейлора вдали от 0. Чего версия с double всё-равно не делает.
Без этого же считать синус (даже Тейлором) можно сильно по-другому.

> Разве не таблицы?
Так можно до того дойти, что любые константы в коде таблицами называть.

Lion007
Я ж выше давал ссылку:
https://randomascii.wordpress.com/2014/10/09/intel-underestimates… -quintillion/
У встроенной FSIN/FCOS/FSINCOS существенные проблемы с точностью на некоторых значениях.

FordPerfectПостоялецwww10 апр. 20183:55#41
eDmkУчастникwww10 апр. 20184:11#42
>любые константы в коде таблицами называть
Ну да. Массив с константами и есть таблица :)
SuslikМодераторwww10 апр. 20185:18#43
eDmk
> 0,49999999999999 - Тоже немного напрягает. Должно быть 0.5
кому должно?

мне кажется, имеет смысл почитать статью от FordPrefect, чтобы сложить хотя бы общее представление о том, как работают числа с плавающей точкой: https://gamedev.ru/code/articles/FloatingPoint

eDmkУчастникwww10 апр. 20185:51#44
>кому должно?
У слова должно 2 ударения и смысла. В моем случае дОлжно. Синоним - полагается.
Страницы: 1 2 3 4 58 Следующая »

/ Форум / Программирование игр / Общее

2001—2018 © GameDev.ru — Разработка игр