3
$\begingroup$

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.

  • 0
    Rather than "up to 16 bits" being not quite right, it's totally wrong. Typical IEEE754 64-bit [double-precision](http://en.wikipedia.org/wiki/Double_precision_floating-point_format) numbers use effectively 53 bits for mantissa, giving ~ 16 decimal digits accuracy. (15.95... = 53 * ln(2)/ln(10).)2011-10-10

1 Answers 1

7

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.

  • 0
    @eudoxos: Thanks for the idea.2011-08-07