Calculation of integer square root

There was a need to check whether the integer is a square, and if so, then calculate the root. And I want to do this in integer arithmetic. It is clear that you can implement the Newton method in integers, but it requires division at every step. But can it be otherwise? Find the square root modulo the power of two, and check if it will be an ordinary square root.



You can restrict yourself to odd numbers: for an even number, if the number of low-order zero digits is odd, then there is no root, but if it is even, then you can shift the number to the right, calculate the root of the odd, and shift back to the left half the original number of zero bits.



For odd N and 2 k , k> 3, if N ≡ 1 mod 8, then there are 4 different roots modulo 2 k , otherwise there are no roots. We need the smallest of the four roots of x. The other three roots are 2 k - x, 2 k-1 + x and 2 k - 2 k-1 - x



I would like something similar to calculating the inverse modulo 2 k - doubling the number of valid bits per iteration.



Suppose we already have a root x 0 of N modulo 2 k : N - x 0 2 = 2 k a

And we want to find x 1 = x 0 + 2 k-1 y such that in N - x 1 2 there are more lower zero bits.

N - (x 0 + 2 k-1 y) 2 = 2 k a - 2 k x 0 * y - 2 2k-2 y 2

Divide by 2 k : a - x 0 * y - 2 k-2 y 2

And equate to 0 modulo 2 k-2 : y = a * x 0 -1 mod 2 k-2

Received x 1 = x 0 + 2 k-1 a * (x 0 -1 mod 2 k-2 )

And finally x 1 = x 0 + (N - x 0 2 ) / 2 * (x 0 -1 mod 2 k-2 )



Of the k bits in the next iteration, you get 2 (k-1) bits. At the same time, we consider the inverse to the root at each iteration.



Test code:



uint8_t sqr16(uint16_t n) { if (n % 8 != 1) return 0; uint16_t sqr = (n + 1) / 2; //4 bit uint16_t inv = 2 - sqr; sqr += inv * (n-sqr*sqr)/2; //6 bit inv *= 2 - sqr * inv; sqr += inv * (n-sqr*sqr)/2; //10 bit //inv *= 2 - sqr * inv; if (sqr & 256) sqr = 0u - sqr; sqr = (uint8_t)sqr; // lowest root if (n == sqr*sqr) return sqr; return 0; }
      
      





Adding a couple of iterations, we get the root from uint_64



All Articles