[Haskell-cafe] PI and series
jerzy.karczmarczuk at info.unicaen.fr
jerzy.karczmarczuk at info.unicaen.fr
Wed Oct 10 12:47:36 EDT 2007
Dan Piponi writes:
> Jules Bean said:
>> If it is true of many Floating instances that (atan 1 * 4) is an
>> accurate way to calculate pi (and it appears to be 'accurate enough' for
>> Float and Double, on my computer) then adding it as a default doesn't
>> appear to do any harm.
>
> ... For the case of power series as an instance of Num, using
> 4*atan 1 gives me the wrong thing as it triggers an infinite
> summation, whereas I'd want pi to simply equal the constant power
> series. Now you could counter that by saying that power series are an
> esoteric case. But apart from code in the libraries, most code I've
> seen that provides an instance of Num is doing something mildly
> esoteric.
For me series are not esoteric. But they can be implemented in a more
cautious way, my own package behaves rationally:
Prelude> :l Series
[1 of 1] Compiling Series ( Series.hs, interpreted )
-- .... nasty messages omitted.
Ok, modules loaded: Series.
*Series> 4*atan 1
3.141592653589793
*Series> 4*atan 1 :: Series Double
Loading package haskell98 ... linking ... done.
[3.141592653589793,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, ...]
*Series> pi
3.141592653589793
*Series> pi :: Series Double
[3.141592653589793,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, ...]
*Series>
The implementation is, with (:-) as the cons-series operator, colon-dash
for those whose browsers will convert it to some nice smiley...:
zeros = 0:-zeros
...
instance Floating a => Floating (Series a)
where
pi = pi:-zeros
...
atan u@(u0:-_) = sInt (atan u0) (sDif u / (1+u*u))
etc. Of course, sDif differentiates and sInt integrates the series, see
2nd volume of Knuth for the algorithms, and papers, either a nice
tutorial of Doug McIlroy, or my old paper for the implementation.
The instance does it all. Default *could* help, and would give the same
answer, but it would be much slower, so I would override it anyway.
Jerzy Karczmarczuk
More information about the Haskell-Cafe
mailing list