root.egi

--
--
-- Algebra
--
--

--
-- Root
--
declare mathfunc rt
declare apply rt n x :=
  if isInteger n
    then
      if n = 1
        then x
        else
          match x as mathValue with
            | #0 -> 0
            | ?isMonomial -> rtMonomial n x
            | poly $xs / poly $ys ->
                let xd := reduce gcdForMathValue xs
                    yd := reduce gcdForMathValue 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: MathValue) (x: MathValue) : MathValue :=
  rtTerm n (numerator x * denominator x ^ (n - 1)) / denominator x

def rtTerm (n: MathValue) (x: MathValue) : MathValue :=
  match x as termExpr with
    | term $a _ ->
      let rtm1 (n: MathValue) : MathValue := 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: MathValue) (x: MathValue) : MathValue :=
  match (n, x) as (mathValue, mathValue) 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: MathValue) (x: MathValue) : MathValue :=
      let f (xs: [(MathValue, MathValue)]) : (MathValue, MathValue) :=
            match xs as assocMultiset mathValue with
              | [] -> (1, 1)
              | ($p, $k) :: $rs ->
                  let (a, b) := f rs
                   in (p ^' i.quotient k n *' a, p ^' (k % n) *' b)
          g (n: MathValue) (x: MathValue) : MathValue :=
            let d := match x as termExpr with
                        | term $m $xs ->
                            gcdForMathValue n (reduce gcdForMathValue (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: MathValue) (x: MathValue) : MathValue :=
  match (n, x) as (integer, integer) with
    | (#2, _) ->
        -- Principal-branch normalization for constant radicands
        -- (design/cas-simplification.md 3.7, option A): a symbol-free
        -- radicand certified negative by interval arithmetic becomes
        -- i * sqrt(-x).  Every surviving sqrt atom then has a positive
        -- radicand, which makes the pair merge sqrt a * sqrt b =
        -- sqrt (a b) sound on principal branches (the formal merge is
        -- off by a sign when both radicands are negative).
        if signOfConst x = "neg"
          then i *' sqrt (- x)
          else sqrtDenest x
    | (_, _) -> 'rt n x

-- Depth-2 denesting (design/cas-simplification.md 3.6; the classic
-- Borodin-Fagin-Hopcroft-Tompa condition):
--   sqrt (a + b sqrt c), with integers a > 0, b /= 0, c > 0 and
--   a^2 - b^2 c = d^2 a perfect square, denests to
--   sqrt ((a+d)/2) + sign(b) * sqrt ((a-d)/2).
-- Example: sqrt (9 - 4 sqrt 5) = sqrt 5 - 2.  Anything outside this
-- shape (or with a non-square a^2 - b^2 c, e.g. sqrt (-10 - 2 sqrt 5))
-- stays symbolic.  The perfect-square test factorizes a^2 - b^2 c, so
-- it is guarded by a size cap.
def sqrtDenest (x: MathValue) : MathValue :=
  match x as mathValue with
    | poly [term $a [], term $b [(apply1 #sqrt $c, #1)]] ->
        if isInteger a && isInteger b && isInteger c && a > 0 && c > 0
          then
            let r := a^2 - b^2 * c
             in if r > 0 && r <= 1000000000000
                  then sqrtDenest' x a b r
                  else 'sqrt x
          else 'sqrt x
    | _ -> 'sqrt x

def sqrtDenest' (x: MathValue) (a: MathValue) (b: MathValue) (r: MathValue) : MathValue :=
  let facs := toAssoc (pF r)
   in if all (\(_, k) -> i.modulo k 2 = 0) facs
        then
          let d := foldl (\acc (p, k) -> acc * p ^ (i.quotient k 2)) 1 facs
              s := if b > 0 then 1 else -1
           in sqrt ((a + d) / 2) + s * sqrt ((a - d) / 2)
        else 'sqrt x

-- sqrt is split into `declare mathfunc` (declares the function name and
-- registers default symbolic behaviour) plus `declare apply` (the
-- algorithmic simplification at application time). Pattern rewrites such
-- as (sqrt $x)^2 = x live separately as `declare rule auto` declarations
-- in lib/math/normalize.egi.
declare mathfunc sqrt
declare apply sqrt x :=
  let m := numerator x
      n := denominator x
   in rt 2 (m *' n) /' n

declare mathfunc rtu
declare apply rtu n :=
  if isInteger n
    then match n as integer with
      | #1 -> 1
      | #2 -> -1
      | #3 -> w
      | #4 -> i
      | _ -> 'rtu n
    else 'rtu n

def rtOfUnity : MathValue -> MathValue := rtu