Il problema di root è che non è possibile rappresentare esattamente la radice quadrata di un non quadrato come un numero in virgola mobile.
Se ξ è il valore esatto e x l'approssimazione (che dovrebbe essere ancora abbastanza buona, in modo che in particolare floor(ξ) = a = floor(x) detiene ancora), quindi la differenza dopo la fase successiva della frazione continua è
ξ' - x' = 1/(ξ - a) - 1/(x - a) = (x - ξ)/((ξ - a)*(x - a)) ≈ (x - ξ)/(ξ - a)^2
Così vediamo che in ogni passo il valore assoluto della differenza tra l'approssimazione e il valore reale aumenta, dal 0 < ξ - a < 1. Ogni volta che si verifica un grande quoziente parziale (ξ - a è vicino a 0), la differenza aumenta di un fattore elevato. Una volta (il valore assoluto di) la differenza è 1 o maggiore, il successivo quoziente parziale calcolato è garantito errato, ma molto probabilmente il primo quoziente parziale errato si verifica prima.
Charlesmentioned l'approssimazione che, con un'approssimazione originale con n cifre corrette, è possibile calcolare circa n quozienti parziali della frazione continua. Questa è una buona regola empirica, ma come abbiamo visto, qualsiasi grande quoziente parziale costa più precisione e quindi riduce il numero di quozienti parziali ottenibili, e di tanto in tanto si ottengono percentuali parziali errate molto prima.
Il caso di √139 è uno con un periodo relativamente lungo, con un paio di grandi quozienti parziali, quindi non è sorprendente che il primo quoziente parziale erroneamente calcolato appare prima del periodo è stata completata (sono piuttosto sorpreso che doesn' t si verificano prima).
Utilizzando l'aritmetica in virgola mobile, non c'è modo di impedirlo.
Ma per il caso di conti quadrati, possiamo evitare questo problema utilizzando solo l'aritmetica dei numeri interi. Dire che si desidera calcolare l'espansione frazione continua di
ξ = (√D + P)/Q
dove Q divide D - P² e D > 1 non è un quadrato perfetto (se la condizione divisibilità non è soddisfatta, è possibile sostituire D con D*Q², P con P*Q e Q con Q²; il tuo caso è P = 0, Q = 1, dove è banalmente soddisfatto). Scrivi i quozienti completi come
ξ_k = (√D + P_k)/Q_k (with ξ_0 = ξ, P_0 = P, Q_0 = Q)
e indicare i quozienti parziali a_k. Poi
ξ_k - a_k = (√D - (a_k*Q_k - P_k))/Q_k
e, con 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],
così Q_{k+1} = (D - P_{k+1}^2)/Q_k — dal P_{k+1}^2 - P_k^2 è un multiplo di Q_k, per induzione Q_{k+1} é un intero ed Q_{k+1} divide D - P_{k+1}^2.
L'espansione in frazione continua di un numero reale ξ è periodico se e solo se ξ è un irrazionale quadratico, e il periodo è completata quando nell'algoritmo sopra, la prima coppia (P_k, Q_k) ripetizioni. Il caso delle radici quadrate pure è particolarmente semplice, il periodo è completato quando il primo Q_k = 1 per un k > 0 e lo P_k, Q_k sono sempre non negativi.
Con R = floor(√D), i quozienti parziali può essere calcolato come
a_k = floor((R + P_k)/Q_k)
modo che il codice per l'algoritmo di cui sopra diventa
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;
}
che calcola facilmente la frazione continua di (ad esempio) √7981 con un periodo lunghezza di 182.
Quali sono i tipi di A e B? –
@PaulR A è doppio B è int –
La mia risposta potrebbe essere utile. Genera e stampa la frazione continua per un "doppio" arbitrario. – Mysticial