[Haskell-cafe] fast integer base-2 log function?
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! Uwe > 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@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] fast integer base-2 log function?
Thanks, guys! It looks at first glance as if the code Thorkil posted is similar to mine (grow comparison number in steps of 2 in the exponent, then binary-search to get the exact exponent), while Stefan's version is more similar to the walk-the-list idea I had in mind. I'll play with both of these when I get a chance... Uwe On Feb 10, 2008 10:44 PM, Thorkil Naur <[EMAIL PROTECTED]> wrote: > Hello, > > If the standard libraries provide such a function, I haven't found it. I > must > admit that I haven't studied your code in detail. I usually do as follows > for > integer logarithms, shamelessly stolen from the Haskell report: > > > -- Integer log base (c.f. Haskell report 14.4): > > > > imLog :: Integer->Integer->Integer > > imLog b x > > = if x < b then > > 0 > > else > > let > > l = 2 * imLog (b*b) x > > doDiv x l = if x < b then l else doDiv (x`div`b) (l+1) > > in > > doDiv (x`div`(b^l)) l > > Best regards > Thorkil > > > On Monday 11 February 2008 07:15, 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! > > > > Uwe > > > > > 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@haskell.org > > http://www.haskell.org/mailman/listinfo/haskell-cafe > > > ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] fast integer base-2 log function?
Hello, If the standard libraries provide such a function, I haven't found it. I must admit that I haven't studied your code in detail. I usually do as follows for integer logarithms, shamelessly stolen from the Haskell report: > -- Integer log base (c.f. Haskell report 14.4): > > imLog :: Integer->Integer->Integer > imLog b x > = if x < b then > 0 > else > let > l = 2 * imLog (b*b) x > doDiv x l = if x < b then l else doDiv (x`div`b) (l+1) > in > doDiv (x`div`(b^l)) l Best regards Thorkil On Monday 11 February 2008 07:15, 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! > > Uwe > > > 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@haskell.org > http://www.haskell.org/mailman/listinfo/haskell-cafe > ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] fast integer base-2 log function?
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 .&. 0x /= 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@haskell.org
> http://www.haskell.org/mailman/listinfo/haskell-cafe
signature.asc
Description: Digital signature
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] fast integer base-2 log function?
Hi, all, a few days ago I had asked about fast integer log-2 routines, and got a couple of answers. I've now had a chance to play with the routines, and here's what I found. Initially, Thorkil's routine taken from the Haskell report was about 30% or so faster than mine. When I replaced the calls to my routine "powi" with native calls to ^, my routine became about 10% faster than Thorkil's routine. (I undoubtedly had some reason for using my own version of powi, but I have no idea anymore what that reason was... :-/ ) I initially thought that that speed difference might be due to the fact that my routine had base 2 hard-wired, whereas his routine is for general bases, but that seems not to be the case: when I modified my version to also do general bases, it stayed pretty much the same. I didn't do enough statistics-gathering to really be absolutely positively certain that my routine is indeed 10.000% faster, but there did seem to be a slight edge in speed there. Here's the latest version, in case anyone's interested. I had previously had effectively a bit-length version; since I made it general base-b I changed it to a log function. > ilogb :: Integer -> Integer -> Integer > ilogb b n | n < 0 = ilogb b (- n) > | n < b = 0 > | otherwise = (up b n 1) - 1 > where up b n a = if n < (b ^ a) > then bin b (quot a 2) a > else up b n (2*a) > bin b lo hi = if (hi - lo) <= 1 > then hi > else let av = quot (lo + hi) 2 > in if n < (b ^ av) > then bin b lo av > else bin b av hi Stefan's routine is, as expected, much much faster still: I tested the first two routines on numbers with 5 million or so bits and they took ~20 seconds of CPU time, whereas I tested Stefan's routine with numbers with 50 million bits, and it took ~11 seconds of CPU time. The limitation of Stefan's routine is of course that it's limited to base 2 -- it is truly a bit-length routine -- and I guess another potential limitation is that it uses GHC extensions, not pure Haskell (at least I had to compile it with -fglasgow-exts). But it's the speed king if those limitations aren't a problem! Uwe ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] fast integer base-2 log function?
On Thu, Feb 14, 2008 at 08:23:29PM -0800, Uwe Hollerbach wrote: > Stefan's routine is, as expected, much much faster still: I tested the > first two routines on numbers with 5 million or so bits and they took > ~20 seconds of CPU time, whereas I tested Stefan's routine with > numbers with 50 million bits, and it took ~11 seconds of CPU time. The > limitation of Stefan's routine is of course that it's limited to base > 2 -- it is truly a bit-length routine -- and I guess another potential > limitation is that it uses GHC extensions, not pure Haskell (at least > I had to compile it with -fglasgow-exts). But it's the speed king if > those limitations aren't a problem! At the risk of stating the obvious, it also hard-codes 32 bit words and certain configurable (but nobody bothers) details of the GMP data representation (big endian, 0 nails). Stefan signature.asc Description: Digital signature ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] fast integer base-2 log function?
Yes, I suspect you are right. I didn't look into that in much detail, although I did try exchanging "(2 ^ 5000)" with "(1 `shiftL` 5000)"; but that didn't make any difference. Uwe On Fri, Feb 15, 2008 at 9:21 AM, Ryan Ingram <[EMAIL PROTECTED]> wrote: > On Thu, Feb 14, 2008 at 8:23 PM, Uwe Hollerbach <[EMAIL PROTECTED]> > wrote: > > Stefan's routine is, as expected, much much faster still: I tested the > > first two routines on numbers with 5 million or so bits and they took > > ~20 seconds of CPU time, whereas I tested Stefan's routine with > > numbers with 50 million bits, and it took ~11 seconds of CPU time. > > This seems wrong to me; that routine should take a small constant > amount of time. I suspect you are measuring the time to construct the > 50-million bit numbers as well. If you constructed a single number > and called this routine on it several times I am sure you would get > far different results, with the first routines taking ~7-11s each and > Stefan's GHC/GMP-magic taking almost nothing. > > -- ryan > ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] fast integer base-2 log function?
On Thu, Feb 14, 2008 at 8:23 PM, Uwe Hollerbach <[EMAIL PROTECTED]> wrote: > Stefan's routine is, as expected, much much faster still: I tested the > first two routines on numbers with 5 million or so bits and they took > ~20 seconds of CPU time, whereas I tested Stefan's routine with > numbers with 50 million bits, and it took ~11 seconds of CPU time. This seems wrong to me; that routine should take a small constant amount of time. I suspect you are measuring the time to construct the 50-million bit numbers as well. If you constructed a single number and called this routine on it several times I am sure you would get far different results, with the first routines taking ~7-11s each and Stefan's GHC/GMP-magic taking almost nothing. -- ryan ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] fast integer base-2 log function?
> {-# 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 .&. 0x /= 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)
>
>
I tried the above but I got wrong results on 2^31..2^32-1 because in the
additions in ilog2i4, the number x was -1 because of sign extension
performed by the shifts all the way (thanks for the ghci debugger). (So,
yes, I do care somewhat about garbage on negatives :)
I modified to the following hoping also to use both on 32 and 64 bit
machines. Have I shot myself in the foot anyway? For example, is there a
guarantee that the most significant limb is non-zero? Is there any
possibility of this or some related function being added to Data.Bits?
{-# LANGUAGE MagicHash #-}
import GHC.Exts
import Data.Bits
limbSize = bitSize (0 :: Int)
ilog2 k x = case limbSize of
64 -> ilog2i0 k (fromIntegral x)
32 -> ilog2i1 k (fromIntegral x)
_ -> undefined
-- experiment with using a LUT here (hint: FFI + static arrays in C)
ilog2i0, ilog2i1, ilog2i2, ilog2i3, ilog2i4, ilog2i5 :: Int -> Word -> Int
ilog2i0 k x | x .&. 0x /= 0 = ilog2i1 (k + 32) (x `shiftR`
32)
| otherwise = ilog2i1 k x
ilog2i1 k x | x .&. 0x /= 0 = ilog2i2 (k + 16) (x `shiftR` 16)
| otherwise = ilog2i2 k x
ilog2i2 k x | x .&. 0xFF00 /= 0 = ilog2i3 (k + 8) (x `shiftR` 8)
| otherwise = ilog2i3 k x
ilog2i3 k x | x .&. 0xF0 /= 0 = ilog2i4 (k + 4) (x `shiftR` 4)
| otherwise = ilog2i4 k x
ilog2i4 k x | x .&. 0xC /= 0= ilog2i5 (k + 2) (x `shiftR` 2)
| otherwise = ilog2i5 k x
ilog2i5 k x | x .&. 0x2 /= 0= k + 1 + fromIntegral (x `shiftR` 1)
| otherwise = k + fromIntegral x
log2i :: Integer -> Int -- actually returns bit length
log2i (J# len adr) = ilog2 0 (I# (indexIntArray# adr (len -# 1#))) + I# (32#
*# (len -# 1#))
log2i (S# sml) = ilog2 0 (I# sml)
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] fast integer base-2 log function?
On Tue, Feb 26, 2008 at 09:33:57PM +, Jens Blanck wrote:
> > {-# 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 .&. 0x /= 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)
> >
> >
> I tried the above but I got wrong results on 2^31..2^32-1 because in the
> additions in ilog2i4, the number x was -1 because of sign extension
> performed by the shifts all the way (thanks for the ghci debugger). (So,
> yes, I do care somewhat about garbage on negatives :)
This is what you get for only testing on 100 and 2^34, I guess :)
If you change all the Int to Word (unsigned) it should work. Should.
> I modified to the following hoping also to use both on 32 and 64 bit
> machines. Have I shot myself in the foot anyway? For example, is there a
> guarantee that the most significant limb is non-zero? Is there any
> possibility of this or some related function being added to Data.Bits?
> [snip code]
It's still not going to be portable because I'm hardwiring the GMP "nail
count" parameter to 0. As for going standard - if you want this,
propose it! I can't think of a sane implementation of Integer that
doesn't support some kind of approximate logarithm.
Stefan
signature.asc
Description: Digital signature
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe
