Generating continued fractions for square roots
Asked Answered
S

8

12

I wrote this code for generating Continued Fraction of a square root N.
But it fails when N = 139.
The output should be {11,1,3,1,3,7,1,1,2,11,2,1,1,7,3,1,3,1,22}
Whilst my code gives me a sequence of 394 terms... of which the first few terms are correct but when it reaches 22 it gives 12!

Can somebody help me with this?

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);
Sonority answered 29/8, 2012 at 16:43 Comment(8)
What are the types of A and B ?Smatter
My answer here might be helpful. It generates and prints out the continued fraction for an arbitrary double.Veridical
The root problem that you are running into is simply a bad case of numerical instability. Generating continued fractions is inherently numerically unstable. The code in that answer I linked to does a work-around.Veridical
@Mysticial: this is exactly what I found also, it does not work very well in practice without modifications because the flooring can cause trouble and checking that the flooring has no effect should be done via tolerance.Zirconium
@Zirconium Is there any test case where it actually fails? Note that my version regenerates the residual at each step by subtracting out the current convergent. So the buildup of round-off is restricted to only two operations. The only case where that floor would fail is if a 2-ULP error pushes it across an integer boundary. (which in that case the next convergent would be huge anyway)Veridical
@Mysticial: I am at work so I did not put it in C++, but in Excel the example of the wiki page en.wikipedia.org/wiki/Continued_fraction for 3.245 fails using floor(), I know that Excel decimal numbers are not doubles though, but my take on it is that you are exactly spot on, the case where it would break is that if the floor returns for example 3 instead of 4 and this is what happens when you do it in Excel with 3.245 (although it may pass in C++ due to the double having a better precision)Zirconium
Oops, I realized that I didn't actually link to my answer. It's here.Veridical
The code in my answer spits out {3, 4, 12, 3, 1} for 3.245. So it looks like it's working correctly.Veridical
S
26

The root problem is that you cannot exactly represent the square root of a non-square as a floating-point number.

If ξ is the exact value and x the approximation (which is supposed to be still quite good, so that in particular floor(ξ) = a = floor(x) still holds), then the difference after the next step of the continued fraction algorithm is

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

Thus we see that in each step the absolute value of the difference between the approximation and the real value increases, since 0 < ξ - a < 1. Every time a large partial quotient occurs (ξ - a is close to 0), the difference increases by a large factor. Once (the absolute value of) the difference is 1 or greater, the next computed partial quotient is guaranteed to be wrong, but very probably the first wrong partial quotient occurs earlier.

Charles mentioned the approximation that with an original approximation with n correct digits, you can compute about n partial quotients of the continued fraction. That is a good rule of thumb, but as we saw, any large partial quotients cost more precision and thus reduce the number of obtainable partial quotients, and occasionally you get wrong partial quotients much earlier.

The case of √139 is one with a relatively long period with a couple of large partial quotients, so it's not surprising that the first wrongly computed partial quotient appears before the period is completed (I'm rather surprised that it doesn't occur earlier).

Using floating-point arithmetic, there's no way to prevent that.

But for the case of quadratic surds, we can avoid that problem by using integer arithmetic only. Say you want to compute the continued fraction expansion of

ξ = (√D + P) / Q

where Q divides D - P² and D > 1 is not a perfect square (if the divisibility condition is not satisfied, you can replace D with D*Q², P with P*Q and Q with ; your case is P = 0, Q = 1, where it is trivially satisfied). Write the complete quotients as

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

and denote the partial quotients a_k. Then

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

and, with 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 — since P_{k+1}^2 - P_k^2 is a multiple of Q_k, by induction Q_{k+1} is an integer and Q_{k+1} divides D - P_{k+1}^2.

The continued fraction expansion of a real number ξ is periodic if and only if ξ is a quadratic surd, and the period is completed when in the above algorithm, the first pair (P_k, Q_k) repeats. The case of pure square roots is particularly simple, the period is completed when first Q_k = 1 for a k > 0, and P_k, Q_k are always nonnegative.

With R = floor(√D), the partial quotients can be calculated as

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

so the code for the above algorithm becomes

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;
}

which easily calculates the continued fraction of (e.g.) √7981 with a period length of 182.

Shakiashaking answered 30/8, 2012 at 1:9 Comment(4)
Thanks a lot .That was a really big help . But can you provide me with a link for paper or something that describes the CF Algorithms with intensive explanationSonority
They are described (with varying depth) in more or less every book whose title is sufficiently similar to "Introduction to Number Theory". From my own experience, I can recommend the classic Hardy/Wright, the book by Hua Loo Keng, and - with a caveat - the book by Don Redmond. (The caveat is that at least in the edition I read, Redmond's book had a lot of typesetting errors, where e.g. 13·27 appeared as 1327, 10 squared as 102 etc, which is often irritating. Apart from that, it's a very good book, and treats CFs in more breadth than the other two.)Shakiashaking
Great solution. Can you please explain why a = (R + P)/Q though? I am not quite getting it.Trela
@ArulxZ The value of the partial quotient is floor((sqrt(D) + P)/Q). Now it's a fact that for a positive integer Q, and any positive real x one has floor(x/Q) == floor(floor(x)/Q), so we have floor((sqrt(D)+P)/Q) == floor(floor(sqrt(D)+P)/Q), but floor(sqrt(D)+P) == floor(sqrt(D)) + P since P is an integer, and with R == floor(sqrt(D)), that becomes floor((R+P)/Q). Since integer division in C++ truncates, and everything is positive, the expression (R+P)/Q in C++ code computes floor((R+P)/Q), as desired.Shakiashaking
F
5

The culprit isn't floor. The culprit is the calculation A= 1.0 / (A - B); Digging deeper, the culprit is the IEEE floating point mechanism your computer uses to represent real numbers. Subtraction and addition lose precision. Repeatedly subtracting as your algorithm is doing repeatedly loses precision.

By the time you have calculated the continued fraction terms {11,1,3,1,3,7,1,1,2,11,2}, your IEEE floating point value of A is good to only six places rather than the fifteen or sixteen one would expect. By the time you get to {11,1,3,1,3,7,1,1,2,11,2,1,1,7,3,1,3,1} your value of A is pure garbage. It has lost all precision.

Footle answered 29/8, 2012 at 18:9 Comment(1)
To a good first approximation (a coincidence that works only in base 10) about n continued fraction terms can be extracted from a number with n digits of precision. If you try to take more you'll get garbage, as you pointed out.Sluiter
T
2

The sqrt function in math is not precise. You can use sympy instead with arbitrarily high precision. Here is a very easy code to calculate continued fractions for any square root or number included in sympy:

from __future__ import division #only needed when working in Python 2.x
import sympy as sp

p=sp.N(sp.sqrt(139), 5000)

n=2000
x=range(n+1)
a=range(n)
x[0]=p

for i in xrange(n):
    a[i] = int(x[i])
    x[i+1]=1/(x[i]-a[i])
    print a[i],

I have set the precision of your number to 5000 and then calculated 2000 continued fraction coefficients in this example code.

Tarpeia answered 29/4, 2017 at 15:14 Comment(1)
Yes, this algorithm works in the general case, but as Daniel Fischer's answer shows, it's possible to determine the exact periodic continued fraction for quadratic numbers without resorting to arbitrary precision arithmetic. Of course, if you want to evaluate that continued fraction to high precision then you do need something better than standard double precision. BTW, there's probably not much point posting Python code to a question tagged C++, OTOH your code doesn't use any Python "tricks" so it should be fairly clear to anyone reading this page.Variety
U
2

In case anyone is trying to solve this is in a language without integers, here is the code from the accepted answer adapted for JavaScript.

Note two ~~ (floor operators) have been added.

export const squareRootContinuedFraction = D =>{
    let R = ~~Math.sqrt(D);
    let f = [];
    f.push(R);
    if (R*R === D) {
        return f;
    }
    let a = R, P = 0, Q = 1;
    do {
        P = a*Q - P;
        Q = ~~((D - P *P)/Q);
        a = ~~((R + P)/Q);
        f.push(a);
    } while (Q != 1);
    return f;
};
Unimpeachable answered 26/12, 2017 at 22:24 Comment(0)
Z
0

I used you algo in a spreadsheet and I get 12 as well, I think you must have made a mistake in your algo, I tried for 253 values, and B did not reach the final value of it.

Can you try to explain a bit more what the algo should do and how it would work ?

I think I got your algo and you made a mistake in your question, it should be 12. For future reference, the algo can be found at this page http://en.wikipedia.org/wiki/Continued_fraction and it is very prone to issues with decimal/numerical computational issues if the inverse value is very close to the integer you are trying to round at.

When doing the prototype under Excel, I could not reproduce the example of the wiki page for 3.245, because at some point Floor() floored the number to 3 instead of 4, so some boundary checking to check for accuracy is required ...

In this case you probably want to add a maximum number of iteration, a tolerance for checking the exit condition (the exit condition should be that A is equal to B btw)

Zirconium answered 29/8, 2012 at 16:51 Comment(0)
I
0

Your code isn't calculating the square root of n. It attempts to calculate the continued fraction of an already calculated √n. I mean it's OK but, had it been correct, your approach is more suitable for general decimal to rational convertion. However for a sqrt function for the regular (simple) continued fractions (where all numerators are ones) the algorithm is slightly different.

However the problem isn't finished. Yes as a rule, the CF coefficients for √n is in the form of a repeating palindrome which ends with the double of the first nonzero coefficient. Such as √31 =[5;1,1,3,5,3,1,1,10,1,1,3,5,3,1,1,10..]. Now there is no simple way to quess the length of the palindrome for every given n. There are some known patterns however they are far from defining a generalized pattern for all n. So stoping the iteration at the end of the first palindrome is a very undeterministic approach. Imagine

          __
√226 =[15;30]

while

          ____________________________________________________
√244 =[15;1,1,1,1,1,2,1,5,1,1,9,1,6,1,9,1,1,5,1,2,1,1,1,1,1,30]

If you decide to stop the iteration at 2*f[0] most of the time you get either a bad approximation like in √226 or an over-calculated one like in √244 case. Besides, once n grows, chasing the end of the palindrome becomes more meaningless since you will never need such a precision.

            ___________________________________________________________________________________________________________________________________________________________________________
√7114 = [84;2,1,9,3,1,10,2,23,1,1,1,1,1,2,1,27,2,1,1,3,1,2,1,1,1,16,4,3,1,3,2,1,6,18,1,1,2,6,11,11,6,2,1,1,18,6,1,2,3,1,3,4,16,1,1,1,2,1,3,1,1,2,27,1,2,1,1,1,1,1,23,2,10,1,3,9,1,2,168]

In this case it would be reasonable to stop the iteration once the necessary precision is obtained. As I have mentioned at the beginning, there are two approaches.

  1. The general Decimal to Rational algorithm obtains simple continued fractions from any decimal number. This will give a CF resolving to exactly to that decimal without any floating point errors. For a detailed explanation of this algorithm you may check a previous answer of mine.
  2. There happens to be a more direct algorithm for √n which is basically the same as 1 but tuned for square roots. In this case you don't provide √n but n.

The idea is as follows. We have to define a general form for the input which involves a square root value at the numerator. Then we try to reach to the same expression in the continued fraction part to be able to iterate.

Let our input be of the form

q + √n
______
   p

for simple square root operation we can assume q is 0 and p is 1. If we can establish this form at the next stage then we can easily iterate on.

Starting with the initial stage where q = 0, p = 1, m is the integer part of √n and 1/x is the floating part, our target is to bring x into (q + √n) / p form;

          1            1             1       (√n + m)        √n + m
√n = m + ___ ⇒ x = _______ ⇒ x = ________ . ________ ⇒ x = ________
          x         √n - m        (√n - m)   (√n + m)        n - m^2

Now √n is at the numerator and we have the form of;

    √n + q
x = ______
       p

where q = m and p = n - m^2. Right at this point you can calculate x and the next m by flooring x. The generalized form of the algorithm becomes;

    √n + q        1                p          p(√n - (q - pm))   p(√n + (pm - q))
x = ______ = m + ___ ⇒ x' = ______________ = _________________ = ________________
       p          x'         √n + (q - pm)     n - (q - pm)^2     n - (q - pm)^2

at this point p is divisible by n - (q - pm)^2. This is now stable and we can extend it as long as we want. Let's make new assignments for q and p as;

q' = pm-q;
p' = (n - q'^2)/p;

     √n + q'
x' = ______
        p'

m' = Math.floor(x')

Note that when p' becomes 1 (n - q'^2 = p) we are at the end of the palindrome. However in order to decide where to stop I am using the same mechanism as described in my toRational algorithm as linked in the alternative 1 above. It basically stops once the JS floating point resolution is reached. The JavaScript code is as follows;

function rationalSqrt(n){
  var nr = Math.sqrt(n),
      m  = Math.floor(nr),
      p  = n-m**2,
      q  = m,
      cs = [m],
      n0 = 1,
      d0 = 0,
      n1 = m,
      d1 = 1,
      n2 = 0,
      d2 = 1;
  if (nr === m) return {n:m,d:1,cs};
  while (Math.abs(nr-n2/d2) > Number.EPSILON){
    m = Math.floor((nr+q)/p);
    q = m*p-q;
    p = (n-q**2)/p;
    cs.push(m);
    n2 = m*n1+n0;
    d2 = m*d1+d0;
    n0 = n1;
    d0 = d1;
    n1 = n2;
    d1 = d2;
  }
  return {n:n2,d:d2,cs};
}

These two algorithms differ slightly.

  1. ~60% of the time they result same fractions.
  2. ~27% of the time toRational gives smaller numerators and denominators within the JS floating point resolution. However surely it's not a better approximation had we have one more digit in JS.
  3. ~13% of the time rationalSqrt (this one) gives smaller numerators and denominators within the JS floating point resolution.
  4. rationalSqrt will result exact coefficients as one would expect from a square root but will be truncated once the resolution is enough.
  5. toRational gives expected coefficients of couse but the last one could be totally unrelated to what you would expect from the square root series.

One such example is;

rationalSqrt(511); //returns
{ n : 882184734
, d : 39025555
, cs: [22,1,1,1,1,6,1,14,4,1,21,1,4,14,1,6,1]
}

while

toRational(Math.sqrt(511));
{ n : 1215746799
, d : 53781472
, cs: [22,1,1,1,1,6,1,14,4,1,21,1,4,14,1,10]
}

Further thoughts:

  • Consider that we are given the RCF coefficients. Can we reverse the rationalSqrt alogrithm to obtain it's(q + √n) / p form? This might come as an interesting task.
  • Practiacally limiting the CF with the JS decimal resolution would yield a similar precision in CF arithmetics and you would be wasting your time. In order to achieve higher square root resolution I advise to replace while (Math.abs(nr-n2/d2) > Number.EPSILON) condition with while (cs.length < 30) if 30 CF coefficients of resolution would be sufficient for your application. That figure is up to you.
Ineffective answered 1/3, 2022 at 10:12 Comment(0)
Y
0

for what it's worth, this is what it looks like for the 1st 6 levels of continued fraction estimation for sqrt(2) :

return \
    (_ < (__ = --__-_) ? \
        (_--+_<__ ? _/(_/(_/(_/(_/(_/(_/ ++_ +_)+_)+_)+_)+_)+_)  \
        :  _+_<__ ? _/(_/(_/(_/(_/(_/    ++_    +_)+_)+_)+_)+_)  \
                  : _/(_/(_/(_/(_/       ++_       +_)+_)+_)+_)) \
        : __>!--_ ? _/(_/(_/(_/          ++_          +_)+_)+_)  \
        :    !__  ? _/(_/(_/             ++_             +_)+_)  \
                  : _/(_/                ++_                +_)) + --_

Here, the variable _ always represent a 1 to the left of the center column, and always a 2 from center column to right edge. It's easier to observe the high-level structure this way without all the numbers and letters getting in the way.

Note that at the bottom right (")) + --_"), the addition of 1 is applicable to all precision depths, since the core fraction here is only estimating the 0.41421356... part.

Ytterbite answered 11/6 at 15:9 Comment(0)
P
-1

I used a Surd Storage type for infinite precision of the square root of n.

(b * \sqrt(n) + d)/c

=

(b * c * sqrt(n) - c * d + a_i * c^2) / (b^2 * n - d^2 - (a_i * c)^2 + 2* a_i * c * d)

The floor value of the sqrt(n) is only used once. After which the remaining iterations are stored as a surd type. This avoids the rounding errors seen in other algorithms and infinite (memory constrained) resolution can be achieved.

a_0 = floor value of sqrt (n)

a_i = (b_i * a_0 + d_i) / c_i

b_i+1 = b_i * c

c_i+1 = (b_i)^2 * n - (d_i)^2 - (a_i * c_i)^2 + 2 * a_i * c_i * d_i

d_i+1 = a_i * (c_i)^2 - c_i * d_i

g = gcd(b_i+1 , c_i+1 , d_i+1)

b_i+1 = b_i+1 / g

c_i+1 = c_i+1 / g

d_i+1 = d_i+1 / g

a_i+1 = (b_i+1 * x + d_i+1) / c_i+1

Then for i=0 to i=Maximum_terms a continued fraction is produced beginning with [a_0;a_1,a_2 ... ,2*a_0]

I terminate the fraction when the a_i term equals 2 times the a_0. This is the point at which the sequence repeats.

The mathematics was done by Electro World and a very nice video on the mathematics can be found here https://youtu.be/GFJsU9QsytM

The source code written in Java with BigInteger is provided below. Hope you like it.

A boolean true is returned if a repeat series is found and false if a repeat sequence is not found for the desired precision.

The precision can be easily modified to suit as per Maximum_terms.

The square root 139 [11;1,3,1,3,7,1,1,2,11,2,1,1,7,3,1,3,1,22] repeat length 18

The square root of 15 [3;1,6] repeat length 2

The square root of 2501 [50;100] repeat length 1

The square root of 10807 [103;1,22,9,2,2,5,4,1,1,1,6,15,1,5,2,1,3,6,34,2,34,6,3,1,2,5,1,15,6,1,1,1,4,5,2,2,9,22,1,206] repeat length 40

A possible two fold speed up would be to look at the palindromic nature of series. In this case 34,2,34. Only half the sequence would need to be determined.

    public static Boolean SquareRootConFrac(BigInteger N) {
BigInteger A,B=BigInteger.ONE,C=B,D=BigInteger.ZERO;
BigInteger A0=N.sqrt(),Bi=B,Ci=C,Di=D,G;
BigInteger TwoA0 = BigInteger.TWO.multiply(A0);
int Frac_Length=0, Maximum_terms=10000; //Precision 10000 terms
String str="";
Boolean Repeat=false, Success=false, Initial_BCD=true;

while(!Repeat) {
    Frac_Length++;                         Success=!(Frac_Length==Maximum_terms);
    A=((B.multiply(A0)).add(D)).divide(C); Repeat=A.equals(TwoA0)||!Success;

    Bi=B.multiply(C);
    Ci=(B.multiply(B).multiply(N)).subtract(D.multiply(D)).subtract(A.multiply(A).multiply(C).multiply(C)).add(BigInteger.TWO.multiply(A).multiply(C).multiply(D));
    Di=(A.multiply(C).multiply(C)).subtract(C.multiply(D));
    G=Bi.gcd(Ci).gcd(Di);
    B=Bi.divide(G);C=Ci.divide(G);D=Di.divide(G);
    
    if(Initial_BCD) {str="["+A+";";System.out.print(str);Initial_BCD=false;}
    else            {str=""+A;System.out.print(str);if(!Repeat){str=",";System.out.print(str);}}
}
str="]";System.out.println(str);
str="repeat length ";System.out.print(str);
if(Success) {str=""+(Frac_Length-1);System.out.println(str);}
else        {str="not found";System.out.println(str);}
return Success;
}
Prospectus answered 5/7, 2021 at 9:43 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.