- Suslik
- Постоялец
Предлагаю на растерзание триг. функции. Пока для начала только 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.755731919494057310-7 x10+2.08767566054833610-9x12-1.147074217758860810-11 x14+4.77945598607155710-14 x16-1.56182665074719310-16 x18+4.107400301161042310-19 x20-8.8342668522528310-22 x22+1.523029374572779410-24x24-1.71589302693384210-27x26
Для COS1
0.9999999999999852-0.49999999999973677x2+0.04166666666627814x4-0.001388888888671486 x6+0.00002480158723868764 x8-2.755731814260099510-7 x10+2.08767450591165410-9x12-1.147065738517490210-11x14+4.779030207369567610-14 x16-1.560360124714008810-16x18+4.07332466111252510-19x20-8.32382064038449710-22 x22+1.077971929355881410-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
Красота!
В дольнешем планирую свести всё в таблицу по скорости и точности. Есть мысли еще ускорить за счет уменьшения интервала. Ну и.т.д.
Пока всё!
asvp
ну здрасте. насколько я понял, это просто первые несколько членов ряда тейлора. ессно, радиус сходимости у такой последовательности далеко не бесконечность, что выдаст твой косинус, к примеру, на аргумент 1000.0f * pi?
>Красота!
Ещё как красота! я насчитал всего-то 28 обращений к памяти на подсчёт одного косинуса! Автор, ты сравнивал свой уберкод со стандартным sincos? как разница?
Бтв что-то мне подсказывает, что если посчитать тот же ряд на процессоре, без доступа к памяти, получится быстрее.
Для тестов не забываем включать floating point model : fast, считать одновременно синус с косинусом с точностью float.
Suslik
Су, забей ^_^
Тут всё и так понятно. fsincos выдаёт синус в формате double с точностью, вполне приемлимой и для научных рассчётов.
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 и т.д. в зависимости от точности. Конвеер проца быстрее освободиться.
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 точности :)
совет на будущее:
попытки оптимизации самых низкоуровневых операций вроде корня обязательно будут либо менее точны, либо будут с урезанной областью определения, либо медленнее стандартных. исключения практически исключены. лучше оптимизируй алгоритм, который эти твои синусы считает - скорее всего, их можно считать реже и меньше, или вовсе без них обойтись :)
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.
Естьествено для точности я пользуюсь стандартными функциями. А для более простых расчетов я вставляю эти функции. Там где мне не значение важно, а знак или четверть или приблизителное значение.
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.
ясно, спасибо, отписался.
Suslik
Насколько я помню, у косинуса, разложенного в ряд Тейлора, точнее - Маклорена - теоретически радиус сходимости как раз таки бесконечность. Только скорость сходимости уменьшается при увеличении аргумента.
asvp
>Нет. Это не ряд Тейлора.
Выходит, что это ряд, полученный при апрокксимации функции косинуса степенным полиномом? Жжошь.
Вообще тригонометрические функции синус и косинус выражаются через тангенс двойного угла, который очень быстро считается. Гораздо быстрее, чем тут.
Zefick
>Насколько я помню, у косинуса, разложенного в ряд Тейлора, точнее - Маклорена - теоретически радиус сходимости как раз таки бесконечность.
понятия "ряд" и "частичная сумма ряда" - почувствуй разницу. вообще, конечно, понятие радиуса сходимости к частичной сумме ряда отношения мало имеет, но, думаю, мы все поняли, что имелось в виду.
asvp
Ээээ.... Не ну ты всё-таки соберись да проверь свой велосипед на точньсть/быстродействие со стандартным sin/cos от FPU.
А так понятно что непонятно если и велика производительность/точность твоего кода - то насколько?
Ответь на этот вопрос - тогда и будет ясно стоит ли овчинка выделки.
Короче - тестируй.
Тема в архиве.