适用于GPU的Math.Sin(double)函数

前言


我需要实时地在视频卡处理器上计算出弧度。

作者并没有设定超越标准函数System.Math.Sin()(C#)的目标,并且没有实现。

工作的结果和我的选择(对于那些不想阅读的人):

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


出版理由


  • HLSL语言中没有用于双精度的标准Sin函数(但这并不准确)
  • Internet上关于此主题的信息很少。

考虑的方法



分析参数


  • Math.Sin的准确性
  • 相对于Math.Sin的速度

除了分析之外,我们还将改善其性能。

泰勒排名


优点:

  • 最高精度
    该函数用于计算Sin值,可用于计算无限精确的 Sin 。 它进行的迭代次数越多,在输出(假设)中获得的值越准确。 在编程实践中,值得考虑的是根据所使用的参数类型(双精度,浮点型,十进制等)的舍入误差。
  • 计算任何角度
    您可以输入任何值作为函数的参数,因此无需监视输入参数。
  • 独立性
    它不需要进行初步计算(如下文所述的功能),并且通常是组装更快的功能的基础。

缺点:

  • 极低的速度(4-10%)
    要使精度接近Math.Sin的精度,需要进行大量的迭代,结果它的工作速度比标准函数慢25倍。
  • 角度越大,精度越低
    在函数中输入的角度越大,要获得与Math.Sin相同的精度,需要进行更多的迭代。

原始外观(速度:4%):

标准函数包括计算阶乘,以及每次迭代提高幂。

图片

修改(速度:10%):

可以在一个周期内减少某些程度的计算(a * = aa;),并且可以预先计算其他阶乘并将其放入数组中,同时不能将符号(+,-,+,...)转换为幂,也可以减少其计算使用以前的值。

结果是以下代码:

罪恶(rad,步骤)
  // <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; } } 


全景(速度:19%):

我们知道角度越小,需要的迭代越少。 我们需要的最小角度为0.25 * PI,即 45度。 考虑到Sin和Cos在45度的范围内,我们可以获得Sin的所有值(从-1到1)(在2 * PI范围内)。 为此,我们将圆(2 * PI)分为8个部分,并为每个部分指示了自己的正弦计算方法。 此外,为了加快计算速度,我们将拒绝使用获得余数(%)的功能(以获取45度区域内角度的位置):

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


多项式


我在Internet上遇到了这种方法,作者需要一种快速的Sin搜索功能,以实现两倍的低准确性(错误<0.000 001),而无需使用预先计算的值的库。

优点:

  • 高速(9-84%)
    最初,未更改的抛出多项式显示的速度是原始Math.Sin的9%,慢了10倍。 由于进行了微小的更改,速度急剧提高到了84%,如果您闭上双眼以保持准确性,这还不错。
  • 无需额外的初步计算和存储
    如果上方和下方都需要组成变量数组以加快计算速度,则此处所有关键系数均由亲切计算,并由作者本人以常量的形式放入公式中。
  • 精度高于Mathf.Sin(浮点)
    为了比较:

    0.84147 1 -Mathf.Sin(1)(统一引擎);
    0.841470984807897-数学正弦(1)(标准C#函数);
    0.8414709 56802368 -sin(1)(GPU,hlsl语言);
    0.84147 1184637935 - Sin_0(1)

缺点:

  • 不普遍
    您无法手动调整精度,因为尚不知道作者使用什么工具来计算此多项式。
  • 怎么了
    为什么作者需要一个不需要任何数组且精度如此之低(相对于两倍)的函数?

原始视图:

罪_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); } 


高级视图:

Sin_0(弧度)
 /// <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; } 


线性插补


此方法基于数​​组中两个记录的结果之间的线性插值。
这些条目分为mem_sin和mem_cos,它们包含标准函数Math.Sin和Math.Cos在输入参数从0到0.25 * PI的段上的预先计算的结果。

角度在0到45度之间的操作原理与泰勒级数的改进版本没有区别,但是同时调用了一个函数,该函数查找-在两个记录之间存在一个角度-并在它们之间找到一个值。

优点:

  • 高速(65%)
    由于插值算法的简单性,其速度达到Math.Sin速度的65%。 我认为速度> 33%是令人满意的。
  • 最高精度
    罕见的拒绝案例:
    0.255835595715180-Math.Sin;
    0.2558355957151 79 - Sin_3
  • 快脚
    我喜欢这个功能,因为它是我自己写的,非常痛苦,超出了要求:速度> 33%,精度在1e-14以上。 我会给她起一个骄傲的名字-VēlōxPes。

缺点:

  • 它需要在内存中的位置
    要工作,必须首先计算两个数组:sin和cos; 每个阵列的重量约为16mb(16 * 2 = 32mb)

原始视图:

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:修复了确定Sin_Lerp(),Cos_Lerp(),Ini_Mem_Sin()和Ini_Mem_Cos()中的索引时的错误。

Source: https://habr.com/ru/post/zh-CN426355/


All Articles