[Math] How to express an irrational as a continued fraction in computer with high precision

computer sciencecontinued-fractionsirrational-numbers

Background
I'm writing a C++ library for continued fraction using MPIR (Multiple Precision Integers and Rationals) library http://www.mpir.org/ due to the limitation of built-in double in C/C++. Indeed, it is only precise up to 16 bits.

Problem
While playing around with the library, I realized one interesting problem that I couldn't figure out how it happened.
I'm using the recursive formula for general irrational number:
$$a_k = \lfloor \alpha_k \rfloor$$
$$\alpha_{k+1} = \dfrac{1}{\alpha_k – a_k}$$

Using normal double type (built-in C/C++), my solution was:

#include <iostream>
#include <vector>
#include <cmath>

const int LENGTH = 32;

void evaluate( const double a, int length ) {
    double ak = a;
    while( length-- ) {
        std::cout << std::floor( ak );      
        ak = 1.0 / ( ak - std::floor( ak ) );
    }
    std::cout << '\n';
}

int main() {
    evaluate( std::sqrt( 6.0 ), LENGTH );
}

As expected, my output was precise up to 16 bits:

*2242424242424242*21482  
Press any key to continue . . .

Next, I tried it with MPIR library with 256 bits precision. Surprisingly, the output is different than ordinary double but not "correct" after 16th bit.

22424242424242423911124644103251171
Press any key to continue . . .

And this is my program,

#include <iostream>
#include <vector>
#include <cmath>
#include <mpir.h>

const int LENGTH = 32;

void evaluate( const double a, int length ) {
    double ak = a;
    while( length-- ) {
        std::cout << std::floor( ak );      
        ak = 1.0 / ( ak - std::floor( ak ) );
    }
    std::cout << '\n';
}

void evaluate_unlimited_bit( const double value, int length ) {
    mpf_t ak;
    mpf_t one;
    mpf_t temp;
    // initialize to 256 bits
    mpf_init2( ak, 256 );
    mpf_init2( temp, 256 );
    mpf_init2( one, 256 );

    // set value
    mpf_set_d( ak, value );
    mpf_set_d( one, 1.0 );
    mpf_set_d( temp, 0.0 );

    while( length-- ) {
        mpf_floor( temp, ak );
        std::cout << mpf_get_si( temp );
        mpf_sub( ak, ak, temp );
        mpf_div( ak, one, ak );
    }
    mpf_clear( ak );
    mpf_clear( one );
    mpf_clear( temp );
    std::cout << "\n\n";
}

int main() {
    evaluate_unlimited_bit( std::sqrt( 6.0 ), LENGTH );
    evaluate( std::sqrt( 6.0 ), LENGTH );
}

As we can see, $\sqrt{6}$ is periodic of the form $[2;\overline{2,4}]$. So I expected my result should be correct up to LENGTH bit. And I have no idea what caused this issue because I think MPIR is a very reliable library. How was the calculation off by many digits?

As a side note, I could have posted at stackoverflow.com; however, I feel MSE forum is somehow more relevant. The built-in latex typesetting is one, as well as many people has helped me over the year in this forum which make me decide to post here because I feel much easier to express my idea. If it's not appropriate, feel free to migrate it to SO. Thank you.

Best Answer

It seems to me that you're feeding your high-precision routine with the ordinary double approximation to $\sqrt{6}$, which is exactly $$x=2.44948974278317788133563226438127458095550537109375.$$ So you would be getting the continued fraction of this number $x$ with high precision, instead of that of $$y=\sqrt{6}= 2.44948974278317809819728407470589139196594748065667\ldots$$ To get what you want, you need to first compute $\sqrt{6}$ to high precision.

Related Question