三角関数 n倍角 Chebyshev多項式 漸化式 FFT

Chebyshev多項式を用いると、三角関数のn倍角を漸化式で表すことができる。三角関数のn倍角を順に求めていく計算では、三角関数を計算するルーチンの呼出は最初の1回だけで、残りは漸化式で逐次的に数値計算することができる。
一松信：初等関数の数値計算、pp.129-132、1974、教育出版
上記文献によれば元ネタは下記文献らしい。
C. W. Clenshaw : A note on the summation of Chebyshev series, Mathematical Tables and other Aids to Computation, Vol.9, No.51, pp.118-120, 1955

Chebyshev多項式と漸化式は、例えば下記を参照。
潮田康夫：チェビシェフの多項式とｎ倍角の公式，数研通信，No.69，pp.8-11，2011→数研出版（PDF）
cosのn倍角はcosの多項式で表されるが、sinのn倍角はcosの多項式とsinの積で表される点に注意。

Chebyshev多項式を用いた漸化式で、三角関数のn倍角を逐次求めるC言語のソースコードは、FFTの一部として下記に記載がある。
奥村晴彦：C言語による標準アルゴリズム事典、pp.346-348、2018、技術評論社
Qiita：技術評論社 Software Technology シリーズ一覧
ソースコードに解説が無く難解だったため、当該部分のみ引用して補足させてもらう。

(1) 三角関数を計算してFFTの回転因子を求める等、下準備をする。
t = sin(PI / n);
dc = 2 * t * t; // = 1 - cos (2 * PI / n) ※0の近傍での桁落ちを避けるため、2 * sin(PI / n) * sin(PI / n) で計算している。
ds = sqrt(dc * (2 - dc)); // = sin (2 * PI / n)
t = 2 * dc;

(2) k倍角（c[k], s[k]）を漸化式で求め、次の(k+1)倍角を求めるための下準備をする。
c[k] = c[k-1] - dc; dc += t * c[k]; // 更新前の dc は -c[k] + c[k-1] の意味合い
s[k] = s[k-1] + ds; ds -= t * s[k]; // 更新前の ds は s[k] - s[k-1] の意味合い

漸化式を用いた計算では桁落ちに注意しなければならない場合もある。
吉田年雄：数値解析の基礎・基本、理工系数学の基礎・基本7、pp.125-127、2005、牧野書店（廃業）
特性方程式の解が初期値となる cos の計算は桁落ちが起きていないか心配であったが、問題ないと解説されている。