C言語なんかにどっぷり浸かっていると、sin()関数とかexp()関数とか数学関数は普通に使えるし、なにも困ることは無いわけです。でもそれって、何方かが数学関数ライブラリを整備してくれているからこそ、なので、オリジナル16ビットCPUコアMoscoviumなんて作って見たりすると、はたと困ってしまうのです。
そこで今回は、「一般的な数学関数くらい自分で作らなきゃね」をテーマとします。
使用するのは、半精度浮動小数演算ユニット(H-FPU)を搭載した、MCOC115FP0440です。デザインデータは以下リンクよりダウンロード可能です。
コンピュータで数学関数を作る、と言ってまず思いつくのはテイラー展開と言う方法です。
円周率を計算した時に使ったのは、arctan (1/x)のグレゴリー級数でした(↓)。
上のまとめに従うと、テイラー級数の一種だった、と言うことになります。
数学関数の値を求める時は、上の様な級数に値を代入して、関数値を求めます。
ところがこれらの級数は、項番nの値(上の数式だと項番k)が無限大まで続くのです。コンピュータでは無限大と言う数値を扱うことが出来ないので、適切なところで打ち切るしかありません。では、どのくらいの値で打ち切るのか？と言うと、「必要な演算精度が得られたところ」で打ち切る必要があります。
今回は半精度浮動小数half型がターゲットで、10進数換算の有効数字は3桁程度しかありません。たったの3桁で良いので、さぞかし少ない項数で十分な精度が得られるのではないか？と期待します。
ところが、実際にはそう上手くは行かない様です。後で関数毎に見てみましょう。
一般的には、テイラー級数第n項の数値がどれだけ小さいか？で、収束速度が異なります。収束の遅い関数は、高次項まで計算する必要があるので計算が遅くなります。
まず、今回作成する関数ライブラリの目標仕様を以下の通りとします。
おそらく、演算速度や演算精度を必要とする用途では、根本的に対処すべき問題が異なるので、このページには到達しないであろうと推察します。
早速、三角関数から作って行きます。三角関数の中で今回実際に作成するのは正弦sin()関数だけです。sin()関数が有れば、余弦cos()関数と正接tan()関数は以下式で求められます。
sin()関数のテイラー展開は以下式となります。
第n項目の数値は、分母に(2n＋1)！←階乗という比較的大きな数値が来るため、小さくなりそうな気がしますが、問題は分子にあるx2n＋1です。sin()関数のxは全ての区間(−∞〜∞)で定義されるので、例えばx＝1,000などの大きな数値だった場合、その2n＋1乗の数値は分母の(2n＋1)！を打ち消すに余りあるほど大きくなってしまいます。
これでは収束どころの話ではありません。xの値をなるべく小さくすることが必要となります。
sin()関数は0〜2πの区間で1周期となり、後は2π毎に同じ波形が繰り返し現れる周期関数です(下図)。
しかも上図の様に、赤色部分、水色部分、黄色部分、緑色部分は左右・上下に反転していますが全て同じ形をしています。この様なグラフの相似特性を使うと、x＝0〜π/2までの範囲の値を求めさえすれば、その値を左右・上下に反転することで全てのxに関してsin (x)を計算することが出来ます。
sin (x)を求めるには、
x’＝x％(2π) …… 2πで割った剰余x’として、
x’’＝x’％(π/2) …… 更にπ/2で割った剰余x’’に関して、
sin (x’’)を求めて、左右・上下に反転すれば良い。
この時、0≦x’’＜π/2
こうすることでx’’の最大値はπ/2≒1.57程度となり、テイラー級数の分母(2n＋1)！と比べて十分小さくなります。half型の精度(有効数字3桁)だと、n＝0〜3までのテイラー級数で、十分な精度を得られることが解りました。
こうして完成したsin()関数をiverilogによるシミュレーションで計算させてみました(H-FPUが必要なので、MCOC115FP0440のデザインを使ってシミュレーションしました)。
sin()関数テストプログラム：mcvm_hfpu_sin.asm
数学関数ライブラリ：mcvm_lib_math.incl
実験サンプル：【 】内はCygwinシェルウインドウからの実行コマンド
【iverilog -f test.lst ; if ( ! $status ) vvp a.out】
：
VCD info: dumpfile test.vcd opened for output.
sin
-6.906250 -0.584961
-6.796875 -0.492920
-6.687500 -0.395020
-6.578125 -0.292480
：
6.578125 0.292480
6.687500 0.395020
6.796875 0.492920
6.906250 0.584961
Test Pass
※シミュレーション用プリントユニットを使っているので、実機(FPGAボード)で実行しても、何も表示されません。表示された値は、「<xの値> <sin (x)の値>」です。
こうして(↑)計算させた値をExcelでグラフにしたのが下図です。
青い線(左軸)がsin (x)です。グラフだと非常にきれいな正弦波に見えますが、赤線の様に期待値からの誤差率(右対数軸)を取って見ると、sin (x)≒0となる付近の誤差が大きいことが解りました。
誤差率は％(＝1/100)ではなく、parts per thousand(＝1/1,000)となっています。有効数字3桁の精度が満たされている場合、誤差率(1/1,000)は1未満の値となります。よって、「目標誤差率(1/1,000)は1未満」です。
改めて赤い線を見てみると、ほとんどの部分は1未満ですが、x＝π、x＝2πとなってsin (x)≒0となる周辺で1以上となって、有効数字3桁を下回っています。
しかも、その誤差率はｘ＝π→2π→3π→…と大きくなるにつれて、大きくなって行く様に見えます(緑矢印)。
以下の様なC言語のプログラムを実行してみましょう。
#include <stdio.h>
#include <math.h>
int main (void)
{
int idx;
float x;
for (idx=1; idx<=65536; idx*=2)
{
x=2.*M_PI*idx;
printf ("sin(2*pi*%d)=%f\n",idx,sin (x));
}
return (0);
}
ずっと、sin (2π×idx)を計算して表示しています。sin()関数は、x＝2π×idx(ただしidx＝整数)では、常にゼロとなります。確かめてみましょう(↓)。
実験サンプル：【 】内はCygwinシェルウインドウからの実行コマンド
【gcc test_sin_large_x.cc ; if ( ! $status ) a.exe】
sin(2*pi*1)=0.000000
sin(2*pi*2)=0.000000
sin(2*pi*4)=0.000001
sin(2*pi*8)=0.000001
sin(2*pi*16)=0.000003
sin(2*pi*32)=0.000006
sin(2*pi*64)=0.000011
sin(2*pi*128)=0.000022
sin(2*pi*256)=0.000045
sin(2*pi*512)=0.000090
sin(2*pi*1024)=0.000179
sin(2*pi*2048)=0.000358
sin(2*pi*4096)=0.000716
sin(2*pi*8192)=0.001432
sin(2*pi*16384)=0.002865
sin(2*pi*32768)=0.005729
sin(2*pi*65536)=0.011458
最初の方こそ、sin (2π×idx)≒0ですが、sin (2π×65,536)≒0.0115と、さすがにゼロとは言い難い値になってきます。実は、C言語の数学関数ライブラリでさえ、こんな感じなのです。
変数xの宣言をfloatからdoubleに変更すると、だいぶマシになります。要は、浮動小数型の仮数部のビット数(解像度)によって、程度の差はあれども、x＞＞2πとなる大きな値に対するsin()関数値の精度はどんどん低下するのです。
仮数部が10ビットしかないhalf型ではこれがより顕著に見えてきて、x＝2πでも、誤差率(1/1,000)が100近くに到達してしまいました(有効数字1桁分くらいの精度しかない！)。
※sin()関数の計算アルゴリズムによる精度低下ではなく、浮動小数型仮数部の解像度(ビット数)による精度低下であることに留意の事。
誤差率(1/1,000)グラフから判断するに、abs (x)＜＜πの範囲であれば、そこそこ使えると思います。ですがちょっと、half型で三角関数の計算は厳しいかもしれません。
指数exp()関数のテイラー展開は、以下の通りです。
exp()関数も、全ての範囲(x＝−∞〜∞)で定義されていて、しかも級数の分子側にxnが存在するため、収束の遅い関数となっています。つまり、必要な精度を確保するためには大きな項番nまで計算しなければなりません。
ここでもxの値をなるべく小さくする工夫が必要です。まず、引数xの整数部分xi、小数部分xfとおきます。
x＝xi＋xf …… ただし、xi＝int (x)、xf＝x−xi
exp (x)＝exp (xi＋xf)＝exp (xi)×exp (xf)
上の式より、exp (xi)は、ネイピア数e≒2.718…の整数べき乗ですので、適切なループを回しながらeのべき乗を計算すれば求めることが出来ます。half型なら、xi＝12程度でexp (12)≒162,755＝INF(無限大)となるので、高々12個前後のテーブル参照にできます。ループして掛け算するよりも高速でしょう。
残るexp (xf)は、0≦xf＜1の範囲です。テイラー級数の第n項では、分子0≦xfn＜1に対して、分母n！なので、級数は比較的速く0に収束して行きます。それでも、有効数字3桁の精度を得るためには、n＝0〜6程度までの級数を計算する必要がありました。
そこで以下の様にして、さらにxfの最大値を小さくする工夫を実施しました。
x＝xi＋xf …… ただし、xi＝int (2×x)/2、xf＝x−xi
0≦xf＜0.5
難点はxiが整数ではなくなってしまうことですが、以下の様にして「整数回の掛け算」に帰着させることが出来ます。
exp (xi)＝exp (int (2×x)/2)＝(√e)int (2×x)
√e(定数)の整数回掛け算で算出できました。
ただしこの場合、ループ回数が×2倍に増加します。折角テイラー級数の項数を減らしても、ループ回数が倍増してしまうと本末転倒になりかねないので良く吟味してください。今回はテーブルを引くことにするので、テーブルサイズが×2倍になるだけのことで、速度低下はありません。
xfの最大値が1→0.5に小さくなったことで、テイラー級数の項番n＝0〜4程度で必要な精度を確保することが出来、2項ほど削減(その分高速化)することが出来ました。
また、exp (x)の計算でx＜0の場合、浮動小数値の加減算で桁落ちが発生するらしいです。桁落ちによる演算精度の低下を回避するため、以下の関係式を使って計算します。
exp (x)＝1/exp (−x) …… x＜0の場合だけ、この関係式を使う
x＜0のケースなので、−x＞0です。exp (−x)をテイラー展開で計算したら、その逆数を取って完了となります。この方が、桁落ちによる誤差を最小化できるようです。
sin()関数と同様にして、シミュレーションで計算した数値をグラフ化したのが以下図です。
sin()関数グラフと同様、exp()関数値(青プロット、左対数軸)と、誤差率(赤線、1/1,000、右対数軸)としてあります。目標は誤差率(1/1,000)で1未満ですが、丁度1前後の値で分布しているようです。
自然対数ln()関数あるいは、log()関数は、前節のexp()関数の逆関数です。収束の遅かったexp()関数とは逆で、log()関数のテイラー展開(下式)の収束は速いです。
※個人的には、log (x)と言う書き方だと常用対数を思い浮かべてしまうので、ln (x)の方が馴染みが良いのですが、C言語のライブラリでは、log (x)が自然対数、log10 (x)が常用対数関数となっているので、これに倣って、本ページではlog (x)と表記します。
それでも、xの値を小さくすることは必要なので、以下の様にして計算します。
有効数字3桁の精度を確保するため、テイラー展開の項番n＝0〜2の計算が必要でした。シミュレーションで求めた数値をグラフ化した物が以下図です。
x＝1−αの付近で少々誤差率(1/1,000)が大きくなっていますが、それ以外の範囲では概ね目標である1未満を確保しています。
逆正接atan()関数あるいは、arctan()関数、tan-1()関数は、正接tan()関数の逆関数です。以下がatan()関数のテイラー展開です。
上図は、円周率の公式を流用しちゃいました。1/x ⇒ xと読み替えて見てください。これ、使ってみたのですが、exp()関数なんて比では無いほど、収束が遅いです。項番n＝0〜64程度まで算出してみたのですが、xの値によっては、全く収束できず、漸く有効数字2桁を確保できた程度でした。
遅すぎ。
そこで、atan()関数については、以下の様な連分数を使用することとします。
連分数については、以下より検索してみてください。
逆三角関数の連分数について検索する(別ウインドウが開きます)
今回、必要な精度を得るため、n＝4までの連分数を計算しました。
尚、この連分数についてはabs (x)≦1の範囲で定義されています。
一方、atan (x)の定義域はxの全範囲(−∞〜∞)なので、abs (x)＞1の範囲は以下の関係式を使って計算します。
※sgn (x)＝(x＜0)？-1：1；です。
シミュレーションで求めた数値をグラフ化すると、以下となります。
概ね、目標である誤差率(1/1,000)で1未満となっており、良好です。
前節までに作成した数学関数があると、その他の数学関数は以下の様にして求めることが可能です。
本ページで作成した数学関数を組み合わせることによって演算可能な関数を記しました。複数の関数を組み合わせることで、演算時間もかかるし精度も低下します。各々のテイラー展開があるので、本来は各関数毎に作成した方が、演算速度でも精度でも有利です。
ライブラリ化すると、それなりの動作検証や評価が必要なので、ここでは作成しません。将来、必要になった時に作ることにします（１＾＾；）。
C言語の数値関数で、使用頻度の高そうな関数を作っておこうと思います。
完成したMoscovium用数学関数ライブラリを以下にリンクしておきます。
数学関数ライブラリ：mcvm_lib_math.incl
収録した関数は以下の通りです。
組み合わせて算出する数学関数📐も含めて、基本的に必要そうなところは網羅してあるので、Moscovium＋H-FPUで何らかの数学関数演算が必要になった場合には、何とか凌げるでしょう。
