Wednesday, April 25, 2007

Overloading Haskell numbers, part 3, Fixed Precision. Fixed precision numbers can be handy to have sometimes. Haskell only provides integral, rational, and floating point in the Prelude, but I want more! I want to define a type of fixed precision numbers, i.e., numbers that have a certain number of decimals. Rational numbers are handy to represent these numbers so here we go.
newtype Fixed = F Rational deriving (Eq, Ord, Enum, Num, Fractional, Real, RealFrac)
This uses an GHC extension that allows you to derive anything from the representation type so we don't have to define all those tedious instances. But wait, where is the fixed precision? What I just defined is just the Rational numbers with a different name. I really want a type which is parametrized over what precision it has. OK, so here we go again:
newtype Fixed e = F Rational deriving (Eq, Ord, Enum, Num, Fractional, Real, RealFrac)
Here e is meant to be a type that specifies the precision. Let's capture those types that we allow as the precision as a type class:
class Epsilon e where
    eps :: e -> Rational
The idea here is that the actual argument to the eps function is not used, it's just a proxy and the type carries all the information. Here's an example:
data Eps1
instance Epsilon Eps1 where
    eps _ = 1
The eps returns 1 for the Eps1 type. So the intended meaning is that the type Fixed Eps1 are numbers within an epsilon of the right answer, and epsilon is 1. Similarly
data Prec50
instance Epsilon Prec50 where
    eps _ = 1e-50
The type Fixed Prec50 has an epsilon of 1e-50, i.e., 50 decimals. Btw, notice that what looks like a floating point literal is actually a rational number; one of the very clever decisions in the original Haskell design. We can also define type constructors like
data EpsDiv10 p
instance (Epsilon e) => Epsilon (EpsDiv10 e) where
    eps e = eps (un e) / 10
       where un :: EpsDiv10 e -> e
             un = undefined
The type Fixed (EpsDiv10 Prec50) has an epsilon of 1e-50/10 = 1e-51. Note how the function un is just a dummy to get the e. If we had used the scoped type variable extension we could have avoided this function. OK, so we can compute with these fixed precision numbers, but where does the epsilon enter? Nowhere, so far. All arithmetic is performed exactly by the rationals, but our first use of the epsilon is in the show function.
instance (Epsilon e) => Show (Fixed e) where
    showsPrec = showSigned showFixed
      where showFixed f@(F x) = showString $ show q ++ "." ++ decimals r (eps $ un f)
              where q :: Integer
             (q, r) = properFraction x
      un :: Fixed e -> e
             un = undefined
             decimals a e | e >= 1 = ""
                          | otherwise = intToDigit b : decimals c (10 * e)
                 where (b, c) = properFraction (10 * a)
To print a fixed point number we first print the integer part. The properFraction method handily returns the integer part and some left over fraction that is <1. We then print decimal after decimal by multiplying the fraction by 10, converting the integer part to a digit, and then recursing. While generating these decimals we also multiply the epsilon by 10 each time and stop printing when epsilon is >= 1, because then we don't need any more decimals. All right, lets put it to the test:
Data.Number.Fixed> 1/3 :: Fixed Prec50
0.33333333333333333333333333333333333333333333333333
Data.Number.Fixed> 3/8 :: Fixed Prec10
0.3750000000
Data.Number.Fixed> (1 :: Fixed Prec10) + (0.5 :: Fixed Prec50)
    Couldn't match expected type `Prec10'
           against inferred type `Prec50'
      Expected type: Fixed Prec10
      Inferred type: Fixed Prec50
As we want, we can't mix different fixed point numbers with different precisions. Sometimes we do want to convert, so let's provide such a function
convertFixed :: Fixed e -> Fixed f
convertFixed (F x) = F x
It looks like this function does nothing, and it's true. It just converts the type, and since the epsilon type is not involved in the number itself it's trivial. Well, basic arithmetic isn't that much fun. What about the transcendental functions? Luckily Jan Skibinski has already implemented all of them. They compute the various function to with an epsilon using continued fractions. The instance declaration looks something like:
instance (Epsilon e) => Floating (Fixed e) where
    sqrt = toFixed1 F.sqrt
    exp = toFixed1 F.exp
    ...

toFixed1 :: (Epsilon e) => (Rational -> Rational -> Rational) -> Fixed e -> Fixed e
toFixed1 f x@(F r) = F (f (eps (un x)) r)
  where un :: Fixed e -> e
        un = undefined
And all the clever code I stole from Jan is in the module F, see link below. Now it's getting more fun.
Data.Number.Fixed> pi :: Fixed Prec50
3.14159265358979323846264338327950288419716939937510
Data.Number.Fixed> exp 10 :: Fixed Prec10
22026.4657881843
Data.Number.Fixed> log 2 :: Fixed Prec50
0.69314718055994530941723212145817656807550013436026
Data.Number.Fixed> exp (sin (sqrt (54/42))) :: Fixed Prec10
2.4745696343
Data.Number.Fixed> sqrt 2 :: Fixed Prec500
1.41421356237309504880168872420969807856967187537694807317667973799073247
8462107038850387534327641572735013846230912297024924836055850737212644121
4970999358314132226659275055927557999505011527820605714701095599716059702
7453459686201472851741864088919860955232923048430871432145083976260362799
5251407989687253396546331808829640620615258352395054745750287759961729835
5752203375318570113543746034084988471603868999706990048150305440277903164
5424782306849293691862158057846311159666871301301561856898723724

Data.Number.Fixed Data.Complex> let i = 0 :+ 1 :: Complex (Fixed Prec50) 
Data.Number.Fixed Data.Complex> exp (i*pi)
(-0.99999999999999999999999999999999999999999999999999) :+ 0.00000000000000000000000000000000000000000000000000
Naturally, we need to package up all these definitions in a module that will keep the Fixed type abstract, and you can also keep the Epsilon class and the various Epsilon types abstract if you provide enough building blocks. Source at http://www.augustsson.net/Darcs/Data.Number/.

2 comments:

  1. Very nice. However, I'm intrigued by the following remark:
    >> Btw, notice that what looks like a floating point literal
    is actually a rational number; one of the very clever decisions in the
    original Haskell design. <<

    This has always seemed a slightly strange design decision to me. For example, in GHC, I can do:
    > toRational 0.6
    5404319552844595%9007199254740992
    which is unexpected and not obviously useful.

    What's the benefit of being able to convert reals to rationals?

    ReplyDelete
  2. First, what type would you have picked for literals that look like floating point (like 1.2e-3)?
    You can't use Float or Double, because that limits the precision at which they are stored. For instance, with my fixed point library I can compute, say, 2*1e-500 and get the right answer. But 1e-500 would be 0 if it was interpreted as a double. So to be able to have floating point literals of arbitrary size and precision we must pick a type that can represent those. And a suitable type in Haskell is Rational.

    Doing something like (toRational (0.6::Double)) is not that useful, I agree. But it is not *not* converting a real number to a rational number. It's converting a floating point number to a rational number. The numbers we can represent with floating point is a finite subset of the rationals. So it makes sense that we can take a floating point number and convert it to the rational number that it actually represents.

    Here's a fun one:

    Prelude> toRational 0.1 - 0.1
    1%180143985094819840

    Show's you why 0.1 is a dangerous thing when using floating point!

    ReplyDelete