円周上の点を有理点で近似したかった

講義のレポートで、 $\LaTeX$ で図を書く時に必要になったので、メモ。
$\TeX$ の貧弱な数値計算環境では、sin, cos を気軽に計算することは、多大なコストを要する。
そこで、出来るだけ簡単な有理点で近似してやろうということになった。

アイデアは簡単で、三角関数の間の関係式を駆使すると、 $\alpha=\tan(\frac{\theta}{2})$ とすると、次が成立する。

$\displaystyle\cos\theta = \frac{1-\alpha^2}{1+\alpha^2}\quad,\qquad\sin\theta=\frac{2\alpha}{1+\alpha^2}$

ここで、 $\alpha=p/q$ なる有理数ならば、

$\displaystyle\cos\theta = \frac{q^2-p^2}{p^2+q^2}\quad,\qquad\sin\theta = \frac{2pq}{p^2+q^2}$

となり、中々きれいな形になることが知られている。

$\tan(\frac{\alpha}{2})$ の近似は、出来るだけ簡単な有理数になって欲しいので、Farey 列を用いる。
Farey 列については、Farey sequence - Wikipedia, the free encyclopedia を参照して下さい。
Haskell で書いてみる：

farey n = farey' n [(0,1),(1,1)]
    where farey' n (h:[]) = [h]
          farey' n ((a,b):t@((c,d):_)) = if b+d<= n
                                         then farey' n $ (a,b):(a+c,b+d):t
                                         else (a,b):farey' n t

ratioapprox n x = rapprox x $ farey n
    where rapprox _ (y:[]) = y
          rapprox x (y:t@((r,s):_)) = if s*x<=r then near x y (r,s)
                                                else rapprox x t
          near a (b,c) (d,e) = case compare (a*c*e-b*e) (c*d-a*c*e) of
                               LT -> (b,c)
                               GT -> (d,e)
                               EQ -> if (c <= e) then (b,c) else (d,e)

ratiotriple n t = restriple
    where (p,q) = ratioapprox n $ tan (t/2)
          restriple = (q*q-p*p,2*p*q,p*p+q*q)

関数 ratiotriple は、精度 n と角度 t を受け取り、 $(\cos t,\sin t)$ に近くにある円周上の有理点に対応するピタゴラスの三数を返す。
精度というのは、 $\tan(\frac{t}{2})$ の近似値を出す時、分母をどれだけ大きく取るかの上限。
色々試してみるとわかるが、必ずしも結果の三数は互いに素でない。
つまり、非自明な約分で結果的に簡単な比になる場合を見落としている可能性はあるが、ご愛嬌。

ghci で動かすと以下のような結果に

*Main> ratiotriple 10 (pi/12)
(63.0,16.0,65.0)

ちなみに

*Main> 180*atan (16/63)/pi
14.250032697803595

なので、まあ、運が悪くなければ悪くはない精度か。

（追記）
しかし、どうせハードコーディングにならざるを得ないので、10進展開でも問題無いということに気付いてしまった。
Haskell とか C++ のテンプレートに慣れ過ぎているせいで、コンパイラがこれ程計算を嫌がる処理系で、上手に息ができない病気になってしまっている。
面白かったから、記事自身は残しておきます。