1. Если вы только начинаете программировать на ассемблере и не знаете с чего начать, тогда попробуйте среду разработки ASM Visual IDE
    (c) на правах рекламы
    Скрыть объявление

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

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

  1. Jin X

    Jin X Active Member

    Публикаций:
    0
    Регистрация:
    15 янв 2009
    Сообщения:
    358
    Адрес:
    Кольца Сатурна
    Я ж замерил x64 на асме (fpu) и на Delphi (sse), вторая оказалась быстрее, см. выше.
    К тому же, я переделывал асм-код под x86, скорость была абсолютно той же, как и на x64, так что разница архитектур тут большой роли не играет.

    Пишем под x86 Sin(1), заходим под отладчиком (IDE-шным), видим:
    Код (Text):
    1. asm
    2.         FLD     tbyte ptr X
    3.         FSIN
    4.         FWAIT
    5. end;
    Пишем под x64 то же самое, видим:
    Код (Text):
    1.  
    2. function Sin(const X: Double): Double;
    3. var
    4.   Q: integer;
    5.   Y,Z: Double;
    6. begin
    7.   if Abs(x) < Pi/4 then
    8.     Result := pSinDouble(X, 0)
    9.   else
    10.   begin
    11.     Q := pRemDouble(X, Y, Z);
    12.     case Q of
    13.       0: Result :=  pSinDouble(Y, Z);
    14.       1: Result :=  pCosDouble(Y, Z);
    15.       2: Result := -pSinDouble(Y, Z);
    16.       3: Result := -pCosDouble(Y, Z);
    17.       else Result := 0; // avoid warning W1035 Return value of function '%s' might be undefined
    18.     end;
    19.   end;
    20. end;
    далее (pRemDouble неинтересен, идём в pCosDouble):
    Код (Text):
    1. {$IF    defined(ARITH_PUREPASCAL_EXT64) or defined(ARITH_X64_SSE)}
    2. function pCosDouble(const x, y: Double) : Double;
    3. const
    4.   CCos : ARRAY[0..5] OF UINT64 =
    5.   ( $BDA8FA6A8A7D84DF,
    6.     $3E21EE9DC12C88AC,
    7.     $BE927E4F7F1EE922,
    8.     $3EFA01A019C8F945,
    9.     $BF56C16C16C15018,
    10.     $3FA555555555554B );
    11. var
    12.   r1, r2, s, t, u, v,
    13.   L, L1, L2,
    14.   D2, D4 : Double;
    15. begin
    16.   D2 := x * x;
    17.   D4 := D2 * D2;
    18.   L1 :=           PDouble(@CCos[0])^;
    19.   L2 :=           PDouble(@CCos[1])^;
    20.   L1 := L1 * D4 + PDouble(@CCos[2])^;
    21.   L2 := L2 * D4 + PDouble(@CCos[3])^;
    22.   L1 := L1 * D4 + PDouble(@CCos[4])^;
    23.   L2 := L2 * D4 + PDouble(@CCos[5])^;
    24.   L := L2 + L1 * D2;
    25.   L := L * D4;
    26.   s := 1.0;
    27.   t := D2 * 0.5;
    28.   u := s - t;
    29.   v := u - s;
    30.   r1 := t + v;
    31.   r2 := x * y;
    32.   r2 := L - r2;
    33.   r2 := r2 - r1;
    34.   Result := u + r2;
    35. end;
    Смотрим дизасм в этом месте:
    https://www.screencast.com/t/pTo3pLs7Urz

    Вот ещё вараинт:
    Код (Text):
    1. {$APPTYPE CONSOLE}
    2. uses Winapi.Windows;
    3. var
    4.   C, i: Integer;
    5.   D, delta: Double;
    6. procedure FpuInit;
    7. asm
    8. FINIT
    9. end;
    10. function FpuSin: Double;
    11. asm
    12.         FLD     qword ptr D
    13.         FSIN
    14.         FWAIT
    15.         FSTP    st
    16. end;
    17. begin
    18.   D := 0;
    19.   delta := Pi / 18000000;
    20.   FpuInit;
    21.   C := GetTickCount;
    22.   for i := 1 to 36000000 do
    23.   begin
    24.     D := D + delta;
    25.     FpuSin();
    26.   end;
    27.   C := GetTickCount - C;
    28.   WriteLn(C);
    29. end.
    Компилим в x64 (чтоб всё по-честному было), получаем: 1125 мсек (против ≈ 850 обычного Sin на sse/Тейлор, тоже x64, исходники выше).
    --- Сообщение объединено, 1 авг 2020 ---
    Х/з, что там внутри и почему это медленнее.
    В исходнике Delphi, по меньшей мере, идёт чередование использования регистров (переменные L1, L2), которые вычисляются параллельно (можно сделать 3 шт вообще, по идее, должно быть ещё быстрее). Может, в харде все вычисления последовательны?
    Кстати, надо ещё на AMD замерить скорость, кстати...
    --- Сообщение объединено, 1 авг 2020 ---
    На AMD разница меньше: 1100 / 800 (вообще говоря, замер sse сильно скачет от 700 до 950 примерно).
    --- Сообщение объединено, 1 авг 2020 ---
    *разница НЕ меньше.
     
    M0rg0t, q2e74 и Indy_ нравится это.
  2. Indy_

    Indy_ Well-Known Member

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

    ~1.3 раза быстрее, это ниочём. Странно конечно почему последовательность инструкций выполняется быстрее чем одна fsin при том же результате. Фишка fpu - там есть стек, на котором идут вычисления, что делает не нужным выгрузку результата в память. Те на последовательности мат функций fpu будет быстрее. Судя по профайлу нет никакого смысла юзать левое(полиномы на sse), тем более для системы это не имеет значения - она выгружает мат блок через xsave.

    > pCosDouble

    Так это чебышев походу, нужно искать по коэффициентам. В точности совпадает с fsin ??
    --- Сообщение объединено, 1 авг 2020 ---
    Jin X,

    > FWAIT

    Эта инструкция доставляет мат фаулт, если его нет работает как nop. Если зациклить L: fwait/jmp L, то ядро будет использовать мат блок и общий профайл потока просядет. Это сам факт использования математики, mmx отображены на fpu, sse в том же блоке и в общем любо обращение к этим блокам включает механизм выгрузки контекста математики, а это долго ибо он большой.
    --- Сообщение объединено, 1 авг 2020 ---
    Jin X,

    На скрине L1 := CCos[0]

    А что это ?
    Какая то константа или может быть перед вычислением синуса вычисляется косинус, откуда ссылки на переменную ?
     
    q2e74 нравится это.
  3. UbIvItS

    UbIvItS Well-Known Member

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

    Indy_ Well-Known Member

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

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

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

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

    UbIvItS Well-Known Member

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

    Indy_ Well-Known Member

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

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

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

    UbIvItS Well-Known Member

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

    Indy_ Well-Known Member

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

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

    Jin X Active Member

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

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

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

    Indy_ Well-Known Member

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

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

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

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

    Jin X Active Member

    Публикаций:
    0
    Регистрация:
    15 янв 2009
    Сообщения:
    358
    Адрес:
    Кольца Сатурна
    По поводу точности.
    Из 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 в 22:13 ---
    Вот же: 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 в 22:28 ---
    Ну не то чтоб ни о чём, но не сильно много.
    Но я думаю, что если алгоритм усовершенствовать, можно добиться лучшего результата (до 2х, думаю, всё же можно довести). А если пожертвовать точностью (в довольно многих случаях большая точность не нужна), но ещё больше.
    --- Сообщение объединено, 2 авг 2020 в 22:44 ---
    А вот, кстати, и фокус-покус.
    Написал на 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 в 22:53 ---
    И ошибок тут 49438 из 36000000 (0,137 %).
    Разница больше, чем на 1, опять же, в 417 случаях :))
    Больше, чем на 16 – в 43-х случаях. Т.е. очень похожая ситуация, только общее число ошибок меньше.
    --- Сообщение объединено, 2 авг 2020 в 22:53 ---
    Можно ещё под шлангом проверить, но уже в лом что-то.
     

    Вложения:

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

    Indy_ Well-Known Member

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

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

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

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