亿加合和智能车制作

标题: C语言简单实现sin/cos函数 [打印本页]

作者: 724847936 时间: 2011-2-17 13:13
标题: C语言简单实现sin/cos函数
在某些情况下我们需要一些更高效的且近似于标准值的 sin 和cos函数。
有时候我们并需要过高的精度，这时 C语言中自带的三角函数（sinf() 和 cosf() f）计算的精度超出了我们所需要的精度要求，所以其效率很低。我们真正需要的是间于精度和效率的一个折中的方案。众所周知的取近似值的方法是：泰勒级数（和著名的马克劳林级数）
代码是：
x - 1/6 x^3 + 1/120 x^5 - 1/5040 x^7 + ...

我们绘制了下图:

绿线是标准的sin函数，红线是4项泰勒级数展开式。这个近似值的效果看起来还不错，但是如果你仔细观察后会发现

它在pi/2之前的效果还是很好的，但是超过了pi/2后就开始快速偏离标准sin。它在pi处的值比原来的0多了0.075.用这个方法来模拟波动忽动忽停，看起来很呆板，这个效果当然不是我们想要的。
我们可以继续增加泰勒级数项的个数来减小误差，但是这将导致我们的公式非常的冗长。用4项的泰勒级数展开式需要我们进行7次乘法和3次加法来完成。泰勒级数不能同时满足我对精度和效率的要求。
刚刚近似如果能满足sin（pi）=0就好了。从上图我还可以发现另一件事：这个曲线看起来很象抛物线！所以我们来寻找一个尽可能和sin接近的抛物线（公式）。抛物线的范式方程是：A + B x + C x^2.这个公式可以让们控制三个自由度。显然我们需要其满足sine(0) = 0, sine(pi/2) = 1 and sine(pi) = 0. 这样我们就得到了3个等式。A + B 0 + C 0^2 = 0
A + B pi/2 + C (pi/2)^2 = 1
A + B pi + C pi^2 = 0

解得：A = 0, B = 4/pi, C = -4/pi^2.我们的抛物线诞生啦！

貌似这个的误差看起来比泰勒级数还要遭。其实不是的！这种方法的最大误差是0.056.（译者：而且这种近似值没有误差积累）而且这个近似值的绘制出的波动是光滑的，而且只需要3次乘法和一次加法。不过它还不够完美。下图是[-pi, pi] 之间的图像：

显然我们至少需要它在1个完整的周期内都符合我们要求。但是我们可以看出，我们缺少的另一半是原抛物线的一个映射。它的公式是：4/pi x + 4/pi^2 x^2。所以我们可以直接这样写：
Code:
if(x > 0) { y = 4/pi x - 4/pi^2 x^2; } else { y = 4/pi x + 4/pi^2 x^2; }
添加一个条件分支不是一个好的方法。它会让程序渐渐的变慢。但是观察一下我们模拟的和标准的图像是多么的接近啊！观察上面两式子，只是中间的一项正负号不同，我的第一个单纯的想法是可以提取x的正负号来消除分支，即使用：x / abs(x)。除法的代价是非常大的，我们来观察一下现在的公式： 4/pi x - x / abs(x) 4/pi^2 x^2。将除法化简后我们得到：4/pi x - 4/pi^2 x abs(x).所以只需要多一步运算我们就得到了我们需要的sin逼近值！下图是结果

接下来我们要考虑cos。有基础的三角公理可以知道：cos(x) = sin(pi/2 + x).把x多加一个pi/2就可以搞定了？事实上它的某一部分不是我们期望得到的。

我们需要做的就是当x > pi/2时“跳回”。这个可以由减去2 pi来实现。
Code:
x += pi/2;
if(x > pi) // Original x > pi/2 { x -= 2 * pi; // Wrap: cos(x) = cos(x - 2 pi)}
y = sine(x);

又出现了一个分支，我们可以用逻辑“与”来消除它，像是这样：
x -= (x > pi) & (2 * pi);
Code:
x -= (x > pi) & (2 * pi);

注意这并不是c的源代码。但是这个应该可以说明它是怎么样运行的。当x > pi是false 时，逻辑“与”（&）运算后得到的是0，也就是（x-=0）大小没有改变，哈哈完美的等价!我会给读这篇文章的读者留一些关于这个练习。虽然cos 比sin需要多一些运算，但是相比之下貌似也没有更好方法可以让程序更快了。现在我们的最大误差是0.056 ，四项泰勒级数展开式每一次都会有一点点误差。再来看看我们sin函数：

它在pi/2之前的效果还是很好的，但是超过了pi/2后就开始快速偏离标准sin。它在pi处的值比原来的0多了0.075.用这个方法来模拟波动忽动忽停，看起来很呆板，这个效果当然不是我们想要的。
我们可以继续增加泰勒级数项的个数来减小误差，但是这将导致我们的公式非常的冗长。用4项的泰勒级数展开式需要我们进行7次乘法和3次加法来完成。泰勒级数不能同时满足我对精度和效率的要求。
刚刚近似如果能满足sin（pi）=0就好了。从上图我还可以发现另一件事：这个曲线看起来很象抛物线！所以我们来寻找一个尽可能和sin接近的抛物线（公式）。抛物线的范式方程是：A + B x + C x^2.这个公式可以让们控制三个自由度。显然我们需要其满足sine(0) = 0, sine(pi/2) = 1 and sine(pi) = 0. 这样我们就得到了3个等式。A + B 0 + C 0^2 = 0
A + B pi/2 + C (pi/2)^2 = 1
A + B pi + C pi^2 = 0

解得：A = 0, B = 4/pi, C = -4/pi^2.我们的抛物线诞生啦！

貌似这个的误差看起来比泰勒级数还要遭。其实不是的！这种方法的最大误差是0.056.（译者：而且这种近似值没有误差积累）而且这个近似值的绘制出的波动是光滑的，而且只需要3次乘法和一次加法。不过它还不够完美。下图是[-pi, pi] 之间的图像：

显然我们至少需要它在1个完整的周期内都符合我们要求。但是我们可以看出，我们缺少的另一半是原抛物线的一个映射。它的公式是：4/pi x + 4/pi^2 x^2。所以我们可以直接这样写：
Code:
if(x > 0) { y = 4/pi x - 4/pi^2 x^2; } else { y = 4/pi x + 4/pi^2 x^2; }
添加一个条件分支不是一个好的方法。它会让程序渐渐的变慢。但是观察一下我们模拟的和标准的图像是多么的接近啊！观察上面两式子，只是中间的一项正负号不同，我的第一个单纯的想法是可以提取x的正负号来消除分支，即使用：x / abs(x)。除法的代价是非常大的，我们来观察一下现在的公式： 4/pi x - x / abs(x) 4/pi^2 x^2。将除法化简后我们得到：4/pi x - 4/pi^2 x abs(x).所以只需要多一步运算我们就得到了我们需要的sin逼近值！下图是结果

接下来我们要考虑cos。有基础的三角公理可以知道：cos(x) = sin(pi/2 + x).把x多加一个pi/2就可以搞定了？事实上它的某一部分不是我们期望得到的。

我们需要做的就是当x > pi/2时“跳回”。这个可以由减去2 pi来实现。
Code:
x += pi/2;
if(x > pi) // Original x > pi/2 { x -= 2 * pi; // Wrap: cos(x) = cos(x - 2 pi)}
y = sine(x);

又出现了一个分支，我们可以用逻辑“与”来消除它，像是这样：
x -= (x > pi) & (2 * pi);
Code:
x -= (x > pi) & (2 * pi);

注意这并不是c的源代码。但是这个应该可以说明它是怎么样运行的。当x > pi是false 时，逻辑“与”（&）运算后得到的是0，也就是（x-=0）大小没有改变，哈哈完美的等价!我会给读这篇文章的读者留一些关于这个练习。虽然cos 比sin需要多一些运算，但是相比之下貌似也没有更好方法可以让程序更快了。现在我们的最大误差是0.056 ，四项泰勒级数展开式每一次都会有一点点误差。再来看看我们sin函数：



现在是不是不能继续提升精准度了呢？当前的版本已经可以满足大多度sin函数的应用了。但是对一些要求更高一些的程序现在做的还够。

仔细观察图像，你会注意到我们的近似值总是比真实值大，当然除了0，pi/2 和 pi。所以我们要做的就是在不改变这些点（0，pi/2 和 pi）的情况下，将函数再“按下去”一些。解决方法是利用抛物线的平方。看起来就像这样：

注意它保持着原来那些关键点，不同的是它比真实的sin函数值更低了。所以我们可以用一个加权的平均值来使两个函数更接近。

Code:
Q (4/pi x - 4/pi^2 x^2) + P (4/pi x - 4/pi^2 x^2)^2

利用Q + P = 1. 你可以灵活的控制绝对误差或相对误差。别急我来告诉你取不同的极限结果时Q，P的值。绝对误差的最佳权值是：Q = 0.775, P = 0.225 ；相对误差的最佳权值是：Q = 0.782，P = 0.218 。让我们来看一下结果的图像。

红线呢？它几乎被绿线完全覆盖了，这足以证明我们的近似十分完美。最大误差是0.001，50倍的提升！这个公式看起来很长，但是括号里面的公式最终得到的值是相同的，也就是说括号里的只需要被计算一次。事实上在原来的基础上只是增加了额外的2次乘法和2次加法就可以得到现在的结果。

先别高兴的太早，我们还要“制造”一个负号出来。我们需要增加一个abs()运算。最终的c代码是：
Code:
float sine(float x)
{
const float B = 4/pi;
const float C = -4/(pi*pi);
float y = B * x + C * x * abs(x);
#ifdef EXTRA_PRECISION // const float Q = 0.775;
const float P = 0.225;
y = P * (y * abs(y) - y) + y;
// Q * y + P * y * abs(y)
#endif }
所以我们仅仅是需要多加5次乘法和3次加法就可以完成了。如果我们忽略abs()这个仍然是比4项泰勒级数展开式快，更精准！Cos只需要相应的变换一下x就可以了。
(译者注：后面是汇编程序，不翻译了)
part2

我选取了最小误差的情况，用as3运行后发现提升了14倍，而且仍然是非常精准。不过你必须直接使用它，不能把它放到一个函数中，因为每调用一次额外的函数调用会削减执行效率，最终你会得到一个比Math.sin() 和 Math.cos()效率更差的结果。还有这里会用到的三角定理：
cos(x) = sin(x + pi/2)
cos(x - pi/2) = sin(x)
下载: fastTrig.as.
可以清楚到对比结果，现在你可以用这个替换Math.sin() 和 Math.cos()了
//always wrap input angle to -PI..PI
if (x< -3.14159265)
x+= 6.28318531;
else
if (x>  3.14159265)
x-= 6.28318531;

//compute sine
if (x< 0)
sin= 1.27323954 * x+ .405284735 * x* x;
else
sin= 1.27323954 * x- 0.405284735 * x* x;

//compute cosine: sin(x + PI/2) = cos(x)
x+= 1.57079632;
if (x>  3.14159265)
x-= 6.28318531;

if (x< 0)
cos= 1.27323954 * x+ 0.405284735 * x* x
else
cos= 1.27323954 * x- 0.405284735 * x* x;
}
High precision sine/cosine (~8x faster)
//always wrap input angle to -PI..PI
if (x< -3.14159265)
x+= 6.28318531;
else
if (x>  3.14159265)
x-= 6.28318531;

//compute sine
if (x< 0)
{
sin= 1.27323954 * x+ .405284735 * x* x;

if (sin< 0)
      sin= .225 * (sin*-sin- sin) + sin;
else
      sin= .225 * (sin* sin- sin) + sin;
}
else
{
sin= 1.27323954 * x- 0.405284735 * x* x;

if (sin< 0)
      sin= .225 * (sin*-sin- sin) + sin;
else
      sin= .225 * (sin* sin- sin) + sin;
}

//compute cosine: sin(x + PI/2) = cos(x)
x+= 1.57079632;
if (x>  3.14159265)
x-= 6.28318531;

if (x< 0)
{
cos= 1.27323954 * x+ 0.405284735 * x* x;

if (cos< 0)
      cos= .225 * (cos*-cos- cos) + cos;
else
      cos= .225 * (cos* cos- cos) + cos;
}
else
{
cos= 1.27323954 * x- 0.405284735 * x* x;

if (cos< 0)
      cos= .225 * (cos*-cos- cos) + cos;
else
      cos= .225 * (cos* cos- cos) + cos;
}

作者: turf456 时间: 2011-2-17 13:20
不错，可惜没图。

作者: rsdlite 时间: 2011-2-17 13:31
图？

作者: zelight 时间: 2012-1-15 14:01
nice

作者: Stanger 时间: 2013-8-27 20:21
大哥，图呢？？？？

作者: Quixote 时间: 2015-1-13 01:00
这个方法实际是很慢的

作者: susutoptom 时间: 2017-11-7 14:33
好歹上个图啊

欢迎光临亿加合和智能车制作 (http://www.znczz.com/)