profiling problems
Damien R. Sullivan
dasulliv@cs.indiana.edu
Tue, 25 Feb 2003 23:48:34 -0500
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))