Square root algorithm in C

Question

For practicing reasons I coded the following algorithm in C:

#include <stdio.h>
#include <stdlib.h>


int main(){

double x=0;

printf("Enter the number: ");

scanf("%lf", &x);

int i = 0;

double v = 0;
double n=0;
int grenze = 12;
double z = 10;
/*
for(i=1; i<(x/2+1); i++){
    v=i;
    if((v*v) <= x){
        n = i;
    }
}

v=n;
*/
for(i=1; i<grenze+1; i++){
    z = z * 0.1;

    while(v*v<x){
        v = v + z;
        if(v*v<x){
            n = v;
        }
    }
    v=n;
}
printf("%.10f\n", n);



}

This works quite well, but for number greater than a certain value (I do not know when this starts) for example 50.000.000.000 the program freezes.

Is there anything I do not see here?

Have you used a debugger to help you? It's designed exactly to help for delving into such problems. — kaylum
– kaylum, Commented Oct 19, 2016 at 23:38
What does that have to do with using a debugger? A debugger (ie, gdb) can run in the terminal. — kaylum
– kaylum, Commented Oct 19, 2016 at 23:41
For a big enough number, v = v + z; is a no-op — v doesn't change. That makes your loop run for a long time. The algorithm is not a good choice. I'm not convinced it will work well for numbers like 1E-70 (and you've shown it won't work for numbers like 1E+70). Note that doubles have a range up to 1E±300 or more on most machines with IEEE floating point math support. — Jonathan Leffler
– Jonathan Leffler, Commented Oct 19, 2016 at 23:42
Floating point numbers have a finite number of decimal digits — around about 16 for double, usually. If you add 1E+16 and 1E-16, then the result is 1E+16; there aren't enough significant digits to store the extra information. You have 5E+10 as the start point; you generate fractions z = 1.0 then 0.1, 0.01, ... 0.000 000 000 01 or thereabouts (what's an order of magnitude between friends?). When you go around adding the smallest values, you don't end up changing anything. — Jonathan Leffler
– Jonathan Leffler, Commented Oct 19, 2016 at 23:52
The canonical document on the subject is What Every Computer Scientist Should Know About Floating-Point Arithmetic or from the ACM. And there should be a Q&A for which this can be a duplicate. I don't have it identified yet. I'll transfer my comments into an answer. — Jonathan Leffler
– Jonathan Leffler, Commented Oct 19, 2016 at 23:58

Jonathan Leffler · Accepted Answer · 2016-10-20 00:21:29Z

For a big enough number v and small enough number z, v = v + z; is a no-op — v doesn't change. That makes your loop run for a long time.

The algorithm is not a good choice — you should look up the Newton-Raphson method. I'm not convinced your algorithm will work well for extreme numbers like 1E-70 (and you've shown it won't work for numbers like 1E+70). Note that doubles have a range up to 1E±300 or more on most machines with IEEE floating point support.

Could you maybe explain why for a certain limit v = v + z is a no-op please?

Floating point numbers have a finite number of decimal digits — around about 16 for double, usually. If you add 1E+16 and 1E-16, then the result is 1E+16; there aren't enough significant digits to store the extra information. You have 5E+10 as the start point; you generate fractions z = 1.0 then 0.1, 0.01, … 0.000 000 000 01 or thereabouts (what's an order of magnitude between friends?). When you go around adding the smallest values, you don't end up changing anything.

The canonical document on the subject is What Every Computer Scientist Should Know About Floating-Point Arithmetic or from the ACM.

Jonathan Leffler · Accepted Answer · 2016-10-20 00:06:57Z

0

If you write a printf statement in the while loop to show you v, you will see that at some point it stops changing (which is bad, because you use it to exit your loop!).

I put printf("DEBUG: entered for loop, v=%lf\n", v); in the while loop and saw this:

DEBUG: entered for loop, v=173205.080757
DEBUG: entered for loop, v=173205.080757
DEBUG: entered for loop, v=173205.080757
DEBUG: entered for loop, v=173205.080757
DEBUG: entered for loop, v=173205.080757
DEBUG: entered for loop, v=173205.080757
DEBUG: entered for loop, v=173205.080757
DEBUG: entered for loop, v=173205.080757

You were adding z which is so small that when v is at 173205 it doesn't actually change the value. You can extend the range by changing the type of v to long double; however, as we only have a certain number of bits to work with, there is still an upper bound where adding a small enough number to it won't change the representable value.

edited Oct 20, 2016 at 0:06

Jonathan Leffler

759k145 gold badges961 silver badges1.3k bronze badges

answered Oct 20, 2016 at 0:04

Greg Schmit

4,6752 gold badges25 silver badges39 bronze badges

2 Comments

Jonathan Leffler Over a year ago

Changing to long double is a temporary band aid — change the source number from 5E10 to 5E50 and long double will be no help. It also won't resolve the problems for tiny source numbers.

Greg Schmit Over a year ago

Very true Jonathan, and I quickly edited it to reflect that reality.

Collectives™ on Stack Overflow

Square root algorithm in C

2 Answers 2

Comments

2 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

Comments

2 Comments

Your Answer

Sign up or log in

Post as a guest

Related