一次计算正弦和余弦
Calculation sine and cosine in one shot
我有一个科学代码,它使用同一参数的正弦和余弦(我基本上需要该参数的复指数)。我想知道是否有可能比分别调用正弦和余弦函数更快。
此外,我只需要大约 0.1% 的精度。那么有什么方法可以找到默认的三角函数并截断幂级数以提高速度吗?
我想到的另一件事是,有什么方法可以执行余数运算,使结果始终为正?在我自己的算法中,我使用了 x=fmod(x,2*pi);
,但如果 x 为负,我需要添加 2pi(较小的域意味着我可以使用较短的幂级数)
编辑:事实证明 LUT 是最好的方法,但是我很高兴我了解了其他近似技术。我还将建议使用显式中点近似。这就是我最终做的:
const int N = 10000;//about 3e-4 error for 1000//3e-5 for 10 000//3e-6 for 100 000
double *cs = new double[N];
double *sn = new double[N];
for(int i =0;i<N;i++){
double A= (i+0.5)*2*pi/N;
cs[i]=cos(A);
sn[i]=sin(A);
}
以下部分近似(中点)sincos(2*pi*(wc2+t[j]*(cotp*t[j]-wc)))
double A=(wc2+t[j]*(cotp*t[j]-wc));
int B =(int)N*(A-floor(A));
re += cs[B]*f[j];
im += sn[B]*f[j];
另一种方法可能是使用切比雪夫分解。您可以使用正交性 属性 来查找系数。针对指数优化,看起来像这样:
double fastsin(double x){
x=x-floor(x/2/pi)*2*pi-pi;//this line can be improved, both inside this
//function and before you input it into the function
double x2 = x*x;
return (((0.00015025063885163012*x2-
0.008034350857376128)*x2+ 0.1659789684145034)*x2-0.9995812174943602)*x;} //7th order chebyshev approx
一种方法是学习如何实现 CORDIC algorithm. It is not difficult and pretty interesting intelectually. This gives you both the cosine and the sine. Wikipedia gives a MATLAB example 应该很容易在 C++ 中适应。
请注意,您可以通过降低参数 n
来提高速度并降低精度。
关于你的第二个问题,已经有人问过 here(在 C 中)。好像没有什么简单的方法。
如果您寻求使用幂级数获得良好(但不高)准确度的快速评估,您应该使用切比雪夫多项式的展开:将系数制成表格(您需要非常少的 0.1% 准确度)并使用这些多项式的递归关系(真的很简单)。
参考文献:
- 制表系数:http://www.ams.org/mcom/1980-34-149/S0025-5718-1980-0551302-5/S0025-5718-1980-0551302-5.pdf
- 切比雪夫展开的计算:https://en.wikipedia.org/wiki/Chebyshev_polynomials
您需要 (a) 在 -pi/2..+pi/2 范围内获取 "reduced" 参数,然后 (b) 处理结果中的符号当参数实际上应该在完整基本区间 -pi..+pi 的 "other" 一半时。这些方面应该不会造成大问题:
- 确定(和"remember"作为整数1或-1)原始角度的符号并继续计算绝对值。
- 使用取模函数缩小到区间0..2PI
- 判断(和"remember"为整数1或-1)是否在"second"一半,如果是则减去pi*3/2,否则减去pi/2.注意:这有效地互换了正弦和余弦(符号除外);在最终评估中考虑到这一点。
这完成了在-pi/2..+pi/2中获取角度的步骤
在使用 Cheb 展开计算正弦和余弦后,应用上面步骤 1 和 3 的 "flags" 以获得值中的正确符号。
只需创建一个查找 table。以下将让您查找 -2PI 和 2PI 之间的任何弧度值的 sin 和 cos。
// LOOK UP TABLE
var LUT_SIN_COS = [];
var N = 14400;
var HALF_N = N >> 1;
var STEP = 4 * Math.PI / N;
var INV_STEP = 1 / STEP;
// BUILD LUT
for(var i=0, r = -2*Math.PI; i < N; i++, r += STEP) {
LUT_SIN_COS[2*i] = Math.sin(r);
LUT_SIN_COS[2*i + 1] = Math.cos(r);
}
您通过以下方式对查找 table 进行索引:
var index = ((r * INV_STEP) + HALF_N) << 1;
var sin = LUT_SIN_COS[index];
var cos = LUT_SIN_COS[index + 1];
这是一个 fiddle,它显示了您可以从不同大小的 LUTS http://jsfiddle.net/77h6tvhj/
中预期的百分比错误
EDIT 这是一个带有 ~benchmark~ 与 float sin 和 cos 对比的 ideone (c++)。 http://ideone.com/SGrFVG 无论 ideone.com 上的基准值是多少,LUT 都要快 5 倍。
您还可以在给定角度和余弦的情况下使用平方根计算正弦。
下例假定角度范围为 0 到 2π:
double c = cos(angle);
double s = sqrt(1.0-c*c);
if(angle>pi)s=-s;
对于单精度浮点数,Microsoft 对正弦使用 11 次多项式逼近,对余弦使用 10 次多项式逼近:XMScalarSinCos。
他们还有更快的版本,XMScalarSinCosEst,使用低阶多项式。
如果您没有使用 Windows,您会在 Boost 许可下找到相同的代码 + 系数 on geometrictools.com。
我有一个科学代码,它使用同一参数的正弦和余弦(我基本上需要该参数的复指数)。我想知道是否有可能比分别调用正弦和余弦函数更快。
此外,我只需要大约 0.1% 的精度。那么有什么方法可以找到默认的三角函数并截断幂级数以提高速度吗?
我想到的另一件事是,有什么方法可以执行余数运算,使结果始终为正?在我自己的算法中,我使用了 x=fmod(x,2*pi);
,但如果 x 为负,我需要添加 2pi(较小的域意味着我可以使用较短的幂级数)
编辑:事实证明 LUT 是最好的方法,但是我很高兴我了解了其他近似技术。我还将建议使用显式中点近似。这就是我最终做的:
const int N = 10000;//about 3e-4 error for 1000//3e-5 for 10 000//3e-6 for 100 000
double *cs = new double[N];
double *sn = new double[N];
for(int i =0;i<N;i++){
double A= (i+0.5)*2*pi/N;
cs[i]=cos(A);
sn[i]=sin(A);
}
以下部分近似(中点)sincos(2*pi*(wc2+t[j]*(cotp*t[j]-wc)))
double A=(wc2+t[j]*(cotp*t[j]-wc));
int B =(int)N*(A-floor(A));
re += cs[B]*f[j];
im += sn[B]*f[j];
另一种方法可能是使用切比雪夫分解。您可以使用正交性 属性 来查找系数。针对指数优化,看起来像这样:
double fastsin(double x){
x=x-floor(x/2/pi)*2*pi-pi;//this line can be improved, both inside this
//function and before you input it into the function
double x2 = x*x;
return (((0.00015025063885163012*x2-
0.008034350857376128)*x2+ 0.1659789684145034)*x2-0.9995812174943602)*x;} //7th order chebyshev approx
一种方法是学习如何实现 CORDIC algorithm. It is not difficult and pretty interesting intelectually. This gives you both the cosine and the sine. Wikipedia gives a MATLAB example 应该很容易在 C++ 中适应。
请注意,您可以通过降低参数 n
来提高速度并降低精度。
关于你的第二个问题,已经有人问过 here(在 C 中)。好像没有什么简单的方法。
如果您寻求使用幂级数获得良好(但不高)准确度的快速评估,您应该使用切比雪夫多项式的展开:将系数制成表格(您需要非常少的 0.1% 准确度)并使用这些多项式的递归关系(真的很简单)。
参考文献:
- 制表系数:http://www.ams.org/mcom/1980-34-149/S0025-5718-1980-0551302-5/S0025-5718-1980-0551302-5.pdf
- 切比雪夫展开的计算:https://en.wikipedia.org/wiki/Chebyshev_polynomials
您需要 (a) 在 -pi/2..+pi/2 范围内获取 "reduced" 参数,然后 (b) 处理结果中的符号当参数实际上应该在完整基本区间 -pi..+pi 的 "other" 一半时。这些方面应该不会造成大问题:
- 确定(和"remember"作为整数1或-1)原始角度的符号并继续计算绝对值。
- 使用取模函数缩小到区间0..2PI
- 判断(和"remember"为整数1或-1)是否在"second"一半,如果是则减去pi*3/2,否则减去pi/2.注意:这有效地互换了正弦和余弦(符号除外);在最终评估中考虑到这一点。
这完成了在-pi/2..+pi/2中获取角度的步骤 在使用 Cheb 展开计算正弦和余弦后,应用上面步骤 1 和 3 的 "flags" 以获得值中的正确符号。
只需创建一个查找 table。以下将让您查找 -2PI 和 2PI 之间的任何弧度值的 sin 和 cos。
// LOOK UP TABLE
var LUT_SIN_COS = [];
var N = 14400;
var HALF_N = N >> 1;
var STEP = 4 * Math.PI / N;
var INV_STEP = 1 / STEP;
// BUILD LUT
for(var i=0, r = -2*Math.PI; i < N; i++, r += STEP) {
LUT_SIN_COS[2*i] = Math.sin(r);
LUT_SIN_COS[2*i + 1] = Math.cos(r);
}
您通过以下方式对查找 table 进行索引:
var index = ((r * INV_STEP) + HALF_N) << 1;
var sin = LUT_SIN_COS[index];
var cos = LUT_SIN_COS[index + 1];
这是一个 fiddle,它显示了您可以从不同大小的 LUTS http://jsfiddle.net/77h6tvhj/
中预期的百分比错误EDIT 这是一个带有 ~benchmark~ 与 float sin 和 cos 对比的 ideone (c++)。 http://ideone.com/SGrFVG 无论 ideone.com 上的基准值是多少,LUT 都要快 5 倍。
您还可以在给定角度和余弦的情况下使用平方根计算正弦。
下例假定角度范围为 0 到 2π:
double c = cos(angle);
double s = sqrt(1.0-c*c);
if(angle>pi)s=-s;
对于单精度浮点数,Microsoft 对正弦使用 11 次多项式逼近,对余弦使用 10 次多项式逼近:XMScalarSinCos。 他们还有更快的版本,XMScalarSinCosEst,使用低阶多项式。
如果您没有使用 Windows,您会在 Boost 许可下找到相同的代码 + 系数 on geometrictools.com。