使用 C 计算 sin 的麦克劳林级数
Calculate maclaurin series for sin using C
我写了一个代码来使用它的麦克劳林级数计算 sin 并且它有效但是当我尝试计算它以获得大的 x 值并尝试通过给出一个大阶数 N(总和的长度)来抵消它时 - 最终它溢出并且没有给我正确的结果。这是代码,我想知道是否有其他方法可以优化它,因此它也适用于大的 x 值(它已经适用于小的 x 值和非常大的 N 值)。
这是代码:
long double calcMaclaurinPolynom(double x, int N){
long double result = 0;
long double atzeretCounter = 2;
int sign = 1;
long double fraction = x;
for (int i = 0; i <= N; i++)
{
result += sign*fraction;
sign = sign*(-1);
fraction = fraction*((x*x) / ((atzeretCounter)*(atzeretCounter + 1)));
atzeretCounter += 2;
}
return result;
}
您可以通过将符号变量合并到 fraction
更新中来避免符号变量,如 (-x*x)
.
使用您的算法,您不会遇到阶乘中整数溢出的问题。
一旦 x*x < (2*k)*(2*k+1)
错误 - 假设精确评估 - 受 abs(fraction)
限制,即系列中下一项的大小。
对于大 x
最大的错误来源是截断 resp。通过取消交替系列的条款放大的浮点错误。对于 k
大约 x/2
,第 k 个项周围的项具有最大的大小,必须被其他大项抵消。
减半平方
在不使用 pi 值的情况下处理大 x
的一种简单方法是使用三角定理,其中
sin(2*x)=2*sin(x)*cos(x)
cos(2*x)=2*cos(x)^2-1=cos(x)^2-sin(x)^2
并首先将 x
减半,同时计算 sin(x/2^n)
和 cos(x/2^n)
的麦克劳林级数,然后使用三角平方(字面平方作为复数 cos(x)+i*sin(x)
)恢复原始参数的值。
cos(x/2^(n-1)) = cos(x/2^n)^2-sin(x/2^n)^2
sin(x/2^(n-1)) = 2*sin(x/2^n)*cos(x/2^n)
然后
cos(x/2^(n-2)) = cos(x/2^(n-1))^2-sin(x/2^(n-1))^2
sin(x/2^(n-2)) = 2*sin(x/2^(n-1))*cos(x/2^(n-1))
等等
见同时计算sin和cos值,然后用
封装
def CosSinForLargerX(x,n):
k=0
while abs(x)>1:
k+=1; x/=2
c,s = getCosSin(x,n)
r2=0
for i in range(k):
s2=s*s; c2=c*c; r2=s2+c2
s = 2*c*s
c = c2-s2
return c/r2,s/r2
主要问题是使用超出收敛范围的级数。
正如 OP 所说 "converted x to radX = (x*PI)/180" 表示 OP 以度数而不是弧度开始,OP 很幸运。找到 my_sin(x)
的第一步是缩小范围。当以度数开始时,减少是准确的。所以在转换为弧度之前减小范围。
long double calcMaclaurinPolynom(double x /* degrees */, int N){
// Reduce to range -360 to 360
// This reduction is exact, no round-off error
x = fmod(x, 360);
// Reduce to range -180 to 180
if (x >= 180) {
x -= 180;
x = -x;
} else if (x <= -180) {
x += 180;
x = -x;
}
// Reduce to range -90 to 90
if (x >= 90) {
x = 180 - x;
} else if (x <= -90) {
x = -180 - x;
}
//now convert to radians.
x = x*PI/180;
// continue with regular code
或者,如果使用 C11,请使用 remquo()
。在 SO 中搜索示例代码。
正如@user3386109 上面评论的那样,不需要"convert back to degrees"。
[编辑]
对于典型的求和级数,首先 对最不重要的项求和可以提高答案的精度。使用 OP 的代码,这可以用
完成
for (int i = N; i >= 0; i--)
或者,不是迭代固定次数,而是循环直到该项对总和没有意义。下面使用递归首先对最低有效项求和。在 -90 到 90 范围内减少范围,迭代次数并不过分。
static double sin_d_helper(double term, double xx, unsigned i) {
if (1.0 + term == 1.0)
return term;
return term - sin_d_helper(term * xx / ((i + 1) * (i + 2)), xx, i + 2);
}
#include <math.h>
double sin_d(double x_degrees) {
// range reduction and d --> r conversion from above
double x_radians = ...
return x_radians * sin_d_helper(1.0, x_radians * x_radians, 1);
}
我写了一个代码来使用它的麦克劳林级数计算 sin 并且它有效但是当我尝试计算它以获得大的 x 值并尝试通过给出一个大阶数 N(总和的长度)来抵消它时 - 最终它溢出并且没有给我正确的结果。这是代码,我想知道是否有其他方法可以优化它,因此它也适用于大的 x 值(它已经适用于小的 x 值和非常大的 N 值)。 这是代码:
long double calcMaclaurinPolynom(double x, int N){
long double result = 0;
long double atzeretCounter = 2;
int sign = 1;
long double fraction = x;
for (int i = 0; i <= N; i++)
{
result += sign*fraction;
sign = sign*(-1);
fraction = fraction*((x*x) / ((atzeretCounter)*(atzeretCounter + 1)));
atzeretCounter += 2;
}
return result;
}
您可以通过将符号变量合并到 fraction
更新中来避免符号变量,如 (-x*x)
.
使用您的算法,您不会遇到阶乘中整数溢出的问题。
一旦 x*x < (2*k)*(2*k+1)
错误 - 假设精确评估 - 受 abs(fraction)
限制,即系列中下一项的大小。
对于大 x
最大的错误来源是截断 resp。通过取消交替系列的条款放大的浮点错误。对于 k
大约 x/2
,第 k 个项周围的项具有最大的大小,必须被其他大项抵消。
减半平方
在不使用 pi 值的情况下处理大 x
的一种简单方法是使用三角定理,其中
sin(2*x)=2*sin(x)*cos(x)
cos(2*x)=2*cos(x)^2-1=cos(x)^2-sin(x)^2
并首先将 x
减半,同时计算 sin(x/2^n)
和 cos(x/2^n)
的麦克劳林级数,然后使用三角平方(字面平方作为复数 cos(x)+i*sin(x)
)恢复原始参数的值。
cos(x/2^(n-1)) = cos(x/2^n)^2-sin(x/2^n)^2
sin(x/2^(n-1)) = 2*sin(x/2^n)*cos(x/2^n)
然后
cos(x/2^(n-2)) = cos(x/2^(n-1))^2-sin(x/2^(n-1))^2
sin(x/2^(n-2)) = 2*sin(x/2^(n-1))*cos(x/2^(n-1))
等等
见同时计算sin和cos值,然后用
封装def CosSinForLargerX(x,n):
k=0
while abs(x)>1:
k+=1; x/=2
c,s = getCosSin(x,n)
r2=0
for i in range(k):
s2=s*s; c2=c*c; r2=s2+c2
s = 2*c*s
c = c2-s2
return c/r2,s/r2
主要问题是使用超出收敛范围的级数。
正如 OP 所说 "converted x to radX = (x*PI)/180" 表示 OP 以度数而不是弧度开始,OP 很幸运。找到 my_sin(x)
的第一步是缩小范围。当以度数开始时,减少是准确的。所以在转换为弧度之前减小范围。
long double calcMaclaurinPolynom(double x /* degrees */, int N){
// Reduce to range -360 to 360
// This reduction is exact, no round-off error
x = fmod(x, 360);
// Reduce to range -180 to 180
if (x >= 180) {
x -= 180;
x = -x;
} else if (x <= -180) {
x += 180;
x = -x;
}
// Reduce to range -90 to 90
if (x >= 90) {
x = 180 - x;
} else if (x <= -90) {
x = -180 - x;
}
//now convert to radians.
x = x*PI/180;
// continue with regular code
或者,如果使用 C11,请使用 remquo()
。在 SO 中搜索示例代码。
正如@user3386109 上面评论的那样,不需要"convert back to degrees"。
[编辑]
对于典型的求和级数,首先 对最不重要的项求和可以提高答案的精度。使用 OP 的代码,这可以用
完成for (int i = N; i >= 0; i--)
或者,不是迭代固定次数,而是循环直到该项对总和没有意义。下面使用递归首先对最低有效项求和。在 -90 到 90 范围内减少范围,迭代次数并不过分。
static double sin_d_helper(double term, double xx, unsigned i) {
if (1.0 + term == 1.0)
return term;
return term - sin_d_helper(term * xx / ((i + 1) * (i + 2)), xx, i + 2);
}
#include <math.h>
double sin_d(double x_degrees) {
// range reduction and d --> r conversion from above
double x_radians = ...
return x_radians * sin_d_helper(1.0, x_radians * x_radians, 1);
}