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

On Tue, Feb 25, 2003 at 09:12:23PM -0800, Hal Daume III wrote:

My guess is that the error is here. There's a bug in ghc where (x `div` 0) and (x `mod` 0) will results in a seg fault. For now (according to one of the Simons, this is fixed in the next version) you should do something

Doesn't help. I didn't expect it to -- my code shouldn't be divving by zero, and anyway the bus error happens *immediately* on startup, whereas it wouldn't call floatRat for half a minute if it were working.

If that doesn't fix the problem, perhaps you can try to whittle the example down to something shorter.

The first code I posted was my shorter attempt. At any rate, it's hard to pull stuff out of this and still have it work. And I've had versions which did compile with profiling when the version before or after wouldn't. This one doesn't, on my Sun. Does on the Linux box... -xx- Damien X-)