egison-5.1.0: hs-src/Language/Egison/Math/CAS.hs
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE PatternSynonyms #-}
{-# LANGUAGE QuasiQuotes #-}
{- |
Module : Language.Egison.Math.CAS
Licence : MIT
This module defines the new CASValue data type for the computer algebra system.
The type structure directly determines the runtime internal representation.
Key design principles:
- Type composition determines normal form (e.g., Poly (Frac Integer) vs Frac (Poly Integer))
- Supports Laurent polynomials (negative exponents allowed in monomials)
- Constructive: coefficients can be recursively nested CASValues
-}
module Language.Egison.Math.CAS
( -- * Core data types
CASValue (..)
, CASTerm (..)
, Monomial
, SymbolExpr (..)
, Id
, makeApplyExpr
-- * Smart constructors
, casInteger
, casFactor
, casPoly
, casFrac
, casTerm
-- * Arithmetic operations
, casPlus
, casMult
, casNegate
, casMinus
, casDivide
, casPower
, casNumerator
, casDenominator
-- * Observed type (Phase 8)
, prettyTypeOf
, casAtomSet
, casDifferentialClosed
-- * Normalization
, casNormalize
, casNormalizePoly
, casReshapeAs
-- * Re-exports for Rewrite.hs (avoids a direct Data import cycle)
, prettyFunctionName
-- * Predicates
, casIsZero
, casIsAtom
-- * Pretty printing
, prettyCAS
-- * GCD operations
, casGcd
, casTermsGcd
-- * Pattern synonyms for CASValue
, pattern CASZero
, pattern CASSingleSymbol
, pattern CASSingleTerm
-- * Pattern matching (control-egison)
, CASM (..)
, CASTermM (..)
, CASSymbolM (..)
, casTerm'
, casTerm'M
, casTermM
, casSymbol
, casSymbolM
, casFunc
, casFuncM
, casApply1
, casApply1M
, casApply2
, casApply2M
, casApply3
, casApply3M
, casApply4
, casApply4M
, casQuote
, casNegQuote
, casNegQuoteM
, casQuoteFunction
, casQuoteFunctionM
, casEqualMonomial
, casEqualMonomialM
, casZero
, casZeroM
, casSingleTerm
, casSingleTermM
) where
import Data.List (sortBy, groupBy, intercalate, intersect, nub)
import Data.Ord (comparing)
import Data.Function (on)
import Data.Ratio ((%), numerator, denominator)
import Control.Egison
import Control.Monad (MonadPlus (..))
import Language.Egison.IExpr (Index (..))
import Language.Egison.Type.Types (Type(..), SymbolSet(..), TypeAtom(..))
import {-# SOURCE #-} Language.Egison.Data (WHNFData, prettyFunctionName)
-- | CASValue represents mathematical values in the CAS.
-- The structure is compositional: each constructor has a well-defined semantics.
data CASValue
= CASInteger Integer
-- ^ Base case: an integer value
| CASFactor SymbolExpr
-- ^ An atomic factor (generated by quote operator ')
-- Represents a symbol or function application that is not yet expanded
| CASPoly [CASTerm]
-- ^ A polynomial (sum of terms). Empty list represents zero.
-- Supports Laurent polynomials: monomial exponents can be negative.
| CASFrac CASValue CASValue
-- ^ A quotient: numerator / denominator
-- Only needed when denominator is non-monomial
deriving (Eq, Show)
-- | CASTerm represents a single term in a polynomial: coefficient × monomial
-- The coefficient is a CASValue, enabling nested polynomial structures.
data CASTerm = CASTerm CASValue Monomial
deriving (Eq, Show)
-- | We choose the definition 'monomials' without its coefficients.
-- ex. 2 x^2 y^3 is *not* a monomial. x^2 t^3 is a monomial.
type Monomial = [(SymbolExpr, Integer)]
-- | Identifier type for symbols
type Id = String
-- | SymbolExpr represents atomic symbolic expressions in the CAS.
-- NOTE: SymbolExpr uses CASValue for function arguments (Apply1-4, Quote, FunctionData).
data SymbolExpr
= Symbol Id String [Index CASValue]
| Apply1 CASValue CASValue
| Apply2 CASValue CASValue CASValue
| Apply3 CASValue CASValue CASValue CASValue
| Apply4 CASValue CASValue CASValue CASValue CASValue
| Quote CASValue -- For backtick quote: `expr
| QuoteFunction WHNFData -- For single quote on functions: 'func
| FunctionData CASValue [CASValue] -- fnname args
-- Manual Eq instance (QuoteFunction comparison uses function name)
instance Eq SymbolExpr where
Symbol id1 s1 js1 == Symbol id2 s2 js2 = id1 == id2 && s1 == s2 && js1 == js2
Apply1 f1 a1 == Apply1 f2 a2 = f1 == f2 && a1 == a2
Apply2 f1 a1 b1 == Apply2 f2 a2 b2 = f1 == f2 && a1 == a2 && b1 == b2
Apply3 f1 a1 b1 c1 == Apply3 f2 a2 b2 c2 = f1 == f2 && a1 == a2 && b1 == b2 && c1 == c2
Apply4 f1 a1 b1 c1 d1 == Apply4 f2 a2 b2 c2 d2 = f1 == f2 && a1 == a2 && b1 == b2 && c1 == c2 && d1 == d2
Quote m1 == Quote m2 = m1 == m2
QuoteFunction whnf1 == QuoteFunction whnf2 =
case (prettyFunctionName whnf1, prettyFunctionName whnf2) of
(Just n1, Just n2) -> n1 == n2
_ -> False -- Anonymous functions are never equal
FunctionData n1 k1 == FunctionData n2 k2 = n1 == n2 && k1 == k2
_ == _ = False
instance Show SymbolExpr where
show = prettySymbolExpr
-- | Pretty print a SymbolExpr
prettySymbolExpr :: SymbolExpr -> String
prettySymbolExpr (Symbol _ (':':':':':':_) []) = "#"
prettySymbolExpr (Symbol _ s []) = s
prettySymbolExpr (Symbol _ s js) = s ++ concatMap showIndex js
where
showIndex (Sup i) = "~" ++ prettyCAS' i
showIndex (Sub i) = "_" ++ prettyCAS' i
showIndex (SupSub i) = "~_" ++ prettyCAS' i
showIndex (DF _ _) = ""
showIndex (User i) = "|" ++ prettyCAS' i
prettySymbolExpr (Apply1 fn a1) = unwords [prettyApplyFn fn, prettyApplyArg a1]
prettySymbolExpr (Apply2 fn a1 a2) = unwords [prettyApplyFn fn, prettyApplyArg a1, prettyApplyArg a2]
prettySymbolExpr (Apply3 fn a1 a2 a3) = unwords [prettyApplyFn fn, prettyApplyArg a1, prettyApplyArg a2, prettyApplyArg a3]
prettySymbolExpr (Apply4 fn a1 a2 a3 a4) = unwords [prettyApplyFn fn, prettyApplyArg a1, prettyApplyArg a2, prettyApplyArg a3, prettyApplyArg a4]
prettySymbolExpr (Quote mExprs) = "`" ++ prettyCAS' mExprs
prettySymbolExpr (QuoteFunction whnf) = "'" ++ maybe "<function>" id (prettyFunctionName whnf)
prettySymbolExpr (FunctionData name args) = unwords (prettyApplyFn name : map prettyApplyArg args)
-- | Pretty print the function slot of an Apply1-4 / FunctionData. The function
-- reference is often stored as a CASPoly wrapper around a single CASFactor
-- (e.g. `'cos` is `CASPoly [CASTerm 1 [(QuoteFunction cos, 1)]]`); print such
-- wrappers transparently as the underlying SymbolExpr to avoid a spurious
-- `('cos) (θ)` rendering.
prettyApplyFn :: CASValue -> String
prettyApplyFn (CASFactor sym) = prettySymbolExpr sym
prettyApplyFn (CASPoly [CASTerm (CASInteger 1) [(sym, 1)]]) = prettySymbolExpr sym
prettyApplyFn v = prettyCAS' v
-- | Pretty print an argument to an Apply or FunctionData. Single-symbol
-- arguments (lifted to CASPoly with coefficient 1) print without parens
-- (e.g. `'cos θ` instead of `'cos (θ)`); composite arguments fall back to
-- the parenthesizing `prettyCAS'`.
prettyApplyArg :: CASValue -> String
prettyApplyArg (CASFactor sym) = prettySymbolExpr sym
prettyApplyArg (CASPoly [CASTerm (CASInteger 1) [(sym, 1)]]) = prettySymbolExpr sym
prettyApplyArg v = prettyCAS' v
-- | Pretty print a CASValue (basic version for SymbolExpr Show instance)
prettyCAS :: CASValue -> String
prettyCAS (CASInteger n) = show n
prettyCAS (CASFactor sym) = prettySymbolExpr sym
prettyCAS (CASPoly []) = "0"
prettyCAS (CASPoly terms) = prettyTerms terms
where
prettyTerms [] = "0"
prettyTerms (t:ts) = prettyTerm t ++ concatMap withSign ts
withSign term@(CASTerm coeff _)
| isNegative coeff = " - " ++ prettyTerm (negateCAST term)
| otherwise = " + " ++ prettyTerm term
prettyTerm (CASTerm coeff []) = prettyCAS coeff
prettyTerm (CASTerm (CASInteger 1) mono) = prettyMono mono
prettyTerm (CASTerm (CASInteger (-1)) mono) = "- " ++ prettyMono mono
-- Use the parenthesizing `prettyCAS'` for the coefficient so that nested
-- polynomial coefficients (e.g. `(2*z + 3) * x`) don't visually flatten
-- into the surrounding sum. Integer/Factor coeffs pass through unwrapped.
prettyTerm (CASTerm coeff mono) = prettyCAS' coeff ++ " * " ++ prettyMono mono
-- For multi-factor monomials, use ` * ` between factors when any is a
-- function-application form (Apply1-4 / FunctionData), since juxtaposition
-- would be ambiguous with the function-call syntax `f x`. Otherwise keep
-- the conventional `x y` juxtaposition for plain symbols. A single factor
-- needs no separator nor wrapping.
prettyMono [single] = prettyPow single
prettyMono mono
| any (isApplyFactor . fst) mono = intercalate " * " (map prettyPow' mono)
| otherwise = unwords (map prettyPow mono)
-- prettyPow' wraps Apply factors in parens for clarity in the explicit-`*`
-- form (e.g. `('cos θ) * r`).
prettyPow' (sym, 1) | isApplyFactor sym = "(" ++ prettySymbolExpr sym ++ ")"
prettyPow' p = prettyPow p
prettyPow (sym, 1) = prettySymbolExpr sym
-- An application form under an exponent needs parens: `g x y^2` would
-- read as g applied to x and y^2, not as (g x y) squared.
prettyPow (sym, n)
| isApplyFactor sym = "(" ++ prettySymbolExpr sym ++ ")^" ++ show n
| otherwise = prettyCAS' (CASFactor sym) ++ "^" ++ show n
isApplyFactor :: SymbolExpr -> Bool
isApplyFactor (Apply1 {}) = True
isApplyFactor (Apply2 {}) = True
isApplyFactor (Apply3 {}) = True
isApplyFactor (Apply4 {}) = True
isApplyFactor (FunctionData {}) = True
isApplyFactor _ = False
isNegative (CASInteger n) = n < 0
isNegative _ = False
negateCAST (CASTerm (CASInteger n) m) = CASTerm (CASInteger (-n)) m
negateCAST t = t
prettyCAS (CASFrac num denom) = prettyCAS' num ++ " / " ++ prettyCAS' denom
prettyCAS' :: CASValue -> String
prettyCAS' v@(CASInteger _) = prettyCAS v
prettyCAS' v@(CASFactor _) = prettyCAS v
prettyCAS' v = "(" ++ prettyCAS v ++ ")"
-- | Compute the observed type of a CASValue and pretty print it.
-- The observed type is the most specific static type that the value
-- inhabits, computed bottom-up from the runtime structure.
-- CASInteger _ → "Integer"
-- CASFactor (Symbol ...)
-- → "Symbol"
-- CASFactor _ → "Factor"
-- CASPoly [] → "Integer" (canonical zero)
-- CASPoly terms → "Poly C [atoms]" where C is the join of term coefficient
-- types and atoms is the sorted list of distinct flat atoms
-- CASFrac n d → "Frac (typeOf n)" if d is integer/poly with single term,
-- otherwise "Frac (Poly typeOf-num [..])"
prettyTypeOf :: CASValue -> String
prettyTypeOf (CASInteger _) = "Integer"
prettyTypeOf (CASFactor (Symbol _ _ _)) = "Symbol"
prettyTypeOf (CASFactor _) = "Factor"
prettyTypeOf (CASPoly []) = "Integer"
prettyTypeOf (CASPoly terms) =
let coeffTypes = map (\(CASTerm c _) -> prettyTypeOf c) terms
coeffType = joinObservedTypes coeffTypes
atoms = collectAtoms terms
atomStr = if null atoms
then "[]"
else "[" ++ commaSep atoms ++ "]"
in "Poly " ++ parenIfApp coeffType ++ " " ++ atomStr
prettyTypeOf (CASFrac n d) =
"Frac " ++ parenIfApp inner
where
nT = prettyTypeOf n
dT = prettyTypeOf d
-- If numerator and denominator share the observed type, that's the inner;
-- otherwise widen to MathValue.
inner = if nT == dT then nT else "MathValue"
-- | Pretty join of observed types of multiple coefficients.
--
-- "Integer" is the bottom of the CAS observed-type lattice (it embeds into
-- every other CAS type), so when one term is observed as "Integer" and the
-- rest as some richer type T, we report T (not the over-broad "MathValue").
-- If two truly distinct non-Integer types appear, we still widen to
-- "MathValue" — a full subtype-aware join over the observed-type strings is
-- left as future work.
joinObservedTypes :: [String] -> String
joinObservedTypes [] = "Integer"
joinObservedTypes ts = case filter (/= "Integer") ts of
[] -> "Integer"
t : ts' | all (== t) ts' -> t
| otherwise -> "MathValue"
-- | Phase 8 differential closure: collect the set of atoms (as their canonical
-- pretty form) appearing in a CASValue's monomials, recursing into nested
-- coefficients. The result is a sorted list of unique atom names.
casAtomSet :: CASValue -> [String]
casAtomSet (CASInteger _) = []
casAtomSet (CASFactor sym) = [prettySymbolExpr sym]
casAtomSet (CASPoly terms) =
let atomNames = concat
[ map (prettySymbolExpr . fst) mono ++ casAtomSet coeff
| CASTerm coeff mono <- terms
]
in unique (sortBy compare atomNames)
where
unique [] = []
unique (x:xs) = x : unique (dropWhile (== x) xs)
casAtomSet (CASFrac n d) =
let combined = casAtomSet n ++ casAtomSet d
in unique (sortBy compare combined)
where
unique [] = []
unique (x:xs) = x : unique (dropWhile (== x) xs)
-- | Check whether differentiation preserved the atom set of the input.
-- Used by the `differentialClosed` primitive: the result is true iff the
-- atom set of `output` is a subset of that of `input` (no new atoms
-- introduced by `∂/∂`).
casDifferentialClosed :: CASValue -> CASValue -> Bool
casDifferentialClosed input output =
let inA = casAtomSet input
outA = casAtomSet output
in all (`elem` inA) outA
-- | Collect distinct atom names from a list of CASTerm monomials.
-- Returns a sorted list of pretty atom forms (`x`, `sin x`, etc.).
collectAtoms :: [CASTerm] -> [String]
collectAtoms terms =
let atomNames = [prettySymbolExpr s | CASTerm _ mono <- terms, (s, _) <- mono]
in unique (sortBy compare atomNames)
where
unique [] = []
unique (x:xs) = x : unique (dropWhile (== x) xs)
-- | Comma-separate strings.
commaSep :: [String] -> String
commaSep [] = ""
commaSep [x] = x
commaSep (x:xs) = x ++ ", " ++ commaSep xs
-- | Wrap in parens if the type printed contains a space (i.e. an application).
parenIfApp :: String -> String
parenIfApp s
| ' ' `elem` s = "(" ++ s ++ ")"
| otherwise = s
-- | Helper function to create Apply constructors based on argument count
makeApplyExpr :: CASValue -> [CASValue] -> SymbolExpr
makeApplyExpr fn [a1] = Apply1 fn a1
makeApplyExpr fn [a1, a2] = Apply2 fn a1 a2
makeApplyExpr fn [a1, a2, a3] = Apply3 fn a1 a2 a3
makeApplyExpr fn [a1, a2, a3, a4] = Apply4 fn a1 a2 a3 a4
makeApplyExpr _ _ = error "makeApplyExpr: unsupported number of arguments (must be 1-4)"
--------------------------------------------------------------------------------
-- Smart Constructors
--------------------------------------------------------------------------------
-- | Create an integer CASValue
casInteger :: Integer -> CASValue
casInteger = CASInteger
-- | Create a factor CASValue from a SymbolExpr
casFactor :: SymbolExpr -> CASValue
casFactor = CASFactor
-- | Create a polynomial CASValue, normalizing the terms
casPoly :: [CASTerm] -> CASValue
casPoly terms = casNormalize (CASPoly terms)
-- | Create a division CASValue, simplifying if possible
casFrac :: CASValue -> CASValue -> CASValue
casFrac num denom = casNormalize (CASFrac num denom)
-- | Create a term
casTerm :: CASValue -> Monomial -> CASTerm
casTerm = CASTerm
--------------------------------------------------------------------------------
-- Predicates
--------------------------------------------------------------------------------
-- | Check if a CASValue is zero
casIsZero :: CASValue -> Bool
casIsZero (CASInteger 0) = True
casIsZero (CASPoly []) = True
casIsZero (CASFrac n _) = casIsZero n
casIsZero _ = False
-- | Check if a CASValue is one
casIsOne :: CASValue -> Bool
casIsOne (CASInteger 1) = True
casIsOne (CASPoly [CASTerm (CASInteger 1) []]) = True
casIsOne _ = False
-- | Check if a CASValue is atomic (no parentheses needed for display)
-- Returns True for atomic values that don't need parentheses for display
casIsAtom :: CASValue -> Bool
casIsAtom (CASInteger _) = True
casIsAtom (CASFactor _) = True
casIsAtom (CASPoly []) = True -- Zero
casIsAtom (CASPoly [CASTerm _ []]) = True -- Integer only
casIsAtom (CASPoly [CASTerm (CASInteger 1) [_]]) = True -- Single symbol with coeff 1
casIsAtom (CASFrac num (CASPoly [CASTerm (CASInteger 1) []])) = casIsAtom num -- n/1 = n
casIsAtom _ = False
--------------------------------------------------------------------------------
-- Arithmetic Operations
--------------------------------------------------------------------------------
-- | Add two CASValues
casPlus :: CASValue -> CASValue -> CASValue
casPlus a b = casNormalize (casPlus' a b)
casPlus' :: CASValue -> CASValue -> CASValue
-- Integer + Integer
casPlus' (CASInteger a) (CASInteger b) = CASInteger (a + b)
-- Poly + Poly
casPlus' (CASPoly ts1) (CASPoly ts2) = CASPoly (ts1 ++ ts2)
-- Integer + Poly: embed integer as polynomial term
casPlus' (CASInteger n) (CASPoly ts) = CASPoly (CASTerm (CASInteger n) [] : ts)
casPlus' (CASPoly ts) (CASInteger n) = CASPoly (CASTerm (CASInteger n) [] : ts)
-- Frac + Frac: cross-multiply and add numerators
casPlus' (CASFrac n1 d1) (CASFrac n2 d2) =
CASFrac (casPlus' (casMult' n1 d2) (casMult' n2 d1)) (casMult' d1 d2)
-- Frac + other: embed other as Frac
casPlus' (CASFrac n d) other = CASFrac (casPlus' n (casMult' other d)) d
casPlus' other (CASFrac n d) = CASFrac (casPlus' (casMult' other d) n) d
-- Factor handling: lift to polynomial before operation
casPlus' (CASFactor sym) other = casPlus' (liftFactorToPoly sym) other
casPlus' other (CASFactor sym) = casPlus' other (liftFactorToPoly sym)
-- | Negate a CASValue
casNegate :: CASValue -> CASValue
casNegate (CASInteger n) = CASInteger (-n)
casNegate (CASPoly terms) = CASPoly (map negateTerm terms)
where
negateTerm (CASTerm coeff mono) = CASTerm (casNegate coeff) mono
casNegate (CASFrac n d) = CASFrac (casNegate n) d
casNegate (CASFactor sym) = CASPoly [CASTerm (CASInteger (-1)) [(sym, 1)]]
-- | Subtract two CASValues
casMinus :: CASValue -> CASValue -> CASValue
casMinus a b = casPlus a (casNegate b)
-- | Multiply two CASValues
casMult :: CASValue -> CASValue -> CASValue
casMult a b = casNormalize (casMult' a b)
casMult' :: CASValue -> CASValue -> CASValue
-- Integer * Integer
casMult' (CASInteger a) (CASInteger b) = CASInteger (a * b)
-- Integer * Poly: scale all coefficients
casMult' (CASInteger n) (CASPoly ts) = CASPoly (map (scaleTerm n) ts)
where
scaleTerm k (CASTerm coeff mono) = CASTerm (casMult' (CASInteger k) coeff) mono
casMult' (CASPoly ts) (CASInteger n) = casMult' (CASInteger n) (CASPoly ts)
-- Poly * Poly: distribute
casMult' (CASPoly []) _ = CASPoly []
casMult' _ (CASPoly []) = CASPoly []
casMult' (CASPoly ts1) (CASPoly ts2) =
CASPoly [multTerms t1 t2 | t1 <- ts1, t2 <- ts2]
where
multTerms (CASTerm c1 m1) (CASTerm c2 m2) =
CASTerm (casMult' c1 c2) (combineMonomials m1 m2)
-- Frac * Frac: multiply numerators and denominators
casMult' (CASFrac n1 d1) (CASFrac n2 d2) =
CASFrac (casMult' n1 n2) (casMult' d1 d2)
-- Frac * other: multiply into numerator
casMult' (CASFrac n d) other = CASFrac (casMult' n other) d
casMult' other (CASFrac n d) = CASFrac (casMult' other n) d
-- Factor handling: lift to polynomial before operation
casMult' (CASFactor sym) other = casMult' (liftFactorToPoly sym) other
casMult' other (CASFactor sym) = casMult' other (liftFactorToPoly sym)
-- | Lift a Factor to a polynomial: sym → 1 * sym^1
liftFactorToPoly :: SymbolExpr -> CASValue
liftFactorToPoly sym = CASPoly [CASTerm (CASInteger 1) [(sym, 1)]]
-- | Combine two monomials by adding exponents of matching symbols
combineMonomials :: Monomial -> Monomial -> Monomial
combineMonomials m1 m2 = foldr insertSymbol m2 m1
where
insertSymbol (sym, expo) mono =
case lookup sym mono of
Just _ -> map (\(s, e) -> if s == sym then (s, e + expo) else (s, e)) mono
Nothing -> (sym, expo) : mono
-- | Divide two CASValues: a / b
casDivide :: CASValue -> CASValue -> CASValue
casDivide a b = casNormalize (CASFrac a b)
-- | Raise a CASValue to an integer power
casPower :: CASValue -> Integer -> CASValue
casPower _ 0 = CASInteger 1
casPower x 1 = x
casPower x n
| n > 0 = casMult x (casPower x (n - 1))
| otherwise = casDivide (CASInteger 1) (casPower x (-n)) -- Negative power
-- | Get the numerator of a CASValue.
-- For tower-fixed level 4 polynomials (Poly with Frac coefficients),
-- compute the LCM of coefficient denominators and return the value
-- multiplied by it (clearing the Fracs from coefficients).
casNumerator :: CASValue -> CASValue
casNumerator (CASFrac num _) = num
casNumerator x@(CASPoly ts) =
let lcmD = polyDenominatorLCM ts
in if lcmD == 1
then x
else casMult x (CASInteger lcmD)
casNumerator x = x
-- | Get the denominator of a CASValue.
-- For tower-fixed level 4 polynomials, the denominator is the LCM of
-- the Frac coefficients' denominators.
casDenominator :: CASValue -> CASValue
casDenominator (CASFrac _ denom) = denom
casDenominator (CASPoly ts) = CASInteger (polyDenominatorLCM ts)
casDenominator _ = CASInteger 1
-- | Compute the LCM of all denominators in a polynomial's Frac coefficients.
-- Returns 1 if there are no Frac coefficients.
polyDenominatorLCM :: [CASTerm] -> Integer
polyDenominatorLCM ts =
let denoms = concatMap termDenoms ts
in if null denoms then 1 else foldl1 lcm denoms
where
termDenoms (CASTerm (CASFrac _ (CASInteger d)) _) = [abs d]
termDenoms _ = []
--------------------------------------------------------------------------------
-- Normalization
--------------------------------------------------------------------------------
-- | Normalize a CASValue
casNormalize :: CASValue -> CASValue
casNormalize (CASInteger n) = CASInteger n
casNormalize (CASFactor sym) = CASFactor sym
casNormalize (CASPoly terms) = casNormalizePoly terms
casNormalize (CASFrac num denom) = casNormalizeFrac num denom
-- | Normalize a polynomial
-- Steps:
-- 1. Fold symbols within each term (x * x^2 → x^3)
-- 2. Remove zero-exponent symbols
-- 3. Fold terms with equal monomials
-- 4. Remove zero-coefficient terms
-- 5. Sort terms in descending order
casNormalizePoly :: [CASTerm] -> CASValue
casNormalizePoly = casNormalizePolyWith FlattenNested
-- | How polynomial normalization treats nested coefficients (Phase
-- gamma-prime of the extensible-tower plan,
-- design/type-cas-tower-implementation.md section 4).
data NestedCoeffPolicy
= FlattenNested
-- ^ The default on every arithmetic path: a coefficient that is itself
-- a CASPoly (the nested canonical form produced by reshape) is
-- distributed out into the outer monomial, so operations always exit
-- in the default flat canonical form and terms coming from nested and
-- flat representations merge (i + (-i) = 0 across representations).
| KeepNested
-- ^ Reshape's final grouping: keep the nested form just constructed.
-- Coefficients are neither re-normalized (they were just produced by a
-- recursive reshape; re-normalizing would flatten deeper nesting) nor
-- distributed. Invariant: nested forms exist only as the direct output
-- of reshape (annotation sites).
deriving (Eq)
-- | Core polynomial normalization, parametrized by the nested-coefficient
-- policy (see 'NestedCoeffPolicy').
casNormalizePolyWith :: NestedCoeffPolicy -> [CASTerm] -> CASValue
casNormalizePolyWith policy terms =
let flatten = policy == FlattenNested
-- Normalize each term's monomial
terms1 = map (normalizeTermMonomialWith flatten) terms
-- Distribute nested coefficients (flat exit form) when asked to
terms1' = if flatten && any hasNestedCoeff terms1
then map (normalizeTermMonomialWith flatten)
(concatMap flattenNestedTerm terms1)
else terms1
-- Fold terms with equal monomials
terms2 = foldTermsWith policy terms1'
-- Remove zero-coefficient terms
terms3 = filter (not . isZeroTerm) terms2
-- Sort in descending order (standard polynomial order)
terms4 = sortTermsDescending terms3
in case terms4 of
[] -> CASInteger 0 -- Empty polynomial is zero
[CASTerm coeff []] | isIntegerCoeff coeff -> extractInteger coeff
ts -> CASPoly ts
where
isZeroTerm (CASTerm coeff _) = casIsZero coeff
isIntegerCoeff (CASInteger _) = True
isIntegerCoeff _ = False
extractInteger (CASInteger n) = CASInteger n
extractInteger _ = error "extractInteger: not an integer"
-- | Does the term carry a nested coefficient that 'flattenNestedTerm'
-- would distribute? Cheap constructor check used as a fast-path guard.
hasNestedCoeff :: CASTerm -> Bool
hasNestedCoeff (CASTerm (CASPoly _) _) = True
hasNestedCoeff (CASTerm (CASFrac (CASPoly _) (CASInteger _)) _) = True
hasNestedCoeff _ = False
-- | Distribute a nested coefficient into flat terms: a CASPoly coefficient
-- is multiplied out into the outer monomial; a CASFrac coefficient with a
-- CASPoly numerator over an integer denominator is distributed likewise.
-- Recursion terminates because each step strictly reduces nesting.
flattenNestedTerm :: CASTerm -> [CASTerm]
flattenNestedTerm t@(CASTerm c mono) = case c of
CASPoly inner ->
concatMap (\(CASTerm ic im) -> flattenNestedTerm (CASTerm ic (im ++ mono))) inner
CASFrac (CASPoly inner) d@(CASInteger _) ->
concatMap (\(CASTerm ic im) ->
flattenNestedTerm (CASTerm (casNormalizeFrac ic d) (im ++ mono))) inner
_ -> [t]
-- | Normalize a term's monomial: combine duplicate symbols, remove zero
-- exponents. When `renormCoeff` is False the coefficient is left untouched
-- (reshape's keep-nested mode — the coefficient was just built by a
-- recursive reshape and re-normalizing it would flatten deeper nesting).
normalizeTermMonomialWith :: Bool -> CASTerm -> CASTerm
normalizeTermMonomialWith renormCoeff (CASTerm coeff mono) =
let -- Fold duplicate symbols
mono1 = foldMonomialSymbols mono
-- Remove zero-exponent symbols
mono2 = filter (\(_, exp) -> exp /= 0) mono1
-- Normalize the coefficient recursively (flatten mode only)
coeff' = if renormCoeff then casNormalize coeff else coeff
in CASTerm coeff' mono2
-- | Fold duplicate symbols in a monomial by adding their exponents
foldMonomialSymbols :: Monomial -> Monomial
foldMonomialSymbols mono =
let grouped = groupBy ((==) `on` fst) (sortBy (comparing (show . fst)) mono)
in concatMap combineGroup grouped
where
combineGroup :: [(SymbolExpr, Integer)] -> [(SymbolExpr, Integer)]
combineGroup [] = []
combineGroup grp@((sym, _):_) = [(sym, sum (map snd grp))]
-- | Fold terms with equal monomials by adding their coefficients; under
-- 'KeepNested' merged coefficients are normalized without flattening, so a
-- nested coefficient produced by reshape survives the grouping.
foldTermsWith :: NestedCoeffPolicy -> [CASTerm] -> [CASTerm]
foldTermsWith policy terms =
let grouped = groupBy equalMonos (sortBy (comparing termMonoKey) terms)
in concatMap combineTerms grouped
where
-- Use show-based key for consistent ordering
termMonoKey (CASTerm _ m) = map (\(s, e) -> (show s, e)) (sortBy (comparing (show . fst)) m)
termMono (CASTerm _ m) = sortBy (comparing (show . fst)) m
equalMonos t1 t2 = termMono t1 == termMono t2
normalizeCoeff c = case (policy, c) of
(KeepNested, CASPoly ts) -> casNormalizePolyWith KeepNested ts
_ -> casNormalize c
combineTerms [] = []
combineTerms grp@((CASTerm _ m):_) =
let mono = sortBy (comparing (show . fst)) m
coeffSum = foldr casPlus' (CASInteger 0) [c | CASTerm c _ <- grp]
in [CASTerm (normalizeCoeff coeffSum) mono]
-- | Sort terms in descending order (highest degree first)
-- Order: by total degree, then lexicographically by symbols
sortTermsDescending :: [CASTerm] -> [CASTerm]
sortTermsDescending = sortBy (flip (comparing termDegree) <> flip (comparing termSymbolsKey))
where
termDegree (CASTerm _ mono) = sum (map snd mono)
-- Use show-based key for consistent lexicographic ordering
termSymbolsKey (CASTerm _ mono) = map (\(s, e) -> (show s, e)) (sortBy (comparing (show . fst)) mono)
-- | Normalize a fraction
-- Steps:
-- 1. If denominator is 1, return numerator
-- 2. If numerator is 0, return 0
-- 3. Simplify using GCD
-- 4. Ensure positive denominator
casNormalizeFrac :: CASValue -> CASValue -> CASValue
casNormalizeFrac num denom =
let num' = casNormalize num
denom' = casNormalize denom
in case (num', denom') of
-- Zero numerator
(n, _) | casIsZero n -> CASInteger 0
-- Denominator is 1
(n, d) | casIsOne d -> n
-- Denominator is -1: negate numerator
(n, CASInteger (-1)) -> casNegate n
-- Integer / Integer: reduce by GCD and normalize sign
(CASInteger n, CASInteger d) ->
let g = gcd n d
-- Normalize sign: ensure positive denominator
sign = if d < 0 then -1 else 1
n' = sign * (n `div` g)
d' = abs (d `div` g)
in if d' == 1
then CASInteger n'
else CASFrac (CASInteger n') (CASInteger d')
-- Poly / Integer (constant denominator): per the type-promotion-tower
-- design (type-cas.md §実行時の型昇格タワー), constant denominators are
-- absorbed into each term's coefficient as a Frac. The result is
-- level 4 (Poly with Frac coefficients) instead of level 5 (Frac of
-- Poly), matching the canonical form for `Poly (Frac Integer) [..]`.
(CASPoly ts1, CASInteger d) | d /= 0 ->
let ts1' = map (\(CASTerm c m) ->
CASTerm (casNormalizeFrac c (CASInteger d)) m) ts1
in casNormalizePoly ts1'
-- Any / single-term-Poly with monomial denominator: per the
-- type-promotion-tower design, monomial denominators are absorbed
-- as negative exponents (Laurent polynomial form, level 3/4) rather
-- than left as level 5 Frac. The numerator is normalized to Poly
-- form (CASInteger n → [CASTerm n []], CASFactor sym → [CASTerm 1
-- [(sym, 1)]]) and each term's exponents are decremented by the
-- denominator's monomial.
(numV, CASPoly [CASTerm denomCoef denomMono]) | not (null denomMono) ->
let numTerms = case numV of
CASPoly ts -> ts
CASInteger n -> [CASTerm (CASInteger n) []]
CASFactor sym -> [CASTerm (CASInteger 1) [(sym, 1)]]
_ -> [CASTerm numV []]
ts' = map (divTermByMonomial denomCoef denomMono) numTerms
in casNormalizePoly ts'
-- Poly / Poly (non-monomial denominator): try to reduce by monomial GCD,
-- then by the univariate polynomial GCD; otherwise keep as level 5 Frac.
-- After a proper polynomial-GCD reduction we re-enter casNormalizeFrac:
-- the reduced pair is coprime (the second pass finds a constant GCD and
-- stops), and a denominator that collapsed to a constant or a monomial
-- is absorbed by the earlier branches.
(CASPoly ts1, CASPoly ts2) ->
let (ts1', ts2') = simplifyPolyDiv ts1 ts2
in case (ts1', ts2') of
(ts1'', [CASTerm (CASInteger 1) []]) -> casNormalizePoly ts1''
_ -> case univariateGcdReduce ts1' ts2' of
Just (ts1'', ts2'') ->
casNormalizeFrac (casNormalizePoly ts1'') (casNormalizePoly ts2'')
Nothing -> case multivariateGcdReduce ts1' ts2' of
Just (ts1'', ts2'') ->
casNormalizeFrac (casNormalizePoly ts1'') (casNormalizePoly ts2'')
Nothing -> CASFrac (casNormalizePoly ts1') (casNormalizePoly ts2')
-- a / (b / c) = (a * c) / b
(n, CASFrac b c) -> casNormalizeFrac (casMult n c) b
-- (a / b) / c = a / (b * c)
(CASFrac a b, c) -> casNormalizeFrac a (casMult b c)
-- Default: no simplification
_ -> CASFrac num' denom'
--------------------------------------------------------------------------------
-- Type-driven reshape (Phase A of the reshape primitive design)
--------------------------------------------------------------------------------
-- | Reshape a CAS value to match the structure implied by the given Type.
--
-- This is the runtime side of the `reshape` primitive: type info comes from
-- a compile-time annotation embedded in the AST, and at evaluation time we
-- structurally rewrite the CASValue to fit that shape.
--
-- The implementation follows the type promotion tower (see type-cas-tower.md):
-- normalize first (which handles tower-level reductions), then recursively
-- adjust coefficients/numerators/denominators to match nested type arguments.
-- For values that cannot be reshaped to the target structure (e.g. a value
-- with free atoms reshaped to TInt), we leave the normalized form as-is —
-- per the "trust the annotation" principle.
casReshapeAs :: Type -> CASValue -> CASValue
casReshapeAs ty v = case ty of
TInt -> casNormalize v
TMathValue -> casNormalize v
TFactor -> casNormalize v
TFrac inner -> reshapeAsFrac inner v
TPoly inner ss -> reshapeAsPoly inner ss v
TTerm inner ss -> reshapeAsTerm inner ss v
_ -> v -- non-CAS types: pass through
-- | Frac inner: keep fraction form (or collapsed Integer when denom=1) and
-- recursively reshape numerator and denominator with the inner type.
reshapeAsFrac :: Type -> CASValue -> CASValue
reshapeAsFrac innerTy v = case casNormalize v of
CASInteger n -> CASInteger n
CASFrac num denom -> casNormalizeFrac (casReshapeAs innerTy num)
(casReshapeAs innerTy denom)
cv -> cv
-- | Poly inner [..]: structural reshape with atom-set separation.
--
-- Atom routing is driven by the whole inner tower — the chain of nested
-- Poly/Term coefficient types, descending through Frac — under the
-- restriction that a nested Poly tower contains AT MOST ONE open symbol
-- set `[..]` (checked statically at annotation and declaration sites,
-- Types.hasAmbiguousOpenTower):
--
-- * The inner tower is all closed (e.g. `Poly (Poly Integer [i]) [x]`,
-- or with an open OUTER set `Poly (Poly Integer [i]) [..]`): atoms
-- listed anywhere in the inner tower go into the coefficient (the deep
-- union, so towers of depth 3+ route correctly), the rest stay at this
-- level.
-- * The inner tower contains the open slot and this level's set is
-- closed (e.g. `Poly (Poly Integer [..]) [i]`): complement split —
-- atoms of this level's closed set stay here, everything else goes
-- into the coefficient.
-- * No routing information at all (base coefficient type such as
-- Integer, or both this level and the inner tower open — excluded by
-- the static check): no separation; just recurse on coefficients
-- (basic widening).
--
-- The inside atoms are folded into the coefficient (which is then
-- recursively reshaped to the inner type, repeating the same routing one
-- level deeper), the outside atoms form the new term's monomial.
-- Note (Phase gamma-prime): the entry `casNormalize v` flattens any nested
-- shape the input may carry, so reshape is a function of the VALUE only —
-- this is what makes the absorption law `casReshapeAs C . casReshapeAs B =
-- casReshapeAs C` hold (D5 coherence). The final grouping then uses the
-- keep-nested normalizer so the structure just built is not re-flattened.
reshapeAsPoly :: Type -> SymbolSet -> CASValue -> CASValue
reshapeAsPoly innerTy outerSS v = case casNormalize v of
CASInteger n -> CASInteger n
CASPoly ts ->
case atomSplit innerTy outerSS of
Just split -> casNormalizePolyWith KeepNested (map (separateTerm innerTy split) ts)
Nothing -> casNormalizePolyWith KeepNested
[CASTerm (casReshapeAs innerTy c) m | CASTerm c m <- ts]
cv -> cv
-- | Term inner [...]: like Poly but expected to be a single-term form.
-- We do not enforce the single-term invariant here; the type checker is
-- responsible for that, and at runtime we just recurse.
reshapeAsTerm :: Type -> SymbolSet -> CASValue -> CASValue
reshapeAsTerm = reshapeAsPoly
-- | Which atoms of a term go into the coefficient at this Poly level.
data AtomSplit
= InsideAtoms [TypeAtom] -- ^ an atom goes inside iff it is listed
| OutsideAtoms [TypeAtom] -- ^ an atom goes inside iff it is NOT listed
-- (complement split for an open inner tower)
-- | Decide the atom split at one Poly level from the inner tower and this
-- level's own symbol set. Returns Nothing when there is no routing
-- information (basic widening).
atomSplit :: Type -> SymbolSet -> Maybe AtomSplit
atomSplit innerTy outerSS =
let (innerClosed, innerOpens) = towerInfo innerTy
in if innerOpens >= 1
then case outerSS of
SymbolSetClosed outerAtoms -> Just (OutsideAtoms outerAtoms)
_ -> Nothing
else if null innerClosed
then Nothing
else Just (InsideAtoms innerClosed)
-- | Collect, over the whole inner tower (nested Poly/Term levels, descending
-- through Frac), the union of closed atom sets and the number of open
-- symbol sets.
towerInfo :: Type -> ([TypeAtom], Int)
towerInfo (TPoly inner ss) = combineSS ss (towerInfo inner)
towerInfo (TTerm inner ss) = combineSS ss (towerInfo inner)
towerInfo (TFrac inner) = towerInfo inner
towerInfo _ = ([], 0)
combineSS :: SymbolSet -> ([TypeAtom], Int) -> ([TypeAtom], Int)
combineSS (SymbolSetClosed atoms) (as, n) = (atoms ++ as, n)
combineSS SymbolSetOpen (as, n) = (as, n + 1)
combineSS (SymbolSetVar _) (as, n) = (as, n)
-- | Split a single CASTerm given the atom split: the inside atoms are
-- multiplied into the coefficient and reshaped to `innTy`; the remaining
-- atoms become the new term's monomial.
separateTerm :: Type -> AtomSplit -> CASTerm -> CASTerm
separateTerm innTy split (CASTerm c mono) =
let (innerMono, outerMono) =
case split of
InsideAtoms atoms -> splitMonomialByAtoms atoms mono
OutsideAtoms atoms -> let (out, inn) = splitMonomialByAtoms atoms mono
in (inn, out)
coeffWithInner =
if null innerMono
then c
else casMult c (CASPoly [CASTerm (CASInteger 1) innerMono])
newCoeff = casReshapeAs innTy coeffWithInner
in CASTerm newCoeff outerMono
-- | Partition a monomial: (atoms-in-inner-set, atoms-not-in-inner-set).
splitMonomialByAtoms :: [TypeAtom] -> Monomial -> (Monomial, Monomial)
splitMonomialByAtoms inAtoms = go [] []
where
go inAcc outAcc [] = (reverse inAcc, reverse outAcc)
go inAcc outAcc (e@(sym, _) : rest)
| symbolInAtomSet sym inAtoms = go (e : inAcc) outAcc rest
| otherwise = go inAcc (e : outAcc) rest
-- | Decide whether a SymbolExpr matches any TypeAtom in the inner set.
-- Matches `Symbol _ name _` against `TANameAtom name` (e.g. atom `i`).
-- Matches `Apply1..4 fn _` against `TAApplyAtom name _` by function name
-- (e.g. atom `sin x` for any single application of `sin`). The argument
-- structure inside the TypeAtom is not currently checked.
symbolInAtomSet :: SymbolExpr -> [TypeAtom] -> Bool
symbolInAtomSet sym = any (matches sym)
where
matches (Symbol _ name _) (TANameAtom n) = name == n
matches (Apply1 fn _) (TAApplyAtom n _) = applyFnName fn == Just n
matches (Apply2 fn _ _) (TAApplyAtom n _) = applyFnName fn == Just n
matches (Apply3 fn _ _ _) (TAApplyAtom n _) = applyFnName fn == Just n
matches (Apply4 fn _ _ _ _) (TAApplyAtom n _) = applyFnName fn == Just n
matches (FunctionData fn _) (TAApplyAtom n _) = applyFnName fn == Just n
matches _ _ = False
applyFnName (CASFactor (Symbol _ n _)) = Just n
applyFnName (CASPoly [CASTerm (CASInteger 1) [(Symbol _ n _, 1)]]) = Just n
applyFnName _ = Nothing
-- | Simplify polynomial division by extracting common monomial GCD
simplifyPolyDiv :: [CASTerm] -> [CASTerm] -> ([CASTerm], [CASTerm])
simplifyPolyDiv [] ts2 = ([], ts2)
simplifyPolyDiv ts1 [] = (ts1, [])
simplifyPolyDiv ts1 ts2 =
let gcdTerm = casTermsGcd (ts1 ++ ts2)
in (map (`divideTermBy` gcdTerm) ts1, map (`divideTermBy` gcdTerm) ts2)
-- | Divide a Term by a monomial denominator (single-term Poly's coef and mono).
-- Used for Laurent absorption: `(coef * mono) / (denomCoef * denomMono)` becomes
-- `(coef / denomCoef) * (mono - denomMono)` where exponents are subtracted.
divTermByMonomial :: CASValue -> Monomial -> CASTerm -> CASTerm
divTermByMonomial denomCoef denomMono (CASTerm c m) =
CASTerm (casNormalizeFrac c denomCoef) (subtractMonomial m denomMono)
-- | Subtract one monomial from another (decrement exponents of shared symbols,
-- introduce negative exponents for symbols only in the divisor). Zero-exponent
-- entries are left in; `normalizeTermMonomial` will filter them out later.
subtractMonomial :: Monomial -> Monomial -> Monomial
subtractMonomial nums denoms = foldr subOne nums denoms
where
subOne (sym, denomExp) acc =
case lookup sym acc of
Just numExp ->
(sym, numExp - denomExp) : filter ((/= sym) . fst) acc
Nothing -> (sym, -denomExp) : acc
--------------------------------------------------------------------------------
-- GCD Operations
--------------------------------------------------------------------------------
-- | GCD of two CASValues (for coefficient reduction)
-- Initial implementation: only handles CASInteger, others return 1
casGcd :: CASValue -> CASValue -> CASValue
casGcd (CASInteger a) (CASInteger b) = CASInteger (gcd a b)
casGcd _ _ = CASInteger 1 -- Fallback: GCD = 1 for other coefficient types
-- | Compute the GCD of a list of terms (coefficient GCD + monomial GCD)
casTermsGcd :: [CASTerm] -> CASTerm
casTermsGcd [] = CASTerm (CASInteger 1) []
casTermsGcd [t] = t
casTermsGcd terms = foldl1 termGcd terms
where
termGcd (CASTerm c1 m1) (CASTerm c2 m2) =
CASTerm (casGcd c1 c2) (monoGcd m1 m2)
-- | Reduce a Poly/Poly fraction by the univariate polynomial GCD over Q,
-- e.g. (x^2 - 1)/(x - 1) -> (x + 1)/1.
--
-- Stage-1 scope (design/type-cas-tower-implementation.md section 7): both
-- term lists must be univariate in the SAME single symbol with positive
-- exponents (the common monomial content has already been divided out by
-- 'simplifyPolyDiv', but per-side Laurent exponents may remain — those
-- bail out) and integer or integer-fraction coefficients. Anything else
-- returns Nothing and the fraction is left untouched.
--
-- The reduced pair is rescaled by a COMMON factor — integer coefficients,
-- joint content 1, positive leading denominator coefficient — so the
-- fraction's value is preserved exactly.
univariateGcdReduce :: [CASTerm] -> [CASTerm] -> Maybe ([CASTerm], [CASTerm])
univariateGcdReduce ts1 ts2
| null ts1 || null ts2 = Nothing
| otherwise = do
sym <- singleCommonSymbol
p1 <- toDense sym ts1
p2 <- toDense sym ts2
let g = polyGcdQ p1 p2
if length g < 2 -- constant gcd: nothing to reduce
then Nothing
else do
q1 <- exactDivQ p1 g
q2 <- exactDivQ p2 g
let l = foldr (lcm . denominator) 1 (q1 ++ q2)
i1 = map (numerator . (* (l % 1))) q1
i2 = map (numerator . (* (l % 1))) q2
c0 = foldr gcd 0 (i1 ++ i2)
c = if c0 == 0 then 1 else c0
d = case i2 of (x : _) | x < 0 -> negate c
_ -> c
return (fromDense sym (map (`div` d) i1), fromDense sym (map (`div` d) i2))
where
-- Euclid on dense rationals is cheap for the degrees CAS code meets;
-- the cutoff only guards against pathological inputs.
maxGcdDegree :: Integer
maxGcdDegree = 200
singleCommonSymbol =
case nub [ s | CASTerm _ m <- ts1 ++ ts2, (s, _) <- m ] of
[s] -> Just s
_ -> Nothing
coefToRational (CASInteger n) = Just (n % 1)
coefToRational (CASFrac (CASInteger a) (CASInteger b))
| b /= 0 = Just (a % b)
coefToRational _ = Nothing
termExponent sym (CASTerm _ m) = case m of
[] -> Just 0
[(s, e)] | s == sym && e > 0 -> Just e
_ -> Nothing
-- Dense, highest-degree-first coefficient list over Q.
toDense sym ts = do
pairs <- mapM (\t@(CASTerm c _) ->
(,) <$> termExponent sym t <*> coefToRational c) ts
let deg = maximum (map fst pairs)
if deg > maxGcdDegree
then Nothing
else Just [ sum [ c | (e, c) <- pairs, e == d ] | d <- [deg, deg-1 .. 0] ]
fromDense sym cs =
[ CASTerm (CASInteger c) (if e == 0 then [] else [(sym, e)])
| (e, c) <- zip [toInteger (length cs) - 1, toInteger (length cs) - 2 .. 0] cs
, c /= 0 ]
trim = dropWhile (== 0)
-- Position-preserving long division (no mid-loop trimming: a zero that
-- appears at the head after cancellation is a zero QUOTIENT coefficient,
-- not a shorter polynomial). One quotient coefficient per step.
polyDivModQ x0 y0 =
let y = trim y0
m = length y
go r | length r < m = ([], r)
| otherwise = case (r, y) of
(rh : rt, yh : yt) ->
let k = rh / yh
r' = zipWith (-) rt
(map (* k) (yt ++ replicate (length r - m) 0))
(qs, rest) = go r'
in (k : qs, rest)
_ -> error "polyDivModQ: division by the zero polynomial"
in go (trim x0)
polyGcdQ a b = go (trim a) (trim b)
where
go x [] = monic x
go x y = go y (trim (snd (polyDivModQ x y)))
monic [] = []
monic xs@(x0:_) = map (/ x0) xs
exactDivQ x y =
let (q, r) = polyDivModQ x y
in if all (== 0) r then Just q else Nothing
-- | Reduce a Poly/Poly fraction by the multivariate polynomial GCD over the
-- rationals (subresultant PRS) — stage 2 of the fraction reduction
-- (design/cas-simplification.md G1).
--
-- Every atom (symbols, symbolic applications such as 'cos θ, quotes,
-- function symbols) is treated uniformly as a variable, so the reducer
-- covers Schwarzschild/T2/thurston-style cancellations with no extra
-- machinery. Fail-open: any input outside the supported shape (Laurent
-- exponents, non-rational coefficients, sizes beyond the guards) returns
-- Nothing and the fraction is left untouched.
--
-- Value preservation: both sides are scaled by ONE common denominator-
-- clearing factor, divided exactly by the same gcd, and finally rescaled by
-- a common integer content and sign, so the fraction's value never changes.
multivariateGcdReduce :: [CASTerm] -> [CASTerm] -> Maybe ([CASTerm], [CASTerm])
multivariateGcdReduce ts1 ts2
| null ts1 || null ts2 = Nothing
| otherwise = do
rs1 <- mapM ratTerm ts1
rs2 <- mapM ratTerm ts2
-- Laurent exponents are out of scope (handled by the monomial layer).
ensure (all (all ((> 0) . snd) . snd) (rs1 ++ rs2))
let atoms = sortBy (comparing show)
(nub [ s | (_, m) <- rs1 ++ rs2, (s, _) <- m ])
-- The univariate reducer owns the single-symbol case.
ensure (length atoms >= 2 && length atoms <= mpMaxAtoms)
let nAtoms = length atoms
scale = foldr (lcm . denominator . fst) 1 (rs1 ++ rs2)
zeroV = replicate nAtoms 0
oneMP = [(1, zeroV)]
toVec m = [ sum [ e | (s, e) <- m', s == a ] | a <- atoms ]
where m' = foldMonomialSymbols m
conv rs = mpNorm [ (numerator (c * (scale % 1)), toVec m) | (c, m) <- rs ]
p = conv rs1
q = conv rs2
-- Descending graded-lexicographic order over the fixed atom list;
-- a proper monomial order, so leading-term exact division works.
mpGrlex a b = compare (sum a, a) (sum b, b)
mpNorm ts =
[ (c, v)
| grp@((_, v) : _) <- groupBy ((==) `on` snd)
(sortBy (flip mpGrlex `on` snd) ts)
, let c = sum (map fst grp)
, c /= 0 ]
mpNeg = map (\(c, v) -> (negate c, v))
mpAdd a b = mpNorm (a ++ b)
mpSub a b = mpAdd a (mpNeg b)
mpMul a b = mpNorm [ (c1 * c2, zipWith (+) v1 v2)
| (c1, v1) <- a, (c2, v2) <- b ]
mpMulTerm (c, v) b = mpNorm [ (c * c2, zipWith (+) v v2) | (c2, v2) <- b ]
mpPow b n = foldr mpMul oneMP (replicate (fromInteger n) b)
mpMaxTotalDeg f = maximum (0 : [ sum v | (_, v) <- f ])
mpDegIn i f = maximum (0 : [ v !! i | (_, v) <- f ])
mpPresent f = [ i | i <- [0 .. nAtoms - 1]
, any (\(_, v) -> v !! i > 0) f ]
mpIntContent f = foldr (gcd . fst) 0 f
-- Positive normalization: strip the integer content, make the
-- leading (grlex) coefficient positive.
mpPosNorm [] = []
mpPosNorm f =
let c0 = mpIntContent f
s = case f of ((c, _) : _) | c < 0 -> -1
_ -> 1
d = s * c0
in map (\(c, v) -> (c `div` d, v)) f
-- Leading-term exact division; Nothing when not exact.
mpDivExact _ [] = Nothing
mpDivExact a0 b@((bc, bv) : _) = go a0 []
where
go [] acc = Just (mpNorm acc)
go a@((ac, av) : _) acc
| all (>= 0) dv && ac `mod` bc == 0 =
let qt = (ac `div` bc, dv)
in go (mpSub a (mpMulTerm qt b)) (qt : acc)
| otherwise = Nothing
where dv = zipWith (-) av bv
zeroAtI i v = take i v ++ [0] ++ drop (i + 1) v
-- Univariate view in atom i: (degree, coefficient poly without i),
-- degrees descending.
mpUniView i f =
[ (v0 !! i, mpNorm [ (c, zeroAtI i v) | (c, v) <- grp ])
| grp@((_, v0) : _) <- groupBy ((==) `on` ((!! i) . snd))
(sortBy (flip compare `on` ((!! i) . snd)) f) ]
mpCoefAt i d f = mpNorm [ (c, zeroAtI i v) | (c, v) <- f, v !! i == d ]
mpLeadCoef i f = mpCoefAt i (mpDegIn i f) f
mpShift i k = map (\(c, v) ->
(c, take i v ++ [v !! i + k] ++ drop (i + 1) v))
-- Content and primitive part with respect to atom i.
mpContentI i f = goC (map snd (mpUniView i f))
where
goC [] = Just oneMP
goC [c0] = Just (mpPosNorm c0)
goC (c0 : cs) = do
rest <- goC cs
if rest == oneMP then Just oneMP else mpGcdM (mpPosNorm c0) rest
-- Knuth's division-free pseudo-remainder of f1 by f2 in atom i:
-- the remainder of lc(f2)^(delta+1) * f1 divided by f2.
mpPrem i f1 f2 =
let dq = mpDegIn i f2
lcq = mpLeadCoef i f2
step r k =
let ck = mpCoefAt i (dq + k) r
r' = mpSub (mpMul lcq r) (mpMul (mpShift i k ck) f2)
in r'
go r k | length r > mpPrsTermCap = Nothing
| k < 0 = Just r
| otherwise = go (step r k) (k - 1)
in go f1 (mpDegIn i f1 - dq)
-- Subresultant PRS on primitive parts; returns the gcd of the
-- primitive parts (a unit when they are coprime in atom i).
mpPrsGcd i f1 f2
| mpDegIn i f1 < mpDegIn i f2 = mpPrsGcd i f2 f1
| mpDegIn i f2 == 0 = Just oneMP
| otherwise = loop f1 f2 oneMP oneMP
where
loop a b g h = do
ensure (length a <= mpPrsTermCap && length b <= mpPrsTermCap)
let delta = mpDegIn i a - mpDegIn i b
r <- mpPrem i a b
if null r
then do c <- mpContentI i b
mpDivExact b c
else if mpDegIn i r == 0
then Just oneMP
else do
b' <- mpDivExact r (mpMul g (mpPow h delta))
let g' = mpLeadCoef i b
h' <- case delta of
0 -> Just h
1 -> Just g'
_ -> mpDivExact (mpPow g' delta)
(mpPow h (delta - 1))
loop b b' g' h'
-- Multivariate gcd, recursing on the set of present atoms.
mpGcdM a b
| null a = Just (mpPosNorm b)
| null b = Just (mpPosNorm a)
| otherwise =
case mpPresent a ++ mpPresent b of
[] -> Just [(gcd (mpIntContent a) (mpIntContent b), zeroV)]
idxs -> do
let i = pickVar (nub idxs)
ca <- mpContentI i a
pa <- mpDivExact a ca
cb <- mpContentI i b
pb <- mpDivExact b cb
c <- mpGcdM ca cb
g <- mpPrsGcd i pa pb
Just (mpPosNorm (mpMul c g))
where
pickVar is = snd (minimum
[ (mpDegIn i a + mpDegIn i b, i) | i <- is ])
-- Deterministic coprimality prefilter: evaluate both sides at a
-- fixed prime point; a common divisor g must satisfy g(pt) | gcd
-- of the evaluations, so gcd 1 certifies that any common divisor
-- is a unit at that point (heuristic skip; fail-open, and
-- deterministic because the points are fixed).
mpEval pt f = sum [ c * product (zipWith (^) pt v) | (c, v) <- f ]
evalCoprime =
any certify [ [3,5,7,11,13,17,19,23], [29,31,37,41,43,47,53,59] ]
where
certify pt =
let a = mpEval (take nAtoms pt) p
b = mpEval (take nAtoms pt) q
in a /= 0 && b /= 0 && gcd a b == 1
ensure (length p <= mpMaxTerms && length q <= mpMaxTerms)
ensure (mpMaxTotalDeg p <= mpMaxDegree && mpMaxTotalDeg q <= mpMaxDegree)
ensure (not (null (mpPresent p `intersect` mpPresent q)))
ensure (not evalCoprime)
g <- mpGcdM p q
ensure (mpMaxTotalDeg g > 0)
pR <- mpDivExact p g
qR <- mpDivExact q g
-- Common rescale of the reduced pair (value-preserving): joint integer
-- content out, denominator's leading coefficient positive.
let cJ = gcd (mpIntContent pR) (mpIntContent qR)
sg = case qR of ((c, _) : _) | c < 0 -> -1
_ -> 1
d = sg * (if cJ == 0 then 1 else cJ)
toTerms f = [ CASTerm (CASInteger (c `div` d))
[ (a, e) | (a, e) <- zip atoms v, e /= 0 ]
| (c, v) <- f ]
return (toTerms pR, toTerms qR)
where
ensure b = if b then Just () else Nothing
ratTerm (CASTerm c m) = (\r -> (r, m)) <$> ratCoef c
ratCoef (CASInteger n) = Just (n % 1)
ratCoef (CASFrac (CASInteger a) (CASInteger b)) | b /= 0 = Just (a % b)
ratCoef _ = Nothing
mpMaxAtoms = 8
mpMaxTerms = 200
mpMaxDegree = 60 :: Integer
mpPrsTermCap = 2000
-- | GCD of two monomials: take minimum exponent for each shared symbol
monoGcd :: Monomial -> Monomial -> Monomial
monoGcd [] _ = []
monoGcd ((sym, expo):rest) mono =
case lookup sym mono of
Just exp' -> (sym, min expo exp') : monoGcd rest mono
Nothing -> monoGcd rest mono
-- | Divide a term by another term (for GCD reduction)
divideTermBy :: CASTerm -> CASTerm -> CASTerm
divideTermBy (CASTerm coeff1 mono1) (CASTerm coeff2 mono2) =
CASTerm (divCoeff coeff1 coeff2) (divMono mono1 mono2)
where
divCoeff (CASInteger a) (CASInteger b) = CASInteger (a `div` b)
divCoeff a _ = a -- Fallback: no division for other types
divMono m [] = m
divMono m ((sym, expo):rest) =
let m' = map (\(s, e) -> if s == sym then (s, e - expo) else (s, e)) m
in divMono m' rest
--------------------------------------------------------------------------------
-- Pattern Synonyms for CASValue
--------------------------------------------------------------------------------
-- | Pattern for zero value
pattern CASZero :: CASValue
pattern CASZero = CASInteger 0
-- | Pattern for a single symbol: sym → 1 * sym^1
pattern CASSingleSymbol :: SymbolExpr -> CASValue
pattern CASSingleSymbol sym = CASPoly [CASTerm (CASInteger 1) [(sym, 1)]]
-- | Pattern for a single term: coeff * mono
pattern CASSingleTerm :: Integer -> Monomial -> CASValue
pattern CASSingleTerm coeff mono = CASPoly [CASTerm (CASInteger coeff) mono]
--------------------------------------------------------------------------------
-- Pattern Matching (control-egison matchers)
--------------------------------------------------------------------------------
-- | Matcher for CASValue
data CASM = CASM
instance Matcher CASM CASValue
-- | Matcher for CASTerm
data CASTermM = CASTermM
instance Matcher CASTermM CASTerm
-- | Matcher for SymbolExpr (CAS version)
data CASSymbolM = CASSymbolM
instance Matcher CASSymbolM SymbolExpr
-- | Match a term and extract its coefficient and monomial
casTerm' :: Pattern (PP CASValue, PP Monomial) CASTermM CASTerm (CASValue, Monomial)
casTerm' _ _ (CASTerm coeff mono) = pure (coeff, mono)
-- | Matcher decomposition for casTerm' pattern
casTerm'M :: CASTermM -> CASTerm -> (CASM, Multiset (CASSymbolM, Eql))
casTerm'M CASTermM _ = (CASM, Multiset (CASSymbolM, Eql))
casTermM :: CASTermM -> CASTerm -> (CASM, Multiset (CASSymbolM, Eql))
casTermM CASTermM _ = (CASM, Multiset (CASSymbolM, Eql))
-- | Match a symbol and extract its name
casSymbol :: Pattern (PP String) CASSymbolM SymbolExpr String
casSymbol _ _ (Symbol _ name []) = pure name
casSymbol _ _ _ = mzero
casSymbolM :: CASSymbolM -> p -> Eql
casSymbolM CASSymbolM _ = Eql
-- | Match a function and extract its name and arguments
casFunc :: Pattern (PP CASValue, PP [CASValue])
CASSymbolM SymbolExpr (CASValue, [CASValue])
casFunc _ _ (FunctionData name args) = pure (name, args)
casFunc _ _ _ = mzero
casFuncM :: CASSymbolM -> SymbolExpr -> (CASM, List CASM)
casFuncM CASSymbolM _ = (CASM, List CASM)
-- | Match Apply1 and extract function name, WHNF, and argument
casApply1 :: Pattern (PP String, PP WHNFData, PP CASValue) CASSymbolM SymbolExpr (String, WHNFData, CASValue)
casApply1 _ _ (Apply1 (CASSingleSymbol (QuoteFunction fnWhnf)) a1) =
case prettyFunctionName fnWhnf of
Just fn -> pure (fn, fnWhnf, a1)
Nothing -> mzero
casApply1 _ _ _ = mzero
casApply1M :: CASSymbolM -> p -> (Eql, Something, CASM)
casApply1M CASSymbolM _ = (Eql, Something, CASM)
-- | Match Apply2 and extract function name, WHNF, and arguments
casApply2 :: Pattern (PP String, PP WHNFData, PP CASValue, PP CASValue) CASSymbolM SymbolExpr (String, WHNFData, CASValue, CASValue)
casApply2 _ _ (Apply2 (CASSingleSymbol (QuoteFunction fnWhnf)) a1 a2) =
case prettyFunctionName fnWhnf of
Just fn -> pure (fn, fnWhnf, a1, a2)
Nothing -> mzero
casApply2 _ _ _ = mzero
casApply2M :: CASSymbolM -> p -> (Eql, Something, CASM, CASM)
casApply2M CASSymbolM _ = (Eql, Something, CASM, CASM)
-- | Match Apply3 and extract function name, WHNF, and arguments
casApply3 :: Pattern (PP String, PP WHNFData, PP CASValue, PP CASValue, PP CASValue) CASSymbolM SymbolExpr (String, WHNFData, CASValue, CASValue, CASValue)
casApply3 _ _ (Apply3 (CASSingleSymbol (QuoteFunction fnWhnf)) a1 a2 a3) =
case prettyFunctionName fnWhnf of
Just fn -> pure (fn, fnWhnf, a1, a2, a3)
Nothing -> mzero
casApply3 _ _ _ = mzero
casApply3M :: CASSymbolM -> p -> (Eql, Something, CASM, CASM, CASM)
casApply3M CASSymbolM _ = (Eql, Something, CASM, CASM, CASM)
-- | Match Apply4 and extract function name, WHNF, and arguments
casApply4 :: Pattern (PP String, PP WHNFData, PP CASValue, PP CASValue, PP CASValue, PP CASValue) CASSymbolM SymbolExpr (String, WHNFData, CASValue, CASValue, CASValue, CASValue)
casApply4 _ _ (Apply4 (CASSingleSymbol (QuoteFunction fnWhnf)) a1 a2 a3 a4) =
case prettyFunctionName fnWhnf of
Just fn -> pure (fn, fnWhnf, a1, a2, a3, a4)
Nothing -> mzero
casApply4 _ _ _ = mzero
casApply4M :: CASSymbolM -> p -> (Eql, Something, CASM, CASM, CASM, CASM)
casApply4M CASSymbolM _ = (Eql, Something, CASM, CASM, CASM, CASM)
-- | Match Quote and extract the inner CASValue
casQuote :: Pattern (PP CASValue) CASSymbolM SymbolExpr CASValue
casQuote _ _ (Quote m) = pure m
casQuote _ _ _ = mzero
-- | Match Quote and extract the negated inner CASValue
casNegQuote :: Pattern (PP CASValue) CASSymbolM SymbolExpr CASValue
casNegQuote _ _ (Quote m) = pure (casNegate m)
casNegQuote _ _ _ = mzero
casNegQuoteM :: CASSymbolM -> p -> CASM
casNegQuoteM CASSymbolM _ = CASM
-- | Match QuoteFunction and extract function name and WHNF
casQuoteFunction :: Pattern (PP String, PP WHNFData) CASSymbolM SymbolExpr (String, WHNFData)
casQuoteFunction _ _ (QuoteFunction whnf) = case prettyFunctionName whnf of
Just name -> pure (name, whnf)
Nothing -> mzero
casQuoteFunction _ _ _ = mzero
casQuoteFunctionM :: CASSymbolM -> p -> Eql
casQuoteFunctionM CASSymbolM _ = Eql
-- | Match equal monomial (checks if two monomials are equal, handling sign)
casEqualMonomial :: Pattern (PP Integer, PP Monomial) (Multiset (CASSymbolM, Eql)) Monomial (Integer, Monomial)
casEqualMonomial (_, VP xs) _ ys = case casIsEqualMonomial xs ys of
Just sgn -> pure (sgn, xs)
Nothing -> mzero
casEqualMonomial _ _ _ = mzero
casEqualMonomialM :: Multiset (CASSymbolM, Eql) -> p -> (Eql, Multiset (CASSymbolM, Eql))
casEqualMonomialM (Multiset (CASSymbolM, Eql)) _ = (Eql, Multiset (CASSymbolM, Eql))
-- | Check if two monomials are equal, returning sign if so
casIsEqualMonomial :: Monomial -> Monomial -> Maybe Integer
casIsEqualMonomial xs ys =
match dfs (xs, ys) (Multiset (CASSymbolM, Eql), Multiset (CASSymbolM, Eql))
[ [mc| ((casQuote $s, $n) : $xss, (casNegQuote #s, #n) : $yss) ->
case casIsEqualMonomial xss yss of
Nothing -> Nothing
Just sgn -> return (if even n then sgn else - sgn) |]
, [mc| (($x, $n) : $xss, (#x, #n) : $yss) -> casIsEqualMonomial xss yss |]
, [mc| ([], []) -> return 1 |]
, [mc| _ -> Nothing |]
]
-- | Match zero CASValue
casZero :: Pattern () CASM CASValue ()
casZero _ _ CASZero = pure ()
casZero _ _ (CASPoly []) = pure () -- Empty polynomial is also zero
casZero _ _ _ = mzero
casZeroM :: CASM -> p -> ()
casZeroM CASM _ = ()
-- | Match a single term in CASValue and extract coefficient, denominator coefficient, and monomial
casSingleTerm :: Pattern (PP Integer, PP Integer, PP Monomial) CASM CASValue (Integer, Integer, Monomial)
casSingleTerm _ _ (CASFrac (CASPoly [CASTerm (CASInteger c) mono]) (CASPoly [CASTerm (CASInteger c2) []])) = pure (c, c2, mono)
casSingleTerm _ _ (CASFrac (CASSingleTerm c mono) (CASInteger c2)) = pure (c, c2, mono)
casSingleTerm _ _ (CASPoly [CASTerm (CASInteger c) mono]) = pure (c, 1, mono)
casSingleTerm _ _ (CASInteger n) = pure (n, 1, []) -- Integer is a single term with empty monomial
casSingleTerm _ _ _ = mzero
casSingleTermM :: CASM -> p -> (Eql, Eql, Multiset (CASSymbolM, Eql))
casSingleTermM CASM _ = (Eql, Eql, Multiset (CASSymbolM, Eql))
-- | ValuePattern instance for CASM
instance ValuePattern CASM CASValue where
value e () CASM v = if e == v then pure () else mzero
-- | ValuePattern instance for CASSymbolM
instance ValuePattern CASSymbolM SymbolExpr where
value e () CASSymbolM v = if e == v then pure () else mzero