--
--
-- Arithmetic Operation
--
--
declare symbol i, w, e, π: MathExpr

def toMathExpr {a} (arg: a) : MathExpr := mathNormalize (toMathExpr' arg)

def (+') : MathExpr -> MathExpr -> MathExpr := i.+
def (-') : MathExpr -> MathExpr -> MathExpr := i.-
def (*') : MathExpr -> MathExpr -> MathExpr := i.*
def (/') : MathExpr -> MathExpr -> MathExpr := i./

def plusForMathExpr (x: MathExpr) (y: MathExpr) : MathExpr :=
  mathNormalize (x +' y)

def minusForMathExpr (x: MathExpr) (y: MathExpr) : MathExpr :=
  mathNormalize (x -' y)

def multForMathExpr (x: MathExpr) (y: MathExpr) : MathExpr :=
  mathNormalize (x *' y)

def divForMathExpr (x: MathExpr) (y: MathExpr) : MathExpr :=
  x /' y

def sum {Num a} (xs: [a]) : a := foldl (+) 0 xs
def sum' (xs: [MathExpr]) : MathExpr := foldl (+') 0 xs

def product {Num a} (xs: [a]) : a := foldl (*) 1 xs
def product' (xs: [MathExpr]) : MathExpr := foldl (*') 1 xs

def power (x: MathExpr) (n: MathExpr) : MathExpr := mathNormalize (power' x n)
def power' (x: MathExpr) (n: MathExpr) : MathExpr := foldl (*') 1 (take n (repeat1 x))

def exp (x: MathExpr) : MathExpr := 'exp x

def (^) (x: MathExpr) (n: MathExpr) : MathExpr :=
  if x = e
    then exp n
    else if isRational n
      then if n >= 0
        then if isInteger n then power x n else '(^) x n
        else 1 / x ^ i.neg n
      else '(^) x n

def (^') (x: MathExpr) (n: MathExpr) : MathExpr :=
  if x = e
    then exp n
    else if isRational n
      then if n >= 0
        then if isInteger n then power' x n else '(^) x n
        else 1 /' x ^' i.neg n
      else '(^) x n

def gcd (x: MathExpr) (y: MathExpr) : MathExpr :=
  match (x, y) as (termExpr, termExpr) with
    | (_, #0) -> x
    | (#0, _) -> y
    | (term $a $xs, term $b $ys) ->
      gcd' (i.abs a) (i.abs b) *' foldl (*') 1 (map (\(s, n) -> s ^' n) (AC.intersect xs ys))

def gcd' (x: Integer) (y: Integer) : Integer :=
  match (x, y) as (integer, integer) with
    | (_, #0) -> x
    | (#0, _) -> y
    | (_, ?(>= x)) -> gcd' (i.modulo y x) x
    | (_, _) -> gcd' y x

def P./ fx gx x :=
  let xs := reverse (coefficients fx x)
      ys := reverse (coefficients gx x)
      (zs, rs) := L./ xs ys
   in ( sum' (map2 (\c n -> c *' x ^' n) (reverse zs) nats0)
      , sum' (map2 (\c n -> c *' x ^' n) (reverse rs) nats0) )