Revised numerical prelude, version 0.02

Dylan Thurston dpt@math.harvard.edu
Tue, 13 Feb 2001 18:32:21 -0500


--dDRMvlgZJXvWKvBx
Content-Type: text/plain; charset=us-ascii
Content-Disposition: inline

Here's a revision of the numerical prelude.  Many thanks to all who
helped.  Changes include:

* Removed "Powerful", replacing it with (^) in Num and (**) in Real.
* Fixed numerous typos
* Removed gcd and co. from Integral
* Added shortcomings & limitation of scope
* Added SmallIntegral, SmallReal
* wrote skeleton VectorSpace, PowerSeries
* Added framework to make it run under hugs.  There are some usability issues.

Any comments welcome!

Best,
	Dylan Thurston

--dDRMvlgZJXvWKvBx
Content-Type: text/plain; charset=us-ascii
Content-Disposition: attachment; filename="NumPrelude.lhs"

Revisiting the Numeric Classes
------------------------------
The Prelude for Haskell 98 offers a well-considered set of numeric
classes which cover the standard numeric types (Integer, Int,
Rational, Float, Double, Complex) quite well.  But they offer limited
extensibility and have a few other flaws.  In this proposal we will
revisit these classes, addressing the following concerns:

(1) The current Prelude defines no semantics for the fundamental
    operations.  For instance, presumably addition should be
    associative (or come as close as feasible), but this is not
    mentioned anywhere.

(2) There are some superfluous superclasses.  For instance, Eq and
    Show are superclasses of Num.  Consider the data type

  data IntegerFunction a = IF (a -> Integer)

    One can reasonably define all the methods of Num for
    IntegerFunction a (satisfying good semantics), but it is
    impossible to define non-bottom instances of Eq and Show.

    In general, superclass relationship should indicate some semantic
    connection between the two classes.

(3) In a few cases, there is a mix of semantic operations and
    representation-specific operations.  toInteger, toRational, and
    the various operations in RealFloating (decodeFloat, ...) are the
    main examples.

(4) In some cases, the hierarchy is not finely-grained enough:
    operations that are often defined independently are lumped
    together.  For instance, in a financial application one might want
    a type "Dollar", or in a graphics application one might want a
    type "Vector".  It is reasonable to add two Vectors or Dollars,
    but not, in general, reasonable to multiply them.  But the
    programmer is currently forced to define a method for (*) when she
    defines a method for (+).

In specifying the semantics of type classes, I will state laws as
follows:
  (a + b) + c === a + (b + c)
The intended meaning is extensional equality: the rest of the program
should behave in the same way if one side is replaced with the
other.  Unfortunately, the laws are frequently violated by standard
instances; the law above, for instance, fails for Float:

  (100000000000000000000.0 + (-100000000000000000000.0)) + 1.0 = 1.0
  100000000000000000000.0 + ((-100000000000000000000.0) + 1.0) = 0.0

Thus these laws should be interpreted as guidelines rather than
absolute rules.  In particular, the compiler is not allowed to use
them.  Unless stated otherwise, default definitions should also be
taken as laws.

This version is fairly conservative.  I have retained the names for
classes with similar functions as far as possible, I have not made
some distinctions that could reasonably be made, and I have tried to
opt for simplicity over generality.

Thanks to Brian Boutel, Joe English, William Lee Irwin II, Marcin
Kowalczyk, and Ken Shan for helpful comments.

Scope & Limitations/TODO:
* It might be desireable to split Ord up into Poset and Ord (a total
 ordering).  This is not addressed here.

* In some cases, this heirarchy is not fine-grained enough.  For
  instance, time spans ("5 minutes") can be added to times ("12:34"),
  but two times are not addable.  ("12:34 + 8:23"??)  As it stands,
  users have to use a different operator for adding time spans to times
  than for adding two time spans.  Similar issues arise for vector space
  et al.  This is a consciously-made tradeoff, but might be changed.

  This becomes most serious when dealing with quantities with units
  like length/distance^2, for which (*) as defined here is useless,
  but Haskell's type system doesn't seem to be strong enough to deal
  with those in any convenient way.

  [One way to see the issue: should
    f x y = iterate (x *) y
  have principal type
    (Num a) => a -> a -> [a]
  or something like
    (Num a, Module a b) => a -> b -> [b]
  ?]

* I stuck with the Haskell 98 names.  In some cases I find them
  lacking.  Given free rein and not worrying about backwards
  compatibility, I might rename the classes as follows:
    Num           --> Ring
    Floating      --> Analytic
    RealFloat     --> RealAnalytic

* I'm not happy with Haskell's current treatment of numeric literals.
  I'm particularly unhappy with their use in pattern matching.  I feel
  like it should be a special case of some more general construction.
  I'd like to make it easier to use a non-standard Prelude, but
  there's a little too much magic.  For instance, the definition of
  round in the Haskell 98 Prelude is

    round x          =  let (n,r) = properFraction x
                             m     = if r < 0 then n - 1 else n + 1
                        in case signum (abs r - 0.5) of
                              -1 -> n
                              0  -> if even n then n else m
                              1  -> m

  I'd like to copy this over to this revised library.  But the numeric
  constants have to be wrapped in explicit calls to fromInteger.
  Worse, the case statement must be rewritten!

> module NumPrelude where
> import qualified Prelude as P
> -- Import some standard Prelude types verbatim verbandum
> import Prelude hiding (
>        Int, Integer, Float, Double, Rational, Num(..), Real(..),
>        Integral(..), Fractional(..), Floating(..), RealFrac(..),
>        RealFloat(..), subtract, even, odd,
>        gcd, lcm, (^), (^^))
>		 
>
> infixr 8  ^, **
> infixl 7  *
> infixl 7 /, `quot`, `rem`, `div`, `mod`
> infixl 6  +, -
>
> class Additive a where
>     (+), (-) :: a -> a -> a
>     negate   :: a -> a
>     zero     :: a
>
>      -- Minimal definition: (+), zero, and (negate or (-1))
>     negate a = zero - a
>     a - b    = a + (negate b)

Additive a encapsulates the notion of a commutative group, specified
by the following laws:

          a + b === b + a
   (a + b) + c) === a + (b + c)
       zero + a === a
 a + (negate a) === 0

Typical examples include integers, dollars, and vectors.

> class (Additive a) => Num a where
>     (*)         :: a -> a -> a
>     one	  :: a
>     fromInteger :: Integer -> a
>     (^)         :: (SmallIntegral b) => a -> b -> a
>
>       -- Minimal definition: (*), (one or fromInteger)
>     fromInteger n | n < 0 = negate (fromInteger (-n))
>     fromInteger n | n >= 0 = reduceRepeated (+) zero one n
>     a ^ n | n < zero = error "Illegal negative exponent"
>           | True  = reduceRepeated (*) one a (toInteger n)
>     one = fromInteger 1

Num encapsulates the mathematical structure of a (not necessarily
commutative) ring, with the laws

  a * (b * c) === (a * b) * c
      one * a === a
      a * one === a
  a * (b + c) === a * b + a * c

Typical examples include integers, matrices, and quaternions.

"reduceRepeat op a n" is an auxiliary function that, for an
associative operation "op", computes the same value as

  reduceRepeated op a0 a n = foldr op a0 (repeat (fromInteger n) a)

but applies "op" O(log n) times and works for large n.  A sample
implementation is below:

> reduceRepeated :: (a -> a -> a) -> a -> a -> Integer -> a
> reduceRepeated op a0 a n
>                   | n == 0 = a0
>                   | even n = reduceRepeated op a0 (op a a) (div n 2)
>                   | True   = reduceRepeated op (op a0 a) (op a a) (div n 2)

> class (Num a) => Integral a where
>     div, mod :: a -> a -> a
>     divMod :: a -> a -> (a,a)
>
>      -- Minimal definition: divMod or (div and mod)
>     div a b = let (d,_) = divMod a b in d
>     mod a b = let (_,m) = divMod a b in m
>     divMod a b = (div a b, mod a b)

Integral corresponds to a commutative ring, where "a mod b" picks a
canonical element of the equivalence class of "a" in the ideal
generated by "b".  Div and mod satisfy the laws

                        a * b === b * a
(a `div` b) * b + (a `mod` b) === a
              (a+k*b) `mod` b === a `mod` b
                    0 `mod` b === 0
                    a `mod` 0 === a

Typical examples of Integral include integers and polynomials over a
field.  Note that for a field, there is a canonical instance defined
by the above rules; e.g.,

  instance Integral Rational where
      divMod a 0 = (a,undefined)
      divMod a b = (0,a/b)

> class (Num a) => Fractional a where
>     (/)          :: a -> a -> a
>     recip        :: a -> a
>     fromRational :: Rational -> a
>
>      -- Minimal definition: recip or (/)
>     recip a = one / a
>     a / b = a * (recip b)
>     fromRational r = fromInteger (numerator r) / fromInteger (denominator r)
>     -- I'd like this next definition to be legal.
>     -- It would only apply if there were an implicit instance for Num a
>     -- through Fractional a.
>  -- a ^ n | n < 0 = reduceRepeated (^) one (recip a) (negate (toInteger n))
>  --       | True  = reduceRepeated (^) one a (toInteger n)

Fractional again corresponds to a commutative ring.  Division is
partially defined and satisfies

           a * b === b * a
   a * (recip a) === one

when it is defined.  To safely call division, the program most take
type-specific action; e.g., the following is appropriate in many
cases:

safeRecip :: (Integral a, Eq a, Fractional a) => a -> Maybe a
safeRecip a = let (q,r) = one `divMod` b in
    if (r == 0) then Just q else Nothing

Typical examples include rationals, the real numbers, and rational
functions (ratios of polynomials).  An instance should not typically
be declared unless most elements are invertible.

> -- Note: I think "Analytic" would be a better name than "Floating".
> class (Fractional a) => Floating a where
>     pi                  :: a
>     exp, log, sqrt      :: a -> a
>     logBase, (**)       :: a -> a -> a
>     sin, cos, tan       :: a -> a
>     asin, acos, atan    :: a -> a
>     sinh, cosh, tanh    :: a -> a
>     asinh, acosh, atanh :: a -> a
> 
>         -- Minimal complete definition:
>         --      pi, exp, log, sin, cos, sinh, cosh
>         --      asinh, acosh, atanh
>     x ** y           =  exp (log x * y)
>     logBase x y      =  log y / log x
>     sqrt x           =  x ** (fromRational 0.5)
>     tan  x           =  sin  x / cos  x
>     tanh x           =  sinh x / cosh x

Floating is the type of numbers supporting various analytic
functions.  Examples include real numbers, complex numbers, and
computable reals represented as a lazy list of rational
approximations.

Note the default declaration for a superclass.  See the comments
below, under "Instance declaractions for superclasses".

The semantics of these operations are rather ill-defined because of
branch cuts, etc.

> class (Num a, Ord a) => Real a where
>     abs    :: a -> a
>     signum :: a -> a
>
>       -- Minimal definition: nothing
>     abs x    = max x (negate x)
>     signum x = case compare x zero of GT -> one
>				        EQ -> zero
>				        LT -> negate one

This is the type of an ordered ring, satisfying the laws

             a * b === b * a
     a + (max b c) === max (a+b) (a+c)
  negate (max b c) === min (negate b) (negate c)
     a * (max b c) === max (a*b) (a*c) where a >= 0

Note that abs is in a rather different place than it is in the Haskell
98 Prelude.  In particular,

  abs :: Complex -> Complex

is not defined.  To me, this seems to have the wrong type anyway;
Complex.magnitude has the correct type.

> class (Real a, Floating a) => RealFrac a where
> -- lifted directly from Haskell 98 Prelude
>     properFraction   :: (Integral b) => a -> (b,a)
>     truncate, round  :: (Integral b) => a -> b
>     ceiling, floor   :: (Integral b) => a -> b
> 
>         -- Minimal complete definition:
>         --      properFraction
>     truncate x   =  m  where (m,_) = properFraction x
>     
>     round x      =  fromInteger (
>                     let (n,r) = properFraction x
>                         m     = if r < zero then n - one else n + one
>                       in case compare (abs r - (fromRational 0.5)) zero of
>                             LT -> n
>                             EQ  -> if even n then n else m
>                             GT  -> m
>                     )
>     
>     ceiling x      =  fromInteger (if r > zero then n + one else n)
>                       where (n,r) = properFraction x
>     
>     floor x        =  fromInteger (if r < zero then n - one else n)
>                       where (n,r) = properFraction x

As an aside, let me note the similarities between "properFraction x"
and "x divMod 1" (if that were defined).  In particular, it might make
sense to unify the rounding modes somehow.

> class (RealFrac a, Floating a) => RealFloat a where
>     atan2            :: a -> a -> a
> {- This needs lots of fromIntegral wrapping.
>     atan2 y x
>       | x>0           =  atan (y/x)
>       | x==0 && y>0   =  pi/2
>       | x<0  && y>0   =  pi + atan (y/x) 
>       |(x<=0 && y<0)  = -atan2 (-y) x
>       | y==0 && x<0   =  pi    -- must be after the previous test on zero y
>       | x==0 && y==0  =  y     -- must be after the other double zero tests
> -}

(Note that I removed the IEEEFloat-specific calls here, so probably
nobody will actually use this default definition.)

> class (Real a, Integral a) => RealIntegral a where
>     quot, rem        :: a -> a -> a   
>     quotRem          :: a -> a -> (a,a)
>
>       -- Minimal definition: nothing required
>     quot a b = let (q,_) = quotRem a b in q
>     rem a b  = let (_,r) = quotRem a b in r
>     quotRem a b = let (d,m) = divMod a b in
>                    if (signum d < (fromInteger 0)) then
>	                   (d+(fromInteger 1),m-b) else (d,m)

Remember that divMod does not specify exactly a `quot` b should be,
mainly because there is no sensible way to define it in general.  For
an instance of RealIntegral a, it is expected that a `quot` b will
round towards minus infinity and a `div` b will round towards 0.

> class (Real a) => SmallReal a where
>     toRational :: a -> Rational
> class (SmallReal a, RealIntegral a) => SmallIntegral a where
>     toInteger :: a -> Integer

These two classes exist to allow convenient conversions, primarily
between the built-in types.  These classes are "small" in the sense
that they can be converted to integers (resp. rationals) without loss
of information.  They should satisfy

    fromInteger . toInteger === id
  fromRational . toRational === id
     toRational . toInteger === toRational

> --- Numerical functions
> subtract         :: (Additive a) => a -> a -> a
> subtract         =  flip (-)
>
> even, odd        :: (Eq a, Integral a) => a -> Bool
> even n           =  n `mod` (one + one) == fromInteger zero
> odd              =  not . even

Additional standard libraries might include IEEEFloat (including the
bulk of the functions in Haskell 98's RealFloat class), VectorSpace,
Ratio, and Lattice.

> -- Support functions so that this whole thing can be tested on top
> -- of a standard prelude.
> -- Alternative: use "newtype".
> type Integer = P.Integer
> type Int = P.Int
> type Float = P.Float
> type Double = P.Double
> type Rational = P.Rational -- This one is lame.

> instance Additive P.Integer where
>     (+)    = (P.+)
>     zero   = 0
>     negate = P.negate
> instance Num P.Integer where
>     (*)    = (P.*)
>     one    = 1
> instance Integral P.Integer where
>     divMod = P.divMod
> instance Real P.Integer
> instance RealIntegral P.Integer
> instance SmallReal P.Integer where
>     toRational = P.toRational
> instance SmallIntegral P.Integer where
>     toInteger = id

> data T a = T a
--dDRMvlgZJXvWKvBx
Content-Type: text/plain; charset=us-ascii
Content-Disposition: attachment; filename="VectorSpace.lhs"

> module VectorSpace where
> import NumPrelude
> import qualified Prelude
>
> -- Is this right?
> infixl 7 *>, <*
>
> class (Num a, Additive b) => Module a b where
>     (*>) :: a -> b -> b

A module over a ring satisfies:

   a *> (b + c) === a *> b + a *> c
   (a * b) *> c === a *> (b *> c)
   (a + b) *> c === a *> c + b *> c

For instance, the following function can be used to define any
Additive as a module over Integer:

> integerMultiply :: (SmallIntegral a, Additive b) => a -> b -> b
> integerMultiply a b = reduceRepeated (+) zero b (toInteger a)

There are no instance declarations by default, since they would
overlap with too many other instances and would be slower than
desired.

> class (Num a, Additive b) => RightModule a b where
>     (<*) :: b -> a -> b

> class (Fractional a, Additive b) => VectorSpace a b

> class (VectorSpace a b) => DivisibleSpace a b where
>     (</>) :: b -> b -> a

DivisibleSpace is used for free one-dimensional vector spaces.  It
satisfies

  (a </> b) *> b = a

Examples include dollars and kilometers.
--dDRMvlgZJXvWKvBx
Content-Type: text/plain; charset=us-ascii
Content-Disposition: attachment; filename="PowerSeries.lhs"

> module PowerSeries where
> import NumPrelude
> import qualified Prelude as P
> import VectorSpace
> import Prelude hiding (
>        Int, Integer, Float, Double, Rational, Num(..), Real(..),
>        Integral(..), Fractional(..), Floating(..), RealFrac(..),
>        RealFloat(..), subtract, even, odd,
>        gcd, lcm, (^), (^^))

Power series, either finite or unbounded.  (zipWith does exactly the
right thing to make it work almost transparently.)

> newtype PowerSeries a = PS [a] deriving (Eq, Ord, Show)
> stripPS (PS l) = l
> truncatePS :: Int -> PowerSeries a -> PowerSeries a
> truncatePS n (PS a) = PS (take n a)

Note that the derived instances only make sense for finite series.

> instance (Additive a) => Additive (PowerSeries a) where
>     negate (PS l) = PS (map negate l)
>     (PS a) + (PS b) = PS (zipWith (+) a b)
>     zero = PS (repeat zero)
>
> instance (Num a) => Num (PowerSeries a) where
>     one = PS (one:repeat zero)
>     fromInteger n = PS (fromInteger n : repeat zero)
>     PS (a:as) * PS (b:bs) = PS ((a*b):stripPS (a *> PS bs + PS as*PS (b:bs)))
>     PS _ * PS _ = PS []
>
> instance (Num a) => Module a (PowerSeries a) where
>     a *> (PS bs) = PS (map (a *) bs)

It would be nice to also provide:

  instance (Module a b) => Module a (PowerSeries b) where
      a *> (PS bs) = PS (map (a *>) bs)

maybe with

  instance (Num a) => Module a a where
      (*>) = (*)

> instance (Integral a) => Integral (PowerSeries a) where
>     divMod a b = (\(x,y)-> (PS x, PS y)) (unzip (aux a b))
>        where aux (PS (a:as)) (PS (b:bs)) =
>                 let (d,m) = divMod a b in
>                 (d,m):aux (PS as - d *> (PS bs)) (PS (b:bs))
>              aux _ _ = []


--dDRMvlgZJXvWKvBx--