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/.
Very nice. However, I'm intrigued by the following remark:
ReplyDelete>> 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?
First, what type would you have picked for literals that look like floating point (like 1.2e-3)?
ReplyDeleteYou 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!