自前のサインテーブル
Boreasでは、サインテーブルを再帰展開で作成しています。
ここでいう再帰とは、直前に求めた値を使って次の値を求めるという程度の意味です。
わりと古くからある手法ですが、一般的なPC上のプログラム開発では開発言語が提供する
三角関数を使う方が楽なので、よっぽどの事情がないかぎり使われません。
●原理
まず、Sin(２Nα）という式を考えます。
（円に内接する正多角形の図を書くとわかりやすいのですが、省略）
αは単位となる固定値の角度で、Nが1,2,3,4…と変化します。
このとき、各Nにおける式の値は、
Sin０α＝０
Sin２α＝Sin０α＋２Sinα・Cos（０α＋α）
Sin４α＝Sin２α＋２Sinα・Cos（２α＋α）
Sin６α＝Sin４α＋２Sinα・Cos（４α＋α）
…
一般式 Sin*1で示されることがわかります。
これらについても、各Nにおける式の値は、
Cos（２α＋α）＝Cos（０α＋α）−２Sinα・Sin２α
Cos（４α＋α）＝Cos（２α＋α）−２Sinα・Sin４α
Cos（６α＋α）＝Cos（４α＋α）−２Sinα・Sin６α
…
一般式 Cos*2を加えて次のSin値(A1)を求めるわけですが、
Sin(A1)とSin(A)の差分は以下のようにして導くことができます。
「Sin(A1)-Sin(A)=差分」としたとき、A=B-c、A1=B+cとおくと、加法定理より
Sin(B+c)-Sin(B-c)=Sin(B)Cos(c)+Cos(B)Sin(c)-( Sin(B)Cos(c)-Cos(B)Sin(c) )
=Sin(B)Cos(c)+Cos(B)Sin(c)-Sin(B)Cos(c)+Cos(B)Sin(c)
=2Cos(B)Sin(c)
というわけで Sin(B+c)=Sin(B-c)+２Cos(B)Sin(c) の解が成り立ちますので、
B-cを分割幅２αのN倍＝２Nα、cを分割幅の半分＝αとおくと、上の一般式になります。
Sin((B-c)+2c)＝Sin(B-c)＋２Sin(c)・Cos((B-c)+c)
Sin(2Nα+2α)＝Sin2Nα＋２Sinα・Cos（2Nα＋α）
Sin((2N+2)α＝Sin2Nα＋２Sinα・Cos((2N+1)α)
（Cosも同様に、加法定理で2つのCosの差を求めると、上の式が導けます）
さて、これらの式を、αを整理して左辺関数のαの係数順に並べてみると
Sin０α＝０
Cos１α＝定数
Sin２α＝Sin０α＋２Sinα・Cos１α
Cos３α＝Cos１α−２Sinα・Sin２α
Sin４α＝Sin２α＋２Sinα・Cos３α
Cos５α＝Cos３α−２Sinα・Sin４α
Sin６α＝Sin４α＋２Sinα・Cos５α
Cos７α＝Cos５α−２Sinα・Sin６α
となります。
3行目以降をみると、1行上と2行上で左辺にあった項が使用されていることが判ります。
また、毎行に２Sinαが使われています。
つまり、初期パラメータとしてCosαとSinαが与えられれば、繰り返し計算することで
残りが求まるということになります。
●スクリプト化
式はSin,Cosそれぞれ２α毎なので、Sin２Nαのみを利用し、NについてのSinテーブルとします。
Cosは位相がSinの90度進みになりますので、求めるのはSinだけで十分となります。
1周256分割の場合を考えると、αは分割幅の半分なので、α＝２π/（２・２５６）
２Sin（２π/５１２）は、だいたい 4/163 となります。
精度として1/2000程度（0〜1を0〜2000で対応）を考えることにします。
Cosαは、0.9999247…という値になるのですが、面倒なのでCos0と同じ扱いにします。
パラメータはそろいましたので、1/4円周分だけ求めるスクリプトを書いてみます。
1/4の部分で精度を保つことができれば、残りの3/4はコピーと符号反転で補完できます。
このスクリプトは手順を正直に書いたものです。HSP wikiにもほぼ同じコードがあります。
dim sintbl,256
zcos=32768000
zsin=0
repeat 65
sintbl.cnt=zsin/16384
zsin += zcos*4/163
zcos -= zsin*4/163
loop
HSP2.6標準では整数型しか扱えませんので精度をもたせるには下駄を履かせる必要がありますが、
下駄が大きいため配列へ格納する際に16384で割っています。
その計算を考慮し、zcosの初期値は2000*16384=32768000としました。
切捨てなども発生する結果、1/4周のSin値は0〜1999になります。
なお、ここで選ばれているパラメータは精度よくSinが得られる組み合わせであり、1周分を計算
しても大きくずれません。-Sinθ＝Sin-θや、Sinθ＝-Sin(π＋θ)を厳密に成立させる必要が
なければ、repeat 256でも問題ないでしょう。
ちなみに、HSP2.Xではオペランド補完と、リテラルを0〜511におさえる手法でサイズを減らす
ことが可能です。（下記のスクリプトはAXファイルにして64バイト分です）
求めた各Sin値は、
dim sintbl,256
zcos=40750000 ;下の計算に合わせて初期値調整
;zsin=0 初期値0であることを見越して初期化省略 6バイト減
repeat 256
sintbl.cnt=zsin/500 ;割る数を511以下に抑える(4バイト減)
zsin + zcos*16/163/ ;← /163/ は/163/163と同じ。4バイト増のところ2バイト増。
zcos - zsin ;*4/163 がなくなったので8バイト減。
loop ; +=,-=を + - に変えたのと合わせて、計20バイト減。
kakudo=kakudo&255:kak90p=kakudo+192&255
sin=sintbl.kakudo:cos=sintbl.kak90p
のようにしてテーブルから直接参照できます。(下駄には注意）
Boreas Ver2.26では分解能を128にしたかったので、1/4周を求めて残りはCopyしています。
(Sin展開だけを見るとサイズは増えますが、他の処理も巻き込んでループを節約すること
を取りました。)
