Aarsh Dave

Sometimes I think I take math.h for granted, so today I'm going to try and create my own square root function.

Since many square roots are irrational, this will be all about trying to get the best approximation that I can.

I'm not really sure where to start, so I'll try to implement the first and only approximation that I learned back in Calc 1, which is using a linear approximation.

A linear approximation

So we're trying to estimate $f(x) = \sqrt{x}$
To take a linear approximation we can find a perfect square close to $x$ , and use the tangent line to the function to approximate.

The tangent line can be found with the derivative $f'\left(x\right)=\frac{1}{2\sqrt{x}}$ Plugging this into a line equation, the tangent line for the square root function at any input $a$ is given by
$y=\frac{1}{2\sqrt{a}}\left(x-a\right)+\sqrt{a}$
Now the only thing left to find is a perfect square $a$ that is close to $x$ I'm sure there's a better way to do this, but for now I'm just going to iterate through numbers until the square is greater than the target number. Theoretically this should be alright as $f(x) = x^2$ grows fairly quickly.

Here's how I'm finding the closest perfect square
(the function returns its square root)

int closestSquare(double x)
{
  int previousSquare;
  int nextSquare = 1;

  while (nextSquare * nextSquare < x)
  {
    previousSquare = nextSquare++;
  }

  return (nextSquare * nextSquare - x) > (x - previousSquare * previousSquare) ? previousSquare : nextSquare;
}

And now we should have everything we need to implement this approximation.

double linear_sqrt(double x)
{
  int cs = closestSquare(x);
  return (1.0 / (2 * cs)) * (x - cs * cs) + cs;
}

Ok so before I try to improve this approximation, lets get a sense of how accurate this is to the sqrt() function in math.h

The error should be less for greater numbers as the square root will grow slower and the tangent line will become a better approximation.

Averaging the error for the first $n$ positive integers we get...

$n$	Average Error
$100$	$0.007871$
$1,000$	$0.002607$
$10,000$	$0.000829$

... which I think is pretty decent for a simple approximation

The second iteration

I'm not really sure of any other types of approximations, so I think the next logical step is to create a quadratic approximation using both the first a second derivatives of the function.

This will use an approximating line that looks like this
$y=-\frac{1}{8\left(\sqrt{a}\right)^{3}}\left(x-a\right)^{2}+\frac{1}{2\sqrt{a}}\left(x-a\right)+\sqrt{a}$
Where this function will have the same first and second derivative as the square root function at any input $a$ .

This is the implementation

double quad_sqrt(double x)
{
  int cs = closestSquare(x);
  return (-1.0 / (8 * power(cs, 3))) * power((x - cs * cs), 2) + (1.0 / (2 * cs)) * (x - cs * cs) + cs;
}

(I created a power function to simplify things, can't be cheating with math.h)

So let's check the error for this function...

$n$	Average Error
$100$	$0.000982$
$1,000$	$0.000135$
$10,000$	$0.000017$

...and it seems to be improving

A general form

Now to improve this to any desired accuracy it seems that we should be able to take derivatives of the function and create better and better approximations...

this sounds a lot like Taylor's Theorem from Calc 2 which means our approximation should be able to be extended in the following way...

$f\left(x\right)≈f\left(a\right)+f'\left(a\right)\left(x-a\right)+\frac{f''\left(a\right)}{2!}\left(x-a\right)^{2}+...+\frac{f^{\left(k\right)}\left(a\right)}{k!}\left(x-a\right)^{k}$
to implement this we need a factorial function and a way to take an arbitrary amount derivatives of the square root function.

The factorial is trivial

int factorial(int n)
{
  return n > 2 ? n * factorial(n - 1) : n;
}

(I'm just using int's for now because I don't expect to plugin in anything too high)

and for the derivative of $f\left(x\right)=\sqrt{x}$ , I'll implement the power rule directly in the loop.

And the completed function looks like

double general_sqrt(double x, int n)
{
  int cs = closestSquare(x);
  double approx = cs;

  double coefficient = 0.5;

  for (int k = 1; k <= n; k++)
  {
    approx += power(x - cs * cs, k) / factorial(k) * coefficient * power(cs, 1 - 2 * k);
    coefficient *= (1 - 2 * k) / 2.0;
  }

  return approx;
}

Here's a table showing the average error for different amounts of iterations:

	1	2	3	4	5
$n=100$	$0.007871$	$0.000982$	$0.000322$	$0.000178$	$0.000122$
$n=1,000$	$0.002607$	$0.000135$	$0.000033$	$0.000018$	$0.000012$
$n=10,000$	$0.000829$	$0.000017$	$0.000003$	$0.000002$	$0.000001$

I think I'll call that a success!

Wrapping Up

Although using Taylor's theorem is fairly accurate, I think its actually too slow compared to modern methods. Modern algorithms as far as I can tell make use of Newton's method after making an approximation to improve said approximation. Newton's method is another interesting bit of calculus but I'll leave that for another time.

Source code if you're interested