feat: improve polyrith by testing for membership in the radical (#7790

)

In @hrmacbeth 's tutorial on `polyrith`, there were examples of problems that it could almost solve, but failed. The goal was not expressible as a linear combination of the hypotheses but a power of the goal was. 

```lean
example (x y z : ℚ) (h : x = y) (h2 : x * y = 0) : x + y*z = 0 :=
sorry
```

Mathematically, `x+y*z` is in the radical of the ideal generated by `x-y, x*y`. There's a "standard trick" for testing membership in the radical without a search for the proper exponent: see e.g. section 2.2 of [arxiv.org/pdf/1007.3615.pdf](https://arxiv.org/pdf/1007.3615.pdf) or 4.2 Prop 8 of Cox, Little, O'Shea.

This PR implements the trick in the Sage call made by `polyrith`. When the power returned is `n > 1`, we use `linear_combination (exp := n)` to check the certificate (#7789 ).

The `polyrith` test infrastructure still needs to be ported from mathlib3. All tests in the test file succeed when they are uncommented. A future PR will restore the old test suite.


Co-authored-by: Rob Lewis <[email protected]>
Co-authored-by: Mario Carneiro <[email protected]>

Mathlib/Tactic/Polyrith.lean

@@ -33,6 +33,11 @@ It then calls a Python file that uses the SageMath API to compute the coefficien
 coefficients are then sent back to Lean, which parses them into pexprs. The information is then
 given to the `linear_combination` tactic, which completes the process by checking the certificate.
 
+In fact, `polyrith` uses Sage to test for membership in the *radical* of the ideal.
+This means it searches for a linear combination of hypotheses that add up to a *power* of the goal.
+When this power is not 1, it uses the `(exp := n)` feature of `linear_combination` to report the
+certificate.
+
 `polyrith` calls an external python script `scripts/polyrith_sage.py`. Because this is not a Lean
 file, changes to this script may not be noticed during Lean compilation if you have already
 generated olean files. If you are modifying this python script, you likely know what you're doing;
@@ -250,14 +255,24 @@ instance : FromJson Poly where
         mon := mon.mul' (.pow' (← fromJson? (← j.getArrVal? 0)) (← fromJson? (← j.getArrVal? 1)))
       pure mon
 
+/-- A schema for the data reported by the Sage calculation -/
+structure SageCoeffAndPower where
+  /-- The function call produces an array of polynomials
+  parallel to the input list of hypotheses. -/
+  coeffs : Array Poly
+  /-- Sage produces an exponent (default 1) in the case where the hypothesess
+  sum to a power of the goal. -/
+  power  : ℕ
+  deriving FromJson, Repr
+
 /-- The result of a sage call in the success case. -/
 structure SageSuccess where
   /-- The script returns a string containing python script to be sent to the remote server,
   when the tracing option is set. -/
   trace : Option String := none
   /-- The main result of the function call is an array of polynomials
-  parallel to the input list of hypotheses. -/
-  data : Option (Array Poly) := none
+  parallel to the input list of hypotheses and an exponent for the goal. -/
+  data : Option SageCoeffAndPower := none
   deriving FromJson, Repr
 
 /-- The result of a sage call in the failure case. -/
@@ -338,7 +353,7 @@ def polyrith (g : MVarId) (only : Bool) (hyps : Array Expr)
     match ← sageOutput (createSageArgs traceOnly α vars hyps' tgt) with
     | .ok { trace, data } =>
       if let some trace := trace then logInfo trace
-      if let some polys := data then
+      if let some {coeffs := polys, power := pow} := data then
         let vars ← liftM <| (← get).atoms.mapM delab
         let p ← Poly.sumM (polys.zip hyps') fun (p, src, _) => do
           let h := .hyp (← delab (match src with | .input i => hyps[i]! | .fvar h => .fvar h))
@@ -348,7 +363,8 @@ def polyrith (g : MVarId) (only : Bool) (hyps : Array Expr)
         let stx := (withRef (← getRef) <| p.toSyntax vars).run
         let tac ←
           if let .const 0 := p then `(tactic| linear_combination)
-          else `(tactic| linear_combination $stx:term)
+          else if pow = 1 then `(tactic| linear_combination $stx:term)
+          else `(tactic| linear_combination (exp := $(quote pow)) $stx:term)
         try
           guard (← Elab.runTactic g tac).1.isEmpty
         catch _ => throwError

scripts/polyrith_sage.py

@@ -19,19 +19,29 @@ def create_query(type: str, n_vars: int, eq_list, goal_type):
     """ Create a query to invoke Sage's `MPolynomial_libsingular.lift`. See
     https://github.com/sagemath/sage/blob/f8df80820dc7321dc9b18c9644c3b8315999670b/src/sage/rings/polynomial/multi_polynomial_libsingular.pyx#L4472-L4518
     for a description of this method. """
-    var_list = ", ".join([f"var{i}" for i in range(n_vars)])
+    var_list = [f"var{i}" for i in range(n_vars)] + ['aux']
     query = f'''
 if {n_vars!r} != 0:
-    P = PolynomialRing({type_str(type)}, 'var', {n_vars!r})
-    [{var_list}] = P.gens()
+    P = PolynomialRing({type_str(type)}, {var_list})
+    [{", ".join(var_list)}] = P.gens()
+    p = P({goal_type})
+    gens = {eq_list} + [1 - p*aux]
+    I = P.ideal(gens)
+    coeffs = P(1).lift(I)
+    power = max(cf.degree(aux) for cf in coeffs)
+    coeffs = [P(cf.subs(aux = 1/p)*p^power) for cf in coeffs[:int(-1)]]
+    print(str(power)+';'+serialize_polynomials(coeffs))
 else:
     # workaround for a Sage shortcoming with `n_vars = 0`,
     # `TypeError: no conversion of this ring to a Singular ring defined`
+    # In this case, there is no need to look for membership in the *radical*;
+    # we just check for membership in the ideal, and return exponent 1 
+    # if coefficients are found.
     P = PolynomialRing({type_str(type)}, 'var', 1)
-p = P({goal_type})
-I = P.ideal({eq_list})
-coeffs = p.lift(I)
-print(serialize_polynomials(coeffs))
+    p = P({goal_type})
+    I = P.ideal({eq_list})
+    coeffs = p.lift(I)
+    print('1;'+serialize_polynomials(coeffs))
 '''
     return query
 
@@ -42,12 +52,17 @@ def __init__(self, ename, evalue, message='Error in Sage communication'):
         self.message = message
         super().__init__(self.message)
 
+def parse_response(resp: str) -> str:
+    exp, data = resp.split(';', 1)
+    return dict(power=int(exp), coeffs=json.loads(data))
+
+
 def evaluate_in_sage(query: str) -> str:
     data = {'code': query}
     headers = {'content-type': 'application/x-www-form-urlencoded'}
     response = requests.post('https://sagecell.sagemath.org/service', data, headers=headers).json()
     if response['success']:
-        return json.loads(response.get('stdout'))
+        return parse_response(response.get('stdout'))
     elif 'execute_reply' in response and 'ename' in response['execute_reply'] and 'evalue' in response['execute_reply']:
         raise EvaluationError(response['execute_reply']['ename'], response['execute_reply']['evalue'])
     else:

test/polyrith.lean

@@ -561,6 +561,26 @@ by polyrith
 --   polyrith,
 -- end
 
+/-
+
+### Examples with exponent
+
+-/
+
+example (x y z : ℚ) (h : x = y) (h2 : x * y = 0) : x + y*z = 0 := by
+  polyrith
+
+example (K : Type)
+    [Field K]
+    [CharZero K]
+    {x y z : K}
+    (h₂ : y ^ 3 + x * (3 * z ^ 2) = 0)
+    (h₁ : x ^ 3 + z * (3 * y ^ 2) = 0)
+    (h₀ : y * (3 * x ^ 2) + z ^ 3 = 0)
+    (h : x ^ 3 * y + y ^ 3 * z + z ^ 3 * x = 0) :
+    x = 0 := by
+  polyrith
+
 /-!
 ### With trace enabled
 Here, the tactic will trace the command that gets sent to sage,