きっかけは、風速計を作った時でした。
「#22-6 風速計を作ろう（２）」
下図の様に脈動する起電圧を評価するため、実効値(RMS)を求めました。
「1.5V」等の浮動小数値を使わず、整数演算が可能な様に、わざわざmV単位で1,500mVなどと工夫してプログラムしたのに、実効値を求めるには、どうしても平方根演算(開平)が必要でした。
最後の最後でsqrt()関数を使わざるを得ず、浮動小数演算を使う羽目になってしまいました。
Arduino DUEで万歩計を作った時(「#21-8 万歩計を作ろう」)も、そうでした。加速度センサから読み出した、A＝(x，y，z)方向の加速度(ベクタ値)から、加速度の大きさ(スケイラ値)だけを知りたくて、本来は以下の計算が必要でした。
abs (A)＝√(x2＋y2＋z2)
(x，y，z)の各成分値は、MPU9250から読み出した16ビット符号付整数でしたので、適切な固定小数演算にて、(x2＋y2＋z2)までは求めましたが、あと一息、平方根を取るには浮動小数演算が必要でした。「万歩計」では開平せずに、そのまま(x2＋y2＋z2)を使って歩行を検知しました。
下図縦軸は本来、平方根を取る必要があるのですが、開平せずに使ったので、単位が重力加速度Gの二乗(G2)という見慣れない値になっています。
本来は、開平してabs (A)で判定すべきだったかと思われます。
浮動小数演算については、H-FPU(Half precision Floating Point Unit)の設計をした時に解説済みなので、繰り返しは避けます。簡単にまとめておくと、以下の通りです。
これをソフトウェアで実現するのは、コードサイズも大きくなるし演算速度も遅くなってしまいます。
C言語のリンカは賢くて、使っていないライブラリの関数は最終的なコードには含まれません。なので、整数演算しか使っていないプログラムでは、浮動小数演算ライブラリはリンクされません(期待値)。逆に、1回でも浮動小数演算を使ってしまうと、浮動小数演算ライブラリがリンクされます。これのサイズが結構大きいのです。
演算速度が遅いだけならまだしも、UNO等メモリの小さいマイコンでは、コードサイズに結構利いて来ます。
ところが、電子工作なんてやっていると、実効値(RMS)って結構使います。ベクタ値の大きさもしかりです(冒頭で言及した通り)。平方根の計算(開平)を、「浮動小数演算に頼らず整数計算できない物かな？」と、考えたのがスタートでした。
それでは、考えてみましょう。
二分探索のアルゴリズムは教科書に載っているので、以下より検索して見てください。
二分探索について調べる(別ウインドウが開きます)
※全然関連の無さそうなAD変換器の話だけど、こちらの解説も参考になるかも？
前節で考案したアルゴリズムをプログラムに落としていくのですが、いきなりArduinoでプログラムを作ってデバッグするのは、コンパイルして→書き込んで→シリアルモニタでコマンド打ち込んで、と、面倒です。そこで、まずはOlive+のスクリプトで整数平方根isqrt()関数のプログラムを作ってデバッグして行きます。
整数開平手続き(procedure)は以下の通りとなります(前節のアルゴリズム通り)。
procedure isqrt x
{
// int (sqrt (x))
var ll=0;
var rr=int (x/2);
var mm;
var y;
// exception
if (x≦3)
{
isqrt=( x≠0 );
return;
}
// binary search
while (ll≦rr)
{
mm=int ((rr + ll)/2);
y=int (x/mm);
// continue
if (y≦mm)
{
rr=mm - ( rr=mm );
}
else
{
ll=mm + ( ll=mm );
}
}
isqrt=mm;
minv isqrt=y;
}
二分探索のアルゴリズム説明は端折りましたが、簡単に上コードの動作を撫でておきます。
解y＝int (√x)の範囲を、下限値ll＝0から上限値rr＝x/2の間にあると設定し、その中間mm＝(rr＋ll)/2の数値(※)が、y＝int (√x)より大きいか、小さいか、判定して行きます。
ll〜rrの範囲が1/2づつ狭まって行って、最後にll＞rrと大小が逆転すると、探索終了となります。
普通、二分探索では、求める解yが0〜xの間にあると仮定してスタートします。この場合の探索回数はlog2 (x)回です。この探索回数の分だけ、除算を実行することになるので、log2 (x)がこのアルゴリズムのコスト💰となります。
つまり、このアルゴリズムでは、log2 (x)に比例して演算時間が増加します。
※本来、mmの計算は、mm＝ll＋(rr−ll)/2として求めないとならないみたいです。理由は、rr＋llの加算結果がMAX_INTを超えた場合に誤動作するからです。今回のアルゴリズムでは、rr＋llがMAX_INTを超えることは無いし、減算1回分遅くなるので、敢えてmm＝(rr＋ll)/2の方を使っています。
二分探索において、探索範囲0〜xに対して最初に探索する値mmは本来、mm＝(rr＋ll)/2＝x/2なのだけど、求める解√xがmm＝x/2より大きくなるのは、x＝1，2，3の3値だけです。
x＝1の時、√x＝1 ＞x/2＝0.5
x＝2の時、√x＝1.414…＞x/2＝1
x＝3の時、√x＝1.732…＞x/2＝1.5
…… これより大きなxでは、解sqrt (x)≦x/2の範囲にある。
x＝4の時、√x＝2 ≦x/2＝2
x＝5の時、√x＝2.236…＜x/2＝2.5
：
従って、x≦3の場合を例外処理して、探索範囲を0〜x/2からスタートすると、探索回数を1回減らすことが可能です。この時の探索回数＝log2 (x)−1回です。
ちょっと厄介なのは、普通にデータベースから条件に合致したデータを探し出す二分探索と違って、探すデータy(＝√x)そのものが未知な状態で、完全一致するデータが見つからない場合でも、条件に合致した最大の整数を探さなければなりません。
一旦テストしたmmのデータは、まだ解である可能性が残っているので、除外できません。
そんなこともあってか普通に二分探索すると、y＝int (√x)が見つかる時と、y＝int (√x)−1が見つかる時があるみたいです。そこで少し工夫して、確実にy＝int (√x)が見つかる様にプログラムしてみました。
新しいll，rrを決定するプロセスで、ll＝mmとrr＝mmの場合分けをして、rr＝llとなるまで探索を続ける部分等が、普通の二分探索では不要な部分です。この処理によって、探索回数が＋1回増加しています。
全体をまとめてみると、探索範囲を0〜x/2に絞ったことで、探索回数−1回の効果と、変則的な二分探索に対応したことによる探索回数＋1回の効果が打ち消し合って、探索回数＝log2 (x)−1＋1＝log2 (x)となります。
例えば整数xに対して開平する時、二分探索(従って除算)の回数はlog2 (x)回です(前節)。一方、逐次探索(シーケンシャルサーチ)では、√x回となります。
下図の通り、log2 (x)−√xのグラフを描いてみると、x＞16では二分探索に軍配が上がります。
x＜16未満では逐次探索の方が速いけど、16ビットとか32ビットとか言う数値の話をしている中の16未満なので、一般的に二分探索の方が速いと考えて差し支えありません。
二分探索で答えを探すと、探索回数、従って除算回数がlog2 (x)回となるので、それが整数開平のコストです。扱う数値の最大値が16ビットか32ビットか？ではなく、開平しようとするxの大きさが、演算速度に関係することに留意しておいてください。数値xが小さければ小さいほど、高速に開平できます。
さて、Arduinoの様にC言語などの環境で浮動小数演算が可能な場合、今回のアルゴリズム(整数演算のみ)で開平する場合と、浮動小数演算(sqrt()関数)で開平する場合、結局どっちが速いのだろう？
今回の整数開平アルゴリズム：
浮動小数開平アルゴリズム：
数値xの値が十分小さい間は、文句なしに整数開平アルゴリズムの方が速いのですが、xが大きくなってくると整数除算回数が増えて、だんだん遅くなっていきます。
結論は一概には出せません。
整数開平isqrtに関しては、整数除算が主たるコストの中身なので、ハードウェア除算器を搭載しているArduino DUE(Cortex-M3)やMoscovium(＋DIVC32)等のMCUでは高速に演算可能ですが、非搭載のUNO(ATmega328)等では、シフト除算ルーチンを使うことになって低速です。
浮動小数演算(sqrt()関数)にしても、FPU持ちのESP32(Xtensa LX6)やMoscovium(＋H-FPU)は高速ですが、DUE(Cortex-M3)やUNO(ATmega328)では、浮動小数演算ライブラリ(ソフトウェア)がリンクされて低速です。
浮動小数sqrt()関数と比べて、どの程度の数値x(大きさ)まで、整数開平isqrtの方が高速なのか、については環境によって異なるとしか言いようがないです。演算速度が最重要な場合は、ご使用になる環境でベンチマークして決める事になると思われます。
※浮動小数演算ライブラリをリンクする環境と比べたら、恐らく、整数開平isqrtの方が、コードサイズは小さくなると思われます。コードサイズの方が重要な場合は間違いなくソリューションになるでしょう。
それでは、作成したisqrt手続きを使って、動作をテストしてみましょう。完成したOlive+スクリプトを以下にリンクします。
完成スクリプト：isqrt.olv
※スクリプト(.olv)の実行には拙作Olive+をWindowsにインストールしてください。
スクリプト(isqrt.olv)を実行すると、以下の様なウインドウが開きます。
変数numに、数値xの値を入力してOKボタンを押すと、isqrt (x)を計算して表示します。以下の様に、期待値としてsqrt (x)の値と、二分探索の様子もレポートされます。
sqrt (1234)=35.1283, isqrt (1234)=35
cnt=1 ll=0 rr=617 y=4 mm=308
cnt=2 ll=0 rr=308 y=8 mm=154
cnt=3 ll=0 rr=154 y=16 mm=77
cnt=4 ll=0 rr=77 y=32 mm=38
cnt=5 ll=0 rr=38 y=64 mm=19
cnt=6 ll=19 rr=38 y=44 mm=28
cnt=7 ll=28 rr=38 y=37 mm=33
cnt=8 ll=33 rr=38 y=35 mm=35
cnt=9 ll=33 rr=35 y=36 mm=34
cnt=10 ll=34 rr=35 y=36 mm=34
cnt=11 ll=35 rr=35 y=35 mm=35
x＝1,234の平方根を計算させました。平方根sqrt (x)≒35.1283に対して、整数の平方根isqrt (x)＝35と、合致しています。
その下に並んだリストは、探索の途中経過です。log2 (1,234)≒10.269なので、11回探索しました(小数点以下は切り上げ)。ll＝0，rr＝x/2からスタートして、だんだん範囲が狭まって行き、最後にll＝rrへ収束する様子が見て取れます。
OKボタンを押すと、またnumの入力を求めるウインドウが開きます。このウインドウでCancelボタンを押すと、スクリプトを終了することができます。
また、numを入力するウインドウでnumに負の数を指定すると、abs (num)回ランダムに発生した数値xに対して、int (sqrt (x))とisqrt (x)を自動計算して、期待値比較します。
例）num＝−1,000を指定した時、1,000個の乱数を開平して、結果を「isqrt_out.log」ファイルに出力します。
以下、ログファイルの内容抜粋です。角かっこ[ ]内が期待値です。期待値と不一致すると、先頭の「○」が「×」に変わります。
○isqrt (48978)=221 [221]
○isqrt (37684)=194 [194]
○isqrt (40408)=201 [201]
：
○isqrt (38313)=195 [195]
○isqrt (8134)=90 [90]
○isqrt (23752)=154 [154]
errors=0/1000
すべてのログを確認しなくても、最後の「errors=0/1000」でエラー(期待値不一致)の数を確認できます。毎回乱数を作り直すので、繰り返す度に異なる数値の開平を検証できます。気の済むまで、自動検証させてみてください。
当方で検証した範囲では、エラー無くisqrt (x)が計算できていました。
前節にて、Olive+のスクリプト(isqrt.olv)版isqrt()関数のデバッグを済ませたので、これをC言語に移植してArduinoでも試してみましょう。以下に完成したプログラムをリンクします。
完成プログラム：test_isqrt.ino
コンパイルしてArduinoへ書き込むと、シリアルモニタから以下のコマンドを受け付けています。
シリアルモニタコマンド仕様：
Arduino UNOを使って検証してみましょう。
※プログラム(test_isqrt.ino)は汎用なので、他のArduinoでも動作するはずです。Arduinoのunsigned intサイズに従ったサイズの整数開平を検証してみてください。例えばArduino DUEで実行すると、32ビット符号無整数の開平が可能なはずです。
実験サンプル：【 】内はシリアルモニタからの送信データ
【isqrt,32768】
OK: isqrt (32768)=181 [181] O
整数x＝32,768の整数平方根を求めました。角カッコ[ ]内は、int (sqrt (x))を計算した期待値です。解は181で正解です。
【self,100】
OK+ isqrt (39765)=199 [199]
OK+ isqrt (4638)=68 [68]
OK+ isqrt (34222)=184 [184]
：
OK+ isqrt (33609)=183 [183]
OK+ isqrt (27997)=167 [167]
OK+ isqrt (4844)=69 [69]
OK: self error=0/100
セルフテストを実行しました。ランダムに整数xを100回生成しながら、期待値との比較を実行して結果を表示しました。期待値と不一致すると、先頭の「OK+」が「ERR+」に変化します。最後に表示した要約「error=0/100」によって、期待値との不一致が0だったことが判ります。
同様に、数1,000個の期待値比較を実施したところ、不一致はありませんでした。
浮動小数演算するsqrt()関数のコードサイズを確認しておきます。プログラム先頭付近に「#define USE_SQRT」と言うマクロが定義されているので、これをコメントしてコンパイルすると、sqrt()関数を使わずにコンパイルします。
※期待値生成に使うsqrt()関数を使わないでコンパイルするので、プログラム自体は正常動作しなくなります。コードサイズ確認用のオプション機能とご理解ください。
「#define USE_SQRT」の有り／無しでコンパイルした結果を以下表にまとめます。
| マイコンボード | USE_SQRT | 差分 | |
| 有り | 無し | ||
| Arduino UNO | 5,338バイト | 4,812バイト | 526バイト |
| Arduino DUE | 32,420バイト | 31,764バイト | 656バイト |
| STM32 | 38,144バイト | 37,496バイト | 648バイト |
| ESP32 | 214,965バイト | 214,365バイト | 600バイト |
コンパイル後に表示されるメッセージから、書き込みに必要なROM(Flash)サイズを拾って表にしました。最後の「差分」がsqrt()関数に要したコードサイズということになります。
何か、ちょっとイメージと違うな……。みんな揃って差分＝500〜600バイトって、どうしてだろう？ESP32はFPU持ちなので、浮動小数演算ルーチンは不要なはず。コンパイル後のアセンブラソースまでは確認していないので、詳細不明です。
※もしかして、float演算を使っていなくても、浮動小数演算ライブラリがリンクされている？？
最後に、オリジナル16ビットCPUコアMoscoviumにも移植しておきます。以下にMoscoviumのアセンブラに移植したプログラムをリンクします。
完成プログラム：mcoc_isqrt.asm
Moscovium版のプログラムは、iverilogを使ったシミュレーション確認のみとしました。アセンブルしたら、適切なMoscoviumマイコンのデザインを使ってシミュレーションして見てください。実行には、除算器DIVC32が必要なので、型名一覧表からDIVC32を搭載したMoscoviumか、Nihonium搭載エディションを使ってください。
実験サンプル：【 】内はCygwinシェルウインドウからの実行コマンド
【iverilog -f test.lst ; if ( ! $status ) vvp a.out】
：
isqrt ( 241)= 15 [ 15]o
isqrt ( 824)= 28 [ 28]o
isqrt ( 2799)= 52 [ 52]o
isqrt ( 6609)= 81 [ 81]o
isqrt (14945)= 122 [ 122]o
isqrt (27375)= 165 [ 165]o
isqrt (42921)= 207 [ 207]o
isqrt (52498)= 229 [ 229]o
isqrt (64169)= 253 [ 253]o
Test Pass
CPI: 72182/ 48455( 0)= 1489/1000
finish: 2888880 ns
整数開平isqrtの結果と、角かっこ[ ]内は期待値です。シミュレーションログから、一部抜粋しました。
行末に「o」の文字が見えますが、期待値比較が合致した印です。不一致した場合、代わりに「x」が表示されます。実行した範囲では、期待値との不一致はありませんでした。問題なさそうです。
isqrtサブルーチン部分のコードサイズは、36ワードでした。ただし、除算命令(divu命令)は使っています。divu命令を使わずに、シフト除算ルーチンを用意すると、もう少しコードサイズが大きくなります。それでも、浮動小数演算ルーチンを作ってsqrt()関数まで実装することを考えたら、ずっと小さいコードサイズです。
isqrtサブルーチンはレジスタサイズに無関係なので、Nihoniumを使って32ビットに拡張(NHビット＝1)すると、そのままのコードで32ビット符号無整数の開平が可能です。そこは、C言語のintみたいで良いでしょ？
整数計算だけで平方根を計算できる可能性を示してみました。何かで活用できると良いのだけど。ちなみにisqrt (5)＝2だけど、×1002倍して、isqrt (50,000)＝223より、√5≒2.23まで求めることができたりします。それでも「富士山麓(√5≒2.236)」にも及ばないのですけどね。
えっと、√xを求めるのに、以下のように考えました。
ですが、よくよく考えてみたら、yを移行する理由がありませんでしたね。
y＝x÷yと変形した場合と大差ないようにも思えますが、本文中で提案した方法(y＝x÷y)だと、割り算が必要なのに対して、今回の方法(y2＝x)では割り算を使わないことが大きな違いです。
一般に、コンピュータ除算は遅いので、除算を使わず、乗算で決着した方が、演算速度は高速となります。
ただし、注意点もあります。
y＝x÷yの方法だと、演算結果x÷yは、必ずxよりも小さくなります。しかし、今回の方法(y2＝x)だと、演算結果y2がxの値を超える場合があります。ただ超えるだけなら良いのですが、もし、超えた結果、整数型(uint16やuint32等)の最大値(MAX_INT)を超えてしまった場合、処理を誤る恐れがあります。
例えば、uint16で処理する場合で、x＝65,535の開平は今回の方法(y2＝x)では求めることができません。一方、本文中の方法(y＝x÷y)ならば、求めることが可能です。
速度を重視する場合、今回の方法(y2＝x)は有効ですが、一般的にすべての数値範囲で整数開平したい場合、本文中の方法(y＝x÷y)しか解となり得ません。
余談でした。
