Fonction Math.Sin (double) pour GPU

Préface


J'avais besoin de calculer l'arc avec une précision accrue sur le processeur de la carte vidéo en temps réel.

L'auteur ne s'est pas fixé pour objectif de dépasser la fonction standard System.Math.Sin () (C #) et ne l'a pas atteint.

Le résultat du travail et mon choix (pour ceux qui ne veulent pas lire):

Sin_3 (rad)
using System; class Math_d { const double PI025 = Math.PI / 4; /// <summary> 2^17 = 131072 (1 ),   10000 ( ),  2^21 = 22097152 (16 )   +-1 ( ) (  ) </summary> const int length_mem = 22097152; const int length_mem_M1 = length_mem - 1; /// <summary>    sin,    . </summary> static double[] mem_sin; /// <summary>    cos,    . </summary> static double[] mem_cos; /// <summary>  ,   sin,    . </summary> public static void Initialise() { Ini_Mem_Sin(); Ini_Mem_Cos(); } /// <summary>       Cos,   . </summary> /// <param name="rad"></param> public static double Sin_3(double rad) { double rad_025; int i; //    //if (rad < 0) { rad = -rad + Math.PI; } i = (int)(rad / PI025); //   rad_025 = rad - PI025 * i; //     ( ) i = i & 7; //      8 //    switch (i) { case 0: return Sin_Lerp(rad_025); case 1: return Cos_Lerp(PI025 - rad_025); case 2: return Cos_Lerp(rad_025); case 3: return Sin_Lerp(PI025 - rad_025); case 4: return -Sin_Lerp(rad_025); case 5: return -Cos_Lerp(PI025 - rad_025); case 6: return -Cos_Lerp(rad_025); case 7: return -Sin_Lerp(PI025 - rad_025); } return 0; } /// <summary>   sin    </summary> static void Ini_Mem_Sin() { double rad; mem_sin = new double[length_mem]; for (int i = 0; i < length_mem; i++) { rad = (i * PI025) / length_mem_M1; mem_sin[i] = Math.Sin(rad); } } /// <summary>   cos    </summary> static void Ini_Mem_Cos() { double rad; mem_cos = new double[length_mem]; for (int i = 0; i < length_mem; i++) { rad = (i * PI025) / length_mem_M1; mem_cos[i] = Math.Cos(rad); } } /// <summary>      sin  0  pi/4. </summary> /// <param name="rad">     0  pi/4. </param> static double Sin_Lerp(double rad) { int i_0; int i_1; double i_0d; double percent; double a; double b; double s; percent = rad / PI025; i_0d = percent * length_mem_M1; i_0 = (int)i_0d; i_1 = i_0 + 1; a = mem_sin[i_0]; b = mem_sin[i_1]; s = i_0d - i_0; return Lerp(a, b, s); } /// <summary>      cos  0  pi/4. </summary> /// <param name="rad">     0  pi/4. </param> static double Cos_Lerp(double rad) { int i_0; int i_1; double i_0d; double percent; double a; double b; double s; percent = rad / PI025; i_0d = percent * length_mem_M1; i_0 = (int)i_0d; i_1 = i_0 + 1; a = mem_cos[i_0]; b = mem_cos[i_1]; s = i_0d - i_0; return Lerp(a, b, s); } /// <summary>      . (return a + s * (b - a)) </summary> /// <param name="a">  . </param> /// <param name="b">  . </param> /// <param name="s">  . 0 = a, 1 = b, 0.5 =   a  b. </param> public static double Lerp(double a, double b, double s) { return a + s * (b - a); } } 


Motifs de publication


  • Il n'y a pas de fonction Sin standard pour double dans le langage HLSL (mais ce n'est pas prĂ©cis)
  • Il y a peu d'informations disponibles sur Internet Ă  ce sujet.

Approches envisagées



Les paramètres analysés


  • PrĂ©cision avec Math.Sin
  • Vitesse par rapport Ă  Math.Sin

En plus de l'analyse, nous améliorerons leurs performances.

Taylor se classe


Avantages:

  • PrĂ©cision maximale
    Cette fonction, utilisée pour calculer la valeur Sin, peut être utilisée pour calculer la valeur Sin infiniment précise . Plus il subit d'itérations, plus la valeur est obtenue en sortie avec précision (dans une hypothèse). Dans la pratique de la programmation, il convient de considérer les erreurs d'arrondi des calculs en fonction des types de paramètres utilisés (double, flottant, décimal et autres).
  • Calcule n'importe quel angle
    Vous pouvez entrer n'importe quelle valeur comme argument de la fonction, il n'est donc pas nécessaire de surveiller les paramètres d'entrée.
  • L'indĂ©pendance
    Il ne nécessite pas de calculs préliminaires (comme les fonctions discutées ci-dessous), et c'est souvent la base sur laquelle des fonctions plus rapides sont assemblées.

Inconvénients:

  • Très faible vitesse (4-10%)
    Il faut beaucoup d'itérations pour rapprocher la précision de la précision de Math.Sin, ce qui fait qu'elle fonctionne 25 fois plus lentement que la fonction standard.
  • Plus l'angle est grand, plus la prĂ©cision est faible
    Plus l'angle entré dans la fonction est grand, plus il faut d'itérations pour atteindre la même précision que Math.Sin.

Aspect original (vitesse: 4%):

La fonction standard implique le calcul des factorielles, ainsi que l'augmentation d'une puissance à chaque itération.

image

Modifié (vitesse: 10%):

Le calcul de certains degrés peut être réduit dans un cycle (a * = aa;), et d'autres factorielles peuvent être pré-calculées et placées dans un tableau, tout en changeant les signes (+, -, +, ...) ne peuvent pas être augmentés à une puissance et leur calcul peut également être réduit en utilisant les valeurs précédentes.

Le résultat est le code suivant:

Péché (rad, étapes)
  // <summary>  ,    Fact </summary> static double[] fact; /// <summary>            . ///   rad,   . ///  ( Math): 4% (fps)  steps = 17 </summary> /// <param name="rad">   .      pi/4. </param> /// <param name="steps">  :  ,   .  pi/4   E-15  8. </param> public static double Sin(double rad, int steps) { double ret; double a; //,     double aa; // *  int i_f; //  int sign; // (  -  +,     = +) ret = 0; sign = -1; aa = rad * rad; a = rad; i_f = 1; //      for (int n = 0; n < steps; n++) { sign *= -1; ret += sign * a / Fact(i_f); a *= aa; i_f += 2; } return ret; } /// <summary>   (n!).  n > fact.Length,  -1. </summary> /// <param name="n"> ,     . </param> public static double Fact(int n) { if (n >= 0 && n < fact.Length) { return fact[n]; } else { Debug.Log("    . n = " + n + ",   = " + fact.Length); return -1; } } /// <summary>  . </summary> static void Init_Fact() { int steps; steps = 46; fact = new double[steps]; fact[0] = 1; for (int n = 1; n < steps; n++) { fact[n] = fact[n - 1] * n; } } 


Vue supérieure (vitesse: 19%):

Nous savons que plus l'angle est petit, moins il faut d'itérations. Le plus petit angle dont nous avons besoin est = 0,25 * PI, c'est-à-dire 45 degrés. En considérant Sin et Cos dans la région de 45 degrés, nous pouvons obtenir toutes les valeurs de -1 à 1 pour Sin (dans la région de 2 * PI). Pour ce faire, nous divisons le cercle (2 * PI) en 8 parties et pour chaque partie nous indiquons notre propre méthode de calcul du sinus. De plus, pour accélérer le calcul, nous refuserons d'utiliser la fonction d'obtention du reste (%) (pour obtenir la position de l'angle à l'intérieur de la zone de 45 degrés):

Sin_2 (rad)
  // <summary>  ,    Fact </summary> static double[] fact; /// <summary>   Sin </summary> /// <param name="rad"></param> public static double Sin_2(double rad) { double rad_025; int i; //rad = rad % PI2; //% -   .  , fps   90  150 (  100 000   ) //rad_025 = rad % PI025; i = (int)(rad / PI025); rad_025 = rad - PI025 * i; i = i & 7; //     8  //    switch (i) { case 0: return Sin(rad_025, 8); case 1: return Cos(PI025 - rad_025, 8); case 2: return Cos(rad_025, 8); case 3: return Sin(PI025 - rad_025, 8); case 4: return -Sin(rad_025, 8); case 5: return -Cos(PI025 - rad_025, 8); case 6: return -Cos(rad_025, 8); case 7: return -Sin(PI025 - rad_025, 8); } return 0; } /// <summary>            . ///   rad,   . ///  ( Math): 10% (fps)  steps = 17 </summary> /// <param name="rad">   .      pi/4. </param> /// <param name="steps">  :  ,   .  pi/4   E-15  8. </param> public static double Sin(double rad, int steps) { double ret; double a; //,     double aa; // *  int i_f; //  int sign; // (  -  +,     = +) ret = 0; sign = -1; aa = rad * rad; a = rad; i_f = 1; //      for (int n = 0; n < steps; n++) { sign *= -1; ret += sign * a / Fact(i_f); a *= aa; i_f += 2; } return ret; } /// <summary>            . ///   rad,   . ///  ( Math): 10% (fps), 26% (test)  steps = 17 </summary> /// <param name="rad">   .      pi/4. </param> /// <param name="steps">  :  ,   .  pi/4   E-15  8. </param> public static double Cos(double rad, int steps) { double ret; double a; double aa; // *  int i_f; //  int sign; // (  -  +,     = +) ret = 0; sign = -1; aa = rad * rad; a = 1; i_f = 0; //      for (int n = 0; n < steps; n++) { sign *= -1; ret += sign * a / Fact(i_f); a *= aa; i_f += 2; } return ret; } /// <summary>   (n!).  n > fact.Length,  -1. </summary> /// <param name="n"> ,     . </param> public static double Fact(int n) { if (n >= 0 && n < fact.Length) { return fact[n]; } else { Debug.Log("    . n = " + n + ",   = " + fact.Length); return -1; } } /// <summary>  . </summary> static void Init_Fact() { int steps; steps = 46; fact = new double[steps]; fact[0] = 1; for (int n = 1; n < steps; n++) { fact[n] = fact[n - 1] * n; } } 


PolynĂ´mes


Je suis tombé sur cette méthode sur Internet, l'auteur avait besoin d'une fonction de recherche rapide Sin pour le double de faible précision (erreur <0,000 001) sans utiliser de bibliothèques de valeurs pré-calculées.

Avantages:

  • Haute vitesse (9-84%)
    Initialement, le polynôme lancé sans modifications a montré une vitesse de 9% du Math.Sin d'origine, ce qui est 10 fois plus lent. Grâce à de petits changements, la vitesse monte brusquement à 84%, ce qui n'est pas mal si vous fermez les yeux à la précision.
  • Aucun calcul prĂ©liminaire supplĂ©mentaire et aucune mĂ©moire requise
    Si ci-dessus et ci-dessous ci-dessous, nous devons composer des tableaux de variables pour accélérer les calculs, alors ici tous les coefficients clés ont été aimablement calculés et placés dans la formule par l'auteur lui-même sous la forme de constantes.
  • PrĂ©cision supĂ©rieure Ă  Mathf.Sin (float)
    A titre de comparaison:

    0,84147 1 - Mathf.Sin (1) (moteur Unity);
    0,841470984807897 - Math.Sin (1) (fonction C # standard);
    0,8414709 56802368 - sin (1) (GPU, langage hlsl);
    0,84147 1184637935 - Sin_0 (1) .

Inconvénients:

  • Pas universel
    Vous ne pouvez pas ajuster la précision manuellement car on ne sait pas quels outils l'auteur a utilisés pour calculer ce polynôme.
  • Pourquoi?
    Pourquoi l'auteur a-t-il eu besoin d'une fonction qui ne nécessite aucun tableau et qui a une précision aussi faible (par rapport au double)?

Vue originale:

Sin_1 (x)
 /// <summary>  ( Math): 9% (fps)</summary> /// <param name="x">     -2*Pi  2*Pi </param> public static double Sin_1(double x) { return 0.9999997192673006 * x - 0.1666657564532464 * Math.Pow(x, 3) + 0.008332803647181511 * Math.Pow(x, 5) - 0.00019830197237204295 * Math.Pow(x, 7) + 2.7444305061093514e-6 * Math.Pow(x, 9) - 2.442176561869478e-8 * Math.Pow(x, 11) + 1.407555708887347e-10 * Math.Pow(x, 13) - 4.240664814288337e-13 * Math.Pow(x, 15); } 


Vue supérieure:

Sin_0 (rad)
 /// <summary>  ( Math): 83% (fps)</summary> /// <param name="rad">     -2*Pi  2*Pi </param> public static double Sin_0(double rad) { double x; double xx; double ret; xx = rad * rad; x = rad; //1 ret = 0.9999997192673006 * x; x *= xx; //3 ret -= 0.1666657564532464 * x; x *= xx; //5 ret += 0.008332803647181511 * x; x *= xx; //7 ret -= 0.00019830197237204295 * x; x *= xx; //9 ret += 2.7444305061093514e-6 * x; x *= xx; //11 ret -= 2.442176561869478e-8 * x; x *= xx; //13 ret += 1.407555708887347e-10 * x; x *= xx; //15 ret -= 4.240664814288337e-13 * x; return ret; } 


Interpolation linéaire


Cette méthode est basée sur une interpolation linéaire entre les résultats de deux enregistrements dans un tableau.
Les entrées sont divisées en mem_sin et mem_cos, elles contiennent les résultats pré-calculés de la fonction standard Math.Sin et Math.Cos sur un segment de paramètres d'entrée de 0 à 0,25 * PI.

Le principe des manipulations avec un angle de 0 à 45 degrés ne diffère pas de la version améliorée de la série Taylor, mais en même temps une fonction est appelée qui trouve - entre deux enregistrements il y a un angle - et trouve une valeur entre eux.

Avantages:

  • Haute vitesse (65%)
    En raison de la simplicité de l'algorithme d'interpolation, la vitesse atteint 65% de la vitesse de Math.Sin. Je pense que la vitesse> 33% est satisfaisante.
  • PrĂ©cision maximale
    Un exemple d'un cas rare de rejet:
    0,2555835595715180 - Math.Sin;
    0,25558355957151 79 - Sin_3 .
  • Pied rapide
    J'adore cette fonction car elle est née à l'agonie, écrite par moi et a dépassé les exigences: vitesse> 33%, précision supérieure à 1e-14. Je vais lui donner un nom fier - Vēlōx Pes.

Inconvénients:

  • Il faut une place en mĂ©moire
    Pour travailler, vous devez d'abord calculer deux tableaux: pour sin et pour cos; chaque tableau pèse ~ 16 Mo (16 * 2 = 32 Mo)

Vue originale:

Sin_3 (rad)
 class Math_d { const double PI025 = Math.PI / 4; /// <summary> 2^17 = 131072 (1 ),   10000 ( ),  2^21 = 22097152 (16 )   +-1 ( ) (  ) </summary> const int length_mem = 22097152; const int length_mem_M1 = length_mem - 1; /// <summary>    sin,    . </summary> static double[] mem_sin; /// <summary>    cos,    . </summary> static double[] mem_cos; /// <summary>  ,   sin,    . </summary> public static void Initialise() { Ini_Mem_Sin(); Ini_Mem_Cos(); } /// <summary>       Cos,   . </summary> /// <param name="rad"></param> public static double Sin_3(double rad) { double rad_025; int i; //    //if (rad < 0) { rad = -rad + Math.PI; } i = (int)(rad / PI025); //   rad_025 = rad - PI025 * i; //     ( ) i = i & 7; //      8 //    switch (i) { case 0: return Sin_Lerp(rad_025); case 1: return Cos_Lerp(PI025 - rad_025); case 2: return Cos_Lerp(rad_025); case 3: return Sin_Lerp(PI025 - rad_025); case 4: return -Sin_Lerp(rad_025); case 5: return -Cos_Lerp(PI025 - rad_025); case 6: return -Cos_Lerp(rad_025); case 7: return -Sin_Lerp(PI025 - rad_025); } return 0; } /// <summary>   sin    </summary> static void Ini_Mem_Sin() { double rad; mem_sin = new double[length_mem]; for (int i = 0; i < length_mem; i++) { rad = (i * PI025) / length_mem_M1; mem_sin[i] = Math.Sin(rad); } } /// <summary>   cos    </summary> static void Ini_Mem_Cos() { double rad; mem_cos = new double[length_mem]; for (int i = 0; i < length_mem; i++) { rad = (i * PI025) / length_mem_M1; mem_cos[i] = Math.Cos(rad); } } /// <summary>      sin  0  pi/4. </summary> /// <param name="rad">     0  pi/4. </param> static double Sin_Lerp(double rad) { int i_0; int i_1; double i_0d; double percent; double a; double b; double s; percent = rad / PI025; i_0d = percent * length_mem_M1; i_0 = (int)i_0d; i_1 = i_0 + 1; a = mem_sin[i_0]; b = mem_sin[i_1]; s = i_0d - i_0; return Lerp(a, b, s); } /// <summary>      cos  0  pi/4. </summary> /// <param name="rad">     0  pi/4. </param> static double Cos_Lerp(double rad) { int i_0; int i_1; double i_0d; double percent; double a; double b; double s; percent = rad / PI025; i_0d = percent * length_mem_M1; i_0 = (int)i_0d; i_1 = i_0 + 1; a = mem_cos[i_0]; b = mem_cos[i_1]; s = i_0d - i_0; return Lerp(a, b, s); } /// <summary>      . (return a + s * (b - a)) </summary> /// <param name="a">  . </param> /// <param name="b">  . </param> /// <param name="s">  . 0 = a, 1 = b, 0.5 =   a  b. </param> public static double Lerp(double a, double b, double s) { return a + s * (b - a); } } 



UPD: correction d'une erreur dans la détermination de l'index dans Sin_Lerp (), Cos_Lerp (), Ini_Mem_Sin () et Ini_Mem_Cos ().

Source: https://habr.com/ru/post/fr426355/


All Articles