ニュートン・ラフソン・ディビジョン（Big Integers）

Question

私はプログラミング練習としてBigIntクラスを作っています。ベース65536の2の補数符号付き整数のベクトルを使用します（32ビット乗算がオーバーフローしないようにします。完全に動作するようにベースを増やします）。

基本的な算術演算はすべて、1つの問題でコード化されています。除算はで痛みを伴います。私は作成できました。 （それは商の各桁のバイナリディビジョンのようなものです...誰かがそれを見たくなければ投稿しません....）

私の遅いアルゴリズムではなく、ニュートン・ラフソンは、（シフトされた）逆数を見つけ、次に掛け算（およびシフト）します。私は基本の周りに私の頭を持っていると思います：式（x1 = x0（2 - x0 *除数））良い初期推測、そしていくつかの反復量の後、xは逆数に収束します。この部分は十分に簡単なようだ...しかし、大きな整数にこの式を適用しようとするとき、私はいくつかの問題に実行しています：

問題1：

を私は整数で働いているので...よく..私は分数を使うことはできません。これはxが常に発散する原因になります（x0 *除数は<でなければなりませんか？）。私の直感は、整数を（正確に）働かせることができる方程式にいくつかの修正を加えなければならないと言いますが、それが何であるかを見つけるのは本当に苦労しています。 （私の数学のスキルが私をここで打ち負かしています....）私はいくつかの同等の等式を見つける必要があると思う。の代わりにがある。d * [base^somePower]？整数で動作する（x1 = x0（2 - x0 * d））のような数式がありますか？

問題2：

私はいくつかの数字の逆数を見つけるために、ニュートンの公式を使用し、その結果は、答えはどうあるべきか、以下のほんの派閥... EXなってしまいます。 （10進数）4の逆数を見つけようとしているとき：

x0 = 0.3 
x1 = 0.24 
x2 = 0.2496 
x3 = 0.24999936 
x4 = 0.2499999999983616 
x5 = 0.24999999999999999999998926258176 

私は、ベース10の数字を表現していた場合、私は25の結果を望む（および2で右シフト製品に覚えておくこと）。 1/3などのいくつかの逆数を使用すると、正確さが十分であることがわかった後に結果を切り捨てることができます。しかし、私はどのように上記の結果から正しい相反を引き出すことができますか？

申し訳ありませんが、これがあまりにも漠然としている場合、またはあまりにも多くを求めている場合は。私はWikipediaとGoogleで見つけられるすべての研究論文を調べましたが、私は壁に頭を向けているように感じています。私は誰でも私に与えることができる任意のヘルプに感謝！

...

編集：私が予想よりもはるかに遅いですが、アルゴリズムの動作を手に入れました。私は数千桁の数字であっても、私の古いアルゴリズムに比べて実際には多くのスピードを失ってしまった...私はまだ何かが欠けている。非常に高速な乗算に問題はありません。 （私は確かにKaratsubaのアルゴリズムを使用しています）。興味がある人々のために

は、ここではニュートン・ラプソンアルゴリズムの私の現在の反復である：

bigint operator/(const bigint& lhs, const bigint& rhs) { 
    if (rhs == 0) throw overflow_error("Divide by zero exception"); 
    bigint dividend = lhs; 
    bigint divisor = rhs; 

    bool negative = 0; 
    if (dividend < 0) { 
     negative = !negative; 
     dividend.invert(); 
    } 
    if (divisor < 0) { 
     negative = !negative; 
     divisor.invert(); 
    } 

    int k = dividend.numBits() + divisor.numBits(); 
    bigint pow2 = 1; 
    pow2 <<= k + 1; 

    bigint x = dividend - divisor; 
    bigint lastx = 0; 
    bigint lastlastx = 0; 
    while (1) { 
     x = (x * (pow2 - x * divisor)) >> k; 
     if (x == lastx || x == lastlastx) break; 
     lastlastx = lastx; 
     lastx = x; 
    } 
    bigint quotient = dividend * x >> k; 
    if (dividend - (quotient * divisor) >= divisor) quotient++; 
    if (negative)quotient.invert(); 
    return quotient; 
} 

そして、ここに高速である私（本当に醜い）古いアルゴリズムです：

bigint operator/(const bigint& lhs, const bigint & rhs) { 
    if (rhs == 0) throw overflow_error("Divide by zero exception"); 
    bigint dividend = lhs; 
    bigint divisor = rhs; 

    bool negative = 0; 
    if (dividend < 0) { 
     negative = !negative; 
     dividend.invert(); 
    } 
    if (divisor < 0) { 
     negative = !negative; 
     divisor.invert(); 
    } 

    bigint remainder = 0; 
    bigint quotient = 0; 
    while (dividend.value.size() > 0) { 
     remainder.value.insert(remainder.value.begin(), dividend.value.at(dividend.value.size() - 1)); 
     remainder.value.push_back(0); 
     remainder.unPad(); 
     dividend.value.pop_back(); 

     if (divisor > remainder) { 
      quotient.value.push_back(0); 
     } else { 
      int count = 0; 
      int i = MSB; 
      bigint value = 0; 
      while (i > 0) { 
       bigint increase = divisor * i; 
       bigint next = value + increase; 
       if (next <= remainder) { 
        value = next; 
        count += i; 
       } 
       i >>= 1; 
      } 
      quotient.value.push_back(count); 
      remainder -= value; 
     } 
    } 

    for (int i = 0; i < quotient.value.size()/2; i++) { 
     int swap = quotient.value.at(i); 
     quotient.value.at(i) = quotient.value.at((quotient.value.size() - 1) - i); 
     quotient.value.at(quotient.value.size() - 1 - i) = swap; 
    } 

    if (negative)quotient.invert(); 
    quotient.unPad(); 
    return quotient; 
} 

Answer 1

まず、O(n^2)の時間分割を妥当な定数で実装することができます。したがって、素朴な乗算よりもはるかに遅くはありません。しかし、KaratsubaのようなアルゴリズムやFFTベースの乗算アルゴリズムを使用すると、実際にはニュートンラフソンを使って除算アルゴリズムをスピードアップすることができます。

xの逆数を計算するためのニュートンラフソン反復は、q[n+1]=q[n]*(2-q[n]*x)です。

floor(2^k/B)を計算するとします。ここで、Bは正の整数です。 WLOG、B≤2^k;それ以外の場合、商は0です。 x=B/2^kのNewton-Raphson反復は、q[n+1]=q[n]*(2-q[n]*B/2^k)を生成する。この種の各反復は唯一の整数乗算とビットシフトを必要とq[n+1]=q[n]*(2^(k+1)-q[n]*B) >> k

として、私たちはそれを並べ替えることができます。それは floor(2^k/B)に収束しますか？必ずしも。しかし、最悪の場合、最終的には floor(2^k/B)と ceiling(2^k/B)の間で交互になります（証明してください！）。だから、もしあなたがこの場合であるかどうかを確認するために、それほど賢明でないテストを使用して、 floor(2^k/B)を抽出することができます。 （この "それほど巧妙なテスト"は、各反復での乗算よりもずっと速くすべきですが、このことを最適化するのが良いでしょう）。

実際には、任意の正の整数A,Bに対してfloor(A/B)を計算するには、floor(2^k/B)を計算すれば十分です。 kを​​とし、floor(A/B)=A*ceiling(2^k/B) >> kを確認します。最後に、このアプローチの単純ではあるが重要な最適化は、Newton-Raphson法の初期反復において乗算を切り捨てる（すなわち、製品の上位ビットのみを計算する）ことである。その理由は、初期の反復の結果は商とは大きく異なり、不正確に実行することは重要ではないということです。時間がM(k)の2つの≤kビット整数を掛けることができると仮定して、O(M(2n))の2つの≤nビットの整数を除算することができます。

Answer 2

ニュートンRaphsonは近似アルゴリズムであり、整数の計算には適していません。あなたは、あなたが見ている種類の問題につながる丸め誤差を得るでしょう。浮動小数点数を使って問題を解決し、指定された桁数で正確な整数が得られるかどうかを調べることができます（次の段落を参照してください）。

2番目の問題点として、精度（小数点以下桁数）精度を求めて、その精度に丸めてください。問題の精度が20桁の場合は、0.25に丸めます。必要な精度の桁が安定するまで反復するだけです。一般に、コンピュータ上の不合理な数字を表すことは、しばしば不正確さをもたらす。

Answer 3

これが正しく表示されていれば、大きな改善がxの良い開始値を選ぶことです。

1/x = pow(2,log2(1/x)) 
1/x = pow(2,-log2(x)) 
1/x >= pow(2,-floor(log2(x))) 

階（LOG2（X））は、単純に、最上位ビットセットの指標であるとして、逆の最上位ビットが、しなければならどこ除数は、あなたが知っているどのように多くの桁数を知ります。