我必须计算以下内容:
float2 y = CONSTANT;
for (int i = 0; i < totalN; i++)
h[i] = cos(y*i);
totalN 是一个很大的数字,所以我想以更有效的方式进行此操作。有什么办法可以改善这一点吗?我怀疑有,因为毕竟,对于 n=1..N,我们知道 cos(n) 的结果是什么,所以也许有一些定理可以让我以更快的方式计算它。我真的很感激任何提示。
提前致谢,
费德里科
我必须计算以下内容:
float2 y = CONSTANT;
for (int i = 0; i < totalN; i++)
h[i] = cos(y*i);
totalN 是一个很大的数字,所以我想以更有效的方式进行此操作。有什么办法可以改善这一点吗?我怀疑有,因为毕竟,对于 n=1..N,我们知道 cos(n) 的结果是什么,所以也许有一些定理可以让我以更快的方式计算它。我真的很感激任何提示。
提前致谢,
费德里科
使用最美丽的数学公式之一,欧拉公式
exp(i*x) = cos(x) + i*sin(x)
,
替换x := n * phi
:
cos(n*phi) = Re( exp(i*n*phi) )
sin(n*phi) = Im( exp(i*n*phi) )
exp(i*n*phi) = exp(i*phi) ^ n
幂^n
是n
重复的乘法。因此,您可以通过从 开始重复复数乘法cos(n*phi)
同时计算 和。sin(n*phi)
exp(i*phi)
(1+i*0)
代码示例:
Python:
from math import *
DEG2RAD = pi/180.0 # conversion factor degrees --> radians
phi = 10*DEG2RAD # constant e.g. 10 degrees
c = cos(phi)+1j*sin(phi) # = exp(1j*phi)
h=1+0j
for i in range(1,10):
h = h*c
print "%d %8.3f"%(i,h.real)
或 C:
#include <stdio.h>
#include <math.h>
// numer of values to calculate:
#define N 10
// conversion factor degrees --> radians:
#define DEG2RAD (3.14159265/180.0)
// e.g. constant is 10 degrees:
#define PHI (10*DEG2RAD)
typedef struct
{
double re,im;
} complex_t;
int main(int argc, char **argv)
{
complex_t c;
complex_t h[N];
int index;
c.re=cos(PHI);
c.im=sin(PHI);
h[0].re=1.0;
h[0].im=0.0;
for(index=1; index<N; index++)
{
// complex multiplication h[index] = h[index-1] * c;
h[index].re=h[index-1].re*c.re - h[index-1].im*c.im;
h[index].im=h[index-1].re*c.im + h[index-1].im*c.re;
printf("%d: %8.3f\n",index,h[index].re);
}
}
我不确定您愿意做出什么样的准确性与性能折衷,但在这些链接上有各种正弦近似技术的广泛讨论:
正弦曲线的乐趣 - http://www.audiomulch.com/~rossb/code/sinusoids/
快速准确的正弦/余弦 - http://www.devmaster.net/forums/showthread.php?t=5784
编辑(我认为这是“Fun with Sinusoids”页面上损坏的“Don Cross”链接):
优化 Trig 计算 - http://groovit.disjunkt.com/analog/time-domain/fasttrig.html
也许最简单的公式是
cos(n+y) = 2cos(n)cos(y) - cos(ny)。
如果您预先计算常数 2*cos(y),则每个值 cos(n+y) 都可以通过一次乘法和一次减法从前 2 个值中计算出来。即,在伪代码中
h[0] = 1.0
h[1] = cos(y)
m = 2*h[1]
for (int i = 2; i < totalN; ++i)
h[i] = m*h[i-1] - h[i-2]
这是一种方法,但它使用了一点内存来进行罪恶。它使用三角恒等式:
cos(a + b) = cos(a)cos(b)-sin(a)sin(b)
sin(a + b) = sin(a)cos(b)+cos(a)sin(b)
然后是代码:
h[0] = 1.0;
double g1 = sin(y);
double glast = g1;
h[1] = cos(y);
for (int i = 2; i < totalN; i++){
h[i] = h[i-1]*h[1]-glast*g1;
glast = glast*h[1]+h[i-1]*g1;
}
如果我没有犯任何错误,那么应该这样做。当然可能存在舍入问题,因此请注意这一点。我在 Python 中实现了这个,它非常准确。
这里有一些很好的答案,但它们都是递归的。使用浮点运算时,递归计算不适用于余弦函数;你总是会得到快速复合的舍入错误。
考虑计算 y = 45 度,总计 N 10 000。最终结果不会是 1。
为了解决 Kirk 的担忧:所有基于 cos 和 sin 递归的解决方案都归结为计算
x(k) = R x(k - 1),
其中 R 是旋转 y 的矩阵,x(0) 是单位向量 (1, 0)。如果 k - 1 的真实结果是 x'(k - 1),而 k 的真实结果是 x'(k),则错误来自 e(k - 1) = x(k - 1) - x' (k - 1) 到 e(k) = R x(k - 1) - R x'(k - 1) = R e(k - 1) 通过线性。由于 R 是所谓的正交矩阵,因此 R e(k - 1) 与 e(k - 1) 具有相同的范数,并且误差增长非常缓慢。(它完全增长的原因是由于四舍五入;R 的计算机表示通常几乎是正交的,但不是完全正交的,因此有必要根据准确性不时使用三角操作重新启动递归需要。这仍然比使用三角运算来计算每个值要快得多。)
您可以使用复数来执行此操作。
如果定义 x = sin(y) + i cos(y),则 cos(y*i) 将是 x^i 的实部。
您可以迭代地计算所有 i。复数乘法是 2 次乘法加上 2 次加法。
知道 cos(n) 并没有帮助——你的数学库已经为你做了这些琐碎的事情。
如果您预先计算cos ( y )和sin(y),并在此过程中跟踪 cos(i y) 和 sin(i*y)。但是,它可能会导致一些精度损失 - 您必须进行检查。
您需要得到的 cos(x) 有多准确?如果你能忍受一些,你可以创建一个查找表,以 2*PI/N 间隔对单位圆进行采样,然后在两个相邻点之间进行插值。将选择 N 以达到某种所需的准确度水平。
我不知道插值实际上是否比计算余弦更便宜。由于它通常在现代 CPU 中以微码形式完成,因此可能不是。