//
//   integer sqrt algorithm.
//
// you should find this file(s) in
// http://www.jjj.de/
//
//   *** isqrt: ***
//
//    ENTRY x: unsigned long
//    EXIT  returns floor(sqrt(x) * pow(2,NN))
//
//    Since the square root never uses more than half the bits
//    of the input, we use the other half of the bits to contain
//    extra bits of precision after the binary point.
//
//    EXAMPLE (for BITSPERLONG=32 and NN=16):
//     then    isqrt(144) = 786432 = 12.0 * 65536
//             isqrt(32)  = 370727 = 5.66 * 65536
//
//    NOTES
//     (1) change NN to 0 if you do not want
//         the answer scaled.  Indeed, if you want NN bits of
//         precision after the binary point, use BITSPERLONG/2+NN.
//         The code assumes that BITSPERLONG is even.
//     (2) This is really better off being written in assembly.
//         These 3 lines are really an "arithmetic shift left"
//         on the double-long value with r in the upper half
//         and x in the lower half.  This operation is typically
//         expressible in only one or two assembly instructions.
//     (3) Unrolling this loop is probably not a bad idea.
//
//    ALGORITHM
//     The calculations are the base-two analogue of the square
//     root algorithm we all learned in grammar school.  Since we're
//     in base 2, there is only one nontrivial trial multiplier.
//
//     Notice that absolutely no multiplications or divisions are performed.
//     This means it'll be fast on a wide range of processors.
//
//------------------------------------------------------------

#include <math.h>

#define BITSPERLONG sizeof(long)*8

#define TOP2BITS(x) (x>>(BITSPERLONG-2))

#define NN 8  // range: 0...BITSPERLONG/2


inline unsigned long isqrt(unsigned long x, unsigned long &r)
{
    int i;
    unsigned long a = 0L;               /* accumulator */
    unsigned long e = 0L;               /* trial product */

    r=0;                                /* remainder */
    for (i=0; i<BITSPERLONG/2+NN; i++)  /* NOTE 1 */
    {
        r <<= 2;                        /* NOTE 2 */
        r +=  TOP2BITS(x);
        x <<= 2;

        a <<= 1;
        e  =  (a<<1)+1;

        if ( r>=e )
        {
            r-=e;
            a++;
        }
    }

    return a;
}
//----------------------------------------------------------------


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

void func(unsigned long i)
{
    unsigned long sq,rm;
    double dsq;

    sq=isqrt(i, rm );
    printf("\n sqrt(%6lu)=%6lu, rem=%6lu ", i, sq>>NN, rm);
    dsq=(double)sq/(1<<NN);
    printf(" ==%-18.12g  =?=%-.12g", dsq, sqrt((double)i));

    return;
}
//----------------------------------------------------------------


main(void)
{

    for (unsigned long i=0; i<37; ++i)
        func(i);


    unsigned long a=2749;

    func(a*a-1);
    func(a*a);
    func(a*a+1);
    func(a*a+2);
    func(a*a+51);

    return 0;
}
//----------------------------------------------------------------