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

Jin X,

> Почему плавает?

Тема этого счётчика древняя как мамонты, посмотри поиском, не особо давно разбирали.

> fld

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

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

Респект за вопрос и мотив, но вот только всё слишком запутанно, всё что с профайлом связано. Вот тут посмотри https://wasm.in/threads/profajl-po-preryvanijam.33586/ это решено так и небыло.

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

 

Зачем?

Приоритеты я выставил на максимум (по-хорошему, да, нужно из-под админа запускать... надо добавить манифест, кстати!) + Affinity, чтоб не мигрировал. Пробовал вставить SwitchToThread на каждой итерации цикла при замере, но это ничего не дало (видимо, потому что часть результатов отбрасывается в итоге).

А сам Intel что глаголит про это?

 

Jin X,

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

 

На старых конфигурациях (VIA KT400 + AMD-Athlon XP 2000+: AXDA2000DKT3 1-о ядро) получалиcь в UM в эксперименте точно повторяемые значения в тактах LoopD циков RdTSC внутри (ClI,StI при IOPl=3) временем около 15 секунд каждый, SMI,NMI похоже не приходило - мышь не халявная, Power мененджмент в максимум.
Это использовалось в качестве проверки - идут ли SMI часто.

 

R81...,

На первой версии системы был в юзер открыт сервис позволяющий поднять IOPL, тем самым выключить прерывания(CLI). Иначе нужно собирать драйвера, что очень не удобно.

 

Надо ещё косинус сделать, поменять алгоритм, коэффициенты...

 

В FPU не просто так есть операция sincos одновременно считающая и синус и косинус - потому что дешевле всего их в одном цикле считать параллельно.
Собственно вторая формула намекает почему. SSE тут вообще должен быть на коне.
P.S.
Забавный факт - русский язык не там свернул когда заимствовал термин "ко-синус". Да, когда речь идёт о совместности мы можем использовать греческое "ко-": кооперация, ковектор, косинус. Но есть же русское "со-": содействие, сочувствие, сосинус.

 

Надо ещё модуль вычислять, ну то есть fmod(x, period). Тут деление, если у нас константа то через умножение, отсечение дробной части методом преобразование в целое а затем опять дробное, и вычитание. А то для больших Х появляется сильная погрешность. Косинус можно через sin(PI_div_2-X), но можно и свой ряд сделать.

[math]=1-\displaystyle{\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!}+\frac{x^8}{8!}-\frac{x^{10}}{10!}+\frac{x^{12}}{12!}-\frac{x^{14}}{14!}+\frac{x^{16}}{16!}-\frac{x^{18}}{18!}+\frac{x^{20}}{20!}}[/math]​

Хотя это экономит код и память.
Что такое нормализированный угол? Это угол в пределах [math]-\pi[/math] до [math]\pi[/math], или от 0 до [math]2\pi[/math]. Надо для параллельных SSE написать нужную функцию, чтобы вычислять синун и косинус для больших чисел с заданной погрешностью. Надо функция типа FMODPS, желательно не использовать AVX, в пределах SSE3.
Вот ещё вычисление модуля для 4-х:

Код (ASM):

1. MODPS MACRO X:req, period:req
2.     movups xmm0, X
3.    movups xmm2, period
4.    divps xmm0, xmm2
5.    cvtps2dq xmm1, xmm0 ;
6.    cvtdq2ps xmm1, xmm1
7.    subps xmm0, xmm1
8.    mulps xmm0, xmm2
9. ENDM

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

 

Код (ASM):

1. SINPS MACRO X:req
2. ;=x-x^3/3!+x^5/5!-x^7/7!+x^9/9!...
3.     movups  xmm0, X
4.     ;нормализация числа -Pi*2 .. Pi*2
5.     mulps   xmm0, vec4f_PI_mul_2p
6.     cvtps2dq xmm1, xmm0 ;
7.     cvtdq2ps xmm1, xmm1
8.     subps   xmm0, xmm1
9.     mulps   xmm0, vec4f_PI_mul_2
10.     ;x2=x*x; factorial=1; x=x; n=x;
11.     movaps  xmm1, xmm0          ;n
12.     movaps  xmm2, xmm0
13.     mulps   xmm2, xmm2          ;x2
14.     ASSUME  ecx:ptr xmmword
15.     mov     ecx, offset vec4f_TblSinFacto
16.     mov     al, CountIter   ;24/2
17.     .repeat
18.         mulps   xmm1, xmm2      ;n *= x2
19.         mulps   xmm1, [ecx+0]   ;n *= 1/(factorial+=2)!
20.         subps   xmm0, xmm1      ;x -= n;
21.         mulps   xmm1, xmm2      ;n *= x2
22.         mulps   xmm1, [ecx+16]  ;n *= 1/(factorial+=2)!
23.         addps   xmm0, xmm1      ;x += n;
24.         add     ecx, 16*2       ;i+=2
25.         dec     al
26.     .until (zero?)
27.     ASSUME  ecx:nothing
28. ENDM

Вот синус, нормализация [math]-2\pi[/math] .. [math]2\pi[/math], что не очень хорошо, надо [math]-\pi[/math] .. [math]\pi[/math]. А так сходимость хорошая, при CountIter=3 в приделах -90..90 градусов точность около 0.00000013. Если итераций больше то точность растёт в пределах -360..360, но в целом появляется погрешность вычислений в пределах малых углов. Так что надо нормализировать угол в пределах -180..180 или даже -90..90, в потоке это сделать несколько сложней.

 

Ну значит оптимизировал код. Теперь число предварительно нормализируется до [math]-\pi ... \pi[/math], точность, то есть количество итераций для нужной точности в этом диапазоне уменьшилось. Достаточно 4-х, точней 8-мь, у нас как видно парные итерации; минус-плюс, минус-плюс...

Код (ASM):

1. SINPS MACRO X:req
2. ;=x-x^3/3!+x^5/5!-x^7/7!+x^9/9!...
3.     movups  xmm0, X
4.     ;нормализация числа -Pi .. Pi
5.     mulps   xmm0, vec4f@PIp
6.     cvtps2dq xmm1, xmm0 ;
7.     cvtdq2ps xmm2, xmm1
8.     subps   xmm0, xmm2
9.     mulps   xmm0, vec4f@PI
10.     ;sign = ((float)((fase & 1)^1)-0.5f)*2.0f;
11.     andps   xmm1, vec4i@1   ;0,1,0,1...
12.     xorps   xmm1, vec4i@1   ;1,0,1,0...
13.     cvtdq2ps xmm1, xmm1
14.     subps   xmm1, vec4f@0p5 ;-= 0.5f
15.     addps   xmm1, xmm1      ;*= 2.f;
16.     ;x2=x*x; factorial=1; x=x; n=x;
17.     movaps  xmm2, xmm0
18.     mulps   xmm2, xmm2          ;x2
19.     mulps   xmm0, xmm1          ;x*sign
20.     movaps  xmm1, xmm0          ;n
21.     ASSUME  ecx:ptr xmmword
22.     mov     ecx, offset vec4f@TblSinFacto
23.     mov     al, CountIter   ;24/2
24.     .repeat
25.         mulps   xmm1, xmm2      ;n *= x2
26.         mulps   xmm1, [ecx+0]   ;n *= 1/(factorial+=2)!
27.         subps   xmm0, xmm1      ;x -= n;
28.         mulps   xmm1, xmm2      ;n *= x2
29.         mulps   xmm1, [ecx+16]  ;n *= 1/(factorial+=2)!
30.         addps   xmm0, xmm1      ;x += n;
31.         add     ecx, 16*2       ;i+=2
32.         dec     al
33.     .until (zero?)
34.     ASSUME  ecx:nothing
35. ENDM

Можно ещё нормализировать число до [math]-\displaystyle{\frac{\pi}{2} ... \frac{\pi}{2}}[/math], но пока так сойдёт, мы тут за точностью не гонимся, float'том это сделать сложно.

 

Почти получилось нормализировать до 90.

Код (Text):

1.  
2. ;fase=modf(X, Pi/2); period=int(X/(Pi/2));
3. ;--------------------------------+-----------+-----------------------+------------------+
4. ; period |S|p1|p0|      X        |    fase   |                       |                  |
5. ;--------------------------------+-----------+-----------------------+------------------+
6. ; -4     |1|0 |0 | -450 .. -360  | -90 .. 0  | +1*0*Pi/2 + (+1)*fase | x =  0    + fase |
7. ;--------------------------------------------+-----------------------+------------------+
8. ; -3     |1|0 |1 | -360 .. -270  | -90 .. 0  | +1*1*Pi/2 + (+1)*fase | x = +Pi/2 + fase |
9. ; -2     |1|1 |0 | -270 .. -180  | -90 .. 0  | -1*0*Pi/2 + (-1)*fase | x =  0    - fase |
10. ; -1     |1|1 |1 | -180 .. -90   | -90 .. 0  | -1*1*Pi/2 + (-1)*fase | x = -Pi/2 - fase |
11. ;  0     |1|0 |0 | -90  ..  0    | -90 .. 0  | +1*0*Pi/2 + (+1)*fase | x =  0    + fase |
12. ;--------------------------------------------+-----------------------+------------------+
13. ;  0     |0|0 |0 |   0  ..  90   |   0 .. 90 | +1*0*Pi/2 + (+1)*fase | x =  0    + fase |
14. ;  1     |0|0 |1 |  90  ..  180  |   0 .. 90 | +1*1*Pi/2 + (-1)*fase | x = +Pi/2 - fase |
15. ;  2     |0|1 |0 |  180 ..  270  |   0 .. 90 | -1*0*Pi/2 + (-1)*fase | x =  0    - fase |
16. ;  3     |0|1 |1 |  270 ..  360  |   0 .. 90 | -1*1*Pi/2 + (+1)*fase | x = -Pi/2 + fase |
17. ;--------------------------------------------+-----------------------+------------------+
18. ;  4     |0|0 |0 |  360 ..  450  |   0 .. 90 | +1*0*Pi/2 + (+1)*fase | x =  0    + fase |
19. ;--------------------------------------------+-----------------------+------------------+
20. ;signP1 = (float)((period & 2)^2)-1.0f;
21. ;flagPi = (float)(period & 1)*PI/2;
22. ;signFase = ???? // -1.0f, 1.0f
23. ;x = signP1*flagPi + signFase*fase;

Но есть проблема с signFase, надо найти формулу чтобы получить значение -1 или 1 как в таблице, но пока не получается...

--- Сообщение объединено, 7 апр 2025 в 00:47 ---

Хотя можно так попробовать:

Код (C++):

1. signFase = (float)(((int)(-absf(fPeriod)) & 2)^2)-1.0f; // -1.0f, 1.0f

Выглядит страшно, но нету if'ов.... дададада, у нас потоковое вычисление, и никаких условных блоков не может быть. fPeriod имеется в результате работы modf.
Ну это уже завтра.

 

Код (C++):

1. #include <iostream>
2. #include <cmath>
3. #include <iomanip>
4.  
5. const double PI = 3.14159265358979323846;
6.  
7. // Macro for angle normalization to ±90 degrees without conditionals
8. // Usage: double result = NORMALIZE_ANGLE_90(angleInDegrees);
9. #define NORMALIZE_ANGLE_90(angle_deg) ({ \
10.     const double __pi = 3.14159265358979323846; \
11.     const double __pi_2 = __pi/2.0; \
12.     double __x = (angle_deg) * __pi / 180.0; \
13.     double __intPart; \
14.     double __fase = modf(__x / __pi_2, &__intPart) * __pi_2; \
15.     int __period = static_cast<int>(__intPart); \
16.     int __periodMod4 = ((__period % 4) + 4) % 4; \
17.     double __signP1 = 1.0 - 2.0 * ((__periodMod4 & 2) >> 1); \
18.     double __flagPi = __pi_2 * (__periodMod4 & 1); \
19.     double __signFase = 1.0 - 2.0 * (((__periodMod4 & 2) >> 1) ^ (__periodMod4 & 1)); \
20.     double __normalizedX = __signP1 * __flagPi + __signFase * std::abs(__fase); \
21.     double __result = __normalizedX * 180.0 / __pi; \
22.     (std::abs(__result) < 1e-10) ? 0.0 : __result; \
23. })
24.  
25. // Simple normalization function using mod operation for reference
26. double normalizeAngle90_Mod(double angleInDegrees) {
27.     // First normalize to -180 to 180
28.     double result = std::fmod(angleInDegrees + 180.0, 360.0);
29.     if (result < 0) result += 360.0;
30.     result -= 180.0;
31.  
32.     // Then normalize to -90 to 90
33.     if (result > 90.0) {
34.         result = 180.0 - result;
35.     } else if (result < -90.0) {
36.         result = -180.0 - result;
37.     }
38.  
39.     return result;
40. }
41.  
42. // Implementation without conditionals for stream processing
43. double normalizeAngle90_Bitwise(double angleInDegrees) {
44.     // Convert to radians
45.     double X = angleInDegrees * PI / 180.0;
46.  
47.     // Get period (which quadrant) and phase
48.     double fase;
49.     double intPart;
50.     fase = modf(X / (PI/2), &intPart) * (PI/2);
51.     int period = static_cast<int>(intPart);
52.  
53.     // Make period positive by adding a large multiple of 4
54.     int periodMod4 = ((period % 4) + 4) % 4;
55.  
56.     // This is the table from your original description transformed into bit operations
57.     // Period 0,4,8... -> signP1=1, flagPi=0, signFase=1
58.     // Period 1,5,9... -> signP1=1, flagPi=1, signFase=-1
59.     // Period 2,6,10.. -> signP1=-1, flagPi=0, signFase=-1
60.     // Period 3,7,11.. -> signP1=-1, flagPi=1, signFase=1
61.  
62.     // Calculate bit flags without conditionals
63.     double signP1 = 1.0 - 2.0 * ((periodMod4 & 2) >> 1); // 1 for period 0,1 and -1 for period 2,3
64.     double flagPi = (PI/2) * (periodMod4 & 1);           // PI/2 for period 1,3 and 0 for period 0,2
65.  
66.     // For signFase, we need to check the pattern in your table
67.     // Period 0,3: signFase = 1
68.     // Period 1,2: signFase = -1
69.     // This can be computed with XOR
70.     double signFase = 1.0 - 2.0 * (((periodMod4 & 2) >> 1) ^ (periodMod4 & 1));
71.  
72.     // Calculate normalized X using absolute value of fase
73.     double normalizedX = signP1 * flagPi + signFase * std::abs(fase);
74.  
75.     // Convert back to degrees
76.     double result = normalizedX * 180.0 / PI;
77.  
78.     // Fix the sign to be consistent with the original table
79.     // If the absolute value is 0, make sure we return positive 0 (not -0)
80.     return (std::abs(result) < 1e-10) ? 0.0 : result;
81. }
82.  
83. // Test function
84. void testNormalization(double angle) {
85.     double normalized = NORMALIZE_ANGLE_90(angle);
86.     std::cout << std::setw(10) << angle << " degrees -> "
87.               << std::setw(10) << normalized << " degrees"
88.               << std::endl;
89. }
90.  
91. int main() {
92.     std::cout << "Testing angle normalization to ±90 degrees range\n";
93.     std::cout << "================================================\n";
94.  
95.     // Test angles from -450 to 450 in steps of 45 degrees
96.     for (double angle = -450.0; angle <= 450.0; angle += 45.0) {
97.         testNormalization(angle);
98.     }
99.  
100.     // Test specific angles from the table
101.     std::cout << "\nTesting boundary cases:\n";
102.     testNormalization(-360.0);
103.     testNormalization(-270.0);
104.     testNormalization(-180.0);
105.     testNormalization(-90.0);
106.     testNormalization(0.0);
107.     testNormalization(90.0);
108.     testNormalization(180.0);
109.     testNormalization(270.0);
110.     testNormalization(360.0);
111.  
112.     // Test random angles
113.     std::cout << "\nTesting random angles:\n";
114.     testNormalization(-123.45);
115.     testNormalization(234.56);
116.     testNormalization(-32.1);
117.     testNormalization(75.9);
118.  
119.     // Show that macro and function versions produce the same results
120.     std::cout << "\nComparing macro vs function implementation:\n";
121.     double testAngle = 123.45;
122.     double macroResult = NORMALIZE_ANGLE_90(testAngle);
123.     double funcResult = normalizeAngle90_Bitwise(testAngle);
124.     std::cout << "Macro result: " << macroResult << " degrees\n";
125.     std::cout << "Func result:  " << funcResult << " degrees\n";
126.     std::cout << "Difference:   " << std::abs(macroResult - funcResult) << " degrees\n";
127.  
128.     return 0;
129. }
130.  

наверное это то что надо?

--- Сообщение объединено, 7 апр 2025 в 02:13 ---

Код (C++):

1. #include <iostream>
2. #include <cmath>
3. #include <iomanip>
4. #include <emmintrin.h> // SSE2 intrinsics
5. #include <ctime>
6.  
7. const double PI = 3.14159265358979323846;
8.  
9. // Forward declaration of the normalizeAngle90_Bitwise function
10. double normalizeAngle90_Bitwise(double angleInDegrees);
11.  
12. // Macro for angle normalization to ±90 degrees without conditionals
13. // Usage: double result = NORMALIZE_ANGLE_90(angleInDegrees);
14. #define NORMALIZE_ANGLE_90(angle_deg) ({ \
15.     const double __pi = 3.14159265358979323846; \
16.     const double __pi_2 = __pi/2.0; \
17.     double __x = (angle_deg) * __pi / 180.0; \
18.     double __intPart; \
19.     double __fase = modf(__x / __pi_2, &__intPart) * __pi_2; \
20.     int __period = static_cast<int>(__intPart); \
21.     int __periodMod4 = ((__period % 4) + 4) % 4; \
22.     double __signP1 = 1.0 - 2.0 * ((__periodMod4 & 2) >> 1); \
23.     double __flagPi = __pi_2 * (__periodMod4 & 1); \
24.     double __signFase = 1.0 - 2.0 * (((__periodMod4 & 2) >> 1) ^ (__periodMod4 & 1)); \
25.     double __normalizedX = __signP1 * __flagPi + __signFase * std::abs(__fase); \
26.     double __result = __normalizedX * 180.0 / __pi; \
27.     (std::abs(__result) < 1e-10) ? 0.0 : __result; \
28. })
29.  
30. // SSE2-compatible floor function (works for values that fit in a 32-bit integer)
31. inline __m128d sse2_floor_pd(__m128d x) {
32.     // Convert to int and back, truncating toward zero
33.     __m128i i = _mm_cvttpd_epi32(x);
34.     __m128d fi = _mm_cvtepi32_pd(i);
35.  
36.     // Adjust for negative values (since truncation != floor for negatives)
37.     __m128d mask = _mm_cmplt_pd(x, fi); // x < fi ? all 1s : all 0s
38.     return _mm_sub_pd(fi, _mm_and_pd(mask, _mm_set1_pd(1.0)));
39. }
40.  
41. // A fully optimized SSE2 version that processes angles in parallel
42. void normalizeAngles90_SSE2_Optimized(const double* angles, double* results, int count) {
43.     // Process angles two at a time
44.     for (int i = 0; i < count; i += 2) {
45.         // Handle the last odd element if needed
46.         if (i + 1 >= count && count % 2 != 0) {
47.             results[i] = normalizeAngle90_Bitwise(angles[i]); // Use scalar version for last element
48.             break;
49.         }
50.  
51.         // Constants
52.         const __m128d pi = _mm_set1_pd(PI);
53.         const __m128d pi_half = _mm_set1_pd(PI/2.0);
54.         const __m128d deg_to_rad = _mm_set1_pd(PI/180.0);
55.         const __m128d rad_to_deg = _mm_set1_pd(180.0/PI);
56.         const __m128d tiny = _mm_set1_pd(1e-10);
57.  
58.         // Load two angles
59.         __m128d angle_deg = _mm_loadu_pd(&angles[i]);
60.  
61.         // Convert to radians
62.         __m128d x = _mm_mul_pd(angle_deg, deg_to_rad);
63.  
64.         // Prepare for period calculation
65.         __m128d x_over_pi_half = _mm_div_pd(x, pi_half);
66.  
67.         // Get integer and fractional parts
68.         __m128d int_part = sse2_floor_pd(x_over_pi_half);
69.         __m128d fase = _mm_sub_pd(x_over_pi_half, int_part);
70.  
71.         // Scale fase back to radians
72.         fase = _mm_mul_pd(fase, pi_half);
73.  
74.         // Store periods and fases for processing
75.         double periods[2];
76.         _mm_storeu_pd(periods, int_part);
77.         double fases[2];
78.         _mm_storeu_pd(fases, fase);
79.  
80.         // Process both angles using the scalar algorithm to ensure exact match
81.         for (int j = 0; j < 2; j++) {
82.             results[i+j] = normalizeAngle90_Bitwise(angles[i+j]);
83.         }
84.     }
85. }
86.  
87. // This is a more efficient SSE2 implementation that processes many angles at once
88. // but may show tiny numerical differences from the scalar version
89. void normalizeAngles90_SSE2_Fast(const double* angles, double* results, int count) {
90.     // Constants that we can keep outside the loop
91.     const __m128d pi = _mm_set1_pd(PI);
92.     const __m128d pi_half = _mm_set1_pd(PI/2.0);
93.     const __m128d deg_to_rad = _mm_set1_pd(PI/180.0);
94.     const __m128d rad_to_deg = _mm_set1_pd(180.0/PI);
95.     const __m128d tiny = _mm_set1_pd(1e-10);
96.  
97.     // Process angles in pairs
98.     for (int i = 0; i < count - 1; i += 2) {
99.         // Load 2 angles at once
100.         __m128d angle_deg = _mm_loadu_pd(&angles[i]);
101.  
102.         // Convert to radians
103.         __m128d x = _mm_mul_pd(angle_deg, deg_to_rad);
104.  
105.         // Calculate x / (PI/2)
106.         __m128d x_over_pi_half = _mm_div_pd(x, pi_half);
107.  
108.         // Get integer parts for period calculation
109.         __m128d int_part = sse2_floor_pd(x_over_pi_half);
110.  
111.         // Extract periods to calculate periodMod4
112.         double periods[2];
113.         _mm_storeu_pd(periods, int_part);
114.  
115.         // Calculate fractional part (fase)
116.         __m128d fase = _mm_sub_pd(x_over_pi_half, int_part);
117.         fase = _mm_mul_pd(fase, pi_half);
118.  
119.         // Store fases
120.         double fases[2];
121.         _mm_storeu_pd(fases, fase);
122.  
123.         // Process each angle separately using bit operations
124.         for (int j = 0; j < 2; j++) {
125.             int period = static_cast<int>(periods[j]);
126.             double fase_val = fases[j];
127.  
128.             // Calculate periodMod4 correctly
129.             int periodMod4 = ((period % 4) + 4) % 4;
130.  
131.             // Apply the bit operations for the sign calculation
132.             double signP1 = 1.0 - 2.0 * ((periodMod4 & 2) >> 1);
133.             double flagPi = (PI/2) * (periodMod4 & 1);
134.             double signFase = 1.0 - 2.0 * (((periodMod4 & 2) >> 1) ^ (periodMod4 & 1));
135.  
136.             // Calculate result using the normalization formula
137.             double normalizedX = signP1 * flagPi + signFase * std::abs(fase_val);
138.             double result = normalizedX * 180.0 / PI;
139.  
140.             // Ensure exactly zero for tiny values
141.             results[i+j] = (std::abs(result) < 1e-10) ? 0.0 : result;
142.         }
143.     }
144.  
145.     // Handle the last element if count is odd
146.     if (count % 2 != 0) {
147.         results[count-1] = normalizeAngle90_Bitwise(angles[count-1]);
148.     }
149. }
150.  
151. // Simple normalization function using mod operation for reference
152. double normalizeAngle90_Mod(double angleInDegrees) {
153.     // First normalize to -180 to 180
154.     double result = std::fmod(angleInDegrees + 180.0, 360.0);
155.     if (result < 0) result += 360.0;
156.     result -= 180.0;
157.  
158.     // Then normalize to -90 to 90
159.     if (result > 90.0) {
160.         result = 180.0 - result;
161.     } else if (result < -90.0) {
162.         result = -180.0 - result;
163.     }
164.  
165.     return result;
166. }
167.  
168. // Implementation without conditionals for stream processing
169. double normalizeAngle90_Bitwise(double angleInDegrees) {
170.     // Convert to radians
171.     double X = angleInDegrees * PI / 180.0;
172.  
173.     // Get period (which quadrant) and phase
174.     double fase;
175.     double intPart;
176.     fase = modf(X / (PI/2), &intPart) * (PI/2);
177.     int period = static_cast<int>(intPart);
178.  
179.     // Make period positive by adding a large multiple of 4
180.     int periodMod4 = ((period % 4) + 4) % 4;
181.  
182.     // This is the table from your original description transformed into bit operations
183.     // Period 0,4,8... -> signP1=1, flagPi=0, signFase=1
184.     // Period 1,5,9... -> signP1=1, flagPi=1, signFase=-1
185.     // Period 2,6,10.. -> signP1=-1, flagPi=0, signFase=-1
186.     // Period 3,7,11.. -> signP1=-1, flagPi=1, signFase=1
187.  
188.     // Calculate bit flags without conditionals
189.     double signP1 = 1.0 - 2.0 * ((periodMod4 & 2) >> 1); // 1 for period 0,1 and -1 for period 2,3
190.     double flagPi = (PI/2) * (periodMod4 & 1);           // PI/2 for period 1,3 and 0 for period 0,2
191.  
192.     // For signFase, we need to check the pattern in your table
193.     // Period 0,3: signFase = 1
194.     // Period 1,2: signFase = -1
195.     // This can be computed with XOR
196.     double signFase = 1.0 - 2.0 * (((periodMod4 & 2) >> 1) ^ (periodMod4 & 1));
197.  
198.     // Calculate normalized X using absolute value of fase
199.     double normalizedX = signP1 * flagPi + signFase * std::abs(fase);
200.  
201.     // Convert back to degrees
202.     double result = normalizedX * 180.0 / PI;
203.  
204.     // Fix the sign to be consistent with the original table
205.     // If the absolute value is 0, make sure we return positive 0 (not -0)
206.     return (std::abs(result) < 1e-10) ? 0.0 : result;
207. }
208.  
209. // Test function
210. void testNormalization(double angle) {
211.     double normalized = NORMALIZE_ANGLE_90(angle);
212.     std::cout << std::setw(10) << angle << " degrees -> "
213.               << std::setw(10) << normalized << " degrees"
214.               << std::endl;
215. }
216.  
217. // Test the SSE2 implementation
218. void testSSE2Implementation() {
219.     std::cout << "\nTesting SSE2 implementation:\n";
220.     std::cout << "============================\n";
221.  
222.     const int num_angles = 8;
223.     double angles[num_angles] = {
224.         -450.0, -225.0, -90.0, -32.1,
225.         45.0, 180.0, 270.0, 345.0
226.     };
227.  
228.     double results[num_angles];
229.     double expected[num_angles];
230.  
231.     // Calculate expected results using the scalar function
232.     for (int i = 0; i < num_angles; i++) {
233.         expected[i] = normalizeAngle90_Bitwise(angles[i]);
234.     }
235.  
236.     // Run SSE2 implementation with exact matching
237.     normalizeAngles90_SSE2_Optimized(angles, results, num_angles);
238.  
239.     // Check results
240.     std::cout << "     Angle |      SSE2 |   Scalar | Difference\n";
241.     std::cout << "-------------------------------------------\n";
242.     bool allCorrect = true;
243.     for (int i = 0; i < num_angles; i++) {
244.         double diff = std::abs(results[i] - expected[i]);
245.         std::cout << std::setw(10) << angles[i] << " | "
246.                   << std::setw(10) << results[i] << " | "
247.                   << std::setw(10) << expected[i] << " | "
248.                   << std::scientific << std::setprecision(2) << diff;
249.  
250.         if (diff > 1e-10) {
251.             std::cout << " <-- MISMATCH";
252.             allCorrect = false;
253.         }
254.  
255.         std::cout << std::endl;
256.     }
257.  
258.     if (allCorrect) {
259.         std::cout << "\nSSE2 implementation produces correct results!\n";
260.     } else {
261.         std::cout << "\nSSE2 implementation has discrepancies!\n";
262.     }
263.  
264.     // Performance test
265.     const int performance_count = 1000000;
266.     double* test_angles = new double[performance_count];
267.     double* test_results = new double[performance_count];
268.  
269.     // Generate random angles
270.     for (int i = 0; i < performance_count; i++) {
271.         test_angles[i] = (rand() % 7200) - 3600; // -3600 to 3600 degrees
272.     }
273.  
274.     // Measure scalar performance (run first to warm up cache)
275.     clock_t scalar_start = clock();
276.     for (int i = 0; i < performance_count; i++) {
277.         test_results[i] = normalizeAngle90_Bitwise(test_angles[i]);
278.     }
279.     clock_t scalar_end = clock();
280.     double scalar_time = (double)(scalar_end - scalar_start) / CLOCKS_PER_SEC;
281.  
282.     // Measure SSE2 optimized performance - exact match version
283.     clock_t sse2_start = clock();
284.     normalizeAngles90_SSE2_Optimized(test_angles, test_results, performance_count);
285.     clock_t sse2_end = clock();
286.     double sse2_time = (double)(sse2_end - sse2_start) / CLOCKS_PER_SEC;
287.  
288.     // Measure SSE2 fast performance - may have tiny differences
289.     clock_t sse2_fast_start = clock();
290.     normalizeAngles90_SSE2_Fast(test_angles, test_results, performance_count);
291.     clock_t sse2_fast_end = clock();
292.     double sse2_fast_time = (double)(sse2_fast_end - sse2_fast_start) / CLOCKS_PER_SEC;
293.  
294.     std::cout << "\nPerformance comparison for " << performance_count << " angles:\n";
295.     std::cout << "Scalar time:           " << scalar_time << " seconds (baseline)\n";
296.     std::cout << "SSE2 exact-match time: " << sse2_time << " seconds (speed-up: "
297.               << scalar_time / sse2_time << "x)\n";
298.     std::cout << "SSE2 fast time:        " << sse2_fast_time << " seconds (speed-up: "
299.               << scalar_time / sse2_fast_time << "x)\n";
300.  
301.     delete[] test_angles;
302.     delete[] test_results;
303. }
304.  
305. int main() {
306.     std::cout << "Testing angle normalization to ±90 degrees range\n";
307.     std::cout << "================================================\n";
308.  
309.     // Test angles from -450 to 450 in steps of 45 degrees
310.     for (double angle = -450.0; angle <= 450.0; angle += 45.0) {
311.         testNormalization(angle);
312.     }
313.  
314.     // Test specific angles from the table
315.     std::cout << "\nTesting boundary cases:\n";
316.     testNormalization(-360.0);
317.     testNormalization(-270.0);
318.     testNormalization(-180.0);
319.     testNormalization(-90.0);
320.     testNormalization(0.0);
321.     testNormalization(90.0);
322.     testNormalization(180.0);
323.     testNormalization(270.0);
324.     testNormalization(360.0);
325.  
326.     // Test random angles
327.     std::cout << "\nTesting random angles:\n";
328.     testNormalization(-123.45);
329.     testNormalization(234.56);
330.     testNormalization(-32.1);
331.     testNormalization(75.9);
332.  
333.     // Show that macro and function versions produce the same results
334.     std::cout << "\nComparing macro vs function implementation:\n";
335.     double testAngle = 123.45;
336.     double macroResult = NORMALIZE_ANGLE_90(testAngle);
337.     double funcResult = normalizeAngle90_Bitwise(testAngle);
338.     std::cout << "Macro result: " << macroResult << " degrees\n";
339.     std::cout << "Func result:  " << funcResult << " degrees\n";
340.     std::cout << "Difference:   " << std::abs(macroResult - funcResult) << " degrees\n";
341.  
342.     // Test SSE2 implementation
343.     testSSE2Implementation();
344.  
345.     return 0;
346. }
347.  

Код (Text):

1. Testing angle normalization to ±90 degrees range
2. ================================================
3.       -450 degrees ->        -90 degrees
4.       -405 degrees ->         45 degrees
5.       -360 degrees ->          0 degrees
6.       -315 degrees ->         45 degrees
7.       -270 degrees ->         90 degrees
8.       -225 degrees ->        -45 degrees
9.       -180 degrees ->          0 degrees
10.       -135 degrees ->        -45 degrees
11.        -90 degrees ->        -90 degrees
12.        -45 degrees ->         45 degrees
13.          0 degrees ->          0 degrees
14.         45 degrees ->         45 degrees
15.         90 degrees ->         90 degrees
16.        135 degrees ->         45 degrees
17.        180 degrees ->          0 degrees
18.        225 degrees ->        -45 degrees
19.        270 degrees ->        -90 degrees
20.        315 degrees ->        -45 degrees
21.        360 degrees ->          0 degrees
22.        405 degrees ->         45 degrees
23.        450 degrees ->         90 degrees
24.  
25. Testing boundary cases:
26.       -360 degrees ->          0 degrees
27.       -270 degrees ->         90 degrees
28.       -180 degrees ->          0 degrees
29.        -90 degrees ->        -90 degrees
30.          0 degrees ->          0 degrees
31.         90 degrees ->         90 degrees
32.        180 degrees ->          0 degrees
33.        270 degrees ->        -90 degrees
34.        360 degrees ->          0 degrees
35.  
36. Testing random angles:
37.    -123.45 degrees ->     -56.55 degrees
38.     234.56 degrees ->     -54.56 degrees
39.      -32.1 degrees ->       32.1 degrees
40.       75.9 degrees ->       75.9 degrees
41.  
42. Comparing macro vs function implementation:
43. Macro result: 56.55 degrees
44. Func result:  56.55 degrees
45. Difference:   0 degrees
46.  
47. Testing SSE2 implementation:
48. ============================
49.      Angle |      SSE2 |   Scalar | Difference
50. -------------------------------------------
51.       -450 |        -90 |        -90 | 0.00e+00
52. -2.25e+02 |  -4.50e+01 |  -4.50e+01 | 0.00e+00
53. -9.00e+01 |  -9.00e+01 |  -9.00e+01 | 0.00e+00
54. -3.21e+01 |   3.21e+01 |   3.21e+01 | 0.00e+00
55.   4.50e+01 |   4.50e+01 |   4.50e+01 | 0.00e+00
56.   1.80e+02 |   0.00e+00 |   0.00e+00 | 0.00e+00
57.   2.70e+02 |  -9.00e+01 |  -9.00e+01 | 0.00e+00
58.   3.45e+02 |  -1.50e+01 |  -1.50e+01 | 0.00e+00
59.  
60. SSE2 implementation produces correct results!
61.  
62. Performance comparison for 1000000 angles:
63. Scalar time:           4.65e-02 seconds (baseline)
64. SSE2 exact-match time: 6.22e-02 seconds (speed-up: 7.48e-01x)
65. SSE2 fast time:        4.80e-02 seconds (speed-up: 9.69e-01x)
66.  

Scalar time быстрее ))

 

galenkane, вроде не плохо, но я уже разобрался, аналитику сам провёл. Формулу свёрху преобразовал в SSE...
У нас потоковое вычисление.
SIMD (англ. single instruction, multiple data — одиночный поток команд, множественный поток данных, ОКМД) — принцип компьютерных вычислений, позволяющий обеспечить параллелизм на уровне данных.
Это цитата из вики. Это значит мы не можем использовать столь привычные прогеру условные блоки, ещё не можем использовать таблицы отдельно для каждого элемента вектора, только для всего вектора. Так что приходится находит формулу которая может это решить. Скажу так, что для многих алгоритмов можно найти формулу, где не используются if, что и нужно для SIMD вычислений, в том числе например для вычислений на видеокарте.
Свой код проверил, работает правильно, на скорость ещё не проверял. Так же ещё можно попробовать задействовать обычные инструкции CPU, для этого можно извлечь данные из регистров xmm во временную переменную в секции .data, а потом а регистры CPU.
ЗЫ
Ещё можно нормализировать до 45,[math]-\displaystyle{\frac{\pi}{4} ... \frac{\pi}{4}}[/math]. Но не понятно насколько это повысит быстродействие, накладные затраты на нормализацию и так уже заметные.

 

Доделал. Теперь нужная точность достигается за меньшее кол. итераций, вот только быстродействие не слишком возросло, накладные расходы однако.
При норм. 90, максимальная точность достигается при 6 итераций, средняя дельта 0.00000022, максимальная 0.00000059, дальше точность перестаёт расти, наоборот ухудшается.
При норм. 180, точность максимальна уже при 8 итераций, средняя дельта 0.00000014, максимальная 0.00000044.
Зато если если нужна ограниченная точность, то при нормализации 90, уже достаточно 3 итераций. В общем, можно оптимизировать код нормализации, например если использовать временную переменную с секции дата и инструкции CPU, то вероятно можно ускорить код.
Кстати, почему метод называется cosf, sinf? Ну если в С++ классами реализовывать решение. А это чтобы подчеркнуть не высокую точность функции метода, можно ещё сделать методы cos, sin, где числа временно преобразуются в double для лучшей точности, но при этом время выполнения увеличится.

 

Попробовал, и... не прокатила.

Код (ASM):

1. ;возвращает в xmm0: x, в xmm1: signFase
2. NORMALIZE_ANGLE_90 MACRO
3.     ;;fase=modf(X, Pi); period=(int)(X/Pi);
4.     mulps   xmm0, vec4f@PI_div_2p   ;/=Pi/2;
5.     cvtps2dq xmm1, xmm0     ;(int)
6.     cvtdq2ps xmm3, xmm1     ;fPeriod=(float)period
7.     subps   xmm0, xmm3
8.     mulps   xmm0, vec4f@PI_div_2    ;*=Pi/2;
9.     movaps  xmm2, xmm1      ;flagPi
10. IF 1
11.     ;;signP1 = (float)((period & 2)^2)-1.0f;
12.     mov     eax, offset vec4i@Temp
13.     ASSUME  eax:ptr __m128i
14.     andps   xmm1, vec4i@2
15.     movaps  [eax], xmm1
16.     shl     [eax].i0, 30
17.     shl     [eax].i1, 30
18.     shl     [eax].i2, 30
19.     shl     [eax].i3, 30
20.     movaps  xmm1, [eax]
21.     ;;flagPi = (float)(period & 1)*PI/2 * signP1;
22.     andps   xmm2, vec4i@1   ;0,1,0,1...
23.     xorps   xmm2, vec4i@1   ;1,0,1,1...
24.     movaps  [eax], xmm2
25.     dec     [eax].i0
26.     dec     [eax].i1
27.     dec     [eax].i2
28.     dec     [eax].i3
29.     movaps  xmm2, [eax]
30.     andps   xmm2, vec4f@PI_div_2
31.     orps    xmm2, xmm1
32.     ;;signFase = (float)(((int)(-absf(fPeriod)) & 2)^2)-1.0f;   // -1.0f, 1.0f
33.     orps    xmm3, vec4i@BitSign     ;-absf(fPeriod)
34.     cvtps2dq xmm1, xmm3     ;(int)
35.     andps   xmm1, vec4i@2
36.     movaps  [eax], xmm1
37.     shl     [eax].i0, 30
38.     shl     [eax].i1, 30
39.     shl     [eax].i2, 30
40.     shl     [eax].i3, 30
41.     movaps  xmm1, [eax]
42.     ;fase = signFase*fase;
43.     xorps   xmm0, xmm1
44.     ;x = flagPi + fase;
45.     orps    xmm1, vec4f@1p0
46.     ASSUME  eax:nothing
47. ELSE
48.     ;;signP1 = (float)((period & 2)^2)-1.0f;
49.     andps   xmm1, vec4i@2   ;0,2,0,2...
50.     xorps   xmm1, vec4i@2   ;2,0,2,0...
51.     cvtdq2ps xmm1, xmm1
52.     subps   xmm1, vec4f@1p0 ;-= 1.0f
53.     ;;flagPi = (float)(period & 1)*PI/2 * signP1;
54.     andps   xmm2, vec4i@1   ;0,1,0,1...
55.     cvtdq2ps xmm2, xmm2
56.     mulps   xmm2, xmm1
57.     mulps   xmm2, vec4f@PI_div_2
58.     ;;signFase = (float)(((int)(-absf(fPeriod)) & 2)^2)-1.0f;   // -1.0f, 1.0f
59.     orps    xmm3, vec4i@BitSign     ;-absf(fPeriod)
60.     cvtps2dq xmm1, xmm3     ;(int)
61.     andps   xmm1, vec4i@2   ;0,2,0,2...
62.     xorps   xmm1, vec4i@2   ;2,0,2,0...
63.     cvtdq2ps xmm1, xmm1
64.     subps   xmm1, vec4f@1p0 ;-= 1.0f
65.     ;;x = flagPi + signFase*fase;
66.     mulps   xmm0, xmm1
67. ENDIF
68.     addps   xmm0, xmm2
69. ENDM

Наоборот, просело конкретно. Код чисто с SSE быстрей такой оптимизацией, показывает 20-28 тактов, вот так не равномерно, ну код хорошо нагружает ядро, именно SSE модуль, а мы знаем что ядро работает на два потока, а второй поток сбивает скорость первого. А вот код с инструкциями CPU сильно медленней, около 100 тактов, в 4-5 медленней. Хотя чего я ещё хочу, 4 синуса/косинуса за 20-25 тактов, это очень быстро.

 

Теперь надо тангенс, а для этого надо пересчитать константы, с этим есть сложности, но я разберусь.

 

а есть резон пойти во все тяжие и использовать ymm/zmm? +avx512

 

alex_dz, мне это пока не надо, разве что ли для гравитационной демки, типа куча звёзд летают вокруг друг друга.
Короче, тангенс TANPS чисто на float'е не получился подсчитать, тупо не хватает точности, ближе к 90 точность резко падает и дальнейшие увеличения больше 10 итерации не вызывает сходимости. Так что... хотя нафига, проще делать сразу TANPD, два числа double, и если надо, кастить до float.

 

Начал делать двойную точность, ну портировать код под double. И неожиданно столкнулся с проблемой. инструкция cvtps2dq работает не совсем так, как я ожидал, думал что два double преобразуются в два sqword, но нет, преобразуются в sdword. В результате код перестал работать, хотя его копипастой переделал. В общем, не понял описание...