less_retarded_wiki/sqrt.md
2024-09-09 14:10:06 +02:00

6.4 KiB

Square Root

Square root (sometimes shortened to sqrt) of number a is such a number b that b^2 = a, for example 3 is a square root of 9 because 3^2 = 9. Finding square root is one of the most basic and important operations in math and programming, e.g. for computing distances, solving quadratic equations etc. Square root is a special case of finding Nth root of a number for N = 2. Square root of a number doesn't have to be a whole number; in fact if the square isn't a whole number, it is always an irrational number (i.e. it can't be expressed as a fraction of two integers, for example square root of two is approximately 1.414...); and it doesn't even have to be a real number (e.g. square root of -1 is i). Strictly speaking there may exist multiple square roots of a number, for example both 5 and -5 are square roots of 25 -- the positive square root is called principal square root; principal square root of x is the same number we get when we raise x to 1/2, and this is what we are usually interested in -- from now on by square root we will implicitly mean principal square root. Programmers write square root of x as sqrt(x) (which should give the same result as raising to 1/2, i.e. pow(x,0.5)), mathematicians write it as:

  _    1/2
\/x = x

Here is the graph of square root function (notice it's a parabola flipped by the diagonal axis, for square root is an inverse function to the function x^2):

      ^ sqrt(x)
      |       :       :       :       :       :
    3 + ~ ~ ~ + ~ ~ ~ + ~ ~ ~ + ~ ~ ~ + ~ ~ ~ + ~
      |       :       :       :       :       :
      |       :       :       :       :       ___....--
    2 + ~ ~ ~ + ~ ~ ~ + ~ ~ ~ + ~ __..---'"""": ~
      |       :       : __..--""""    :       :
      |       :  __.--'"      :       :       :
    1 + ~ ~ _--'" ~ ~ + ~ ~ ~ + ~ ~ ~ + ~ ~ ~ + ~
      |  _-"  :       :       :       :       :
      | /     :       :       :       :       :
  ----+"------|-------|-------|-------|-------|----> x
      |0      1       2       3       4       5
      |

TODO

Programming

Square root is a very common operation and oftentimes has to be VERY fast -- it's used a lot for example in computer graphics where it may need to be computed several million times per second. For this reason there oftentimes exist special instructions and otherwise hardware accelerated options, you will very likely have such function around on most computers -- in C math standard library there's the sqrt floating point function that will probably be very fast. But let's now consider you want to program your own square root, which may happen for example when dealing with embedded computers.

If we really need extreme speed, we may always use a look up table with precomputed values. Of course a table with ALL the values would be very big, but remember we can make a much smaller table (where each item spans a bigger range) that will provide just a quick estimate from which you'll make just a few extra iterations towards the correct answer.

Within desired precision square root can be relatively quickly computed iteratively by binary search. Here is a simple C function computing integer square root this way:

{ I checked the code works for all values with 32 bit integer (remember that C specification allows integers to be as small as 16 bit though, remember to adjust the constants if you suspect you might hit this, overflows may happen). On my computer I measured this to be about 2 (with compiler optimization) to 4 (without) times slower than using the hardware accelerated float point stdlib function. ~drummyfish }

unsigned int sqrtInt(unsigned int x)
{
  unsigned int m, l = 0,
    r = x < 8300000 ? x / 128 + 40 : 65535; // upper bound est. for 32 bit int
    //r = x < 15000 ? x / 64 + 20 : 255;    // <-- for 16 bit int use this

  while (1)
  {
    if (l > r)
      break;

    m = (l + r) / 2;

    if (m * m > x)
      r = m - 1;
    else
      l = m + 1;
  }

  return r;
}

But now let's take a look at maybe even a better algorithm -- it only gives an approximation, however it's pretty accurate and we may modify it to give a precise value by simply adjusting the estimate at the end. The advantage is that the approximation only uses bit shifts, no multiplication! This can be crucial on some simple platforms where multiplication may be expensive. Furthermore the algorithm is pretty simple, it works like this: given input number x, we may imagine we have a rectangle of size 1 times x; now we can try to make it into square by doubling the first side and halving the other. If we get equal sides, the side length is the square root. Of course we don't always iterate to the same side sizes -- if the first side gets bigger than the other, we stop and simply average them -- and that's it! Here's the code:

{ I remember I came up with a simple form of this algorithm when I was still in elementary school, however it's pretty obvious so it probably already exists, I didn't bother checking. Measuring the speed of this the pure approximation is very fast, I measured it basically at the same speed as the stdlib float square root; the exact version is of course considerably slower -- anyway for lower input values I measured it even faster than the binary search, however for higher values not anymore. ~drummyfish }

unsigned int sqrtIntApprox(unsigned int x)
{
  unsigned int a = 1;

  while (x > a)
  {
    a <<= 1;
    x >>= 1;
  }

  return (a + x) >> 1;
}

unsigned int sqrtIntExact(unsigned int x)
{
  unsigned int a = 1, b = x;

  while (b > a)
  {
    a <<= 1;
    b >>= 1;
  }

  a = (a + b) >> 1;

  while (a * a > x) // iterate to exact value
    a--;

  return a;
}

TODO: Heron's method?

The following is a non-iterative approximation of integer square root in C that has acceptable accuracy to about 1 million (maximum error from 1000 to 1000000 is about 7%): { Painstakingly made by me. This one was even faster than the stdlib function! ~drummyfish }

int32_t sqrtApprox(int32_t x)
{
  return
    (x < 1024) ?
    (-2400 / (x + 120) + x / 64 + 20) :
      ((x < 93580) ?
       (-1000000 / (x + 8000) + x / 512 + 142) :
       (-75000000 / (x + 160000) + x / 2048 + 565));
}