> {-# OPTIONS_GHC -fglasgow-exts #-}

Standard Forms

> module Simplify       (  standard1, standard2, standard3,
>                          summ, prod, (^^^),
>                          reduce, compact, distribute, lower
>                       )
> where
> import Ratio
> import List           (  groupBy  )
> import Array
>
> import Sort           (  mergeSortBy, mergeLists  )
>
> import Function       (  Function (..), Prim (..), constant,
>                          summands, factors, esplit,
>                          nat, zero, one, mone  )

> infixr 1 ...
> infixr 5 ^^^
> infixr 5 #
> infix  5 -->

 >, `pow`
 > infix  4 -=>

> f # g  = g . f
> sortBy = mergeSortBy

----------------------------------------------------------------------
Simplification
----------------------------------------------------------------------

> fold (...) summ prod (^^^)    =  fold
>     where
>     fold (f :.: g)    =  fold f ... fold g
>     fold (Summ fs)    =  summ [ fold f | f<-fs ]
>     fold (Prod fs)    =  prod [ fold f | f<-fs ]
>     fold (f :^: g)    =  fold f ^^^ fold g
>     fold f            =  f

----------------------------------------------------------------------
Note: maximum of 10 iterations
----------------------------------------------------------------------

> simplify red  =  iterate red # take 10 # limit
>
> reduce        =  fold (...) summ prod (^^^)

----------------------------------------------------------------------
Reduction rules

Composition
----------------------------------------------------------------------

> f         ... Id              =  f
> Ratio n   ... f               =  Ratio n
> Const c   ... f               =  Const c
> Id        ... f               =  f
> Prim f    ... (Prim g :.: h)
>    | f  `isInverse` g         =  h
> (f :.: g) ... h               =  f ... (g ... h)
> (Summ fs) ... h               =  Summ [ f ... h | f<-fs ]
> (Prod fs) ... h               =  Prod [ f ... h | f<-fs ]
> (f :^: g) ... h               =  (f ... h) ^^^ (g ... h)
> f         ... g               =  f :.: g              -- catch all

> f `isInverse` g     = (f,g) `elem` inverses || (g,f) `elem` inverses
> inverses = [ ( Sin,     Arcsin ),
>              ( Cos,     Arccos ),
>              ( Tan,     Arctan ),
>              ( Cot,     Arccot ),
>              ( Exp,     Ln     ) ]

----------------------------------------------------------------------
Addition.
----------------------------------------------------------------------

> summ []               =  zero
> summ [f]              =  f
> summ fs               =  Summ (peval (mergeLists [ summands f | f<-fs ]))
>     where
>     peval (Ratio m:Ratio n:fs)        =  peval (Ratio (m + n):fs)
>     peval fs                          =  fs

----------------------------------------------------------------------
Multiplication.
----------------------------------------------------------------------

> prod []               =  one
> prod [f]              =  f
> prod fs
>     | zero `elem` gs  =  zero
>     | otherwise       =  Prod (peval gs)
>     where
>     gs                                =  mergeLists [ factors f | f<-fs ]
>     peval (Ratio m:Ratio n:fs)        =  peval (Ratio (m * n):fs)
>     peval fs                          =  fs

----------------------------------------------------------------------
Exponentiation.
----------------------------------------------------------------------

> Ratio m   ^^^ Ratio r                      -- Ratio (n :% 1)
>     | denominator r == 1
>     , n <- numerator r
>                     = case n >= 0 of
>                          True        ->  Ratio (m ^ n)
>                          False       ->  Ratio (recip m ^ negate n)
> --Ratio n   ^^^ f     | n == 0        =  zero         -- f > 0 !!
> f         ^^^ Ratio n | n == 0        =  one          -- f /= 0 ??
> Ratio n   ^^^ f       | n == 1        =  one
> f         ^^^ Ratio n | n == 1        =  f
> (Prod fs) ^^^ h                       =  prod [ f ^^^ h | f<-fs ]
> (f :^: g) ^^^ h                       =  f ^^^ prod [g, h]
> f         ^^^ g                       =  f :^: g              -- catch all

----------------------------------------------------------------------
Collection of Partial Terms
----------------------------------------------------------------------

> compact       :: Function -> Function
> compact       =  fold (...) scompact pcompact (^^^)
>
> scompact      =  map fsplit # sortBy le2 # groupBy eq2 # map melt # summ
>     where
>     melt []           =  zero
>     melt ((c, f):cfs) =  prod [summ (c : [ c | (c, _)<-cfs ]), f]
>
> fsplit (Prod (m@(Ratio _):fs))        =  (m, Prod fs)
> fsplit f                              =  (one, f)
>
> pcompact      =  map esplit # sortBy le1 # groupBy eq1 # map melt # prod
>     where
>     melt []           =  one
>     melt ((b, e):bes) =  b ^^^ summ (e : [ e | (_, e)<-bes ])

----------------------------------------------------------------------
Note: this is important, that the factors will simplify (`^^^` instead of ':^:'),
otherwise fsplit will not work (the Pattern 'Prod (m@(Ratio _):fs)' won't match
----------------------------------------------------------------------

> eq1 (a, _) (b, _)     =  a == b
> le1 (a, _) (b, _)     =  a <= b
>
> eq2 (_, a) (_, b)     =  a == b
> le2 (_, a) (_, b)     =  a <= b

> standard1             =  simplify (reduce # compact)

----------------------------------------------------------------------

> standard2             =  simplify (reduce # compact # distribute)

> distribute            =  fold (...) summ mult (^^^)

> mult fs               =  Summ [ Prod gs | gs<-distr [ summands f | f<-fs ] ]

> distr []              =  [[]]
> distr (x:xs)          =  [ a:as | a<-x, as<-distr xs ]

----------------------------------------------------------------------
Standardform 3
----------------------------------------------------------------------

> standard3 n           =  simplify (reduce # compact # lower n)

> isWhole  r            = denominator r == 1

> lower n               =  fold (...) summ prod pow
>     where
>     (Summ (f:fs)) `pow` (Ratio r)     -- (i :% 1)
>         | isWhole r , i <- denominator r
>         , i > 0 &&  i <= n    =  binom i f (Summ fs)
>         | isWhole r , i <- denominator r
>         , i < 0 && -i <= n    =  binom (-i) f (Summ fs) ^^^ mone
>     f `pow` g                 =  f ^^^ g              -- catch all

> a --> b = (a,b)

> binom n f g   =  summ [
>                      prod [nat (c!(n,k)), f^^^nat k, g^^^nat (n-k)]
>                  | k<-[0..n] ]
>     where c   =  array ((0,0),(n,n)) (
>                  [ (i,0) --> 1 | i<-[0..n] ]
>               ++ [ (i,i) --> 1 | i<-[1..n] ]
>               ++ [ (i,j) --> c!(i-1,j) + c!(i-1,j-1) | i<-[2..n],j<-[1..i-1] ])

----------------------------------------------------------------------
Helper functions
----------------------------------------------------------------------

Walk down the list until too consecutive elements are the same

> limit []              =  error "Simplify.limit: empty list"
> limit [a]             =  a
> limit (a:x@(b:_))
>     | a == b          =  a
>     | otherwise       =  limit x