1

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.

WinterWeb
  • 31
  • 1
  • 5
  • Your two values (r1 and r2) are within the 'machine tolerance' of being equal to each other. Print out the value of `nextafter(r2, INFINITY)` and you will get the *exact* value of r1. Thus, the two values are either side of the actual value of sqrt(2). – Adrian Mole Dec 14 '22 at 10:29
  • 1
    `sqrt` can be approximated in several ways. And the results you got are very similar anyway. What did you expect ? – wohlstad Dec 14 '22 at 10:30
  • 1
    Although I don't understand the downvotes on this question, it *could* be considered a duplicate of [Is floating point math broken?](https://stackoverflow.com/q/588004/10871073) – Adrian Mole Dec 14 '22 at 10:31
  • Floating-point math is often inexact. It is also often difficult.It can be hard to get a result that's accurate down to the last bit or digit, and even harder to guarantee a result that's perfectly rounded — that is, closer to the actual value than any other result could be. – Steve Summit Dec 14 '22 at 11:34
  • The fact that you've succeeded on the first point — your result is one of the two that's closest to the true value — is cause for considerable celebration. Choosing between those two, to find the one that's closer to the true value — would require some considerably fancier techniques (that I for one have no idea how to do). – Steve Summit Dec 14 '22 at 11:47

2 Answers2

3

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;
}
Costantino Grana
  • 3,132
  • 1
  • 15
  • 35
0

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

WinterWeb
  • 31
  • 1
  • 5