12

平方根 N の連分数を生成するためにこのコードを書きまし
たが、N = 139
の場合は失敗します{11,1,3,1,3,7,1,1,2,11,2,1,1,7,3,1,3,1,22}
。 22 に達すると 12 になります。

誰かがこれで私を助けることができますか?

vector <int> f;
int B;double A;
A = sqrt(N*1.0);
B = floor(A);
f.push_back(B);                 
while (B != 2 * f[0])) {
    A = 1.0 / (A - B);
    B =floor(A);                            
    f.push_back(B);     
}
f.push_back(B);
4

7 に答える 7

24

根本的な問題は、非平方の平方根を浮動小数点数として正確に表現できないことです。

ξが正確な値であり近似値である場合x(これはまだ非常に適切であると考えられているため、特にfloor(ξ) = a = floor(x)まだ保持されます)、連分数アルゴリズムの次のステップの後の差は次のとおりです。

ξ' - x' = 1/(ξ - a) - 1/(x - a) = (x - ξ) / ((ξ - a)*(x - a)) ≈ (x - ξ) / (ξ - a)^2

したがって、各ステップで、近似値と実際の値の差の絶対値が増加することがわかります0 < ξ - a < 1。ξ - a大きな部分商が発生する (が 0 に近い)たびに、差は大きな係数で増加します。差 (の絶対値) が 1 以上になると、次に計算される部分商が間違っていることが保証されますが、おそらく最初に間違った部分商が先に発生します。

チャールズ は、正しい数字を使用した元の近似を使用して、連分数の部分商について計算できるという近似について言及しました。これは良い経験則ですが、これまで見てきたように、部分商が大きいと精度が高くなるため、取得できる部分商の数が減り、部分商を誤ってしまうことがあります。nn

の場合は√139、いくつかの大きな部分商を伴う比較的長い期間のケースであるため、最初に誤って計算された部分商が期間が完了する前に表示されることは驚くべきことではありません (それがより早く発生しないことにむしろ驚いています)。

浮動小数点演算を使用すると、それを防ぐ方法はありません。

しかし、二次スルドの場合、整数演算のみを使用することでその問題を回避できます。の連分数展開を計算したいとします。

ξ = (√D + P) / Q

ここでQは分割され、完全な正方形D - P²でD > 1はありません (分割可能条件が満たされない場合は、 、 および に置き換えることができDますD*Q²。ケースはでPあり、自明に満たされます)。完全な商を次のように書きます。P*QQQ²P = 0, Q = 1

ξ_k = (√D + P_k) / Q_k (with ξ_0 = ξ, P_0 = P, Q_0 = Q)

と は部分商を示しa_kます。それで

ξ_k - a_k = (√D - (a_k*Q_k - P_k)) / Q_k

そしてP_{k+1} = a_k*Q_k - P_k、

ξ_{k+1} = 1/(ξ_k - a_k) = Q_k / (√D - P_{k+1}) = (√D + P_{k+1}) / [(D - P_{k+1}^2) / Q_k],

so Q_{k+1} = (D - P_{k+1}^2) / Q_k—P_{k+1}^2 - P_k^2は の倍数であるQ_kため、帰納法によりQ_{k+1}は整数であり、 をQ_{k+1}割りD - P_{k+1}^2ます。

実数の連続分数展開は、 が 2 次スルドであるξ場合にのみ周期的ξであり、上記のアルゴリズムで最初のペア(P_k, Q_k)が繰り返されるときに周期が完了します。純粋な平方根の場合は特に単純です。ピリオドはQ_k = 1a の最初のときに完了しk > 0、P_k, Q_k常に非負です。

ではR = floor(√D)、部分商は次のように計算できます。

a_k = floor((R + P_k) / Q_k)

したがって、上記のアルゴリズムのコードは次のようになります

std::vector<unsigned long> sqrtCF(unsigned long D) {
    // sqrt(D) may be slightly off for large D.
    // If large D are expected, a correction for R is needed.
    unsigned long R = floor(sqrt(D));
    std::vector<unsigned long> f;
    f.push_back(R);
    if (R*R == D) {
        // Oops, a square
        return f;
    }
    unsigned long a = R, P = 0, Q = 1;
    do {
        P = a*Q - P;
        Q = (D - P*P)/Q;
        a = (R + P)/Q;
        f.push_back(a);
    }while(Q != 1);
    return f;
}

√7981これは、182 の周期長を持つ(eg) の連続分数を簡単に計算します。

于 2012-08-30T01:09:02.607 に答える
5

犯人はそうではありませんfloor。原因は計算A= 1.0 / (A - B);より深く掘り下げると、原因はコンピュータが実数を表すために使用する IEEE 浮動小数点メカニズムです。減算と加算は精度を失います。アルゴリズムが繰り返し行っているように繰り返し減算すると、精度が失われます。

連分数項 {11,1,3,1,3,7,1,1,2,11,2} を計算するまでに、A の IEEE 浮動小数点値は、 15 人か 16 人が予想するだろう。{11,1,3,1,3,7,1,1,2,11,2,1,1,7,3,1,3,1} に到達するまでに、A の値はまったくのゴミです. それはすべての精度を失いました。

于 2012-08-29T18:09:32.777 に答える
0

スプレッドシートであなたのアルゴを使用しましたが、12も取得しました。あなたはアルゴを間違えたに違いないと思います。私は、253の値を試しましたが、Bはその最終値に達しませんでした。

アルゴが何をすべきか、そしてそれがどのように機能するかをもう少し説明してみてください。

私はあなたのアルゴを手に入れました、そしてあなたはあなたの質問で間違いを犯したと思います、それは12であるはずです。将来の参考のために、アルゴはこのページhttp://en.wikipedia.org/wiki/Continued_fractionで見つけることができ、それは非常に傾向があります逆の値が丸めようとしている整数に非常に近い場合の小数/数値の計算上の問題。

Excelでプロトタイプを作成するとき、3.245のwikiページの例を再現できませんでした。これは、ある時点で、Floor（）が数値を4ではなく3にフロア化したため、精度をチェックするための境界チェックが必要になるためです...

この場合、おそらく最大反復回数、終了条件をチェックするための許容値を追加する必要があります（終了条件はAがBに等しいことである必要があります）

于 2012-08-29T16:51:06.780 に答える