Algorithms for calculating square roots
平方根アルゴリズムは、正の 実数 の 非負 平方根を計算します。 自然数 の平方根は 、 完全平方数 を除き、すべて 無理数 であるため、 平方根は通常、有限の精度でしか計算できません。これらの アルゴリズムは 通常、精度が徐々に上がっていく一連の 近似値 を構築します。 S {\displaystyle {\sqrt {S}}} S {\displaystyle S}
ほとんどの平方根計算法は反復的です。適切な初期推定値 を選択した後 、何らかの終了基準が満たされるまで 反復的な改良が行われます。改良法の一つに、 ニュートン法 の特殊なケースであるヘロン法があります 。除算が乗算よりもはるかにコストがかかる場合は、代わりに逆平方根を計算する方が望ましい場合があります。 S {\displaystyle {\sqrt {S}}}
平方根を桁ごとに計算する方法やテイラー級数を用いる方法もあります。平方根の有理近似値は、連分数展開を用いて計算できます。
採用される方法は、必要な精度、利用可能なツール、そして計算能力によって異なります。これらの方法は、暗算に適したもの、通常は少なくとも紙と鉛筆を必要とするもの、そしてデジタル電子計算機やその他の計算装置で実行されるプログラムとして実装されたものに大別されます。アルゴリズムは、収束性(指定された精度を達成するために必要な反復回数)、個々の演算(例えば除算)または反復の計算複雑性、そして誤差伝播(最終結果の精度)を考慮する場合があります。
紙と鉛筆による合成除算や級数展開など、いくつかの方法では初期値は必要ありません。一部の応用では、 整数の平方根 、つまり平方根を最も近い整数に丸めたり切り捨てたりした値が必要になります(この場合は、手順を修正する必要があります)。
歴史 平方根(特に 2の平方根 )を求める手順は、少なくとも紀元前17世紀の古代バビロニア時代から知られていました。 バビロニアの数学者たちは、1の後に60進法で3 桁の 数字を使って2の平方根を計算していました が、その正確な方法は分かっていません。彼らは斜辺を (例えば、高さが ロッド、幅がロッド の門の対角線を)近似する方法を知っていました。そして、同様の手法を用いて の近似値を求めていた可能性があります。 a 2 + b 2 ≈ a + b 2 2 a {\displaystyle {\sqrt {a^{2}+b^{2}}}\approx a+{\frac {b^{2}}{2a}}} 41 60 + 15 3600 {\displaystyle {\frac {41}{60}}+{\frac {15}{3600}}} 40 60 {\displaystyle {\frac {40}{60}}} 10 60 {\displaystyle {\frac {10}{60}}} 2 . {\displaystyle {\sqrt {2}}.}
1世紀エジプトのヘロン法は、平方根を計算するための最初の確認可能なアルゴリズムでした。
現代の分析手法は、ルネサンス初期に西ヨーロッパに アラビア数字 が導入されてから開発され始めました。 [ 要出典 ]
現在、ほぼすべてのコンピューティング デバイスには、プログラミング言語構造、コンパイラの組み込み関数またはライブラリ関数、あるいは説明した手順のいずれかに基づくハードウェア演算子として、高速かつ正確な平方根関数が備わっています。
初期見積もり 多くの反復平方根アルゴリズムでは、初期の シード値 が 必要です。シードは 0 でない正の数値で、1 から ( 平方根を求める数)までの範囲でなければなりません。平方根はその範囲内になければならないからです。シードが平方根から大きく離れている場合、アルゴリズムはより多くの反復を必要とします。 (または ) で初期化すると、 平方根の大きさのオーダーを取得するだけで、おおよそ 回の反復が無駄になります。したがって、精度は限られるかもしれませんが計算しやすい大まかな推定値を持つことは有用です。一般に、初期推定値が優れているほど、収束は速くなります。ニュートン法では、平方根よりいくらか大きいシードは、平方根よりいくらか小さいシードよりもわずかに速く収束します。 S {\displaystyle S} x 0 = 1 {\displaystyle x_{0}=1} S {\displaystyle S} 1 2 | log 2 S | {\displaystyle {\tfrac {1}{2}}\vert \log _{2}S\vert }
一般に、推定は根を含むことが分かっている任意の区間( など )に従います。推定は、 区間 に対する関数近似の特定の値です。より良い推定値を得るには、区間のより厳しい境界を得るか、 に対するより良い関数近似を見つける必要があります 。後者は通常、近似においてより高次の多項式を使用することを意味しますが、すべての近似が多項式であるとは限りません。一般的な推定方法には、スカラー、線形、双曲線、対数などがあります。暗算または紙と鉛筆による推定には、通常、10進法が使用されます。コンピュータによる推定には、2進法がより適しています。推定では、 数値が科学的記数法で表されるため、
指数と 仮数は通常別々に扱われます。 [ x 0 , S / x 0 ] {\displaystyle [x_{0},S/x_{0}]} f ( x ) = x {\displaystyle f(x)={\sqrt {x}}} f ( x ) {\displaystyle f(x)}
小数点推定 通常、数値は 科学的記数法 で と 表され 、 n は整数で、可能な平方根の範囲は です 。 S {\displaystyle S} a × 10 2 n {\displaystyle a\times 10^{2n}} 1 ≤ a < 100 {\displaystyle 1\leq a<100} a × 10 n {\displaystyle {\sqrt {a}}\times 10^{n}} 1 ≤ a < 10 {\displaystyle 1\leq {\sqrt {a}}<10}
スカラー推定 スカラー法では、範囲を区間に分割し、各区間の推定値は単一のスカラー値で表されます。範囲を単一の区間とみなす場合、算術平均(5.5)または幾何平均( )の時間は 妥当な推定値となります。これらの絶対誤差と相対誤差は異なります。一般的に、単一のスカラー値は非常に不正確です。より正確な推定値を得るには、範囲を2つ以上の区間に分割しますが、スカラー推定値は本質的に精度が低くなります。 10 ≈ 3.16 {\displaystyle {\sqrt {10}}\approx 3.16} 10 n {\displaystyle 10^{n}}
2つの区間を幾何学的に分割した場合、平方根は 次のように推定できる [注1] S = a × 10 n {\displaystyle {\sqrt {S}}={\sqrt {a}}\times 10^{n}} S ≈ { 2 ⋅ 10 n if a < 10 , 6 ⋅ 10 n if a ≥ 10. {\displaystyle {\sqrt {S}}\approx {\begin{cases}2\cdot 10^{n}&{\text{if }}a<10,\\6\cdot 10^{n}&{\text{if }}a\geq 10.\end{cases}}}
この推定値は、a = 100 で最大絶対誤差が になり 、a = 1 で最大相対誤差が 100% になります。 4 ⋅ 10 n {\displaystyle 4\cdot 10^{n}}
たとえば、 を として 因数 分解すると 、推定値は となり 、絶対誤差は 246、相対誤差はほぼ 70% になります。 S = 125348 {\displaystyle S=125348} 12.5348 × 10 4 {\displaystyle 12.5348\times 10^{4}} S ≈ 6 ⋅ 10 2 = 600 {\displaystyle {\sqrt {S}}\approx 6\cdot 10^{2}=600} 125348 = 354.0 {\displaystyle {\sqrt {125348}}=354.0}
線形推定 より良い推定値、そして標準的な方法は、小さな弧上の 関数の線形近似です。上記のように、底のべき乗を数から因数分解し 、区間を まで縮約した場合 、弧を張る 割線 、または弧上のどこかにある接線を近似値として使用できますが、弧と交差する最小二乗回帰直線の方がより正確です。 y = x 2 {\displaystyle y=x^{2}} S {\displaystyle S} [ 1 , 100 ] {\displaystyle [1,100]}
最小二乗回帰直線は、推定値と関数の値の平均差を最小化します。その式は です 。並べ替え、 。計算を容易にするために係数を四捨五入します。 y = 8.7 x − 10 {\displaystyle y=8.7x-10} x = 0.115 y + 1.15 {\displaystyle x=0.115y+1.15} S ≈ ( a / 10 + 1.2 ) ⋅ 10 n {\displaystyle {\sqrt {S}}\approx (a/10+1.2)\cdot 10^{n}}
これは、区間 における関数 y=x 2 の単一部分線形近似によって得られる 平均的に 最良の推定値である 。a=100 において最大絶対誤差は1.2、S=1および10において最大相対誤差は30%である。 [注 2] [ 1 , 100 ] {\displaystyle [1,100]}
10で割るには、 の指数から1を引くか 、比喩的に小数点を1桁左に移動します。この定式化では、任意の加法定数1に小さな増分を加えることで十分な推定値が得られるため、正確な数値を覚えておくことは負担になりません。範囲を1本の直線で結んだ近似値(四捨五入の有無にかかわらず)は、有効桁精度1桁未満です。相対誤差は1/2 2 より大きいため 、得られる情報は2ビット未満です。範囲が2桁と、この種の推定値としては非常に大きいため、精度は大幅に制限されます。 a {\displaystyle a} [ 1 , 100 ] {\displaystyle [1,100]}
区分線形近似(複数の線分で、それぞれが元の円弧の一部を近似する)を用いることで、はるかに正確な推定値を得ることができます。使用する線分の数が多いほど、近似値は向上します。最も一般的な方法は接線を用いることですが、重要なのは円弧をどのように分割するか、そして接線をどこに配置するかです。y = 1から y = 100 までの円弧を 分割 する効果的な方法 は、幾何学的に分割することです。2つの区間の境界は、元の区間の境界である 1×100 の平方根、つまり [1, 2 √ 100 ] と [ 2 √ 100 ,100] となります。 3つの区間の場合、境界 は 100 の立方根、すなわち[1, 3√100 ] 、[ 3√100 ,( 3√100)2]、[(3√100) 2,100 ] など です 。2 つ の 区間の場合、 2√100 = 10となり 、非常に便利な数 です 。接線は簡単に導出でき、x = √1 * √10 および x = √10 * √10 にあります 。これらの式は、 および です 。これを反転すると、平方根は、およびです 。 したがって 、 の 場合 : y = 3.56 x − 3.16 {\displaystyle y=3.56x-3.16} y = 11.2 x − 31.6 {\displaystyle y=11.2x-31.6} x = 0.28 y + 0.89 {\displaystyle x=0.28y+0.89} x = .089 y + 2.8 {\displaystyle x=.089y+2.8} S = a ⋅ 10 2 n {\displaystyle S=a\cdot 10^{2n}} S ≈ { ( 0.28 a + 0.89 ) ⋅ 10 n if a < 10 , ( .089 a + 2.8 ) ⋅ 10 n if a ≥ 10. {\displaystyle {\sqrt {S}}\approx {\begin{cases}(0.28a+0.89)\cdot 10^{n}&{\text{if }}a<10,\\(.089a+2.8)\cdot 10^{n}&{\text{if }}a\geq 10.\end{cases}}}
絶対誤差の最大値は区間の最大値である a =10および100で発生し、それぞれ0.54および1.7です。相対誤差の最大値は区間の両端である a =1、10、100で発生し、どちらの場合も17%です。17%、つまり0.17は1/10より大きいため、この方法では小数点以下の精度しか得られません。
双曲線推定 場合によっては、双曲線推定が有効な場合があります。双曲線も凸曲線であり、 直線よりも y = x 2 の円弧に沿っている可能性が高いためです。双曲線推定は浮動小数点除算を必要とするため、計算が複雑になります。区間 x 2 の近似最適双曲線近似は、y = 190/(10−x) − 20です。転置すると、平方根は x = 10 − 190/(y+20)です。したがって 、 [ 1 , 100 ] {\displaystyle [1,100]} S = a ⋅ 10 2 n {\displaystyle S=a\cdot 10^{2n}} S ≈ ( 10 − 190 a + 20 ) ⋅ 10 n {\displaystyle {\sqrt {S}}\approx \left(10-{\frac {190}{a+20}}\right)\cdot 10^{n}}
除算の精度は小数点1桁までで十分です。なぜなら、全体的な推定値はその精度しかなく、暗算で行うことができるからです。この双曲線推定値は、スカラー推定値や線形推定値よりも平均的に優れています。絶対誤差は a = 100で最大1.58、相対誤差は a = 10で最大となり、推定値3.67は3.16のルートより16.0%高くなります。ニュートン・ラプソン法を推定値10から開始して実行すれば、双曲線推定値と一致する3.66に到達するまでに2回の反復が必要です。75のようなより一般的なケースでは、双曲線推定値8.00はわずか7.6%低いだけであり、より正確な結果を得るには75から5回のニュートン・ラプソン法の反復が必要になります。
算術推定 区分線形近似に似ていますが、代数方程式ではなく算術演算のみを使用する方法では、掛け算表を逆に使用します。1から100までの数の平方根は1から10までの範囲にあります。したがって、25が完全な平方数(5 × 5)、36が完全な平方数(6 × 6)であることが分かっている場合、25以上36未満の数の平方根は5で始まります。他の平方数の間の数についても同様です。この方法では最初の桁は正しく計算されますが、1桁まで正確ではありません。例えば、35の平方根の最初の桁は5ですが、35の平方根はほぼ6です。
より良い方法は、範囲を正方形の中間で区切ることです。つまり、25から36の中間値である30.5までの数を5と推定し、30.5より大きく36までの数を6と推定します。 [注3] この手順では、九九から2つの積の中間の境界値を求めるのに少し計算をするだけで済みます。境界値の参照表を以下に示します。
1つの 最も近い正方形 k = a {\displaystyle k={\sqrt {a}}} EST(東部基準時。 1 1 (= 1 2 ) 1 2.5 4 (= 2 2 ) 2 6.5 9 (= 3 2 ) 3 12.5 16 (= 4 2 ) 4 20.5 25 (= 5 2 ) 5 30.5 36 (= 6 2 ) 6 42.5 49 (= 7 2 ) 7 56.5 64 (= 8 2 ) 8 72.5 81 (= 9 2 ) 9 90.5 100 (= 10 2 ) 10 100
最後の演算は推定値 k に10の累乗を2で割った値を掛けることです。 つまり、 S = a ⋅ 10 2 n {\displaystyle S=a\cdot 10^{2n}} S ≈ k ⋅ 10 n {\displaystyle {\sqrt {S}}\approx k\cdot 10^{n}}
この方法では、最適な最初の桁に丸められるため、暗黙的に 1 桁の精度が得られます。
この方法は、ほとんどの場合、被演算子を囲む最も近い正方形間の補間によって3桁まで拡張できます。 の場合 、 はおおよそ k に分数を加えた値となり、 a と k 2 の差を2つの正方形の差で割った値となります。 ここで k 2 ≤ a < ( k + 1 ) 2 {\displaystyle k^{2}\leq a<(k+1)^{2}} a {\displaystyle {\sqrt {a}}} a ≈ k + R {\displaystyle {\sqrt {a}}\approx k+R} R = ( a − k 2 ) ( k + 1 ) 2 − k 2 {\displaystyle R={\frac {(a-k^{2})}{(k+1)^{2}-k^{2}}}}
最終的な演算は、上記のように、結果に 10 の累乗を 2 で割って掛けることです。 S = a ⋅ 10 n ≈ ( k + R ) ⋅ 10 n {\displaystyle {\sqrt {S}}={\sqrt {a}}\cdot 10^{n}\approx (k+R)\cdot 10^{n}}
k は小数で、 R は小数に変換する必要がある分数です。通常、分子は1桁、分母は1桁か2桁なので、小数への変換は暗算で行うことができます。
例: 75 の平方根を求めます。 75 = 75 × 10 2 · 0 なので、 a は 75、 n は 0 です。九九から、仮数の平方根は 8 の小数点 何位でなければなりません。なぜなら、 a は 8×8 = 64 と 9×9 = 81 の間にある ため、 k は 8 です。 何位は R の小数表現です 。 分数 R は、分子が 75 − k 2 = 11、分母が 81 − k 2 = 17 です。 11/17 は、 12/18 = 2/3 = .67 より少し小さいので、.66 と推測します (ここでは推測しても問題ありません。誤差は非常に小さいです)。最終的な推定値は 8 + .66 = 8.66 です。 √ 75 を 3 桁有効にすると 8.66 なので、推定値は 3 桁有効で適切です。この方法を使用した推定値はすべて正確になるわけではありませんが、近い値になります。
バイナリ推定 コンピュータ内部で二進法 で作業する場合、 と 表すと 、平方根は 次 のように推定できる。 S {\displaystyle S} a × 2 2 n {\displaystyle a\times 2^{2n}} 0.1 2 ≤ a < 10 2 {\displaystyle 0.1_{2}\leq a<10_{2}} S = a × 2 n {\displaystyle {\sqrt {S}}={\sqrt {a}}\times 2^{n}} S ≈ ( 0.485 + 0.485 ⋅ a ) ⋅ 2 n {\displaystyle {\sqrt {S}}\approx (0.485+0.485\cdot a)\cdot 2^{n}}
これは、3桁の係数を持つ最小二乗回帰直線です。 最大絶対誤差は =2 で0.0408、最大相対誤差は =1 で3.0%です。係数が2のべき乗であるため、計算上都合の良い丸め推定値は以下のとおりです。 a {\displaystyle {\sqrt {a}}} a {\displaystyle a} a {\displaystyle a}
S ≈ ( 0.5 + 0.5 ⋅ a ) ⋅ 2 n {\displaystyle {\sqrt {S}}\approx (0.5+0.5\cdot a)\cdot 2^{n}} [注4] 最大絶対誤差は2で0.086、最大相対誤差は a = 0.5 および a = 2.0 で6.1%となります。
の場合 、2進近似は となるため 、 推定値の絶対誤差は19、相対誤差は5.3%となります。相対誤差は1/2 4 よりわずかに小さいため、推定値は4ビット以上まで良好です。 S = 125348 = 1 1110 1001 1010 0100 2 = 1.1110 1001 1010 0100 2 × 2 16 {\displaystyle S=125348=1\;1110\;1001\;1010\;0100_{2}=1.1110\;1001\;1010\;0100_{2}\times 2^{16}\,} S ≈ ( 0.5 + 0.5 ⋅ a ) ⋅ 2 8 = 1.0111 0100 1101 0010 2 ⋅ 1 0000 0000 2 = 1.456 ⋅ 256 = 372.8 {\displaystyle {\sqrt {S}}\approx (0.5+0.5\cdot a)\cdot 2^{8}=1.0111\;0100\;1101\;0010_{2}\cdot 1\;0000\;0000_{2}=1.456\cdot 256=372.8} 125348 = 354.0 {\displaystyle {\sqrt {125348}}=354.0}
8ビットまでのの推定値は、 の上位8ビットをテーブル参照することで得られます 。ただし、ほとんどの浮動小数点表現では上位ビットは暗黙的に指定され、8ビットの下位ビットは丸められる必要があることに注意してください。テーブルは、256バイトの事前計算された8ビット平方根値です。例えば、 1.8515625 10を表すインデックス 11101101 2 の場合、エントリは 10101110 2で、これは 1.8515625 10 の平方根である 1.359375 10 を表します (8ビット精度(2桁以上の小数点以下桁数))。 a {\displaystyle a} a {\displaystyle a}
ヘロン法 近似値を求める 最初の明示的な アルゴリズムは 、 1世紀 ギリシャの 数学者 アレクサンドリアのヘロンにちなんで ヘロン法 と呼ばれ、彼は 西暦60年の 著書 『メトリカ』 でこの方法を説明しました 。 バビロニア法 とも呼ばれます(斜辺を近似するバビロニア法と混同しないでください)が、この方法が バビロニア人 に知られていたという証拠はありません 。 S {\displaystyle \ {\sqrt {S~}}\ }
正の実数 が与えられたとき 、 x 0 > 0 を任意の正の初期推定値とする。ヘロン法は、 所望の精度が達成されるまで繰り返し計算を行う。この式で定義される数列は、 次のように収束する。 S {\displaystyle S} x n + 1 = 1 2 ( x n + S x n ) , {\displaystyle x_{n+1}={\frac {1}{2}}\left(x_{n}+{\frac {S}{x_{n}}}\right),} ( x 0 , x 1 , x 2 , x 3 , … ) {\displaystyle \ {\bigl (}\ x_{0},\ x_{1},\ x_{2},\ x_{3},\ \ldots \ {\bigr )}\ } lim n → ∞ x n = S . {\displaystyle \ \lim _{n\to \infty }x_{n}={\sqrt {S~}}~.}
これはニュートン法を 用いてを解くこと と同等である 。このアルゴリズムは 二次収束性が あり、 の正しい桁数は 反復ごとにほぼ倍増する。 x 2 − S = 0 {\displaystyle x^{2}-S=0} x n {\displaystyle x_{n}}
導出 基本的な考え方は、 正の実数の平方根に対して が過大評価されている場合 、 は過小評価され、その逆もまた同様であるため、これら2つの数値の平均はより正確な近似値を提供すると合理的に期待できるというものです。(この主張の正式な証明は、 算術平均と幾何平均の不等式に依存しており、 平方根 に関する記事で述べたように、この平均は常に平方根に対して過大評価され、 収束が保証されることを示しています。) x {\displaystyle \ x\ } S {\displaystyle \ S\ } S x {\displaystyle \ {\tfrac {\ S\ }{x}}\ }
より正確には、 が の初期推定値であり 、 が推定値の誤差である場合 、二項式を次のように展開することができます。
そして、誤差項を解くと、 x {\displaystyle \ x\ } S {\displaystyle \ {\sqrt {S~}}\ } ε {\displaystyle \ \varepsilon \ } S = ( x + ε ) 2 , {\displaystyle \ S=\left(x+\varepsilon \right)^{2}\ ,} ( x + ε ) 2 = x 2 + 2 x ε + ε 2 {\displaystyle \ {\bigl (}\ x+\varepsilon \ {\bigr )}^{2}=x^{2}+2x\varepsilon +\varepsilon ^{2}}
ε = S − x 2 2 x + ε ≈ S − x 2 2 x , {\displaystyle \varepsilon ={\frac {\ S-x^{2}\ }{\ 2x+\varepsilon \ }}\approx {\frac {\ S-x^{2}\ }{2x}}\ ,} もし仮に ε ≪ x {\displaystyle \ \varepsilon \ll x~} したがって、誤差を補正し、以前の推定値を次のように更新することができます
。計算された誤差は正確ではないため、これは実際の答えではありませんが、次回の補正で使用する新しい推定値となります。この更新プロセスは、所望の精度が得られるまで繰り返されます。 x + ε ≈ x + S − x 2 2 x = S + x 2 2 x = S x + x 2 ≡ x r e v i s e d . {\displaystyle \ x+\varepsilon \ \approx \ x+{\frac {\ S-x^{2}\ }{2x}}\ =\ {\frac {\ S+x^{2}\ }{2x}}\ =\ {\frac {\ {\frac {S}{\ x\ }}+x\ }{2}}\ \equiv \ x_{\mathsf {revised}}~.}
このアルゴリズムはp 進数 でも同様に機能します が、実数の平方根を p 進平方根と同一視するためには使用できません。たとえば、この方法を使用すると、実数では +3 に収束しますが、 2 進数で
は -3 に収束する有理数のシーケンスを構築できます。
小数点 以下の インポート Decimal 、 localcontext 、 getcontext NumberType = int | float | Decimal デフ sqrt_Heron ( s : 数値型 、 精度 : int | なし = なし 、 推測 : NumberType | なし = なし ) -> 小数点 : 「」 任意の精度でヘロン・ニュートン法を使用して sqrt(s) を計算します。 パラメータ ---------- s : int | float | 小数点 平方根を計算する非負数。 精度: int、オプション 有効桁数。デフォルトは現在の小数点コンテキストです。 サポートされる最小精度は 2 です。 (精度 = 1 は、丸めの異常を回避するために許可されません。) 推測: int | float | Decimal、オプション 初期推測値。デフォルトは s / 2 です。 返品 ------- 小数点 指定された精度に丸められた sqrt(s) の近似値。 「」 s == 0 の 場合 : 小数点 ( 0 ) を返す s = 小数 ( s ) s < 0 の 場合 : ValueError を発生させます ( 「sqrt(s) は負の数に対して定義されていません。」 ) 精度 が None の 場合 : precision = getcontext () . prec # 指定されていない場合は現在のグローバルコンテキストを使用する # 最小精度を暗黙的に強制する 精度 == 1 の 場合 : 精度 = 2 guess = Decimal ( s / 2 if guess が None else guess ) guard = 25 # 内部の安定性のための一時的な追加桁 最大反復回数 = 10_000 # ローカルコンテキスト: 精度の変化を分離する localcontext ()を ctx として 使用 : ctx . prec = 精度 + ガード 推測 = ( 推測 + s / 推測 ) / 2 _が 範囲 内 ( max_iter ) の場合: 次の推測 = ( 推測 + s / 推測 ) / 2 # 改善が十分に小さければ停止する 推測値 - next_guess < Decimal ( f "1e- { precision } " ) の場合 : 壊す 推測 = 次の推測 それ以外 : ArithmeticError が発生します ( f "Heron法は { max_iter } 回の反復で収束しませんでした" ) # 目標精度に丸める(ガードを外す) ctx . prec = 精度 戻り値 + 次の推測
計算例 次の例は、 sqrt_Heronさまざまな入力で関数
を実行する方法を示しています。
print ( f "1) { sqrt_Heron ( 125348 , 精度 = 7 , 推測値 = 600 ) } " ) print ( f "2) { sqrt_Heron ( Decimal ( '3.1415926535897932384626433832795028841971693993' )) } " ) print ( f "3) { sqrt_Heron ( 2 , 1_000_157 ) } " ) print ( f "4) { sqrt_Heron ( 2 , 10_000_005 , 1.414 ) } " ) print ( f "5) { sqrt_Heron ( 2 , 100_000_000 , 1 ) } " ) これにより、次の出力が生成されます。
1) 354.0452 2) 1.772453850905516027298167483 3) 1.4142135623730950488016887242 ... 269732025731849141493880004856742892 4) 1.4142135623730950488016887242 ... ... 872480508054123572727872131589714262 5) 1.4142135623730950488016887242 ... ... ... 023678977744844723443287604232894971 広告1)
を7 桁の有効数字 で計算する 手順は次のとおりです。 s {\displaystyle {\sqrt {s\,}}} s = 125348 {\displaystyle s=125348} x 0 = 6 ⋅ 10 2 = 600 x 1 = 1 2 ( x 0 + S x 0 ) = 1 2 ( 600 456666666 + 125348 600 4566666 ) ≈ 404.456666666 x 2 = 1 2 ( x 1 + S x 1 ) = 1 2 ( 404.456666666 + 125348 404.456666666 ) ≈ 357.186837334 x 3 = 1 2 ( x 2 + S x 2 ) = 1 2 ( 357.186837334 + 125348 357.186837334 ) ≈ 354.059011038 x 4 = 1 2 ( x 3 + S x 3 ) = 1 2 ( 354.059011038 + 125348 354.059011038 ) ≈ 354.045195124 x 5 = 1 2 ( x 4 + S x 4 ) = 1 2 ( 354.045195124 + 125348 354.045195124 ) ≈ 354.045194855 {\displaystyle {\begin{alignedat}{5}x_{0}&=6\cdot 10^{2}=600\\x_{1}&={\frac {1}{2}}\left(x_{0}+{\frac {S}{x_{0}}}\right)&&={\frac {1}{2}}\left(600{\phantom {456666666}}+{\frac {125348}{600}}{\phantom {4566666}}\right)&&\approx 404.456666666\\x_{2}&={\frac {1}{2}}\left(x_{1}+{\frac {S}{x_{1}}}\right)&&={\frac {1}{2}}\left(404.456666666+{\frac {125348}{404.456666666}}\right)&&\approx 357.186837334\\x_{3}&={\frac {1}{2}}\left(x_{2}+{\frac {S}{x_{2}}}\right)&&={\frac {1}{2}}\left(357.186837334+{\frac {125348}{357.186837334}}\right)&&\approx 354.059011038\\x_{4}&={\frac {1}{2}}\left(x_{3}+{\frac {S}{x_{3}}}\right)&&={\frac {1}{2}}\left(354.059011038+{\frac {125348}{354.059011038}}\right)&&\approx 354.045195124\\x_{5}&={\frac {1}{2}}\left(x_{4}+{\frac {S}{x_{4}}}\right)&&={\frac {1}{2}}\left(354.045195124+{\frac {125348}{354.045195124}}\right)&&\approx 354.045194855\end{alignedat}}}
したがって、 有効数字は 7 桁(四捨五入)になります。 125348 ≈ 354.0452 {\displaystyle {\sqrt {\,125348\,}}\approx 354.0452}
2) デフォルトの精度 での計算(6回の反復ステップ) 。 [注 5] ⌊ π × 10 46 ⌋ × 10 − 46 {\displaystyle {\sqrt {\left\lfloor \pi \times 10^{46}\right\rfloor \times 10^{-46}}}}
3) 1,000,157桁まで の計算(22回の反復ステップ) 。 2 {\displaystyle {\sqrt {2}}}
4) 10,000,005桁まで の計算(23回の反復ステップ) 。 2 {\displaystyle {\sqrt {2}}}
広告 5) 1億桁まで の計算(28回の反復ステップ) 。 [ 引用が必要 ] 2 {\displaystyle {\sqrt {2}}}
妥当な初期推測であれば、それほど多くの反復は必要ないようです。
注記
49行目と54行目の説明 ヘロン法には次のような特性があります。
x n < S ⟹ x n + 1 > S x n = S ⟹ x n + 1 = x n x n > S ⟹ x n + 1 < x n {\displaystyle {\begin{array}{rcl}x_{n}<{\sqrt {S}}&\implies &x_{n+1}>{\sqrt {S}}\\x_{n}={\sqrt {S}}&\implies &x_{n+1}=x_{n}\\x_{n}>{\sqrt {S}}&\implies &x_{n+1}<x_{n}\end{array}}} 簡単に言うと、反復処理によって より大きい値が生成されると ( の場合は直ちに 、 の場合は 1 ステップ後に発生します )、後続の推定値はすべて を上回った まま、毎回小さくなっていきます。つまり、シーケンスは に向かって「滑り落ち」、 収束します。 S {\displaystyle {\sqrt {S}}} x 0 > S {\displaystyle x_{0}>{\sqrt {S}}} x 0 < S {\displaystyle x_{0}<{\sqrt {S}}} S {\displaystyle {\sqrt {S}}} S {\displaystyle {\sqrt {S}}}
プログラムの49行目では、 guessは値 に設定されています 。そしてコードの54行目では、 は負の値にはなりません。 ≥ S {\displaystyle \geq {\sqrt {S}}} δ n = x n − x n + 1 {\displaystyle \delta _{n}=x_{n}-x_{n+1}}
停止基準の正当性 連続推定値の差を利用することで、
δ n = x n − x n + 1 {\displaystyle \delta _{n}=x_{n}-x_{n+1}} 、 停止基準として、この方法は近似値の列が 真の値に向かって収束することを保証する 。連続する差が十分に小さくなると、与えられた目標が達成される。重要な洞察は、 絶対誤差 が x n {\displaystyle x_{n}} S {\displaystyle {\sqrt {S}}} δ n {\displaystyle \delta _{n}}
ε n = x n − S {\displaystyle \varepsilon _{n}=x_{n}-{\sqrt {S}}} 、 は、逐次改善の大きさに直接関係している。具体的には、線形収束または二次収束する反復法の場合、 次のような 定数が存在する。 δ n {\displaystyle \delta _{n}} C < 1 {\displaystyle C<1}
ε n + 1 = C ⋅ δ n {\displaystyle \varepsilon _{n+1}=C\cdot \delta _{n}} 。 この関係は 、が減少するにつれて絶対誤差 も小さくなることを意味します。したがって、が所定の閾値を下回った時点で反復処理を停止することで 、実際の誤差が最大でもその閾値内に収まることが保証されます。 δ n {\displaystyle \delta _{n}} ε n {\displaystyle \varepsilon _{n}} δ n {\displaystyle \delta _{n}}
収束 100 の平方根を求めるヘロン法の収束速度を、 初期推定値を変えて比較した片対数グラフです。負の推定値は負の平方根に収束し、正の推定値は正の平方根に収束します。平方根に近い値ほど収束が速く、近似値はすべて過大評価されていることに注意してください。SVGファイルでは、グラフにマウスポインターを合わせると、その点が表示されます。 と仮定すると、 任意の自然数に対して、 相対 誤差 が
次のよう に定義され 、 x 0 > 0 a n d S > 0 . {\displaystyle \ x_{0}>0~~{\mathsf {and}}~~S>0~.} n : x n > 0 . {\displaystyle \ n:x_{n}>0~.} x n {\displaystyle \ x_{n}\ } ε n = x n S − 1 > − 1 {\displaystyle \ \varepsilon _{n}={\frac {~x_{n}\ }{\ {\sqrt {S~}}\ }}-1>-1\ } x n = S ⋅ ( 1 + ε n ) . {\displaystyle \ x_{n}={\sqrt {S~}}\cdot \left(1+\varepsilon _{n}\right)~.}
すると、 ε n + 1 = ε n 2 2 ( 1 + ε n ) ≥ 0 . {\displaystyle \ \varepsilon _{n+1}={\frac {\varepsilon _{n}^{2}}{2(1+\varepsilon _{n})}}\geq 0~.}
したがって 、収束が保証され、 二次方程式 となります。 ε n + 2 ≤ min { ε n + 1 2 2 , ε n + 1 2 } {\displaystyle \ \varepsilon _{n+2}\leq \min \left\{\ {\frac {\ \varepsilon _{n+1}^{2}\ }{2}},{\frac {\ \varepsilon _{n+1}\ }{2}}\ \right\}\ }
収束の最悪のケース 上記の大まかな推定をバビロニア方式で使用する場合、最も精度の低いケースは昇順で次のようになります。 S = 1 ; x 0 = 2 ; x 1 = 1.250 ; ε 1 = 0.250 . S = 10 ; x 0 = 2 ; x 1 = 3.500 ; ε 1 < 0.107 . S = 10 ; x 0 = 6 ; x 1 = 3.833 ; ε 1 < 0.213 . S = 100 ; x 0 = 6 ; x 1 = 11.333 ; ε 1 < 0.134 . {\displaystyle {\begin{aligned}S&=\ 1\ ;&x_{0}&=\ 2\ ;&x_{1}&=\ 1.250\ ;&\varepsilon _{1}&=\ 0.250~.\\S&=\ 10\ ;&x_{0}&=\ 2\ ;&x_{1}&=\ 3.500\ ;&\varepsilon _{1}&<\ 0.107~.\\S&=\ 10\ ;&x_{0}&=\ 6\ ;&x_{1}&=\ 3.833\ ;&\varepsilon _{1}&<\ 0.213~.\\S&=\ 100\ ;&x_{0}&=\ 6\ ;&x_{1}&=\ 11.333\ ;&\varepsilon _{1}&<\ 0.134~.\end{aligned}}}
いずれにせよ、 ε 1 ≤ 2 − 2 . ε 2 < 2 − 5 < 10 − 1 . ε 3 < 2 − 11 < 10 − 3 . ε 4 < 2 − 23 < 10 − 6 . ε 5 < 2 − 47 < 10 − 14 . ε 6 < 2 − 95 < 10 − 28 . ε 7 < 2 − 191 < 10 − 57 . ε 8 < 2 − 383 < 10 − 115 . {\displaystyle {\begin{aligned}\varepsilon _{1}&\leq 2^{-2}.\\\varepsilon _{2}&<2^{-5}<10^{-1}~.\\\varepsilon _{3}&<2^{-11}<10^{-3}~.\\\varepsilon _{4}&<2^{-23}<10^{-6}~.\\\varepsilon _{5}&<2^{-47}<10^{-14}~.\\\varepsilon _{6}&<2^{-95}<10^{-28}~.\\\varepsilon _{7}&<2^{-191}<10^{-57}~.\\\varepsilon _{8}&<2^{-383}<10^{-115}~.\end{aligned}}}
丸め誤差は収束を遅くします。 大きな丸め誤差を避けるため、計算する値の望ましい精度よりも少なくとも1桁多く残しておくことをお勧めします。 x n {\displaystyle \ x_{n}\ }
ハレー法
上記のプログラムでは、推定値(強調表示された行49と51)は、
推測 = ( 推測 + s / 推測 ) / 2 # ... 次の推測 = ( 推測 + s / 推測 ) / 2 は
推測 *= ( 推測 * 推測 + 3 * s ) / ( 3 * 推測 * 推測 + s ) # ... next_guess = guess * ( guess * guess + 3 * s ) / ( 3 * guess * guess + s ) この関数は ハレー法の sqrt_Heron実装に変換され 、
x n + 1 = x n ⋅ x n 2 + 3 S 3 x n 2 + S {\displaystyle x_{n+1}=x_{n}\cdot {\frac {x_{n}^{2}+3S}{3x_{n}^{2}+S}}} ハレー法では推定値が必ずしも一方向に動くわけではないので、54行目は
abs ( 推測値 - 次の推測値 ) < Decimal ( f "1e- { 精度 } " ) の場合 : ハレー法は反復ごとに 収束速度が 速く(根への収束速度は3乗で、2乗よりも速い ) 、反復ごとに5回の乗算(除算は3回の乗算としてカウント)が必要になります。5つの例の計算はそれぞれ4、14、15、19回の反復ステップで完了します。一方、ヘロン法では除算が1回、つまり乗算が3回で済むため、長期的にはヘロン法の方がわずかに優れています。
バクシャリ法 平方根の近似値を求めるこの方法は、 バクシャーリ写本 と呼ばれる古代 インドの 写本に記載されています。これは代数的にはヘロン法の2回の反復に等価であり、したがって4乗収束します。つまり、反復ごとに近似値の正しい桁数はほぼ4倍になります。 現代の記法を用いた元の表現は次のとおりです。 を計算するには 、 を の初期近似値とします 。次に、次のように繰り返します。 S {\displaystyle {\sqrt {S}}} x 0 2 {\displaystyle x_{0}^{2}} S {\displaystyle S} a n = S − x n 2 2 x n , x n + 1 = x n + a n , x n + 2 = x n + 1 − a n 2 2 x n + 1 . {\displaystyle {\begin{aligned}a_{n}&={\frac {S-x_{n}^{2}}{2x_{n}}},\\x_{n+1}&=x_{n}+a_{n},\\x_{n+2}&=x_{n+1}-{\frac {a_{n}^{2}}{2x_{n+1}}}.\end{aligned}}}
と の値は 、 ヘロン法で計算される値と全く同じです。これを確認するには、ヘロン法の2番目のステップで と を計算します。 と の定義を用いて 分子 を次のように並べ替えます。 x n + 1 {\displaystyle x_{n+1}} x n + 2 {\displaystyle x_{n+2}} x n + 2 = x n + 1 2 + S 2 x n + 1 = x n + 1 + S − x n + 1 2 2 x n + 1 {\displaystyle x_{n+2}={\frac {x_{n+1}^{2}+S}{2x_{n+1}}}=x_{n+1}+{\frac {S-x_{n+1}^{2}}{2x_{n+1}}}} x n + 1 {\displaystyle x_{n+1}} a n {\displaystyle a_{n}} S − x n + 1 2 = S − ( x n + a n ) 2 = S − x n 2 − 2 x n a n − a n 2 = S − x n 2 − ( S − x n 2 ) − a n 2 = − a n 2 . {\displaystyle {\begin{aligned}S-x_{n+1}^{2}&=S-(x_{n}+a_{n})^{2}\\&=S-x_{n}^{2}-2x_{n}a_{n}-a_{n}^{2}\\&=S-x_{n}^{2}-(S-x_{n}^{2})-a_{n}^{2}\\&=-a_{n}^{2}.\end{aligned}}}
これを用いて、整数から始めて平方根の有理近似値を構築することができます。が に近い 整数で 、 が絶対値が最小となる差である場合、最初の反復は次のように表すことができます。 x 0 = N {\displaystyle x_{0}=N} N 2 {\displaystyle N^{2}} S {\displaystyle S} d = S − N 2 {\displaystyle d=S-N^{2}} S ≈ N + d 2 N − d 2 8 N 3 + 4 N d = 8 N 4 + 8 N 2 d + d 2 8 N 3 + 4 N d = N 4 + 6 N 2 S + S 2 4 N 3 + 4 N S = N 2 ( N 2 + 6 S ) + S 2 4 N ( N 2 + S ) . {\displaystyle {\sqrt {S}}\approx N+{\frac {d}{2N}}-{\frac {d^{2}}{8N^{3}+4Nd}}={\frac {8N^{4}+8N^{2}d+d^{2}}{8N^{3}+4Nd}}={\frac {N^{4}+6N^{2}S+S^{2}}{4N^{3}+4NS}}={\frac {N^{2}(N^{2}+6S)+S^{2}}{4N(N^{2}+S)}}.}
バクシャリ法は分数根を含む任意の根の計算に一般化できる。
バクシャリ法の後半部分は、ヘロン反復法のより簡略化された形式として繰り返し使用できると考える人もいるかもしれない。例えば、 である。 しかし、これは 数値的に不安定で ある。元の入力値 を参照しない場合 、精度は の元の計算の精度によって制限され 、すぐに不十分になってしまう。 a n + 1 = − a n 2 2 x n + 1 , x n + 2 = x n + 1 + a n + 1 , a n + 2 = − a n + 1 2 2 x n + 2 , x n + 3 = x n + 2 + a n + 2 , etc. {\displaystyle {\begin{aligned}a_{n+1}&={\frac {-a_{n}^{2}}{2x_{n+1}}},&x_{n+2}&=x_{n+1}+a_{n+1},\\a_{n+2}&={\frac {-a_{n+1}^{2}}{2x_{n+2}}},&x_{n+3}&=x_{n+2}+a_{n+2},{\text{ etc.}}\end{aligned}}} S {\displaystyle S} a n {\displaystyle a_{n}}
例 ヘロン法の例と 同じ例を用いると、最初の反復では次のようになる。 S = 125348 {\displaystyle S=125348} x 0 = 600 a 0 = 125348 − 600 2 2 × 600 = − 195.5433 ≈ − 200 x 1 = 600 + ( − 200 ) = − 400 x 2 = 400 − ( − 200 ) 2 2 × 400 = − 350 {\displaystyle {\begin{alignedat}{3}x_{0}&=600\\[1ex]a_{0}&={\frac {125348-600^{2}}{2\times 600}}&&=-195.5433\approx -200\\[1ex]x_{1}&=600+(-200)&&={\phantom {-}}400\\[1ex]x_{2}&=400-{\frac {(-200)^{2}}{2\times 400}}&&={\phantom {-}}350\end{alignedat}}}
同様に、 2 回目の反復では が得られます。 ヘロン法とは異なり、 の式で は の誤差が修正されないため、 は 8 桁で計算する必要があります 。 a 2 = 125348 − 350 2 2 × 350 = 00 4.06857 x 3 = 350 + 4.06857 = 354.06857 x 4 = 354.06857 − 4.06857 2 2 × 354.06857 = 354.045194 {\displaystyle {\begin{alignedat}{3}a_{2}&={\frac {125348-350^{2}}{2\times 350}}&&={\phantom {00}}4.06857\\[1ex]x_{3}&=350+4.06857&&=354.06857\\[1ex]x_{4}&=354.06857-{\frac {4.06857^{2}}{2\times 354.06857}}&&=354.045194\end{alignedat}}} x 3 {\displaystyle x_{3}} x 4 {\displaystyle x_{4}} x 3 {\displaystyle x_{3}}
桁ごとの計算 これは、数列中の平方根の各桁を求める手法です。この手法は、 1600年頃に発表された フランソワ・ヴィエトの研究に由来しています。 [9]この手法は 二項定理 に基づいており 、本質的には逆アルゴリズムで を解くものです 。バビロニア法よりも時間がかかりますが、いくつかの利点があります。 ( x + y ) 2 = x 2 + 2 x y + y 2 {\displaystyle (x+y)^{2}=x^{2}+2xy+y^{2}}
手動で計算する方が簡単になります。 見つかったルートのすべての数字は正しいことがわかっているため、後で変更する必要はありません。 平方根の展開が終了している場合、アルゴリズムは最後の桁が見つかった時点で終了します。したがって、与えられた整数が 平方数 であるかどうかを確認するために使用できます。 このアルゴリズムはどの 基数 でも機能しますが、当然、その進行方法は選択された基数によって異なります。 デメリットは次のとおりです。
上位のルートでは管理できなくなります。 この方法では、不正確な推測や部分計算は許容されません。このようなエラーがあると、近似エラーを自己修正する ニュートン法 とは異なり、結果の次の桁がすべて間違ってしまいます。 桁ごとの計算は理論上は十分効率的ですが、ソフトウェア実装にはコストがかかりすぎます。反復処理ごとに大きな数値が処理されるため、より多くのメモリを必要としますが、正解の桁は1桁分しか進みません。そのため、アルゴリズムは桁が増えるごとに時間がかかります。 ネイピアの骨には、 このアルゴリズムの実行を支援する機能が含まれています。n乗根 シフト アルゴリズムはこの手法の一般化です。
基本原則 まず、数S の平方根を求めるケースを考えてみましょう。これは、10を底とする2桁の数 XY の平方です 。ここで、 X は十の位、 Y は一の位です。具体的には、 S は 3桁または4桁の小数点で構成されます。 S = ( 10 X + Y ) 2 = 100 X 2 + 20 X Y + Y 2 . {\displaystyle S=\left(10X+Y\right)^{2}=100X^{2}+20XY+Y^{2}.}
桁ごとのアルゴリズムを開始するために、 S の桁を右から始めて 2 桁の 2 つのグループに分割します。つまり、最初のグループは 1 桁または 2 桁になります。次に、 Xの値を、 X 2 が最初のグループ以下となる 最大の桁として決定します。次に、最初のグループと X 2 の差を計算し、それに 2 番目のグループを連結して 2 回目の反復を開始します。これは、 S から減算することと同じで 、 が残ります。 S' を 10 で割り、さらに 2X で割って 整数部分を保持し、 Y を推測します。 2X を 仮の Y と連結し 、それに Y を掛けます。推測が正しければ、これは次を計算することと同じです。したがって、 S' と結果 の差である剰余は0 になります。結果が S' より大きい場合は 、推測値を 1 下げて、余りが 0 になるまで再試行します。これは答えが完全な平方根 XY である単純なケースなので、アルゴリズムはここで停止します。 100 X 2 {\displaystyle 100X^{2}} S ′ = 20 X Y + Y 2 {\displaystyle S'=20XY+Y^{2}} ( 10 ( 2 X ) + Y ) Y = 20 X Y + Y 2 = S ′ , {\displaystyle (10(2X)+Y)Y=20XY+Y^{2}=S',}
同じ考え方は、任意の平方根計算にも適用できます。Sの平方根を n 個 の正の数 の和として表すと、 S = ( a 1 + a 2 + a 3 + ⋯ + a n ) 2 . {\displaystyle S=\left(a_{1}+a_{2}+a_{3}+\dots +a_{n}\right)^{2}.}
基本的な恒等式を繰り返し適用することで、
右辺の項は次のように展開できる。 ( x + y ) 2 = x 2 + 2 x y + y 2 , {\displaystyle (x+y)^{2}=x^{2}+2xy+y^{2},} ( a 1 + a 2 + a 3 + ⋯ + a n ) 2 = a 1 2 + 2 a 1 a 2 + a 2 2 + 2 ( a 1 + a 2 ) a 3 + a 3 2 + ⋯ + a n − 1 2 + 2 ( ∑ i = 1 n − 1 a i ) a n + a n 2 = a 1 2 + [ 2 a 1 + a 2 ] a 2 + [ 2 ( a 1 + a 2 ) + a 3 ] a 3 + ⋯ + [ 2 ( ∑ i = 1 n − 1 a i ) + a n ] a n . {\displaystyle {\begin{aligned}&(a_{1}+a_{2}+a_{3}+\dotsb +a_{n})^{2}\\=&\,a_{1}^{2}+2a_{1}a_{2}+a_{2}^{2}+2(a_{1}+a_{2})a_{3}+a_{3}^{2}+\dots +a_{n-1}^{2}+2\left(\sum _{i=1}^{n-1}a_{i}\right)a_{n}+a_{n}^{2}\\=&\,a_{1}^{2}+[2a_{1}+a_{2}]a_{2}+[2(a_{1}+a_{2})+a_{3}]a_{3}+\dots +\left[2\left(\sum _{i=1}^{n-1}a_{i}\right)+a_{n}\right]a_{n}.\end{aligned}}}
この式を使うと、s の値を順に推測することで平方根を求めることができます 。数値が 既に推測されていると仮定すると、 上記の和の右辺の m 番目の項は次のように与えられます。ここ で、 はこれまでに求められた平方根のおおよその値です。ここで、新しい推測は
それぞれ、 以降のすべての項の合計 、つまり剰余であり、 すべての に対して となるよう に初期化されます。正確な平方根が求められている 場合 は、s の合計が 平方根の適切な近似値となり、 は 近似誤差となります。 a i {\displaystyle a_{i}} a 1 , … , a m − 1 {\displaystyle a_{1},\ldots ,a_{m-1}} Y m = [ 2 P m − 1 + a m ] a m , {\displaystyle Y_{m}=\left[2P_{m-1}+a_{m}\right]a_{m},} P m − 1 = ∑ i = 1 m − 1 a i {\textstyle P_{m-1}=\sum _{i=1}^{m-1}a_{i}} a m {\displaystyle a_{m}} X m = X m − 1 − Y m , {\displaystyle X_{m}=X_{m-1}-Y_{m},} X m {\displaystyle X_{m}} Y m {\displaystyle Y_{m}} X m ≥ 0 {\displaystyle X_{m}\geq 0} 1 ≤ m ≤ n , {\displaystyle 1\leq m\leq n,} X 0 = S . {\displaystyle X_{0}=S.} X n = 0 , {\displaystyle X_{n}=0,} a i {\displaystyle a_{i}} X n {\displaystyle X_{n}}
例えば、10進法では、 は プレースホルダー、係数
となります 。平方根計算の任意のm番目の段階では、これまでに求められた近似根 と総和項は 次のように与えられます。 S = ( a 1 ⋅ 10 n − 1 + a 2 ⋅ 10 n − 2 + ⋯ + a n − 1 ⋅ 10 + a n ) 2 , {\displaystyle S=\left(a_{1}\cdot 10^{n-1}+a_{2}\cdot 10^{n-2}+\cdots +a_{n-1}\cdot 10+a_{n}\right)^{2},} 10 n − i {\displaystyle 10^{n-i}} a i ∈ { 0 , 1 , 2 , … , 9 } {\displaystyle a_{i}\in \{0,1,2,\ldots ,9\}} P m − 1 {\displaystyle P_{m-1}} Y m {\displaystyle Y_{m}} P m − 1 = ∑ i = 1 m − 1 a i ⋅ 10 n − i = 10 n − m + 1 ∑ i = 1 m − 1 a i ⋅ 10 m − i − 1 , {\displaystyle P_{m-1}=\sum _{i=1}^{m-1}a_{i}\cdot 10^{n-i}=10^{n-m+1}\sum _{i=1}^{m-1}a_{i}\cdot 10^{m-i-1},} Y m = [ 2 P m − 1 + a m ⋅ 10 n − m ] a m ⋅ 10 n − m = [ 20 ∑ i = 1 m − 1 a i ⋅ 10 m − i − 1 + a m ] a m ⋅ 10 2 ( n − m ) . {\displaystyle Y_{m}=\left[2P_{m-1}+a_{m}\cdot 10^{n-m}\right]a_{m}\cdot 10^{n-m}=\left[20\sum _{i=1}^{m-1}a_{i}\cdot 10^{m-i-1}+a_{m}\right]a_{m}\cdot 10^{2(n-m)}.}
ここで の位の値は 10の偶数乗なので、任意のm番目の段階において、剰余 の最上位桁のペア( その 第1項は )のみを処理すればよい。以下のセクションでは、この手順を体系化する。 Y m {\displaystyle Y_{m}} X m − 1 {\displaystyle X_{m-1}} Y m {\displaystyle Y_{m}}
同様の方法を使用して、10 進数以外の記数法で平方根を計算できることは明らかです。たとえば、2 進数で桁ごとの平方根を求めることは、 の値が {0,1} のより小さな 2 進数セットから検索されるため、非常に効率的です。これにより、各段階で の値が の場合 または の場合の いずれ かになる ため、計算が高速になります 。 の可能な選択肢が 2 つしかないという事実は、計算の m 番目の段階で の値を決定するプロセスも容易にします。これは、 について かどうかを確認するだけでよいためです。 この条件が満たされている場合は を取り、 そうでない場合は をとります。 また、2 による乗算が左ビットシフトによって行われるという事実も、計算に役立ちます。 a i {\displaystyle a_{i}} Y m {\displaystyle Y_{m}} Y m = 0 {\displaystyle Y_{m}=0} a m = 0 {\displaystyle a_{m}=0} Y m = 2 P m − 1 + 1 {\displaystyle Y_{m}=2P_{m-1}+1} a m = 1 {\displaystyle a_{m}=1} a m {\displaystyle a_{m}} a m {\displaystyle a_{m}} Y m ≤ X m − 1 {\displaystyle Y_{m}\leq X_{m-1}} a m = 1. {\displaystyle a_{m}=1.} a m = 1 {\displaystyle a_{m}=1} a m = 0. {\displaystyle a_{m}=0.}
10進数(10進数) 元の数を小数で書きましょう。数字は 長除法の アルゴリズムに似ており、長除法と同様に、根号は上の行に書きます。次に、小数点から始めて左右に数字を2つに分けます。根号の小数点は平方数の小数点の上に来ます。根号の1桁は平方数の各桁の上に来ます。
左端の数字のペアから始めて、各ペアに対して次の手順を実行します。
左から始めて、まだ使われていない最上位(左端)の2桁の数字を下ろし(すべての桁が使われている場合は「00」と書きます)、前のステップの余りの右側に書きます(最初のステップでは余りはありません)。つまり、余りに100を掛け、2桁を足します。これが 現在の値 c になります。 次のようにp 、 y 、 x を求めます 。 これまでに求めた根の一部を p と し 、小数点を無視します。(最初のステップでは p = 0です。) となる 最大の数字 x を決定します。新しい変数 y = x (20 p + x )を使用します。 x ( 20 p + x ) ≤ c {\displaystyle x(20p+x)\leq c} 注: 20 p + x は単に p を 2 倍して、右側に 数字 xを追加したものです。 注: x は、 c /(20· p ) を推測し、 y を試算して 、必要に応じて x を 上下に調整することで求めることができます。 その数字を 根号の次の数字、つまり先ほど下ろした平方根の2つの数字の上に置きます。つまり、次の p は元の p に10を掛けたものに x を 加えたものになります 。 x {\displaystyle x} c から y を引いて 新しい剰余を作成します。 余りがゼロで、下げる桁がもうない場合、アルゴリズムは終了します。そうでない場合は、ステップ1に戻って次の反復処理に進みます。
例 152.2756 の平方根を求めます。
1 2. 3 4 / \/ 01 52.27 56 01 1*1 <= 1 < 2*2 x=1 01 y = x*x = 1*1 = 1 00 52 22*2 <= 52 < 23*3 x=2 00 44 y = (20+x)*x = 22*2 = 44 08 27 243*3 <= 827 < 244*4 x=3 07 29 y = (240+x)*x = 243*3 = 729 98 56 2464*4 <= 9856 < 2465*5 x=4 98 56 y = (2460+x)*x = 2464*4 = 9856 00 00 アルゴリズム終了: 答え=12.34
二進法(基数2) このセクションでは、上記の桁ごとの計算セクションの形式論を使用しますが、わずかに 異なり、 または 各 について とします 。 から まで すべて の を反復し、 値が決定された すべての の合計で ある近似解 を構築します。 が または に 等しい かどうか を判断するために 、 とします 。 (つまり、 を含む近似解の平方が 対象の平方を超えない)場合は 、そうでない場合は 、 です 。各ステップで 平方することを避けるために 、差を保存し、 で を 設定することで増分的に更新します 。最初は、 で 最大の を 設定します 。 N 2 = ( a n + ⋯ + a 0 ) 2 {\displaystyle N^{2}=(a_{n}+\dotsb +a_{0})^{2}} a m = 2 m {\displaystyle a_{m}=2^{m}} a m = 0 {\displaystyle a_{m}=0} 2 m {\displaystyle 2^{m}} 2 n {\displaystyle 2^{n}} 2 0 {\displaystyle 2^{0}} P m = a n + a n − 1 + … + a m {\displaystyle P_{m}=a_{n}+a_{n-1}+\ldots +a_{m}} a i {\displaystyle a_{i}} a m {\displaystyle a_{m}} 2 m {\displaystyle 2^{m}} 0 {\displaystyle 0} P m = P m + 1 + 2 m {\displaystyle P_{m}=P_{m+1}+2^{m}} P m 2 ≤ N 2 {\displaystyle P_{m}^{2}\leq N^{2}} 2 m {\displaystyle 2^{m}} a m = 2 m {\displaystyle a_{m}=2^{m}} a m = 0 {\displaystyle a_{m}=0} P m = P m + 1 {\displaystyle P_{m}=P_{m+1}} P m {\displaystyle P_{m}} X m = N 2 − P m 2 {\displaystyle X_{m}=N^{2}-P_{m}^{2}} X m = X m + 1 − Y m {\displaystyle X_{m}=X_{m+1}-Y_{m}} Y m = P m 2 − P m + 1 2 = 2 P m + 1 a m + a m 2 {\displaystyle Y_{m}=P_{m}^{2}-P_{m+1}^{2}=2P_{m+1}a_{m}+a_{m}^{2}} a n = P n = 2 n {\displaystyle a_{n}=P_{n}=2^{n}} n {\displaystyle n} ( 2 n ) 2 = 4 n ≤ N 2 {\displaystyle (2^{n})^{2}=4^{n}\leq N^{2}}
追加の最適化として、 がゼロでない場合 の 2 つの項で ある と を 別々の変数 、に格納します 。 P m + 1 2 m + 1 {\displaystyle P_{m+1}2^{m+1}} ( 2 m ) 2 {\displaystyle (2^{m})^{2}} Y m {\displaystyle Y_{m}} a m {\displaystyle a_{m}} c m {\displaystyle c_{m}} d m {\displaystyle d_{m}} c m = P m + 1 2 m + 1 {\displaystyle c_{m}=P_{m+1}2^{m+1}} d m = ( 2 m ) 2 {\displaystyle d_{m}=(2^{m})^{2}} Y m = { c m + d m if a m = 2 m 0 if a m = 0 {\displaystyle Y_{m}={\begin{cases}c_{m}+d_{m}&{\text{if }}a_{m}=2^{m}\\0&{\text{if }}a_{m}=0\end{cases}}}
c m {\displaystyle c_{m}} 各ステップで 効率的に更新できます。 d m {\displaystyle d_{m}} c m − 1 = P m 2 m = ( P m + 1 + a m ) 2 m = P m + 1 2 m + a m 2 m = { c m / 2 + d m if a m = 2 m c m / 2 if a m = 0 {\displaystyle c_{m-1}=P_{m}2^{m}=(P_{m+1}+a_{m})2^{m}=P_{m+1}2^{m}+a_{m}2^{m}={\begin{cases}c_{m}/2+d_{m}&{\text{if }}a_{m}=2^{m}\\c_{m}/2&{\text{if }}a_{m}=0\end{cases}}} d m − 1 = d m 4 {\displaystyle d_{m-1}={\frac {d_{m}}{4}}}
注意してください: これは以下の関数で返される最終結果です。 c − 1 = P 0 2 0 = P 0 = N , {\displaystyle c_{-1}=P_{0}2^{0}=P_{0}=N,}
実装 Python プログラムは を計算します 。このアルゴリズムは 整数平方根を 桁ごと(ビットごと) に計算する 方法 です。 isqrt ( n ) = ⌊ n ⌋ . {\displaystyle \operatorname {isqrt} (n)=\lfloor {\sqrt {n}}\rfloor .}
def isqrt ( x : int ) -> int : assert x >= 0 , "sqrt 入力は非負である必要があります" op : int = x # X_(n+1) res : int = 0 # c_n # d_n は 4 の最大の累乗 <= n から始まります one : int = 1 while one <= op : one <<= 2 # これで 'one' が 4 の最大の累乗になります <= x one >>= 2 # dₙ … d₀ の 場合、 one != 0 の 場合 : op >= res + one の場合 : # X_(m+1) ≥ Y_m の場合、a_m = 2^m op -= res + one # X_m = X_(m+1) - Y_m res += 2 * one # c_m = c_m + 2*d_m res //= 2 # c_(m-1) = c_m / 2 one //= 4 # d_(m-1) = d_m / 4 # c_(-1) 戻り 値 2進数や10進数、その他の基数を用いたより高速なアルゴリズムは、ルックアップテーブルを使用することで実現できます。これは、実質的に、 より多くのストレージスペースと引き換えに実行時間を短縮することになります 。
指数関数的恒等式 ポケット電卓に は通常、 指数関数 と 自然対数 を計算するための優れたルーチンが実装されており、対数( )と指数関数 ( ) の性質を用いて得られる恒等式を用いて S の平方根を計算します 。 分数の分母は n乗根に対応します。上記の例では分母が2であるため、この式は平方根を求めることを示しています。 対数表 や 計算尺 を用いて平方根を計算する場合も、同じ恒等式が使用されます 。 ln x n = n ln x {\displaystyle \ln x^{n}=n\ln x} e ln x = x {\displaystyle e^{\ln x}=x} S = e 1 2 ln S . {\displaystyle {\sqrt {S}}=e^{{\frac {1}{2}}\ln S}.}
2変数反復法 この方法は の平方根を求めるのに適用でき 、 は に対して最もよく収束します 。しかし、これはコンピュータベースの計算では実際の制限にはなりません。2 を基数とする浮動小数点数および固定小数点数表現では、指数を変更するかシフトすることにより 、それぞれ 4 の整数乗、したがって 対応する 2 の累乗を乗算するのは簡単です。したがって、 は の範囲に移動できます 。さらに、以下の方法では一般的な除算は使用せず、2 の累乗による加算、減算、乗算、除算のみを使用します。これらも簡単に実装できます。この方法の欠点は、バビロニア法などの単一変数反復法とは対照的に、数値誤差が蓄積することです。 0 < S < 3 {\displaystyle 0<S<3\,\!} S ≈ 1 {\displaystyle S\approx 1} S {\displaystyle S\,\!} S {\displaystyle {\sqrt {S}}} S {\displaystyle S\,\!} 1 2 ≤ S < 2 {\textstyle {\tfrac {1}{2}}\leq S<2}
このメソッドの初期化ステップは、 反復ステップが 次に (while ) を読み取りながら行われます。 a 0 = S c 0 = S − 1 {\displaystyle {\begin{aligned}a_{0}&=S\\c_{0}&=S-1\end{aligned}}} a n + 1 = a n − a n c n / 2 c n + 1 = c n 2 ( c n − 3 ) / 4 {\displaystyle {\begin{aligned}a_{n+1}&=a_{n}-a_{n}c_{n}/2\\c_{n+1}&=c_{n}^{2}(c_{n}-3)/4\end{aligned}}} a n → S {\displaystyle a_{n}\to {\sqrt {S}}} c n → 0 {\displaystyle c_{n}\to 0}
の収束は2 次式です。 したがって の収束も 2 次式です。 c n {\displaystyle c_{n}\,\!} a n {\displaystyle a_{n}\,\!}
この方法の証明は比較的簡単です。まず、 の反復定義を と書き直します。 すると 、 が であることが帰納法によって簡単に証明できます
。 したがって、 が 目的の結果に収束することは、 が 0 に 収束することによって保証され 、これは から導かれます 。 c n {\displaystyle c_{n}} 1 + c n + 1 = ( 1 + c n ) ( 1 − 1 2 c n ) 2 . {\displaystyle 1+c_{n+1}=(1+c_{n})(1-{\tfrac {1}{2}}c_{n})^{2}.} S ( 1 + c n ) = a n 2 {\displaystyle S(1+c_{n})=a_{n}^{2}} a n {\displaystyle a_{n}\,\!} S {\displaystyle {\sqrt {S}}} c n {\displaystyle c_{n}\,\!} − 1 < c 0 < 2 {\displaystyle -1<c_{0}<2\,\!}
この手法は1950年頃にMVウィルクス 、 DJホイーラー 、 S.ギル によって、世界初の電子計算機の一つである EDSAC で使用するために開発されました 。 この手法は後に一般化され、平方根以外の計算も可能になりました。
逆平方根の反復法 以下は、 S の逆平方根を求めるための反復法です 。逆平方根が求められたら、 単純な乗算でを求めます 。これらの反復処理では乗算のみが行われ、除算は行われません。そのため、 バビロニア法 よりも高速です。ただし、安定していません。初期値が逆平方根に近くない場合、反復処理は逆平方根に収束するのではなく、そこから発散してしまいます。したがって、これらの方法を適用する前に、バビロニア法の反復処理を大まかな推定値に対して実行することが有利な場合があります。 1 / S {\displaystyle 1/{\sqrt {S}}} S {\displaystyle {\sqrt {S}}} S = S ⋅ ( 1 / S ) {\displaystyle {\sqrt {S}}=S\cdot (1/{\sqrt {S}})}
ニュートン法を この方程式に 適用すると 、ステップごとに 3 回の乗算を使用して二次収束する方法が生成されます。 ( 1 / x 2 ) − S = 0 {\displaystyle (1/x^{2})-S=0} x n + 1 = x n 2 ⋅ ( 3 − S ⋅ x n 2 ) = x n ⋅ ( 3 2 − S 2 ⋅ x n 2 ) . {\displaystyle x_{n+1}={\frac {x_{n}}{2}}\cdot (3-S\cdot x_{n}^{2})=x_{n}\cdot \left({\frac {3}{2}}-{\frac {S}{2}}\cdot x_{n}^{2}\right).} もう1つの反復は、 ハレー法(ハウス ホルダー法 の2次の方法)によって得られる 。これは 3次収束する が、反復ごとに5回の乗算が必要となる。 [ 注6] y n = S ⋅ x n 2 , {\displaystyle y_{n}=S\cdot x_{n}^{2},} x n + 1 = x n 8 ⋅ ( 15 − y n ⋅ ( 10 − 3 ⋅ y n ) ) = x n ⋅ ( 15 8 − y n ⋅ ( 10 8 − 3 8 ⋅ y n ) ) . {\displaystyle x_{n+1}={\frac {x_{n}}{8}}\cdot (15-y_{n}\cdot (10-3\cdot y_{n}))=x_{n}\cdot \left({\frac {15}{8}}-y_{n}\cdot \left({\frac {10}{8}}-{\frac {3}{8}}\cdot y_{n}\right)\right).} 固定小数点演算 を行う場合 、3倍と8除算はシフトと加算を用いて実装できます。浮動小数点を使用する場合、ハレー法は、 他のすべての定数を事前に計算して調整することで、反復ごとに4回の乗算に削減できます 。 3 / 8 S {\textstyle {\sqrt {3/8}}S} y n = 3 8 S ⋅ x n 2 , {\displaystyle y_{n}={\sqrt {\frac {3}{8}}}S\cdot x_{n}^{2},} x n + 1 = x n ⋅ ( 15 8 − y n ⋅ ( 25 6 − y n ) ) . {\displaystyle x_{n+1}=x_{n}\cdot \left({\frac {15}{8}}-y_{n}\cdot \left({\sqrt {\frac {25}{6}}}-y_{n}\right)\right).}
ゴールドシュミットのアルゴリズム ゴールドシュミットのアルゴリズムは、ロバート・エリオット・ゴールドシュミットにちなんで名付けられた ゴールドシュミット除算 の拡張であり、 [15] [16] 平方根の計算に使用できます。一部のコンピュータでは、ゴールドシュミットのアルゴリズムを使用して、 とを同時に計算します 。ゴールドシュミットのアルゴリズムは、 融合積和 命令とパイプライン型浮動小数点ユニットまたは2つの独立した浮動小数点ユニットを備えたコンピュータ上で、ニュートン・ラプソン反復 法よりも高速に計算を実行します。 S {\displaystyle {\sqrt {S}}} 1 / S {\displaystyle 1/{\sqrt {S}}} S {\displaystyle {\sqrt {S}}}
ゴールドシュミットのアルゴリズムの最初の書き方は
b 0 = S {\displaystyle b_{0}=S} Y 0 ≈ 1 / S {\displaystyle Y_{0}\approx 1/{\sqrt {S}}} (通常はテーブル検索を使用) y 0 = Y 0 {\displaystyle y_{0}=Y_{0}} x 0 = S y 0 {\displaystyle x_{0}=Sy_{0}} が1に十分近づくまで 、あるいは一定回数の反復 処理を繰り返します
。反復処理は と に収束します。 計算から と の いずれかを省略することも可能です。 また
、両方が必要な場合は 、各反復処理で を計算するのではなく、最後に を使用します。 b n + 1 = b n Y n 2 Y n + 1 = 1 2 ( 3 − b n + 1 ) x n + 1 = x n Y n + 1 y n + 1 = y n Y n + 1 {\displaystyle {\begin{aligned}b_{n+1}&=b_{n}Y_{n}^{2}\\Y_{n+1}&={\tfrac {1}{2}}(3-b_{n+1})\\x_{n+1}&=x_{n}Y_{n+1}\\y_{n+1}&=y_{n}Y_{n+1}\end{aligned}}} b i {\displaystyle b_{i}} lim n → ∞ x n = S , {\displaystyle \lim _{n\to \infty }x_{n}={\sqrt {S}},} lim n → ∞ y n = 1 / S . {\displaystyle \lim _{n\to \infty }y_{n}=1/{\sqrt {S}}.} x n {\displaystyle x_{n}} y n {\displaystyle y_{n}} x n = S y n {\displaystyle x_{n}=Sy_{n}}
2番目の形式は、 融合乗算加算 演算を使用して始まります。
y 0 ≈ 1 / S {\displaystyle y_{0}\approx 1/{\sqrt {S}}} (通常はテーブル検索を使用) x 0 = S y 0 {\displaystyle x_{0}=Sy_{0}} h 0 = 1 2 y 0 {\displaystyle h_{0}={\tfrac {1}{2}}y_{0}} が0に十分近づくまで 、あるいは一定回数の反復 を繰り返す
。これは収束し 、 r n = 0.5 − x n h n x n + 1 = x n + x n r n h n + 1 = h n + h n r n {\displaystyle {\begin{aligned}r_{n}&=0.5-x_{n}h_{n}\\x_{n+1}&=x_{n}+x_{n}r_{n}\\h_{n+1}&=h_{n}+h_{n}r_{n}\end{aligned}}} r i {\displaystyle r_{i}} lim n → ∞ x n = S , {\displaystyle \lim _{n\to \infty }x_{n}={\sqrt {S}},} lim n → ∞ 2 h n = 1 / S . {\displaystyle \lim _{n\to \infty }2h_{n}=1/{\sqrt {S}}.}
テイラー級数 N が の近似値である 場合、 平方根 関数の テイラー級数 を使用することでより良い近似値を求めることができます 。 S {\displaystyle {\sqrt {S}}} N 2 + d = N ∑ n = 0 ∞ ( − 1 ) n ( 2 n ) ! ( 1 − 2 n ) n ! 2 4 n d n N 2 n = N ( 1 + d 2 N 2 − d 2 8 N 4 + d 3 16 N 6 − 5 d 4 128 N 8 + ⋯ ) {\displaystyle {\sqrt {N^{2}+d}}=N\sum _{n=0}^{\infty }{\frac {(-1)^{n}(2n)!}{(1-2n)n!^{2}4^{n}}}{\frac {d^{n}}{N^{2n}}}=N\left(1+{\frac {d}{2N^{2}}}-{\frac {d^{2}}{8N^{4}}}+{\frac {d^{3}}{16N^{6}}}-{\frac {5d^{4}}{128N^{8}}}+\cdots \right)}
反復法であるため、 収束の次数は 使用される項の数に等しくなります。項が2つの場合、バビロニア法と同等です。項が3つの場合、各反復はバクシャーリ近似とほぼ同じ演算数を必要としますが、収束速度はより遅くなります。 [ 要出典 ] したがって、これは特に効率的な計算方法ではありません。収束速度を最大化するには、 Nを 可能な限り小さくするよう に選択します。 | d | N 2 {\displaystyle {\frac {|d|}{N^{2}}}\,}
連分数展開 実数の連分数 表現 は、その 10 進数または 2 進数展開の代わりに使用できます。この表現には、有理数 (完全な平方数ではない) の平方根が、10 進数表記法で有理数が繰り返し展開されるのと同様に、周期的な繰り返し展開を持つという特性があります。
二次無理数( a 、 b 、 c が整数である 形の数 )、特に整数の平方根は、 周期的な連分数 を持つ。平方根の数値を求めるのではなく、 連分数 展開、つまり有理数近似を求めることが必要な場合もある。 平方根を求める正の数を Sとする。a を 初期推定値、 r を剰余項とすると、次のように書ける。 となるので、 S の平方根は 次のように
表せる。 a + b c {\displaystyle {\frac {a+{\sqrt {b}}}{c}}} S = a 2 + r . {\displaystyle S=a^{2}+r.} S − a 2 = ( S + a ) ( S − a ) = r {\displaystyle S-a^{2}=({\sqrt {S}}+a)({\sqrt {S}}-a)=r} S = a + r a + S . {\displaystyle {\sqrt {S}}=a+{\frac {r}{a+{\sqrt {S}}}}.}
この式を 分数の分母に適用すると、次のようになります。 S {\displaystyle {\sqrt {S}}} S = a + r a + ( a + r a + S ) = a + r 2 a + r a + S . {\displaystyle {\sqrt {S}}=a+{\frac {r}{a+(a+{\frac {r}{a+{\sqrt {S}}}})}}=a+{\frac {r}{2a+{\frac {r}{a+{\sqrt {S}}}}}}.}
コンパクト記法
連分数の分子・分母展開(上記)は、記述もテキストフォーマットシステムへの埋め込みも面倒です。そこで数学者は、次 のような
いくつかの代替表記法を考案しました ( 一般化連分数 § 表記法 を 参照) 。 S = a + r 2 a + r 2 a + r 2 a + ⋯ {\displaystyle {\sqrt {S}}=a+{\frac {r}{2a+}}\,{\frac {r}{2a+}}\,{\frac {r}{2a+}}\cdots }
全体を通すと 、さらに簡潔な表記は次のようになります。 [注 7] 連分数の繰り返し(非完全平方の平方根はすべて繰り返します)の場合、繰り返し部分は1回だけ表され、上線は上線部分の終わりのない繰り返しを示します。 [注 8] r = 1 {\displaystyle r=1} [ a ; 2 a , 2 a , 2 a , ⋯ ] {\displaystyle [a;2a,2a,2a,\cdots ]} [ a ; 2 a ¯ ] {\displaystyle [a;{\overline {2a}}]}
√2 の 場合 、 の値は 1なので、その表現は次のようになります。 a {\displaystyle a} [ 1 ; 2 ¯ ] {\displaystyle [1;{\overline {2}}]}
この方法を続けると、平方根の 一般化された連分数は 次のように表せる。 S = a + r 2 a + r 2 a + r 2 a + ⋱ {\displaystyle {\sqrt {S}}=a+{\cfrac {r}{2a+{\cfrac {r}{2a+{\cfrac {r}{2a+\ddots }}}}}}}
このような分数 を評価して根を求める最初のステップは、根を求める数と選択した分母の数を数値的に代入することです。例えば、標準形では は1であり、√2 の場合は は 1 なので 、分母が3つの場合の数値連分数は次のようになります。 r {\displaystyle r} a {\displaystyle a} 2 ≈ 1 + 1 2 + 1 2 + 1 2 {\displaystyle {\sqrt {2}}\approx 1+{\cfrac {1}{2+{\cfrac {1}{2+{\cfrac {1}{2}}}}}}}
ステップ2では、連分数を一つずつ下から上へ約分し、分子と分母が整数である有理分数にします。約分は次のように行われます(最初の3つの分母を取ります)。 1 + 1 2 + 1 2 + 1 2 = 1 + 1 2 + 1 5 2 = 1 + 1 2 + 2 5 = 1 + 1 12 5 = 1 + 5 12 = 17 12 {\displaystyle {\begin{aligned}1+{\cfrac {1}{2+{\cfrac {1}{2+{\cfrac {1}{2}}}}}}&=1+{\cfrac {1}{2+{\cfrac {1}{\frac {5}{2}}}}}\\&=1+{\cfrac {1}{2+{\cfrac {2}{5}}}}=1+{\cfrac {1}{\frac {12}{5}}}\\&=1+{\cfrac {5}{12}}={\frac {17}{12}}\end{aligned}}}
最後に (手順 3)、有理分数の分子を分母で割って、根のおおよその値を取得します。 これは 3 桁の精度に丸められます。 17 ÷ 12 = 1.42 {\displaystyle 17\div 12=1.42}
√2 の実際の値は 、 有効数字3桁で1.41です。相対誤差は0.17%なので、有理分数はほぼ3桁の精度で計算できます。分母の数を増やすほど、より正確な近似値が得られます。分母を4つにすると分数 となり 、ほぼ4桁の精度で計算できます。 41 29 = 1.4137 {\displaystyle {\frac {41}{29}}=1.4137}
以下は、平方根、その単純な連分数、および 分母 99 までの
最初の項 ( 収束項と呼ばれる) の例です。
√S ~小数点 連分数 収束 √2 1.41421 [ 1 ; 2 ¯ ] {\displaystyle [1;{\overline {2}}]} 3 2 , 7 5 , 17 12 , 41 29 , 99 70 {\displaystyle {\frac {3}{2}},{\frac {7}{5}},{\frac {17}{12}},{\frac {41}{29}},{\frac {99}{70}}} √3 1.73205 [ 1 ; 1 , 2 ¯ ] {\displaystyle [1;{\overline {1,2}}]} 2 1 , 5 3 , 7 4 , 19 11 , 26 15 , 71 41 , 97 56 {\displaystyle {\frac {2}{1}},{\frac {5}{3}},{\frac {7}{4}},{\frac {19}{11}},{\frac {26}{15}},{\frac {71}{41}},{\frac {97}{56}}} √5 2.23607 [ 2 ; 4 ¯ ] {\displaystyle [2;{\overline {4}}]} 9 4 , 38 17 , 161 72 {\displaystyle {\frac {9}{4}},{\frac {38}{17}},{\frac {161}{72}}} √6 2.44949 [ 2 ; 2 , 4 ¯ ] {\displaystyle [2;{\overline {2,4}}]} 5 2 , 22 9 , 49 20 , 218 89 {\displaystyle {\frac {5}{2}},{\frac {22}{9}},{\frac {49}{20}},{\frac {218}{89}}} √10 3.16228 [ 3 ; 6 ¯ ] {\displaystyle [3;{\overline {6}}]} 19 6 , 117 37 {\displaystyle {\frac {19}{6}},{\frac {117}{37}}} π {\displaystyle {\sqrt {\pi }}} 1.77245 [ 1 ; 1 , 3 , 2 , 1 , 1 , 6... ] {\displaystyle [1;1,3,2,1,1,6...]} 2 1 , 7 4 , 16 9 , 23 13 , 39 22 {\displaystyle {\frac {2}{1}},{\frac {7}{4}},{\frac {16}{9}},{\frac {23}{13}},{\frac {39}{22}}} e {\displaystyle {\sqrt {e}}} 1.64872 [ 1 ; 1 , 1 , 1 , 5 , 1 , 1... ] {\displaystyle [1;1,1,1,5,1,1...]} 2 1 , 3 2 , 5 3 , 28 17 , 33 20 , 61 37 {\displaystyle {\frac {2}{1}},{\frac {3}{2}},{\frac {5}{3}},{\frac {28}{17}},{\frac {33}{20}},{\frac {61}{37}}} ϕ {\displaystyle {\sqrt {\phi }}} 1.27202 [ 1 ; 3 , 1 , 2 , 11 , 3 , 7... ] {\displaystyle [1;3,1,2,11,3,7...]} 4 3 , 5 4 , 14 11 {\displaystyle {\frac {4}{3}},{\frac {5}{4}},{\frac {14}{11}}}
一般に、有理分数の分母が大きいほど、近似値は良くなります。また、連分数を切り捨てると、その分母が分母以下の分数の根に最も近似する有理分数が得られることも示されています。例えば、分母が70以下の分数は、99 /70ほど √2に 近似 できません。
浮動小数点表現に依存する近似値 数値は 浮動小数点 形式で表され、これは 科学的記数法 とも呼ばれます 。その平方根は で 、立方根や対数にも同様の公式が適用されます。一見すると、これは簡潔さの向上にはならないように見えますが、近似値のみが必要な場合を考えてみましょう 。その場合、 は1桁まで適切です。次に、いくつかのべき乗 p が奇数になることを認識してください。つまり、3141.59 = 3.14159 × 10です。 m × b p {\displaystyle m\times b^{p}} m × b p / 2 {\displaystyle {\sqrt {m}}\times b^{p/2}} b p / 2 {\displaystyle b^{p/2}} 3 底の分数乗を扱うのではなく、仮数部に底を掛け、その乗数から1を引いて偶数にします。調整された表現は31.4159 × 10 2 なので平方根は √31.4159 × 10 となる。 1 .
調整仮数の整数部を取る場合、値は1から99までしか取れません。これは、99個の事前計算された平方根の表へのインデックスとして使用し、推定値を完成させることができます。16進数のコンピュータではより大きな表が必要になりますが、2進数のコンピュータでは3つのエントリしか必要ありません。調整仮数の整数部の可能なビットは01(べき乗が偶数なのでシフトは発生しません。 正規化された 浮動小数点数は常に上位桁が0以外の値を持つことを覚えておいてください)または、べき乗が奇数の場合は10または11で、これらは 元の仮数の最初の 2 ビットです。したがって、6.25 = 110.01 は 2 進数で 1.1001 × 2 2 に正規化されて偶数累乗になり、仮数部のペアビットは 01 になります。一方、.625 = 0.101 は 2 進数で 1.01 × 2 −1 の奇数累乗に正規化されるため、調整は 10.1 × 2 −2 になり、ペアビットは 10 になります。累乗の下位ビットがペア仮数の上位ビットに反映されることに注意してください。偶数累乗の下位ビットは 0 で、調整された仮数は 0 で始まりますが、奇数累乗の下位ビットは 1 で、調整された仮数は 1 で始まります。したがって、累乗が半分になると、下位ビットがシフトされてペア仮数部の最初のビットになるかのようになります。
3つのエントリしかないテーブルは、仮数部のビットを追加することで拡張できます。しかし、コンピュータでは、テーブルに補間を計算するよりも、同等の結果をもたらすより単純な計算を見つける方が多くの場合効果的です。すべては、表現形式の正確な詳細と、数値の各部分にアクセスして操作するためにどのような演算が利用可能かに依存します。例えば、 Fortran にはべき乗を求める関数があります EXPONENT(x)。良好な初期近似値を求めるために費やされた労力は、粗悪な近似値を求めるために必要となるであろう改良プロセスの追加の反復を回避することで回収されます。これらの反復は少ないため(1回の反復で除算、加算、2分の1が実行される)、制約は厳格です。
多くのコンピュータは IEEE (または十分に類似した)表現に従っており、ニュートン法の出発点として平方根の近似値を非常に高速に得ることができます。ここで説明する手法は、浮動小数点形式(2を底とする)が2を底とする対数を近似するという事実に基づいています。つまり、 log 2 ( m × 2 p ) = p + log 2 ( m ) {\displaystyle \log _{2}(m\times 2^{p})=p+\log _{2}(m)}
したがって、IEEE形式の32ビット単精度浮動小数点数(特に、表現された形式では累乗に127の バイアス が加えられている)の場合、そのバイナリ表現を32ビット整数として解釈し、それを でスケーリングし、127のバイアスを除去する ことで、近似対数を得ることができます。 2 − 23 {\displaystyle 2^{-23}} x int ⋅ 2 − 23 − 127 ≈ log 2 ( x ) . {\displaystyle x_{\text{int}}\cdot 2^{-23}-127\approx \log _{2}(x).}
例えば、1.0 は 16進 数 0x3F800000 で表されますが、これを 整数にすると となります。上記の式を用いると 、 から期待される が得られます 。同様に、1.5 (0x3FC00000) から 0.5 が得られます。 1065353216 = 127 ⋅ 2 23 {\displaystyle 1065353216=127\cdot 2^{23}} 1065353216 ⋅ 2 − 23 − 127 = 0 {\displaystyle 1065353216\cdot 2^{-23}-127=0} log 2 ( 1.0 ) {\displaystyle \log _{2}(1.0)}
平方根を求めるには、対数を2で割り、その値を元の値に戻します。次のプログラムはこの考え方を示しています。指数の最下位ビットは意図的に仮数部に伝播するようにしています。このプログラムの手順を正当化する一つの方法は、 指数バイアスを 、 仮数部に明示的に格納されているビット数を と仮定し、次の式を示すことです。 b {\displaystyle b} n {\displaystyle n} ( ( 1 2 ( x int / 2 n − b ) ) + b ) ⋅ 2 n = 1 2 ( x int − 2 n ) + ( 1 2 ( b + 1 ) ) ⋅ 2 n . {\displaystyle \left(\left({\tfrac {1}{2}}\left(x_{\text{int}}/2^{n}-b\right)\right)+b\right)\cdot 2^{n}={\tfrac {1}{2}}\left(x_{\text{int}}-2^{n}\right)+\left({\tfrac {1}{2}}\left(b+1\right)\right)\cdot 2^{n}.}
/* float は IEEE 754 単精度浮動小数点形式であると想定します */ #include <stdint.h> ユニオン FloatUInt { float f ; uint32_t i ; } float sqrtapprox ( float z ) { union FloatUInt val = { z }; // ビット パターンを維持しながら型を変換します /* * 次のコードを正当化するには、次の式を証明してください。 * * ((((val.i / 2^m) - b) / 2) + b) * 2^m = ((val.i - 2^m) / 2) + ((b + 1) / 2) * 2^m) * * ここで 、 * * b = 指数バイアス * m = 仮数ビット数 */ val . i -= 1 << 23 ; // 2^m を減算します。 val . i >>= 1 ; // 2 で割ります。 val . i += 1 << 29 ; // ((b + 1) / 2) * 2^m を加算します。 // float として再度解釈します return val . f ; } 上記の関数の中核を成す3つの数学演算は、1行で表現できます。さらに調整を加えることで、最大相対誤差を低減できます。したがって、キャストを除く3つの演算は、次のように書き直すことができます。
val.i = ( 1 << 29 ) + ( val.i >> 1 ) - ( 1 << 22 ) + a ; ここで、 a は近似誤差を調整するためのバイアスです。例えば、 a = 0 の場合、2の偶数乗(例:1.0)では正確な結果が得られますが、他の数値では結果がわずかに大きくなります(例:2.0の場合、1.414ではなく1.5となり、誤差は6%になります)。a = −0x4B0D2 の場合 、 最大相対誤差は±3.5%に最小化されます。
近似値を方程式に対する ニュートン法の 初期推定値として使用する場合は 、次のセクションに示す逆数形式が推奨されます。 ( 1 / x 2 ) − S = 0 {\displaystyle (1/x^{2})-S=0}
平方根の逆数 上記のルーチンの派生版を以下に掲載します。これは平方根の 逆数 を計算するために使用できます。 つまり、グレッグ・ウォルシュによって書かれたものです。整数シフト近似により相対誤差は4%未満となり、 次の行で ニュートン法を1回反復することで誤差はさらに0.15%まで減少しました。 コンピュータグラフィックスにおいて、これはベクトルを正規化する非常に効率的な方法です。 x − 1 / 2 {\displaystyle x^{-1/2}}
ユニオン FloatInt { float x ; int i ; }; float inv_sqrt ( float x ) { float xhalf = 0.5f * x ; union FloatInt u ; u . x = x ; u . i = 0x5f375a86 - ( u . i >> 1 ); // 次の行は、精度を上げるために任意の回数繰り返すことができます u . x = u . x * ( 1.5f - xhalf * u . x * u . x ); return u . x ; } 一部のVLSIハードウェアでは、2次多項式推定とそれに続くゴールドシュミット 反復法を使用して逆平方根を実装しています 。
負の平方または複素平方 S < 0の場合 、その主平方根は S = | S | i . {\displaystyle {\sqrt {S}}={\sqrt {\vert S\vert }}\,\,i\,.}
S = a + bi ( a と b は実数、 b ≠ 0)とすると 、その主平方根は S = | S | + a 2 + sgn ( b ) | S | − a 2 i . {\displaystyle {\sqrt {S}}={\sqrt {\frac {\vert S\vert +a}{2}}}\,+\,\operatorname {sgn}(b){\sqrt {\frac {\vert S\vert -a}{2}}}\,\,i\,.}
これはルートを二乗することで検証できる。 ここで | S | = a 2 + b 2 {\displaystyle \vert S\vert ={\sqrt {a^{2}+b^{2}}}}
はS の 法 である。 複素数 の主平方根は、 実部が非負である根として定義される。
参照
注記 ^ 係数 2 と 6 は、指定された桁数で可能な最小値と最大値の幾何平均を近似するため使用 さ れ ます 。 1 ⋅ 10 = 10 4 ≈ 1.78 {\displaystyle {\sqrt {{\sqrt {1}}\cdot {\sqrt {10}}}}={\sqrt[{4}]{10}}\approx 1.78\,} 10 ⋅ 100 = 1000 4 ≈ 5.62 {\displaystyle {\sqrt {{\sqrt {10}}\cdot {\sqrt {100}}}}={\sqrt[{4}]{1000}}\approx 5.62\,} ^ 四捨五入していない推定値の最大絶対誤差は100で2.65、最大相対誤差はy=1、10、100で26.5%である。 ^ 数字が2つの平方数のちょうど中間、例えば30.5の場合は、大きい方の数字を推測します。この場合は6です。 ちなみに これは、 y = 1における y = x 2 の接線の方程式です 。 ^ デフォルトの精度は1から decimal.MAX_PREC^ 上記参照。 ^ 参照: 連分数#表記法 ^ 参照: 周期連分数
参考文献 ^ Herrero Piñeyro, PJ; Linero Bas, A.; Massa Esteve, MR; Mellado Romero, A. (2023). 「Vièteの研究に基づくn根近似値に関する問題」 (PDF) . MATerials MATemàtics . 5 : 1– 27. 2025年 11月15日 閲覧 。 8ページをご覧ください。 ^ Goldschmidt, Robert E. (1964). 収束による除算の応用 (PDF) (論文). 理学修士論文. MIT OCLC 34136725. 2015年12月10日時点のオリジナルより アーカイブ (PDF) . 2015年 9月15日 閲覧 。 ^ 「著者」. IBM Journal of Research and Development . 11 : 125–127 . 1967. doi :10.1147/rd.111.0125. 2018年7月18日時点の オリジナル よりアーカイブ。
参考文献 アブラモウィッツ、ミルトン 、 ステガン 、アイリーン・A. 編 (1970) [1964]. 数式、グラフ、数表付き数学関数ハンドブック . 応用数学シリーズ55(第9刷). ワシントンD.C.: ドーバー出版 . ISBN 978-0-486-61272-0 LCCN 64-60036。MR 0167642 。 ベイリー、デイビッド、 ボルウェイン、ジョナサン (2012). 「古代インドの平方根:法医学的古数学の演習」 (PDF) . アメリカ数学月刊誌 . 第119巻第8号. 646–657頁 . 2017 年 9月14日 閲覧 . キャンベル=ケリー、マーティン (2009年9月)「コンピューティングの起源」 サイエンティフィック・アメリカン誌 、 301 (3): 62–69 . Bibcode :2009SciAm.301c..62C. doi :10.1038/scientificamerican0909-62. JSTOR 26001527. PMID :19708529. クック、ロジャー(2008年) 『古典代数学:その性質、起源、そして用途 』ジョン・ワイリー・アンド・サンズ社、 ISBN 978-0-470-25952-8 。 ファウラー、デイビッド;ロブソン、エレノア (1998). 「古バビロニア数学における平方根近似:YBC 7289の文脈」 (PDF) . Historia Mathematica . 25 (4): 376. doi : 10.1006/hmat.1998.2209 . ガワー, ジョン・C. (1958). 「根抽出のための反復法に関するノート」. コンピュータジャーナル . 1 (3): 142–143 . doi : 10.1093/comjnl/1.3.142 . ガイ、マーティン (1985). 「ウー氏のアバカスアルゴリズムによる高速整数平方根」. ケント大学カンタベリー校 (UKC). 2012年3月6日時点のオリジナルよりアーカイブ。 2025年 10月5日 閲覧 。 ヒース、トーマス (1921年)『ギリシャ数学史』第2巻、オックスフォード:クラレンドン・プレス、323-324頁。 ジャクソン、テレンス (2011年7月1日). 「95.42 自然数の無理平方根 ― 幾何学的アプローチ」 . The Mathematical Gazette . 95 (533): 327– 330. doi :10.1017/S0025557200003193. ISSN 0025-5572. S2CID 123995083. Johnson, SG (2015年2月4日). 「ニュートン法による平方根」 (PDF) . MITコース18.335: 数値解析入門 . 2025年 10月12日 閲覧. ロモント、クリス (2003).「高速逆平方根」 (PDF) . Markstein, Peter (2004年11月). Goldschmidtアルゴリズムを用いたソフトウェアによる除算と平方根 (PDF) . 第6回実数とコンピュータに関する会議. ドイツ 、ダグストゥール. CiteSeerX 10.1.1.85.9648 . ロバート・ネミロフ、ジェリー・ボンネル (1994). 「2の平方根を100万桁まで」 NASA Astronomy Picture of the Day . 2025年 10月21日 閲覧 。 ロバート・ネミロフ、ジェリー・ボンネル (1994a). 「2から1000万桁までの平方根」 NASA Astronomy Picture of the Day . 2025年 10月21日 閲覧 。 Piñeiro, José-Alejandro; Díaz Bruguera, Javier (2002年12月). 「逆数、除算、平方根、逆平方根の高速倍精度計算」 IEEE Transactions on Computers . 51 (12): 1377– 1388. Bibcode :2002ITCmp..51.1377P. doi :10.1109/TC.2002.1146704. サルディーナ、マニー (2007).「(折り畳み)連分数を用いた根の一般的な抽出法」サリー(英国) Simply Curious (2018年6月5日). 「バクシャーリー写本に屈する」. Simply Curious ブログ. 2020年 12月21日 閲覧 。 スタイナーソン、アーン。ダン、コービット。ヘンドリー、マシュー (2003)。 「整数平方根関数」。 Wilkes, MV ; Wheeler, DJ ; Gill, S. (1951). 『電子デジタル計算機のためのプログラムの作成法』Oxford: Addison-Wesley. p. 262. OCLC 475783493.
外部リンク ワイスタイン、エリック・W. 「平方根アルゴリズム」 。MathWorld 。 Jarvis, AF「減算による平方根」 (PDF) . afjarvis.staff.shef.ac.uk . 2022年4月10日時点のオリジナル (PDF) からアーカイブ。 ラドヴィッチ、アンドリヤ。 「整数平方根アルゴリズム」。 アンドリジャル.com 。 2025 年 11 月 7 日 に取得 。 エグバート, ウィリアム・E. (1977年5月). 「パーソナル電卓アルゴリズム I: 平方根」 (PDF) . ヒューレット・パッカード・ジャーナル : 22–24 . 2025年 11月7日 閲覧 。 「平方根を学ぶための計算機」. calculatorsquareroot.com . 2025年 11月7日 閲覧 。