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.