Continued fraction code error

Question

 #include <iostream>
 #include <math.h>

 int main()
 {
   double result = sqrt(4732);
   int intResult = (int)result;

   result = result - (double)intResult;
   double first = result;
   while(1)
   {   
     result = 1/result;
     intResult = (int)result;
     result = result - intResult;
     std::cout<<intResult<<std::endl;
     double absDiff = result > first ? (result-first):(first-result);
     if(absDiff < 0.000001)
       break;
   }   
   return 0;
}

I am trying to calculate continued fraction of square root of 4732. Here is the wiki description for the continued fraction. https://en.wikipedia.org/wiki/Continued_fraction

The correct answer is [68;1,3,1,3,45,1,1,2,11,15,5,34,5,15,11,2,1,1,45,3,1,3,1,136]. My code goes in an infinite loop. Here is the output from the first few iterations. [68;1,3,1,3,45,1,1,2,11,15,5,126,.. Note that the output starts diverging from here.

I have checked this code for other numbers like 9999, 13, 12, etc. and it is giving the correct answer. Can anybody point me to the problem in the code? Thanks in advance.

Update: I used float128 as well. But it does not solve the issue. It looks like I may have to use some another algorithm to avoid losing precision.

@OliverCharlesworth I am currently debugging it using gdb. But I am hoping that I have missed out something trivial. So had posted this question. — Sandeep, Aug 27 '16 at 22:17
Perhaps [Boost.Multiprecision](http://www.boost.org/doc/libs/1_61_0/libs/multiprecision/doc/html/index.html). — Igor R., Aug 27 '16 at 22:17
Possible duplicate of [Is floating point math broken?](http://stackoverflow.com/questions/588004/is-floating-point-math-broken) — Richard Critten, Aug 27 '16 at 22:30
Don't use floating-point arithmetic when you need exact resutls. You can use [this algorithm](https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Algorithm) using integers. — interjay, Aug 28 '16 at 11:38

score 1 · Accepted Answer · answered Aug 28 '16 at 21:36

Thanks to interjay for pointing me to how to solve the problem. Here is the code that gave the correct results for me. The key here is not to use floating point arithmetic as it would propagate error even if you use 128 bits of precision.

#include <iostream> 
#include <math.h>

int main()
{
  int num = 4732;
  double result = sqrt(num);
  int intResult = (int)result;

  int a = intResult;
  int b = 1;
  int first_a = a;
  int first_b = b;

  do  
  {   
    double temp = (double)b/ (result - (double)a);
    int I = (int)temp;
    int new_a = (I*((num - (a*a))/b)) - a;
    int new_b = (num - (a*a))/b;
    a = new_a;
    b = new_b;
    std::cout<<I<<std::endl;
  }   
  while((first_a != a) || (first_b != b));

  return 0;
}

Continued fraction code error

1 Answers1