[Haskell] for large x, log (x::Integer) :: Double

[Axel Simon]

On Sun, Jul 04, 2004 at 12:32:39PM +0200, Dylan Thurston wrote:
> On Wed, Jun 30, 2004 at 03:07:00PM -0700, Greg Buchholz wrote:
> > -- Inspired from Mr. Howard Oakley.  Might not qualify as "good",
> > -- but with this function I get log10(x)=849.114419903382
> > ... 
> 
> For those who aren't aware: working with logs base 2 internally will
> be very much faster than logs base 10, since the numbers are stored
> internally in a base-2 representation.  (Note that 'show' converts to
> base 10, which involves a large number of divisions in the easy
> algorithm.)

You could actually have a constant-time log function for Integer by 
looking at the internal representation. A gmp integer consists of an array 
of "limbs" (usually the size of your machine's register), the used bytes 
"size" within that allocated space and the sign bit in the MSB of "size". 
All integers are thus stored as unsigned quantities, hence it should be 
sufficient to look at size and the highest limb in the array to determine 
how many bits the number uses. I once wrote a marshaling function from 
Integer to C gmp integers which I attach. It should help as starting 
point. (The conversion from C to Haskell might be flawed, I always had 
problems reading them back in and never found out where. So beware).

Axel.

Run through hsc2hs:
                                                                                
import GHC.Base
import GHC.IOBase
import Foreign
import Foreign.C
 
-- A dummy data type for the external mpz_struct.
data MpzStruct = MpzStruct
                                                                                
-- To use the guts of an integer in C land, we need to copy the data of the
-- number from the GCed Haskell land to safer C land.
withInteger :: Integer -> (Ptr MpzStruct -> IO a) -> IO a
withInteger (J## size## barr##) act = do
  -- Make space for a mpz_t struct.
  allocaBytes #{const sizeof(__mpz_struct)} $ \mpzPtr -> do
    let bytesInLimbs = I## (sizeofByteArray## barr##)
    let limbs@(I## limbs##) = bytesInLimbs `div` #{const sizeof(mp_limb_t)}
    -- Make space for a limb as big as the one of the Integer.
    allocaBytes bytesInLimbs $ \limbPtr@(Ptr addr##) -> do
      -- Copy the limb of the Integer to the newly created mpz_t.
      let
        copy## :: Int## -> State## RealWorld -> (## State## RealWorld, () ##)
        copy## 0## s = (## s, () ##)
        copy## i## s =
          case i## -## 1## of { i'## ->
            case indexWord#{const 8*sizeof(mp_limb_t)}Array##
              barr## i'## of { val## ->
              case writeWord#{const 8*sizeof(mp_limb_t)}OffAddr##
                addr## i'## val## s of { s' ->
                copy## i'## s'
              }
            }
          }
      IO $ \s -> copy## limbs## s
      -- Fill the structure
      #{poke __mpz_struct, _mp_alloc} (castPtr mpzPtr) limbs
      #{poke __mpz_struct, _mp_size} (castPtr mpzPtr) (I## size##)
      #{poke __mpz_struct, _mp_d} (castPtr mpzPtr) limbPtr
      act mpzPtr
withInteger small act = withInteger (toBig small) act
                                                                                
-- To transfer an integer to Haskell land we allocate an mpz_t
-- structure, initialize it, send it to the C function,
-- then copy its content to Haskell and finally free it.
extractInteger :: (Ptr MpzStruct -> IO a) -> IO (Integer, a)
extractInteger act =
  -- Make space for a mpz_t struct.
  allocaBytes #{const sizeof(__mpz_struct)} $ \mpzPtr ->
    bracket_ (mpz_init mpzPtr) (mpz_clear mpzPtr) $ do
      res <- act mpzPtr
      -- Extract the contents of the mpz_t to a Haskell integer.
      (I## alloc##) <- #{peek __mpz_struct, _mp_alloc} (castPtr mpzPtr)
      (I## size##) <- #{peek __mpz_struct, _mp_size} (castPtr mpzPtr)
      (Ptr addr##) <- #{peek __mpz_struct, _mp_d} (castPtr mpzPtr)
      val <- IO $ \s## ->
        case newByteArray## (alloc## *## #{const sizeof(mp_limb_t)}##) s## of
          { (## s1##, mutarr## ##) ->
            let
              copy## :: Int## -> State## RealWorld -> State## RealWorld
              copy## 0## s## = s##
              copy## i## s##=
                let
                  i'##  = i## -## 1##
                  val## = indexWord#{const 8*sizeof(mp_limb_t)}OffAddr##
                            addr## i'##
                  s'##  = writeWord#{const 8*sizeof(mp_limb_t)}Array##
                            mutarr## i'## val## s##
                in copy## i'## s'##             in
              case copy## alloc## s1## of { s2## ->
                case unsafeFreezeByteArray## mutarr## s2## of {
                  (## s3##, barr## ##) -> (## s3##, J## size## barr## ##)
                }
              }
           }
      return (val, res)
                                                                                
foreign import ccall unsafe "__gmpz_init" mpz_init :: Ptr MpzStruct -> IO ()
foreign import ccall unsafe "__gmpz_clear" mpz_clear :: Ptr MpzStruct -> IO ()