整数平方根

数論では、非負整数nの整数平方根（isqrt）は、nの平方根以下の最大の整数である非負整数mである。 $${\displaystyle \operatorname {isqrt} (n)=\lfloor {\sqrt {n}}\rfloor .}$$

例えば、${\displaystyle \operatorname {isqrt} (27)=\lfloor {\sqrt {27}}\rfloor =\lfloor 5.19615242270663...\rfloor =5.}$

およびを非負の整数 とします。${\displaystyle y}$${\displaystyle k}$

（の10進表現）を計算するアルゴリズムは、完全な平方数ではない各入力に対して永遠に実行されます。^{[注 1]}${\displaystyle {\sqrt {y}}}$ ${\displaystyle y}$

計算アルゴリズムは永久に実行されるわけではありません。それでも、任意の精度で計算を行うことができます。 ${\displaystyle \lfloor {\sqrt {y}}\rfloor }$ ${\displaystyle {\sqrt {y}}}$${\displaystyle k}$

いずれかを選択して計算します。 ${\displaystyle k}$${\textstyle \lfloor {\sqrt {y\times 100^{k}}}\rfloor }$

たとえば（設定）： ${\displaystyle y=2}$$${\displaystyle {\begin{aligned}&k=0:\lfloor {\sqrt {2\times 100^{0}}}\rfloor =\lfloor {\sqrt {2}}\rfloor =1\\&k=1:\lfloor {\sqrt {2\times 100^{1}}}\rfloor =\lfloor {\sqrt {200}}\rfloor =14\\&k=2:\lfloor {\sqrt {2\times 100^{2}}}\rfloor =\lfloor {\sqrt {20000}}\rfloor =141\\&k=3:\lfloor {\sqrt {2\times 100^{3}}}\rfloor =\lfloor {\sqrt {2000000}}\rfloor =1414\\&\vdots \\&k=8:\lfloor {\sqrt {2\times 100^{8}}}\rfloor =\lfloor {\sqrt {20000000000000000}}\rfloor =141421356\\&\vdots \\\end{aligned}}}$$

結果を比較する${\displaystyle {\sqrt {2}}=1.41421356237309504880168872420969807856967187537694...}$

入力を乗算するとk桁の精度が得られるようです。 ^{[注 2]}${\displaystyle 100^{k}}$

の（全体の）10進表現を計算するには、各パスごとに 係数を増やしながら、無限回実行することができます。${\displaystyle {\sqrt {y}}}$${\displaystyle \operatorname {isqrt} (y)}$${\displaystyle y}$${\displaystyle 100}$

次のプログラム（）では、手順がすでに定義されており、議論のために、すべての変数が無制限の大きさの整数を保持できると仮定します。 ${\displaystyle \operatorname {sqrtForever} }$${\displaystyle \operatorname {isqrt} (y)}$

次に、の10進表現全体を出力します。^{[注 3]}${\displaystyle \operatorname {sqrtForever} (y)}$${\displaystyle {\sqrt {y}}}$

// 停止せずにsqrt(y)を出力します
void sqrtForever (符号なし整数y )   
{
	符号なし整数結果= isqrt ( y );    
	printf ( "%d." , result ); // 結果を出力し、その後に小数点を付けます 	

	while ( true ) // 永久に繰り返す ... 	
	{
		y = y * 100 ; // 理論的な例: オーバーフローは無視されます    	
		結果= isqrt ( y );  
		printf ( "%d" , result % 10 ); // 結果の最後の桁を出力   	
	}
}

結論としては、 を計算するアルゴリズムは、を計算するアルゴリズムisqrt()と計算上同等であるということです。 sqrt()

非負整数の平方根は次のように定義される。 ${\displaystyle y}$$${\displaystyle \lfloor {\sqrt {y}}\rfloor =x:x^{2}\leq y<(x+1)^{2},x\in \mathbb {N} }$$

たとえば、なぜなら。 ${\displaystyle \operatorname {isqrt} (27)=\lfloor {\sqrt {27}}\rfloor =5}$${\displaystyle 6^{2}>27{\text{ および }}5^{2}\ngtr 27}$

次の C プログラムは簡単な実装です。

// 整数の平方根
// (線形検索、昇順を使用) 
unsigned int isqrt ( unsigned int y ) { // 初期過小評価、L <= isqrt(y) unsigned int L = 0 ;    

	
	    

	(( L + 1 ) * ( L + 1 ) <= y ) L = L + 1である一方、         
		    

	Lを返す; } 

// 整数の平方根
// (線形検索、降順を使用) 
unsigned int isqrt ( unsigned int y ) { // 初期過大評価、isqrt(y) <= R unsigned int R = y ;    

	
	    

	一方、( R * R > y ) R = R - 1です。     
		    

	Rを返します; } 

上記のプログラム（線形探索、昇順）では、等価性を使用して乗算を加算に置き換えることができます。 $${\displaystyle (L+1)^{2}=L^{2}+2L+1=L^{2}+1+\sum _{i=1}^{L}2.}$$

// 整数の平方根
// (線形検索、昇順) 加算を使用
符号なし整数isqrt (符号なし整数y )    
{
	符号なし整数L = 0 ;    
	符号なし整数a = 1 ;    
	符号なし整数d = 3 ;    

	( a <= y )であるとき   
	{
		a = a + d ; // (a + 1) ^ 2    	
		d = d + 2 ;    
		L = L + 1 ;    
	}

	Lを返します。 
}

線形探索は、の最小値に達するまですべての値を順番にチェックします。 ${\displaystyle x}$${\displaystyle x^{2}>y}$

代わりにバイナリ検索を使用することで高速化が実現されます。次の C プログラムは実装です。

// 整数平方根（二分探索を使用）
符号なし整数isqrt (符号なし整数y )    
{
	符号なし整数L = 0 ;    
	符号なし整数M ;  
	符号なし整数R = y + 1 ;      

    ただし、( L != R - 1 )     
    {
        M = ( L + R ) / 2 ;      

		( M * M <= y )の場合     
			L = M ;  
		それ以外
			R = M ;  
	}

    Lを返します。 
}

数値例

例えば、二分探索法を使って計算すると、次のシーケンス が得られる。${\displaystyle \operatorname {isqrt} (2000000)}$${\displaystyle [左,右]}$$${\displaystyle {\begin{aligned}&[0,2000001]\rightarrow [0,1000000]\rightarrow [0,500000]\rightarrow [0,250000]\rightarrow [0,125000]\rightarrow [0,62500]\rightarrow [0,31250]\rightarrow [0,15625]\\&\rightarrow [0,7812]\rightarrow [0,3906]\rightarrow [0,1953]\rightarrow [976,1953]\rightarrow [976,1464]\rightarrow [1220,1464]\rightarrow [1342,1464]\rightarrow [1403,1464]\\&\rightarrow [1403,1433]\rightarrow [1403,1418]\rightarrow [1410,1418]\rightarrow [1414,1418]\rightarrow [1414,1416]\rightarrow [1414,1415]\end{aligned}}}$$

この計算には21回の反復ステップが必要ですが、線形探索（昇順、から開始）では${\displaystyle 0}$1414段。

とを計算する一つの方法は、ニュートン法の特殊なケースであるヘロン法を使って方程式の解を求めることで、反復式は次のようになる。 ${\displaystyle {\sqrt {n}}}$${\displaystyle \operatorname {isqrt} (n)}$${\displaystyle x^{2}-n=0}$$${\displaystyle x_{k+1}={\frac {1}{2}}\!\left(x_{k}+{\frac {n}{x_{k}}}\right),\quad k\geq 0,\quad x_{0}>0.}$$

この数列は として二次収束します。 ${\displaystyle \{x_{k}\}}$ ${\displaystyle {\sqrt {n}}}$${\displaystyle k\to \infty }$

上記のアルゴリズムで 停止基準 が保証する最大の可能数であることを証明することができます^{[引用が必要] 。}${\displaystyle c=1}$$${\displaystyle |x_{k+1}-x_{k}|<c}$$${\displaystyle \lfloor x_{k+1}\rfloor =\lfloor {\sqrt {n}}\rfloor }$

すべての有理数を正確に表現できない数値形式(浮動小数点など) を使用する実装では、丸め誤差を防ぐために 1 未満の停止定数を使用する必要があります。

は多くの に対して無理数ですが、 が有理数である場合、数列には有理数項のみが含まれます。したがって、この方法を使用すると、 を計算するために有理数のフィールドから出る必要はなく、この事実には理論的な利点がいくつかあります。 ${\displaystyle {\sqrt {n}}}$${\displaystyle n}$${\displaystyle \{x_{k}\}}$${\displaystyle x_{0}}$${\displaystyle \operatorname {isqrt} (n)}$

非常に大きな整数nを計算する場合、両方の除算演算にユークリッド除算の商を使用することができます。これには、中間値に整数のみを使用するという利点があり、大きな数の浮動小数点表現を使用する必要はありません。これは、反復式を使用するのと同じです。 ${\displaystyle \lfloor {\sqrt {n}}\rfloor }$$${\displaystyle x_{k+1}=\left\lfloor {\frac {1}{2}}\!\left(x_{k}+\left\lfloor {\frac {n}{x_{k}}}\right\rfloor \right)\right\rfloor ,\quad k\geq 0,\quad x_{0}>0,\quad x_{0}\in \mathbb {Z} .}$$

という事実を利用して $${\displaystyle \left\lfloor {\frac {1}{2}}\!\left(x_{k}+\left\lfloor {\frac {n}{x_{k}}}\right\rfloor \right)\right\rfloor =\left\lfloor {\frac {1}{2}}\!\left(x_{k}+{\frac {n}{x_{k}}}\right)\right\rfloor ,}$$

これは有限回数の反復で 到達できることを示すことができます。${\displaystyle \lfloor {\sqrt {n}}\rfloor }$

元のバージョンでは、に対して、に対してとなります。したがって、整数バージョンでは、最終解に達するまでおよび となります。最終解 に対して、 およびとなるため、停止基準は となります。 ${\displaystyle x_{k}\geq {\sqrt {n}}}$${\displaystyle k\geq 1}$${\displaystyle x_{k}>x_{k+1}}$${\displaystyle x_{k}>{\sqrt {n}}}$${\displaystyle \lfloor x_{k}\rfloor \geq \lfloor {\sqrt {n}}\rfloor }$${\displaystyle x_{k}\geq \lfloor x_{k}\rfloor >x_{k+1}\geq \lfloor x_{k+1}\rfloor }$${\displaystyle x_{s}}$${\displaystyle x_{s}}$${\displaystyle \lfloor {\sqrt {n}}\rfloor \leq \lfloor x_{s}\rfloor \leq {\sqrt {n}}}$${\displaystyle \lfloor x_{s+1}\rfloor \geq \lfloor x_{s}\rfloor }$${\displaystyle \lfloor x_{k+1}\rfloor \geq \lfloor x_{k}\rfloor }$

しかし、は必ずしも上記の反復式の不動点ではありません。 実際、が完全な平方でない場合に限り、 が不動点であることが示されます。 が完全な平方である場合、シーケンスは収束する代わりに 、 との間の周期 2 サイクルで終了します。${\displaystyle \lfloor {\sqrt {n}}\rfloor }$${\displaystyle \lfloor {\sqrt {n}}\rfloor }$${\displaystyle n+1}$${\displaystyle n+1}$${\displaystyle \lfloor {\sqrt {n}}\rfloor }$${\displaystyle \lfloor {\sqrt {n}}\rfloor +1}$

// 整数の平方根
符号なし整数int_sqrt (符号なし整数s )    
{
	// ゼロはゼロになる
    // 1つは1つを生み出す
	s < = 1の場合    
		sを返します。 

    // 初期見積もり（高すぎるはずです）
	符号なし整数x0 = s / 2 ;      

	// アップデート
	符号なし整数x1 = ( x0 + s / x0 ) / 2 ;          

	while ( x1 < x0 ) // 境界チェック   	
	{
		x0 = x1 ;  
		x1 = ( x0 + s / x0 ) / 2 ;        
	}		
	x0を返します。 
}

たとえば、上記のアルゴリズムを使用して 2000000の整数平方根を計算すると、次のシーケンス が得られます。合計で 13 回の反復ステップが必要です。Heron 法は解に 2 乗的に収束しますが、最初は反復ごとに 1 ビット未満の精度しか得られません。つまり、初期推定値の選択がアルゴリズムのパフォーマンスにとって重要であるということです。 $${\displaystyle {\begin{aligned}&1000000\rightarrow 500001\rightarrow 250002\rightarrow 125004\rightarrow 62509\rightarrow 31270\rightarrow 15666\rightarrow 7896\\&\rightarrow 4074\rightarrow 2282\rightarrow 1579\rightarrow 1422\rightarrow 1414\rightarrow 1414\end{aligned}}}$$

2 進対数の整数部分またはビット長の高速計算が利用可能な場合 (たとえば、C++20std::bit_widthなど)、より大きい2 の 最小の累乗である から開始する方がよいでしょう 。2000000の整数平方根の例では、、となり、結果のシーケンスは となります。 この場合、必要な反復ステップは 4 つだけです。 $${\displaystyle x_{0}=2^{\lfloor (\log _{2}n)/2\rfloor +1},}$$${\displaystyle {\sqrt {n}}}$${\displaystyle \lfloor \log _{2}n\rfloor =20}$${\displaystyle x_{0}=2^{11}=2048}$$${\displaystyle 2048\rightarrow 1512\rightarrow 1417\rightarrow 1414\rightarrow 1414.}$$

従来の紙とペンを使った平方根の計算アルゴリズムは、上位の桁から下位の桁へと計算を進め、各新しい桁で平方根が得られる最大の数字を選択するというものです。1 の位の後で停止すると、計算結果は整数の平方根になります。 ${\displaystyle {\sqrt {n}}}$${\displaystyle \leq n}$

2 を基数として操作する場合、桁の選択は 0 (「小さい候補」) と 1 (「大きい候補」) の間で簡略化され、桁の操作はバイナリ シフト演算で表現できます。 は*乗算、<<は左シフト、 は>>論理右シフトであるため、任意の自然数の整数平方根を求める再帰アルゴリズムは次のようになります。

定義 integer_sqrt ( n :  int )  ->  int :
    assert  n  >=  0 、 「sqrt は非負の入力に対してのみ機能します」
     n  <  2の場合:
        戻る n

    # 再帰呼び出し:
    小さい方の値 = 整数の平方根( n  >>  2 )  <<  1
    大きいカンデラ = 小さいカンデラ +  1
     Large_cand  *  large_cand  >  nの場合:
         小さいcandを返す
    それ以外：
         large_candを返す


# 同等:
定義 integer_sqrt_iter ( n :  int )  ->  int :
    assert  n  >=  0 、 「sqrt は非負の入力に対してのみ機能します」
     n  <  2の場合:
        戻る n

    # シフト量を検索します。[[find first set]]も参照してください。
    # シフト = ceil(log2(n) * 0.5) * 2 = ceil(ffs(n) * 0.5) * 2
    シフト =  2
    while  ( n  >>  shift )  !=  0の場合:
        シフト +=  2

    # ビット設定ループを展開します。
    結果 =  0
     シフト >=  0の場合:
        結果 = 結果 <<  1
        大きい_cand  =  (
            結果 +  1
        )   # 最後のビットは常に 0 なので、結果 ^ 1 (xor) と同じです。
         large_cand  *  large_cand  <=  n  >> シフトの場合:
            結果 =  large_cand
        シフト -=  2

     結果を返す

従来のペンと紙による桁ごとのアルゴリズムの表現には、上記のコードには存在しないさまざまな最適化が含まれています。特に、前の桁の2乗を事前に減算して一般的な乗算ステップを不要にするトリックがあります。例については、平方根の計算方法§ 2進数システム（基数2）を参照してください。 ^[1]

カラツバ平方根アルゴリズムは、バーニケル・ツィーグラーのカラツバ除算とカラツバ乗算を使用した場合、「50〜1,000,000桁」の大きな整数に対する高速アルゴリズムです。^[2]

64 ビットの符号なし整数のアルゴリズムの例を以下に示します。アルゴリズム:

1. u64_isqrt内の入力を正規化します。
2. 正規化された入力を必要とするu64_normalized_isqrt_remを呼び出します。
3. 正規化された入力ビットの最上位半分を使用してu32_normalized_isqrt_rem を呼び出します。最上位ビットは同じままなので、すでに正規化されています。
4. ビット数が十分に小さいときに、より高速なアルゴリズムが見つかるまで再帰的に続行します。
5. u64_normalized_isqrt_rem は、返された整数の平方根と剰余を取り、指定された正規化されたu64の正しい結果を生成します。
6. 次に、 u64_isqrt は結果を非正規化します。

/// `u64` で Karatsuba 平方根を実行します。
pub fn u64_isqrt ( mut n : u64 ) -> u64 {     
    n <= u32の場合:: MAXをu64として{      
        // `n` が `u32` に収まる場合は、`u32` 関数で処理します。
        u32_isqrt ( nをu32として)をu64として返します。     
    }それ以外{  
        // 正規化シフトはカラツバ平方根を満たす
        // アルゴリズムの前提条件「a₃ ≥ b/4」ここでa₃は最も
        // `n` のビットの重要な4分の1、b は
        // その 4 分の 1 のビットで表すことができる値。
        //
        // b/4は2番目に重要なものを除いてすべて0になります
        // 2進数のビット（010...0）。a₃は少なくともb/4でなければならないので、a₃の
        // 最上位ビットまたはその隣のビットは1でなければなりません。a₃の
        // 最上位ビットは`n`の最上位ビットであり、
        // `n` にも同じことが適用されます。
        //
        // 偶数ビットシフトする理由は、
        // 偶数ビットは平方根をシフトしたものを生成します
        // 正規化シフトの半分だけ左に移動します:
        //
        // 平方根(n << (2 * p))
        // sqrt(2.pow(2 * p) * n)
        // sqrt(2.pow(2 * p)) * sqrt(n)
        // 2.pow(p) * sqrt(n)
        // 平方根(n) << p
        //
        // 奇数ビットシフトすると醜い sqrt(2) が残る
        // 乗算されます。
        const EVEN_MAKING_BITMASK : u32 = ! 1 ;   
        normalization_shift = n . leading_zeros () & EVEN_MAKING_BITMASKとします。     
        n <<=正規化シフト;  

        sと_をu64_normalized_isqrt_rem ( n )とします。    

        非正規化シフトを正規化シフト/ 2とします。     
        sを返します>> denormalization_shift ;   
    }
}

/// 正規化された `u64` に対してカラツバ平方根を実行し、平方数を返す。
/// 根と剰余。
fn  u64_normalized_isqrt_rem ( n : u64 )  -> ( u64 , u64 ) {  
    定数HALF_BITS : u32 = u64 :: BITS >> 1 ;     
    const QUARTER_BITS : u32 = u64 :: BITS >> 2 ;     
    定数LOWER_HALF_1_BITS : u64 = ( 1 << HALF_BITS ) - 1 ;       

    デバッグアサート! (
        n . leading_zeros () <= 1 、  
        「入力は正規化されていません: {n} の先頭には 0 または 1 ではなく、{} 個のゼロ ビットがあります。」
        n .先頭ゼロ()
    );

    hi = ( n >> HALF_BITS )をu32とします。       
    lo = n & LOWER_HALF_1_BITSとします。     

    s_primeとr_primeをu32_normalized_isqrt_rem ( hi )とします。    

    分子を( ( r_prime as u64 ) << QUARTER_BITS ) | ( lo >> QUARTER_BITS )とします。           
    分母= ( s_prime as u64 ) << 1とします。       

    q =分子/分母とします。     
    u =分子%分母とします。     

    mut s = ( s_prime << QUARTER_BITS )をu64 + qとします。          
    mut r = ( u << QUARTER_BITS ) | ( lo & (( 1 << QUARTER_BITS ) - 1 ));とします。              
    q_squared = q * qとします。     
    r < q_squared {の場合    
        r += 2 * s - 1 ;      
        s -= 1 ;  
    }
    r -= q_squared ;  

    戻り値( s , r );  
}

一部のプログラミング言語では、一般的なケースに加えて、整数の平方根計算専用の明示的な操作が用意されていたり、この目的のためにライブラリによって拡張されたりします。

• 平方根を計算する方法