[Haskell-cafe] Re: Is logBase right?

[Daniel Fischer]

Am Sonntag 23 August 2009 15:27:48 schrieb Steve:
> On Sun, 2009-08-23 at 15:12 +0400, Eugene Kirpichov wrote:
> > > There is *not* the same problem in Python:
> > > $ python
> > > Python 2.6.2 (r262:71600, Jul  9 2009, 23:16:53)
> > > [GCC 4.4.0 20090506 (Red Hat 4.4.0-4)] on linux2
> > > Type "help", "copyright", "credits" or "license" for more information.
> > >
> > >>>> import math
> > >>>> math.log10(1000)
> > >
> > > 3.0
> > >
> > >>> import math
> > >>> math.log(1000,10)
> >
> > 2.9999999999999996
>
> That is surprising. I think its a Python bug - the results should be
> consistent.

On the other hand, it is also desirable to have log(a,b) == log(a)/log(b).
If you have a special implementation of log10, you must decide which consistency you want 
to break, log(a,10) == log10(a) or log(a,10) == log(a)/log(10).

>
> > > Recent work in Python 3 (and Python 2.6) has improved the handling of
> > > floating point numbers, and addresses exactly the problem that Roberto
> > > has raised.
> > >
> > > I see no reason why Haskell could not improve its handling of floating
> > > point numbers by using similar techniques.
> >
> > You mean introducing a "log10" function into the definition of the
> > Floating class ? That might be a proposal for Haskell Prime.
>
> I was not thinking of the log10 function. I was thinking of the changes
> mentioned in
> http://docs.python.org/3.1/whatsnew/3.1.html
> where it says
> "Python now uses David Gay’s algorithm for finding the shortest floating
> point representation that doesn’t change its value. This should help
> mitigate some of the confusion surrounding binary floating point
> numbers."

Note however that that only concerns the output, not the underlying value.
It doesn't change the fact that log(1000)/log(10) (== log(1000,10)) < 3 in Python, too.

>
> Also, I had a problem using floating point in Python where
>
> >>> round(697.04157958254996, 10)
>
> gave
> 697.04157958259998
> which is closer to
> 697.0415795826
> than the desired result of
> 697.0415795825

Well, 10^10*(697.04157958254996 :: Double) == 6970415795825.5,
that gets rounded to 6970415795826 under both, round-half-up and round-half-to-even, and 
(6970415795826/10^10 :: Double) == 697.04157958259998

Though ghci displays it as 697.0415795826, unlike Python.
>
> Its been fixed in the latest versions of Python:
> >>> round(697.04157958254996, 10)
>
> 697.0415795825
>
>
> I'm not sure what the equivalent in Haskell is. Is there a function for
> rounding to a number of decimal digits ?
> I came up with this:
>
> roundN :: Double -> Int -> Double
> roundN n ndigits = fromIntegral (round $ n * m) / m
>   where m = 10 ^ ndigits
>
> ghci> roundN 697.04157958254996 10
> 697.0415795826
>
> which is not the desired result.

Prelude Numeric> let roundN :: Double -> Int -> Double; roundN x d = let { m = 10^(d+1); i 
= floor (m*x); r = i `mod` 10; j = if r < 5 then i-r else i-r+10 } in fromInteger j/m

roundN :: Double -> Int -> Double
roundN x d = fromInteger j / m
      where
        m = 10^(d+1)
        i = floor (m*x)
        r = i `mod` 10
        j | r < 5       = i-r
          | otherwise   = i+10-r

ghci> roundN 697.04157958254996 10
697.0415795825

I haven't treated negative input, but, more importantly, floating point representation 
being what it is, there are cases where this doesn't do what you want/expect either.

>
> Steve