possible bug

[Andrew J Bromage]

G'day.

On Sun, Jun 22, 2003 at 02:26:00PM -0700, John Mann wrote:

> Now, if I change the number of terms to be about 10
> times as many, by changing the program to:
[...]
> I get the LESS-accurate result:

Your problem is entirely due to floating point round-off error.

I can think of two solutions to your problem off the top of my head.
The most idiomatic Haskell solution is to work in rationals and convert
to Double at the end.  This also gives you the most accurate answer
because there is no rounding error (until the last minute, of course).

The most "numeric analyst-like" solution is to be more careful about
the order that you add the terms.  This, for example, is a simple
modification of your code which always ensure that the two smallest
numbers-to-be-added are added first:

	import List

	-- This is like your "pie" function with two modifications:
	--     1) It returns a list of terms
	--     2) It doesn't have the special case for the last
	--        term.  This is dealt with in the case expression
	--        instead.

	pieSeries :: Integer -> [Double]
	pieSeries n
	    = case ps' n of
		(t:ts) -> t * 0.5 : ts
		[]     -> []
	    where
		ps' n
		    | n==1 = [4]
		    | otherwise =  4/fromInteger n :
				  -4/(fromInteger n - 2) :
				  ps' (n-4)

	-- Note that we need to sort by absolute value.  The
	-- simplest way is to use the Schwartzian transform.
	-- (See http://haskell.org/hawiki/SchwartzianTransform)

	pie :: Integer -> Double
	pie n
	    = sumTerms (sort [ (abs x, x) | x <- pieSeries n ])
	    where
		sumTerms [(_,t)] = t
		sumTerms ((_,t1):(_,t2):ts)
		    = let t = t1+t2
		      in sumTerms (insert (abs t, t) ts)

Cheers,
Andrew Bromage