profiling problems

[Damien R. Sullivan]

On Tue, Feb 25, 2003 at 06:08:42PM -0800, Hal Daume III wrote:
> Have you removed the .o and .hi files before compiling with -prof?  These
> are not compatible.

Well, I got this code to run, although I thought I couldn't earlier from
clean.  But the longer code I'm really concerned with still bus errors.
Compiled with profiling on a Linux box it seems to run okay, so perhaps it's a
profiling-Sun interaction?


import Numeric
import Ratio
import System
-- less general; wants n to be power of 2, for faster logs
-- version 2: memoing
-- version 3: trying out sophisticated calculations

-- TODO cache comb; try calculating H(n-1) - ln n directly; try grouping terms
-- so fewer GCD calls are made; try figuring series for top and bottom
-- directly so no GCD calls

logtop 1 = 1
logtop n = 9*n*logtop (n-2) + oddfact (n-2)

loddfact :: [Integer]
loddfact = 1 : zipWith (*) [3,5 .. ] loddfact

oddfact :: Integer -> Integer
oddfact n = loddfact!!ni
    where 
    ni = nint `div` 2
    nint = fromInteger n

log_two :: Int -> Rational
log_two lim = 2* logtop n % (3^lim * oddfact n)
    where n = fromIntegral lim

log_tup :: Int -> Int -> (Rational, Rational)
log_tup lim b =
    (l + br*lo_err, l + br*hi_err)
    where 
    br = toRational b
    l = br*log_two lim 
    -- TODO double check these
    lo_err = (2)/(n+1)*(d/(1+c))^(lim+1)
    hi_err = (2)/(n+1)*d^(lim+1)
-- original logs
--     lo_err = (-2)/(n+1)*(d/(1+c))^(lim+1)
--     hi_err = -lo_err
    -- lo_err = (-1)/(n+1)*x^(lim+1)
    -- c = 0
    n = toRational lim
    d = 1%3
    c = d

floatRat :: Rational -> Int -> Integer
floatRat r lim =
    num*10^lim `div` den
    where
    num = numerator r
    den = denominator r

bernoulli :: Integer -> Rational
bernoulli n = btab!!ni
    where
    ni = fromInteger n
    btab :: [Rational]
    -- TODO use zip bhelp [2..]; should be clearer, but check output
    -- might not work actually, With gives the function...
    btab = 1 : -1/2 : zipWith bhelp [2 .. ] btab
    bhelp n _ | odd n = 0
              | otherwise = -1/fromIntegral(n+1) *sumbn n

lfact :: [Integer]
lfact = 1 : zipWith (*) [1 .. ] lfact

fact :: Integer -> Integer
fact 0 = 1
fact n = lfact!!ni
    where 
    ni = fromInteger n
    
-- 22% time spent here; improve?  Probably needs smarts, not memo.  Not clear
-- if smarts would speed it up now, test separately
-- yeah, memo worn't work; comb (n+1) i is only called once
-- comb m n = comb m (n-1) * (m-n+1)/n ?  Check.  Guess this does call for
-- memo.
comb :: Integer -> Integer -> Integer
comb m 0 = 1
comb m n = (numerator(toRational (fact m) / toRational (fact n * fact (m-n))))
         
sumbn :: Integer -> Rational
sumbn n = sum [ bernoulli i * fromIntegral(comb (n+1) i) | i <- [0 .. n-1] ]

harmonic :: Integer -> Rational
harmonic n = harmtop n % fact n

harmtop 1 = 1
harmtop n = n * harmtop (n-1) + fact (n-1)  

-- assuming one bound of the integral is zero, else use the abs line.  Doesn't
-- seem to make too much difference to accuracy.
gamma_tup :: Int -> Integer -> Int -> (Rational, Rational)
gamma_tup lim m b = 
    (base - log_hi - integ_hi, base - log_lo - integ_lo)
--    (base - log_hi - abs integ, base - log_lo + abs integ)
    where
    base = harmonic (n-1) - bigsum m n
    (log_lo, log_hi) = log_tup lim b
    n = 2^b
    -- bern of 4k is negative; bern of 4k+2 is positive; abs value of integral
    -- is capped by Bm/(m*n^m)
    integ_hi = if m `mod` 4 == 0 then 0 else integ
    integ_lo = if m `mod` 4 == 0 then integ else 0
    integ = integral m n

bigsum :: Integer -> Integer -> Rational
bigsum m n = sum [(-1)^(i-1)*bernoulli i /fromInteger (i*n^i) | i <- 1:[2,4..m] ]

integral :: Integer -> Integer -> Rational
integral m n = bernoulli m/toRational m / toRational n^^m

main :: IO ()
main = 
     do
        args <- System.getArgs
        let (m, b) = (read (args!!0), read (args!!1))
        let lim :: Int
            lim = read (args!!2)
        let (glo, ghi) = gamma_tup lim m b
        putStrLn (show glo)
        putStrLn (show ghi)
        putStrLn (show (floatRat glo lim))
        putStrLn (show (floatRat ghi lim))
        putStrLn (show (floatRat ghi lim - floatRat glo lim))