Make restrictions more robust

Author: Andrew Luo

sage_acsv/asymptotics.py

@@ -129,15 +129,18 @@
 from sage.rings.asymptotic.asymptotic_ring import AsymptoticRing
 from sage.rings.ideal import Ideal
 from sage.rings.imaginary_unit import I
+from sage.rings.fraction_field import FractionField
 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
 from sage.rings.qqbar import AA, QQbar
 from sage.rings.rational_field import QQ
+from sage.schemes.affine.affine_space import AffineSpace
 from sage.symbolic.constants import pi
 from sage.symbolic.ring import SR
 
 from sage_acsv.kronecker import _kronecker_representation
 from sage_acsv.helpers import (
     ACSVException,
+    Restriction,
     is_contributing,
     compute_newton_series,
     rational_function_reduce,
@@ -148,7 +151,7 @@
 from sage_acsv.debug import Timer, acsv_logger
 from sage_acsv.settings import ACSVSettings
 from sage_acsv.whitney import whitney_stratification
-from sage_acsv.groebner import compute_primary_decomposition, compute_saturation
+from sage_acsv.groebner import compute_primary_decomposition, compute_saturation, compute_radical
 
 
 # we need to monkeypatch a function from the asymptotics module such that creating
@@ -778,18 +781,28 @@ def compute_asymptotic_contribution(
     unit, all_factors, all_multiplicities = _get_factorization(R(H))
     unit, all_factors, all_multiplicities = SR(unit), [SR(f) for f in all_factors], [SR(f) for f in all_multiplicities]
 
+    fvs = list(set([v for ri in r for v in SR(ri).variables()] + [v for cpi in cp for v in SR(cpi).variables()]))
+    Rf = PolynomialRing(QQ, len(fvs), fvs)
+
     asymptotics = []
+    restrictions = []
 
     # Since critical point might be symbolic, we need to consider all flats it may lie on
-    for factor_multiplicities in Combinations(zip(all_factors, all_multiplicities)):
-        if not factor_multiplicities:
+    for combo in Combinations(list(range(len(all_factors)))):
+        if not combo:
             continue
-        factors, multiplicities = zip(*factor_multiplicities)
+        factors = [all_factors[i] for i in combo]
+        multiplicities = [all_multiplicities[i] for i in combo] 
+        
+        valid_set = Restriction(Rf)
+        valid_set.setminus(
+            [all_factors[j] for j in range(len(all_factors)) if j not in combo]
+        )
+
         # Determine whether it's possible for the critical point vanish in this
         # combination of factors
         if  any([f.subs(subs_dict) != 0 for f in factors]):
-            # TODO - check if some parameters allow this to hold in certain conditions
-            continue
+            valid_set.intersection([f.subs(subs_dict) for f in factors])
 
         s = len(factors)
         normals = matrix(
@@ -799,6 +812,12 @@ def compute_asymptotic_contribution(
             raise ACSVException(
                 "Not a transverse intersection. Cannot deal with this case."
             )
+        else:
+            try:
+                Id_minors = [Rf(f) for f in normals.minors(s)]
+                valid_set.setminus(Id_minors)
+            except:
+                valid_set = None
 
         # Step 2: Find the locally parametrizing coordinates of the point pt
         # Since we have d variables and s factors, there should be d-s of these
@@ -815,6 +834,11 @@ def compute_asymptotic_contribution(
             if Jac.determinant() != 0:
                 break
 
+            try:
+                valid_set = valid_set.setminus([Jac.determinant()])
+            except:
+                valid_set = None
+
             acsv_logger.info("Variables do not parametrize, shuffling")
             vs_r_cp = list(zip(vs, r, cp))
             shuffle(vs_r_cp)  # shuffle mutates the list
@@ -823,16 +847,17 @@ def compute_asymptotic_contribution(
             raise ACSVException("Cannot find parametrizing set.")
 
         expansion, constant_squared, base, exponent, s = _compute_asm_quantity(
-            G, H, vs, cp, r, subs_dict, unit ,factors, multiplicities, 1
+            G, H, vs, cp, r, subs_dict, unit,factors, multiplicities, 1
         )
         constant = constant_squared.sqrt()
 
         n = SR.var("n")
         asymptotics.append(base**n * n**exponent * (pi ** (s - d)).sqrt() * constant * expansion)
+        restrictions.append(valid_set)
 
     # TODO: If input parameters are rational, do some extra checking of validity
 
-    return asymptotics
+    return restrictions, asymptotics
 
 def _compute_asm_quantity(G, H, vs, cp, r, subs_dict, unit, factors, multiplicities, expansion_precision):
     s = len(factors)

sage_acsv/helpers.py

@@ -19,12 +19,16 @@
     MonomialGrowthGroup,
 )
 from sage.rings.ideal import Ideal
+from sage.rings.fraction_field import FractionField
+from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
 from sage.rings.qqbar import AA, AlgebraicNumber, QQbar
 from sage.rings.rational_field import QQ
+from sage.schemes.affine.affine_space import AffineSpace
 from sage.symbolic.expression import Expression
 from sage.symbolic.ring import SymbolicRing, SR
 from sage.symbolic.operators import add_vararg
 from sage.symbolic.constants import pi
+from sage_acsv.groebner import compute_radical, compute_saturation
 
 
 @dataclass
@@ -57,7 +61,6 @@ class Term:
     def __lt__(self, other):
         return (self.base, self.power) < (other.base, other.power)
 
-
 def collapse_zero_part(algebraic_number: AlgebraicNumber) -> AlgebraicNumber:
     if algebraic_number.real().is_zero():
         algebraic_number = QQbar(algebraic_number.imag()) * QQbar(-1).sqrt()
@@ -567,6 +570,83 @@ def get_expansion_terms(
 
     return decomposed_terms
 
+class Restriction():
+    def __init__(self, base_ring):
+        if len(base_ring.gens()) == 0:
+            self.base_ring = SR
+            self.base_field = SR
+            self.inclusion = []
+            self.exclusion = []
+            return
+
+        self.base_ring = base_ring
+        self.base_field = FractionField(base_ring)
+        self.inclusion = Ideal(base_ring.zero())
+        self.exclusion = Ideal(base_ring.one())
+
+    def intersection(self, eqns):
+        if self.base_ring != SR:
+            try:
+                eqns = [self.base_field(f) for f in eqns]
+                denoms = [f.denominator() for f in eqns]
+                Id = compute_radical(Ideal([self.base_ring.zero()] + [self.base_ring(f * prod(denoms)) for f in eqns]))
+                self.inclusion += Id
+            except:
+                self.inclusion = eqns + [SR(f) for f in self.inclusion.gens()]
+                self.exclusion = [SR(f) for f in self.exclusion.gens()]
+                self.base_ring = SR
+                self.base_field = SR
+        else:
+            self.inclusion += eqns
+
+        return self
+    
+    def setminus(self, eqns):
+        if self.base_ring != SR:
+            try:
+                eqns = [self.base_field(f) for f in eqns]
+                denoms = [f.denominator() for f in eqns]
+                Id = compute_radical(Ideal([self.base_ring.zero()] + [self.base_ring(f * prod(denoms)) for f in eqns]))
+                self.exclusion = self.exclusion.intersection(Id)
+            except:
+                self.inclusion = [SR(f) for f in self.inclusion.gens()]
+                self.exclusion = eqns + [SR(f) for f in self.exclusion.gens()]
+                self.base_ring = SR
+                self.base_field = SR
+        else:
+            self.exclusion += eqns
+
+        return self
+    
+    def simplify(self):
+        if self.base_ring == SR:
+            return
+        self.inclusion = compute_saturation(self.inclusion, self.exclusion)
+        self.inclusion = compute_radical(self.inclusion)
+        self.exclusion = compute_radical(self.exclusion)
+
+    def contains(self, other : Restriction):
+        if self.base_ring == SR:
+            return False
+        
+        self.simplify()
+        other.simplify()
+        for f in self.inclusion.gens():
+            if f not in other.inclusion:
+                return False
+        
+        exclusion_set = compute_radical(self.exclusion + other.inclusion)
+        for f in other.exclusion.gens():
+            if f not in exclusion_set:
+                return False
+        
+        return True
+    
+    def is_empty(self):
+        return self.contains(Restriction(self.base_ring).intersection([1]))
+    
+    def __str__(self):
+        return f"Quasi-Affine variety defined by X - Y where X is \n {self.inclusion} and Y is \n {self.exclusion}"
 
 class ACSVException(Exception):
     def __init__(self, message, retry=False):