Войти
ПрограммированиеФорумОбщее

Быстрые и точные тригонометрические функции

Advanced: Тема повышенной сложности или важная.

Страницы: 1 2 3 4 5 6 Следующая »
#0
18:40, 24 мар. 2009

Предлагаю на растерзание триг. функции. Пока для начала только COS.

F64 COS0(const F64& v)
{
  PBL_ALIGN(16) F64 Sq=v*v;
  return (((((((((((((-1.715893026933842e-27 * Sq)+
                    1.5230293745727794e-24)  * Sq-
                    8.83426685225283e-22)  * Sq+
                    4.1074003011610423e-19)  * Sq-
                    1.561826650747193e-16) * Sq+
                    4.779455986071557e-14) * Sq-
                    1.1470742177588608e-11) * Sq+
                    2.087675660548336e-9)  * Sq-
                    2.7557319194940573e-7)  * Sq+
                    0.00002480158730015702) * Sq-
                    0.0013888888888846677)  * Sq+
                    0.04166666666666017)    * Sq-
                    0.4999999999999963)    * Sq+
                    1.0000000000000024;
}

F64 COS1(const F64& v)
{
  PBL_ALIGN(16) F64 Sq=v*v;
  return ((((((((((((1.0779719293558814e-24 * Sq)-
                    8.323820640384497e-22)  * Sq+
                    4.073324661112525e-19)  * Sq-
                    1.5603601247140088e-16) * Sq+
                    4.7790302073695676e-14) * Sq-
                    1.1470657385174902e-11) * Sq+
                    2.087674505911654e-9)  * Sq-
                    2.7557318142600995e-7)  * Sq+
                    0.00002480158723868764) * Sq-
                    0.001388888888671486)  * Sq+
                    0.04166666666627814)    * Sq-
                    0.49999999999973677)    * Sq+
                    0.9999999999999852;
}

Разница в функциях скорость и точность. Точность порядка <=6*10^-14.
Полиномы:
Для COS0
1.0000000000000024-0.4999999999999963*x2+0.04166666666666017x4-0.0013888888888846677x6+0.00002480158730015702x8-2.755731919494057310-7 x10+2.08767566054833610-9x12-1.147074217758860810-11 x14+4.77945598607155710-14 x16-1.56182665074719310-16 x18+4.107400301161042310-19 x20-8.8342668522528310-22 x22+1.523029374572779410-24x24-1.71589302693384210-27x26

Для COS1
0.9999999999999852-0.49999999999973677x2+0.04166666666627814x4-0.001388888888671486 x6+0.00002480158723868764 x8-2.755731814260099510-7 x10+2.08767450591165410-9x12-1.147065738517490210-11x14+4.779030207369567610-14 x16-1.560360124714008810-16x18+4.07332466111252510-19x20-8.32382064038449710-22 x22+1.077971929355881410-24x24

Для более быстрых и менее точных функций требуется пересчет полиномов.
Да, чуть не забыл. Прелесть всей этой фигни, в том, что интервал v от 0 до 2PI. Хотя и при интервале от -2PI до 2PI тоже работают.

Вот код на asme для COS1 после оптимизации. (VS2008, Intel C++, SSE2)
10002CDC  movsd      xmm0,mmword ptr [a]
10002CE4  mulsd      xmm0,mmword ptr [a]
10002CEC  movsd      mmword ptr [ebp-0D8h],xmm0
10002CF4  movsd      xmm1,mmword ptr [___xt_z+108h (1002A370h)]
10002CFC  mulsd      xmm1,mmword ptr [ebp-0D8h]
10002D04  subsd      xmm1,mmword ptr [___xt_z+110h (1002A378h)]
10002D0C  mulsd      xmm1,mmword ptr [ebp-0D8h]
10002D14  addsd      xmm1,mmword ptr [___xt_z+118h (1002A380h)]
10002D1C  mulsd      xmm1,mmword ptr [ebp-0D8h]
10002D24  subsd      xmm1,mmword ptr [___xt_z+120h (1002A388h)]
10002D2C  mulsd      xmm1,mmword ptr [ebp-0D8h]
10002D34  addsd      xmm1,mmword ptr [___xt_z+128h (1002A390h)]
10002D3C  mulsd      xmm1,mmword ptr [ebp-0D8h]
10002D44  subsd      xmm1,mmword ptr [___xt_z+130h (1002A398h)]
10002D4C  mulsd      xmm1,mmword ptr [ebp-0D8h]
10002D54  addsd      xmm1,mmword ptr [___xt_z+138h (1002A3A0h)]
10002D5C  mulsd      xmm1,mmword ptr [ebp-0D8h]
10002D64  subsd      xmm1,mmword ptr [___xt_z+140h (1002A3A8h)]
10002D6C  mulsd      xmm1,mmword ptr [ebp-0D8h]
10002D74  addsd      xmm1,mmword ptr [___xt_z+148h (1002A3B0h)]
10002D7C  mulsd      xmm1,mmword ptr [ebp-0D8h]
10002D84  subsd      xmm1,mmword ptr [___xt_z+150h (1002A3B8h)]
10002D8C  mulsd      xmm1,mmword ptr [ebp-0D8h]
10002D94  addsd      xmm1,mmword ptr [___xt_z+158h (1002A3C0h)]
10002D9C  mulsd      xmm1,mmword ptr [ebp-0D8h]
10002DA4  subsd      xmm1,mmword ptr [___xt_z+160h (1002A3C8h)]
10002DAC  mulsd      xmm1,mmword ptr [ebp-0D8h]
10002DB4  addsd      xmm1,mmword ptr [___xt_z+168h (1002A3D0h)]
10002DBC  mov        eax,dword ptr [this]
10002DC2  movsd      mmword ptr [eax+180h],xmm1

Красота!
В дольнешем планирую свести всё в таблицу по скорости и точности. Есть мысли еще ускорить за счет уменьшения интервала. Ну и.т.д.
Пока всё!


#1
18:58, 24 мар. 2009

asvp
ну здрасте. насколько я понял, это просто первые несколько членов ряда тейлора. ессно, радиус сходимости у такой последовательности далеко не бесконечность, что выдаст твой косинус, к примеру, на аргумент 1000.0f * pi?
>Красота!
Ещё как красота! я насчитал всего-то 28 обращений к памяти на подсчёт одного косинуса! Автор, ты сравнивал свой уберкод со стандартным sincos? как разница?

#2
19:02, 24 мар. 2009

Бтв что-то мне подсказывает, что если посчитать тот же ряд на процессоре, без доступа к памяти, получится быстрее.

#3
19:04, 24 мар. 2009

Для тестов не забываем включать floating point model : fast, считать одновременно синус с косинусом с точностью float.

#4
19:08, 24 мар. 2009

Suslik
Су, забей ^_^
Тут всё и так понятно. fsincos выдаёт синус в формате double с точностью, вполне приемлимой и для научных рассчётов.

#5
19:17, 24 мар. 2009

Suslik
1. Нет. Это не ряд Тейлора.
Я же написал: "интервал v от 0 до 2PI. Хотя и при интервале от -2PI до 2PI тоже работают".
Ряд просчитан именно для этого интервала. А преобразование интервала 1000.0*PI к 2PI до 2PI делай сам.

2. Float -32-x битный. При расчетах офигенная погрешность и потеря точность + накопление ошибок. Для float тоже планирую, но нам полиномы будут гораздо меньше. Точность предтавления чисел не позволит.

3. "Бтв что-то мне подсказывает, что если посчитать тот же ряд на процессоре, без доступа к памяти, получится быстрее." Фразу не понял. Это как без доступа к памяти?

4. Насчет fsincos. Согласен. Но эта функция выполняется на FPU. А эти функции SSE2. Т.о. ты часть можешь выполнять через fsincos, если тебе нужна точность, а часть, если тебе точность не важна, SINCOS0 или SINCOS1, SINCOS2 и т.д. в зависимости от точности. Конвеер проца быстрее освободиться.

#6
19:26, 24 мар. 2009

asvp
>Нет. Это не ряд Тейлора.
o'rly?
cos(x) = 1 - x2 / 2! + x4 / 4! + o(x5) - это, должно быть, называется теперь "ряд asvp'а"? :)

>Ряд просчитан именно для этого интервала. А преобразование интервала 1000.0*PI к 2PI до 2PI делай сам.
зачем мне делать сомнительной скорости преобразования, когда есть гарантированно быстрое точное, сделанное давным-давно за меня <math.h>::cosf(float) ? а вот то, что преобразование из интервала [-inf +inf] -> [0.0f, 2.0f * pi] удастся сделать достаточно быстро - совершенно не очевидно.

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

>Но эта функция выполняется на FPU. А эти функции SSE2.
может быть хватит изобретать сферические велосипеды в вакууме? просто замени в любой здоровой программе стандартный sinf на свой UBERSIN99, да сравни скорость, если не веришь синтетическим тестам. обещаю, math.h разрабатывали не дураки, и если что-то было возможно без ограничения общности ускорить, это это уже, наверное, ускорили.

>Т.о. ты часть можешь выполнять через fsincos, если тебе нужна точность
то считаешь через hsincos, который считает в long double точности :)

#7
19:27, 24 мар. 2009

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

#8
19:46, 24 мар. 2009

Suslik

Asm функции SINCOS подставляемая автоматически при Release. (Это стандартаня мат. библиотек Intel C++). А я уж думаю они не зря не пишут просто fsincos.

10010680  push        ebp 
10010681  mov        ebp,esp
10010683  sub        esp,68h
10010686  pextrw      eax,xmm0,3
1001068B  and        eax,7FFFh
10010690  sub        eax,3030h
10010695  cmp        eax,10C5h
1001069A  ja          ___libm_sse2_sincos+16Ch (100107ECh)
100106A0  unpcklpd    xmm0,xmm0
100106A4  movapd      xmm1,xmmword ptr ds:[1002A8E0h]
100106AC  mulpd      xmm1,xmm0
100106B0  movapd      xmm2,xmmword ptr ds:[1002A8F0h]
100106B8  cvtsd2si    edx,xmm1
100106BC  addpd      xmm1,xmm2
100106C0  movapd      xmm3,xmmword ptr ds:[1002A900h]
100106C8  subpd      xmm1,xmm2
100106CC  movapd      xmm2,xmmword ptr ds:[1002A910h]
100106D4  mulpd      xmm3,xmm1
100106D8  add        edx,1C7600h
100106DE  movapd      xmm4,xmm0
100106E2  and        edx,3Fh
100106E5  movapd      xmm5,xmmword ptr ds:[1002A8D0h]
100106ED  lea        eax,ds:[100298A0h]
100106F3  shl        edx,6
100106F6  add        eax,edx
100106F8  mulpd      xmm2,xmm1
100106FC  subpd      xmm0,xmm3
10010700  mulpd      xmm1,xmmword ptr ds:[1002A920h]
10010708  subpd      xmm4,xmm3
1001070C  movapd      xmm7,xmmword ptr [eax+10h]
10010711  movapd      xmm3,xmm4
10010715  subpd      xmm4,xmm2
10010719  mulpd      xmm5,xmm0
1001071D  subpd      xmm0,xmm2
10010721  movapd      xmm6,xmmword ptr ds:[1002A8B0h]
10010729  mulpd      xmm7,xmm4
1001072D  subpd      xmm3,xmm4
10010731  mulpd      xmm5,xmm0
10010735  mulpd      xmm0,xmm0
10010739  subpd      xmm3,xmm2
1001073D  movapd      xmm2,xmmword ptr [eax]
10010741  subpd      xmm1,xmm3
10010745  movapd      xmm3,xmmword ptr [eax+30h]
1001074A  addpd      xmm2,xmm3
1001074E  subpd      xmm7,xmm2
10010752  mulpd      xmm1,xmm7
10010756  movapd      xmm7,xmmword ptr [eax+10h]
1001075B  mulpd      xmm2,xmm4
1001075F  mulpd      xmm6,xmm0
10010763  mulpd      xmm3,xmm4
10010767  mulpd      xmm2,xmm0
1001076B  mulpd      xmm7,xmm0
1001076F  mulpd      xmm0,xmm0
10010773  addpd      xmm5,xmmword ptr ds:[1002A8C0h]
1001077B  mulpd      xmm4,xmmword ptr [eax]
1001077F  addpd      xmm6,xmmword ptr ds:[1002A8A0h]
10010787  mulpd      xmm5,xmm0
1001078B  movapd      xmm0,xmm3
1001078F  addpd      xmm3,xmmword ptr [eax+10h]
10010794  addpd      xmm6,xmm5
10010798  movapd      xmm5,xmm6
1001079C  unpckhpd    xmm6,xmm6
100107A0  unpcklpd    xmm5,xmm5
100107A4  mulpd      xmm6,xmm7
100107A8  mulpd      xmm2,xmm5
100107AC  movapd      xmm7,xmm4
100107B0  addpd      xmm4,xmm3
100107B4  movapd      xmm5,xmmword ptr [eax+10h]
100107B9  subpd      xmm5,xmm3
100107BD  subpd      xmm3,xmm4
100107C1  addpd      xmm1,xmmword ptr [eax+20h]
100107C6  addpd      xmm5,xmm0
100107CA  addpd      xmm3,xmm7
100107CE  addpd      xmm1,xmm5
100107D2  addpd      xmm1,xmm3
100107D6  addpd      xmm1,xmm2
100107DA  addpd      xmm1,xmm6
100107DE  addpd      xmm1,xmm4
100107E2  movapd      xmm0,xmm1
100107E6  unpckhpd    xmm1,xmm1
100107EA  jmp        ___libm_sse2_sincos+1E0h (10010860h)
100107EC  jg          ___libm_sse2_sincos+1B4h (10010834h)
100107EE  movapd      xmm1,xmm0
100107F2  pextrw      eax,xmm0,3
100107F7  or          eax,8000h
100107FC  pinsrw      xmm1,eax,3
10010801  addsd      xmm1,mmword ptr ds:[1002A960h]
10010809  cmp        eax,8010h
1001080E  jge        ___libm_sse2_sincos+19Ah (1001081Ah)
10010810  mulsd      xmm0,mmword ptr ds:[1002A950h]
10010818  jmp        ___libm_sse2_sincos+1E0h (10010860h)
1001081A  movsd      xmm3,mmword ptr ds:[1002A930h]
10010822  mulsd      xmm3,xmm0
10010826  subsd      xmm0,xmm3
1001082A  mulsd      xmm0,mmword ptr ds:[1002A940h]
10010832  jmp        ___libm_sse2_sincos+1E0h (10010860h)
10010834  sub        esp,20h
10010837  movsd      mmword ptr [esp],xmm0
1001083C  lea        eax,[esp+20h]
10010840  mov        dword ptr [esp+8],eax
10010844  mov        eax,3
10010849  mov        dword ptr [esp+0Ch],eax
1001084D  call        ___libm_sincos_huge (10011170h)
10010852  add        esp,20h
10010855  movq        xmm0,mmword ptr [esp+8]
1001085B  movq        xmm1,mmword ptr [esp]
10010860  mov        esp,ebp
10010862  pop        ebp 
10010863  ret

Для инфы.
FSINCOS выполняется 170-250 тактов. Следующую такющу команду он может выполнять только через 140 таков.

Я так полагаю ты вообще не грамотен в математике. Этот ряд основан на "Экономичной рациональной аппроксимации + регрессия"

<math.h>::cosf(float) работает с 32-битными вещественными числами. Я посмотрю на твои расчеты, когда тебе необходимо будет вычислять метры в масштабе тыс. км.

>"открой секрет - где тебе не хватило точности стандартного дабл-косинуса, что ты заметил накопление ошибки в нём?"
То ты про cosf, то ты "дабл-косинуса". ???

Я велосипед не изобретаю. Во всех книгах по аппроксимации и полиномам я не нашел формулу для интервала от 0 до 2Pi. Там везде от 0 до Pi/2 или от 0 до Pi/4.
Естьествено для точности я пользуюсь стандартными функциями. А для более простых расчетов я вставляю эти функции. Там где мне не значение важно, а знак или четверть или приблизителное значение.

#9
19:56, 24 мар. 2009

0.9984083758964231- 0.4969997328713499 x^2 + 0.040636129463344 x^4 -
0.0012687183602229071 x^6 + 0.000018410314478490932 x^8 -
1.0811929725975097*10^-7 x^10

Вот полином с точностью ~0.004
в нем 6 умножений и 5 слож. и вычит.
MULSD выпоняется за 5 тактов, следующая команада может выполняться через 1-2 такт.
ADDPD выпоняется за 3 тактов, следующая команада может выполняться через 1-2 такт.
Итого пусть хоть 6*5+3*5=45 тактов. Да пусть хоть 100. Мне наплевать, когда мне всего лишь нужно определить четверть угла без приведения к интервалу от 0 до Pi/2 или того хуже Pi/4.

#10
20:05, 24 мар. 2009

ясно, спасибо, отписался.

#11
20:05, 24 мар. 2009

Suslik
  Насколько я помню, у косинуса, разложенного в ряд Тейлора, точнее - Маклорена - теоретически радиус сходимости как раз таки бесконечность. Только скорость сходимости уменьшается при увеличении аргумента.

asvp
>Нет. Это не ряд Тейлора.
  Выходит, что это ряд, полученный при апрокксимации функции косинуса степенным полиномом? Жжошь.

  Вообще тригонометрические функции синус и косинус выражаются через тангенс двойного угла, который очень быстро считается. Гораздо быстрее, чем тут.

#12
20:08, 24 мар. 2009

Zefick
>Насколько я помню, у косинуса, разложенного в ряд Тейлора, точнее - Маклорена - теоретически радиус сходимости как раз таки бесконечность.
понятия "ряд" и "частичная сумма ряда" - почувствуй разницу. вообще, конечно, понятие радиуса сходимости к частичной сумме ряда отношения мало имеет, но, думаю, мы все поняли, что имелось в виду.

#13
20:13, 24 мар. 2009
Изображение
#14
20:19, 24 мар. 2009

asvp

Ээээ.... Не ну ты всё-таки соберись да проверь свой велосипед на точньсть/быстродействие со стандартным sin/cos от FPU.
А так понятно что непонятно если и велика производительность/точность твоего кода - то насколько?
Ответь на этот вопрос - тогда и будет ясно стоит ли овчинка выделки.

Короче - тестируй.

Страницы: 1 2 3 4 5 6 Следующая »
ПрограммированиеФорумОбщее

Тема в архиве.