C - Square root function not working properly

Question

I've written a function called "mysqrt(double r)" which calculates the square root of "r". However, the result is different from the one given by the standard sqrt() function from the math.h library.

The code is the following:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
 
double mysqrt (double r);
 
int main()
{
   double r1, r2;
 
   r1 = sqrt (2.0);
   r2 = mysqrt (2.0);
 
   if (r1 == r2)
   {
      printf ("OK! Result = %.40lf\n", r1);
   }
   else
   {
      printf ("!? r1 = %.40lf, r2 = %.40lf\n", r1, r2);
   }
 
   return EXIT_SUCCESS;
}
 
double mysqrt (double r)
{
   double n1, n2;
 
   n1 = r;
   n2 = 1.0;
   while (n1 - n2 > 0)
   {
      n1 = (n1 + n2) / 2.0;
      n2 = r / n1;
   }
   return n1;
}

And the output is: !? r1 = 1.4142135623730951454746218587388284504414, r2 = 1.4142135623730949234300169337075203657150

I've tried changing the types from double to long double, but it didn't work.

Answer 1

I re-wrote the mysqrt function like this:

double mysqrt (double r)
{
   double n1, n2;
   n1 = r;
   n2 = 1.0;
   while (absolute_value(n1 - n2) > EPSILON)
   {
      n1 = (n1 + n2) / (double) 2;
      n2 = r / n1;
   }
   return n1;
}

double absolute_value (double n)
{
    if (n < 0)
    {
        n *= -1;
    }
    return n;
}

with

#define EPSILON 1e-15

By doing so the output becomes:

OK! Result = 1.4142135623730951454746218587388284504414

Answer 2

You want to approximate the irrational number sqrt(2) with a IEEE-754 binary64 representation. You have two possibilities:

approx1 double: 0 01111111111 0110101000001001111001100110011111110011101111001100
approx2 double: 0 01111111111 0110101000001001111001100110011111110011101111001101

What are those numbers? Let's convert them back to decimal:

approx1 (exact): 1.41421356237309492343001693370752036571502685546875
sqrt(2)        : 1.4142135623730950488016887242096980785696718753769480731766797379...
approx2 (exact): 1.4142135623730951454746218587388284504413604736328125

They are both pretty close. In fact they are the two closest numbers you can obtain with binary64. Which is better?

abs(sqrt(2)-approx1): 0.0000000000000001253716717905021777128546450199081980731766797379...
abs(sqrt(2)-approx2): 0.0000000000000000966729331345291303718716885982558644268233202620...

It seems that approx2 is slightly closer to the real value than approx1. For most purposes they are pretty equivalent.

I've used this website for the conversion and WolframAlpha for the comparisons.

EDIT

Floating point math is hard. Notice this:

int main(void)
{
   double r1, r2;
 
   r1 = sqrt (2.0);
   r2 = mysqrt (2.0);
   printf ("sqrt(2):\nr1 = %.40lf\nr2 = %.40lf\n", r1, r2);

   r1 *= r1;
   r2 *= r2;
   printf ("Square them:\nr1 = %.40lf\nr2 = %.40lf\n", r1, r2);

   r1 = fabs(r1 - 2.0);
   r2 = fabs(r2 - 2.0);
   printf ("Subtract 2 (abs):\nr1 = %.40lf\nr2 = %.40lf\n", r1, r2);

   return EXIT_SUCCESS;
}

The output is:

sqrt(2):
r1 = 1.4142135623730951454746218587388284504414
r2 = 1.4142135623730949234300169337075203657150
Square them:
r1 = 2.0000000000000004440892098500626161694527
r2 = 1.9999999999999995559107901499373838305473
Subtract 2 (abs):
r1 = 0.0000000000000004440892098500626161694527
r2 = 0.0000000000000004440892098500626161694527

So, the squared version of the two approximations of sqrt(2) are equally distant from 2. The same happens for sqrt(8), for example.

EDIT 2

Floating point math is really hard. Your function enters an infinite loop for some inputs. Try for example 7.0. You can fix it as follows:

double mysqrt(double r)
{
    double n1 = r;
    double n2 = 1.0;
    double old = n2;
    while (old != n1 && n1 - n2 > 0) {
        old = n1;
        n1 = (n1 + n2) / 2.0;
        n2 = r / n1;
    }
    return n1;
}