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.