Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[Merged by Bors] - feat: improve polyrith by testing for membership in the radical #7790

Closed
wants to merge 15 commits into from
Closed
24 changes: 20 additions & 4 deletions Mathlib/Tactic/Polyrith.lean
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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 : ℕ
robertylewis marked this conversation as resolved.
Show resolved Hide resolved
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
robertylewis marked this conversation as resolved.
Show resolved Hide resolved
deriving FromJson, Repr

/-- The result of a sage call in the failure case. -/
Expand Down Expand Up @@ -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))
Expand All @@ -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
Expand Down
31 changes: 23 additions & 8 deletions scripts/polyrith_sage.py
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Comment on lines +26 to +33
Copy link
Member

@alexjbest alexjbest Nov 3, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't really understand why this is done with this trick of checking if 1 is in the ideal <hyps, 1/goal>, I would expect singular/sage has a way to ask for a lift to the radical ideal directly which I would hope would be easier to maintain in the long run (especially if we have aspirations of making this work for non-char 0 rings properly). Did you consider this approach?
I'm not completely against merging this, but the whole thing is starting to look a bit scary ;)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

After experimenting a bit I believe that using the builtin radical support, while probably a better algorithm underneath, is not very convenient to use with the way things are set up here, so I'm happy with this I guess

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, I remember looking for this when I first implemented this (granted, a year and a half ago) and not turning up anything useful. Sage/Singular will check for membership in the radical just fine but extracting the power and coefficient is a pain. My impression is that what I do in this PR is a "standard" approach and that a native Sage implementation would do roughly the same, although I'm somewhat out of my element here.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When I read the sage docstring for radical earlier it said:

             From the Singular manual: A combination of the algorithms
            of Krick/Logar and Kemper is used. Works also in positive
            characteristic (Kempers algorithm).

but I have no idea how accurate that is!

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Right -- but there's space between the membership check $p \in \sqrt{\langle H_1, \ldots, H_n \rangle}$, and explicitly producing $k$ and $Q_i$ such that $p^k = \sum Q_iH_i$. Sage's radical makes the former easy, but I don't see how to get from there to the latter. Maybe there's a way, I'd love to hear it.

(It sounds like you came to a similar conclusion, just leaving this comment for reference later!)

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))
robertylewis marked this conversation as resolved.
Show resolved Hide resolved
'''
return query

Expand All @@ -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:
Expand Down
20 changes: 20 additions & 0 deletions test/polyrith.lean
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
Loading