[Haskell-cafe] fast integer base-2 log function?

[Stefan O'Rear]

On Sun, Feb 10, 2008 at 10:15:58PM -0800, Uwe Hollerbach wrote:
> Hello, haskellers,
> 
> Is there a fast integer base-2 log function anywhere in the standard
> libraries? I wandered through the index, but didn't find anything that
> looked right. I need something that's more robust than logBase, it
> needs to handle numbers with a few to many thousands of digits. I
> found a thread from a couple of years ago that suggested there was no
> such routine, and that simply doing "length (show n)" might be the
> best. That seems kind of... less than elegant. I've come up with a
> routine, shown below, that seems reasonably fast (a few seconds of CPU
> time for a million-bit number, likely adequate for my purposes), but
> it seems that something with privileged access to the innards of an
> Integer ought to be even much faster -- it's just a simple walk along
> a list (array?) after all. Any pointers? Thanks!

Even easier.


{-# LANGUAGE MagicHash #-}
import GHC.Exts
import Data.Bits

-- experiment with using a LUT here (hint: FFI + static arrays in C)
ilog2i0, ilog2i1, ilog2i2, ilog2i3, ilog2i4 :: Int -> Int -> Int
ilog2i0 k x | x .&. 0xFFFF0000 /= 0 = ilog2i1 (k + 16) (x `shiftR` 16)
            | otherwise             = ilog2i1 k x
ilog2i1 k x | x .&. 0xFF00 /= 0     = ilog2i2 (k + 8)  (x `shiftR` 8)
            | otherwise             = ilog2i2 k x
ilog2i2 k x | x .&. 0xF0 /= 0       = ilog2i3 (k + 4)  (x `shiftR` 4)
            | otherwise             = ilog2i3 k x
ilog2i3 k x | x .&. 0xC /= 0        = ilog2i4 (k + 2)  (x `shiftR` 2)
            | otherwise             = ilog2i4 k x
ilog2i4 k x | x .&. 0x2 /= 0        = k + 1 + (x `shiftR` 1)
            | otherwise             = k + x

log2i :: Integer -> Int  -- actually returns bit length, and returns garbage on negatives, but do you care?
log2i (J# len adr) = ilog2i0 0 (I# (indexIntArray# adr (len -# 1#))) + I# (32# *# (len -# 1#))
log2i (S# sml) = ilog2i0 0 (I# sml)

> > powi :: Integer -> Integer -> Integer
> > powi b e | e == 0    = 1
> >          | e < 0     = error "negative exponent in powi"
> >          | even e    = powi (b*b) (e `quot` 2)
> >          | otherwise = b * (powi b (e - 1))
> 
> > ilog2 :: Integer -> Integer
> > ilog2 n | n < 0      = ilog2 (- n)
> >         | n < 2      = 1
> >         | otherwise  = up n (1 :: Integer)
> >   where up n a = if n < (powi 2 a)
> >                     then bin (quot a 2) a
> >                     else up n (2*a)
> >         bin lo hi = if (hi - lo) <= 1
> >                        then hi
> >                        else let av = quot (lo + hi) 2
> >                             in if n < (powi 2 av)
> >                                   then bin lo av
> >                                   else bin av hi
> 
> (This was all properly aligned when I cut'n'pasted; proportional fonts
> might be messing it up here.)
> _______________________________________________
> Haskell-Cafe mailing list
> Haskell-Cafe at haskell.org
> http://www.haskell.org/mailman/listinfo/haskell-cafe
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 189 bytes
Desc: Digital signature
Url : http://www.haskell.org/pipermail/haskell-cafe/attachments/20080210/19b63cfc/attachment.bin