Портирование FPU на SSE

пч фпу кривой и его лишь для обратной совместимости держат. видать, не стали мутить приличную хард реализацию, чтоб сохранить место на плашке под "широкие" функи проца

 

UbIvItS,

Интересна сама последовательность в хард, те какой ряд использован чебышев ?

Выше был пример и сказали что результат совпадает с хард, но я не верю. Вероятно товарищ тупо округлил значения и они сошлись. Тем более там дельфи, сомнительно там всё - вот выше не ясно что за синус, может вначале вычисляется косинус, сохраняется в те переменные.

Та либа Интел по математике - почему используется fpu но не используются непосредственно готовые инструкции тот же fsin не понимаю.

 

хз, что там пользовали (можЬ в манах Ынтеля написано). но в общем и целом для харда лучше кордик юзать == схема более компактная и можно хорошо балансировать погрешность.

 

UbIvItS,

> харда лучше кордик юзать

Что такое кордик" ?

 

https://wasm.in/threads/portirovanie-fpu-na-sse.33827/page-2#post-420627

 

UbIvItS,

Кордик" какой то алго ссылка на вики. Всё равно это полином, их много эффективность определяет как быстро сходится ряд. Чебышев быстрее чем Тейлор.

 

Я убрал fwait, результат тот же. По крайней мере, изменения чисел не заметны на глаз.

Нет, это Тейлор, см. ниже.

Смотри исходник, там же всё есть. Это делители:
1/4!, -1/6!, 1/8!, -1/10!, 1/12!, 1/14!, а 1/2! = 0.5 вшито в исходник. Я так понял, сделано в бинарном виде для повышения точности.

 

Jin X,

Где сурец посмотреть выложи сюда, я дельфи не юзаю так что хз. Замечательно если можно посмотреть алго в сурке.

 

По поводу точности.
Из 36000000 вычислений 1139430 несовпадений fpu и sse (3,165%), суммарно ошибка составляет 8.08899921015606E-0021.
Под несовпадением я имею в виду разницу значений хотя бы в один бит.

Больше, чем на единицу отличается 417 значений (если реинтерпрет-кастить их к uint64),
больше, чем на 2 - 241,
больше, чем на 3 - 175,
больше, чем на 4 - 144,
больше, чем на 8 - 87,
больше, чем на 16 - 43,
больше, чем на 32 - 22,
больше, чем на 64 - 11,
больше, чем на 128 - 6,
больше, чем на 256 - 3,
больше, чем на 512..8192 - 2,
больше, чем на 16384 - 1,
больше, чем на 32768 - 0.

Причём, все значения, которые отличаются больше, чем на 1 (417 штук которые) – при аргументе близком к Пи (и кратным Пи, кроме 0).
В частности, в диапазоне от 3.14156263286791 до 3.14162110139778.
Т.е. в диапазоне от 0 до почти Пи ошибка составляет не более, чем на 1 (при рассмотрении чисел как uint64).
Максимальная разница – при X = 3.14159265253101, получается fpu = 1.05878562617649E-0009 , sse = 1.05878562618053E-0009 (разница составляет 4.04449960507803E-0021, т.е. половину от всех ошибок).

--- Сообщение объединено, 2 авг 2020 ---

Вот же: https://wasm.in/threads/portirovanie-fpu-na-sse.33827/page-3#post-420686

Вот ещё pRemDouble до кучи:

Код (Text):

1. function pRemDouble(D: Double; var X, Y: Double): integer;
2. var
3.   PiOf2H, PiOf2M, PiOf2L : Double;
4.   lowestPart: double;
5. const
6.   PiOf2Hi  : UInt64 = $3FF921FB54442D18; // Pi/2 : first 53bits
7.   PiOf2Mi  : UInt64 = $3C91A62633145C04; // Pi/2 : Middle 53bits
8.   PiOf2Lo  : UInt64 = $396707344A409382; // Pi/2 : Lowest 53bits
9.   ExpOffset1 =  54; // $3FF - $3C9
10.   ExpOffset2 = 105; // $3FF - $396
11.  
12.   procedure adddouble(const x: double; var rh, rl : double); inline;
13.   var
14.     hh, hl: double;
15.     temp : double;
16.   begin
17.     hh := x + rh;
18.     temp := hh - x;
19.     hl := (x - (hh - temp)) + (rh - temp);
20.     hl := hl + rl;
21.     rh := hh + hl;
22.     rl := hl - (rh - hh);
23.   end;
24.  
25.   procedure AddDouble2(const xh, xl: Double; var rh, rl : Double); inline;
26.   var
27.     hh, hl: Double;
28.     temp : double;
29.   begin
30.     hh := xh + rh;
31.     temp := hh - xh;
32.     hl := (xh - (hh - temp)) + (rh - temp);
33.     hl := hl + xl + rl;
34.     rh := hh + hl;
35.     rl := hl - (rh - hh);
36.   end;
37.  
38.   // 3*Double precision reduction
39.   procedure ReductionPi2n(n: integer);
40.   var
41.     X0,
42.     T0, T1, T2, T3: Double;
43.   begin
44.     PDoubleRec(@PiOf2H).Exp := n + $3FF;
45.     PDoubleRec(@PiOf2M).Exp := n + $3FF - ExpOffset1;
46.     PDoubleRec(@PiOf2L).Exp := n + $3FF - ExpOffset2;
47.  
48.     T0 := -PiOf2L;
49.     T1 := 0;
50.     AddDouble( lowestPart, T0, T1);
51.     AddDouble(-PiOf2M,     T0, T1);
52.     AddDouble( Y,          T0, T1); // T0&1 - partial sum
53.  
54.     T2 := T0;
55.     T3 := T1;
56.     AddDouble(-PiOf2H, T0, T1);
57.     AddDouble( X,      T0, T1); // T0&1 - X&Y
58.  
59.     X0 := X;
60.     X := T0;
61.     Y := T1;
62.  
63.     T0 := -T0;
64.     T1 := -T1;
65.     AddDouble( X0,     T0, T1);
66.     AddDouble(-PiOf2H, T0, T1);
67.     AddDouble2(T2,T3,  T0, T1);
68.  
69.     lowestPart := T0;
70.   end;
71.  
72. const
73.   nPi2High  : UInt64 = $BFF921FB54400000; // Pi/2 : Upper 32bit.
74.   nPi2Low   : UInt64 = $BDD0B4611A600000; // Pi/2 : Lower 32bit
75.   nPi2Small : UInt64 = $BBA3198A2E037073; // Pi - Pi2High - Pi2Low
76.   TwoOfPi  : UInt64 = $3FE45F306DC9C883; // 2 / Pi
77.   TwoPow22 : UInt64 = $4150000000000000; // 2 ^ 22
78.   TwoPow54 : UInt64 = $4350000000000000; // 2 ^ 54
79. var
80.   D2 : Double;
81.   Q: Int64;
82.   n: integer;
83.   hh, hl,
84.   pHi, pLow, pSmall,
85.   X2, Y2,
86.   temp : double;
87. begin
88.   Result := 0;
89.   X := D;
90.   Y := 0;
91.  
92.   D2 := Abs(D);
93.   if D2 <= Pi / 4 then Exit
94.   else if D2 <= 5 * Pi/4 then
95.   begin
96.     if D2 <= 3 * Pi/4 then Q := 1
97.     else Q := 2;
98.     if D < 0 then Q := -Q;
99.  
100.     pHi := Q * PDouble(@PiOf2Hi)^;
101.     pLow := Q * PDouble(@PiOf2Mi)^;
102.     pSmall := Q * PDouble(@PiOf2Lo)^;
103.  
104.     X2 := D - pHi;
105.     X := X2 - pLow;
106.     Y := - pLow - (X - X2);
107.     Y := Y - pSmall;
108.   end
109.   else if D2 < PDouble(@TwoPow22)^ then
110.   begin
111.     // Reduction for small/mid number (up to 2^22)
112.     Q := Trunc((D2 - Pi/4) * PDouble(@TwoOfPi)^) + 1;
113.  
114.     if D < 0 then Q := -Q;
115.     pHi := Q * PDouble(@nPi2High)^;
116.     pLow := Q * PDouble(@nPi2Low)^;
117.     pSmall := Q * PDouble(@nPi2Small)^;
118.  
119.     X2 := D + pHi;
120.  
121.     hh := pLow + X2;
122.     temp := hh - pLow;
123.     hl := (pLow - (hh - temp)) + (X2 - temp);
124.     X2 := hh + hl;
125.     Y2 := hl - (X2 - hh);
126.  
127.     hh := pSmall + X2;
128.     temp := hh - pSmall;
129.     hl := (pSmall - (hh - temp)) + (X2 - temp);
130.     hl := hl + Y2;
131.     X2 := hh + hl;
132.     Y2 := hl - (X2 - hh);
133.  
134.     X := X2;
135.     Y := Y2;
136.   end
137.   else
138.   begin
139.     // Reduction for large numbber between 2^22 to 2^53
140.     n := PDoubleRec(@D2).Exp - $3FF;
141.     PUInt64(@PiOf2H)^ := PiOf2Hi;
142.     PUInt64(@PiOf2M)^ := PiOf2Mi;
143.     PUInt64(@PiOf2L)^ := PiOf2Lo;
144.  
145.     X := D2;
146.     Y := 0;
147.     lowestPart := 0;
148.  
149.     while n >= 3 do
150.     begin
151.       PDoubleRec(@PiOf2H).Exp := n + $3FF;
152.       if PiOf2H < X then
153.         ReductionPi2n(n);
154.       Dec(n);
155.     end;
156.  
157.     Q := 0;
158.     while n >= 0 do
159.     begin
160.       PDoubleRec(@PiOf2H).Exp := n + $3FF;
161.       if PiOf2H < X then
162.       begin
163.         Q := Q or (1 shl n);
164.         ReductionPi2n(n);
165.       end;
166.       Dec(n);
167.     end;
168.  
169.     if X > Pi/4 then
170.     begin
171.       Q := Q + 1;
172.       ReductionPi2n(0);
173.     end;
174.     if D < 0 then
175.     begin
176.       X := -X;
177.       Y := -Y;
178.       Q := -Q;
179.     end;
180.   end;
181.   Result := Q and 3;
182. end;

--- Сообщение объединено, 2 авг 2020 ---

Ну не то чтоб ни о чём, но не сильно много.
Но я думаю, что если алгоритм усовершенствовать, можно добиться лучшего результата (до 2х, думаю, всё же можно довести). А если пожертвовать точностью (в довольно многих случаях большая точность не нужна), но ещё больше.

--- Сообщение объединено, 2 авг 2020 ---

А вот, кстати, и фокус-покус.
Написал на C++

Код (C++):

1. #include <iostream>
2. #include <windows.h>
3.  
4. volatile double y;
5.  
6. int main()
7. {
8.   double x = 0, delta = 1.74532925199432957692369e-7;
9.   DWORD c = GetTickCount();
10.   for (int i = 0; i < 36000000; ++i) {
11.     x = x + delta;
12.     y = std::sin(x);
13.   }
14.   c = GetTickCount() - c;
15.   std::cout << c;
16.   return 0;
17. }

и

Код (C++):

1. #include <iostream>
2. #include <windows.h>
3.  
4. volatile double y;
5.  
6. double finit()
7. {
8.   __asm {
9.         finit
10.   }
11. }
12.  
13. double fsin(double x)
14. {
15.   __asm {
16.         fld     qword ptr x
17.         fsin
18.         fstp    st
19.   }
20. }
21.  
22. int main()
23. {
24.   double x = 0, delta = 1.74532925199432957692369e-7;
25.   finit();
26.   DWORD c = GetTickCount();
27.   for (int i = 0; i < 36000000; ++i) {
28.     x = x + delta;
29.     y = fsin(x);
30.   }
31.   c = GetTickCount() - c;
32.   std::cout << c;
33.   return 0;
34. }

Компилю VC++
Первый код – 640 (x86) и 480 (x64).
Второй – порядка 1170 (x86, в x64 асм не поддерживается).
Так что, тут разница 1.8 - 2.4х.

--- Сообщение объединено, 2 авг 2020 ---

И ошибок тут 49438 из 36000000 (0,137 %).
Разница больше, чем на 1, опять же, в 417 случаях )
Больше, чем на 16 – в 43-х случаях. Т.е. очень похожая ситуация, только общее число ошибок меньше.

--- Сообщение объединено, 2 авг 2020 ---

Можно ещё под шлангом проверить, но уже в лом что-то.

 

Jin X,

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

> Но я думаю, что если алгоритм усовершенствовать, можно добиться лучшего результата

Сними профайл на последовательности стековых fpu вычислений, вероятно на порядок будет быстрее, чем полиномы. Тк не будет паразитных выборок в память, а числа длинные в памяти. Я не проверял, но уверен что тайминг резко пойдёт на спад с каждой +1 операцией fpu. И если что то вычислить на пол сотни операций fpu(сложную математику), то от этих полиномов не останется ничего. Примитивная операция +30% - это никто не заметит.

 

Не совсем по теме.
Переделка методов класса Fmatrix в движке XRay с помощью интринсиков.

Спойлер

Код (C++):

1.  
2. // Multiply RES = A[4x3]*B[4x3] (no projection), faster than ordinary multiply
3. ICF Fmatrix& mul_43 (const Fmatrix &A,const Fmatrix &B)
4. {
5. [indetn]// VERIFY ((this!=&A)&&(this!=&B));
6. /* m[0][0] = A.m[0][0] * B.m[0][0] + A.m[1][0] * B.m[0][1] + A.m[2][0] * B.m[0][2];
7. m[0][1] = A.m[0][1] * B.m[0][0] + A.m[1][1] * B.m[0][1] + A.m[2][1] * B.m[0][2];
8. m[0][2] = A.m[0][2] * B.m[0][0] + A.m[1][2] * B.m[0][1] + A.m[2][2] * B.m[0][2];
9. m[0][3] = 0;
10. m[1][0] = A.m[0][0] * B.m[1][0] + A.m[1][0] * B.m[1][1] + A.m[2][0] * B.m[1][2];
11. m[1][1] = A.m[0][1] * B.m[1][0] + A.m[1][1] * B.m[1][1] + A.m[2][1] * B.m[1][2];
12. m[1][2] = A.m[0][2] * B.m[1][0] + A.m[1][2] * B.m[1][1] + A.m[2][2] * B.m[1][2];
13. m[1][3] = 0;
14. m[2][0] = A.m[0][0] * B.m[2][0] + A.m[1][0] * B.m[2][1] + A.m[2][0] * B.m[2][2];
15. m[2][1] = A.m[0][1] * B.m[2][0] + A.m[1][1] * B.m[2][1] + A.m[2][1] * B.m[2][2];
16. m[2][2] = A.m[0][2] * B.m[2][0] + A.m[1][2] * B.m[2][1] + A.m[2][2] * B.m[2][2];
17. m[2][3] = 0;
18. m[3][0] = A.m[0][0] * B.m[3][0] + A.m[1][0] * B.m[3][1] + A.m[2][0] * B.m[3][2] + A.m[3][0];
19. m[3][1] = A.m[0][1] * B.m[3][0] + A.m[1][1] * B.m[3][1] + A.m[2][1] * B.m[3][2] + A.m[3][1];
20. m[3][2] = A.m[0][2] * B.m[3][0] + A.m[1][2] * B.m[3][1] + A.m[2][2] * B.m[3][2] + A.m[3][2];
21. m[3][3] = 1;*/
22. register __m128 Ai = A.i.m;
23. register __m128 Aj = A.j.m;
24. register __m128 Ak = A.k.m;
25. i.m = ADDPS(ADDPS(MULPS(SETPS1(B.i.x), Ai), MULPS(SETPS1(B.i.y), Aj)), MULPS(SETPS1(B.i.z), Ak));
26. j.m = ADDPS(ADDPS(MULPS(SETPS1(B.j.x), Ai), MULPS(SETPS1(B.j.y), Aj)), MULPS(SETPS1(B.j.z), Ak));
27. k.m = ADDPS(ADDPS(MULPS(SETPS1(B.k.x), Ai), MULPS(SETPS1(B.k.y), Aj)), MULPS(SETPS1(B.k.z), Ak));
28. c.m = ADDPS(ADDPS(ADDPS(MULPS(SETPS1(B.c.x), Ai), MULPS(SETPS1(B.c.y), Aj)), MULPS(SETPS1(B.c.z), Ak)), A.c.m);
29. VERIFY(c._==1.f);
30. return *this;[/indent]
31. }
32.  

Просто демонстрация оптимизации, код быстрей оригинала в 2.3 раза и где так же компактней, а значит быстрей будет загружаться.
PS
Дефайны!

Спойлер

Код (C++):

1.  
2. #define MINPS _mm_min_ps
3. #define MINSS _mm_min_ss
4. #define MAXPS _mm_max_ps
5. #define MAXSS _mm_max_ss
6. #define ADDPS _mm_add_ps
7. #define ADDSS _mm_add_ss
8. #define SUBPS _mm_sub_ps
9. #define SUBSS _mm_sub_ss
10. #define MULPS _mm_mul_ps
11. #define MULSS _mm_mul_ss
12. #define DIVPS _mm_div_ps
13. #define DIVSS _mm_div_ss
14. #define XORPS _mm_xor_ps
15. #define ANDPS _mm_and_ps
16. #define ANDNOTPS _mm_andnot_ps
17. #define ORPS _mm_or_ps
18. #define SETPS1 _mm_set_ps1
19. #define MOVSS _mm_set_ss
20. #define SHUFFLEPS _mm_shuffle_ps
21. #define SQRTPS _mm_sqrt_ps
22. #define SQRTSS _mm_sqrt_ss
23. #define SETPS3(Z) _mm_shuffle_ps(_mm_set_ss(Z),_mm_set_ss(Z),_MM_SHUFFLE(3,0,0,0))
24. #define ABSPS(A) _mm_andnot_ps(SIGNMASK3, A)
25.  

--- Сообщение объединено, 8 дек 2020 ---

А вот макро инструкция SINPS

Код (ASM):

1.  
2. SINPS MACRO X:req
3.     ;=x-x^3/3!+x^5/5!-x^7/7!
4.    movups xmm0, X
5.    movaps xmm2, xmm0
6.    movaps xmm3, xmm0
7.    ;1
8.    mulps xmm3, xmm3 ;=x*x
9.    mulps xmm2, xmm3 ;=x^3
10.    movaps xmm1, xmm2
11.    mulps xmm1, trigoFactorial3
12.    subps xmm0, xmm1
13.    ;2
14.    mulps xmm2, xmm3 ;=x^5
15.    movaps xmm1, xmm2
16.    mulps xmm1, trigoFactorial5
17.    addps xmm0, xmm1
18.    ;3
19.    mulps xmm2, xmm3 ;=x^7
20.    mulps xmm2, trigoFactorial7
21.    subps xmm0, xmm2
22. ENDM
23.  

Проверил, работает! На сколько быстро пока не проверял, думаю быстрей fsin.
ЗЫ
Константы такие, да я их секцию конст поместил, чтобы ида мне готовый код с интринсиками выдала!

Код (ASM):

1. .const
2. align 16
3. trigoFactorial3 real4 0.16666667, 0.16666667, 0.16666667, 0.16666667 ;=1/3!
4. trigoFactorial5 real4 0.0083333333, 0.0083333333, 0.0083333333, 0.0083333333 ;=1/5!
5. trigoFactorial7 real4 1.9841270e-4, 1.9841270e-4, 1.9841270e-4, 1.9841270e-4 ;=1/7!

 

Код (ASM):

1. SINPS MACRO X:req
2.   ;=x-x^3/3!+x^5/5!-x^7/7!
3.   movups xmm0, X  ; и здесь бы выровнять
4.   movaps xmm2, xmm0
5.   movaps xmm3, xmm0
6.   mulps xmm3, xmm3 ;=x*x
7.   ;1
8.   mulps xmm2, xmm3 ;=x^3
9.   mulps xmm2, 1div6 ;=x^3/2/3
10.   subps xmm0, xmm2
11.   ;2
12.   mulps xmm2, xmm3 ;=x^5/2/3
13.   mulps xmm2, 1div20 ;=x^5/2/3/4/5
14.   addps xmm0, xmm2
15.   ;3
16.   mulps xmm2, xmm3 ;=x^7/2/3/4/5
17.   mulps xmm2, 1div42 ;=x^7/2/3/4/5/6/7
18.   subps xmm0, xmm2
19. ENDM

 

НетРегистрации, прикольно. Я тоже думал что с константами возможна дополнительная оптимизация, но не смог математически доказать правильность.
Да, оказалось надо больше итераций, а так же ограничить х от -PI до PI, на что нужен дополнительный код, так что тут не всё просто. А с другой стороны, что нам запрещает предпринять специальные методы, чтобы угол ограничить -180..180, а если что пошло не так, сработает ассерт.

 

Intro, а как там с точностью результата?
Можно ещё сделать макросы SINPS_HP для более точного расчёта (взять ещё 4 члена полинома, до 1/15!). А также SINPD, SINPD_HP.

--- Сообщение объединено, 22 дек 2020 ---

По скорости SINPS у меня получился в 25 раз быстрее

--- Сообщение объединено, 22 дек 2020 ---

(чем fsin)

--- Сообщение объединено, 22 дек 2020 ---

Вот такой код работает у меня в 1.7x быстрее (я его решил назвать LP, т.к. медленнее):

Код (ASM):

1. ; fasm
2. MACRO SINPS_LP X*
3. {
4.   ;=x-x^3/3!+x^5/5!-x^7/7!
5.   movups xmm0, X
6.   movaps xmm1, xword [_1div6] ;=1/2/3
7.   movaps xmm2, xmm0 ;=x
8.   movaps xmm3, xmm0
9.   mulps xmm3, xmm3 ;=x^2
10.   movaps xmm4, xmm3
11.   mulps xmm4, xmm4 ;=x^4
12.   ;1
13.   mulps xmm3, xmm0 ;=x^3
14.   mulps xmm1, xmm3 ;=x^3/2/3
15.   ;2
16.   mulps xmm2, xmm4 ;=x^5
17.   mulps xmm2, xword [_1div120] ;=x^5/2/3/4/5
18.   addps xmm0, xmm2 ;x+x^5/5!
19.   ;3
20.   mulps xmm3, xmm4 ;=x^7
21.   mulps xmm3, xword [_1div5040] ;=x^7/2/3/4/5/6/7
22.   addps xmm1, xmm3 ;=x^3/3!+x^7/7!
23.  
24.   subps xmm0, xmm1
25. }
26. ...
27. _1div6 dd 4 dup (0.1666666666667)
28. _1div120 dd 4 dup (8.333333333333e-3)
29. _1div5040 dd 4 dup (1.984126984127e-4)
30.  

И более точные вычисления (медленнее в 2.1x, чем предыдущий SINPS_LP и в 1.3x медленнее, чем код от НетРегистрации):

Код (Text):

1. ; fasm
2. MACRO SINPS X*
3. {
4.   ;=x-x^3/3!+x^5/5!-x^7/7!+x^9/9!-x^11/11!+x^13/13!-x^15/15!
5.   movups xmm0, xword [X]  ; и здесь бы выровнять
6.   movaps xmm1, xword [_1div6] ;=1/2/3
7.   movaps xmm2, xmm0 ;=x
8.   movaps xmm3, xmm0
9.   mulps xmm3, xmm3 ;=x^2
10.   movaps xmm4, xmm3
11.   mulps xmm4, xmm4 ;=x^4
12.   ;1
13.   mulps xmm3, xmm0 ;=x^3
14.   mulps xmm1, xmm3 ;=x^3/2/3
15.   ;2
16.   mulps xmm2, xmm4 ;=x^5
17.   mulps xmm2, xword [_1div120] ;=x^5/2/3/4/5
18.   addps xmm0, xmm2 ;x+x^5/5!
19.   ;3
20.   mulps xmm3, xmm4 ;=x^7
21.   mulps xmm3, xword [_1div5040] ;=x^7/2/3/4/5/6/7
22.   addps xmm1, xmm3 ;=x^3/3!+x^7/7!
23.   ;4
24.   mulps xmm2, xmm4 ;=x^9/2/3/4/5
25.   mulps xmm2, xword [_1div3024] ;=x^9/2/3/4/5/6/7/8/9
26.   addps xmm0, xmm2 ;x+x^5/5!+x^9/9!
27.   ;5
28.   mulps xmm3, xmm4 ;=x^11/2/3/4/5/6/7
29.   mulps xmm3, xword [_1div7920] ;=x^11/2/3/4/5/6/7/8/9/10/11
30.   addps xmm1, xmm3 ;=x^3/3!+x^7/7!+x^11/11!
31.   ;6
32.   mulps xmm2, xmm4 ;=x^13/2/3/4/5/6/7/8/9
33.   mulps xmm2, xword [_1div17160] ;=x^13/2/3/4/5/6/7/8/9/10/11/12/13
34.   addps xmm0, xmm2 ;x+x^5/5!+x^9/9!+x^13/13!
35.   ;7
36.   mulps xmm3, xmm4 ;=x^11/2/3/4/5/6/7/8/9/10/11
37.   mulps xmm3, xword [_1div32760] ;=x^11/2/3/4/5/6/7/8/9/10/11/12/13/14/15
38.   addps xmm1, xmm3 ;=x^3/3!+x^7/7!+x^11/11!+x^15/15!
39.  
40.   subps xmm0, xmm1
41. }
42. ...
43. _1div6 dd 4 dup (0.1666666666667)
44. _1div120 dd 4 dup (8.333333333333e-3)
45. _1div5040 dd 4 dup (1.984126984127e-4)
46. _1div3024 dd 4 dup (3.306878306878e-4)
47. _1div7920 dd 4 dup (1.262626262626e-4)
48. _1div17160 dd 4 dup (5.827505827506e-5)
49. _1div32760 dd 4 dup (3.052503052503e-5)
50.  

Все вычисления верные, проверено.
Но погрешности всё равно есть.
Можно подумать ещё как сделать вычисления более точными и более быстрыми.
К примеру, можно свернуть формулу вот так:

--- Сообщение объединено, 22 дек 2020 ---

Блин, отправил случайно...
Ладно, без примера, уже плохо соображаю в 3 часа ночи.

Вот тестовый код.

 

Обычная точность ~3байта, внутренняя SSE мне неизвестна, 1/15! ну вы поняли, что может быть при суммировании.
Соответственно можно лишь попытаться улучшить точность, начиная суммирование с конца.

 

Можно ещё через умножение сделать.
Т.е. не складывать x^5/5! + x^9/9! + x^13/13!, а (((x^4/13!+1/9!)*x^4+1/5!)*x^5. Как, собственно, это в Delphi примерно и делается.

По поводу 25x – это я наврал, неправильно затестил.
fsin медленнее в сравнении с моим sinps в 2.9x.
Вот тестовый код вместе с ним.

Спойлер: И вывод программы

Код (Text):

1. Success: both RDTSCP instruction and invariant TSC are supported.
2. Testing SINPS (unoptimized)...... 60 ticks
3. Testing SINPS.................... 36 ticks
4. Testing SINPS_LP (unoptimized)... 27 ticks
5. Testing SINPS_LP................. 18 ticks
6. Testing FSIN..................... 105 ticks
7. Press a key to exit...

 

О как!
Вынес fld1 в начало кода (после finit), скорость fsin увеличилась до 66 тактов.
И вообще, все цифры чуть уменьшились... почему-то

Спойлер: Вывод

Код (Text):

1. Testing SINPS (unoptimized)...... 55 ticks
2. Testing SINPS.................... 31 ticks
3. Testing SINPS_LP (unoptimized)... 25 ticks
4. Testing SINPS_LP................. 13 ticks
5. Testing FSIN..................... 67 ticks

Хотя, если оставить fld1 + fsin + fstp st0, либо добавить ещё одну пару fld1 + fstp st0, разницы почти не будет.
Интересно, почему так?

 

Jin X,

Так потому что это ошибка профайлера(системных таймеров), у него разрешение на много порядков ниже, а есчо он периодически вытисняется как задача. Замер должен быть в цикле на очень большом числе итераций. Иначе будет рандом.

--- Сообщение объединено, 22 дек 2020 ---

> Success: both RDTSCP instruction and invariant TSC are supported.

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

 

Так, там 256 циклов + 16 на разогрев. И из 256 половина отбрасывается (64 самых быстрых и 64 самых медленных). Минус холостой, разумеется. И я этот тест запускал много раз. Показывает одинаковые значения с разбегом 1-2%.
А сейчас сделал 65536 (с разогревом 256), разбег стал ещё меньше (чаще всего числа всегда одинаковые). Но проблема осталась: когда fld/ftsp внутри функции, все тесты показывают чуть большее значение, а когда снаружи, все показывают чуть меньше. Хотя, они только в одном тесте присутствуют. Причём, все функции выравнены по 16 байтам, а этот тест последний, т.е. тестовые процедуры сдвинуться не могут.

--- Сообщение объединено, 23 дек 2020 ---

А вот сейчас я ещё и сами процедуры замера выравнил по 16 байтам.
И результаты стали одинаковыми, кроме последнего теста. Как и должно быть.

--- Сообщение объединено, 23 дек 2020 ---

Вот.

--- Сообщение объединено, 23 дек 2020 ---

Почему плавает?
У меня разница 4 такта (хотя понимаю, что может быть и больше).
Но учитывая, что он ожидает завершение предыдущих операций (не надо вызывать cpuid, который действительно медленный), то его использование наоборот быстрее получается.

Вообще, хотелось бы, чтобы спецы проинспектировали функции чтения счётчиков и поправили, если что не так:

Код (ASM):

1.  
2. ; Read TSC via RDTSC [for internal use]
3. ; Returns: rax = current TSC counter value, ecx = 0 (processor id detection is not supported)
4. ; Changes ebx !!!
5. if used SpeedTestRDTSC
6. SpeedTestRDTSC:
7.         xor    eax,eax        ; cpuid execution time may vary depending on eax value
8.         cpuid            ; serialization
9.         xor    ecx,ecx        ; processor id (not supported)
10.         rdtsc
11.         shl    rdx,32
12.         or    rax,rdx
13.         mfence
14.         ret
15. end if ; used SpeedTestRDTSC
16.  
17. ; Read TSC via RDTSCP [for internal use]
18. ; Returns: rax = current TSC counter value, ecx = processor id
19. if used SpeedTestRDTSCP
20. SpeedTestRDTSCP:
21.         mfence
22.         rdtscp
23.         shl    rdx,32
24.         or    rax,rdx
25.         mfence
26. SpeedTestEmptyFunc:
27.         ret
28. end if ; used SpeedTestRDTSCP

--- Сообщение объединено, 23 дек 2020 ---

Блин, mfence перед rdtscp тут не должно быть, случайно написал :-/