LabVIEW による浮動小数点表現参考プログラム
(最終改訂 2018/12/12)
書庫ファイル内容
LV_IEEE754.zipには本稿(拡張子.pdf)、本稿掲載図(フォルダFigures:拡張子.png)、LabVIEW による浮動小数点表現参考プログラム3本および有限なデータ長での演算と関連の深い数値積分参考プ
ログラム 2 本のLabVIEW2017 版(情報科学研究教育センター(以下「センター」と表記)での実行
にはランタイムエンジン不要)および LabVIEW2013 版(ランタイムエンジンのダウンロードに National Instruments社(以下「NI」と表記)のアカウント登録が不要)の実行プログラム(拡張子.exe)、 設定ファイル(拡張子.ini)、実行時ライブラリ(サブフォルダ data:拡張子.dll)がフォルダLV2017 およびフォルダ LV2013 に格納されている。実行プログラムの場所に同名の設定ファイルおよび data フォルダを置いて開く(p.12備考参照)。
1. LV_IEEE754DP
LabVIEW のフロントパネル上に数値制御器、数値表示器を作成するときの既定のデータ型 DBL
(DouBLe-precision)の内容はIEEE754規格で表現された64b浮動小数点数である。歴史的に「倍精 度浮動小数点数」と称するが現在はこれが標準であり、むしろ科学技術計算では範囲も精度も不十分な 32bの「単精度」(LabVIEWのSGL:SinGLe-precision)が「『短』精度」と呼ばれるべきものである。
LV_IEEE754DPはLabVIEWの数式評価VIにより数式文字列で入力した1変数関数(註参照)を 評価して関数値の64ビット内部表現を表示し、また色分けされた64個のブールボタン(符号:赤、指 数:黄、仮数:緑)または16個の16進数字メニューリング(一方の操作は他方に連動)で内部表現を 指定して相互転送(関数値⇒内部表現設定制御器、内部表現⇒変数値設定制御器)することでIEEE754 表現を対話的・視覚的に確認し、浮動小数点表現についての理解を深めることを目的とするものである。
註:ここでは変数としてxを用いているが、1変数である限り変数名は任意の英小文字1字(これに数 字1字が続いてもよい:備考参照)が使える。使用できる関数は、例えば次のURLにリストがある。
http://www.mathmachines.net/Construction/SSSP/LabVIEWFormulaStrings.pdf
図1に起動時のフロントパネル、図2(p.2)~図5(p.5)に制御器と表示器の内容を示す。
図1 LV_IEEE754DP起動時の既定の画面(6!=720を表示している)
数式文字列の設定
関数を指定する数式文字列は文字列制御器(キャプション ”arithmetic expression (single variable at most)”:起動時既定は”gamma(x+1)”)に直接入力するかまたは、① ”MEMORY” ボタンON(起 動時既定)の状態でメニューリングの項目を指定して”↓RCL”をクリックしてプリセットメモリから転 送、もしくは②”READ”(註参照)で読込んだテキストファイルの各行で項目を構成したメニューリン グ(起動時既定はLV_IEEE754.zip同梱のsample.txt を読込んで構成した内容)の数式が登録され た項目を指定して” ⇨ ENTER”をクリックして転送する方法がある(図2左)。
註:フロントパネル左上の ”READ”ボタンをクリックして関数、16 進データを記述したテキストファ イルを読込むことができる。VIのプロパティで「自動エラー処理を有効」のチェックを外して実行ファ イルをビルドしており、ファイルダイアログをキャンセルしてもエラーとはならない(ログのファイル 出力についても同様である)。テキストファイルには、1行に数式文字列フォーマットによる1個の数式 または$で始まる1個の16進数値(左詰:データが16桁を超える場合は先頭から16桁が使われ16桁 に満たない場合は右に0が補われる)を記述し、何れの形式も%に続けて注釈を付加することができる。
数式として評価できない文字列が入力された場合、文字列制御器の上に赤 LED(ラベル”error”)を 点灯表示し、フロントパネル左下の数式評価結果のまとめと”PRINT”ボタン(「ブールパネルの設定・
転送とログ出力」の項を参照)はグレーアウトされる。
文字列制御器右の押しボタンSW ”identity function”は恒等関数の指定で、ON(文字列制御器の背景 色が緑に変る)の状態では文字列制御器の書込みは禁止され、”↓RCL”ボタンとメニューリング項目が 数式の場合に” ⇨ ENTER”ボタンが無効化され、他の関数に変更することはできない(図2右)。
図2 関数指定に関係する制御器
起動時のプリセットメモリには、「表現可能な最大の数は?」、「正規化された絶対値最小の数は?」、
「階乗はいくつまで『整数』として表せる?」、「2から√を何回繰返したら1になる?」「π、eなどの 身近な定数の内部表現は?」等の浮動小数点表現に関する素朴な疑問に答える次の 22 個がプリセット されている。
プリセットメモリの内容は変更可能で、”↑STO”の押下で数式文字列がメニューリングの表示項目に 上書き(押しボタンSW ”append” OFF)または項目の最後尾に追加(押しボタンSW “append”ON:
起動時既定)される。
0/0 : 非数(NaN)となる例 1/0 : 正の無限大(Inf)となる例
2^53-1 : 整数として正確に表現される最大の数 (2^53-1)*2^971 : 表現できる最大の数
2^-1022 : 53ビット精度で表現できる最小の正数 2^-1074 : 表現できる最小の正数
exp(x) : 指数関数 pi(x) : πの倍数 ln(x) : 自然対数
log(x) : 常用対数 sqrt(x) : 正の平方根
100^0.2 : 星の1等級の光度比(4dB)
gamma(x+1) : 非負の整数のとき階乗
exp(-x^2/2)/sqrt(pi(2)) : 正規分布の確率密度 exp(-exp(x)) : 二重指数関数
exp(exp(x)) : 二重指数関数
exp(ln(2)*(2^(-x))) : 繰返し平方根
(x+(2/x))/2 : Newton-Raphson法による2の平方根 2^x : 2の冪乗
10^x : 10の冪乗 16^x : 16の冪乗 1/x : 逆数
変数値の設定
関数に設定した数式文字列が変数を含むときフロントパネルの関数指定制御器の左下に変数値の表 示器と帰還設定の制御器(図3下)が表示され、押しボタンSW ”feedback mode”(起動時既定OFF)
が OFF のときフロントパネルの関数指定制御器の上に変数値設定用の制御器(図 3 上)が表示され
る。”feedback mode” OFFのとき変数値は数値制御器または文字列制御器で設定し、制御器を垂直スラ
イドSW ”formula”(起動時既定OFF)で選択する。
図3 関数指定の数式が変数を含むときの変数値指定
数値スライド(ラベル”numeric”)のスケール両端の値は数値スライド上の数値ボックス(ラベ ル”left”、”right”)で指定し、左端より右端の値が小さければスケールは左右反転、両端共に正の場合は 押しボタンSW ”logarithmic”が表示され、ONで対数スケールとなる。数値スライド右の数値ボックス には数値スライドでは設定できないスケール範囲外の値(例:-inf)も入力可能である。
文字列制御器(ラベル”formula”)には数値文字列(数式文字列では指数記号で小文字は使えない:
備考参照)、変数を含まない数式文字列(例:sqrt(2))または64b内部表現を”$”に続く16個以内の 16 進数字(テキストファイルからの読込みと同様の左詰で下位の 0 は省略可能で、例えば$C の値は
$C000000000000000の-2を表す)で直接入力する。
数値スライド右上の押しボタンSW ”force to integer”をON(起動時既定)にすると垂直スライドSW により選択された制御器の設定値を直近の整数に丸めた値が使用される。
フロントパネル右上の押しボタンSW “⇦ argument”(起動時既定OFF)をONにすると(関数評価 値をブールパネルに転送する押しボタンSW ”⇧ Boole / hex”とは排他的で一方がONのとき他方は無効 化される:p.5図5参照)ONにした直後とブールパネル設定値に変化があったときSW ”formula”で選 択された制御器に(”formula”ONのときは$xxx…xxの形で)転送される(ブールパネル設定値に変更 が無い場合数値制御器または文字列制御器の操作は有効)。
"feedback mode” ONの帰還モードでは変数値の設定に関係するフロントパネル上段の制御器は非表 示となり、"feedback mode”ボタンの右に表示される” ↷ feedback”ボタンをクリックして放したときク リック前の関数値が変数値に渡される(p.3図3下左)。ただし関数評価で使用する変数値として渡すの みで非表示となっている変数値設定の制御器の内容は変化しない。
メニューリングで選択した項目が16進データの場合にはモードに無関係に"ENTER"のクリックでブ ールパネルに送られるが、帰還モードでは押しボタンSW "set argument"SWが表示され、ON(起動 時既定)のとき"ENTER"のクリックで当該データが変数値にも渡される(p.3図3下右)。
関数値の表示
評価された関数値はフロントパネル中央部に64b内部表現を8×8の2Dブール配列(row-major order、
左上が符号を示すMSB)、フロントパネル下部に64b全体、11b指数部、52b仮数部(註参照)をそれ ぞれのビット位置を縦に揃えた1D ブール配列(左が上位)で表示し、10進数値、16 進数値、仮数値 等をその周囲に表示する(図4左下)。メモリが潤沢で64b浮動小数点演算回路が標準である現在、敢 えて32b表現を使用する必然性は失われているが、32bでは不十分であることを示すため、2Dブール 配列の下に32b に短縮した場合の値とあふれ情報を表示している(例を図4右上に示す)。各LEDは 精度短縮による指数あふれ、非正規化の情報であり、64b表現で既に無限大、非数、零であった場合に は点灯させていない。
註:52b仮数左の赤LED(ラベル”implicit MSB”)はデータ中には無く指数部(1023げたばき)が1~2046 の正規化数(仮数 53b)で 1が仮定される暗黙の仮数 MSB(隠しビット)である。無限大と非数では 意味を持たないが本VIではIntel FPUでの表現に揃えて±InfとNaNの場合にも1を表示している。
図4 関数値の表示内容(左下)と32ビットに短縮時の指数あふれと非正規化の例(右上)
関数指定の数式文字列で ”gamma(x+1)”(xが非負整数のときxの階乗)、変数値入力で数値制御器 を選択し、”force to integer” をONとしたLV_IEEE754DP起動時既定(p.1図1参照)は、整数とし て表現できる階乗の最大値を求める設定である。
数値スライドを操作して18! = 6402373705728000までは”non-zero integer” LEDが点灯し、整数と して正確に表現できることを確認できる。既定の変数値6のとき8×8ビット2D配列上の数値表示器
(ラベル”numeric”)には 719.999…99と表示されているが、これは書式編集の2進10 進変換誤差に よるものでビット内容から実際の内部データは正しく6!=720の整数値を表していることが分る。
指数値が0~52の範囲にある(「浮動」する真の1の位が暗黙のMSBを含めた仮数内にある)とき、
仮数部LED配列の右に”non-zero integer” LEDと”blink true binary place”ボタンが表示され、1の位 が暗黙の MSB を含めた仮数部内にあり、かつ仮数部内の 1 未満のビットが全て 0 のとき”non-zero integer” が点灯する。”blink true binary place”をONにすると暗黙MSBを含む仮数部の1の位のLED が点滅する。図4は1の位(偶数であり値は0)が点灯しているループのキャプチャ画面である。
ブールパネルの設定・転送とログ出力
フロントパネル右上のブールパネルはブールボタン(ビットが持つ意味により色分けされた 64 個:
ボタン上の数字はLSBを0とするビット位置)または16進数字メニューリング(8×2の2D配列:上 が上位、横2個で1バイト)で内部表現を1ビットまたは4ビット単位で指定するもので、転送指定の 押しボタンSWの設定により、設定した内部表現が表す値を関数値(内容はp.4図4左下)として表示
【恒等関数指定 ”identity function” および転送指定 “⇦ argument”を ON とする】すること、逆に関 数値をブールパネルに表示【転送指定 ”⇧ Boole / hex”をONとする】することができる(図5左)。転 送指定の2個の押しボタンSWの機能は排他的であるため、一方がONのとき他方を無効化しており、
同時にONにすることはできない。
図5 ブールパネル(左)とログ出力(右)
フロントパネル左下の文字列表示器(ラベル”current value”)には関数値のまとめが#で始まる次の 形式で表示される。(i) 変数を含まない数式では、数式、コロン、16進表現、10進表現、(ii) 変数を含 む数式では、これに ”/ x”で区切って変数値(16進表現、10進表現)が続き、(iii) 恒等関数を設定し た場合、16進表現、コロン、基数2の仮数(10進数値)・指数表示、10進表現となる。
文字列表示器下の押しボタン SW ”PRINT”のクリック(帰還モードでは” ↷ feedback”ボタンのクリ ックにより”PRINT”のクリック無しで強制的に)でログに蓄積され文字列表示器(キャプション”log text (to be written)”)に表示される。押しボタンSW ”include comment”がONのときは、文字列制御器(ラ ベル”comment”)に入力した文字列が関数値のまとめに%に続けて追加される。(i) (ii) (iii)それぞれの例 は次の通りである。
# (2^53-1)*2^971 : 7FEFFFFFFFFFFFFF (1.7976931348623157E+308) % no variable
# gamma(x+1) : 4086800000000000 (7.2000000000000000E+2) / x : 4018000000000000 (6.0000000000000000E+0) % single variable
# 400921FB54442D18 : 1.5707963267948966*2^1 (3.1415926535897931E+0) % identity function
ログは文字列表示器(キャプション”log text (to be written)”)上の押しボタンSW ”WRITE”のクリ ックでテキストファイル(既定のファイル名は保存日時の yyyymmddHHMMSS.txt)に保存される。
ファイルダイアログをキャンセルしてもエラーとはならないのは”READ”と同様である。
LV_IEEE754.zipに同梱の20181127093506.txtは、帰還モードの実行例としてNewton-Raphson 法により6を初期値として2の正の平方根(負の初期値を設定すれば負の平方根が得られる)を求めた 例で次の操作で保存したものである。①起動時既定の状態(ログは空、変数値6)でメニューリングの 項目(x+(2/x))/2を指定し”↓RCL”をクリック、② "feedback mode” をONにし”PRINT”をクリッ ク(帰還前の関数値まとめをログに送る)、③ ”feedback converged”が点灯するまで ” ↷ feedback”を7 回クリック、④”WRITE”をクリックしプロンプトに対し既定のファイル名で保存。
2 次の収束をする Newton-Raphson 法では適切な初期値に対しては毎回精度がビット数で倍になる が、初期値が解と大きく異なる場合は相対誤差が1/2になるだけである。最も極端な場合として、上記 の②の操作で ”PRINT”をクリックする前に”READ”ボタン下のメニューリング項目で表現可能な最大 値”$7FEFFFFFFFFFFFFF % maximum positive normalized”を選択し”ENTER”をクリックして 初期値に選ぶと ” ↷ feedback”のクリックで関数値の指数値が1変化するだけ(収束するまでには1000 回以上の反復操作が必要)であることが分る。
2. LV_IEEE754DP_compact
LabVIEW概要掲載のDP_propertyは計算モードとデータ設定モードを分けていたLV_IEEE754DP の旧版(2018/02/14解説掲載)をほぼ同等の機能(テキストファイルの読込みとログの保存はクリップ ボードのコピー・ペーストで代替)でフロントパネルサイズを小さく構成したもので、これにプリセッ ト数式の項目(2、10、16 の冪乗と逆数を追加)、選択方式(文字列制御器をメニューリングに変更)、 16進文字列の扱い(左詰)を変更してLV_IEEE754DP新版と揃えたものがLV_IEEE754DP_compact である。それ以外はDP_property と同じで、プログラムの構成詳細についてはLabVIEW 概要を参照 されたい。http://www.ns.kogakuin.ac.jp/~ct13050/johogaku/LV_sample.zip
起動時既定の関数もLV_IEEE754DPと同じgamma(x+1)であるが、既定の変数値を真の1の位が仮 数LSBとなる18、点滅表示の既定をONとしている。図6は18!の真の1の位(値は0)が点灯して いるループのキャプチャ画面である。
図6 LV_IEEE754DP_compact起動時の既定の画面(18!のビット位置0が点滅)
“result ⇒ clipboard”のクリックで、”panel” ONではブールパネルの設定値、”panel OFF”では数式、
変数値(変数がある場合)、関数値がこの順でコピー対象文字列(”CLEAR”のクリックで空になる)に 連結されクリップボードにコピーされる。”clipboard ⇒ formula etc.”でペーストできる文字列は、0 個または1個の数式文字列(何れのモードでも”formula”に読込まれる)とこれに続く(空白0x20を挟 んでよい)0個または1個の数値($に続く左詰16進数字または#に続くLabVIEWの数値入力で許さ れる文字列:LabVIEWでは”1.2 e 34”の様に数値文字列中に空白を入れることはできない)である。
数値の読込先は、”panel”ONでブールパネル、”panel”OFFで変数値(数式が変数を含まないときは表 示されないが格納はされている)である。
例として起動時の既定の設定で ”CLEAR”をクリック後、クリップボードから”sqrt(x)$4” をペー ストして結果をコピーしたクリップボードの内容は次の様になる。
Formula mode : sqrt(x)
-- argument 4000000000000000 / 2.0000000000000000E+0 -- value 3FF6A09E667F3BCD / 1.4142135623730951E+0
3. LV_FP_comparison
それまで機種により異なっていた浮動小数点表現を統一した IEEE754 規格は過去の各表現の特長
(数値以外にinfiniteとindefiniteを導入:CDC6600、基数2では正規化数の仮数MSBに1を仮定し
データ中に含めない:DEC PDP-11、複数の型がある場合に長いデータでは仮数部だけではなく指数部 も拡大:Burroughs B6700)を統合したもので、LV_FP_comparisonは、DBL、SGLに拡張精度EXT
(EXTended-precision)を加えたLabVIEW の3種類の実数型データを2種類の参考表現を含めて比 較するデモVIである。
フロントパネル構成はLV_IEEE754DP_compactからブールパネルを除き(単一モード)、各表現(げ たばき指数値の表示は割愛)を横に並べたもので起動時既定の画面を図 7に示す。10の素因数に2以 外の5が含まれることから0.1の様な10進の有限小数が2進では無限循環小数となることはプログラ ミング指導の最初に指摘されるが、起動時の設定(関数:1/x、変数値:5)では 0.2=1/5 が 4 桁(24
≡1 mod 5)で循環する様子、xを変化させて0.04=1/25の循環節20桁(220≡1 mod 25)の3周期分、
1/17の10進循環節(2進8桁10進16桁)をEXTでは確認することができる。
図7 LV_FP_comparison起動時の既定の画面(1/5は2進表現では4桁で循環する)
LV_FP_comparisonの仕様は次の通りである。
● プリセット数式にEXT用の(2^64-1)*2^16320、2^-16382、2^-16445を追加。
● 1の位の点滅指定ボタン(全表現に対し1個)はEXTの指数値が0~63のとき表示し、整数表示LED の表示と点灯は各表現で個別に行う。LabVIEWの数式評価VIはDBL型専用であるため 2018/02/14 解説掲載の旧版では三角・逆三角関数、指数・対数関数、平方根、逆数のごく限られた関数のみを評価 したが、本稿掲載版ではEXT型、SGL型にそれぞれ対応させた数式評価を用いている。
● LabVIEWの数値入力では18 桁を超えた部分は無視されるため変数値入力の制御器に文字列制御器
を加え、”1.2 e 34”の様なLabVIEWでは許されない空白を含む数値文字列も入力できる。$に続く 16進文字列は左詰で20桁に満たない場合は右に0を補い20桁を超える部分は無視する。文字列制御 器が選択されているとき、関数値を変数値に帰還する ” ↷ feedback”のクリックで 16 進文字列を格納 する。
● LabVIEWの数値出力でも10進18桁を超える指定は無視されるため(変換過程の対数・指数関数の
誤差を見込んだ制限と思われ、例えばEXT型のacos(-1)の値は有効数字桁数指定で3.1415…7932、
小数点以下桁数指定で3.1415…79324までしか表示しない)、本VIのEXT型の10進表示では絶対値 が10-27以上1028未満の場合(527までの5の正整数冪は1の位がEXTの64b仮数内にあり、仮数表 現が同じ10の冪は1の位が仮数内に無い場合も浮動小数点区間代表値と一致する)に限り対数・指数 関数を用いない10進変換により3.1415926535897932385 E0の様に精度一杯の20桁を出力する(範 囲外ではLabVIEWの出力と同様に18桁に制限している)。
● ペーストできる文字列の形式はLV_IEEE754DP_compactと同様であるが、#に続く10進数値の指 数記号の前後に空白が許される。コピーされる文字列には各表現の値に加え、”negate”SWの設定、EXT 型のLSBの扱い、CDC6600の設定(後述)の情報が含まれる。
● 変則的なフォーマット(註1および図8参照)のEXT型ではデータ長16B全体ではなく数値情報を 格納する上位10Bのみを表示する。EXT型ではLabVIEW内部表現に含まれない仮数MSBのLEDを 赤ではなく黄、キャプションを”implicit MSB”ではなく”fraction MSB”と表示している。
註1:EXT型データの内容はプラットフォームにより異なるが、Intel CPUで実行する場合8087以来
の80b FPレジスタ(指数15b、仮数64b)の上位16b(符号と指数)が全て0(+0または正の非正規 化数)の場合はレジスタの80bがそのまま16Bデータの上位10Bとなる(2Dブール配列の下に”fraction MSB included”と表示:図8左)。レジスタ上位16bに0でないビットがある場合、DBL/SGLと同様 に仮数MSB(ビット位置63重み1)を外に出し上位16bと残る63bの仮数データを格納する(80b内 部表現のLSBは無効な0であり2Dブール配列右下のLSB要素に"LSB disabled"を重ねている:図8
右LabVIEWでは同じ場所に表示されたオブジェクトは時間的に後に作成されたものを優先する)。
● 参考としてメインフレーム代表機IBM System/360と元祖スーパーコンピュータCDC6600の浮動小 数点表現(EXT型の関数値をフォーマット変換したもの)を右に並べている(註2参照)。参考表現は 当該機種による実際の演算結果ではないため、関数値と変数値が等しい場合に点灯する ”converged”
LEDを置いていない。回路の基本が3b(当時は16進数字ではなく8進数字を使用した)で語長60b の CDC では、他の表現で”hexadecimal value”キャプションの文字列を押しボタン SW(キャプショ ン ”oct hex”:既定はOFF)により”octal”と”hexadecimal”を選択可能(既定は8進)としている。ま たCDCには通常のrounded命令の他にunrounded命令があり、2Dブール配列右下に切替SW(既定 はONのrounded)を置いている。
註2:IBM System/360は回路構成をそれまでの3b基本から4b基本としアドレス単位を語長から8b にした。仮数基数が2ではなく24 =16であるため、正規化数の精度は2進で見れば53b~56bと変動 する(1の位は4b単位で移動し仮数最上位に最大3bの0を生じる)。CDC6600は科学技術計算に特化 し整数演算に浮動小数点演算ユニットを流用するため仮数の1の位がMSBの左ではなくLSBの右にあ り絶対値の範囲が1より大きい側に偏っている。正規化数の逆数は0になることはあっても溢れること が無いため後継のCDC7600で逆数乗算による除算回路を採用できた。IEEE754もnビット指数のバ イアスを2n-1ではなく2n-1-1として絶対値の範囲を1より大きい側にシフトしているがシフトが軽微 なため正規化数の逆数が非正規化数になることはあっても0になることはない。符号絶対値表現ではな く1の補数表現を用いるCDC6600では押しボタンSW ”negate”の押下で全ビットが反転する。
仮数符号すなわち全体の値の符号が符号絶対値表現である IEEE754 表現では符号反転により MSB のみが反転するが、非正規化数の表現が変則的なEXT型では符号反転で仮数部分が1bずれる。図8に この様子を示す。
図8 拡張精度(EXT)の変則的なフォーマット(sqrt(.5)*2^-16382の符号反転例)
関数で10^xを指定し、”force to integer”をONにして数値スライドでxを1から増して行くとき、
DBLでは15まで、SGLでは7までそれぞれ仮数精度で整数を表現できる範囲では全て10^xの正しい 結果が得られる。DBL、SGL では必要ならば FP レジスタの下位桁を利用可能であるのに対し精度未 満のデータが無い EXT では 2 以上の整数で全て真の値とは異なっている(x=19 で”non-zero integer”LEDが点灯するが、これは仮数LSBが1の奇数1019+1であり単に整数の表現であるというだ けで正しい値ではない)が10^xをexp(x*ln(10))で求めた場合に通常発生する誤差であって演算精 度自体が悪いのではない。
4. LV_num_int
数値積分の台形則とSimpson則の比較では、区間[0,1]における1/(1+x2)の積分を取上げることが多い が、この被積分関数と区間の組合せは Simpson 則の精度が良くなる特殊な例であって公式の違いによ る精度の比較としては適切でない。これは区間[-1,1]の積分(偶関数では [0,1]での積分も同様)におけ るSimpson則の誤差特性関数の零点が区間分割数を増したときに±i に収束し被積分関数1/(1+x2)の極 を相殺するためである。一方、有限な端の無い積分では関数値の重みに分点の偶数奇数による差を付け ない台形則の精度が圧倒的に良くなる(「Simpson則が悪くなる」のではない)。LV_num_intは被積分 関数と積分区間を数式文字列で設定し、条件による両積分公式の得失を確認するデモVIである。
LV_num_intのメニューリング(キャプション ”integrand”)には、ユーザ定義の ”custom”の他、図 9 に示す 3 種の被積分関数が固定(一部のパラメータは自由)されている。図9 左は上述の Simpson 則がよい結果を与える例で、区間を[0.5,1.5]に変更して区間[0,1]が特殊であることを確かめることがで きる。
図9 台形則(上段)とSimpson則(下段)による数値積分結果の相対誤差の比較
図9中は有限な端の無い積分の代表例である exp(-x2)の区間 [0,∞) における積分である(偶関数で あり区間 (-∞,∞)の積分と同様)。exp(-x2) は急速に減少し DBL 型の演算精度では積分区間の上限∞
を6とすれば十分である(erfc(6)≒2.15×10-17)。周期関数の1周期積分も有限な端の無い例で、図9 右は整数次数のBessel関数を求める例である。
メニューリングで “custom”を選択すると被積分関数を数式文字列で設定する文字列制御器(ラベ ル ”f(x)”:起動時既定はx^3)と押しボタンSW ”reference”(起動時設定ON)が表示され、”reference”ON では相対誤差表示の基準値を文字列制御器(キャプション ”primitive function F(x) or const. r”:起動
時既定はx^4/4)に原始関数または定積分の値を数式文字列で設定する。
p.10図10左はSimpson則で誤差の無い3次式(起動時既定)、図10中は誤差が標準的な収束をす る被積分関数、図10右は区間の端で微分できないためSimpson則の優位性が失われる(誤差の収束は 両公式共に区間分割数の-3/2乗に比例)被積分関数の例である。単に数値積分の結果を知る目的であれ ば押しボタンSW ”reference”をOFFにする(相対誤差の数値表示器は非表示となる)。
図10 “custom”を選択して相対誤差が特徴的な被積分関数を設定した例
5. LV_DE
±∞で急速に減少する関数の無限区間での積分が台形則で精度よく求められることを利用して被積 分関数に変数変換を行い台形則を適用する二重指数関数型数値積分公式(森・高橋によるDE公式)は 科学技術計算ソフトウェアMathematicaのNIntegrateにも実装されている。DE公式による積分では 区間両端付近の関数値の寄与が変数変換により圧縮されるため図 10右の様な端点で微分ができない場 合であっても十分な精度が得られ、更に被積分関数が有限区間の端点で±∞に発散するが端点との差を εとして発散が1/εよりも遅く広義積分が存在する場合もDE公式では求めることができる。ただし、「公 式では」であって、データ長が有限な演算では端点付近の寄与を正しく反映させる配慮が必要であり、
DE公式による数値積分デモVIの LV_DE(2018/02/14解説掲載の旧LV_DErev)で確認できる。DE 公式の具体的な変換の式は被積分関数の性質、原積分の積分区間により異なり下記の文献に詳しい。
http://www.kurims.kyoto-u.ac.jp/~okamoto/paper/Publ_RIMS_DE/41-4-38.pdf 起動時のLV_DEのフロントパネル画面を図11に示す。
図11 LV_DEの起動時既定画面 広義積分の例
LV_DEの仕様は次の通りである。
● 数式文字列で被積分関数と積分区間の端点(±Infはそれぞれ1/0、-1/0)を入力する。
● 入力した積分区間に応じて変数変換の式(上記文献の式4.8、4.16~4.18)が選択される。
● 区間端点を除く数式文字列が数式として評価できない場合 ”error” LED を当該文字列制御器の上に 点灯表示する。
● 積分区間の種類と端点の順序を表示器クラスタ(ラベル ”interval”)に表示する。
● 変数変換後の積分区間は(-∞,∞)であるが実際には前項の台形則によるexp(-x2)の積分に見る様に 大きな値を必要としない。起動時既定は分点の絶対値上限を4、区間の刻み幅を0.15としている。
● 相対誤差を表示する場合は押しボタンSW ”reference”(起動時既定ON)をONにし原始関数(端点 で不定形の極限となる場合は極限値が 0 となる積分定数を選び端点に対応する 0 強制の押しボタン SW ”F(a)=0”、”F(b)=0”をONにする)または定積分を数式文字列で入力する。
● 積分区間が有限の場合は式4.8により変換する。被積分関数が端点で発散する広義積分では変換過程 であふれを生じ変換後の積分区間を十分に確保できない(既定の絶対値上限では結果が±Inf となる)。 有限区間の場合に表示される押しボタンSW ”alternative”をON(端点毎に独立)にして端点の近傍(起 動時の既定は10-5)では端点との差を引数とする関数形を指定(関数を入力する文字列制御器のキャプ ションは端点の大小関係により変化)して変数変換過程を別処理(変換の式自体は同じ)で行いあふれ を防止する(p.10図11、図12左)。
● 積分区間が半無限の場合は式 4.16 により変換する。被積分関数の絶対値が指数関数的に減少する場 合は、この変換では|x|→∞で|f(x)|を過剰に減少させ十分な精度が得られないため、半無限である 場合に表示される押しボタン SW ”exponentially decreasing”(起動時既定は OFF)をON にして式 4.17の変換を選択する(図12中)。
● 積分区間が無限である場合は式4.18により変換する(図12右)。
起動時既定の被積分関数1/√(1-x2) は端点±1の近傍で 1-x2の桁落ちのため Infとなる。変数変換 後の積分区間を桁落ちが生じない範囲に限定したのでは端点付近の寄与を反映できないため、そのまま 公式に従って積分すると高々10進10桁程度の精度しか得られない。起動時の既定は1/√(1-x2) を±1 の近傍でεの関数として与え(註参照)区間 (-1,1) の広義積分を演算精度相応で得ている(p.10図11)。
註:-1の近傍ではx+1をε、+1の近傍では1-xをεとおいて1/√(1-x2) をεで級数展開した2次の項 まで取った数式((3*e/8+1)*e/4+1)/sqrt(2*e)(eはεの代り:両端点で同じ関数形になる)を指 定している。近傍を10-5とした場合に1/√(1-x2) の個別の数値としてはDBL型の精度に不足するが、
総和全体への寄与を考えた精度では十分なため3次以上の項を含める必要は無い。
図12に各種の積分区間での積分の例を示す。何れも既定より少し広い区間刻み0.2(分点数41)で それぞれ演算精度相応の結果を得られることが分る。図12左はln(x)の区間 (0, 1]での広義積分の例 で ”alternative” ONのeの関数形は被積分関数と同じ(端点が0)になる。端点での発散が遅いため 同じ区間刻みでOFF(既定の上限では積分値が-Infとなる)にして絶対値上限をあふれない3に下げ ても相対誤差は-5.88418203051333E-15とON の2.22044604925031E-16から極端に増えるこ とはない。
図12中はexp(-x)の[0,∞)での積分で ”exponentially decreasing”をONにしている。OFFにする と相対誤差は1.67237917276797E-5とONの2.22044604925031E-16から極端に悪化する。図 12右は1/(1+x2)の積分区間を無限区間に広げ端点の大小関係を逆にしてinterval orderを”reverse”と表 示させた例(積分の値は-π)である。積分区間の端点の一方を 0にすると(偶関数の被積分関数に対 し本質的に同じ積分)半無限区間の変換が選ばれ相対誤差は2.87211708194235E-7と無限区間での -3.33066907387547E-16からこれも極端に悪化し、被積分関数の性質に合った積分区間を選ぶ必要 があることが分る。
図12 各種の積分区間での積分例:有限区間(左)、半無限区間(中)、無限区間(右)
備考
ランタイムエンジンについて
LabVIEW の実行ファイル(スタンドアロンアプリケーション)を開くには当該バージョンの
LabVIEW本体またはランタイムエンジンが必要である。センターのPCにはLabVIEW2017がインス
トールされており、LV2017フォルダ中の実行ファイルをそのまま開くことができるがLV2013フォル ダ中の実行ファイルは開けない。
個人所有PCでLabVIEW2017版実行ファイルを開くにはバーチャルカフェテリア(本学の学生は特
別な利用申請無しに利用可能)にログインするかまたは以下のページでNIユーザアカウントを登録し てダウンロードしたLabVIEW2017ランタイムエンジンをインストールする。
ダウンロードページ:http://www.ni.com/download/labview-run-time-engine-2017-sp1/7191/en/
ファイル名:LVRTE2017SP1_f3Patchstd.exe、ファイルサイズ:363.29 MiB
NIユーザアカウントの登録に抵抗のある人は次の場所からLabVIEW2013ランタイムエンジンをダ ウンロード(登録不要)、インストールしてLV2013フォルダ中の実行ファイルを開く。
http://ftp.ni.com/support/softlib/labview/labview_runtime/2013/Windows/LVRTE2013std.exe ファイル名:LVRTE2013std.exe、ファイルサイズ:257 MiB
実行時ライブラリとフォント指定について
数式評価VIを使用する実行ファイルを開くには版に対応した実行時ライブラリlvanlys.dll(LV2017、
LV2013それぞれのdataサブフォルダ中にあり両者に互換性は無い)が必要で、それぞれのサブフォル
ダdataを実行ファイルと同じ場所に置いて開く。
LabVIEW では、フォントを VI の作成時に明示的に指定することもできるが、通常実行ファイルを
作成する場合には(指定フォントがインストールされていない環境では表示できなくなるため)OS 既 定のフォントを使用しており、作成時のフォントと実行時のフォントが異なる場合、文字列が枠に収ま らない、配列が斜めにずれるなどのレイアウトの乱れを生じる。以下の4行が書かれた同梱の設定ファ イル(VIと同じファイル名で拡張子が.ini)を実行ファイルと同じ場所に置いて開くこと。
LV2017版設定ファイル
[VIファイル名(拡張子無し)] appFont="メイリオ" 17 dialogFont="メイリオ" 17 systemFont="メイリオ" 17 LV2013版設定ファイル
[VIファイル名(拡張子無し)] appFont="MS UI Gothic" 12 dialogFont="MS UI Gothic" 12 systemFont="MS UI Gothic" 12
数式文字列について
数式評価VIによる解釈では以下の点に注意が必要である。
● 数値制御器のボックスに入力する数値文字列では指数記号として大文字のE 以外に小文字の e を使 えるが空白は許されない(”1.2 e 34”はエラーとなる)。一方、文字列制御器に入力する数式文字列 では数値文字列の空白は許されるがeが変数とみなされるため指数記号として使えるのは英大文字Eの みである。また数値制御器で入力できるInf、-Inf、NaN(小文字で入力した場合は大文字に変換され る)を数式文字列では評価できないためそれぞれ1/0、-1/0、0/0(NaNと0/0は共にquiet NaNで あるが両者の内部表現は異なる)で代替する。
● 型が浮動小数点数の数値制御器のボックスに入力した文字列 ”-0”は負号付0(MSBのみ1で他の
ビットが 0)と解釈されるが、数式文字列では通常の 0(全ビットが 0)と解釈される。浮動小数点数
に対して数値関数の反転(negate: )もフォーミュラノードの負号も常にMSBを反転するが、引数 0を入力した数式 ”-x”の結果のMSBは0のままで0の符号反転ができないことに注意する。