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 #include #include 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 #include #include #include 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.