-1

Playing around with formulas in C, I realized that I found a formula for calculating the square root of a number. I would like to know if such an algorithm already exists, or if it is widely known to scholarly mathematicians. I'm sending the code so you guys take a look. I tested it in C++ Builder with TimeSpan and it is almost as fast as the language's standard sqrt function, which is written in assembly. If you can take a look and see if this is interesting I would appreciate it. It's for a school assignment.

Ps: for most numbers it gets the precision of the sqrt function to about 20 iterations.

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

int findInitial(double number){

    int i,n,res;
    n = trunc(number);
    res = 0;
    i = 1;
    while (1){
        if ((i * i) >= n) {
           res = i - 1;
           break;
        }
        i++;
    }
    return res;
}

int main(int argc, char *argv[])
{
    int i = 0;
    double number = 23;
    int initial = findInitial(number);
    double f = number;
    double e;
    double temp = 1;
    
    printf("%.18f\n",sqrt(number));

    while ((temp < -0.000000000000000001) ^ (temp > 0.000000000000000001)){
        e = f - (f * f - number)/(f - initial);
        if (temp == ((f - e) * -1)) {
          break;
        }
        temp = f - e;
        f = e;
        i++;
        printf("%d - %.18f\n",i,f*-1);
    }

    system("pause");    
    return 0;
}

Ps2: I had to create a conditional because in the case of number 23, the variable temp oscillated from negative to positive and never reached the desired precision.

vanzuita
  • 37
  • 5
  • 1
    this looks like a basic variant of Newton's method that's not guaranteed to converge – Marcus Müller Dec 12 '21 at 22:31
  • `int n = trunc(number);` is UB when `trunc(number)` is much more than `INT_MAX`. – chux - Reinstate Monica Dec 12 '21 at 22:32
  • You want to check https://en.wikipedia.org/wiki/Newton%27s_method#Square_root ; your code is a bit ummm obfuscated, but what you really do is "next iterations's `f` is last iteration's `f` minus f²-number), divided by a factor proportional to the derivative of x²-a. This is really Newton's method with unclean notation. – Marcus Müller Dec 12 '21 at 22:33
  • 1
    this doesn't necessarily converge. Newton converges for all strictly positive numbers. The lack of guaranteed convergense is due to the bad choice of divisor; should be `2*f`instead of `f-initial`. – Marcus Müller Dec 12 '21 at 22:35
  • You have indeed re-invented Newton's method. For a quadratic equation it will work nicely, unless there is a double root. Your stopping test could be more elegant if you used `std::numeric_limits`. Maybe you can eliminate the oscillation test by starting always to the right of the root. – Victor Eijkhout Dec 12 '21 at 22:38
  • 1
    @MarcusMüller Thanks. I hadn't spotted that. So it's an "over-relaxed" Newton's method in many cases. – Victor Eijkhout Dec 12 '21 at 22:39
  • @VictorEijkhout `std::numeric_limits` is [tag:C++], but this is [tag:C]. – Marcus Müller Dec 12 '21 at 22:40
  • _"almost as fast as the language's standard sqrt function"_ - that seems unlikely, the FPU has a `FSQRT` instruction - the standard library implementation operation is performed in hardware. You might usefully post a question about your test method, it is clearly flawed. If you measured the exact code posted in your question, then what you were measuring is dominated by the `printf()` call which will be many orders of magnitude longer than either your algorithm or `sqrt()`. – Clifford Dec 12 '21 at 23:32
  • @Clifford It wasn't this code that I tested, it was one that I tested in the C++ Builder compiler using the methods I researched on the Embarcadero website. (I took the printfs out) but I really don't know if my way of testing is correct. I simply put a variable of type TStopwatch, initialized it with the StartNew() method and after the function I gave a Stop(). I initialized a variable of type TTimeSpan, took the Elapsed time and called the Ticks method. Gave a similar number. Maybe the sqrt() function is so fast that it can't calculate right. – vanzuita Dec 12 '21 at 23:51
  • Testing a single call is unlikely to give an accurate impression. You should perhaps call the code iteratively, say 10000 calls. The resolution of TTimeSpan.ticks() is 100ns which is significant cf to a single instruction cycle. Generally if you make an assertion about a test you have performed, that is the code you should post. With respect to sqrt() being written in assembler; the machine code generated from C compilation is often just as fast (if optimised); there is nothing magically fast about assembly code, its speed is in the hands of the author, and the compiler optimized is expert. – Clifford Dec 13 '21 at 07:59
  • Your `find initial` is much much slower than binary based sqrt ... so you most likely did not measure it properly ... there is no way this is faster than FPU instruction... see [Power by squaring for negative exponents](https://stackoverflow.com/a/30962495/2521214) and try [tbeg(),tend() from here](https://stackoverflow.com/a/70142282/2521214) for the time measurement – Spektre Dec 13 '21 at 11:53

1 Answers1

4

Aside from this not working for all legal inputs (so it's, strictly speaking, not a method for finding the square root of numbers in general, just for a few of the possible numbers), this is basically Newton's method for approximating the square root.

Using Newton's method for this is like the standard example for a numerical math lecture, so Wikipedia has a concrete section on Newton's Method for finding Square Roots.

Your method is not quite Newton's method, because the convergence of that works for all inputs, by the fortune of dividing by the derivative of the function you want to find a null of (which is f*f-number in your case).

Marcus Müller
  • 34,677
  • 4
  • 53
  • 94
  • Thanks for the explanations. I took a look at the proposed wiki and it's definitely Newton's method with a bad derivative. Also Newton's method converges faster. I don't know if it converges for all the numbers either, at least the ones I tested do. – vanzuita Dec 12 '21 at 23:07
  • @vanzuita that's the good thing, the math *guarantees* it converges :) f²-number fulfills the conditions on the continuity, absence of further nulls, and bounding of derivative that it does for all strictly positive numbers. :) – Marcus Müller Dec 13 '21 at 11:48