packages feed

egison-5.0.0: lib/math/algebra/root.egi

--
--
-- Algebra
--
--

--
-- Root
--
def rt (n: MathExpr) (x: MathExpr) : MathExpr :=
  if isInteger n
    then
      if n = 1
        then x
        else
          match x as mathExpr with
            | #0 -> 0
            | ?isMonomial -> rtMonomial n x
            | poly $xs / poly $ys ->
                let xd := reduce gcd xs
                    yd := reduce gcd ys
                    d := rtMonomial n (xd / yd)
                 in d *' rt'' n (sum' (map (/' xd) xs) /' sum' (map (/' yd) ys))
            | _ -> rt'' n x
    else rt'' n x

def rtMonomial (n: MathExpr) (x: MathExpr) : MathExpr :=
  rtTerm n (numerator x * denominator x ^ (n - 1)) / denominator x

def rtTerm (n: MathExpr) (x: MathExpr) : MathExpr :=
  match x as termExpr with
    | term $a _ ->
      let rtm1 (n: MathExpr) : MathExpr := match n as integer with
                    | #1 -> -1
                    | #2 -> i
                    | ?isOdd -> -1
                    | _ -> undefined
       in if a < 0 then rtm1 n *' rtPositiveTerm n (- x) else rtPositiveTerm n x

def rtPositiveTerm (n: MathExpr) (x: MathExpr) : MathExpr :=
  match (n, x) as (mathExpr, mathExpr) with
    | (#3, $a * #i * $r) -> (- i) * rt 3 (a *' r)
    | (_, $a * (apply1 #sqrt $b) * $r) -> rt (n * 2) (a ^' 2 *' b) *' rt n r
    | (_, $a * (apply2 #rt $n' $b) * $r) -> rt (n * n') (a ^' n' *' b) *' rt n r
    | (_, _) -> rtPositiveTerm1 n x
  where
    rtPositiveTerm1 (n: MathExpr) (x: MathExpr) : MathExpr :=
      let f (xs: [(MathExpr, MathExpr)]) : (MathExpr, MathExpr) :=
            match xs as assocMultiset mathExpr with
              | [] -> (1, 1)
              | ($p, $k) :: $rs ->
                  let (a, b) := f rs
                   in (p ^' i.quotient k n *' a, p ^' (k % n) *' b)
          g (n: MathExpr) (x: MathExpr) : MathExpr :=
            let d := match x as termExpr with
                        | term $m $xs ->
                            gcd n (reduce gcd (map snd (toAssoc (pF m) ++ xs)))
             in rt'' (n / d) (rt d x)
          in match x as termExpr with
            | term $m $xs ->
                match f (toAssoc (pF (abs m)) ++ xs) as (integer, integer) with
                  | ($a, #1) -> a
                  | ($a, $b) -> a *' g n b

def rt'' (n: MathExpr) (x: MathExpr) : MathExpr :=
  match (n, x) as (integer, integer) with
    | (#2, _) -> 'sqrt x
    | (_, _) -> 'rt n x

def sqrt (x: MathExpr) : MathExpr :=
  let m := numerator x
      n := denominator x
   in rt 2 (m *' n) /' n

def rtOfUnity : MathExpr -> MathExpr := rtu

def rtu (n: MathExpr) : MathExpr :=
  if isInteger n
    then match n as integer with
      | #1 -> 1
      | #2 -> -1
      | #3 -> w
      | #4 -> i
      | _ -> 'rtu n
    else 'rtu n