sage: Fix incompatibility with sage 9.4

[real-or-random]

allexprs is already the product all numerators. Don't take it's
numerator again.

Fixes #1067.

To be honest, I have no idea what I'm doing here.

[jonasnick]

I get the same errors here as with the change_ring "workaround":

Formula secp256k1_gej_add_var:
  add:
    branch 4: FAIL, <map object at 0x7f417d48b820> fails (assuming )

Formula secp256k1_gej_add_zinv_var:
  add:
    branch 4: FAIL, <map object at 0x7f417d67d420> fails (assuming )

[real-or-random]

@jonasnick Oh I missed this... I believe this a real of case of incompleteness (something in sage changed?) and now our prover isn't complete enough to prove the formula. But the fix is easy. If you prints() to the relevant part of prove_nonzero like this...

  print("nonzero", nonzero)
  for (f, n) in allexprs.factor():
    print("factor", f)

... you get

nonzero {by, ax, bx, Az, bx - ax, ay, Bz}
factor -bx + ax

So it put bx - ax in the set of nonzero expressions but then fails to see that then -bx + ax is nonzero as well. I added a commit to address this.

Disclaimer: I don't really know what this code is doing in the end, so please take this with a grain of salt. @sipa should check if I'm right.

[sipa]

I wonder if this is generic enough, and there couldn't be arbitrary constant factor differences. If the results aren't normalized in a way that result in unique signs, other factors might differ too. What about:

  allexprs = reduce(lambda a,b: numerator(a)*numerator(b), exprs, 1)
  for (f, n) in allexprs.factor():
    if not any((f % nz) == 0 and (nz % f) == 0 for nz in nonzero):
      ok = False

[real-or-random]

Both the elements of nonzero and of allexprs.factor() are prime factors obtained by factorization, so I think they can differ at most by their unit() (https://doc.sagemath.org/html/en/reference/structure/sage/structure/factorization.html#sage.structure.factorization.Factorization.unit), which is -1 or 1.

So an element of nonzero or allexprs.factor() cannot have an additional factor c with c != -1 and c != 1 because then it wouldn't be a prime factor. But there could be a constant prime element c in allexprs.factor(), which is then trivially nonzero but the prover fails this.

So if we want to allow for this, I'd suggest (which does not on the % operator for which I can't find docs on):

  for (f, n) in allexprs.factor():
    if f not in nonzero and -f not in nonzero and not f.is_constant():
      ok = False

Docs are here:
https://doc.sagemath.org/html/en/reference/polynomial_rings/sage/rings/polynomial/multi_polynomial_element.html#sage.rings.polynomial.multi_polynomial_element.MPolynomial_polydict.is_constant

[sipa]

Are you sure? These are (multivariate) polynomial rings, and all constant terms are units for that.

Specifically, I'm concerned about the case where you'd have factors a + 2*b or b + 1/2*a, depending on which variable is considered the one to make monic.

[real-or-random]

Are you sure? These are (multivariate) polynomial rings, and all constant terms are units for that.

I'm not sure at all.

Specifically, I'm concerned about the case where you'd have factors a + 2b or b + 1/2a, depending on which variable is considered the one to make monic.

Hm yeah, I'd imagine this is the same everytime in the same ring class in sage but I don't know.

Would you then change every similar if within this function?

[sipa]

Hm yeah, I'd imagine this is the same everytime in the same ring class in sage but I don't know.

Indeed, but the same reasoning should IMO apply to the just the factor 1 or -1. If there is a canonical ordering of factors, and e.g. factor() always returns the lexicographically smallest list of factors, there should never be any ambiguity at all. The fact that the sign differs makes me suspicious, so I'd rather be overly cautious (the code still works fine with that change).

EDIT: I realize the above is incorrect. Just because individual inputs to factor() are canonical doesn't imply the same is true for all factors individually. Still, I'm unsure if 1 and -1 are the only possible differences.

Would you then change every similar if within this function?

Yeah. Perhaps introduce a unit_ratio_between(a, b) function or so to abstract it and make it a bit more readable?

[jonasnick]

I did some independent experiments and the only factorizations I was able to get Sage (9.4) to produce had the coefficients of the non-unit factors be integers with gcd 1. This would make them unique if it weren't for this +,- inconsistency.

I considered that this is a regression in Sage >9.2 with order="invlex", because replacing it with order="lex" makes the script work without the -f not in nonzero fix. But no.

sage: a, b, c = QQ["a,b,c"].gens()
....: f = 1 / 2 * a - b
....: f.factor()
....: 
(1/2) * (a - 2*b)
sage: f2 = (1 / 2 * a - b) * (a - c)
....: f2.factor()
(1/2) * (-a + c) * (-a + 2*b)

I'd be fine with leaving the PR as is. Another solution would be to wrap the factor() method such that it behaves as we would expect, i.e. the smallest variable name in a factor gets a positive coefficient (and our other assumptions are asserted).

[jonasnick]

nit: perhaps note in the commit description that this commit is required to make the script pass with sage 9.4.

[sipa]

$ sage
┌────────────────────────────────────────────────────────────────────┐
│ SageMath version 9.2, Release Date: 2020-10-24                     │
│ Using Python 3.9.7. Type "help()" for help.                        │
└────────────────────────────────────────────────────────────────────┘
sage: R.<ax,bx,ay,by,Az,Bz,Ai,Bi> = PolynomialRing(QQ,8,order='invlex')
sage: R
Multivariate Polynomial Ring in ax, bx, ay, by, Az, Bz, Ai, Bi over Rational Field
sage: R(10 * ax + bx / 2).factor()
(1/2) * (bx + 20*ax)
sage: R(20 * ax + bx).factor()
bx + 20*ax
sage: R(20 * ax + bx).factor().unit()
1
sage: R(10 * ax + bx/2).factor().unit()
1/2
sage: list(R(20 * ax + bx).factor())
[(bx + 20*ax, 1)]
sage: list(R(10 * ax + bx / 2).factor())
[(bx + 20*ax, 1)]
sage: R((-2 * (bx - ax)) ^ 1).factor()
(2) * (-bx + ax)
sage: R((-2 * (bx - ax)) ^ 2).factor()
(4) * (-bx + ax)^2
sage: R((-2 * (bx - ax)) ^ 3).factor()
(8) * (-bx + ax)^3
sage: R((1/10 * ax + 1/15*bx)*(ax/7+bx)).factor()
(1/210) * (2*bx + 3*ax) * (7*bx + ax)

[sipa]

Is this (f not in nonzero) and (-f not in nonzero) needed in other places where this pattern occurs below?

[real-or-random]

Another solution would be to wrap the factor() method such that it behaves as we would expect, i.e. the smallest variable name in a factor gets a positive coefficient (and our other assumptions are asserted).

That's a clean idea. I did this.

Is this (f not in nonzero) and (-f not in nonzero) needed in other places where this pattern occurs below?

Not for the current group laws but I added the wrapper to the other places where we call factor().

[real-or-random]

I added sage to CI. Ready from my side (if CI passes...)

[real-or-random]

Upstream report: https://trac.sagemath.org/ticket/33283

[jonasnick]

ACK bed969b

[real-or-random]

This is buggy for --exhaustive, let me fix this.

[real-or-random]

Ok done. Also cherry-picked a commit from https://github.com/real-or-random/secp256k1/commits/202202-gej_add_ge that I discovered when working on #1033.

[jonasnick]

The commit looks fine but is there a specific reason why it was added?

[real-or-random]

I added it because it fixes a bug in the group prover that prevents us to look into other group formulas. In particular, it stopped me from playing around with the formula here because it may return a literal 1: https://github.com/real-or-random/secp256k1/blob/dc53e24a1b63ec07a9ab086280cab1355cd3a881/sage/prove_group_implementations.sage#L209

[jonasnick]

Ah, thanks that makes sense.

[sipa]

ACK ebb1bee

[jonasnick]

ACK ebb1bee