Реализация sinpi () и cospi () с использованием стандартной математической библиотеки C

Функция sinpi(x) вычисляет sin (πx), а функция cospi(x) вычисляет cos (πx), где умножение с π неявно внутри функций. Эти функции были первоначально введены в стандартную математическую библиотеку C как расширение Sun Microsystems в конце 1980-х годов . IEEE Std 754 ™ -2008 определяет эквивалентные функции sinPi и cosPi в разделе 9.

Существуют многочисленные вычисления, где sin (πx) и cos (πx) встречаются естественным образом. Очень простой пример – преобразование Box-Muller (GEP Box и Mervin E. Muller, «Замечание о генерации случайных нормальных девиаций»). Annals of Mathematical Statistics , т. 29, № 2, стр. 610-611 ), который, учитывая две независимые случайные величины U₁ и U₂ с равномерным распределением, создает независимые случайные величины Z₁ и Z₂ со стандартным нормальным распределением:

 Z₁ = √(-2 ln U₁) cos (2 π U₂) Z₂ = √(-2 ln U₁) sin (2 π U₂) 

Другим примером является вычисление синуса и косинуса для аргументов степени, как при этом вычислении расстояния большого круга с использованием формулы Хаверсина:

 /* This function computes the great-circle distance of two points on earth using the Haversine formula, assuming spherical shape of the planet. A well-known numerical issue with the formula is reduced accuracy in the case of near antipodal points. lat1, lon1 latitude and longitude of first point, in degrees [-90,+90] lat2, lon2 latitude and longitude of second point, in degrees [-180,+180] radius radius of the earth in user-defined units, eg 6378.2 km or 3963.2 miles returns: distance of the two points, in the same units as radius Reference: http://en.wikipedia.org/wiki/Great-circle_distance */ double haversine (double lat1, double lon1, double lat2, double lon2, double radius) { double dlat, dlon, c1, c2, d1, d2, a, c, t; c1 = cospi (lat1 / 180.0); c2 = cospi (lat2 / 180.0); dlat = lat2 - lat1; dlon = lon2 - lon1; d1 = sinpi (dlat / 360.0); d2 = sinpi (dlon / 360.0); t = d2 * d2 * c1 * c2; a = d1 * d1 + t; c = 2.0 * asin (fmin (1.0, sqrt (a))); return radius * c; } 

Для C ++ библиотека Boost предоставляет функции sin_pi и cos_pi , а некоторые поставщики предлагают функции sinpi и cospi как расширения в системных библиотеках. Например, Apple добавила __sinpi , __cospi и соответствующие версии с одной точностью __sinpif , __cospif в iOS 7 и OS X 10.9 ( презентация , слайд 101). Но для многих других платформ нет реализации, доступной для программ C.

По сравнению с традиционным подходом, использующим, например, sin (M_PI * x) и cos (M_PI * x) , использование sinpi и cospi повышает точность за счет уменьшения ошибки округления посредством внутреннего умножения на π, а также дает преимущества производительности из-за большого количества более простое сокращение аргументов.

Как можно использовать стандартную математическую библиотеку C для реализации sinpi() и cospi() в разумно эффективном и стандартном соответствии?

    Для простоты я сосредоточусь на sincospi() , который одновременно обеспечивает результаты синуса и косинуса. sinpi и cospi затем могут быть сконструированы как функции-обертки, которые отбрасывают ненужные данные. Во многих приложениях обработка флагов с плавающей запятой (см. fenv.h ) не требуется, и нам не нужно сообщать об ошибках errno большую часть времени, поэтому я их не буду опускать.

    Основная алгоритмическая структура проста. Поскольку очень большие аргументы всегда являются целыми числами и поэтому, таким образом, кратно 2π, их значения синуса и косинуса хорошо известны. Другие аргументы складываются в диапазон [-¼, + ¼] при записи квадрантной информации. Полиномиальные минимаксные аппроксимации используются для вычисления синуса и косинуса на интервале первичных приближений. Наконец, данные квадранта используются для сопоставления предварительных результатов с конечным результатом путем циклического обмена результатами и смены знака.

    Правильная обработка специальных операндов (в частности, -0, бесконечности и NaN) требует от компилятора применения только оптимизаций, соответствующих правилам IEEE-754. Он не может преобразовывать x*0.0 в 0.0 (это неверно для -0, бесконечностей и NaN), а также не может оптимизировать 0.0-x в -x поскольку отрицание – это операция на уровне бит в соответствии с разделом 5.5.1 IEEE- 754 (что дает разные результаты для нhive и NaN). Большинство компиляторов предложит флаг, который обеспечивает использование «безопасных» преобразований, например -fp-model=precise для компилятора Intel C / C ++.

    Еще одно оговорка относится к использованию функции nearbyint при сокращении аргументов. Как и rint , эта функция задается для округления в соответствии с текущим режимом округления. Когда fenv.h не используется, по умолчанию режим округления fenv.h до «ближайшего или ровного». Когда он используется, существует риск, что действует режим направленного округления. Это можно было бы обойти с помощью round , который всегда обеспечивает режим округления «от округлой до ближайшей, галстуки от нуля» независимо от текущего режима округления. Однако эта функция будет работать медленнее, поскольку она не поддерживается эквивалентной машинной инструкцией на большинстве архитектур процессоров.

    Замечание о производительности: приведенный ниже код C99 в значительной степени зависит от использования fma() , который реализует операцию fma() добавления . На большинстве современных аппаратных архитектур это напрямую поддерживается соответствующей аппаратной инструкцией. Если это не так, код может испытывать значительное замедление из-за обычно медленной эмуляции FMA.

      #include  #include  /* Writes result sine result sin(πa) to the location pointed to by sp Writes result cosine result cos(πa) to the location pointed to by cp In extensive testing, no errors > 0.97 ulp were found in either the sine or cosine results, suggesting the results returned are faithfully rounded. */ void my_sincospi (double a, double *sp, double *cp) { double c, r, s, t, az; int64_t i; az = a * 0.0; // must be evaluated with IEEE-754 semantics /* for |a| >= 2**53, cospi(a) = 1.0, but cospi(Inf) = NaN */ a = (fabs (a) < 9.0071992547409920e+15) ? a : az; // 0x1.0p53 /* reduce argument to primary approximation interval (-0.25, 0.25) */ r = nearbyint (a + a); // must use IEEE-754 "to nearest" rounding i = (int64_t)r; t = fma (-0.5, r, a); /* compute core approximations */ s = t * t; /* Approximate cos(pi*x) for x in [-0.25,0.25] */ r = -1.0369917389758117e-4; r = fma (r, s, 1.9294935641298806e-3); r = fma (r, s, -2.5806887942825395e-2); r = fma (r, s, 2.3533063028328211e-1); r = fma (r, s, -1.3352627688538006e+0); r = fma (r, s, 4.0587121264167623e+0); r = fma (r, s, -4.9348022005446790e+0); c = fma (r, s, 1.0000000000000000e+0); /* Approximate sin(pi*x) for x in [-0.25,0.25] */ r = 4.6151442520157035e-4; r = fma (r, s, -7.3700183130883555e-3); r = fma (r, s, 8.2145868949323936e-2); r = fma (r, s, -5.9926452893214921e-1); r = fma (r, s, 2.5501640398732688e+0); r = fma (r, s, -5.1677127800499516e+0); s = s * t; r = r * s; s = fma (t, 3.1415926535897931e+0, r); /* map results according to quadrant */ if (i & 2) { s = 0.0 - s; // must be evaluated with IEEE-754 semantics c = 0.0 - c; // must be evaluated with IEEE-754 semantics } if (i & 1) { t = 0.0 - s; // must be evaluated with IEEE-754 semantics s = c; c = t; } /* IEEE-754: sinPi(+n) is +0 and sinPi(-n) is -0 for positive integers n */ if (a == floor (a)) s = az; *sp = s; *cp = c; } 

    Версия с одной точностью отличается в основном только в основных приближениях. Использование исчерпывающего тестирования позволяет точно определить границы ошибок.

     #include  #include  /* Writes result sine result sin(πa) to the location pointed to by sp Writes result cosine result cos(πa) to the location pointed to by cp In exhaustive testing, the maximum error in sine results was 0.96677 ulp, the maximum error in cosine results was 0.96563 ulp, meaning results are faithfully rounded. */ void my_sincospif (float a, float *sp, float *cp) { float az, t, c, r, s; int32_t i; az = a * 0.0f; // must be evaluated with IEEE-754 semantics /* for |a| > 2**24, cospi(a) = 1.0f, but cospi(Inf) = NaN */ a = (fabsf (a) < 0x1.0p24f) ? a : az; r = nearbyintf (a + a); // must use IEEE-754 "to nearest" rounding i = (int32_t)r; t = fmaf (-0.5f, r, a); /* compute core approximations */ s = t * t; /* Approximate cos(pi*x) for x in [-0.25,0.25] */ r = 0x1.d9e000p-3f; r = fmaf (r, s, -0x1.55c400p+0f); r = fmaf (r, s, 0x1.03c1cep+2f); r = fmaf (r, s, -0x1.3bd3ccp+2f); c = fmaf (r, s, 0x1.000000p+0f); /* Approximate sin(pi*x) for x in [-0.25,0.25] */ r = -0x1.310000p-1f; r = fmaf (r, s, 0x1.46737ep+1f); r = fmaf (r, s, -0x1.4abbfep+2f); r = (t * s) * r; s = fmaf (t, 0x1.921fb6p+1f, r); if (i & 2) { s = 0.0f - s; // must be evaluated with IEEE-754 semantics c = 0.0f - c; // must be evaluated with IEEE-754 semantics } if (i & 1) { t = 0.0f - s; // must be evaluated with IEEE-754 semantics s = c; c = t; } /* IEEE-754: sinPi(+n) is +0 and sinPi(-n) is -0 for positive integers n */ if (a == floorf (a)) s = az; *sp = s; *cp = c; }