Fast approximation of a Square Root?

Question

I'm playing around with some calculations in a program I'm writing. Although my problem may seem obscure and perhaps avoidable, it's not, and even if it is, I'd rather not avoid it for now :-D

The problem is that I have a huge number (thousands of digits, maybe millions at some point) which I need to find an approximate SQRT for. I'm using System.Numerics.BigInteger, I'm need the approximation to be less than or equal to the real SQRT. For example, if the number provided was 100, I'd be happy with 7, 8, 9 or 10, but not 11.

Current Code:

public static BigInteger IntegerSqrt(BigInteger n)
{
    var oldValue = new BigInteger(0);
    BigInteger newValue = n;

    while (BigInteger.Abs(newValue - oldValue) >= 1)
    {
        oldValue = newValue;
        newValue = (oldValue + (n / oldValue)) / 2;
    }

    return newValue;
}

Whilst this gives me the right answer (so far as I can tell), it's crazy slow as it's trying to get a precise answer.

If getting the cube-root, or anything else like that would be faster, I'd be happy with that. Also, yes, I know finding square roots quickly can make you lots of money, not really what I'm trying to do though, I just want a quick approximation...

Any help you could give me would be great!

Edit - Gabe

This is the updated IntegerSqrt function I am using which doesn't appear to work any quicker.

public static BigInteger IntegerSqrt(BigInteger n)
{
    var oldValue = n.RoughSqrt();
    BigInteger newValue = n;
    while (BigInteger.Abs(newValue - oldValue) >= 1)
    {
        oldValue = newValue;
        newValue = (oldValue + (n / oldValue)) / 2;
    }

    return newValue;
}

Edit 2 Is this what you were thinking? - I ran this for 30 minutes on a sample large number (50k digits) and it looped around 100 times without completing. I feel as though I'm missing something...

public static BigInteger IntegerSqrt(BigInteger n)
{
    var oldValue = new BigInteger(0);
    BigInteger newValue = n.RoughSqrt();
    int i = 0;
    while (BigInteger.Abs(newValue - oldValue) >= 1)
    {
        oldValue = newValue;
        newValue = (oldValue + (n / oldValue)) / 2;
        i++;
    }
    return newValue;
}

Answer 1

Using my answer here to calculate Log₂(N), I prepared a test case.

Although your code is closer to the real square root, my tests shows millions of times speed up when the input number is getting larger.

I think you will have to decide between more precise result vs. thousands/millions of times speed-up,

Random rnd = new Random();
var bi = new BigInteger(Enumerable.Range(0, 1000).Select(x => (byte)rnd.Next(16))
                                                 .ToArray());


var sqrt1 = bi.Sqrt().ToString();
var sqrt2 = IntegerSqrt(bi).ToString();

var t1 = Measure(10 * 200000, () => bi.Sqrt());
var t2 = Measure(10, () => IntegerSqrt(bi));

The extension method for BigInteger.Sqrt

public static class BigIntegerExtensions
{
    static int[] PreCalc = new int[] { 8, 7, 6, 6, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
    static int Log2(this BigInteger bi)
    {
        byte[] buf = bi.ToByteArray();
        int len = buf.Length;
        return len * 8 - PreCalc[buf[len - 1]] - 1;
    }

    public static BigInteger Sqrt(this BigInteger bi)
    {
        return new BigInteger(1) << (Log2(bi) / 2);
    }
}

Code to measure the speed

long Measure(int n,Action action)
{
    action(); 
    var sw = Stopwatch.StartNew();
    for (int i = 0; i < n; i++)
    {
        action();
    }
    return sw.ElapsedMilliseconds;
}

RESULTS:

In my computer, both code results ~in 8 seconds, but my posted code runs 10*200000 times, comparing your code runs only 10 times. (Yes, hard to believe, for 8000 bit numbers the difference is 200,000 times)

When comparing with 16000 bit numbers the difference raises to 1,000,000 times...

Answer 2

Here's the rough approximation that will get you to within a factor of 2 of your square root:

System.Numerics.BigInteger RoughSqrt(System.Numerics.BigInteger x)
{
    var bytes = x.ToByteArray();    // get binary representation
    var bits = (bytes.Length - 1) * 8;  // get # bits in all but msb
    // add # bits in msb
    for (var msb = bytes[bytes.Length - 1]; msb != 0; msb >>= 1)
        bits++;
    var sqrtBits = bits / 2;    // # bits in the sqrt
    var sqrtBytes = sqrtBits / 8 + 1;   // # bytes in the sqrt
    // avoid making a negative number by adding an extra 0-byte if the high bit is set
    var sqrtArray = new byte[sqrtBytes + (sqrtBits % 8 == 7 ? 1 : 0)];
    // set the msb
    sqrtArray[sqrtBytes - 1] = (byte)(1 << (sqrtBits % 8));
    // make new BigInteger from array of bytes
    return new System.Numerics.BigInteger(sqrtArray);
}

Instead of starting with n as your approximation (newValue = n;), you would start off by calling this (newValue = RoughSqrt(n);). This will get you to your answer with way fewer iterations.

If you want to take the Nth root instead, you just have to change the rough approximation to use 1/Nth the bits and change Newton's method to use the derivative of x^N:

static BigInteger RoughRoot(BigInteger x, int root)
{
    var bytes = x.ToByteArray();    // get binary representation
    var bits = (bytes.Length - 1) * 8;  // get # bits in all but msb
    // add # bits in msb
    for (var msb = bytes[bytes.Length - 1]; msb != 0; msb >>= 1)
        bits++;
    var rtBits = bits / root + 1;   // # bits in the root
    var rtBytes = rtBits / 8 + 1;   // # bytes in the root
    // avoid making a negative number by adding an extra 0-byte if the high bit is set
    var rtArray = new byte[rtBytes + (rtBits % 8 == 7 ? 1 : 0)];
    // set the msb
    rtArray[rtBytes - 1] = (byte)(1 << (rtBits % 8));
    // make new BigInteger from array of bytes
    return new BigInteger(rtArray);
}

public static BigInteger IntegerRoot(BigInteger n, int root)
{
    var oldValue = new BigInteger(0);
    var newValue = RoughRoot(n, root);
    int i = 0;
    // I limited iterations to 100, but you may want way less
    while (BigInteger.Abs(newValue - oldValue) >= 1 && i < 100)
    {
        oldValue = newValue;
        newValue = ((root - 1) * oldValue 
                    + (n / BigInteger.Pow(oldValue, root - 1))) / root;
        i++;
    }
    return newValue;
}

Answer 3

Firstly, write your number in the approximate form

a = n.10^2m

Where n is an integer and as large as you can make it subject to a system sqrt being available for its size. m is also an integer. The dot means multiply.

You should be able to get m from the number of digits in the initial integer and the number of digits in n.

Note that this will never be greater than the given number (as we have set many of the less significant digits to 0): so satisfying the requirement of the answer being a lower bound.

The square root is

a^1/2 = n^1/2.10^m

This is fast to compute: the first part using the system sqrt; the second part is a trivial multiplication; well within the scope of BigInteger.

The accuracy of this approach approaches floating point precision and its speed is independent of the size of the number.

Answer 4

√ N ≈ ½(N/A + A)

That is Newtons Square Root Approximation - found using...Google :-P

EDIT

The function below returns

7.9025671444237E+89 <-approx sqrt 7.9019961626505E+89 <-actual sqrt

for

624415433545435435343754652464526256453554543543254325247654435839457324987585439578943795847392574398025432658763408563498574398275943729854796875432457385798342759843263456484364

Here is what I came up with (proud of myself for the logic - not so much for the coding :-P)

Ok so this is PHP (was fun to do this :-D) but should be easy to convert.

Looks scary but it isn't

1. find first 3/4 digits (depending on string length).

2. Find the approx sqrt of those digits using closest squares technique (e.g. to find approx sqrt of 45 we know nearest whole number squares (perfect squares - forget term) are 6 (36) and 7 (49) - and that 45 is about 70% of the way between 36 and 49 - so 6.7 is good approximation)

3. Stick this through Newtons square approximation formula.

What I have written is a very ugly way of doing it - with all the steps left in - but it actually makes it easier to see the thought pattern (I hope).

function approxSqrRt($number){

    $numLen = strlen($number);//get the length of the number - this is used to work out whether to get the first 3 or 4 digits from the number.

    $testNum = 0;
    $zerosToAdd = 0;


    //find out if even or odd number of digits as 400 and 40,000 yield same start 20 -> 200 BUT 4000 yields 63.245
    if($numLen % 2 == 0){
        $testNum = substr($number, 0, 4);
        $zerosToAdd = ($numLen - 4) / 2; 
    }else{
        $testNum = substr($number, 0, 3);
        $zerosToAdd = ($numLen - 3) / 2;
    }

    $square1 = 0;

    $x = 0;
    $z = false;
     //this should be made static array and look up nearest value but for this purpose it will do.
    // It is designed to find the nearest 2 squares - probably a much better way but was typing as I thought.
    while ($z == false){
        $x++;
        if(($x*$x) >= $testNum){

            $square1 = $x - 1;

            $z = true;


        }



    }

    $square2 = $square1 + 1;



    $ratio = 0;
    //ratio - we are simply qorking out if our square is closer to $square1 or $square2 and how far it is in between them - this will yield a much more accurate approximation.
    //
    if($testNum != ($square1 * $square1)){
        if($testNum != ($square2 * $square2)){

            $topOfEquation = $testNum;
            $BottomOfEquation = $square2 * $square2;
            $ratio = ($topOfEquation / $BottomOfEquation);


        }else{
            $ratio = 1;
        }
    }else{
        $ratio = 0;
    }

    //get the final approximation before adding the relevant number of 0's.
    $testNumSqrRt = $square1 + $ratio;

    //add the 0's on to make the number the correct length.
    $testNumberFinal = $testNumSqrRt  *  pow(10, $zerosToAdd);

    //run through Newtons Approximation.
    $approxSqrtFinal = (($number / $testNumberFinal) + $testNumberFinal) / 2;

   return $approxSqrtFinal;

}



$num = "624415433484364";

echo approxSqrRt($num) . "<br/>";
echo sqrt($num);

Answer 5

Looking for an efficient integer square root algorithm for ARM Thumb2

Anyway the algorithm must use shifts, +/- operations instead division.

Also if you implement the routine in assembly language, you might get some extra optimization. For using assembly in C#, please see:

http://www.codeproject.com/Articles/1392/Using-Unmanaged-code-and-assembler-in-C

Answer 6

        static public double SQRTByGuess(double num)
        {
            // Assume that the number is less than 1,000,000. so that the maximum of SQRT would be 1000.
            // Lets assume the result is 1000. If you want you can increase this limit
            double min = 0, max = 1000;
            double answer = 0;
            double test = 0;
            if (num < 0)
            {
                Console.WriteLine("Negative numbers are not allowed");
                return -1;
            }
            else if (num == 0)
                return 0;

            while (true)
            {
                test = (min + max) / 2;
                answer = test * test;
                if (num > answer)
                {
                    // min needs be moved
                    min = test;
                }
                else if (num < answer)
                {
                    // max needs be moved
                    max = test;
                }
                if (num == answer)
                    break;
                if (num > (answer - 0.0001) &&
                    num < (answer + 0.0001))
                    break;
            }
            return test;
        }

Reference: http://www.softwareandfinance.com/CSharp/Find_SQRT_Approximation.html