Implémentation de sinpi () et cospi () à l’aide de la bibliothèque mathématique C standard

La fonction sinpi(x) calcule sin (πx) et la fonction cospi(x) calcule cos (πx), où la multiplication avec π est implicite à l’intérieur des fonctions. Ces fonctions ont été initialement introduites dans la bibliothèque mathématique standard C en tant qu’extension de Sun Microsystems à la fin des années 1980 . IEEE Std 754 ™ -2008 spécifie les fonctions équivalentes sinPi et cosPi dans la section 9.

Il existe de nombreux calculs où sin (πx) et cos (πx) se produisent naturellement. Un exemple très simple est la transformation de Box-Muller (GEP Box et Mervin E. Muller, “Note sur la génération de déviations normales aléatoires”. Annals of Mathematical Statistics , vol. 29, n ° 2, p. 610 – 611 ), qui, à partir de deux variables aléatoires indépendantes U et U₂ à dissortingbution uniforme, produit des variables aléatoires indépendantes Z et Z₂ à dissortingbution normale normale:

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

Un autre exemple est le calcul des arguments sinus et cosinus pour les degrés, comme dans ce calcul de la distance de grand cercle utilisant la formule de Haversine:

 /* 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; } 

Pour C ++, la bibliothèque Boost fournit sin_pi et cos_pi , et certains fournisseurs proposent des fonctionnalités sinpi et cospi comme extensions dans les bibliothèques système. Par exemple, Apple a ajouté __sinpi , __cospi et les versions à simple précision correspondantes __sinpif , __cospif à iOS 7 et OS X 10.9 ( présentation , diapositive 101). Mais pour de nombreuses autres plates-formes, aucune mise en œuvre n’est facilement accessible aux programmes C.

Comparé à une approche traditionnelle qui utilise par exemple sin (M_PI * x) et cos (M_PI * x) , l’utilisation de sinpi et cospi améliore la précision en réduisant les erreurs d’arrondi via la multiplication interne avec π, et offre également des avantages en cospi performances du fait de la réduction des arguments plus simple.

Comment utiliser la bibliothèque de maths standard C pour implémenter les fonctionnalités sinpi() et cospi() de manière relativement efficace et conforme aux normes?

Pour des raisons de simplicité, je vais me concentrer sur sincospi() , qui fournit simultanément les résultats sinus et cosinus. sinpi et cospi peuvent ensuite être construits en tant que fonctions wrapper qui éliminent les données inutiles. Dans de nombreuses applications, la gestion des indicateurs de virgule flottante (voir fenv.h ) n’est pas requirejse, et nous n’avons pas besoin de signaler la plupart du temps les erreurs errno . C’est pourquoi je les omettrai.

La structure algorithmique de base est simple. Comme les très grands arguments sont toujours des entiers pairs, et donc des multiples de 2π, leurs valeurs de sinus et de cosinus sont bien connues. Les autres arguments sont repliés dans l’intervalle [-¼, + ¼] lors de l’enregistrement des informations sur le quadrant. Les approximations polynomiales de minimax sont utilisées pour calculer les sinus et les cosinus sur l’intervalle d’approximation principal. Enfin, les données du quadrant sont utilisées pour mapper les résultats préliminaires au résultat final par échange cyclique des résultats et changement de signe.

Le traitement correct des opérandes spéciaux (en particulier -0, infinities et NaN) nécessite que le compilateur applique uniquement les optimisations conformes aux règles IEEE-754. Il ne peut pas transformer x*0.0 en 0.0 (ceci n’est pas correct pour -0, les infinités et les NaN) et il ne peut pas optimiser 0.0-x en -x car la négation est une opération au niveau du bit conformément à la section 5.5.1 de IEEE- 754 (donnant des résultats différents pour les zéros et les NaN). La plupart des compilateurs proposent un indicateur qui impose l’utilisation de transformations “sûres”, par exemple -fp-model=precise pour le compilateur Intel C / C ++.

Une mise en garde supplémentaire concerne l’utilisation de la fonction nearbyint lors de la réduction des arguments. Comme rint , cette fonction est définie pour arrondir en fonction du mode d’arrondi actuel. Lorsque fenv.h n’est pas utilisé, le mode d’arrondi est défini par défaut sur “arrondi au plus proche ou même”. Lorsqu’il est utilisé, un mode d’arrondi dirigé risque d’être appliqué. Cela pourrait être contourné par l’utilisation de round , qui fournit toujours le mode d’arrondi “arrondi au plus proche, égalité loin de zéro”, indépendamment du mode d’arrondi actuel. Toutefois, cette fonction aura tendance à être plus lente car elle n’est pas prise en charge par une instruction machine équivalente sur la plupart des architectures de processeur.

Remarque sur les performances: le code C99 ci-dessous repose largement sur l’utilisation de fma() , qui implémente une opération de fusion multiplic-add . Sur la plupart des architectures matérielles modernes, ceci est directement pris en charge par une instruction matérielle correspondante. Si ce n’est pas le cas, le code risque de connaître un ralentissement important en raison de l’émulation généralement lente de 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; } 

La version à simple précision ne diffère fondamentalement que dans les approximations principales. L'utilisation de tests exhaustifs permet de déterminer avec précision les limites d'erreur.

 #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; }