Implementation of the square root for real time application

Question

For the real time control software I need square root calculation. I have heard that the sqrt function from the standard library isn't good choice due to the unpredictable number of iterations used during calculation. Based on that information I have decided to implement the sqrt function in my own based on the Newton tangent method. So I have exploited the iteration formula

$$x_{k+1} = -\frac{f(x_k)}{f'(x_k)} + x_k,$$ where $$f(x) = x^2 - a$$ and $$a\in\mathbb R_0^+, a_{min} = 0.1, a_{max} = 400.0$$.

Based on the above mentioned formulas I have

$$x_{k+1} = \frac{1}{2}\cdot\left(x_k + \frac{a}{x_k}\right)$$

I have chosen $$x(0) = 2$$ and based on my observations four iterations give good results in comparison with the sqrt from the standard library.

#define ITERATIONS 4

double squareRoot(double a)
{
    double xk  = 2.0; // x(k)
    
    uint8_t k;
    for(k = 0; k < ITERATIONS; k++)
    {
        xk  = 0.5*(xk + a/xk);
    }
    
    return xk;
}

Despite that fact I have some doubts regarding the choice of the initial value and number of iterations.

Answer 1

I have heard that the sqrt function from the standard library isn't good choice due to the unpredictable number of iterations used during calculation.

The goal of standard functions do not generally include a uniform time requirement.

Far more often, a precise and correct as able solution is sought. Speed is of secondary concern.

---

my observations four iterations give good results

Below is a test harness.

I found the numeric results disappointing for general double usage. Notice the "relative difference" of 1.0 (very bad).

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

// Insert OP's code here

void test(double x) {
  double y0 = sqrt(x);
  double y1 = squareRoot(x);
  if (y0 != y1) {
    double diff = y1 - y0;
    double rel = diff / (y1 + y0);
    printf(
        "x:%13g sqrt(x):%23.17g squareRoo(x):%23.17g relative difference:%13g\n",
        x, y0, y1, rel);
  }
}

int main() {
  double a[] = {0.0, 1.0, 2.0, 42, DBL_TRUE_MIN, DBL_MIN, sqrt(DBL_MIN), sqrt(
      DBL_MAX), DBL_MAX};
  size_t n = sizeof a / sizeof a[0];
  for (size_t i = 0; i < n; i++) {
    test(a[i]);

  }
  puts("Done");
}

Output

x:            0 sqrt(x):                      0 squareRoo(x):                  0.125 relative difference:            1
x:            1 sqrt(x):                      1 squareRoo(x):     1.0000000464611474 relative difference:  2.32306e-08
x:            2 sqrt(x):     1.4142135623730951 squareRoo(x):     1.4142135623746899 relative difference:   5.6382e-13
x:           42 sqrt(x):     6.4807406984078604 squareRoo(x):     6.4812185874674091 relative difference:  3.68686e-05
x: 4.94066e-324 sqrt(x):2.2227587494850775e-162 squareRoo(x):                  0.125 relative difference:            1
x: 2.22507e-308 sqrt(x):1.4916681462400413e-154 squareRoo(x):                  0.125 relative difference:            1
x: 1.49167e-154 sqrt(x):  1.221338669755462e-77 squareRoo(x):                  0.125 relative difference:            1
x: 1.34078e+154 sqrt(x): 1.1579208923731618e+77 squareRoo(x):4.1899399781070611e+152 relative difference:            1
x: 1.79769e+308 sqrt(x):1.3407807929942596e+154 squareRoo(x):5.6177910464447366e+306 relative difference:            1
Done

---

Design problem

Newton tangent method rapidly converges to a good answer once it is is near the correct result. Trouble with OP's approach is the slow convergence when the result is not near 2.0.

Instead, better to first get the result exponent in range.

Below, with exponent halving and 1 more iteration that OP's, give very good result for all double > 0. Notice the "relative difference" of 1e-16 (very good).

double squareRoot_improved(double a) {
  int expo;
  frexp(a, &expo);
  double xk = ldexp(1, expo/2);

  uint8_t k;
  for (k = 0; k < 5; k++) {
    xk = 0.5 * (xk + a / xk);
  }

  return xk;
}

Results

x:            0 sqrt(x):                      0 squareRoo(x):                0.03125 relative difference:            1
x:            2 sqrt(x):     1.4142135623730951 squareRoo(x):     1.4142135623730949 relative difference: -7.85046e-17
x: 4.94066e-324 sqrt(x):2.2227587494850775e-162 squareRoo(x):  2.22275874948508e-162 relative difference:  5.55112e-16
x: 2.22507e-308 sqrt(x):1.4916681462400413e-154 squareRoo(x): 1.491668146240043e-154 relative difference:  5.55112e-16
x: 1.49167e-154 sqrt(x):  1.221338669755462e-77 squareRoo(x): 1.2213386697554618e-77 relative difference: -7.85046e-17
x: 1.34078e+154 sqrt(x): 1.1579208923731618e+77 squareRoo(x):  1.157920892373162e+77 relative difference:  5.55112e-17
x: 1.79769e+308 sqrt(x):1.3407807929942596e+154 squareRoo(x):1.3407807929942597e+154 relative difference:  5.55112e-17
Done

---

See also Fast inverse square root.