Algorithme Fast Arc Cos?

J’ai ma propre fonction cos très rapide:

float sine(float x) { const float B = 4/pi; const float C = -4/(pi*pi); float y = B * x + C * x * abs(x); // const float Q = 0.775; const float P = 0.225; y = P * (y * abs(y) - y) + y; // Q * y + P * y * abs(y) return y; } float cosine(float x) { return sine(x + (pi / 2)); } 

Mais maintenant, quand je profile, je vois que acos () tue le processeur. Je n’ai pas besoin de précision intense. Quel est un moyen rapide de calculer acos (x) Merci.

Une approximation cubique simple, le polynôme de Lagrange pour x ∈ {-1, -½, 0, ½, 1}, est:

 double acos(x) { return (-0.69813170079773212 * x * x - 0.87266462599716477) * x + 1.5707963267948966; } 

Il a une erreur maximale d’environ 0,18 rad.

Vous avez de la mémoire disponible? Une table de recherche (avec interpolation, si nécessaire) sera la plus rapide.

nVidia dispose d’excellentes ressources qui montrent comment approximer des fonctions mathématiques très coûteuses, telles que: acos asin atan2 etc etc …

Ces algorithmes produisent de bons résultats lorsque la rapidité d’exécution est plus importante (dans des limites raisonnables) que la précision. Voici leur fonction acos:

 // Absolute error <= 6.7e-5 float acos(float x) { float negate = float(x < 0); x = abs(x); float ret = -0.0187293; ret = ret * x; ret = ret + 0.0742610; ret = ret * x; ret = ret - 0.2121144; ret = ret * x; ret = ret + 1.5707288; ret = ret * sqrt(1.0-x); ret = ret - 2 * negate * ret; return negate * 3.14159265358979 + ret; } 

Et voici les résultats pour le calcul des acos (0.5):

 nVidia: result: 1.0471513828611643 math.h: result: 1.0471975511965976 

C'est assez proche! Selon votre degré de précision requirejs, cela pourrait être une bonne option pour vous.

J’ai le mien. C’est assez précis et assez rapide. Cela fonctionne à partir d’un théorème que j’ai construit autour de la convergence quartique. C’est vraiment intéressant, et vous pouvez voir l’équation et à quelle vitesse cela peut faire converger mon approximation de log naturelle ici: https://www.desmos.com/calculator/yb04qt8jx4

Voici mon code arccos:

 function acos(x) local a=1.43+0.59*xa=(a+(2+2*x)/a)/2 local b=1.65-1.41*xb=(b+(2-2*x)/b)/2 local c=0.88-0.77*xc=(c+(2-a)/c)/2 return (8*(c+(2-a)/c)-(b+(2-2*x)/b))/6 end 

Une grande partie de cela est juste approximation de la racine carrée. Cela fonctionne très bien aussi, sauf si vous vous approchez trop de la racine carrée de 0. L’erreur moyenne (à l’exclusion de x = 0,99 à 1) est 0,0003. Le problème, cependant, est qu’à 0,99, il commence à aller chier, et à x = 1, la différence de précision devient 0,05. Bien sûr, cela pourrait être résolu en faisant plus d’itérations sur les racines carrées (lol nope) ou, juste un peu comme, si x> 0.99, utilisez un ensemble différent de linéarisations de racines carrées, mais cela rend le code très long et moche. .

Si vous ne vous souciez pas autant de la précision, vous pouvez simplement faire une itération par racine carrée, ce qui devrait vous garder quelque part dans la plage de 0,0162 ou quelque chose au niveau de la précision:

 function acos(x) local a=1.43+0.59*xa=(a+(2+2*x)/a)/2 local b=1.65-1.41*xb=(b+(2-2*x)/b)/2 local c=0.88-0.77*xc=(c+(2-a)/c)/2 return 8/3*cb/3 end 

Si cela vous convient, vous pouvez utiliser du code racine carrée préexistant. Cela éliminera l’équation un peu folle à x = 1:

 function acos(x) local a = math.sqrt(2+2*x) local b = math.sqrt(2-2*x) local c = math.sqrt(2-a) return 8/3*db/3 end 

Franchement, cependant, si vous êtes pressé par le temps, souvenez-vous que vous pouvez linéariser arccos en 3.14159-1.57079x et simplement:

 function acos(x) return 1.57079-1.57079*x end 

Quoi qu’il en soit, si vous voulez voir une liste de mes équations d’approximation arccos, vous pouvez aller à https://www.desmos.com/calculator/tcaty2sv8l Je sais que mes approximations ne sont pas les meilleures pour certaines choses, mais si vous Si vous faites quelque chose pour lequel mes approximations seraient utiles, veuillez les utiliser, mais essayez de me donner du crédit.

Une autre approche consiste à utiliser des nombres complexes. De la formule de de Moivre ,

x = cos (π / 2 * x) + ⅈ * sin (π / 2 * x)

Soit θ = π / 2 * x. Alors x = 2θ / π, donc

  • sin (θ) = ℑ (^ 2θ / π )
  • cos (θ) = ℜ (^ 2θ / π )

Comment pouvez-vous calculer les puissances de sans péché et cos? Commencez avec une table précalculée pour des puissances de 2:

  • 4 = 1
  • 2 = -1
  • 1 = ⅈ
  • 1/2 = 0,7071067811865476 + 0,7071067811865475 *
  • 1/4 = 0,9238795325112867 + 0,3826834323650898 *
  • 1/8 = 0,9807852804032304 + 0,19509032201612825 * ⅈ
  • 1/16 = 0,9951847266721969 + 0.0980171403295606 *
  • 1/32 = 0,9987954562051724 + 0.049067674327418015 * ⅈ
  • 1/64 = 0,9996988186962042 + 0,024541228522912288 *
  • 1/128 = 0,9999247018391445 + 0.012271538285719925 *
  • 1/256 = 0,9999811752826011 + 0.006135884649154475 *

Pour calculer des valeurs arbitraires de ⅈ x , approcher l’exposant sous forme de fraction binary, puis multiplier les valeurs correspondantes de la table.

Par exemple, pour trouver sin et cos de 72 ° = 0.8π / 2:

0,8 205/256 = ⅈ 0b11001101 = ⅈ 1/2 * 1/4 * 1/32 * 1/64 * 1/256
= 0.3078496400415349 + 0.9514350209690084 * ⅈ

  • sin (72 °) ≈ 0.9514350209690084 (la valeur “exacte” est 0.9510565162951535)
  • cos (72 °) ≈ 0,3078496400415349 (la valeur “exacte” est 0,30901699437494745).

Pour trouver asin et acos, vous pouvez utiliser cette table avec la méthode Bisection:

Par exemple, pour trouver asin (0.6) (le plus petit angle d’un sortingangle 3-4-5):

  • 0 = 1 + 0 *. Le péché est trop petit, alors augmentez x de 1/2.
  • 1/2 = 0,7071067811865476 + 0,7071067811865475 *. Le péché est trop grand, diminuez donc x de 1/4.
  • 1/4 = 0,9238795325112867 + 0,3826834323650898 *. Le péché est trop petit, alors augmentez x de 1/8.
  • 3/8 = 0,8314696123025452 + 0,5555702330196022 *. Le péché est encore trop petit, alors augmentez x de 1/16.
  • 7/16 = 0,773010453362737 + 0,6343932841636455 *. Le péché est trop grand, diminuez donc x de 1/32.
  • 13/32 = 0,8032075314806449 + 0,5956993044924334 *.

Chaque fois que vous augmentez x, multipliez par la puissance correspondante de. Chaque fois que vous diminuez x, divisez par la puissance correspondante de.

Si nous nous arrêtons ici, nous obtenons des acos (0.6) 13/32 * π / 2 = 0.6381360077604268 (La valeur “exacte” est 0.6435011087932844.)

La précision dépend bien entendu du nombre d’itérations. Pour une approximation rapide et approximative, utilisez 10 itérations. Pour une “précision intense”, utilisez 50 à 60 itérations.

Vous pouvez approximer le cosinus inverse avec un polynôme comme suggéré par dan04 , mais un polynôme est une approximation plutôt mauvaise proche de -1 et 1 où la dérivée du cosinus inverse va à l’infini. Lorsque vous augmentez le degré de polynôme, vous obtenez rapidement des rendements décroissants, et il est toujours difficile d’obtenir une bonne approximation autour des points d’extrémité. Une fonction rationnelle (le quotient de deux polynômes) peut donner une bien meilleure approximation dans ce cas.

 acos(x) ≈ π/2 + (ax + bx³) / (1 + cx² + dx⁴) 

 a = -0.939115566365855 b = 0.9217841528914573 c = -1.2845906244690837 d = 0.295624144969963174 

a une erreur absolue maximale de 0,017 radian (0,96 degré) sur l’intervalle (-1, 1). Voici un graphique (le cosinus inverse en noir, l’approximation polynomiale cubique en rouge, la fonction ci-dessus en bleu) pour la comparaison:

Tracé de acos (x) (noir), une approximation polynomiale cubique (rouge) et la fonction ci-dessus (bleu)

Les coefficients ci-dessus ont été choisis pour minimiser l’erreur absolue maximale sur l’ensemble du domaine. Si vous souhaitez autoriser une erreur plus importante aux points d’extrémité, l’erreur sur l’intervalle (-0,98, 0,98) peut être considérablement réduite. Un numérateur de degré 5 et un dénominateur de degré 2 sont à peu près aussi rapides que la fonction ci-dessus, mais légèrement moins précis. Au désortingment des performances, vous pouvez augmenter la précision en utilisant des polynômes de degré supérieur.

Remarque à propos des performances: le calcul des deux polynômes est toujours très économique et vous pouvez utiliser des instructions de multiplication-addition fusionnées. La division n’est pas si mauvaise, car vous pouvez utiliser l’approximation réciproque matérielle et une multiplication. L’erreur dans l’approximation réciproque est négligeable par rapport à l’erreur dans l’approximation acos. Sur un Skylake i7 à 2,6 GHz, cette approximation peut générer environ 8 cosinus inverses tous les 6 cycles avec AVX. (Il s’agit du débit, la latence est supérieure à 6 cycles.)

Voici un site Web formidable avec de nombreuses options: https://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/arcsin/onlyelem.html

Personnellement, je suis allé à l’approximation du quotient de Chebyshev-Pade avec le code suivant:

 double arccos(double x) { const double pi = 3.141592653; return pi / 2 - (.5689111419 - .2644381021*x - .4212611542*(2*x - 1)*(2*x - 1) + .1475622352*(2*x - 1)*(2*x - 1)*(2*x - 1)) / (2.006022274 - 2.343685222*x + .3316406750*(2*x - 1)*(2*x - 1) + .02607135626*(2*x - 1)*(2*x - 1)*(2*x - 1)); } 

Une implémentation rapide d’arccosine, avec une précision d’environ 0,5 degré, peut être basée sur l’ observation que pour x in [0,1], acos (x) ≈ √ (2 * (1-x)). Un facteur d’échelle supplémentaire améliore la précision près de zéro. Le facteur optimal peut être trouvé par une simple recherche binary. Les arguments négatifs sont traités selon acos (-x) = π – acos (x).

 #include  #include  #include  #include  #include  // Approximate acos(a) with relative error < 5.15e-3 // This uses an idea from Robert Harley's posting in comp.arch.arithmetic on 1996/07/12 // https://groups.google.com/forum/#!original/comp.arch.arithmetic/wqCPkCCXqWs/T9qCkHtGE2YJ float fast_acos (float a) { const float PI = 3.14159265f; const float C = 0.10501094f; float r, s, t, u; t = (a < 0) ? (-a) : a; // handle negative arguments u = 1.0f - t; s = sqrtf (u + u); r = C * u * s + s; // or fmaf (C * u, s, s) if FMA support in hardware if (a < 0) r = PI - r; // handle negative arguments return r; } float uint_as_float (uint32_t a) { float r; memcpy (&r, &a, sizeof(r)); return r; } int main (void) { double maxrelerr = 0.0; uint32_t a = 0; do { float x = uint_as_float (a); float r = fast_acos (x); double xx = (double)x; double res = (double)r; double ref = acos (xx); double relerr = (res - ref) / ref; if (fabs (relerr) > maxrelerr) { maxrelerr = fabs (relerr); printf ("xx=% 15.8e res=% 15.8e ref=% 15.8e rel.err=% 15.8e\n", xx, res, ref, relerr); } a++; } while (a); printf ("maximum relative error = %15.8e\n", maxrelerr); return EXIT_SUCCESS; } 

La sortie de l’échafaud de test ci-dessus devrait ressembler à ceci:

 xx= 0.00000000e+000 res= 1.56272149e+000 ref= 1.57079633e+000 rel.err=-5.14060021e-003 xx= 2.98023259e-008 res= 1.56272137e+000 ref= 1.57079630e+000 rel.err=-5.14065723e-003 xx= 8.94069672e-008 res= 1.56272125e+000 ref= 1.57079624e+000 rel.err=-5.14069537e-003 xx=-2.98023259e-008 res= 1.57887137e+000 ref= 1.57079636e+000 rel.err= 5.14071269e-003 xx=-8.94069672e-008 res= 1.57887149e+000 ref= 1.57079642e+000 rel.err= 5.14075044e-003 maximum relative error = 5.14075044e-003