Journal of Chemical Software, Vol.1, No.2, p.99

ＵＢＡＳＩＣの化学計算への応用

成沢芳男、宮前雄一

Ⅰ．緒言

　ＵＢＡＳＩＣの化学計算への応用について述べる。ＵＢＡＳＩＣは木田祐司１）によって開発された、巨大整数計算のための多倍長計算用ＢＡＳＩＣで、２６００桁の整数計算および実数計算、更には１３００桁の複素数計算と有理数計算ができる。その他に整数論のための関数や複素数をサポートした初等関数も扱えるので、数学の分野では種々応用されているが、残念ながら化学分野への応用はこれまでない。
　化学計算においては現在最も広く利用されている言語はＦＯＲＴＲＡＮであり、それに代わる言語はまだ開発されていない。われわれはＵＢＡＳＩＣをＦＯＲＴＲＡＮに代わるものとして提案しようという考えではなく、個別の計算で非常に高精度を必要とする場合に、ＵＢＡＳＩＣが大変に威力を発揮することをいくつかの計算例によって示そうと思う。
　ＵＢＡＳＩＣはパソコンで２５００桁のπの計算をわずか６秒で実行する。これはこれまでの言語では達成できない能力で、今後その高精度演算機能は化学の色々な分野で応用されるであろう。πの計算に限って考えれば、いわゆるスーパーコンピューターによる性能テストで日米欧の技術競争があって、数万桁の計算を正確に行うことが報告されている。しかしスーパーコンピューターは誰でも気軽に使えるものではなく、またπの計算をすることが目的ではない。化学計算でよく遭遇する問題は最小２乗近似であろう。例えば多項式近似により最小２乗を行なおうとすると多項式のベキに相当する次数のマトリックス計算をする必要が生じる。近似の次数を上げれば一般に精度は良くなると考えられるが、膨大な計算のため計算誤差が大きくなり、逆に精度は下がる。しかしこれまでそのような計算誤差を検証する高精度演算ソフトがなく、化学の分野で系統的な検討が行われてこなかった。本論文では最小２乗近似の計算誤差を検討するためにＢＡＳＩＣ、ＦＯＲＴＲＡＮ及びＵＢＡＳＩＣのプログラムを作成し、演算誤差について論じる。

Ⅱ．序論

　著者はすでに水の密度の多項式近似と測容器の補正２）の論文で、ＢＡＳＩＣによる多項式近似の計算を行った。水の密度の温度依存は４℃で0.999973の最大値を取り、０℃から３０．５℃の範囲で0.999841～0.995494の変化しかなく、比較的次数の小さい多項式で近似でき、なおかつ有効桁数も小数点以下７桁程度であるため当時の８ビットパソコンで速度、精度とも充分であった。しかし更に複雑な曲線を最小２乗近似するためには、速度、精度とも高度なものが要求され、一挙に３２ビットパソコンによる高精度計算が必要となるが、高次元マトリックス計算においてどれだけの演算誤差が生じるのかと言う点に関しては全く研究例がない。そこでＵＢＡＳＩＣを用い、他言語による計算結果と比較し、その有用性について論じたい。

Ⅲ．理論

１．最小２乗法３），４）
　最小２乗近似は区間（ａ，ｂ）で与えられた関数ｇ（ｘ）に対して、適当な条件を満足する関数の中から、２乗誤差

を最小にするｆ（ｘ）を選定する方法である。
　区間（ａ，ｂ）で定義された関数ｇ（ｘ）がその区間内に指定された選点

と、その上の値

とのみによって与えられている場合、適当な条件を満足するある範囲の関数からｆ（ｘ）を選んで、２乗誤差

を最小にすればよい。
　ｆ（ｘ）がある既知関数系｛ｆｋ（ｘ）｝の１次結合として

で与えられる場合には、２乗誤差Ｅを最小にする条件方程式は

である。

２．Ｇａｕｓｓの消去法５），６），７）
　一般にｎ元連立１次方程式

を解くには、(2-1)の両辺をａ１１で割る。

ただし　ａ１２'＝ａ１２／ａ１１ ,ａ１ｎ'＝ａ１ｎ／ａ１１ ,ｂ１'＝ｂ１／ａ１１である。これにａ２１を掛けて(2-2)から引く。

ただし　
　　　　である。以下同様にして

までを作る。一般にａｉｊ' およびｂｉ' の計算式は

である。
これでｘ１が消去された。　　　　　　　　　　消去の第１段
次に(2-2')をａ２２'で割り、そのａｉ２'倍を式(2-i')から引く（ｉ＝３，・・・・，ｎ）。
その際の係数の計算式は

である（ｊ＝２，・・・・，ｎ　，またｉ＝３，・・・・，ｎ）。これでｘ２が消去される。
　　　　　　　　　　　　　　　　　　　　　　消去の第２段
あとは同様にして第ｋ段の計算は次のようになる。

ただし右肩の（ｋ－１）は「第ｋ－１段の結果」を表す。このような処理をｋ＝ｎ－１まで実行すれば、ｘ１，ｘ２，・・・・，ｘｎ－１は消去されて、ｘｎだけを未知数とする方程式ができる。

そして得られた結果を順に前の式に代入していけば、残りの未知数が全部求まる。

３．曲線のあてはめ８）
　離散的な点で与えられたデータに滑らかな曲線をあてはめることを考える。データは

とし、これから近似式ｙ＝ｆ（ｘ）を作ることが目的である。一般にｎ個のパラメータを含む関数ｆ（ｘ，ｐ１，ｐ２，・・・・，ｐｎ）を作り

になるようにｐ１，ｐ２，・・・・，ｐｎを定めるのが１つの解答である。ここでは議論を簡単にするために、ｆがｐ１，ｐ２，・・・・，ｐｎに関して線形１次の形

で表されるとする。ＵＢＡＳＩＣによる数値計算の演算精度を論ずるので複雑な非線形モデルは考慮しない。このようにすれば(3-2)は

という連立１次方程式になるから、これを解けばｐ１,・・・・,ｐｎが求められる。
　実験データ(3-1)は計測誤差を含むので、パラメーターの数をｎより少なくして、(3-2)の代わりに

が最小という条件でパラメーターを決める方が、自然な結果が得られる。式(3-5)を満たすｐ１，・・・・，ｐｎを求めるには

を連立方程式として解けば良い。ｆ（ｘ）が式(3-3)の形の場合には、これは

ただし

となる。式(3-7)は正規方程式と呼ばれる。

Ⅳ．計算

　計算に用いたＢＡＳＩＣ及びＦＯＲＴＲＡＮのプログラムは水の密度を求める多項式近似の計算に用いたもの２）を一部修正して用いた。ＢＡＳＩＣは高精度演算を、ＦＯＲＴＲＡＮはＭＳ－ＦＯＲＴＲＡＮで倍精度演算、Ｓｕｎ－ＦＯＲＴＲＡＮで４倍精度演算を行った。ＵＢＡＳＩＣはＢＡＳＩＣプログラムを改良して用いた。
　０℃から３０．５℃まで 0.5℃刻みの水の密度データ９）を用いＢＡＳＩＣ、ＦＯＲＴＲＡＮ及びＵＢＡＳＩＣプログラムによる各次元ごとの計算結果とデータとの差の２乗和を求めた。用いたパソコンは演算コプロセッサーを内蔵したＮＥＣ　ＰＣ９８０１ＲＡ及びＰＣ９８０１ＦＡであり、またＦＯＲＴＲＡＮ計算はＭＳ－ＦＯＲＴＲＡＮ　Ｖｅｒ．５．１及びＳｕｎＦＯＲＴＲＡＮ１．４を用いてコンパイル及びリンクをして行なった。

Ⅴ．結果と考察

　ＢＡＳＩＣ、ＦＯＲＴＲＡＮ及びＵＢＡＳＩＣプログラムによる各次元ごとの計算結果と水の密度データとの差の２乗和を計算した。０℃から３０．５℃までの水の密度は変化が小さく、５次の多項式で近似すれば十分であって、水の密度データの最適化では各プログラムの差は顕著に示されない。そこで０℃から３０．５℃まで 0.5℃刻みの６２ヶのデータを複雑な変化のあるものに変え、より高次の多項式でなければ近似できないようにして各プログラム相互の差異を検討した。すなわち温度をｘ、密度をｙとし、新しい縦軸の値をＹとする。プログラムの内容を変えずにデータだけを変更し、また計算結果をグラフィック表示する関係から次式

　　　　Ｙ＝ 0.80ｙ＋ sin ｘ × 0.0003 ＋ 0.199 (Ａ)

を計算して、水の密度のデータと同様にｘの値が０から３０．５まで 0.5刻みで６２ヶのデータ点を持つデータファイルを作成し、ｘとＹの関係曲線に最もよくあてはまる多項式近似を行った。各プログラムによる水の密度データおよび(Ａ)式のデータについての誤差の２乗和を計算し、結果を Table１および２に示す。ただしＢＡＳＩＣの計算はＳｕｎ－ＦＯＲＴＲＡＮ及びＵＢＡＳＩＣのものと１２次元まで全く同じであったが、１３次元以降はOVERFLOWのメッセージが表示され計算不能になったので、表には示さなかった。表中の DIMENSION（次元）はマトリックスの次元を示し、多項式の最高ベキ数はこれより１だけ減じたものである。表から明かな通りＵＢＡＳＩＣは多倍長計算ソフトの特長を示し、いずれのデータでも５０次元まで計算誤差が漸減した。５１次元以上の計算は実行していないが、次元数の限界はコンピューターのメモリ環境によって異なり、ＵＢＡＳＩＣの限界ではない。
　ＵＢＡＳＩＣについて(Ａ)式のデータを１０次元及び２０次元で近似したときの最適化の状況をFig.１に示す。１０次元では曲線のあてはめはかなり不満足であるが、２０次元では最適化はかなり満足すべきものであった。
Table １のＦＯＲＴＲＡＮ及びＵＢＡＳＩＣの誤差の２乗和を次元に対してプロットし、水の密度に対する結果をFig.2(a)に、（Ａ）式に対する結果をFig.2(b)に示す。Fig.2(a) の結果が示す通り、比較的低次で多項式近似できる場合は各プログラムによる差は顕著に現れない。ただしＭＳ－ＦＯＲＴＲＡＮの場合は３０次元と３９次元のところに特異的に精度が悪くなるところがあるが、この点に関しては後に議論する。Fig.2(b) の結果はＭＳ－ＦＯＲＴＲＡＮもＵＢＡＳＩＣやＳｕｎ－ＦＯＲＴＲＡＮと同じく１１次元までは全く同じであるが、１２次元以降は誤差が減少しないし、バラツキがひどい。これはＭＳ－ＦＯＲＴＲＡＮがDOUBLE PRECISION型の計算をするために、精度的に適さないことを示す。なお(Ａ)式のデータについてＵＢＡＳＩＣ及びＳｕｎ－ＦＯＲＴＲＡＮの誤差のカーブは単調減少ではなく、階段状になったが、データが複雑になったため最適化するのにより高次の多項式近似が必要であることを示している。
　なおFig.2(a)及び(b)のＭＳ－ＦＯＲＴＲＡＮの誤差の２乗和が３０次元と３９次元で特異的に大きくなることに関して原因は明かでないが、理由はいくつか考えられる。即ち１：パソコンを用いていること、２：演算コプロセッサーを使っていること、３：ＭＳ－ＦＯＲＴＲＡＮのコンパイラー及びリンカーを使っていること及び４：プログラムが正しく書かれていないこと。そこでEthernetを通してＰＣ９８０１ＦＡを富士通Ｓ４／４７０ワークステーションにアクセスし、Allied Telesisの通信ソフトTCP/IPを用いてプログラムおよびデータを転送し、ワークステーション上で各次元毎の水の密度データ及び（Ａ）式に対する誤差の２乗和を計算した。なお実数計算はすべてQUADRUPLE PRECISION 型で行った。Fig.2から明かなように、ＭＳ－ＦＯＲＴＲＡＮはDOUBLE PRECISION型で、精度上問題があったが、ワークステーションによるQUADRUPLE PRECISION の計算ではこの点が改善された。従って、ＭＳ－ＦＯＲＴＲＡＮの計算結果に異常が生じるのは上記４項のプログラムが正しく書かれていないことであるという理由は除かれる。また上記１～３項は１つに纏めることができ、パソコンによる倍精度演算では精度上問題があることを示しており、それ以上原因を明らかにするにはＭＳ－ＦＯＲＴＲＡＮの計算過程を詳細に検討する必要がある。それは今回のＵＢＡＳＩＣに関する性能検査の問題ではないので、これ以上の議論はしない。
　再びFig.１の話に戻って、Fig.１のグラフはＵＢＡＳＩＣで計算した結果ではあるが、 Table ２より１次元から２３次元の範囲においてはＵＢＡＳＩＣとＳｕｎ－ＦＯＲＴＲＡＮの計算結果は同じであるので、この１０次元と２０次元についての結果を表したグラフは、もしＳｕｎ－ＦＯＲＴＲＡＮで計算した結果を同様にグラフ表示しても同じものが得られると言える。すなわち、(Ａ)式のデータではＵＢＡＳＩＣとＳｕｎ－ＦＯＲＴＲＡＮの間の顕著な差はグラフ上では見られない。しかしながら、（Ａ）式のデータでは２０次元でほぼ満足な近似が出来ているが、Ｓｕｎ－ＦＯＲＴＲＡＮとＵＢＡＳＩＣが同じ誤差の２乗和を与えているのは２３次元までであり、より複雑なデータを用いれば精度の差が明らかになると考え、水の密度のデータを（Ａ）式より更に複雑にした次式

　　　　Ｙ＝ 0.80ｙ＋ sin ２ｘ × 0.0003 ＋ 0.199 (Ｂ)

を計算して、（Ａ）式のデータと同様に６２ヶのデータ点を持つデータファイルを作成し、ＭＳ－ＦＯＲＴＲＡＮ、Ｓｕｎ－ＦＯＲＴＲＡＮ、およびＵＢＡＳＩＣの３種類の言語で多項式近似および誤差の２乗和の計算を行い、結果をTable ３に、またこれらの次元に対するプロットをFig.3に示す。Fig.3より明らかなように、（Ｂ）式のデータを用いた場合には、ＵＢＡＳＩＣで計算した誤差の２乗和は２７次元から急激に減少しているが、Ｓｕｎ－ＦＯＲＴＲＡＮの結果はあたかもFig.2(b)のＭＳ－ＦＯＲＴＲＡＮのような、誤差の２乗和が減少しないバラツキのひどい結果となった。
　ＵＢＡＳＩＣについて（Ｂ）式のデータを２０次元及び３５次元で近似したときの最適化の状況をFig.4に示す。（Ｂ）式のデータでは、２０次元では満足な曲線のあてはめができていないが、３５次元ではかなり満足な結果が得られた。Ｓｕｎ－ＦＯＲＴＲＡＮでＵＢＡＳＩＣと同じ結果が出せる限界が２３次元であり、Fig.4の２０次元の近似曲線はほぼＳｕｎ－ＦＯＲＴＲＡＮプログラムの近似できる限界に相当するので、このグラフからもＵＢＡＳＩＣの精度的な優秀性が確認できた。
　従来、誤差の２乗和を計算して評価をするための高精度計算プログラムがなかったので、結果を比較することができなかったが、今回ＵＢＡＳＩＣを用いて、広く使われている基本的な言語の計算精度を論じることができた。

Table 1 水の密度データ使用時の誤差の２乗和

Table 2 （Ａ）式によるデータ使用時の誤差の２乗和

Table 3 （Ｂ）式によるデータ使用時の誤差の２乗和

Fig.1. Simulation curves with UBASIC program
:data points calculated from equation(A) (see text);
:dotted lines are simulation curves calculated with
10th and 20th dimension.

Fig.4. Simulation curves with UBASIC program
:data points calculated from equation(B) (see text);
:dotted lines are simulation curves calculated with 20th and 35th dimension.

REFERENCES

1) 木田祐司, "ＵＢＡＳＩＣ８６ Ver8.1 ユーザーズマニュアル", 日本評論社 (1990).
2) 成澤芳男・中埜邦夫, 分析化学, 32, T25 (1982).
3) 数学ハンドブック編集委員会編, "理工学のための数学ハンドブック", p.326, 丸善 (1960).
4) 一松信・竹之内脩編, "改訂増補新数学事典", 増補p.41, 大阪書籍 (1979).
5) 戸川隼人著, "計算機のための数値計算", p.18, サイエンス社 (1976).
6) 戸川隼人著, "ＢＡＳＩＣによる線形代数", p.46, 共立出版 (1985).
7) 矢野健太郎監修，宮本敏雄訳編, "数学ハンドブック", p.834, 森北出版 (1985).
（I.N.ブロシュテイン・K.A.セメンジャーエフ原著、G.グロッシェ・V.ツィーグラー・D.ツィーグラー編集）
8) Ref. 5, p.39.
9) K. Sch fer, G. Beggerow, "Landolt-B rnstein", 6Aufl, ⅡBand, 1Tl, s.36(1971),
　(Springer Verlag, Berlin�%�Heidelberg�%�New York).

Return