C言語による最新アルゴリズム事典 FFT 三角関数表の計算方法を調べた

# FFT向け三角関数表の作成
以下は 奥村 晴彦 著 "[C言語による最新アルゴリズム事典](https://www.amazon.co.jp/C%E8%A8%80%E8%AA%9E%E3%81%AB%E3%82%88%E3%82%8B%E6%9C%80%E6%96%B0%E3%82%A2%E3%83%AB%E3%82%B4%E3%83%AA%E3%82%BA%E3%83%A0%E4%BA%8B%E5%85%B8-%E3%82%BD%E3%83%95%E3%83%88%E3%82%A6%E3%82%A7%E3%82%A2%E3%83%86%E3%82%AF%E3%83%8E%E3%83%AD%E3%82%B8%E3%83%BC-%E5%A5%A5%E6%9D%91-%E6%99%B4%E5%BD%A6/dp/4874084141)"に記載のあるFFT計算のコードの一部です。
FFTを実行する際に利用する三角関数表(sinテーブル)を作成しています。
sinを一度求めたら、それを使って次々にcos,sinのn倍角の値を計算しています。
が、書籍中に特に説明はなく、なぜこれで計算できるのかすぐにはわかりません。
（FFTのソースコードは[ここ](https://github.com/okumuralab/algo-c/blob/master/src/fft.c)にあります。）

```C
/*
  関数{\tt fft()}の下請けとして三角関数表を作る.
*/
static void make_sintbl(int n, double sintbl[])
{
    int i, n2, n4, n8;
    double c, s, dc, ds, t;

    n2 = n / 2;  n4 = n / 4;  n8 = n / 8;
    t = sin(PI / n);
    dc = 2 * t * t;  ds = sqrt(dc * (2 - dc));
    t = 2 * dc;  c = sintbl[n4] = 1;  s = sintbl[0] = 0;
    for (i = 1; i < n8; i++) {
        c -= dc;  dc += t * c;
        s += ds;  ds -= t * s;
        sintbl[i] = s;  sintbl[n4 - i] = c;
    }
    if (n8 != 0) sintbl[n8] = sqrt(0.5);
    for (i = 0; i < n4; i++) {
        sintbl[n2 - i] = sintbl[i];
        sintbl[i + n2] = - sintbl[i];
    }
}
```

### 三角関数を計算している中心的な部分
上のコードの最も肝心なところを以下に抜粋しました。

```C
    t = sin(PI / n);
    dc = 2 * t * t;  ds = sqrt(dc * (2 - dc));
    t = 2 * dc;  c = sintbl[n4] = 1;  s = sintbl[0] = 0;
    for (i = 1; i < n8; i++) {
        c -= dc;  dc += t * c; // cosに関する計算
        s += ds;  ds -= t * s; // sinに関する計算
        sintbl[i] = s;  sintbl[n4 - i] = c; // 求めたcos,sinの値をテーブルへ設定
    }
```

上の部分で$\theta =2\pi/N $ として
変数 c に$\cos(n\theta)$、変数 s に$\sin(n\theta)$ が計算されているようです。（実際計算できている）
これについて以下で示しました。

# チェビシェフ多項式
$\cos(\theta\)$ のn倍角の値は[チェビシェフ多項式](https://ja.wikipedia.org/wiki/%E3%83%81%E3%82%A7%E3%83%93%E3%82%B7%E3%82%A7%E3%83%95%E5%A4%9A%E9%A0%85%E5%BC%8F)で表せる。

```math
\begin{align}
\cos(1t) &= \cos(t) \\
\cos(2t) &= \cos(t)\cos(t) - \sin(t)\sin(t)=2(\cos(t))^2 - 1\\
\cos(3t) &= 4(\cos(t))^3-3\cos(t) \\
\vdots&
\end{align}
```

以下の式で$x = \cos(t)$ とすれば、

```math
\begin{align}
T_0(x)& = 1 \\
T_1(x)& = x \\
T_2(x)& = 2x^2-1 \\
T_3(x)& = 4x^3-3x \\
\vdots&
\end{align}
```

$T_n(x)$を漸化式で表せば

```math
T_{n+1}(x) = 2xT_n(x)-T_{n-1}(x)
```

と表せる。これを利用する。

# cosのn倍角について
$k = 2\sin^2(\theta/2)$とすると 
$\cos(\theta) = 1 - 2\sin^2(\theta/2) = 1 - k$
また、
$ T_n(\cos\theta) $ を $T_n$と書けばチェビシェフ多項式を使って

```math
\begin{align}
T_0 &= 1 \\
T_1 &= 1 - k = \cos(\theta) \\
T_2 &= 2(1 - k)T_1 - T_0 \\
    &= 2T_1 - 2kT_1 - T_0 = T_1 - 2kT_1 + T_1 - T_0 \\
    &= T_1 - 2kT_1 - k \\
T_3 &= 2(1 - k)T_2 - T_1 \\
    &= 2T_2 - 2kT_2 - T_1 = T_2 - 2kT_2 + T_2 - T_1 \\
    &= T_2 - 2kT_2 - 2kT_1 - k \\
\vdots&
\end{align}
```

以上のように書ける。これを見ると以下の式が成り立っているように見える。

```math
T_{n+1} = T_n - \sum_{i=1}^{n} 2kT_i - k  \tag{1}
```
### 式(1)がn>=1で成り立つことの証明
式$(1)$は、$n=1$のとき成り立つ。
$n=p$のとき成り立つとすると、$n=p+1$は、

```math
\begin{align}
T_{p+2} &= 2(1 - k)T_{p+1} - T_p \\
        &= 2T_{p+1} - 2kT_{p+1} - T_p \\
        &= T_{p+1} - 2kT_{p+1} + T_{p+1} - T_{p} \\
        &= T_{p+1} - 2kT_{p+1} + (T_p - \sum_{i=1}^{p} 2kT_i - k) - T_p \\
        &= T_{p+1} - \sum_{i=1}^{p+1} 2kT_i - k
\end{align}
```

つまり、$n=p+1$で成り立つ。

## プログラム中のcos計算との関係

```C
c -= dc;  dc += t * c;
```

```math
T_{n+1} = T_n - \sum_{i=1}^{n} 2kT_i - k  \tag{1}
```

上記のコードの、変数 c は式$(1)$に従って計算されていて、$T_n$に相当する。
変数 t は$2k$に相当し、変数dcは

```math
\sum_{i=1}^{n} 2kT_i + k
```

に相当することがわかる。

# sinのn倍角について
sinについてもcosと同様に求める。
$k = 2\sin^2(\theta/2) $とし、
$l = \sqrt{k(2 - k)}$
とすると、

```math
\begin{align}
l &= \sqrt{k(2 - k)} = \sqrt{2\sin^2(\theta/2)(2-2\sin^2(\theta/2)} \\
&=\sqrt{4\sin^2(\theta/2)-4\sin^4(\theta/2)} \\
&=2\sin(\theta/2)\cos(\theta/2) \\
&=\sin(\theta) \\
&=\cos(\pi/2-\theta)
\end{align}
```

ここで、チェビシェフ多項式$T_n$を使って$S_n$を

```math
S_n = T_{N/4-n}
```

として、$\pi/2$から右回りに$\theta$角ごとの$\cos$を求める。
これはつまり、$\sin(n\theta)$を求めることになる。

```math
\begin{align}
S_0 &= \cos(\pi/2) = 0 \\
S_1 &= \cos(\pi/2 - \theta) = \sin(\theta) = l \\
S_2 &= 2(1 - k)S_1 - S_0 \\
&= 2S_1 - 2kS_1 - S_0 = S_1-2kS_1+S_1-S_0 \\
&= S_1 - 2kS_1 + l \\
S_3&= 2(1-k)S_2-S1\\
&=2S_2-2kS_2-S_1 =S_2-2kS_2+S_2-S_1  \\
&= S_2-2kS_2-2kS_1+l \\
\vdots&
\end{align}
```

これを見ると以下の式が成り立っているように見える。

```math
S_{n+1} = S_n + \sum_{i=1}^{n} -2kS_i + l  \tag{2}
```

### 式(2)がn>=1で成り立つことの証明
式$(2)$は$n=1$のときに成り立つ。$n=p$のときに成り立つとすると、
$n=p+1$について

```math
\begin{align}
S_{p+2} &= 2(1 - k)S_{p+1} - S_p \\
        &= 2S_{p+1} - 2kS_{p+1} - S_p \\
        &= S_{p+1} - 2kS_{p+1} + S_{p+1} - S_{p} \\
        &= S_{p+1} - 2kS_{p+1} + (S_p + \sum_{i=1}^{p} -2kS_i + l) - S_p \\
        &= S_{p+1} + \sum_{i=1}^{p+1} -2kS_i + l
\end{align}
```

つまり、$n=p+1$で成り立つ。

## プログラム中のsin計算との関係

```C
s += ds;  ds -= t * s;
```

```math
S_{n+1} = S_n + \sum_{i=1}^{n} -2kS_i + l  \tag{2}
```

上記のコードの、変数 s は式$(2)$に従って計算されていて、$S_n$に相当する。
変数 t は$2k$に相当し、変数 ds は式$(2)$で考えれば、

```math
\sum_{i=1}^{n} -2kS_i + l
```

に相当することがわかる。

以上で三角関数テーブルの計算方法がわかった。