计算序列的余弦

Question

我必须计算以下内容：

float2 y = CONSTANT;
for (int i = 0; i < totalN; i++)
   h[i] = cos(y*i);

totalN 是一个很大的数字，所以我想以更有效的方式进行计算。有什么办法可以改善这一点吗？我怀疑是有的，因为毕竟我们知道 n=1..N 时 cos(n) 的结果是什么，所以也许有一些定理可以让我以更快的方式计算这个结果。我真的很感激任何提示。

预先感谢，

费德里科

Answer 1

使用最美丽的数学公式之一，欧拉公式
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 次重复乘法。
因此，您可以通过重复复数乘以 exp(i*phi) 来计算 cos(n*phi) 并同时计算 sin(n*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);
  }
} 

Answer 2

我不确定您愿意在精度与性能之间做出什么样的妥协，但在这些链接中对各种正弦曲线近似技术进行了广泛的讨论：

正弦曲线的乐趣 - http://www.audiomulch.com/~rossb/code/sinusoids/
快速准确的正弦/余弦 - http://www.devmaster.net/forums /showthread.php?t=5784

编辑（我认为这是“Fun with Sinusoids”页面上损坏的“Don Cross”链接）：

优化三角函数计算 - http://groovit.disjunkt.com/analog/time-domain/fasttrig.html

Answer 3

也许最简单的公式是

cos(n+y) = 2cos(n)cos(y) - cos(ny)。

如果预先计算常数 2*cos(y)，则可以通过一次乘法和一次减法根据前 2 个值计算每个值 cos(n+y)。
即，用伪代码

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]

Answer 4

这是一个方法，但它使用了一点内存。它使用三角恒等式：

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实现了这个，而且非常准确。

Answer 5

这里有一些很好的答案，但它们都是递归的。使用浮点运算时，递归计算不适用于余弦函数；您总是会遇到舍入误差，这些误差很快就会复合。

考虑计算 y = 45 度，totalN 10 000。最终结果不会是 1。

Answer 6

为了解决 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) = Re (k - 1) 通过线性关系。由于R是所谓的正交矩阵，因此R e(k - 1) 与e(k - 1) 具有相同的范数，并且误差增长非常缓慢。 （它增长的原因是由于舍入；R 的计算机表示一般来说几乎是正交的，但不完全正交，因此有必要根据精度不时地使用三角运算重新启动递归这仍然比使用三角运算来计算每个值要快得多。）

Answer 7

您可以使用复数来做到这一点。

如果定义 x = sin(y) + i cos(y)，cos(y*i) 将是 x^i 的实部。

您可以迭代计算所有 i 。复数乘法是 2 次乘法加两次加法。

Answer 8

知道 cos(n) 并没有帮助——你的数学库已经为你做了这些琐碎的事情。

知道 cos((i+1)y)=cos(iy+y)=cos(iy)cos(y)-sin(iy)<如果您预先计算 cos(y) 和 sin(y)，并一路跟踪 cos(iy) 和 sin(i*y)，em>sin(y) 会有所帮助。不过，这可能会导致一些精度损失 - 您必须进行检查。

Answer 9

您需要得到的 cos(x) 有多准确？如果您可以忍受一些，您可以创建一个查找表，以 2*PI/N 间隔对单位圆进行采样，然后在两个相邻点之间进行插值。将选择 N 以达到某种所需的精度水平。

我不知道插值是否实际上比计算余弦成本更低。由于它通常在现代 CPU 中以微代码完成，因此可能不是。