(Approximate) integer square root

Posted on 2020-09-04 by Oleg Grenrus

Quoting Wikipedia article: In number theory, the integer square root (intSqrt) of a positive integer n is the positive integer m which is the greatest integer less than or equal to the square root of n ,

\mathsf{intSqrt}\, n = \left\lfloor \sqrt{n} \right\rfloor

How to compute it in Haskell? The Wikipedia article mentions Newton’s method, but doesn’t discuss how to make the initial guess.

In base-4.8 (GHC-7.10) we got countLeadingZeros function, which can be used to get good initial guess.

Recall that finite machine integers look like

n = 0b0......01.....
      ^^^^^^^^       -- @countLeadingZeros n@ bits
              ^^^^^^ -- @b = finiteBitSize n - countLeadingZeros n@ bits 

We have an efficient way to get ``significant bits” count b , which can be used to approximate the number.

2^b \le n \le 2^{b+1}, \qquad n > 0

It is also easy to approximate the square root of numbers like 2^b :

\sqrt{2^b} = 2^{\frac{b}{2}} \approx 2^{\left\lfloor \frac{b}{2} \right\rfloor}

We can use this approximation as the initial guess, and write simple implementation of intSqrt:

module IntSqrt where

import Data.Bits

intSqrt :: Int -> Int
intSqrt 0 = 0
intSqrt 1 = 1
intSqrt n = case compare n 0 of
    LT -> 0           -- whatever :)
    EQ -> 0
    GT -> iter guess  -- only single iteration
  where
    iter :: Int -> Int
    iter 0 = 0
    iter x = shiftR (x + n `div` x) 1 -- shifting is dividing

    guess :: Int
    guess = shiftL 1 (shiftR (finiteBitSize n - countLeadingZeros n) 1)

Note, I do only single iteration1. Is it enough? My need is to calculate square roots of small numbers. We can test quite a large range exhaustively. Lets define a correctness predicate:

correct :: Int -> Int -> Bool
correct n x = sq x <= n && n < sq (x + 1) where sq y = y * y

Out of hundred numbers

correct100 = length
    [ (n,x) | n <- [ 0..99 ], let x = intSqrt n, correct n x ]

the computed intSqrt is correct for 89! Which are the incorrect ones?

incorrect100 =
    [ (8,3)
    , (24,5)
    , (32,6), (33,6), (34,6), (35,6)
    , (48,7)
    , (80,9)
    , (96,10), (97,10), (98,10), (99,10)
    ]

The numbers which are close to perfect square ( 8 + 1 = 3^2 , 24 + 1 = 5^2 , …) are over estimated.

If we take bigger range, say 0...99999 then with single iteration 23860 numbers are correct, with two iterations 96659.

For my usecase (mangling the size of QuickCheck generators) this is good enough, small deviations are very well acceptable. Bit fiddling FTW!


  1. Like infamous Fast inverse square root algorithm, which also uses only single iteration, because the initial guess is very good,↩︎

Site proudly generated by Hakyll