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

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.

-
A picky remark: "up to 16 bits" is not quite right, since each number in the continued fraction can be arbitrarily large, so you can't represent each of them with just one bit. (By the way, not that it matters to your question, but your first example ends with 2,4,2,2,1,48,2 and not 2,4,2,2,1,4,8,2 as it might seem from the output.) –  Hans Lundmark Aug 7 '11 at 8:41
Rather than "up to 16 bits" being not quite right, it's totally wrong. Typical IEEE754 64-bit double-precision numbers use effectively 53 bits for mantissa, giving ~ 16 decimal digits accuracy. (15.95... = 53 * ln(2)/ln(10).) –  jwpat7 Oct 10 '11 at 20:55

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.
@Chan: "you need to first compute $\sqrt6$ to high precision." - luckily enough, that is easily done: start with your double precision approximation, and then polish with Newton-Raphson: $x_{i+1}=x_i-\frac{x_i^2-6}{2x_i}$. –  Ｊ. Ｍ. Aug 7 '11 at 9:22
@Hans Lundmark: Great thanks. Interesting point, I never thought $\sqrt{6}$ would cause the issue though. –  Chan Aug 7 '11 at 17:48