root.egi

--
--
-- Algebra
--
--

--
-- Root
--
def rt n x :=
  if isInteger n
    then 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 x :=
  rtTerm n (numerator x * denominator x ^ (n - 1)) / denominator x

def rtTerm n x :=
  match x as termExpr with
    | term $a _ ->
      let rtm1 n := 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 x :=
  match (n, x) as (mathExpr, mathExpr) with
    | (#3, $a * #i * $r) -> (- i) * rt 3 (a *' r)
    | (_, $a * #sqrt $b * $r) -> rt (n * 2) (a ^' 2 *' b) *' rt n r
    | (_, $a * #rt $n' $b * $r) -> rt (n * n') (a ^' n' *' b) *' rt n r
    | (_, _) -> rtPositiveTerm1 n x
  where
    rtPositiveTerm1 n x :=
      let f xs :=
            match xs as assocMultiset mathExpr with
              | [] -> (1, 1)
              | ncons $p $k $rs ->
                let (a, b) := f rs
                 in (p ^' quotient k n *' a, p ^' (k % n) *' b)
          g n x :=
            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 x :=
  match (n, x) as (integer, integer) with
    | (#2, _) -> `sqrt x
    | (_, _) -> `rt n x

def sqrt x :=
  if isScalar x
    then let m := numerator x
             n := denominator x
          in rt 2 (m * n) / n
    else b.sqrt x

def rtOfUnity := rtu

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