#include <math.h>

/// Integer square root of x.
/// Returns s such that  s*s <= x < (s+1)*(s+1)
/// adapted from: https://web.archive.org/web/20120306040058/http://medialab.freaknet.org/martin/src/sqrt/sqrt.c
/// see also: https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Binary_numeral_system_(base_2)
/// should work for the full range of x
unsigned sqrti(unsigned long x)
{
  #if __SSE2__
    // fast path: if x not too large, convert to floating point.
    //if (x < (1UL<<24))  return sqrtf(x);	// implicit round towards zero (floor)
    //if (x < (1UL<<52))  return sqrt(x);	// implicit round towards zero (floor)
  #endif

    unsigned long s = 0;
    unsigned long bit = 1;

  #if !defined(__GNUC__) && !defined(__clang__)
    while(bit <= x/4) bit*=4;   // largest power of 4 <= x  (except for x=0)
  #else
    // faster than above if __builtin_clzll() is available:
    if (x<=1) return x;
    unsigned z = 63 - __builtin_clzll(x);
    z = z &~ 1;		// round down to even
    bit = 1UL << z;	// largest power of 4 <= x
  #endif

    do {
        unsigned long t = s+bit;
	s /= 2;
	if (x >= t) {
	    s += bit;
	    x -= t;
	}
    } while (bit /= 4);
    return s;     //  s*s <= x < (s+1)*(s+1)
}


// Here are some tests to check if everything is alright:

#include <stdio.h>

int main()
{
    // for the first integers:
    for (long i=0; i<16000000; i++) {
	if (sqrti(i) != (unsigned) sqrt(i))	printf("FAIL for i=%lu : %u != floor(%g)\n",i, sqrti(i),sqrt(i));
    }
    // for powers of 2 up to 2**63:
    for (unsigned long i=(1UL<<63); i>0; i>>=1) {
	if (sqrti(i) != (unsigned) sqrt(i))	printf("FAIL for i=%lu : %u != floor(%g)\n",i, sqrti(i),sqrt(i));
    }

    // this shows rounding problems for large integers (>2**52) where sqrt() does not lead the correct result due to rounding issues.
    for (unsigned long i=((1UL<<63)-1)<<1; i>=0; i>>=1) {
	if (sqrti(i) != (unsigned) sqrt(i))	printf("note: sqrti returns correct value while sqrt fails for i=%lu : %u != floor(%.18g)\n",i, sqrti(i),sqrt(i));
	if (i==0) break;
    }
}