(57)【要約】【課題】逆数または平方根逆数の近似計算をニュートンラフソン法によって行う際、同程度の大きさの浮動小数点数の減算処理に伴ってプロセッサで生じるビットシフト処理に起因するストールを回避でき、計算速度を向上できる演算装置を提供する。【解決手段】演算装置１は、初期近似精度を持つ初期近似逆数に対して近似計算を行う。まず、中間値演算部１５が、ニュートンラフソン法の１段近似計算により、中間近似精度を持つ中間近似逆数を計算し、そして、目標値演算部１７が、ニュートンラフソン法の多段近似計算によって目標精度まで精度を高める。中間近似精度は、プロセッサでの浮動小数点数の減算処理に伴う仮数部のビットシフトによるストールが生じない上限ビットシフト量に対応する非ストール誤差を残すように設定されている。この中間近似精度が得られるように、中間値演算部１５では、故意に誤差を拡大した１段近似計算が行われる。【選択図】図２
2008287555.tif 000002
【特許請求の範囲】
【請求項１】
逆数を求めるべき入力値に対応し所定の初期近似精度を有する逆数である初期近似逆数に対して、ニュートンラフソン法による近似計算を行って、予め設定された目標精度を有する目標逆数を算出する演算装置であって、
浮動小数点数を処理する演算機能を有するプロセッサと、逆数の近似計算プログラムを記憶する記憶部とを有し、前記記憶部に記憶された前記近似計算プログラムを前記プロセッサにて実行することによって中間値演算部および目標値演算部が前記演算装置に備えられ、
前記中間値演算部は、前記初期近似逆数に対して、ニュートンラフソン法による１段近似計算を少なくとも１回行って、前記初期近似精度と前記目標精度の間に設定された所定の中間近似精度を有する中間近似逆数を算出する中間値演算処理を行い、
前記目標値演算部は、前記中間近似逆数に対して、前記中間近似精度から前記目標精度まで精度が上がるようにニュートンラフソン法による多段近似計算を行って、前記目標逆数を算出する目標値演算処理を行い、
前記中間近似精度は、前記プロセッサにおける浮動小数点数の減算処理に伴う仮数部のビットシフトによるストールが生じない上限ビットシフト量に対応する非ストール誤差が残る精度に設定され、前記中間値演算部は、前記非ストール誤差が近似計算結果に残るように誤差を拡大した誤差拡大１段近似計算を行うことにより、前記中間近似精度を越えないように近似精度を制限した前記中間近似逆数を算出することを特徴とする演算装置。
【請求項２】
前記中間値演算部は、前記入力値をａとし、近似計算前後の近似逆数をｘｉ、ｘｉ＋１としたときに、前記誤差拡大１段近似計算として、ｘｉ＋１＝ｘｉ×（（２＋α）−a×ｘｉ）を求める処理を行い、誤差成分αが前記非ストール誤差に応じて設定されていることを特徴とする請求項１に記載の演算装置。
【請求項３】
前記プロセッサがＩＥＥＥ７５４により規定される浮動小数点数を処理し、前記初期近似精度が８ｂｉｔ精度であり、前記目標近似精度が５３ｂｉｔ精度であり、前記中間近似精度が１４ｂｉｔ精度であることを特徴とする請求項２に記載の演算装置。
【請求項４】
前記誤差成分αが、２−１４であることを特徴とする請求項３に記載の演算装置。
【請求項５】
前記目標値演算部は、前記中間近似逆数をｘｉとし、εｉ＝１−ａ×ｘｉとしたときに、前記多段近似計算として、ｘｉ×（１＋εｉ×（１＋εｉ×（１＋εｉ×（１＋εｉ））））を求める処理を行うことを特徴とする請求項１〜４のいずれかに記載の演算装置。
【請求項６】
浮動小数点数を処理する演算機能を有するプロセッサと、平方根逆数の近似計算プログラムを記憶する記憶部とを有し、前記記憶部に記憶された前記近似計算プログラムを前記プロセッサにて実行することにより、平方根逆数を求めるべき入力値に対応し所定の初期近似精度を有する平方根逆数である初期近似平方根逆数に対して、ニュートンラフソン法による近似計算を行って、予め設定された目標精度を有する目標平方根逆数を算出する演算装置であって、
前記初期近似平方根逆数に対して、ニュートンラフソン法による１段近似計算を少なくとも１回行って、前記初期近似精度と前記目標精度の間に設定された所定の中間近似精度を有する中間近似平方根逆数を算出する中間値演算部と、
前記中間近似平方根逆数に対して、前記中間近似精度から前記目標精度まで精度が上がるようにニュートンラフソン法による多段近似計算を行って、前記目標平方根逆数を算出する目標値演算部と、
を備え、
前記中間近似精度は、前記プロセッサにおける浮動小数点数の減算処理に伴う仮数部のビットシフトによるストールが生じない上限ビットシフト量に対応する非ストール誤差が残る精度に設定され、前記中間値演算部は、前記非ストール誤差が近似計算結果に残るように誤差を拡大した誤差拡大１段近似計算を行うことにより、前記中間近似精度を越えないように近似精度を制限した前記中間近似平方根逆数を算出することを特徴とする演算装置。
【請求項７】
前記中間値演算部は、誤差拡大を行わない１段近似計算を１回と前記誤差拡大１段近似計算を１回行うように構成され、前記入力値をａとし、近似計算前後の近似平方根逆数をｘｉ、ｘｉ＋１としたときに、誤差拡大を行わない１段近似計算として、ｘｉ＋１＝ｘｉ＋ｘｉ×（０．５−０．５×a×ｘｉ）を求める処理を行い、前記誤差拡大１段近似計算として、ｘｉ＋１＝ｘｉ＋ｘｉ×（（０．５＋β）−０．５×a×ｘｉ）を求める処理を行い、誤差成分βが前記非ストール誤差に応じて設定されていることを特徴とする請求項６に記載の演算装置。
【請求項８】
前記プロセッサがＩＥＥＥ７５４により規定される浮動小数点数を処理し、前記初期近似精度が５ｂｉｔ精度であり、前記目標近似精度が５３ｂｉｔ精度であり、前記中間近似精度が１４ｂｉｔ精度であることを特徴とする請求項７に記載の演算装置。
【請求項９】
前記誤差成分βが、２−１５であることを特徴とする請求項８に記載の演算装置。
【請求項１０】
前記目標値演算部は、前記中間近似逆数をｘｉとし、εｉ＝１−ａ×ｘｉ２としたときに、多段近似計算として、ｘｉ×（１＋（1/2）×εｉ＋（3/8）×εｉ２＋（5/16）×εｉ３）を求める処理を行うことを特徴とする請求項６〜９のいずれかに記載の演算装置。
【請求項１１】
浮動小数点数を処理する演算機能を有するプロセッサを備えた演算装置に記憶され、前記プロセッサにより実行される逆数の近似計算プログラムであって、逆数を求めるべき入力値に対応し所定の初期近似精度を有する逆数である初期近似逆数に対して、ニュートンラフソン法による近似計算を行って、予め設定された目標精度を有する目標逆数を算出する処理を前記プロセッサに実行させる逆数の近似計算プログラムにおいて、
前記初期近似逆数に対して、ニュートンラフソン法による１段近似計算を少なくとも１回行って、前記初期近似精度と前記目標精度の間に設定された所定の中間近似精度を有する中間近似逆数を算出する中間値演算処理と、
前記中間近似逆数に対して、前記中間近似精度から前記目標精度まで精度が上がるようにニュートンラフソン法による多段近似計算を行って、前記目標逆数を算出する目標値演算処理とを前記プロセッサに実行させ、
前記中間近似精度は、前記プロセッサにおける浮動小数点数の減算処理に伴う仮数部のビットシフトによるストールが生じない上限ビットシフト量に対応する非ストール誤差が残る精度に設定され、前記中間値演算処理は、前記非ストール誤差が近似計算結果に残るように誤差を拡大した誤差拡大１段近似計算を行うことにより、前記中間近似精度を越えないように近似精度を制限した前記中間近似逆数を算出する処理であることを特徴とする逆数の近似計算プログラム。
【請求項１２】
浮動小数点数を処理する演算機能を有するプロセッサを備えた演算装置にて実行される逆数の近似計算方法であって、逆数を求めるべき入力値に対応し所定の初期近似精度を有する逆数である初期近似逆数に対して、ニュートンラフソン法による近似計算を行って、予め設定された目標精度を有する目標逆数を算出する処理を前記演算装置に実行させる逆数の近似計算方法において、
前記初期近似逆数に対して、ニュートンラフソン法による１段近似計算を少なくとも１回行って、前記初期近似精度と前記目標精度の間に設定された所定の中間近似精度を有する中間近似逆数を算出する中間値演算処理と、
前記中間近似逆数に対して、前記中間近似精度から前記目標精度まで精度が上がるようにニュートンラフソン法による多段近似計算を行って、前記目標逆数を算出する目標値演算処理とを前記プロセッサに実行させ、
前記中間近似精度は、前記プロセッサにおける浮動小数点数の減算処理に伴う仮数部のビットシフトによるストールが生じない上限ビットシフト量に対応する非ストール誤差が残る精度に設定され、前記中間値演算処理は、前記非ストール誤差が近似計算結果に残るように誤差を拡大した誤差拡大１段近似計算を行うことにより、前記中間近似精度を越えないように近似精度を制限した前記中間近似逆数を算出する処理であることを特徴とする逆数の近似計算方法。
【請求項１３】
浮動小数点数を処理する演算機能を有するプロセッサを備えた演算装置に記憶され、前記プロセッサにより実行される平方根逆数の近似計算プログラムであって、平方根逆数を求めるべき入力値に対応し所定の初期近似精度を有する平方根逆数である初期近似平方根逆数に対して、ニュートンラフソン法による近似計算を行って、予め設定された目標精度を有する目標平方根逆数を算出する処理を前記プロセッサに実行させる平方根逆数の近似計算プログラムにおいて、
前記初期近似平方根逆数に対して、ニュートンラフソン法による１段近似計算を少なくとも１回行って、前記初期近似精度と前記目標精度の間に設定された所定の中間近似精度を有する中間近似平方根逆数を算出する中間値演算処理と、
前記中間近似平方根逆数に対して、前記中間近似精度から前記目標精度まで精度が上がるようにニュートンラフソン法による多段近似計算を行って、前記目標平方根逆数を算出する目標値演算処理とを前記プロセッサに実行させ、
前記中間近似精度は、前記プロセッサにおける浮動小数点数の減算処理に伴う仮数部のビットシフトによるストールが生じない上限ビットシフト量に対応する非ストール誤差が残る精度に設定され、前記中間値演算処理は、前記非ストール誤差が近似計算結果に残るように誤差を拡大した誤差拡大１段近似計算を行うことにより、前記中間近似精度を越えないように近似精度を制限した前記中間近似平方根逆数を算出する処理であることを特徴とする平方根逆数の近似計算プログラム。
【請求項１４】
浮動小数点数を処理する演算機能を有するプロセッサを備えた演算装置にて実行される平方根逆数の近似計算方法であって、平方根逆数を求めるべき入力値に対応し所定の初期近似精度を有する平方根逆数である初期近似平方根逆数に対して、ニュートンラフソン法による近似計算を行って、予め設定された目標精度を有する目標平方根逆数を算出する処理を前記演算装置に実行させる平方根逆数の近似計算方法において、
前記初期近似平方根逆数に対して、ニュートンラフソン法による１段近似計算を少なくとも１回行って、前記初期近似精度と前記目標精度の間に設定された所定の中間近似精度を有する中間近似平方根逆数を算出する中間値演算処理と、
前記中間近似平方根逆数に対して、前記中間近似精度から前記目標精度まで精度が上がるようにニュートンラフソン法による多段近似計算を行って、前記目標平方根逆数を算出する目標値演算処理とを前記プロセッサに実行させ、
前記中間近似精度は、前記プロセッサにおける浮動小数点数の減算処理に伴う仮数部のビットシフトによるストールが生じない上限ビットシフト量に対応する非ストール誤差が残る精度に設定され、前記中間値演算処理は、前記非ストール誤差が近似計算結果に残るように誤差を拡大した誤差拡大１段近似計算を行うことにより、前記中間近似精度を越えないように近似精度を制限した前記中間近似平方根逆数を算出する処理であることを特徴とする平方根逆数の近似計算方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、逆数または平方根逆数の近似計算をニュートンラフソン法によって行う演算装置に関し、特に、近似計算にて同程度の大きさの浮動小数点数の減算処理に伴ってプロセッサで生じる長いビットシフト処理に起因するストールを回避できる演算装置に関する。
【背景技術】
【０００２】
逆数および平方根逆数（平方根の逆数、以下同じ）を高精度で計算する機能を持った演算装置は、各種の分野で有用に用いられている。例えば、分子生物学の分野では、クーロンポテンシャル（エネルギ）およびクーロン力を計算するのに、逆数および平方根逆数が用いられる。そして、膨大な情報を処理するために高い演算速度が求められる。
【０００３】
従来より、逆数および平方根逆数を計算するために、ニュートンラフソン法が一般に利用されている（例えば特許文献１参照）。以下、ニュートンラフソン法の原理について説明する。
【０００４】
ニュートンラフソン法は、ｆ（ｘ）＝０を満たすｘを求める手法であり、下の式（１）による近似計算が行われる。
ｘｉ＋１＝ｘｉ−ｆ（ｘｉ）／ｆ’（ｘｉ）・・・式（１）
ｘに十分近いｘｉに対して式（１）の計算を１ステップ行うと、精度が倍になる。そこで、目標精度に到達するまで式（１）の計算が繰り返される。
【０００５】
次に、ニュートンラフソン法によって逆数および平方根逆数を計算する方法を説明する。
「逆数」
数値ａの逆数１／ａを求める場合、ｆ（ｘ）＝ａ−１／ｘとする。
ｆ’（ｘ）＝１／ｘ２
ｘｉ＋１＝ｘｉ−（ａ−１／ｘｉ）／（１／ｘｉ２）
ｘｉ＋１＝ｘｉ×（２−ａ×ｘｉ）    ・・・式（２）
【０００６】
式（２）を繰り返すことにより目標精度の逆数が計算される。ここで、式（２）の計算は、ニュートンラフソン法の１ステップに相当する。そこで、式（２）の計算を「１段近似計算」と呼ぶ。
【０００７】
ニュートンラフソン法は以下のように展開される。εｉ＝１−ａ×ｘｉとおくと、ニュートンラフソン法の２ステップが、式（３）で示すように、差分のないかたちで表される。
ｘｉ＋２＝ｘｉ×（１＋εｉ×（１＋εｉ×（１＋εｉ）））  ・・・式（３）
【０００８】
式（３）は、ニュートンラフソン法の２ステップに相当するので、以下、「多段近似計算」と呼ぶ。式（３）をさらに変形し、（１＋εｉ）の項を同じかたちで挿入すると、ニュートンラフソン法のステップ数を増大するのと同じ効果が得られる。このような計算も、「多段近似計算」に含まれる。
【０００９】
「平方根逆数」
数値ａの平方根逆数１／ａ１／２を求める場合、ｆ（ｘ）＝ａ−１／ｘ２とする。
ｆ’（ｘ）＝２／ｘ３
ｘｉ＋１＝（３×ｘｉ−ａ×ｘｉ３）／２    ・・・式（４）
【００１０】
式（４）の計算がニュートンラフソン法の１ステップであり、すなわち１段近似計算である。次に、εｉ＝１−ａ×ｘ２とおくと、ニュートンラフソン法の２ステップが、式（５）で表される。式（５）が多段近似計算に相当する。
ｘｉ＋２＝ｘｉ（１＋（1/2）εｉ＋（3/8）εｉ２＋（5/16）εｉ３＋（1/16）εｉ４）                                                            ・・・式（５）
【特許文献１】特開平１０−５５２６６号公報
【発明の開示】
【発明が解決しようとする課題】
【００１１】
従来の演算装置では、逆数および平方根逆数が上述したニュートンラフソン法の原理に従って計算される。ニュートンラフソン法は、理論的には、逆数および平方根逆数を高精度に計算できる好適な手法である。しかし、実際のプロセッサでの計算にニュートンラフソン法を単純に適用すると、以下に説明するように、プロセッサでの計算時間が長大になる場合がある。
【００１２】
プロセッサでは周知のように数値が浮動小数点数で表現される。プロセッサとしては、例えば、ＩＢＭ（登録商標、以下、同じ）のＰＯＷＥＲ４およびＰＯＷＥＲ５などが挙げられる。また、浮動小数点は例えばＩＥＥＥ７５４に規定されている。図１１は、ＩＥＥＥ７５４の倍精度の浮動小数点数を表しており、この場合、数値が、符号、指数部および仮数部で構成される。「仮数部」は５２ビットの値であり、仮数「１.ｘｉｘ２ｘ３・・・ｘ５２」のうちの小数点以下の部分のみを表している。
【００１３】
ニュートンラフソン法にて逆数を計算する場合に、精度よく計算しようとすると、同程度の大きさの２つの数値の差分を計算する必要が生じてくる。しかしながら、浮動小数点数を扱うプロセッサにて同程度の２つの数値の差分を計算すると、プロセッサ内部で長いビットシフトが行われることになる。
【００１４】
例として、図１２に示すように、同程度の２つの数値が引き算され、両数値が図示の仮数部を有していたとする。この場合、仮数部が近いので、仮数部の差を見ると、前側に値０のビットが並ぶ。ＩＥＥＥ７５４の正規化数の表現では、図１１に示したように、仮数は、「１．ｘｉｘ２ｘ３・・・ｘ５２」と規定されており、すなわち、最初の桁は「１．」に固定されており、その後に小数点以下の５２ｂｉｔの値が続くように規定されている。そのため、仮数部の差が図１２のようになる場合、最初の「１」が現れるまで、すなわち、差分にて値１を持つ先頭のビットが仮数の前端に達するまで、ビットシフトを行う必要がある。２つの数値が近いと、仮数部の差分にならぶ値０のビットが増え、そのためにビットシフト量も大きくなる。
【００１５】
このように、ニュートンラフソン法にて逆数を高精度に計算しようとすると、同程度の大きさの数値の差分を計算することが必要になり、その結果、プロセッサ内部で長いビットシフトが行われてしまう。そして、長いビットシフトは、以降の浮動小数点演算命令の実行をストールさせてしまい、これにより計算速度が低下してしまう。
【００１６】
ここでは、逆数を求める場合に関して従来の問題を説明したが、平方根逆数を求める場合にも同様に計算速度の向上が望まれる。
【００１７】
本発明は上記背景の下でなされたものであり、その目的は、逆数または平方根逆数の近似計算をニュートンラフソン法によって行う際、同程度の大きさの浮動小数点数の減算処理に伴ってプロセッサで生じるビットシフト処理に起因するストールを回避でき、計算速度を向上できる演算装置を提供することにある。
【課題を解決するための手段】
【００１８】
本発明の演算装置は、逆数を求めるべき入力値に対応し所定の初期近似精度を有する逆数である初期近似逆数に対して、ニュートンラフソン法による近似計算を行って、予め設定された目標精度を有する目標逆数を算出する演算装置であって、浮動小数点数を処理する演算機能を有するプロセッサと、逆数の近似計算プログラムを記憶する記憶部とを有し、前記記憶部に記憶された前記近似計算プログラムを前記プロセッサにて実行することによって中間値演算部および目標値演算部が前記演算装置に備えられ、前記中間値演算部は、前記初期近似逆数に対して、ニュートンラフソン法による１段近似計算を少なくとも１回行って、前記初期近似精度と前記目標精度の間に設定された所定の中間近似精度を有する中間近似逆数を算出する中間値演算処理を行い、前記目標値演算部は、前記中間近似逆数に対して、前記中間近似精度から前記目標精度まで精度が上がるようにニュートンラフソン法による多段近似計算を行って、前記目標逆数を算出する目標値演算処理を行い、前記中間近似精度は、前記プロセッサにおける浮動小数点数の減算処理に伴う仮数部のビットシフトによるストールが生じない上限ビットシフト量に対応する非ストール誤差が残る精度に設定され、前記中間値演算部は、前記非ストール誤差が近似計算結果に残るように誤差を拡大した誤差拡大１段近似計算を行うことにより、前記中間近似精度を越えないように近似精度を制限した前記中間近似逆数を算出する。
【００１９】
上記のように、本発明によれば、中間近似精度を持つ中間近似逆数が誤差拡大１段近似計算により算出され、それから、目標精度を持つ目標逆数が中間近似逆数から多段近似計算により算出される。中間近似逆数では近似精度が中間近似精度に制限されており、プロセッサの浮動小数点数の減算処理に伴う仮数部のビットシフトによるストールが生じない上限ビットシフト量に対応する非ストール誤差が残っている。したがって、差が小さすぎる２つの値の減算処理を行うことなく、近似計算を行える。これにより、ストールを回避でき、計算速度を向上できる。
【００２０】
また、中間値演算部は、前記入力値をａとし、近似計算前後の近似逆数をｘｉ、ｘｉ＋１としたときに、前記誤差拡大１段近似計算として、ｘｉ＋１＝ｘｉ×（（２＋α）−a×ｘｉ）を求める処理を行い、誤差成分αが前記非ストール誤差に応じて設定されてよい。
【００２１】
また、前記プロセッサがＩＥＥＥ７５４により規定される浮動小数点数を処理し、前記初期近似精度が８ｂｉｔ精度であり、前記目標近似精度が５３ｂｉｔ精度であり、前記中間近似精度が１４ｂｉｔ精度であってよい。また、前記誤差成分αが、２−１４であってよい。
【００２２】
また、前記目標値演算部は、前記中間近似逆数をｘｉとし、εｉ＝１−ａ×ｘｉとしたときに、前記多段近似計算として、ｘｉ×（１＋εｉ×（１＋εｉ×（１＋εｉ×（１＋εｉ））））を求める処理を行ってよい。
【００２３】
また、本発明の別の態様に係る演算装置は平方根逆数を計算する。この演算装置は、浮動小数点数を処理する演算機能を有するプロセッサと、平方根逆数の近似計算プログラムを記憶する記憶部とを有し、前記記憶部に記憶された前記近似計算プログラムを前記プロセッサにて実行することにより、平方根逆数を求めるべき入力値に対応し所定の初期近似精度を有する平方根逆数である初期近似平方根逆数に対して、ニュートンラフソン法による近似計算を行って、予め設定された目標精度を有する目標平方根逆数を算出する演算装置であって、前記初期近似平方根逆数に対して、ニュートンラフソン法による１段近似計算を少なくとも１回行って、前記初期近似精度と前記目標精度の間に設定された所定の中間近似精度を有する中間近似平方根逆数を算出する中間値演算部と、前記中間近似平方根逆数に対して、前記中間近似精度から前記目標精度まで精度が上がるようにニュートンラフソン法による多段近似計算を行って、前記目標平方根逆数を算出する目標値演算部と、を備え、前記中間近似精度は、前記プロセッサにおける浮動小数点数の減算処理に伴う仮数部のビットシフトによるストールが生じない上限ビットシフト量に対応する非ストール誤差が残る精度に設定され、前記中間値演算部は、前記非ストール誤差が近似計算結果に残るように誤差を拡大した誤差拡大１段近似計算を行うことにより、前記中間近似精度を越えないように近似精度を制限した前記中間近似平方根逆数を算出する。
【００２４】
上記のように、本発明によれば、中間近似精度を持つ中間近似平方根逆数が誤差拡大１段近似計算により算出され、それから、目標精度を持つ目標平方根逆数が中間近似平方根逆数から多段近似計算により算出される。中間近似平方根逆数では近似精度が中間近似精度に制限されており、プロセッサの浮動小数点数の減算処理に伴う仮数部のビットシフトによるストールが生じない上限ビットシフト量に対応する非ストール誤差が残っている。したがって、差が小さすぎる２つの値の減算処理を行うことなく、近似計算を行える。これにより、ストールを回避でき、計算速度を向上できる。
【００２５】
また、前記中間値演算部は、誤差拡大を行わない１段近似計算を１回と前記誤差拡大１段近似計算を１回行うように構成され、前記入力値をａとし、近似計算前後の近似平方根逆数をｘｉ、ｘｉ＋１としたときに、誤差拡大を行わない１段近似計算として、ｘｉ＋１＝ｘｉ＋ｘｉ×（０．５−０．５×a×ｘｉ）を求める処理を行い、前記誤差拡大１段近似計算として、ｘｉ＋１＝ｘｉ＋ｘｉ×（（０．５＋β）−０．５×a×ｘｉ）を求める処理を行い、誤差成分βが前記非ストール誤差に応じて設定されてよい。
【００２６】
また、前記プロセッサがＩＥＥＥ７５４により規定される浮動小数点数を処理し、前記初期近似精度が５ｂｉｔ精度であり、前記目標近似精度が５３ｂｉｔ精度であり、前記中間近似精度が１４ｂｉｔ精度であってよい。また、前記誤差成分βが、２−１５であってよい。
【００２７】
また、前記目標値演算部は、前記中間近似逆数をｘｉとし、εｉ＝１−ａ×ｘｉ２としたときに、多段近似計算として、ｘｉ×（１＋（1/2）×εｉ＋（3/8）×εｉ２＋（5/16）×εｉ３）を求める処理を行ってよい。
【００２８】
本発明は上記の演算装置の態様に限定されない。本発明の別の態様は、例えば、下記に挙げるように、逆数の近似計算プログラムおよび近似計算方法であり、また、平方根逆数の近似計算プログラムおよび近似計算方法である。上記の演算装置の態様における各種の付加的な構成は、これらプログラムおよび方法の態様にも、各態様に適したかたちで適用されてよい。
【００２９】
本発明において、逆数の近似計算プログラムの態様は、浮動小数点数を処理する演算機能を有するプロセッサを備えた演算装置に記憶され、前記プロセッサにより実行される逆数の近似計算プログラムであって、逆数を求めるべき入力値に対応し所定の初期近似精度を有する逆数である初期近似逆数に対して、ニュートンラフソン法による近似計算を行って、予め設定された目標精度を有する目標逆数を算出する処理を前記プロセッサに実行させる逆数の近似計算プログラムにおいて、前記初期近似逆数に対して、ニュートンラフソン法による１段近似計算を少なくとも１回行って、前記初期近似精度と前記目標精度の間に設定された所定の中間近似精度を有する中間近似逆数を算出する中間値演算処理と、前記中間近似逆数に対して、前記中間近似精度から前記目標精度まで精度が上がるようにニュートンラフソン法による多段近似計算を行って、前記目標逆数を算出する目標値演算処理とを前記プロセッサに実行させ、前記中間近似精度は、前記プロセッサにおける浮動小数点数の減算処理に伴う仮数部のビットシフトによるストールが生じない上限ビットシフト量に対応する非ストール誤差が残る精度に設定され、前記中間値演算処理は、前記非ストール誤差が近似計算結果に残るように誤差を拡大した誤差拡大１段近似計算を行うことにより、前記中間近似精度を越えないように近似精度を制限した前記中間近似逆数を算出する処理である。
【００３０】
また、本発明における逆数の近似計算方法の態様は、浮動小数点数を処理する演算機能を有するプロセッサを備えた演算装置にて実行される逆数の近似計算方法であって、逆数を求めるべき入力値に対応し所定の初期近似精度を有する逆数である初期近似逆数に対して、ニュートンラフソン法による近似計算を行って、予め設定された目標精度を有する目標逆数を算出する処理を前記演算装置に実行させる逆数の近似計算方法において、前記初期近似逆数に対して、ニュートンラフソン法による１段近似計算を少なくとも１回行って、前記初期近似精度と前記目標精度の間に設定された所定の中間近似精度を有する中間近似逆数を算出する中間値演算処理と、前記中間近似逆数に対して、前記中間近似精度から前記目標精度まで精度が上がるようにニュートンラフソン法による多段近似計算を行って、前記目標逆数を算出する目標値演算処理とを前記プロセッサに実行させ、前記中間近似精度は、前記プロセッサにおける浮動小数点数の減算処理に伴う仮数部のビットシフトによるストールが生じない上限ビットシフト量に対応する非ストール誤差が残る精度に設定され、前記中間値演算処理は、前記非ストール誤差が近似計算結果に残るように誤差を拡大した誤差拡大１段近似計算を行うことにより、前記中間近似精度を越えないように近似精度を制限した前記中間近似逆数を算出する処理である。
【００３１】
また、本発明における平方根逆数の近似計算プログラムの態様は、浮動小数点数を処理する演算機能を有するプロセッサを備えた演算装置に記憶され、前記プロセッサにより実行される平方根逆数の近似計算プログラムであって、平方根逆数を求めるべき入力値に対応し所定の初期近似精度を有する平方根逆数である初期近似平方根逆数に対して、ニュートンラフソン法による近似計算を行って、予め設定された目標精度を有する目標平方根逆数を算出する処理を前記プロセッサに実行させる平方根逆数の近似計算プログラムにおいて、前記初期近似平方根逆数に対して、ニュートンラフソン法による１段近似計算を少なくとも１回行って、前記初期近似精度と前記目標精度の間に設定された所定の中間近似精度を有する中間近似平方根逆数を算出する中間値演算処理と、前記中間近似平方根逆数に対して、前記中間近似精度から前記目標精度まで精度が上がるようにニュートンラフソン法による多段近似計算を行って、前記目標平方根逆数を算出する目標値演算処理とを前記プロセッサに実行させ、前記中間近似精度は、前記プロセッサにおける浮動小数点数の減算処理に伴う仮数部のビットシフトによるストールが生じない上限ビットシフト量に対応する非ストール誤差が残る精度に設定され、前記中間値演算処理は、前記非ストール誤差が近似計算結果に残るように誤差を拡大した誤差拡大１段近似計算を行うことにより、前記中間近似精度を越えないように近似精度を制限した前記中間近似平方根逆数を算出する処理である。
【００３２】
また、本発明における平方根逆数の近似計算方法の態様は、浮動小数点数を処理する演算機能を有するプロセッサを備えた演算装置にて実行される平方根逆数の近似計算方法であって、平方根逆数を求めるべき入力値に対応し所定の初期近似精度を有する平方根逆数である初期近似平方根逆数に対して、ニュートンラフソン法による近似計算を行って、予め設定された目標精度を有する目標平方根逆数を算出する処理を前記演算装置に実行させる平方根逆数の近似計算方法において、前記初期近似平方根逆数に対して、ニュートンラフソン法による１段近似計算を少なくとも１回行って、前記初期近似精度と前記目標精度の間に設定された所定の中間近似精度を有する中間近似平方根逆数を算出する中間値演算処理と、前記中間近似平方根逆数に対して、前記中間近似精度から前記目標精度まで精度が上がるようにニュートンラフソン法による多段近似計算を行って、前記目標平方根逆数を算出する目標値演算処理とを前記プロセッサに実行させ、前記中間近似精度は、前記プロセッサにおける浮動小数点数の減算処理に伴う仮数部のビットシフトによるストールが生じない上限ビットシフト量に対応する非ストール誤差が残る精度に設定され、前記中間値演算処理は、前記非ストール誤差が近似計算結果に残るように誤差を拡大した誤差拡大１段近似計算を行うことにより、前記中間近似精度を越えないように近似精度を制限した前記中間近似平方根逆数を算出する処理である。
【発明の効果】
【００３３】
上記のように、本発明は、逆数および平方根逆数の近似計算におけるストールを回避でき、計算を高速に行うことができる演算装置、近似計算プログラムおよび方法を提供できる。
【発明を実施するための最良の形態】
【００３４】
以下、本発明の実施の形態に係る演算装置について、図面を用いて説明する。
【００３５】
「第１の実施の形態」
図１は、第１の実施の形態に係る演算装置を示している。第１の実施の形態では、逆数が計算される。
【００３６】
図１において、演算装置１は、コンピュータであり、プロセッサ３と記憶部５とを備えている。プロセッサ３は、浮動小数点数を処理する四則演算機能を有している。例えばプロセッサ３がＩＢＭのＰＯＷＥＲ５であり、浮動小数点数はＩＥＥＥ７５４の倍精度（図１１参照）で表される。記憶部５はメモリであり、本発明に係る逆数の近似計算プログラムを記憶している。記憶部５に記憶された近似計算プログラムがプロセッサ３で実行される。これにより、プロセッサ３に入力値ａが入力され、入力値ａの逆数１／ａが計算されて、逆数１／ａがプロセッサ３から出力される。
【００３７】
図２は、演算装置１の構成を示す機能ブロック図である。図２の各構成は、図１のプロセッサ３および記憶部５により、すなわち、記憶部５の近似計算プログラムをプロセッサ３にて実行することにより実現される。
【００３８】
図２に示すように、演算装置１は、入力部１１、初期値取得部１３、中間値演算部１５、目標値演算部１７および出力部１９で構成されている。入力部１１は、逆数を求めるべき入力値ａを入力する。
【００３９】
初期値取得部１３は、入力値ａに対応する初期近似逆数を取得する。初期近似逆数は、所定の初期近似精度を有する逆数であり、初期近似精度は比較的低い値でよい。初期値取得部１３としては典型的には以下の構成が考えられる。初期値取得部１３としては、逆数テーブル（数値とその逆数を関連づけるテーブル）が予め記憶部５に記憶される。そして、入力値ａに対応する逆数がテーブルから読み出される。精度が低くてよいので、逆数テーブルは小さくてよい。入力値から、逆数テーブルの精度に対応する上位桁部分が取り出され、上位桁部分に対応する逆数がテーブルから読み出されてよい。また、この初期値取得部１３の機能の一部として、逆数テーブルから読み出された低精度の逆数に近似計算が施されて、精度がある程度まで高められてもよい（ここでの近似計算がニュートンラフソン法であってもよい）。
【００４０】
初期値取得部１３は、プロセッサ３に一般に備えられた精度が比較的低い逆数取得機能により好適に実現される。一般のプロセッサには、このような基本的な逆数取得機能についてのコマンドが設定されている。このコマンドが逆数計算のプログラムに組み込まれ、同コマンドに応答して逆数が呼び出されてよい。つまり、プログラム中の逆数呼出のコマンドとそれに対応するプロセッサの逆数呼出機能により初期値取得部１３が好適に実現される。例えば、本実施の形態の例として取り上げているＩＢＭのＰＯＷＥＲ５では、コマンドｆｒｅ（ａ）によって、８ｂｉｔ精度の逆数が呼び出される。浮動小数点の仮数が５３桁であることを考慮すると、初期の逆数の精度は低いといえる。
【００４１】
上記のようにして本実施の形態の演算装置１では入力値ａに対応する初期近似逆数が取得される。初期近似逆数が有する初期近似精度は上述したように比較的低い。この初期近似逆数を対象として本実施の形態の演算装置１では下記の構成、すなわち、中間値演算部１５および目標値演算部１７によりニュートンラフソン法の近似計算が行われ、近似精度が所定の目標精度まで高められる。
【００４２】
中間値演算部１５は、初期値取得部１３で取得された初期近似逆数に対して、ニュートンラフソン法による１段近似計算を少なくとも１回行って、所定の中間近似精度を有する中間近似逆数を算出する。中間近似精度は、初期近似精度と目標精度の間の精度に設定されており、そして、以下に説明するように、プロセッサでのストールを回避できる適切な精度に設定されている。
【００４３】
プロセッサでは、図１１および図１２を用いて説明したように、同程度の２つの数値の減算処理が行われると、仮数部のビットシフトが発生する。減算対象の２つの数値が近いほど、ビットシフト量が大きくなる。したがって、本実施の形態のような演算装置１の場合、中間近似精度が高すぎると、次の目標値演算部１７での近似計算における減算処理にてビットシフト量が大きくなりすぎて、以降の浮動小数点演算命令の実行をストールさせてしまう。ストールが生じると、プロセッサでは、数サイクルの間、命令が発行されなくなる。
【００４４】
そこで、中間近似精度は、上限ビットシフト量に対応する誤差（以下、非ストール誤差という）が残る精度に設定されている。上限ビットシフト量は、上記のビットシフトによるストールが生じない範囲に予め設定された所定のビットシフト量である。本実施の形態では、後述にて具体例を用いて示すように、ストール発生の閾になるビットシフト量が特定される。閾値よりビットシフト量が大きいとストールが発生する。そこで、ストールを確実に防ぐために、閾値より更に１ビット小さいシフト量が上限ビットシフト量に定められる。非ストール誤差は、入力値ａの逆数１／ａ（真の逆数）に対する相対的な誤差であり、相対誤差ということができる。中間値演算部１５は、近似計算結果に非ストール誤差を残すために、通常の１段近似計算ではなく、本実施の形態の特徴である誤差拡大１段近似計算を行う。誤差拡大１段近似計算とは、非ストール誤差が近似計算結果に残るように故意に誤差を拡大したニュートンラフソン法の１ステップの近似計算である。これにより、中間値演算部１５は、中間近似精度を越えないように近似精度を制限した中間近似逆数を算出できる。
【００４５】
目標値演算部１７は、中間値演算部で求められた中間近似逆数に対して、中間近似精度から目標精度まで精度が上がるようにニュートンラフソン法による多段近似計算を行って、目標逆数を算出する。目標逆数は、目標精度を有する逆数である。この目標逆数が、計算結果の逆数１／ａとして出力部１９から出力される。
【００４６】
次に、図３を参照し、上記の中間値演算部１５および目標値演算部１７による処理について、具体例を用いてさらに説明する。この具体例では、プロセッサ３が前出のＰＯＷＥＲ５であり、浮動小数点はＩＥＥＥ７５４に従うとする。この場合、プロセッサ３の標準装備の逆数呼出機能により、８ｂｉｔ精度の逆数が得られる。すなわち、初期近似精度は８ｂｉｔである。また、浮動小数点の精度は５３ｂｉｔであり、そこで、目標精度は、浮動小数点で表現できる最大限の精度である５３ｂｉｔとする。
【００４７】
また、上述した本実施の形態で例示するプロセッサ３では、ストールが生じるビットシフト量の閾値が、１５ビットである。それより大きいビットシフト量が減算処理で生じると、ストールが発生する。より詳細には、数サイクルの期間は新しい命令が発行されなくなるといった種類のストールが発生し、これにより性能が低下する。そこで、本実施の形態では、ストールを確実に防ぐために、上限ビットシフト量が、閾値（１５ビット）より１ビット少ない１４ビットに設定される。そして、上限ビットシフト量に対応する非ストール誤差が残るように、上限ビットシフト量の１４ビットに対応して、中間近似精度が１４ビットに設定される。
【００４８】
中間値演算部１５は、上記の中間近似精度を達成するために、下記の式（６）の誤差拡大１段近似計算を行う。
ｘｉ＋１＝ｘｉ×（（２＋α）−ａ×ｘｉ）    ・・・式（６）
上記の式（６）を通常のニュートンラフソン法の１段近似計算（式（２））と比較すると、誤差成分αが挿入されている。α＝２−１４＝０．００００６１０３５１５６２５とする。この誤差成分αを挿入することにより、誤差が故意に拡大されて、１４ｂｉｔ精度の近似結果が得られる。
【００４９】
次に、目標値演算部１７では、下記の式（７）の多段近似計算が行われる。
ｘｉ＋２＝ｘｉ×（１＋εｉ×（１＋εｉ×（１＋εｉ×（１＋εｉ））））
εｉ＝１−ａ×ｘｉ                                      ・・・式（７）
この多段近似計算式により、１４ｂｉｔ精度の中間近似逆数から、５３ｂｉｔ精度の目標逆数が得られる。
【００５０】
上記の中間値演算部１５および目標値演算部１７の処理についてさらに説明する。通常のニュートンラフソン法の１段近似計算（式（２））が行われると、近似精度は倍になる。この例では、初期近似精度が８ｂｉｔなので、通常の１段近似計算では精度が１６ｂｉｔになる。これでは、精度が高くなりすぎる。すなわち、次の多段近似計算（式（７））にて、εｉ（＝１−ａ×ｘｉ）を計算するときに、大きなビットシフトが生じ、ストールが発生する。
【００５１】
一方、本発明では、式（６）の誤差拡大１段近似計算にて、誤差成分α＝２−１４が挿入されている。したがって、中間近似逆数の精度が、１６ｂｉｔまで高くならず、１４ｂｉｔで留まる。これにより、次のステップにてεｉ（＝１−ａ×ｘｉ）を計算するときに、多大なビットシフトが発生せず、ストールが回避される。
【００５２】
また、式（７）の多段近似計算を、前出の式（３）と比較すると、多少異なっていることが分かる。式（３）は、ニュートンラフソン法の２ステップ分の近似計算の基本的な式である。式（３）と式（７）を比較すると、式（７）では、（１＋εｉ）の項が１つ多く挿入されている。したがって、式（７）は、２ステップよりも大きい近似計算、すなわち、２ステップ分より精度が上がる近似計算を行っていることを意味する。式（７）をこのように設定している理由は下記の通りである。
【００５３】
既に説明したように、中間近似精度は１４ｂｉｔ精度に達している。ニュートンラフソン法の１ステップでは精度が倍であり、２ステップでは精度が４倍である。したがって、２ステップの近似計算を行えば、十分に目標精度（＝５３ｂｉｔ）が得られると考えられる。それにも拘わらず、本実施の形態では、より高い精度が得られるように、上記の如く多段近似計算が調整されており、（１＋εｉ）の項が増やされている。これは、確実に目標精度が得られるようにするためである。プロセッサ３にて上記のような計算を行うと、プロセッサ内部の機構により精度がばらつくことがあり、通常の２ステップ分の近似計算を採用したのでは、検算を実施すると５３ｂｉｔが得られないことがある。この点を考慮し、本実施の形態では、多段近似計算が調整されて、精度が高められている。なお、上記のように本実施の形態では多段近似計算の精度を２ステップ分より上げているが、前出の式（７）では式の左辺をｘｉ＋２と記載し、表現を前出の式（３）と揃えている。
【００５４】
以上に、演算装置１の構成について説明した。次に、図４を参照し、演算装置１の動作を説明する。図４のフローチャートは、本実施の形態による逆数の近似計算方法を示すものであり、また、同近似計算方法をプロセッサ３に実行させる逆数の近似計算プログラムを示している。
【００５５】
図４に示すように、演算装置１では、入力部１１により、逆数を求めるべき入力値ａが入力され（Ｓ１）、初期値取得部１３により入力値ａに対応する初期近似逆数が取得される（Ｓ３）。初期近似逆数は、比較的低い初期近似精度（８ｂｉｔ精度）を有している。
【００５６】
そして、中間値演算部１５にて、初期近似逆数に対して、誤差拡大１段近似計算（式（６））が施されて、中間近似逆数が算出される（Ｓ５）。中間近似逆数は、中間近似精度（１４ｂｉｔ精度）を有している。さらに、目標値演算部１７にて、中間近似逆数に対して、多段近似計算（式（７））が施されて、目標逆数が計算される（Ｓ７）。目標逆数は目標精度（５３ｂｉｔ精度）を有している。目標逆数が、最終的な逆数１／ａとして出力部１９から出力される（Ｓ９）。
【００５７】
図５は、図４のフローチャートに対応する具体的なプログラムの好適な例を示している。図５のステップＳ１１（ｘ＝ｆｒｅ（ａ））は、図４のステップ３に対応し、初期近似逆数を取得するステップである。「ｆｒｅ（ａ）」は、ａの逆数を呼び出すコマンドであり、このコマンドに応答して８ｂｉｔ精度の逆数が返ってくる。なお、同コマンドは、ＩＢＭのプロセッサ「ＰＯＷＥＲ５」に設定されている。同「ＰＯＷＥＲ４」での同様のコマンドは、「ｆｒｅｓ（ａ）」である。このようにコマンドは適用対象のプロセッサに従うことはもちろんである。
【００５８】
次に、図５のステップＳ１３は、図４のステップ５に対応し、式（６）の誤差拡大１段近似計算に対応している。図示のように、誤差成分２−１４（＝０．００００６１０３５１５６２５）が２に加算されており、１４ｂｉｔ精度の解（逆数）が得られる。
【００５９】
また、図５のステップＳ１５は、図４のステップＳ７に対応し、式（７）の多段近似計算に対応している。ステップＳ１５の１行目で、ｙ＝１．０−ｘ×ａが計算される。ｘは中間近似逆数であり、ｙは式（７）のεｉに相当する。この計算は図示のようにちょうど誤差計算に相当している。ｘの精度が良すぎると、ｘ×ａが１に近くなりすぎて、大きなビットシフトとそれによるストールが発生する。これに対して、本実施の形態では、ｘの精度が１４ｂｉｔに抑えられているので、ストールが回避される。したがって、図５の８個の計算命令は、パイプラインストールを起こすことなく実行される。
【００６０】
図６および図７は、比較例として別のプログラムを示している。これらプログラムは、通常のニュートンラフソン法による単純な実装を行った場合に相当する。図６の例では、通常の１段近似計算が３回行われる。この場合、７つの命令で計算がすむ。しかし、ｘ＝ｘ×ｙという積の計算精度にプロセッサ上での制約があり、最終の精度が本実施の形態の手法に劣る。また、図７の例では、精度は高いが、同程度の数値の減算が行われて、パイプラインストールが生じるため、性能が悪く、本実施の形態の手法と比べると計算速度が約半分になる。これらと比べると、図５に示す本発明の手法は、高精度かつ高速に逆数を算出することができる。
【００６１】
以上に本発明の第１の実施の形態について説明した。上記のように、本発明によれば、中間近似精度を持つ中間近似逆数が誤差拡大１段近似計算により算出され、それから、目標精度を持つ目標逆数が中間近似逆数から多段近似計算により算出される。中間近似逆数では近似精度が中間近似精度に制限されており、プロセッサの浮動小数点数の減算処理に伴う仮数部のビットシフトによるストールが生じない上限ビットシフト量に対応する非ストール誤差が残っている。したがって、差が小さすぎる２つの値の減算処理を行うことなく、近似計算を行える。これにより、ストールを回避でき、計算速度を向上できる。
【００６２】
また、本発明によれば、入力値をａとし、近似計算前後の近似逆数をｘｉ、ｘｉ＋１としたときに、誤差拡大１段近似計算として、式（６）に示したように、ｘｉ＋１＝ｘｉ×（（２＋α）−a×ｘｉ）が求められ、誤差成分αが非ストール誤差に応じて設定されている。これにより、非ストール誤差が近似計算結果に残る誤差拡大１段近似計算を好適に実現できる。
【００６３】
また、本発明によれば、プロセッサにより、ＩＥＥＥ７５４により規定される浮動小数点数が処理される。初期近似精度が８ｂｉｔ精度であり、目標近似精度が５３ｂｉｔ精度であり、中間近似精度が１４ｂｉｔ精度である。これにより、ＩＥＥＥ７５４により規定される浮動小数点数を処理するプロセッサにおいて、ストールを回避できる範囲で中間近似精度を高くすることができ、その後の多段近似計算をできるだけ短くでき、計算速度を好適に向上できる。
【００６４】
ここで、本発明を適用せず、１段近似計算は行わないとする。そして、最初から、初期近似逆数を出発点として、目標精度まで達するように、多段近似計算を行ったとする。このような処理も理論的には可能かもしれない。また、比較的低い初期近似精度の逆数を用いた減算処理が行われればよいので、本発明が問題とするような減算での大きなビットシフトも生じない。しかしながら、この場合、多段近似計算で出発点の精度が低いことは相当に不利である。十分な精度を得るためには、相当な数の（１＋εｉ）の項を、多段近似計算の式に挿入しなければならず、そのために計算時間が長くなってしまう。このような構成と比較しても、本発明は逆数を高速に計算でき、有利である。
【００６５】
また、本発明によれば、誤差成分αが、２−１４である。これにより、非ストール誤差を残した中間近似逆数を好適に得られる。
【００６６】
また、本発明によれば、中間近似逆数をｘｉとし、式（７）に示したように、εｉ＝１−ａ×ｘｉとしたときに、多段近似計算として、ｘｉ×（１＋εｉ×（１＋εｉ×（１＋εｉ×（１＋εｉ））））が求められる。これにより、目標精度に到達できる多段近似計算を好適に行える。
【００６７】
「第２の実施の形態」
次に、本発明に係る第２の実施の形態について説明する。第２の実施の形態に係る演算装置は、第１の実施の形態と同様にプロセッサおよび記憶部を備えている。上述した第１の実施の形態では、逆数が計算されたのに対して、第２の実施の形態では、平方根の逆数（以下、平方根逆数という）が算出される。すなわち、平方根逆数の近似計算プログラムが記憶部に記憶され、同プログラムがプロセッサにて実行されて、入力値ａから平方根逆数１／ａ１／２が計算される。以下の説明では、第１の実施の形態と共通する事項の説明は適当に省略する。
【００６８】
図８は、第２の実施の形態における演算装置２１を示している。第１の実施の形態と同様、本実施の形態でも、演算装置２１は、入力部３１、初期値取得部３３、中間値演算部３５、目標値演算部３７および出力部３９で構成されている。ただし、初期値取得部３３、中間値演算部３５および目標値演算部３７の機能は、対応する第１の実施の形態の構成とは異なっている。
【００６９】
初期値取得部３３は、入力値ａに対応する初期近似平方根逆数を取得する。初期値取得部３３は、第１の実施の形態と同様に、平方根逆数のテーブルを備えてよい。さらに、初期値取得部３３は、第１の実施の形態と同様に、プロセッサ３に予め備えられた平方根逆数の取得機能によって好適に実現される。本実施の形態で取り上げるＩＢＭのＰＯＷＥＲ５の場合、平方根逆数を読み出すコマンドは、ｆｒｓｑｒｔｅ（ａ）であり、このコマンドを平方根逆数の計算プログラムに組み込むことにより初期値取得部３３が実現される。そして、同プロセッサでは、５ｂｉｔ精度の平方根逆数が得られ、すなわち、初期近似精度が５ｂｉｔ精度になる。
【００７０】
次に、中間値演算部３５について説明する。中間値演算部３５は、初期近似平方根逆数から、中間近似精度を持つ中間近似平方根逆数を計算する。中間値演算部３５では、第１の実施の形態と同様に、ニュートンラフソン法の１段近似計算が行われる。ただし、第１の実施の形態では、１回の１段近似計算が行われ、この１回の１段近似計算が、誤差が挿入された誤差拡大１段近似計算であった。これに対して、第２の実施の形態では、２回の１段近似計算が行われる。そして、１回目は、通常の１段近似計算が行われ、２回目は誤差拡大近似計算が行われる。
【００７１】
中間値演算部３５の処理について、第１の実施の形態と同様に具体例を用いて説明する。１回目の通常の１段近似計算は、下の式（８）で表される。
ｘｉ＋１＝ｘｉ＋ｘｉ×（０．５−０．５×a×ｘｉ）    ・・・式（８）
【００７２】
式（８）は、前出の式（４）と実質的に同じであり、つまり、前出のニュートンラフソン法において平方根逆数を求めるための通常の１段近似計算に相当する。これに対して、２回目には、下の式（９）で表される誤差拡大１段近似計算が行われる。
ｘｉ＋１＝ｘｉ＋ｘｉ×（（０．５＋β）−０．５×a×ｘｉ）    ・・・式（９）
この２回目の誤差拡大１段近似計算で、中間近似精度が得られるように、誤差成分βが設定されている。具体的には、初期近似精度は上述の通り５ｂｉｔである。中間近似精度は、第１の実施の形態と同様に、ストールを回避するために１４ｂｉｔ精度に設定される。この場合に、β＝２−１５＝０．００００３０５１７６が好適である。この設定により、まず、１回目の通常の１段近似計算（式（８））にて近似精度が５ｂｉｔから１０ｂｉｔに上がり、２回目の誤差拡大１段近似計算にて近似精度が１０ｂｉｔから１４ｂｉｔに上がる。
【００７３】
次に、目標値演算部３７について説明する。目標値演算部３７は、中間近似平方根逆数に対して、中間近似精度から目標精度まで精度が上がるようにニュートンラフソン法による多段近似計算を行って、目標平方根逆数を算出する。この場合、多段近似計算として、下の式（１０）の計算が好適に行われる。式（１０）で算出された目標平方根逆数が、計算結果の平方根逆数１／ａ１／２として出力部３９から出力される。
ｘｉ＋２＝ｘｉ（１＋（1/2）εｉ＋（3/8）εｉ２＋（5/16）εｉ３）
εｉ＝１−ａ×ｘｉ２                                        ・・・式（１０）
【００７４】
式（１０）を、前出の式（５）と比較する。式（５）は、ニュートンラフソン法の原理に沿った２ステップ分の近似計算であった。式（５）と比べると、式（１０）では、３次（係数５／１６）の項までの計算が行われる。４次の項は省略されている。４次の項が計算されなくても、目標精度（５３ｂｉｔ精度）が達成されるからである。
【００７５】
式（１０）において、ｘｉが中間近似平方根逆数である。ｘｉの精度が高すぎた場合、εｉ（＝１−ａ×ｘｉ２）を計算するときに、大きなビットシフトが生じ、ストールが発生する。これに対して、本実施の形態では、式（９）の誤差拡大１段近似計算にてｘｉの精度が抑えられており、これによりストールが回避される。
【００７６】
以上に、演算装置２１の構成について説明した。次に、図９を参照し、演算装置２１の動作を説明する。図９のフローチャートは、本実施の形態による平方根逆数の近似計算方法を示すものであり、また、同近似計算方法をプロセッサ３に実行させる平方根逆数の近似計算プログラムを示している。
【００７７】
図９に示すように、演算装置１では、入力部３１により、平方根逆数を求めるべき入力値ａが入力され（Ｓ２１）、初期値取得部３３により入力値ａに対応する初期近似平方根逆数が取得される（Ｓ２３）。初期近似平方根逆数は、比較的低い初期近似精度（５ｂｉｔ精度）を有している。
【００７８】
そして、中間値演算部３５にて、初期近似平方根逆数に対して、１回の１段近似計算（式（８））と、１回の誤差拡大１段近似計算（式（９））が順次施されて、中間近似平方根逆数が算出される（Ｓ２５）。中間近似平方根逆数は、中間近似精度（１４ｂｉｔ精度）を有している。さらに、目標値演算部３７にて、中間近似平方根逆数に対して、多段近似計算（式（１０））が施されて、目標精度（５３ｂｉｔ精度）を有する目標平方根逆数が計算される（Ｓ２７）。目標平方根逆数が、最終的な平方根逆数１／ａ１／２として出力部３９から出力される（Ｓ２９）。
【００７９】
図１０は、図９のフローチャートに対応する具体的なプログラムの好適な例を示している。図１０のステップＳ３１（ｘ＝ｆｒｓｑｒｔｅ（ａ））は、図９のステップ２３に対応し、初期近似平方根逆数を取得するステップである。「ｆｒｓｑｒｔｅ（ａ）」は、ａの平方根逆数を呼び出すコマンドであり、このコマンドに応答して５ｂｉｔ精度の平方根逆数が返ってくる。なお、同コマンドは、ＩＢＭのプロセッサ「ＰＯＷＥＲ５」に設定されている。第１の実施の形態と同様、適用対象のプロセッサに従うコマンドが用いられることはもちろんである。
【００８０】
次に、図１０のステップＳ３３、Ｓ３５は、図９のステップＳ２５に対応し、式（８）の通常の１段近似計算と式（９）の誤差拡大１段近似計算にそれぞれ対応している。ステップＳ３５では、図示のように、誤差成分２−１５（＝０．００００３０５１７６）が０．５に加算されている。
【００８１】
また、図１０のステップＳ３７は、図９のステップＳ２７に対応し、式（１０）の多段近似計算に対応している。ステップＳ３７の１行目で、ｙ＝ｘ×ｘが計算され、２行目で、ｙ＝０．５−ａ×ｙが計算される。また、前出のステップＳ３３の最初に、ａ×＝０．５と記述され、ａ＝０．５×ａになっている。したがって、２行目の計算は、ｙ＝０．５−０．５×ａ×ｘ２＝０．５（１−ａ×ｘ２）であり、ここでεｉ（＝１−ａ×ｘｉ２）が計算される。そして、この計算の際に、ｘの精度が予め抑えられているので、プロセッサ３でのパイプラインストールの発生を回避することができる。
【００８２】
また、図１０のプログラム例では、ステップＳ３７の最後に、ｘ＝（ｆｌｏａｔ）ｘという行が付加されており、ｘ（中間近似平方根逆数）の下位に０クリア処理を行って、値が丸められる。このステップは、次のステップ（多段近似計算）においてεｉを計算する際の精度を上げるために設けられている。この（ｆｌｏａｔ）の処理により、２４ｂｉｔ以下の値が０になる。その結果、次のｘ２の計算精度が上がり、その次のεｉの計算精度が上がる。εｉの計算精度を上げることで、最終の計算精度を上げることができる。このような理由で、ステップ３７に（ｆｌｏａｔ）の処理が追加されている。
【００８３】
以上に本発明の第２の実施の形態について説明した。上記のように、本発明によれば、中間近似精度を持つ中間近似平方根逆数が誤差拡大１段近似計算により算出され、それから、目標精度を持つ目標平方根逆数が中間近似平方根逆数から多段近似計算により算出される。中間近似平方根逆数では近似精度が中間近似精度に制限されており、プロセッサの浮動小数点数の減算処理に伴う仮数部のビットシフトによるストールが生じない上限ビットシフト量に対応する非ストール誤差が残っている。したがって、差が小さすぎる２つの値の減算処理を行うことなく、近似計算を行える。これにより、ストールを回避でき、計算速度を向上できる。
【００８４】
また、本発明によれば、中間値演算処理では、誤差拡大を行わない１段近似計算が１回行われ、そして、誤差拡大１段近似計算が１回行われる。入力値をａとし、近似計算前後の近似平方根逆数をｘｉ、ｘｉ＋１としたときに、誤差拡大を行わない１段近似計算として、ｘｉ＋１＝ｘｉ＋ｘｉ×（０．５−０．５×a×ｘｉ）が求められ、誤差拡大１段近似計算として、ｘｉ＋１＝ｘｉ＋ｘｉ×（（０．５＋β）−０．５×a×ｘｉ）が求められ、誤差成分βが非ストール誤差に応じて設定されている。これにより、非ストール誤差が近似計算結果に残る誤差拡大１段近似計算を好適に実現できる。
【００８５】
また、本発明によれば、プロセッサにより、ＩＥＥＥ７５４により規定される浮動小数点数が処理される。初期近似精度が５ｂｉｔ精度であり、目標近似精度が５３ｂｉｔ精度であり、中間近似精度が１４ｂｉｔ精度である。これにより、ＩＥＥＥ７５４により規定される浮動小数点数を処理するプロセッサにおいて、ストールを回避できる範囲で中間近似精度を高くすることができ、その後の多段近似計算をできるだけ短くでき、計算速度を好適に向上できる。
【００８６】
また、本発明によれば、誤差成分βが、２−１５である。これにより、非ストール誤差を残した中間近似平方根逆数を好適に得られる。
【００８７】
また、本発明によれば、中間近似平方根逆数をｘｉとし、式（１０）に示したように、εｉ＝１−ａ×ｘｉ２としたときに、多段近似計算として、ｘｉ×（１＋（1/2）×εｉ＋（3/8）×εｉ２＋（5/16）×εｉ３）が求められる。これにより、目標精度に到達できる多段近似計算を好適に行える。
【００８８】
以上に本発明の好適な実施の形態を説明した。しかし、本発明は上述の実施の形態に限定されず、当業者が本発明の範囲内で上述の実施の形態を変形可能なことはもちろんである。
【産業上の利用可能性】
【００８９】
以上のように、本発明にかかる演算装置は、プロセッサにおける逆数および平方根逆数の計算を高速化できるという効果を有し、バイオインフォマティクスにおけるクーロン力やクーロンポテンシャルの計算をはじめとして、各種の分野における各種の計算において有用である。
【図面の簡単な説明】
【００９０】
【図１】本発明の第１の実施の形態に係る演算装置を示す図である。
【図２】第１の実施の形態に係る演算装置の構成を示す機能ブロック図である。
【図３】演算装置の中間値演算部および目標値演算部による処理について具体例を用いて説明する図である。
【図４】第１の実施の形態に係る演算装置の動作を示すフローチャートである。
【図５】図４のフローチャートに対応する具体的なプログラムの好適な例を示す図である。
【図６】図５に対する比較例としての別のプログラムを示す図である。
【図７】図５に対する比較例としての別のプログラムを示す図である。
【図８】本発明の第２の実施の形態に係る演算装置を示す機能ブロック図である。
【図９】第２の実施の形態に係る演算装置の動作を示すフローチャートである。
【図１０】図９のフローチャートに対応する具体的なプログラムの好適な例を示す図である。
【図１１】ＩＥＥＥ７５４に規定される浮動小数点数を示す図である。
【図１２】同程度の２つの浮動小数点数を減算したときにプロセッサで発生するビットシフトを示す図である。
【符号の説明】
【００９１】
１  演算装置
３  プロセッサ
５  記憶部
１１  入力部
１３  初期値取得部
１５  中間値演算部
１７  目標値演算部
１９  出力部
【図１】
2008287555.tif 000003
【図２】
2008287555.tif 000004
【図３】
2008287555.tif 000005
【図４】
2008287555.tif 000006
【図５】
2008287555.tif 000007
【図６】
2008287555.tif 000008
【図７】
2008287555.tif 000009
【図８】
2008287555.tif 000010
【図９】
2008287555.tif 000011
【図１０】
2008287555.tif 000012
【図１１】
2008287555.tif 000013
【図１２】
2008287555.tif 000014