How to express an irrational as a continued fraction in computer with high precision?

Question

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.

Answer 1

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.