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))