您很可能从学校知道以下比率：

$$

$$\ sin（\ alpha + \ beta）= \ sin \ alpha \ times \ cos \ beta + \ cos \ alpha \ times \ sin \ beta \\ \ cos（\ alpha + \ beta）= \ cos \ alpha \ times \ cos \ beta-\ sin \ alpha \ times \ sin \ beta$$

当您在童年时期第一次熟悉此公式时，很可能首先要感到痛苦是由于必须记住该公式。 这非常糟糕，因为实际上您不需要记住该公式-在纸上旋转三角形时会显示该公式。 实际上，写下这个公式时我会做同样的事情。 这篇文章的中间部分将使这种解释变得显而易见。 但是现在，为了把所有的乐趣留给以后再推迟说“ Eureka！”的时间，让我们考虑为什么我们甚至应该考虑这个公式。

参赛作品

三角函数sin()和cos()在计算机图形学cos()可能是最流行的，因为它们是以参数方式描述任何圆形的基础。 在可能的应用场合中，当计算傅立叶变换时会生成圆或旋转的体积对象，在水平面上生成波的程序生成，用于软件声音合成器的生成器等等。 在所有这些情况下，都在循环内调用sin()和cos() ，如下所示：

 for(int n=0; n < num; n++) { const float t = 2.0f*PI*(float)n/(float)num; const float s = sinf(t); const float c = cosf(t); // do something with s and c ... } 

我们开始以增量方式重写周期（请参见下面的代码），因此我们更容易想象在周期为t n该周期的迭代n ，下一个迭代n+1将考虑t+f sin()和cos() 。 换句话说，为我们计算了sin(t)和cos(t) ，我们需要计算sin(t+f)和cos(t+f) ：

 const float f = 2.0f*PI/(float)num; const float t = 0.0f; for( int n=0; n < num; n++ ) { const float s = sinf(t); const float c = cosf(t); // do something with s and c ... t += f; } 

无论如何获取t及其值的范围是什么（在上面的示例中- $$$[0; 2 \ pi]$ ） 唯一令我们困扰的是，存在一个循环，该循环不断调用sin()和cos() ，并带有一个以恒定步长递增的参数（在这种情况下， $$$\ frac {2 \ pi} {\ text {num}}$ ） 本文是关于如何优化此代码的速度，从而完全可以执行相同的计算，而无需使用sin()或cos()函数（在内循环中），甚至不使用更快的sincos()函数。

但是，如果您看一下文章中的第一个公式，我们将看到 $$$f = \ alpha$ 和 $$$t = \ beta$ 我们可以改写成

 sin(t+f) = sin(f)*cos(t) + cos(f)*sin(t) cos(t+f) = cos(f)*cos(t) - sin(f)*sin(t) 

或者换句话说

 new_s = sin(f)*old_c + cos(f)*old_s new_c = cos(f)*old_c - sin(f)*old_s 

由于f是常数，因此sin(f)和cos(f)也是。 分别称它们为a和b ：

 new_s = b*old_c + a*old_s new_c = a*old_c - b*old_s 

可以在我们的代码中直接使用该表达式，然后得到一个循环，在该循环中调用了昂贵的（实际上没有）三角函数！

 const float f = 2.0f*PI/(float)num; const float a = cosf(f); const float b = sinf(f); float s = 0.0f; float c = 1.0f; for( int n=0; n < num; n++ ) { // do something with s and c ... const float ns = b*c + a*s; const float nc = a*c - b*s; c = nc; s = ns; } 

口译

到目前为止，我们一直在盲目地玩数学，而不是了解实际情况。 让我们这样重写内部循环：

$$

$$s_ {n + 1} = s_na + c_nb \\ c_ {n + 1} = c_na-s_nb$$

可能有些人已经注意到，这是在二维空间中旋转对象的公式。 如果您仍然不明白这一点，那么矩阵形式可能会为您提供帮助。

$$

$$\ left（\ begin {matrix} s_ {n + 1} \\ c_ {n + 1} \ end {matrix} \ right）= \ left（\ begin {matrix} a＆b \\ -b＆a \ end {matrix} \ right）\ cdot \ left（\ begin {matrix} s_ {n} \\ c_ {n} \ end {matrix} \ right）$$

实际上，可以将sin(t)和cos(t)分组为长度为1的矢量，并将其绘制为屏幕上的一个点。 将此向量称为x 。 然后 $$$x = \ {\ sin \ beta，\ cos \ beta \}$ 。 所以表达的向量形式是 $$$x_ {n + 1} = Rx_n$ 在哪里 $$$R = \左（\开始{matrix} a和b \\-b＆a \结束{matrix} \ right）$ 。 我们看到我们的循环在每次迭代中执行一个小的二维旋转，以便在循环执行过程中x绕圆旋转。 每次迭代x旋转 $$$\ frac {2 \ pi} {\ text {num}}$ 弧度
所以基本上

$$

$$\ sin（\ alpha + \ beta）= \ sin \ alpha \ times \ cos \ beta + \ cos \ alpha \ times \ sin \ beta \\ \ cos（\ alpha + \ beta）= \ cos \ alpha \ times \ cos \ beta-\ sin \ alpha \ times \ sin \ beta$$

这是点运动公式 $$$x = \ {\ cos \ alpha，\ sin \ alpha \}$ 周长以 $$$\ beta$ 弧度 为此，我们将构造两个旋转轴之一，例如， $$$u = \ {\ cos \ beta，\ sin \ beta \}$ 。 旋转的第一部分是投影。 $$$x$ 在 $$$u$ 。 由于 $$$x$ 和 $$$u$ 归一化（长度为1），则投影为其标量积。 因此 $$$s = x \ cdot u = \ sin \ alpha \ cdot \ cos \ beta + \ cos \ alpha \ cdot \ sin \ beta$ ，当然第二个成分是反投影，可以通过投影到垂直轴上找到它， $$$v$ 。 我们可以通过扩展坐标来创建此向量 $$$u$ 并在第一个坐标处将符号更改为相反的符号： $$$c = x \ cdot v = \ cos \ alpha \ cdot \ cos \ beta + \ sin \ alpha \ cdot \ sin \ beta$

注意事项

通常，您应该能够一遍又一遍地执行这些微小的旋转。 其实 $$$| R | = \左| \开始{matrix} a和b \\-b＆a \结束{matrix} \ right | = a ^ 2 + b ^ 2 = \ sin ^ 2 \ alpha + \ cos ^ 2 \ alpha = 1$ ，这表示矩阵 $$$R$ 不会增加或减少其应用的空间，这意味着 $$$x$ 将移动一个完美的圆圈。 但是，由于计算机不准确， $$$x$ 将以螺旋形运动，并最终与旋转圆心重合。 我没有这样的问题，但我认为它们可能以非常大的数量出现，即 小转弯角。

例子

在Kindercrasher （2008年以来的4096字节演示）（KDPV上的截图）中，一群球体向音乐发出脉动。 为此，我计算了声音的傅立叶变换。 有一些算法可以实时执行此操作，例如FFT 。 但是，我的代码应该可以容纳几千字节，因此我决定采用另一种方式。 我没有实现FFT，而是通过简单的定义编写了DFT 。 您可以在Wikipedia上查看。

$$

$$X_k = \ sum_ {n = 0} ^ {N-1} {x_ne ^ {-\ frac {2 \ pi i} {N} kn}} \ quad k = 0,1，\ ldots，N-1$$

我的函数还使用16位立体声音频缓冲区x ，并返回声音y频谱的前128个频率。 查看内部循环的组织方式，该循环运行4096次：不是对sin()或cos()函数的单个调用，尽管在其他实现中，这些调用也会如此。

 #include <math.h> void iqDFT12( float *y, const short *x ) { for( int i=0; i<128; i++ ) { const float wi = (float)i*(2.0f*3.1415927f/4096.0f); const float sii = sinf( wi ); const float coi = cosf( wi ); float co = 1.0f; float si = 0.0f; float acco = 0.0f; float acsi = 0.0f; for( int j=0; j<4096; j++ ) { const float f = (float)(x[2*j+0]+x[2*j+1]); const float oco = co; acco += co*f; co = co*coi - si*sii; acsi += si*f; si = si*coi + oco*sii; } y[i] = sqrtf(acco*acco+acsi*acsi)*(1.0f/32767.0f); } }