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

Тема в разделе "WASM.ASSEMBLER", создана пользователем Mikl___, 14 июл 2020.

  1. UbIvItS

    UbIvItS Well-Known Member

    Публикаций:
    0
    Регистрация:
    5 янв 2007
    Сообщения:
    5.442
    пч фпу кривой и его лишь для обратной совместимости держат. видать, не стали мутить приличную хард реализацию, чтоб сохранить место на плашке под "широкие" функи проца :)
     
  2. Indy_

    Indy_ Well-Known Member

    Публикаций:
    4
    Регистрация:
    29 апр 2011
    Сообщения:
    4.787
    UbIvItS,

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

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

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

    UbIvItS Well-Known Member

    Публикаций:
    0
    Регистрация:
    5 янв 2007
    Сообщения:
    5.442
    хз, что там пользовали (можЬ в манах Ынтеля написано). но в общем и целом для харда лучше кордик юзать == схема более компактная и можно хорошо балансировать погрешность.
     
  4. Indy_

    Indy_ Well-Known Member

    Публикаций:
    4
    Регистрация:
    29 апр 2011
    Сообщения:
    4.787
    UbIvItS,

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

    Что такое кордик" ?
     
  5. UbIvItS

    UbIvItS Well-Known Member

    Публикаций:
    0
    Регистрация:
    5 янв 2007
    Сообщения:
    5.442
  6. Indy_

    Indy_ Well-Known Member

    Публикаций:
    4
    Регистрация:
    29 апр 2011
    Сообщения:
    4.787
    UbIvItS,

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

    Jin X Active Member

    Публикаций:
    0
    Регистрация:
    15 янв 2009
    Сообщения:
    365
    Адрес:
    Кольца Сатурна
    Я убрал fwait, результат тот же. По крайней мере, изменения чисел не заметны на глаз.

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

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

    Indy_ Well-Known Member

    Публикаций:
    4
    Регистрация:
    29 апр 2011
    Сообщения:
    4.787
    Jin X,

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

    TermoSINteZ Синоби даоса Команда форума

    Публикаций:
    1
    Регистрация:
    11 июн 2004
    Сообщения:
    3.495
    Адрес:
    Russia
    maalchemist, чистый синус лучше не замерять. не будет адекватного сравнению. напишите формулу где юзается синус. чтоб были другие инструкции но желательно чтоб отличалась только реализация синуса. Тогда будет более валидный тест. Но тут сами решайте нужно ли вам это )
    просто мерить чистый синус - не валидно как с точки зрения общих алгоритмов работы конвеера процессора, так и в плане практических задач.
     
  10. Jin X

    Jin X Active Member

    Публикаций:
    0
    Регистрация:
    15 янв 2009
    Сообщения:
    365
    Адрес:
    Кольца Сатурна
    По поводу точности.
    Из 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 ---
    Можно ещё под шлангом проверить, но уже в лом что-то.
     

    Вложения:

    • System.zip
      Размер файла:
      222,6 КБ
      Просмотров:
      286
    Последнее редактирование: 2 авг 2020
  11. Indy_

    Indy_ Well-Known Member

    Публикаций:
    4
    Регистрация:
    29 апр 2011
    Сообщения:
    4.787
    Jin X,

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

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

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

    Intro Active Member

    Публикаций:
    0
    Регистрация:
    29 авг 2009
    Сообщения:
    347
    Не совсем по теме.
    Переделка методов класса 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!
     
    Jin X нравится это.
  13. НетРегистрации

    НетРегистрации Member

    Публикаций:
    0
    Регистрация:
    1 фев 2020
    Сообщения:
    58
    Код (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
     
    Jin X нравится это.
  14. Intro

    Intro Active Member

    Публикаций:
    0
    Регистрация:
    29 авг 2009
    Сообщения:
    347
    НетРегистрации, прикольно. Я тоже думал что с константами возможна дополнительная оптимизация, но не смог математически доказать правильность.
    Да, оказалось надо больше итераций, а так же ограничить х от -PI до PI, на что нужен дополнительный код, так что тут не всё просто. А с другой стороны, что нам запрещает предпринять специальные методы, чтобы угол ограничить -180..180, а если что пошло не так, сработает ассерт.
     
  15. Jin X

    Jin X Active Member

    Публикаций:
    0
    Регистрация:
    15 янв 2009
    Сообщения:
    365
    Адрес:
    Кольца Сатурна
    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 часа ночи.

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

    Вложения:

    • sinps.asm
      Размер файла:
      12,6 КБ
      Просмотров:
      119
  16. R81...

    R81... Member

    Публикаций:
    0
    Регистрация:
    1 фев 2020
    Сообщения:
    61
    Обычная точность ~3байта, внутренняя SSE мне неизвестна, 1/15! ну вы поняли, что может быть при суммировании.
    Соответственно можно лишь попытаться улучшить точность, начиная суммирование с конца.
     
  17. Jin X

    Jin X Active Member

    Публикаций:
    0
    Регистрация:
    15 янв 2009
    Сообщения:
    365
    Адрес:
    Кольца Сатурна
    Можно ещё через умножение сделать.
    Т.е. не складывать 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...
     

    Вложения:

    • sinps.asm
      Размер файла:
      12,9 КБ
      Просмотров:
      123
    Pavia и R81... нравится это.
  18. Jin X

    Jin X Active Member

    Публикаций:
    0
    Регистрация:
    15 янв 2009
    Сообщения:
    365
    Адрес:
    Кольца Сатурна
    О как!
    Вынес fld1 в начало кода (после finit), скорость fsin увеличилась до 66 тактов.
    И вообще, все цифры чуть уменьшились... почему-то :dntknw:
    Код (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, разницы почти не будет.
    Интересно, почему так?
     
  19. Indy_

    Indy_ Well-Known Member

    Публикаций:
    4
    Регистрация:
    29 апр 2011
    Сообщения:
    4.787
    Jin X,

    Так потому что это ошибка профайлера(системных таймеров), у него разрешение на много порядков ниже, а есчо он периодически вытисняется как задача. Замер должен быть в цикле на очень большом числе итераций. Иначе будет рандом.
    --- Сообщение объединено, 22 дек 2020 ---
    > Success: both RDTSCP instruction and invariant TSC are supported.

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

    Jin X Active Member

    Публикаций:
    0
    Регистрация:
    15 янв 2009
    Сообщения:
    365
    Адрес:
    Кольца Сатурна
    Так, там 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 тут не должно быть, случайно написал :-/
     

    Вложения:

    • sinps.asm
      Размер файла:
      12,9 КБ
      Просмотров:
      87
    Indy_ нравится это.