Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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.

share|improve this question
    
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
add comment

1 Answer 1

up vote 6 down vote accepted

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.

share|improve this answer
1  
@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}$. –  J. M. Aug 7 '11 at 9:22
2  
@J.M. Or use function for computing arbitrary-precision square root from MPIR itself. –  eudoxos Aug 7 '11 at 13:54
1  
@eudoxos: ...which does Newton-Raphson internally. :D –  J. M. Aug 7 '11 at 13:56
    
@J.M: you've got the point mathematically. As programmer however, I would prefer to concentrate on the problem in question. –  eudoxos Aug 7 '11 at 14:00
    
@Hans Lundmark: Great thanks. Interesting point, I never thought $\sqrt{6}$ would cause the issue though. –  Chan Aug 7 '11 at 17:48
show 1 more comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.