#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 xunsignedsqrti(unsignedlongx){#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)#endifunsignedlongs=0;unsignedlongbit=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)returnx;unsignedz=63-__builtin_clzll(x);z=z&~1;// round down to evenbit=1UL<<z;// largest power of 4 <= x#endifdo{unsignedlongt=s+bit;s/=2;if(x>=t){s+=bit;x-=t;}}while(bit/=4);returns;// s*s <= x < (s+1)*(s+1)}// Here are some tests to check if everything is alright:#include<stdio.h>intmain(){// for the first integers:for(longi=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(unsignedlongi=(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(unsignedlongi=((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;}}
Comments (0)
HTTPSSSH
You can clone a snippet to your computer for local editing.
Learn more.