読む 編集 履歴表示

Fast inverse square root

ウィキペディアから、無料の百科事典

高速逆平方根（こうそくぎゃくへいほうこん、英: fast inverse square root）は、IEEE 754形式の32ビット単精度浮動小数点数 x について、平方根の逆数 1/x の近似値を高速に求めるアルゴリズムである[注釈 1]。浮動小数点数のビット列を整数とみなし、1ビットの右シフトと定数0x5F3759DFからの減算によって初期近似値を得たうえで、ニュートン法によって補正する。1999年に発売されたファーストパーソン・シューティングゲーム『Quake III Arena』のソースコードに含まれていたことで広く知られるようになった[1.2][3]。

3次元コンピュータグラフィックスで多用されるベクトルの正規化を、浮動小数点の除算や平方根の計算を用いずに高速に行えるため、グラフィックスの計算に用いられた[4.1][5.1]。作成者は長らく不明とされ、ジョン・カーマックの作と考えられることも多かったが[3]、2006年の調査で、1980年代後半にアーデントコンピュータのグレッグ・ウォルシュがクリーブ・モラーの影響を受けて作成したものとされた[6][7]。モラー自身は、この技法を1986年ごろのウィリアム・カハンとK・C・ンのコードから学んだと述べている[7]。

このコードが知られるようになると、定数の由来や最適な値が研究され、クリス・ロモント（wikidata）らによって、ニュートン法を適用した後の誤差がわずかに小さい定数0x5F375A86が示されたほか、補正式の係数まで含めて最適化した改良版も提案されている[8.1][2.2]。一方、平方根の逆数の近似値を求める命令を備えたプロセッサが普及したことで、汎用の計算における有用性は低下しているが、そうした命令を持たない組み込み機器などへの応用が期待されている[2.2][9.1]。

背景

3次元コンピュータグラフィックスでは、ベクトルを長さ1に正規化する処理が至る所で必要となる[5.1][10.1]。ベクトル 𝐯=(v1,v2,v3) の正規化は

𝐯̂=𝐯‖𝐯‖=𝐯⋅1v12+v22+v32

で与えられるから、各成分の2乗和の平方根の逆数が求まれば、あとは乗算だけで済む[1.3]。マシュー・ロバートソンは、『Quake III Arena』のような当時のゲームのグラフィックスは、2012年当時より遅いプロセッサで多くの計算を実時間で行う必要があったため、平方根の逆数を高速に求める必要があったと説明している[1.3]。また、浮動小数点数の平方根の逆数そのものは、グラフィックス以外にも画像処理やデジタル信号処理、通信、科学技術計算などで広く使われている[5.2]。

素朴には平方根を計算してから除算すればよいが、除算や平方根は、加算・減算・乗算に比べて回路の設計が複雑で、処理にも時間がかかる演算である[5.2]。高速逆平方根は、整数演算と浮動小数点の乗算・減算だけで 1/x の近似値を得られるため、このような計算の高速化に用いられた[4.1]。マイク・デイは、1990年代のゲームのグラフィックスエンジンで、照明計算に用いるベクトルの正規化を高速に行えたことから、このアルゴリズムが重宝されたとしている[2.2]。平方根の逆数を反復で近似する場合、初期値は多くの場合、メモリを消費するルックアップテーブルから求められるが[9.2]、高速逆平方根は表を使わずに初期値を作る。このコードを使っていたゲイリー・タロリ（wikidata）は、計算量が少ないわりに精度がよく、表も必要としない点を高く評価し、整数演算が速く浮動小数点演算が遅かった当時を回想している[3]。2003年にロモントが Visual C++ .NET で行った計測では、(float)(1.0/sqrt(x)) による素朴な計算よりおよそ4倍高速であった[8.2]。

アルゴリズム

『Quake III Arena』の実装

以下は『Quake III Arena』のソースコード（code/game/q_math.c）に含まれるC言語の関数 Q_rsqrt である。プリプロセッサ指令による部分は省略したが、コメントは原文のままである[11]。

float Q_rsqrt( float number )
{
	long i;
	float x2, y;
	const float threehalfs = 1.5F;

	x2 = number * 0.5F;
	y  = number;
	i  = * ( long * ) &y;                       // evil floating point bit level hacking
	i  = 0x5f3759df - ( i >> 1 );               // what the fuck?
	y  = * ( float * ) &i;
	y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//	y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

	return y;
}

同じファイルの関数 VectorNormalizeFast は、ベクトルの成分の2乗和をこの関数に渡し、得られた値を各成分に掛けることでベクトルを正規化している[11]。

この関数の処理は次のとおりである[8.3][1.4]。

  1. 入力 x の半分の値 x2 を後のために保存する。
  2. x のビット列を、変換を行わずにそのまま整数 i として読み出す[注釈 2]。
  3. i を1ビット右にシフトし（整数として2で割ることに相当する）、その結果を定数0x5F3759DFから引く。
  4. 得られたビット列を再び浮動小数点数として読み出し、初期近似値 y とする。
  5. ニュートン法の反復を1回行って近似値を改善する。2回目の反復はコメントアウトされている。

浮動小数点数の整数としての解釈

正の単精度浮動小数点数 x のビット列を整数として読んだ値 Ix（青）と、二進対数を拡大・平行移動した Llog2x+L(B−σ)（灰）の比較。図の σ≈0.0430 は、対数の近似誤差の最大値を最小にする値である[4.3]。

正の正規化数である単精度浮動小数点数 x は、仮数の小数部分 mx∈[0,1) と指数 ex を用いて

x=(1+mx)2ex

と表される。IEEE 754の単精度形式では、先頭の1ビットが符号、続く8ビットが指数部、残りの23ビットが仮数部の小数部分であり、指数部にはバイアス B=127 を加えた値 Ex=ex+B が、仮数部には整数 Mx=Lmx（L=223）が格納される。したがって、このビット列を符号なし整数として読んだ値 Ix は

Ix=LEx+Mx=L(ex+B+mx)

となる[5.1][4.4]。

一方、x の二進対数は log2x=ex+log2(1+mx) である。0≤mx<1 の範囲では log2(1+mx)≈mx と近似でき、誤差を補正する小さな定数 σ を加えて log2(1+mx)≈mx+σ とおくと

Ix≈L(log2x+B−σ)

が得られる。つまり、浮動小数点数のビット列を整数として読んだ値は、その数の二進対数を拡大・平行移動したものの区分線形な近似になっている[4.4][12][10.2]。ジム・ブリンは1997年の記事で、この性質を「正の数だけを扱うなら、浮動小数点数のビットパターンを整数として解釈したものは、対数関数の区分線形近似を与える」と紹介している[10.2]。ジャン＝ミシェル・ミュレール（英語版）は、これが1962年にJ・N・ミッチェルが乗除算の近似のために用いた二進対数の近似と同じものであると指摘している[14.1]。

初期近似値と定数0x5F3759DF

求める値 y=1/x は log2y=−12log2x を満たす。両辺を上の近似で整数表現に置き換えると

Iy≈32L(B−σ)−12Ix

となる。右辺第2項の 12Ix は1ビットの右シフト i >> 1 で、第1項は定数 R=32L(B−σ) で計算できる。これがコードの i = 0x5f3759df - (i >> 1); に対応する[4.5]。この定数 R は、マジックナンバー、マジック定数などと呼ばれる[8.4][5.1]。

σ=0 とすると R は0x5F400000となり、0x5F3759DFに対応する σ の値はおよそ0.0450466である[4.6]。ロモントは、入力の指数部が偶数の場合に初期近似値の相対誤差を1/8未満に抑えうる可能性がある指数部の値は190（0xBE）しかないことを示し、奇数の場合にも同じ値を用いている[8.5]。定数の上位9ビット（符号ビットの0と、指数部に当たる0xBE）は、指数を −12 倍したうえでバイアスを付け直すための部分であり[15.1]、残りの仮数部に当たるビットは近似の誤差を小さくするように選ばれている[8.6]。

実際には、右シフトによって指数部の最下位ビットが仮数部の最上位ビットへ移り、また減算で仮数部が桁借りを起こすことがあるため、初期近似値は入力の指数が偶数か奇数か、仮数部が大きいか小さいかによって異なる一次式で表される区分線形関数となる。ロモントやロバートソンは、こうした場合分けを詳しく解析している[8.7][1.5]。ロバートソンは、実際に起こりうるのは「指数が偶数で仮数部が小さい」「指数が偶数で仮数部が大きい」「指数が奇数」の3通りだけで、ロモントやデイヴィッド・エバリー（wikidata）が扱った4通り目（指数が奇数で桁借りが起きない場合）は定数の仮数部が小さいために起こらないと指摘した。そのうえで、この3通りの式がすべての正の単精度浮動小数点数について実際の計算結果と一致することを、全数検査で確かめている[1.6]。

ニュートン法による補正

初期近似値は、関数

f(y)=1y2−x

の正の零点をニュートン法で求めることで改善される。反復式は

yn+1=yn−f(yn)f′(yn)=yn(32−x2yn2)

となり、除算を含まない[8.4][1.7]。コードの y = y * ( threehalfs - ( x2 * y * y ) ); がこれに当たる。この反復は、初期値が 0<y0<3/x の範囲にあれば 1/x に収束し、初期値が真の値に近ければ収束は2次（反復ごとに正しい桁数がおよそ倍になる）である[14.2]。

初期値の相対誤差を δ（yn=(1+δ)/x）とすると、反復後の相対誤差は

δ′=−32δ2−12δ3

となり、δ の符号によらず0以下になる[5.3][9.3]。

計算例

x = 16の場合、ビット列は整数として0x41800000であり、1ビット右シフトすると0x20C00000となる。これを0x5F3759DFから引くと0x3E7759DFとなり、浮動小数点数として読むと約0.24155であり[8.8]、真の値0.25に対する相対誤差は約−3.4%である。ここからニュートン法の反復を1回行うと約0.24958となり、相対誤差は約−0.17%に縮まる。各ビット列の内訳は次のとおりである。

x = 16 の場合のビット列
値符号指数部仮数部16進表記
入力 x010000011000000000000000000000000x41800000
右シフト後001000001100000000000000000000000x20C00000
定数010111110011011101011001110111110x5F3759DF
減算結果001111100111011101011001110111110x3E7759DF

右シフトによって入力の指数部の最下位ビット（1）が仮数部の最上位へ移り、減算では仮数部で桁借りが起きて指数部が1減る。結果の指数部は124（真の指数 −3）、仮数は約1.9324となり、1.9324×2−3≈0.24155 が初期近似値となる[8.8][注釈 3]。

精度

入力 x（0.01〜10000）に対する 1/x（緑）と、ニュートン法を1回適用した高速逆平方根の結果との差（赤）。両軸とも対数目盛である。

ロモントは、すべての単精度浮動小数点数について相対誤差の最大値を調べ、定数0x5F3759DFの場合、初期近似値では約3.44%、ニュートン法1回の後では約0.175%（0.00175228）、2回の後では約0.000466%であることを示した。ロモントは、2回目の反復を加えても所要時間の増加はごくわずかだが、精度は大きく向上したと述べている[8.9]。一方、ロバートソンの2012年の計測では、2回目の反復を有効にすると速度上の利点がおよそ半分になった[1.8]。ロバートソンも、正の正規化数全体での最大相対誤差を0.0017522874と求めている[1.9]。ただし、0に非常に近い非正規化数では精度が低下する[1.10][13.1]。

1/x の曲線の形状のため、この反復を1回行うと、近似値は初期値の大小によらず真の値以下になる。そのため Q_rsqrt の結果は常にわずかに小さめとなる[1.11][5.3]。ロバートソンによれば、このため初期値は真の値よりわずかに小さい方が、大きい場合よりも反復後の誤差が小さくなる[1.11]。

また、正規化数の範囲では、入力を4倍すると結果はニュートン法の反復回数によらずちょうど半分になる。真の値 1/x もちょうど半分になるので、相対誤差は x と 4x とで等しい。したがって、1以上4未満の入力について誤差を調べれば、すべての正規化数について調べたことになる[13.3][12]。

歴史

起源

ウィリアム・カハンとK・C・ンは1986年5月、IEEE 754の倍精度浮動小数点数の平方根をソフトウェアで正しく丸めて計算する方法を述べた論文の草稿を書いた。この草稿は未公刊であるが[2.2]、サン・マイクロシステムズの数学ライブラリfdlibmに含まれる平方根関数のソースコードに、コメントとして写しが収められている。草稿には、除算を避けるための方法として、倍精度数の上位32ビットを整数とみなして1ビット右シフトし、定数0x5FE80000から引いたうえで表による補正を加え、1/x の約7.8ビットの近似値を得てから反復で精度を高める手順が記されている。同じ草稿には、上位32ビットを1ビット右シフトして定数0x1FF80000を加え、表で補正して平方根そのものの約8ビットの近似値を得る、ニュートン法向けの手順も並べて記されている[16]。

MATLABの作者でマスワークスの創業者の一人でもある[6]クリーブ・モラーは、2012年に自身のブログへのコメントで、この技法はバークレーでカハンとンが1986年ごろに書いたコードから学んだと述べている。同じコメントでモラーは、この技法がグラフィックス界隈で一種のカルト的な地位を得ていたとは知らなかったとも述べている[7]。1980年代後半、モラーはアーデントコンピュータで、同社の共同創業者であるグレッグ・ウォルシュと共に働いていた[6][7]。ウォルシュは1970年代から計算機業界に携わり、スタンフォード大学在学中にゼロックスのパロアルト研究所で初期のWYSIWYGワープロの開発にも関わった人物で、のちにはタロリと共にAccel Graphicsでも働いている。アーデントは独自のベクトルプロセッサやIntel i860を用いた並列グラフィックスミニコンピュータを開発していたが、後に出資元のクボタの意向で競合のステラー（Stellar）と合併し、スターデント（Stardent）となった[6]。同社のグラフィックスミニコンピュータ「Titan」は公称の性能を出すのに苦労しており、ウォルシュは、ベクトル演算ハードウェアを利用できないソフトウェアを高速化する手段として 1/x の高速な計算法を作成した。ウォルシュは、ニュートン法による近似を研究していたモラーとの協働からその着想を得たといい、同時期に立方根の逆数を求める関数も同様の（より複雑な）手法で書いている[6]。

流布

シリコングラフィックスや3dfxに在籍したタロリは、アーデントに出資していたクボタのコンサルティングも行っており、3dfxではこのコードを自身の用途に使い、調整していた[6]。タロリは、後述するBeyond3Dのリス・ソムフェルトの問い合わせに対する2005年の電子メールで、このコードを10年以上前に見かけて自ら導出し直し、定数の値をいろいろと変えて試したことを覚えていると述べ、自身の功績ではないが、多用したことでこのコードの普及と長寿に貢献しただろうとしている[3]。ヘンリー・S・ウォーレン（wikidata）は2013年の『Hacker's Delight』第2版で、最初期の利用者の一人としてタロリを挙げ、タロリはこのコードが1995年かそれ以前にさかのぼると述べているとしている[13.4]。ソムフェルトは、少なくともブライアン・フック（wikidata）を通じて3dfxからid Softwareへつながる経路をたどることができると述べている[6]。

『Quake III Arena』は1999年12月に発売された[17]。ソースコードが公開される前から、この関数はインターネット上で出回っており、2002年1月9日にはニュースグループcomp.graphics.algorithmsに、Quake 3のものとされるコードが投稿された[15.2]。同時期のゲームには、似た考え方を用いた別の実装もある。アクティビジョンの『インターステート'76（英語版）』（1997年）を逆アセンブルして解析したシェーン・ピーラーによれば、同作は倍精度の平方根の逆数を求めるのに、指数部を整数の減算と右シフトで、仮数部の上位8ビットを256項目の表から作って初期値とし、ニュートン法を1回適用してから1.00001を掛けていた。精度は最悪でも有効数字4桁程度だったという[18]。

ロモントもゲーム開発者向けフォーラムgamedev.netでこのコードを知り、2003年2月に解析を発表した[8.10]。ウォーレンは、2000年代初頭にこの関数がプログラマの間で話題になったと記している[13.1]。

作者探し

コードが『Quake III Arena』に含まれていたことから、作者は同作の開発者であるカーマックだと考えられることが多く、x86アセンブリによる最適化で知られるマイケル・アブラッシュ（英語版）の名も挙がった。ソムフェルトが2004年4月にカーマックへ問い合わせたところ、カーマックは自身の作ではなく、アブラッシュでもないだろうと答え、テリエ・マティセンの名を挙げた[3]。

2005年8月19日、id Softwareは『Quake III Arena』のソースコードをGNU General Public Licenseの下で公開し[6][17][19]、関数の作者についての議論が再燃した[3][1.3]。同月にソムフェルトが問い合わせたマティセンも、自身の作であることを否定した。マティセンは2000年ごろ、水分子の計算化学シミュレーションを高速化するため、48ビットの精度を保つ逆平方根を、表引きとニュートン法による別の方法で書いたことがあったという。またマティセンは、このコードの書きぶりはMITのHAKMEM（英語版）文書を思わせるとも述べている[3]。ソムフェルトはタロリに行き着いたところでいったん調査を締めくくり、個人サイトで発表していたその経緯を、2006年11月29日にBeyond3Dへ記事として再掲した[3]。この記事がスラッシュドットで取り上げられたことをきっかけにウォルシュ本人から連絡があり、同年12月19日の続報で、モラーを着想の源、ウォルシュを作者とする結論が示された[6]。2006年以降の文献での扱いは一様ではない。チャールズ・マケニリーは2007年に、この続報によってウォルシュが最も有力な起源とされたと紹介し[4.1]、デイも2023年に、ソムフェルトの調査の結論としてウォルシュとモラーによる開発を紹介している[2.2]。一方、ロバートソンは2012年の卒業論文で原作者を依然として不明とし[1.12]、ウォーレンも2013年の著書で、作者は歴史の中に埋もれてしまったようだとしている[13.4]。

定数の解析と改良

このコードが知られるようになってから、定数0x5F3759DFがどのように選ばれたのか、より良い定数があるのかが研究されてきた。それに先立つ例として、ブリンの記事がある。

  • ブリンは1997年の『IEEE Computer Graphics and Applications』誌の記事で、ベクトルの正規化に役立つものとして同種の逆平方根の近似を示した。その定数は1.0のビット表現と、それを1ビット右シフトした値との和、すなわち0x5F400000であり[10.1][14.2]、初期近似値の最大相対誤差は約8.9%と粗い[14.2]。ブリンはニュートン法を1回組み込んだ版も示したうえで、補正式の係数を1.5と0.5から1.47と0.47に変えると精度が2倍になると述べている[10.3]。
  • エバリーは2002年1月、この定数を仮数部に対する一次関数 0.966215−M/4 による近似として説明した[15.3]。ロモントは、この説明は誤りだが示唆に富むとし、実際の初期近似値は区分線形であると指摘した[8.11]。エバリーは、ロモントの論文を受けて2010年に自身の解説を改訂し、扱う場合分けを増やしている[1.13]。
  • ロモントは2003年、初期近似値の最大相対誤差を最小にする定数として0x5F37642Fを導いた。しかしこの定数は、ニュートン法を1回適用した後では元の定数より誤差が大きかった。そこで数値的な探索を行い、反復1回後と2回後のいずれでも元の定数より誤差の小さい0x5F375A86を見いだした。ロモントは、元の定数が導出によるものか試行錯誤によるものかを問うている[8.1]。
  • マケニリーは2007年、前述の対数近似に基づいて定数を導出して誤差を解析し、反復1回後の最大相対誤差を最小にする定数として0x5F375A86を得た[4.7]。マケニリーの試したどの方法でも0x5F3759DFそのものは導けなかったが、σ = 0 と σ = τ（log2(1+m)−m の最大値、約0.086）の間を重みを付けた二分法で探索する過程でこの値が現れたことから、マケニリーは、元の定数は許容誤差に達するまで二分法で求められたのかもしれないと推測している[4.8]。マケニリーは、減算を先に行ってから1ビット右シフトする変形（定数は0xBE6EB50C）も代案として示している[4.2]。
  • ロバートソンは2012年のニューブランズウィック大学（英語版）の卒業論文で、ニュートン法適用後の誤差を直接最小化することで、0x5F375A86を解析的に導いた[1.14]。
  • ウォーレンは『Hacker's Delight』第2版で、ニュートン法を省けば整数命令2つだけで最大相対誤差約3.5%の近似が得られること（定数は0x5F37642F）や、定数を0x5F375A82とした版を示し、補正式の係数1.5を1.5008908に変えると反復1回後の相対誤差を±0.000892の範囲に抑えられること、反復を2回行う場合の最適な定数は0x5F37599Eであることなどを述べている。また、0.5f*x の乗算を指数部から1を引く整数演算に置き換える案も示しているが、その場合は約 2.34×10−38 未満の正規化数で結果が不正確になり、非正規化数ではNaNになるとしている[13.1]。
  • レオニード・モロズらは、このアルゴリズムの各段階を初めて数学的に厳密に記述したとし、相対誤差を最小にする定数として、初期近似値に対する0x5F37642Fと、反復1回後・2回後に対する0x5F375A86を解析的に導いた。0x5F375A86の場合の最大相対誤差は、反復1回後で約 1.75118×10−3、2回後で約 4.60×10−6 である[5.4]。モロズらはさらに、相対誤差ではなく絶対誤差を最小にする定数として、初期近似値に対する0x5F3863F7、反復1回後に対する0x5F37E75A、反復2回後に対する0x5F37ADD5も導いており、最適な定数は誤差の尺度と反復回数によって異なることを示した[5.5]。

倍精度浮動小数点数についても同様の定数が求められている。ロモントは初期近似値に最適な定数として0x5FE6EC85E7DE30DAを示した[8.12]。ニュートン法1回後の誤差を最小にする定数として、マケニリーは0x5FE6EB50C7B537AAを[4.9]、ロバートソンは0x5FE6EB50C7B537A9を求めている[1.15]。ロバートソンは、この方法の精度は浮動小数点数の幅に依存しないため、倍精度や四倍精度に適用しても反復1回後の最大相対誤差は単精度とほぼ同じ（約0.175%）になると指摘し、四倍精度に対する定数として0x5FFE6EB50C7B537A9CD9F02E504FCFBFも示している[1.16]。

主な定数（単精度）と最大相対誤差
定数由来最大相対誤差
初期近似値反復1回後
0x5F400000σ=0 の場合。ブリンの記事のコードと同じ約8.9%[14.2]約1.2%[注釈 4]
0x5F3759DF『Quake III Arena』約3.44%[8.9][14.2]約0.1752%[8.9]
0x5F37642F初期近似値の誤差を最小化（ロモント、モロズら）約3.42%[8.9][14.2]約0.1776%[8.9]
0x5F375A86反復後の誤差を最小化（ロモント、マケニリー、ロバートソン、モロズら）約3.44%[8.9]約0.1751%[8.9][5.6]

定数だけでなくニュートン法の補正式の係数も変えることで、計算量を増やさずに精度を高める試みもある。ヤン・カドレツは2010年、定数と補正式の2つの係数を同時に最適化し、

y←0.703952253y(2.38924456−xy2)（定数は0x5F1FFFF9）

によって、反復1回後の最大相対誤差を約 6.50×10−4 と、元の約2.7分の1に減らした[12][2.2]。デイによれば、3つの定数をすべて調整して誤差をおよそ2.7分の1にする試みは、2008年にピザーが先に行っていた[2.2]。ツェザリ・ヴァルチク、モロズ、ヤン・チェシリンスキ（wikidata）は2021年、補正式を改めた変種を解析的に導き、反復2回の場合、定数を0x5F375A86とした原型と比べて、計算量が同じ変種で約6.8倍、乗算が1回多い変種で約12倍（いずれも単精度）精度が高くなることを示した[9.4]。この論文は、非正規化数への拡張や、倍精度およびそれ以上の精度への適用も扱っている[9.5]。

同じ考え方は平方根の逆数以外のべき乗にも適用できる。ブリンの1997年の記事は、逆数・平方根・平方根の逆数の近似に加え、対数関数や2のべき乗の近似も同じ方法で示している[10.4]。ミュレールによれば、平方根を整数として1ビット右シフトして 127×222 を加えるブリンの近似の最大相対誤差は0.0607で、加える定数をやや小さい532369100にすると0.03476まで下がる[14.2]。マケニリーは、平方根や立方根、一般のべき乗 xa に対応する初期近似の式と反復式を示したうえで、多くの場合は整数演算の利点を生かせないと述べている[4.2]。デイは2023年に、この手法を任意の有理数乗に一般化し、定数や補正式の係数の最適値を自動的に構成する手順を示している[2.2]。

ミュレールは2020年の概説で、浮動小数点数のビット列を操作するこの種の手法について、ソフトウェアで数パーセント程度の粗い近似で足りる場合には勝る方法はなかなかないとする一方、それだけで得られる精度は5〜6ビット程度にとどまり、適用できる関数も2を底とする指数関数・対数関数や小さなべき乗などに限られると評している[14.3]。

性能と現在の利用

ゲーム開発者のエラン・ラスキンが2009年に2.66 GHzのCore 2で計測したところ、1要素あたりの所要時間と平均相対誤差は、高速逆平方根が3.54ナノ秒・0.0990%、SSEのRSQRTSS命令が0.85ナノ秒・0.0094%、RSQRTSSにニュートン法を1回加えたものが2.74ナノ秒で誤差はほぼ0であった。一方、平方根そのものを求める場合には、高速逆平方根の結果に x を掛ける方法（4.33ナノ秒）が、x87のFSQRT命令（24ナノ秒）やSSEのSQRTSS命令（11.9ナノ秒）より速かった[20]。

ロバートソンが2012年に当時のパソコン用のプロセッサ数種で計測したところ、GCCで最適化してコンパイルした 1.0f/sqrtf() と比べて、高速逆平方根はおよそ3倍から4倍高速であった[1.8]。一方モロズらは、この手法はインテルのプロセッサではよく機能するが、たとえばPowerPCでは浮動小数点レジスタと整数レジスタとの間の受け渡しに伴う遅延があり、ハードウェアが備える代替手段の方が効果的な場合があると指摘している[5.1]。

マケニリーは2007年の時点で、GPUの多くやプロセッサの拡張命令セットが、同様の、より複雑な方法をハードウェアで実装しているとしている[4.1]。x86のSSEには単精度の平方根の逆数の近似値を求める命令RSQRTSSがあり、その相対誤差は 1.5×2−12 以下と規定されている[21]。AMDの3DNow!にも、単精度の平方根の逆数を15ビットの精度で推定する命令PFRSQRTがある[22]。デイは、多くのハードウェアが平方根の逆数の計算を高速化する命令を備えるようになったことで、元の技法の有用性は低下したとしている[2.2]。モロズらは、高速逆平方根にはソフトウェアだけでなくハードウェアでの実装例もあることを挙げている[5.1]。ヴァルチクらは、高速逆平方根の手法はなお多くのソフトウェアやハードウェアの実装で有用であるとする一方[9.2]、汎用の計算ではSSEやAVXの命令の方が高速かつ高精度であることを認めている。そのうえで、自分たちの改良版は、こうした命令を持たないマイクロコントローラやFPGAなどの組み込みシステムや、ハードウェアに倍精度の命令がない場合などに応用できると期待している[9.6]。

また、2007年にFPGA上に実装された組み込み向けの頂点シェーダには、ニュートン法の初期値を求めるための命令rsqとして 0x5F3759DF − (a >> 1) が設けられ、その遅延は2サイクルであった（除算命令divは27サイクル）[23.1]。また、Z80などの8ビットCPU向けの開発環境z88dkの浮動小数点ライブラリ（math32）は、2026年10月時点の説明で、平方根（逆平方根を経由する）の初期値を「Quakeのマジックナンバー法」（定数はロモントの0x5F375A86）で求めたうえでニュートン法を3回適用するとしており、ゲームなどでは反復を減らして精度と速度を引き換えにできるとしている[24]。

脚注

注釈

  1. ↑ ここでいう「逆」（inverse）は逆関数ではなく乗法についての逆元、すなわち逆数を意味する。マシュー・ロバートソンは、平方根の逆関数は2乗であるから「逆数平方根」（reciprocal square root）と呼ぶ方が正確だと指摘しており[1.1]、マイク・デイも「逆」の語が曖昧であるとして「高速逆数平方根」（Fast Reciprocal Square Root）の呼称を用いている[2.1]。
  2. ↑ Q_rsqrt は、ポインタの型変換によって float 型のオブジェクトを long 型として読み出している[11]。後の改良版や解説には、共用体を用いて同じ読み替えを行うものもある[12][4.2][13.1]。ヘンリー・S・ウォーレンは、共用体を用いた自身のコードについて、公式に認められたC言語の書き方ではないが、ほとんどのコンパイラで動作すると注記している[13.2]。
  3. ↑ ロモントが示しているのは、初期近似値のビット列0x3E7759DFとその値（0.241553）までである[8.8]。入力とシフト後のビット列、相対誤差、反復後の値は、定義と上記の反復式に従って単精度で計算したものである。
  4. ↑ 補正式の係数を通常どおり1.5と0.5とした場合[1.17]。ブリンが示した係数1.47と0.47を用いると約0.6%となる[1.17]。

出典

  1. Robertson, Matthew (2012年4月24日). A Brief History of InvSqrt [InvSqrt小史] (PDF) (B.Sc. thesis) (英語). University of New Brunswick. 2026年10月14日時点のオリジナルよりアーカイブ (PDF). 2026年11月4日閲覧.
    1. ↑ p. 3
    2. ↑ pp. 1–3
    3. 1 2 3 p. 2
    4. ↑ pp. 9–12
    5. ↑ pp. 11–19
    6. ↑ pp. 13, 17–20, 37
    7. ↑ p. 12
    8. 1 2 p. 21
    9. ↑ p. 5
    10. ↑ p. 6
    11. 1 2 p. 28
    12. ↑ pp. iii, 3
    13. ↑ p. 37
    14. ↑ pp. 28–33
    15. ↑ p. 33
    16. ↑ pp. 32–33
    17. 1 2 p. 35
  2. Day, Mike (2023年7月28日). “Generalising the Fast Reciprocal Square Root Algorithm” [高速逆数平方根アルゴリズムの一般化] (英語). arXiv:2307.15600 [math.NA].
    1. ↑ p. 2
    2. 1 2 3 4 5 6 7 8 9 p. 1
  3. 1 2 3 4 5 6 7 8 Sommefeldt, Rys (2006年11月29日). “Origin of Quake3's Fast InvSqrt()” [Quake3の高速InvSqrt()の起源]. Beyond3D (英語). 2026年5月30日時点のオリジナルよりアーカイブ. 2026年11月4日閲覧.
  4. McEniry, Charles (2007). The Mathematics Behind the Fast Inverse Square Root Function Code [高速逆平方根関数のコードの背後にある数学] (PDF) (Report) (英語). 2009年12月29日時点のオリジナル (PDF)よりアーカイブ. 2026年11月5日閲覧.
    1. 1 2 3 4 p. 1
    2. 1 2 3 p. 16
    3. ↑ pp. 6–7
    4. 1 2 pp. 2–3
    5. ↑ p. 3
    6. ↑ pp. 3, 12
    7. ↑ pp. 2–15
    8. ↑ pp. 6, 16
    9. ↑ p. 15
  5. Moroz, Leonid V.; Walczyk, Cezary J.; Hrynchyshyn, Andriy; Holimath, Vijay; Cieśliński, Jan L. (2018年1月). “Fast calculation of inverse square root with the use of magic constant – analytical approach” [マジック定数を用いた逆平方根の高速計算：解析的アプローチ]. Applied Mathematics and Computation (英語). 316: 245–255. arXiv:1603.04483. doi:10.1016/j.amc.2017.08.025.
    1. 1 2 3 4 5 6 p. 246
    2. 1 2 p. 245
    3. 1 2 p. 250
    4. ↑ pp. 245, 250–251
    5. ↑ pp. 245, 253–254
    6. ↑ p. 251
  6. 1 2 3 4 5 6 7 8 9 Sommefeldt, Rys (2006年12月19日). “Origin of Quake3's Fast InvSqrt() - Part Two” [Quake3の高速InvSqrt()の起源 第2部]. Beyond3D (英語). 2026年5月10日時点のオリジナルよりアーカイブ. 2026年11月4日閲覧.
  7. 1 2 3 4 Moler, Cleve (2012年6月19日). “Symplectic Spacewar” [シンプレクティック・スペースウォー]. Cleve's Corner (英語). MathWorks. 2014年10月15日時点のオリジナルよりアーカイブ. 2026年11月4日閲覧. 2012年6月27日付の著者によるコメント。
  8. Lomont, Chris (2003年2月). “Fast Inverse Square Root” [高速逆平方根] (PDF) (英語). 2026年10月19日時点のオリジナルよりアーカイブ (PDF). 2026年11月4日閲覧.
    1. 1 2 pp. 10–11
    2. ↑ p. 1
    3. ↑ pp. 1, 3
    4. 1 2 p. 3
    5. ↑ pp. 5, 7
    6. ↑ pp. 3–5
    7. ↑ pp. 3–8
    8. 1 2 3 p. 8
    9. 1 2 3 4 5 6 7 p. 10
    10. ↑ pp. 1, 12
    11. ↑ pp. 1–3
    12. ↑ p. 12
  9. Walczyk, Cezary J.; Moroz, Leonid V.; Cieśliński, Jan L. (2021年1月9日). “Improving the Accuracy of the Fast Inverse Square Root by Modifying Newton–Raphson Corrections” [ニュートン・ラフソン補正の修正による高速逆平方根の精度向上]. Entropy (英語). 23 (1). Article 86. doi:10.3390/e23010086.
    1. ↑ p. 15
    2. 1 2 p. 2
    3. ↑ p. 3
    4. ↑ pp. 3, 13–14
    5. ↑ p. 1
    6. ↑ pp. 1–2, 15
  10. Blinn, James F. (1997). “Floating-Point Tricks” [浮動小数点のトリック]. IEEE Computer Graphics and Applications (英語). 17 (4): 80–84. doi:10.1109/38.595279.
    1. 1 2 p. 83
    2. 1 2 p. 81
    3. ↑ pp. 83–84
    4. ↑ pp. 81–84
  11. 1 2 3 “Quake-III-Arena/code/game/q_math.c”. GitHub (英語). id Software. 2026年9月6日時点のオリジナルよりアーカイブ. 2026年11月4日閲覧.
  12. 1 2 3 4 Kadlec, Jan (2010). “Improving the fast inverse square root” [高速逆平方根の改良]. Řrřlog (英語). 2025年12月19日時点のオリジナルよりアーカイブ. 2026年11月4日閲覧.
  13. Warren, Henry S. Jr. (2013). Hacker's Delight [ハッカーのたのしみ] (英語) (2nd ed.). Addison-Wesley. ISBN 978-0-321-84268-8. 第1版の邦訳は『ハッカーのたのしみ』（エスアイビー・アクセス、2004年）。
    1. 1 2 3 4 p. 383
    2. ↑ p. 384
    3. ↑ pp. 384–385
    4. 1 2 p. 478
  14. Muller, Jean-Michel (2020年12月). “Elementary Functions and Approximate Computing” [初等関数と近似計算]. Proceedings of the IEEE (著者最終稿) (英語). 108 (12): 2136–2149. doi:10.1109/JPROC.2020.2991885. 2025年2月10日時点のオリジナルよりアーカイブ. 2026年11月4日閲覧. 節番号は著者最終稿による。
    1. ↑ 第II節、第VI節、参考文献[15]
    2. 1 2 3 4 5 6 7 第VI節
    3. ↑ 第VI節、第VII節
  15. Eberly, David (2002年1月26日). “Fast Inverse Square Root” [高速逆平方根] (PDF) (英語). Magic Software. 2003年4月26日時点のオリジナルよりアーカイブ. 2026年11月4日閲覧.
    1. ↑ pp. 1–2
    2. ↑ p. 1
    3. ↑ p. 2
  16. ↑ “e_sqrt.c”. fdlibm (英語). Netlib. 2026年6月9日時点のオリジナルよりアーカイブ. 2026年11月4日閲覧. ファイル末尾のコメント「Other methods (use floating-point arithmetic)」。
  17. 1 2 Howarth, Robert (2005年8月19日). “Quake 3 Arena Source Code” [Quake 3 Arenaのソースコード]. IGN (英語). 2026年3月6日時点のオリジナルよりアーカイブ. 2026年11月4日閲覧.
  18. ↑ Peelar, Shane (2021年6月1日). “Fast reciprocal square root... in 1997?!” [高速逆平方根が1997年に？]. Shane Peelar's Blog (英語). 2026年1月4日時点のオリジナルよりアーカイブ. 2026年11月5日閲覧.
  19. ↑ 朝倉哲也「往年の名作「Quake III Arena」のソースコードが公開に」『4Gamer.net』2005年8月24日。2014年10月13日時点のオリジナルよりアーカイブ。2026年11月4日閲覧。
  20. ↑ Ruskin, Elan (2009年10月16日). “Timing square root” [平方根の計算時間の計測]. Some Assembly Required (英語). 2021年4月15日時点のオリジナルよりアーカイブ. 2026年11月5日閲覧.
  21. ↑ “Intel 64 and IA-32 Architectures Software Developer's Manual, Volume 2 (2A, 2B, 2C, & 2D): Instruction Set Reference, A-Z” [インテル64およびIA-32アーキテクチャー ソフトウェア開発者マニュアル 第2巻：命令セット・リファレンス] (PDF) (英語). Intel. 2026年9月. Order Number 325383-093US. 2026年10月13日時点のオリジナルよりアーカイブ (PDF). 2026年11月4日閲覧. Vol. 2B, p. 4-596「RSQRTSS—Compute Reciprocal of Square Root of Scalar Single Precision Floating-Point Value」。
  22. ↑ “3DNow! Technology Manual” [3DNow!テクノロジー・マニュアル] (PDF) (英語). Advanced Micro Devices. 2000. Publication 21928G/0 (March 2000). 2026年3月12日時点のオリジナルよりアーカイブ (PDF). 2026年11月5日閲覧. p. 47「PFRSQRT」の項。
  23. Middendorf, Lars; Mühlbauer, Felix; Umlauf, Georg; Bobda, Christophe (2007). “Embedded Vertex Shader in FPGA” [FPGAによる組み込み頂点シェーダ]. Embedded System Design: Topics, Techniques and Trends. IFIP Advances in Information and Communication Technology (英語). Springer. pp. 155–164. doi:10.1007/978-0-387-72258-0_14. ISBN 978-0-387-72257-3.
    1. ↑ p. 161
  24. ↑ “z88dk/libsrc/math/float/math32/readme.md”. GitHub (英語). z88dk. 2026年10月4日. 2026年11月5日時点のオリジナルよりアーカイブ. 2026年11月5日閲覧.

関連項目

関連記事

浮動小数点数