packages feed

egison-5.1.0: lib/math/normalize.egi

--
--
-- Term Rewriting
--
--

-- Built-in normalization. Now reduced to the residual Haskell-side rewriter
-- (`symbolNormalize` -> `casRewriteSymbol` -> `casRewriteDd`); all other
-- rewrite rules live as `declare rule auto` declarations below.
--
-- The user-facing `mathNormalize` is defined as an alias that desugar
-- overrides per `declare rule auto`. The override pattern is:
--
--   def mathNormalize := \v -> iterateRulesCAS [autoRule.0, ...] (mathNormalizeBuiltin v)
--
-- so user-declared auto rules apply on each `mathNormalize` call (which is
-- triggered by every `+`/`*`/`/` operation on MathValue, see arithmetic.egi).
def mathNormalizeBuiltin (x: MathValue) : MathValue := symbolNormalize x

def mathNormalize (x: MathValue) : MathValue := mathNormalizeBuiltin x

--
-- Built-in auto rules (migrated from Math/Rewrite.hs).
-- Each `declare rule auto` extends mathNormalize so the rule fires on every
-- arithmetic operation (+, -, *, /, ^).
--

-- Imaginary unit: i^2 = -1.
-- Replaces casRewriteI in Math/Rewrite.hs. The auto-rule engine handles
-- arbitrary powers (i^3, i^4, ...) by iterating the rule on each term.
declare rule auto term i^2 = -1

-- Cube root of unity: the primitive cube root w satisfies w^2 + w + 1 = 0.
-- Declared as an ideal (G3 of design/cas-simplification.md): the Groebner
-- basis of the generator is computed once and registered as term-level
-- rewrite rules.  The single generated rule w^2 -> -1 - w subsumes the
-- previously hand-written pair (w^3 = w * w^2 reduces in two steps), so
-- the completeness of the rule set is a theorem instead of a convention.
declare ideal [w^2 + w + 1]

-- Logarithm identities (replaces casRewriteLog in Math/Rewrite.hs).
-- log 1 = 0 and log e = 1 are also handled at the lib `log` function entry,
-- but are kept here so the symbolic 'log form normalises consistently.
-- log (e^n) reaches the rule as log (exp n) because (^) rewrites e^n to exp n
-- in lib/math/common/arithmetic.egi.
declare rule auto term log 1 = 0
declare rule auto term log e = 1
declare rule auto term log (exp $n) = n

-- Exponential identities (replaces casRewriteExp in Math/Rewrite.hs).
-- The multi-factor cases ((exp x)^n, exp x * exp y) are now expressible
-- because lib/math/expression.egi was extended with apply1-4 patterns on
-- the multExpr matcher (so the "rest of factors" slot accepts apply1).
declare rule auto term exp 0 = 1
declare rule auto term exp 1 = e
declare rule auto term exp ($n * i * π) = (-1)^n
-- The structural exp rules ((exp x)^n -> exp (n x) and the multi-factor
-- product merge) are implemented in Haskell (Math/Rewrite.hs,
-- casRewriteExp), for the same per-term match-cost reason as the sqrt
-- rules; the value rules above stay here.

-- Power identities (replaces casRewritePower in Math/Rewrite.hs).
-- (x^y)^n = x^(n*y) for n >= 2; n=1 returns the term unchanged.
-- The pattern matches a single-factor term containing one apply2 #(^) factor.
declare rule auto term term $c ((apply2 #(^) $x $y, $n) :: []) =
  if n >= 2
    then c *' '(^) x (n *' y)
    else c *' '(^) x y

-- x^y * x^z = x^(y+z) when both factors are apply2 #(^) with the same base.
-- Exact 2-factor variant; multi-factor cases stay in casRewritePower for now.
declare rule auto term term $c ((apply2 #(^) $x $y, #1) :: (apply2 #(^) #x $z, #1) :: []) =
  c *' '(^) x (y +' z)

-- nth-root power reduction (replaces casRewriteRt in Math/Rewrite.hs).
-- (rt n x)^k for k >= n -> (rt n x)^(k mod n) * x^(k div n).
-- The isInteger guard skips the rule for symbolic n.
declare rule auto term term $c ((apply2 #rt $n $x, $k) :: []) =
  if isInteger n
    then if k >= n
      then c *' x^'(i.quotient k n) *' ('rt n x)^'(i.modulo k n)
      else c *' ('rt n x)^'k
    else c *' ('rt n x)^'k

-- nth-root-of-unity power reduction (replaces casRewriteRtu in
-- Math/Rewrite.hs).
-- - k >= n: reduce exponent via mod n.
-- - k = n-1: apply minimal-polynomial reduction
--   (rtu n)^(n-1) = -1 - rtu n - (rtu n)^2 - ... - (rtu n)^(n-2).
-- - otherwise: keep as-is.
-- Note: the original casRewriteRtu g stage was buggy (foldr casMinus gave
-- alternating signs); the foldl form below is the correct minimal-polynomial.
declare rule auto term term $c ((apply1 #rtu $n, $k) :: []) =
  if isInteger n
    then if k >= n
      then c *' ('rtu n)^'(i.modulo k n)
      else if k = n - 1
        then c *' (foldl (+') (-1) (map (\j -> (-1) *' ('rtu n) ^' j) (between 1 (n - 2))))
        else c *' ('rtu n)^'k
    else c *' ('rtu n)^'k

-- Square root power reduction and pair merging are implemented in
-- Haskell (Math/Rewrite.hs, casRewriteSqrt): the declare-rule versions
-- paid a pattern-match attempt on every term of every sqrt-carrying
-- value per normalization, making arithmetic on such values ~20x
-- slower (thurston.egi's bottleneck; design/cas-simplification.md G6).

-- Note: FunctionData same-shape term merging stays in casRewriteDd
-- (Math/Rewrite.hs). The poly-level declare rule version
--   $a * ($f & func $g $args) * $mr + $b * (func #g #args) * #mr + $rest
--     = (a + b) * f * mr + rest
-- works correctly but the multi-term + same-binding constraint is too
-- expensive for complex differential-form computations like
-- riemann-curvature-tensor-of-S2xS3 (>120s vs few-seconds with Haskell).

-- abs of a manifestly non-negative monomial reduces to the monomial itself.
-- A monomial is non-negative if all its symbol exponents are even AND the
-- coefficient is a non-negative rational. This handles common cases like
-- `abs(sin²θ * r⁴)` (showing up in spherical Laplacians, hodge-spherical
-- etc.) by stripping the abs wrapper.
--
-- The LHS pattern triggers on any abs factor; the body re-matches the inner
-- argument as a single Term and checks the conditions before stripping.
declare rule auto term term $c ((apply1 #abs $a, $n) :: $rr) =
  match a as termExpr with
    | term $aCoeff $aMs ->
        if isRational aCoeff && aCoeff >= 0
            && all (\(_, k) -> i.modulo k 2 = 0) aMs
          then c *' (foldl (*') aCoeff (map (\(p, k) -> p ^' k) aMs)) ^' n
                 *' foldl (*') 1 (map (\(p, k) -> p ^' k) rr)
          else c *' ('abs a) ^' n
                 *' foldl (*') 1 (map (\(p, k) -> p ^' k) rr)
    | _ -> c *' ('abs a) ^' n
             *' foldl (*') 1 (map (\(p, k) -> p ^' k) rr)

--
-- sin/cos pythagorean identity (replaces lib's rewriteRuleForSinAndCos)
--

-- a*mr + (-a)*cos(x)^2*mr -> a*sin(x)^2*mr (since 1 - cos²x = sin²x)
declare rule auto poly $a * $mr + #(- a) * (apply1 #cos $x) ^ #2 * #mr + $pr =
  a *' (sin x)^2 *' mr +' pr

-- a*cos(x)^2*mr + b*sin(x)^2*mr -> a*mr + (b-a)*sin(x)^2*mr
-- (sin²+cos²=1 form: a*cos² + a*sin² + (b-a)*sin² = a + (b-a)*sin²)
declare rule auto poly $a * (apply1 #cos $x) ^ #2 * $mr + $b * (apply1 #sin #x) ^ #2 * #mr + $pr =
  a *' mr +' (b -' a) *' (sin x)^2 *' mr +' pr