公開日：2021/07/03

16th Season

#16-9　数学関数を整備しよう

　C言語なんかにどっぷり浸かっていると、sin()関数とかexp()関数とか数学関数は普通に使えるし、なにも困ることは無いわけです。でもそれって、何方かが数学関数ライブラリを整備してくれているからこそ、なので、オリジナル16ビットCPUなんて作って見たりすると、はたと困ってしまうのです。

　そこで今回は、「一般的な数学関数くらい自分で作らなきゃね」をテーマとします。

　使用するのは、半精度浮動小数点演算ユニット(HFPU)を搭載した、MCOC115FP0440です。デザインデータは以下リンクよりダウンロード可能です。

テイラー展開

　コンピュータで数学関数を作る、と言ってまず思いつくのはテイラー展開と言う方法です。

• テイラー展開の代わりにテイラー級数などと呼ぶ場合もあります。
• 0を中心としたテイラー級数は、マクローリン級数とも呼ばれます。
• グレゴリー級数と呼ばれることもあるみたいです。
• 全部、同様な級数に言及した数学者の名前らしいです。
• 何らか、呼び分ける指標があるのかどうかについては不明です。
• 本ページでは、テイラー展開(級数)と呼ぶことにします。

　円周率を計算した時に使ったのは、arctan (1/x)のグレゴリー級数でした(↓)。



　上のまとめに従うと、テイラー級数の一種だった、と言うことになります。


　数学関数の値を求める時は、上の様な級数に値を代入して、関数値を求めます。
　ところがこれらの級数は、項番nの値(上の数式だと項番k)が無限大まで続くのです。コンピュータでは無限大と言う数値を扱うことが出来ないので、適切なところで打ち切るしかありません。では、どのくらいの値で打ち切るのか？と言うと、「必要な演算精度が得られたところ」で打ち切る必要があります。
　今回は半精度浮動小数点数half型がターゲットで、10進数換算の有効数字は3桁程度しかありません。たったの3桁で良いので、さぞかし少ない項数で十分な精度が得られるのではないか？と期待します。
　ところが、実際にはそう上手くは行かない様です。後で関数毎に見てみましょう。

　一般的には、テイラー級数第n項の数値がどれだけ小さいか？で、収束速度が異なります。収束の遅い関数は、高次項まで計算する必要があるので計算が遅くなります。

用意する関数

　まず、今回作成する関数ライブラリの目標仕様を以下の通りとします。

• 最低限度必要な関数を揃える。
• 具体的用途も無い状態で、かつ、有効数字3桁程度のhalf型数学関数であるため、片っ端から全部そろえたところで労力の無駄遣い。
• 作るからには、ある程度評価も必要だし、デバッグも要るので、その手間を幾分か減らす。
• 今回作成する数学関数：　正弦sin()関数、自然対数log()関数、指数exp()関数、逆正接atan()関数
• その他、数値演算関数：　round()関数、floor()関数、ceil()関数、abs()関数
• 演算速度は不問とする。
• 演算速度を求める場合、そもそも、テイラー級数ではなく、別な高速アルゴリズムが用いられている模様。
• 速度が必要なら、Cのライブラリ関数のインプリを真似したら良い。GPLあるいはLGPLなどで利用可能なはず。
• 筆者の様な素人ではなく、その道のプロたちが最適なアルゴリズムでインプリしていると考えられる。
• ここでは、テイラー級数を基本とする。
• 演算精度は可能な限り確保するが、精度についても不問とする。
• と言っても、作成する関数については、可能な限り精度を追求する。
• 後ほど触れるが、例えば、平方根sqrt (x)＝exp (log (x)/2)で計算することになる。
• 各々、exp()関数とlog()関数の精度は追及するが、それら2つを組み合わせて演算することになるsqrt()関数の精度は、exp()関数、log()関数のそれよりも低下するのは必定。
• それについては不問に付す、と言う事。
• ここでは、10進数換算有効数字3桁程度のhalf型で、「値を演算できれば良い」と割り切る。

　おそらく、演算速度や演算精度を必要とする用途では、根本的に対処すべき問題が異なるので、このページには到達しないであろうと推察します。

三角関数

　早速、三角関数から作って行きます。三角関数の中で今回実際に作成するのは正弦sin()関数だけです。sin()関数が有れば、余弦cos()関数と正接tan()関数は以下式で求められます。





　sin()関数のテイラー展開は以下式となります。



　第n項目の数値は、分母に(2n＋1)！←階乗という比較的大きな数値が来るため、小さくなりそうな気がしますが、問題は分子にあるx^2n＋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によるシミュレーションで計算させてみました(HFPUが必要なので、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＝−∞〜∞)で定義されていて、しかも級数の分子側にxⁿが存在するため、収束の遅い関数となっています。つまり、必要な精度を確保するためには大きな項番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≦xfⁿ＜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の値を小さくすることは必要なので、以下の様にして計算します。

• 定義域はx＝0〜∞、ただしx＞0。
• x＝xf×2^xeという、half型の各要素、仮数部xf，指数部xeに分解して計算する。
• half型では、1≦xf＜2、−14≦xe≦15。
• log (x)＝log (xf×2^xe)＝log (xf)＋xe×log (2)
• log (2)は定数。そこに2の指数xeを掛け算して足し込む。
• log (xf)に関しては、1≦xf＜2と、範囲も狭くて値も小さいので高速に収束できる。

　有効数字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未満となっており、良好です。

その他の数学関数

　前節までに作成した数学関数があると、その他の数学関数は以下の様にして求めることが可能です。

• べき乗関数pow (x，y)：
• 
• z＝x^yと置いて両辺の自然対数を取ると、log (z)＝y×log (x)。
• 更に両辺ネイピア数のべき乗を取って、exp (log (z))＝z＝exp (y×log (x))。
• 平方根関数sqrt (x)：
• 
• √x＝sqrt (x)＝x^1/2＝exp (log (x)/2)。
• 一般に割り算より掛け算の方が速いので、「2で割る」より「0.5を掛ける」方が速い。
• ※Moscovium＋HFPUでも、hmul命令＝3〜4[サイクル]、hdiv命令＝6〜7[サイクル]と、掛け算の方が速い。
• 余弦関数cos (x)：
• 
• 正接関数tan (x)：
• 
• 常用対数関数log₁₀ (x)：
• 
• 自然対数として求めた値log (x)をlog (10)で割り算する。log (10)は定数。
• 一般に割り算より掛け算の方が速いので、1/log (10)を掛けた方が速い。1/log (10)も定数。
• 逆正弦関数asin (x)：
• 
• 逆余弦関数acos (x)：
• 

　本ページで作成した数学関数を組み合わせることによって演算可能な関数を記しました。複数の関数を組み合わせることで、演算時間もかかるし精度も低下します。各々のテイラー展開があるので、本来は各関数毎に作成した方が、演算速度でも精度でも有利です。
　ライブラリ化すると、それなりの動作検証や評価が必要なので、ここでは作成しません。将来、必要になった時に作ることにします（１＾＾；）。

その他の数値演算関数

　C言語の数値関数で、使用頻度の高そうな関数を作っておこうと思います。

• floor (x)：　xの値を超えない最大の整数
• half型→uint16型変換を実行するhuint命令がありますが、huint命令では、half型Xに対して絶対値abs (X)の小数点以下を切り捨てた整数値が返ります。
• つまり、符号付浮動小数点数値として見ると、ゼロ0方向への丸め(正値は切り捨て、負値は切り上げ＋絶対値)です。
• これはfloor (x)ではありません。floor (x)とするには、以下の様な手順となります。
• // r0=floor (r1)
  huint　r2,r1　　// 一旦整数化して
  hhalf　r0,r2　　// half型に戻す
  hmvsg　r0,r1　　// 符号も復元
  hsub　r1,r0　　// 元の数値と比べて
  mov　r1,r1
  bpl　pcnt+6
  addi　r2,1　　// 大きかったら、r1の値を超えない様に調整
  hhalf　r0,r2
  hmvsg　r0,r1
  
• 上の処理でr1＜r0となるケースは、元の数値r1＜0だったケースなので、＋1で調整します。
• r1＝HALF_MAX_N＝−65,504の1つ上の数値は−65,472ですが、r2＝−65,504＋1＝−65,503をhhalfした時点で、r1＝half (65,472)になるため、ループは不要で、上のコードで充足します。
• ceil (x)：　xの値を下回らない最小の整数
• floor (x)の応用で求められます。
• // r0=ceil (r1)
  huint　r2,r1
  hhalf　r0,r2
  hmvsg　r0,r1
  hsub　r0,r1　　// floor()とは減算方向が逆
  mov　r0,r0
  bpl　pcnt+2
  addi　r2,1　　// 小さかったら、r1の値を下回らない様に調整
  hhalf　r0,r2
  hmvsg　r0,r1
  
• 上の処理でr0＜r1となるケースは、元の数値r1＞0だったケースなので、＋1で調整します。
• abs (x)：　xの絶対値
• 符号ビットを0にすれば良いだけなので、色々な方法で実装可能です。今回は以下としました。
• // r0=abs (r1)
  ldbiu　r0,0
  hmvsg　r1,r0
  mov　r0,r1
  
• round (x)：　xの小数点以下を四捨五入
• C言語などで使う場合には、単にfloor (x＋0.5)として算出してしまいます。
• ですが、本当はこれではダメです。負の整数値xに対してfloor (x＋0.5)を実行すると、負数xに対して0.5を加えたことで減算が発生します。桁落ち分のボロウによって、例えば−65,504＋0.5 → −65,472となってしまいました。
• これは浮動小数点演算の特徴なので、程度の差こそあれ、float型やdouble型でも発生する可能性があります(未確認)。
• そこでまず、xが整数かどうかを調べ、整数だった場合にはxをそのまま返すように例外処理を追加しました。
• // integer?
  hfrac　r2,r1
  mov　r2,r2
  bne　pcnt+4
  mov　r0,r1
  rtn
  
• hfrac命令で整数かどうかを確かめて、整数だったらそのまま返します。
• half型の絶対値最小値HALF_MIN≒6.10×10^-5でも、hfrac命令によって非ゼロ値が返るので、上のコードで整数かどうかの判定が可能です。

完成した数学関数ライブラリ

　完成したMoscovium用数学関数ライブラリを以下にリンクしておきます。

　　数学関数ライブラリ：mcvm_lib_math.incl

　収録した関数は以下の通りです。

• // math_atan: half r0=atan (half r1), disturbed r1~r3
• // math_exp: half r0=exp (half r1), disturbed r1~r5
• // math_log: half r0=log (half r1), disturbed r1~r4
• // math_sin: half r0=sin (half r1), disturbed r1~r4
• // math_abs: half r0=abs (half r1), disturbed r1
• // math_ceil: half r0=ceil (half r1) ; disturbed r2
• // math_floor: half r0=floor (half r1), disturbed r1~r3
• // math_round: half r0=round (half r1), disturbed r1,r2

　組み合わせて算出する数学関数📐も含めて、基本的に必要そうなところは網羅してあるので、Moscovium＋HFPUで何らかの数学関数演算が必要になった場合には、何とか凌げるでしょう。