Hugs Humor

[Jerzy Karczmarczuk]

Jon Fairbairn comments //Steffen Mazanek//:

>>Prelude> 0.1::Rational
>>13421773 % 134217728
>>Prelude> 13421773/134217728
>>0.1
>>
>>I do not know how this fraction is calculated, but
>>it does not fit my expectations :-)
>  
> Remember that internally arithmetic is binary, and that 0.1
> can't be expressed exactly as a floating point number. I
> think that's the explanation.
> 
>>Ok, ok, it is no bug...
> 
> No, I think it is a bug: 0.1 ought to be equivalent to
> fromRational (1%10), but Hugs isn't intended for numerical
> work. GHCi gets the right answer.


This is less a bug than a Nessie monster which haunts Hugs
some centuries already, and on Internet the issue has been
discussed at least 4 times. The old, experimental Gofer
Prelude numeric functions were sometimes abominable, since
Mark Jones concentrated on other things, and nobody really
complained, people were busy with other stuff as well.

But I can't understand why this continues until now. The
obvious technique to convert floats to rationals is the
continued fraction expansion which gives the simplified
answer fast.

I don't understand the remark that the internal arithmetic is
binary. Sure, it is, so what? Why Ross Paterson underlies
this as well? He concludes:

 > The real fix would be to keep the literals as Rationals, but
 > this would be too expensive in the Hugs setting.

Andrew Bromage says some words about errors and representation.
I think that the problem can (perhaps should) be dealt with at
a higher level. What's wrong with conversion functions like
those below. First, convert a float to a lazy list of coeffs.
of a regular continued fraction:

tocfrac x =
  let n = floor x
      y = x-fromInteger n
  in n : if y==0.0 then [] else tocfrac (recip y)


and then reconstruct using Euler sequences as described in
Knuth or perhaps an optimized method, the one I cite is not
very efficient. It gives a list of (N,D) -- *all* rational
approximants of the original float.


continuant l@(_:ql) = zip (tail (euseq l)) (euseq ql)
  where
    euseq [] = []
    euseq (x:q) = eus 1 x q
    eus p0 p1 (a:cq) = p0 : eus p1 (a*p1+p0) cq
    eus p0 p1 [] = [p0,p1]



Now, test it:

pp=3.141592653589793
r=take 10 (continuant (tocfrac pp))

You should get

[(3,1),(22,7),(333,106),(355,113),(103993,33102),(104348,33215), ...
etc;, anyway all that is already inexact...

For 0.1 one gets [(0,1),(1,10)]

Jerzy Karczmarczuk