计算序列的余弦值

Question

我需要计算以下内容：

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

totalN是一个很大的数，所以我希望能够更高效地进行计算。有没有什么方法可以提高效率？我怀疑肯定有这种方法，毕竟我们知道当n=1..N时cos(n)的结果，所以也许有某个定理可以让我更快地计算出结果。如果有任何提示，我将不胜感激。

谢谢您提前的帮助！

Federico

- Federico

---

y * i 是以弧度还是角度为单位？如果是角度，您可以使用：cos a = -1 * cos (a - 180)。如果是弧度，则使用：cos a = -1 * cos (a - pi)。y是否是一个不错的常数，只需要计算少量迭代（即需要计算的余弦值少于totalN）？ - shoover

y * i 是以弧度为单位的；问题是我需要找出是否可以利用余弦函数的周期性质。我认为我需要检查这个区间 y*[1, totalN] 是否在 [0, pi] 内，或者是否更大，如果更大，我需要找出由于周期性质而重复的点。 - Federico

除非y是pi的分数（如pi/10），否则cos的周期性可能没有帮助。 - Justin Peel

Answer 1

使用数学中最美的公式之一Euler's formula，进行替换x := n * phi: exp(i*x) = cos(x) + i*sin(x),

得到以下结果： cos(n*phi) = Re( exp(i*n*phi) ) sin(n*phi) = Im( exp(i*n*phi) ) 又因为exp(i*n*phi) = exp(i*phi) ^ 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);
  }
} 

- Curd

---

不错的解决方案，但在大量迭代后浮点误差怎么办？这是不可避免的。 - Kirk Broadhurst

尝试一下：使用我找到的Python示例，即使进行大量迭代（例如数万次），与内置cos函数计算出的值之间的差异仍保持在1E-10的范围内。 - Curd

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

- datageist

---

在我看来，这里最后一个链接中的最后一节方法似乎比这里的任何其他答案都要快。 - Justin Peel

Answer 3

也许最简单的公式是

cos(n+y) = 2cos(n)cos(y) - cos(n-y).

如果您预先计算常量2*cos(y)，那么每个值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]

- Accipitridae

---

那很漂亮，非常好的解决方案，+1 - xxxxxxx

+1非常好的解决方案，但是在迭代后很容易出现浮点误差，特别是如果totalN很大。 - Kirk Broadhurst

@Kirk Broadhurst：是的，你提出了一个重要的观点。在发布上述算法之前，我进行了一次小测试。看起来错误只会线性累积。（例如，在使用double时，经过一百万次迭代后，误差约为10^{-11}）。需要更好地分析误差。在执行一定数量的迭代后，计算一对正确的余弦值也是有意义的。 - Accipitridae

为了避免边界错误，只需每1000次迭代左右从头开始计算余弦。 - Keith Randall

Answer 4

这里有一种方法，但是它需要一点内存来存储sin值。它使用三角恒等式：

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实现的，它非常准确。

- Justin Peel

Answer 5

这里有一些不错的答案，但它们都是递归的。当使用浮点算术计算余弦函数时，递归计算将无法工作；您将不可避免地得到舍入误差，这些误差会很快累积。

考虑计算 y = 45 度，totalN 为 10,000。您最终不会得到 1 作为结果。

- Kirk Broadhurst

---

+1 你说得很有道理。然而，使用y=45会得到一个糟糕的基准。至少在我进行的小型测试中，一些舍入误差相互抵消。也就是说，当y=45时，我得到了10^-15的误差，而其他y值则会产生大约10^-13的误差。 - Accipitridae

Answer 6

为了解决柯克的疑虑：所有基于余弦和正弦递归的解决方案都可以归结为计算

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的计算机表示通常几乎是正交的，但不完全正交，因此根据所需的精度需要定期使用三角函数操作重新启动递归。这仍然比使用三角函数操作计算每个值要快得多。)

- user287792

Answer 7

你可以使用复数来实现这个功能。

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

你可以通过迭代计算所有的 i。复数乘法是 2 次乘法加上两次加法。

- user180326

---

但计算N个指数函数并不比计算N个余弦函数更快。至少没有真正的理由说明它应该更快。 - AVB

@AB：我添加了一个迭代计算它们的方法，将当前值乘以x。每个条目只需要进行一次复杂的乘法运算。 - user180326

Answer 8

知道cos(n)并没有帮助——你的数学库已经为你处理了这些微不足道的事情。

但是，如果你预先计算cos(y)和sin(y)，并在途中跟踪cos(i*y)和sin(i*y)，那么知道cos((i+1)y)=cos(iy+y)=cos(iy)cos(y)-sin(iy)sin(y)可能会有所帮助。不过这可能会导致一些精度损失，你需要进行检查。

- AVB

Answer 9

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

我不知道的是，插值是否比计算余弦更费时。由于现代CPU通常在微码中执行此操作，因此可能不是这样。

- Larry

---

上一次我检查，fcos（汇编语言）需要约35个周期，访问未缓存的内存需要几百个周期。以两种方式计时，并查看哪种更快（表格缓存或不缓存）。 - Arthur Kalliokoski

@Arthur：如果你看一下OP的代码，他想要计算cos然后更新一个数组h[i]，或许为了以后使用。所以真正的问题很可能是插值和计算之间的选择，就像Larry提到的那样，如果可以避免插值（通过选择足够大的N），这可能正是OP正在寻找的。这个解决方案还有一个优点，就是不会出现浮点累积误差。 - Aryabhatta