Можноли как то ускорить код

[Missir]

Компилировал в gcc с O2, профилировал с помощью perf
слабое место судя по всему здесь, на него уходит 60% времени всей программы

   .............................
    g=sqrt(p/3); r=3*q/(g*p);     
    t=g*(0.51304+2.71*r-0.04115123*r*r+0.00739962*r*r*r);
    v=t*(t*t+p)+q; g=3*t*t+p; q=v/g-3*t*v*v/(g*g*g)+a*12;   
    p=q-a; e=q*p+b; h=p*p-4*e; r=sqrt(h);  g=(-p-r)/2; s=(-p+r)/2;
        
    return sqrt(q)+sqrt(g)+sqrt(s);
}

perf выдаёт аннотированный дезасемблированный код, но то ли я в нём не разобрался, то ли он некорректный.
Например, судя по нему, некоторые subsd оказываются очень ресурсоёмкими, а рядом стоящие sqrtd наоборот.
Есть подозрение, что на самом деле основная нагрузка это корни, просто pref неправильно определяет.

Пробовал заменять sqrt() на sqrtf() но толку от этого нет, только появляются дополнительные преобразования double <-> float.

Даже написал кастомный корень:

double sqrta(double x)
{                        
uint64_t y, z, r, mts, u=67108864;
const uint64_t mask_exp=9218868437227405312, b1=4503599627370496,
                   b511=2301339409586323456, mask=4503599627370495;
               
y = *((uint64_t*) &x); mts=((y & mask) | b1);

z=1470631420174624+((3951755717*(mts>>21))>>11)
  -((1405815*(mts>>32)*(mts>>32))>>10)
  +((55937*(mts>>37)*(mts>>37)*(mts>>37))>>12)
  -((2577*(mts>>41)*(mts>>41)*(mts>>41)*(mts>>41))>>9)
  +((806*(mts>>43)*(mts>>43)*(mts>>43)*(mts>>43)*(mts>>43))>>12);

if ((y & b1) != b1) mts=(((z>>21)*3037000500)>>10); else mts=z;  

y+=b1; y>>=1; y+=b511; y&=mask_exp; y+=(mts & mask);

 return *((double*)&y);
 }

Не сказать что он медленный, но похоже всё таки медленнее чем в <math.h>, ну и точность у него только 6 знаков, хотя мне бы этого хватило.
Не знаю, можно ли его как то допилить?

[mihaild]

Очевидный первый вопрос: у вас эти вычисления, вероятно, делаются для очень многих входов. Вычисления для следующего входа всегда требуют всего предыдущего, или большую часть времени несколько вызовов подряд независимы? Если второе, то можно попробовать переписать так чтобы компилятор начал активнее SSE использовать.

Ещё из очевидного попробуйте добавить флаги -ffast-math -march=native (даже если они не подходят для боевого применения, могут подсказать, куда смотреть).

[Missir]

mihaild в сообщении #1637462 писал(а):

Очевидный первый вопрос: у вас эти вычисления, вероятно, делаются для очень многих входов. Вычисления для следующего входа всегда требуют всего предыдущего

нет, это не рекуррентные вычисления, программа вызывает эту функцию, даёт ей данные и она обсчитывает, на выходе одно число double. Но как это использовать не понимаю, данные всё время разные. Если Вы о распараллеливании, то это делается элементарно, но у меня и нет с этим затруднений. У меня задача оптимизировать этот единственный вызов.
Ключ -march=native поставил, изменений никаких, -ffast-math уже был. Переход с О3 на О2 ни на что не влияет.
Сейчас попробовал профилировать в MSWC, там результаты более адекватные. Одиночный корень занимает 6.2%, примерно столько же многочлен (переписывал его по схеме Горнера, но изменений никаких). Строчка с делениями 5.7%, и ещё одна, примерно также, странно но в ней только вычитания и умножения. Тройной корень забирает только 8.2% а не 20 как можно было бы ожидать, видимо компилятор использует sse. В общем как я и думал, корни бурут 20% от всей программы + ещё 3 строки и выходит около 40%.
Можно это хозяйство вообще как то оптимизировать? Может перейти во float или использовать целочисленную арифметику? Это может что то дать, или не стоит даже пытаться?

[mihaild]

Missir в сообщении #1637464 писал(а):

Если Вы о распараллеливании, то это делается элементарно, но у меня и нет с этим затруднений

Не совсем, я именно про SSE (и AVX).

Missir в сообщении #1637464 писал(а):

переписывал его по схеме Горнера, но изменений никаких

Это логично, компилятор с такими оптимизациями должен сам справиться.

Missir в сообщении #1637464 писал(а):

Может перейти во float или использовать целочисленную арифметику?

float если вы упираетесь в CPU и не используется SIMD не поможет.

[photon]

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

[GAA]

photon в сообщении #1637741 писал(а):

Если большой, то можно попробовать заменить корень на вычисление величины обратной корню по известному алгоритму из Quake и затем делить единицу на эту величину

А для каких это процессоров? В начальном сообщении указаны инструкции SSE. Т.е. что-то очень старинное рассматривать не надо. Написать код работающий быстрее аппаратного разве возможно? Можно хотя бы грубые оценки в тактах?

upd
RSQRTSS xmm1, xmm2/m128
вычисляет приближенное значение обратной величины квадратного корня для младшего из 4-х упакованных коротких вещественных чисел из операнда-источника. Вычисленное значение помещается в младшее 32-битное поле операнда-назначения.
Относительная ошибка аппроксимации:$1.5\cdot 2^{-12}$.

Instruction Latency 4/5, Throughput 1
(Latency зависит от модели процессора)
upd 2

Intel® 64 and IA-32 Architectures Optimization Reference Manual Volume 2, 2023 писал(а):

The latency and throughput of SQRTPS/SQRTSS can vary with input values. For certain values, hardware can complete quickly, throughput may be as low as ~ 6 cycles. Similarly, latency for certain input values may be as low as less than 10 cycles

-- Wed 01.05.2024 21:18:17 --

А по теме. Были праздники. Можно было бы и посмотреть текст, но ТС его не выложил. В любом случае, если компилятор не слишком убогий, то обычно за счёт ручной оптимизации можно выиграть проценты (или доли процента).

-- Wed 01.05.2024 22:05:22 --

photon в сообщении #1637741 писал(а):

и затем делить единицу на эту величину

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

[GAA]

По поводу корня уже было: Как на компьютере вычисляются корни и др мат. функции?

[photon]

GAA в сообщении #1637756 писал(а):

А для каких это процессоров? В начальном сообщении указаны инструкции SSE. Т.е. что-то очень старинное рассматривать не надо

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

[GAA]

Компиляция под 64. Пусть соглашение о передаче аргументов функции и возврата результата следующее. Single-параметр передаётся функции через младшие 32 бита xmm0, функция возвращает Single-результат в младших 32 битах xmm0.

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

1) быстрый обратный квадратный корень (без одного шага метода Ньютона — Рафсона)
y = number;
i = * ( long * ) &y;
i = 0x5f3759df - ( i >> 1 );
y = * ( float * ) &i;
"подстрочный перевод" (без попыток оптимизации)

Используется синтаксис ASM

 movaps  xmm1, xmm0
 movd    xmm2, MagicConst
 psrld   xmm1, 1
 psubd   xmm2, xmm1
 mulss   xmm0, xmm2 

Здесь 32-битная переменная MagicConst содержит 0x5f3759df. Назовём эту функцию asqrt.

2) Использование rsqrtss

Используется синтаксис ASM

 rsqrtss xmm1, xmm0
 mulss   xmm0, xmm1 

Назовём эту функцию rsqrtss.

Для очень грубой оценки погрешности вычислим абсолютную и относительную погрешности для точек 0, 0.1, 0.2, … 10. Так как по величине погрешности этих двух функций сильно отличаются, то выведем на две оси ординат.
Абсолютная погрешность:

Вложение:

absolute_error.PNG

Относительная погрешность:

Вложение:

relative_error.PNG

И абсолютная, и относительная погрешности asqrt чудовищно проигрывают rsqrtss.

Теперь по скорости. Если не разбираться в заметке, на которую Вы привели ссылку, то в ней найдём

PatientZero в заметке Полезен ли сегодня быстрый обратный квадратный корень из Quake III? писал(а):

Результаты довольно схожи на всех машинах, время методов приблизительно ранжировано в порядке rsqrtps <= appr < exact <= appr_nr. Применение метода appr_nr или медленнее, или схоже с методом exact, то есть в данном случае не даёт реальной выгоды.

В тексте заметки код, основанный на rsqrtps, обозначен через rsqrtps; код, основанный на быстром обратном квадратном корне (без одного шага метода Ньютона — Рафсона), — через appr; код с итерацией методом Ньютона — Рафсона — через appr_nr.
Т.е. и для нулевого приближения, и для случая более точного вычисления «аппаратный вариант» выигрывает у быстрого обратного квадратного корня.

В заметке отмечается, что если требуется обеспечить одинаковые результаты и на Intel, и на AMD, то приближённый аппаратный обратный корень (rsqrtps/rsqrtss) не выйдет применять. Может это где-то для отладки и нужно. Однако учитывая точность быстрого обратного квадратного корня найти кроме графики его применение очень трудно. (А для графики, вроде, используются GPU).

Тестить такие короткие функции крайне трудно. И совершенно не видно причины: для практических расчётов (rsqrtps/rsqrtss) + одна итерация Ньютона — Рафсона.

-- Fri 03.05.2024 18:17:34 --

Такие же выводы (rsqrtps/rsqrtss + одна итерация Ньютона — Рафсона) можно найти в заметке, указанной в Википедии: Timing square root.

-- Fri 03.05.2024 18:22:35 --

Возвращаясь к вопросу начального сообщения треда.

Missir в сообщении #1637460 писал(а):

sqrt(q)+sqrt(g)+sqrt(s)

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

[Missir]

photon в сообщении #1637741 писал(а):

по известному алгоритму из Quake и затем делить единицу на эту величину

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

double fast_sqrt(double x)
{ 
  double y=0;
  asm ("movq %1, %%rax; shr  $1, %%rax; movq $6910469321099104000, %%rbx; "
       "subq %%rax, %%rbx; movq %%rbx, %0;"
       :"=r"(y)
       :"r"(x)
       :"%rax","%rbx"
       );
       y=y*(3-(x*y*y));
       y=y*(12-(x*y*y));
       y=y*(768-(x*y*y));
       
  return y*x/8192;
}

даёт точность почти double, если убрать последнюю итерацию - получается почти float, если оставить только одну, то получается 3 верных знака. Но но скорости это хозяйство отстаёт от sqrt. Если ограничиться точностью float и вставить тело функции прямо в код, то работает в 1.5 раза быстрее, что тоже несущественно.

-- 04.05.2024, 18:50 --

photon в сообщении #1637741 писал(а):

можно сделать LUT и считать результат линейной интерполяцией

я не знаю что это такое, но делал вот такую штуку, если Вы об этом:

   const uint64_t root[64]=
   { 0,69827414112268,138604612414727,206377602297819,
       273189127023833,339078981026126,404084287098361,
       468239740956932,531577827751080,594129014356004,
       655921920679557,716983472715828,777339039667879,
       837012557120638,896026637960214,954402672497414,
       1012160919052590,1069320586089555,1125899906842624,
       1181916207258746,1237385967972342,1292324880941192,
       1346747901294976,1400669294881971,1454102681942250,
       1507061077286234,1559556927314419,1611602144176602,
       1663208137336207,1714385842776648,1765145750061507,
       1815497927438199,1865452045155277,1964202921218288,
       2061468567839650,2157314049294893,2251799813685248,
       2344982138838073,2436913524311641,2527643037258634,
       2617216618617683,2705677355056896,2793065721238116,
       2879419796267360,2964775457615668,3049166555311868,
       3132625068806179,3215181248566246,3296863744183466,
       3377699720527872,3457714963286720,3536933975049166,
       3615380062951923,3693075418774505,3770041192264134,
       3846297558376920,3921863779041097,3996758259978076,
       4070998603056240,4144601654599387,4217583550025434,
       4289959755150440,4361745104457490,4432953836598623 };

  
   asm ("movq %1, %%rax;"
        "movq %%rax, %%rbx;"
        "xorq %%rcx, %%rcx; incq %%rcx; rcr  $2, %%rcx;"  // знаковый бит
        "andq %%rcx, %%rbx; subq %%rbx, %%rax;"
        "shr  $2, %%rcx; decq %%rcx; movq %%rcx, %%rdx; incq %%rcx;"
        "shr  $9, %%rcx; addq %%rcx, %%rax;"
        "decq %%rcx; notq %%rcx; andq %%rcx, %%rdx; shr $1, %%rax;"
        "addq %%rdx, %%rax; addq %%rbx, %%rax;"
        "movq %%rax, %%rbx; andq %%rcx, %%rbx;"            // порядок
        "notq %%rcx; andq %%rcx, %%rax;"
        "shr  $38, %%rax; xor %%al, %%al; shr  $5, %%rax;" // 8-битный адрес мантиссы        
        "movq %2, %%rcx; add %%rax, %%rcx; xorq (%%rcx), %%rbx;" // из массива  
        "movq %%rbx, %0;"
                
         :"=r"(br)                        
         :"r"(ar), "r"(root)            
         :"%rax", "%rbx", "%rcx", "%rdx"   
        );
       tmp=ar/br; br+=tmp; br/=2;

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

[photon]

Что ж: факир был пьян и фокус не удался

[Missir]

Вообще, нужно найти максимум из множества $sqrt(x)+sqrt(y)+sqrt(z)$, можно было бы сначала использовать
$x+y+z < (sqrt(x)+sqrt(y)+sqrt(z))^2<3(x+y+z)$ для отсева бесперспективных вариантов, но $x+y+z$ различаются незначительно, так что отделить ничего не получается

[Missir]

mihaild в сообщении #1637535 писал(а):

Не совсем, я именно про SSE (и AVX)

сейчас до меня дошло как можно сделать - промежуточные значения 4-х соседних шагов цикла переводить в float и грузить в регистры sse, а потом вместе обрабатывать в пакетном режиме. Сложностей с этим особо не вижу. Единственный вопрос,
mihaild, не подскажете, как мне в завершении всех расчётом найти максимальный элемент 4-х элементного sse регистра? (собственно в этом и вся цель).

[mihaild]

Missir в сообщении #1638048 писал(а):

mihaild, не подскажете, как мне в завершении всех расчётом найти максимальный элемент 4-х элементного sse регистра?

Т.е. у Вас есть 4 float32 числа в xmm регистре, и Вы хотите найти максимум? Такие операции называются горизонтальными, но я не уверен, что горизонтальный максимум существует (см. https://www.intel.com/content/www/us/en ... index.html), и они обычно медленные.
А не получится поменять порядок - сначала посчитать 4 четверки, а потом в них всех искать максимумы одновременно (это за 3 вызова _mm_max_ps даст результат; т.к. там latency сильно больше, чем throughtput, а вычисления сильно зависимы, то лучше, если получится, считать сразу для 8 четверок).

[GAA]

Пусть 4 single находятся в регистре xmm0, доступно два регистра (пусть xmm0 и xmm1), a результат нужно вернуть в младших 32 xmm0.

Используется синтаксис ASM

                                              // xmm0 = (A, B, C, D)
  pshufd xmm1, xmm0, 01001110b                // xmm1 = (C, D, A, B)
  maxps  xmm0, xmm1                          // xmm0 = (max (A, C), max (B, D), max(A, C), max(B, D)
  pshufd xmm1, xmm0, 10110001b              // xmm1 = (max (B, D), max (A, C), max(B, D), max(A, C))
  maxps  xmm0, xmm1