Почему sinf() так медленно при компиляции моего кода в Linux/ARMv7?

Question

Тот же C-код с математикой sin() работает намного медленнее, когда компилирую свою программу на ARMv7 по сравнению с тем, что я получаю в Windows. Я компилирую с -O2 -Wall -Wextra -mfpu=neon -mtune=cortex-a9 -march=armv7 -std=c++11, а мой gcc - gcc (Ubuntu/Linaro 4.8.2-19ubuntu1) 4.8.2.

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

Моя программа создает несколько волновых таблиц при запуске, и в то время как она начинается практически мгновенно на Windows, она занимает около 25-30 секунд, чтобы начать на Linux/ARM ...

Вот код, который показывает, где sinf() функция используется, которая замедляет все.

for (int n = 0; n < 73; ++n) 
{ 
    // Max number of harmonics 
    int hrm = int(16000.f/twf[n]); 

    // Set vectors 
    basic_wf.assign(wavelength[n], 0); 

    for (int i = 0; i < wavelength[n]; ++i) 
    { 
     // Add harmonics 
     for (int h = 1; h < hrm; ++h) 
     { 
      const float harm = 0.14f * (sinf((float)i * FACTOR * twf[n] * (float)h)/(float)h); 
      if (h % 2 == 0) basic_wf[i] -= harm; // add even negative harmonic 
      else basic_wf[i] += harm;    // add odd positive harmonic 
     } 
    } 
} 

Здесь я заполняю 73 стола с пилообразной формой сигнала, добавляя необходимое количество гармоник для каждой частоты. Чем ниже высота ноты, тем выше количество гармоник (фактические вычисления sin()). Это почти мгновенно запускается в Windows ... занимает всю жизнь на моем Linux-сервере.

Answer 1

Вы могли бы сделать то, что уже сделали Napier и Co., чтобы вычислить таблицы логарифмов - или, точнее, таблицы для полномочий 1.000001 или аналогичные.

Если вам нужен вектор значений sin(k*w) затем вычислить c1000=cos(1000*w) и s1000=sin(1000*w) установите

c[0] = 1; s[0] = 0; 
c[1000]=c1000; s[1000] = s1000; 

, а затем итеративно

c[1000*(k+1)] = c1000*c[1000*k]-s1000*s[1000*k]; 
s[1000*(k+1)] = c1000*s[1000*k]+s1000*c[1000*k]; 

, а затем заполнить пробелы с помощью c1=cos(w) и s1=sin(w) снова используя тригонометрические тождества , 1000 шагов вперед или если вы хотите 500 вперед и 500 назад. Этот многоуровневый подход должен держать накопление ошибок с плавающей точкой достаточно малыми.

На «больших» процессорах не должно быть большой разницы, стоимость 2 умножения и оценка sincos должны быть сопоставимы. В вашем случае должно быть некоторое преимущество в мультипликативной процедуре.

Answer 2

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

Вы упомянули в комментарии, что операнды до sinf() включают множитель π. Это означает, что вы действительно хотите использовать sinpif(), где sinpi (x) = sin (x * π). Функция sinpi была введена в стандарте IEEE-754 (2008) с плавающей запятой, но пока не дошла до языковых стандартов. Однако несколько цепей инструментов предлагают его как расширение. Преимущество sinpi() состоит в том, что он требует только очень простого сокращения аргументов независимо от величины аргументов, что может значительно сократить время извлечения. Это приводит к повышению производительности. Так как умножение на π неявно, оно также может обеспечить улучшенную точность по дискретному подходу с использованием sinf().

Я показываю примерную реализацию C99 sinpif() ниже. Обратите внимание, что этот код в значительной степени зависит от стандартной математической функции fmaf() для достижения высокой скорости обработки и отличной точности. Если ваш процессор не имеет аппаратной поддержки для плавной операции многократного добавления (FMA), эта функция будет выполняться очень медленно, так как правильная эмуляция fmaf() нетривиальна. Так как код написан модульным способом, вы хотите настроить компилятор для применения максимального количества функций вложения или добавить соответствующие атрибуты вложения для всех составных функций.

Как вы указали, что ваша аппаратная платформа не предлагает встроенную поддержку FMA, вы можете заменить каждый fmaf(a,b,c) на (a*b+c) с некоторой потерей точности. Согласно моим тестам, максимальная ошибка ulp увеличивается до 1,71364 ulps. Это все еще очень хорошо, но my_sinf() больше не будет точно округлено больше в этом случае, что обычно считается желательным свойством.

/* Argument reduction for sinpi, cospi, sincospi. Reduces to [-0.25, +0.25] */ 
float trig_red_pi_f (float a, int *i) 
{ 
    float r; 
    r = rintf (a + a); 
    *i = (int)r; 
    r = a - 0.5f * r; 
    return r; 
} 

/* Approximate cos(pi*x) for x in [-0.25,0.25]. Maximum ulp error = 0.87440 */ 
float cospif_poly (float s) 
{ 
    float r; 
    r =    0x1.d98dcep-3f; // 2.31227502e-1f 
    r = fmaf (r, s, -0x1.55c4e8p+0f); // -1.33503580e+0f 
    r = fmaf (r, s, 0x1.03c1d4p+2f); // 4.05870533e+0f 
    r = fmaf (r, s, -0x1.3bd3ccp+2f); // -4.93480206e+0f 
    r = fmaf (r, s, 0x1.000000p+0f); // 1.00000000e+0f 
    return r; 
} 

/* Approximate sin(pi*x) for x in [-0.25,0.25]. Maximum ulp error = 0.96441 */ 
float sinpif_poly (float a, float s) 
{ 
    float r; 
    r =    -0x1.2dc6f8p-1f; // -5.89408636e-1f 
    r = fmaf (r, s, 0x1.46602ep+1f); // 2.54981017e+0f 
    r = fmaf (r, s, -0x1.4abbc0p+2f); // -5.16770935e+0f 
    r = r * s; 
    r = fmaf (r, a, -0x1.777a5cp-24f * a); // PI_lo // -8.74227766e-8f 
    r = fmaf (a, 0x1.921fb6p+1f, r);  // PI_hi // 3.14159274e+0f 
    return r; 
} 

/* Compute sin(pi*x) and cos(pi*x) based on quadrant */ 
float sinpif_cospif_core (float a, int i) 
{ 
    float r, s; 
    s = a * a; 
    r = (i & 1) ? cospif_poly (s) : sinpif_poly (a, s); 
    if (i & 2) { 
     r = 0.0f - r; // don't change "sign" of NaNs or create negative zeros 
    } 
    return r; 
} 

/* maximum ulp error = 0.96411 */ 
float my_sinpif (float a) 
{ 
    float r; 
    int i; 
    r = trig_red_pi_f (a, &i); 
    r = sinpif_cospif_core (r, i); 
    /* IEEE-754: sinPi(+n) is +0 and sinPi(-n) is -0 for positive integers n */ 
    r = (a == truncf (a)) ? (a * 0.0f) : r; 
    return r; 
}