Optimize Psi2()

[embray]

We greatly optimize the Psi2() function.

Before:

sage: from sage.schemes.elliptic_curves.isogeny_small_degree import Psi2
sage: timeit('Psi2.f(71)')
5 loops, best of 3: 9.36 s per loop

After:

sage: from sage.schemes.elliptic_curves.isogeny_small_degree import Psi2
sage: timeit('Psi2.f(71)')
5 loops, best of 3: 875 ms per loop

(note the use of Psi2.f() to bypass caching)

---

This also provides a workaround to a test that caused the Python interpreter to segfault on Windows. The specific call that caused the segfault is:

Psi2(71)

But in brief, the crash occurs because we have (in the original code) a large element of Univariate Quotient Polynomial Ring in v over Multivariate Polynomial Ring in x, u over Rational Field with modulus v^2 - u^4 + 10*u^3 + 3*u^2 - 4*u + 8. This is then being cast simply to a multivariate rational polynomial over x, u, v. Because there is no direct coercion between these types this involves eval()-ing the polynomial as a Python expression.

The problem is that this expression can become too large for the stack--specifically in Python's symbol visibility analysis, a step it performs just before compiling an expression to bytecode. The implementation of that recurses into binary expressions, leading to a stack overflow for such a large expression. This issue has come up once before in #16014 where it was worked around by a rewrite of the code. This particular test worked on Linux where the default stack size is 8MB, but it crashed on Windows where the typical stack is just 1MB.

The optimization gives smaller expressions to eval.

A more general fix to this problem would be nice though--I'm writing up some thoughts I have on it in a separate post.

Component: elliptic curves

Author: Erik Bray, Jeroen Demeyer

Branch/Commit: bdce5dc

Reviewer: Jeroen Demeyer, Erik Bray

Issue created by migration from https://trac.sagemath.org/ticket/21388

[jdemeyer]

Changed branch from u/embray/psi2 to u/jdemeyer/psi2

[jdemeyer]

Changed commit from 9b22d75 to 10ff2ec

[jdemeyer]

comment:5

In fact, we don't need the variable x in the ring R at all...

Hang on...

---

New commits:

10ff2ec

Workaround to prevent this calculation to blow the stack

[embray]

comment:6

Is the idea that you would just take the multiplication by x outside?

[sagetrac-git]

Changed commit from 10ff2ec to bdce5dc

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

bdce5dc

Optimize Psi2()

[jdemeyer]

Changed author from Erik Bray to Erik Bray, Jeroen Demeyer

[jdemeyer]

comment:9

This also happens to speed up the function by a factor of 10.

Please review...

[embray]

comment:10

Some of this makes sense to me (I was already irked by the use of += []) but the rest is a little beyond my understanding of how Sage works (though I'd appreciate a brief explanation :)

That said it's no surprise to me that more could be done with this than what I did.

[jdemeyer]

comment:11

Replying to @embray:

Some of this makes sense to me (I was already irked by the use of += []) but the rest is a little beyond my understanding of how Sage works (though I'd appreciate a brief explanation :)

0. Some general code cleanup and reformatting, in particular adding some spaces.

1. To construct the rational number 1/x, I use QQ((1,x)) instead of 1/QQ(x).

2. I remove raise ValueError("%s must be one of %s."%(l,hyperelliptic_primes)) since that exception is already raised by _hyperelliptic_isogeny_data.

3. I rename y -> v to avoid confusion (the first v and the second v are really the same thing conceptually).

4. I remove the variable x from the first R since it's never used. This is probably where the main speed-up comes from.

5. I replace s = [1] by s = [K(1)] to ensure that all elements of s have K as parent.

6. I replace your R((-1)**i*s[i]*x**(d-i)) by (-1)**i * x**(d-i) * R(s[i]). Since s[i] does not involve the variable x at all, I can indeed first convert s[i] to R and then multiply with x.

7. I replace s[i] by s[i].lift() before conversion. This means "forgetting" the quotient and converting s[i] to a true polynomial. Put differently, it's the conversion from K to L. Then there is a coercion from L to R (but not from K to R).

[embray]

comment:12

Thanks, that all makes sense. In particular I didn't realize what lift() did.

[vbraun]

comment:14

Reviewer name..

[jdemeyer]

Reviewer: Jeroen Demeyer, Erik Bray

[vbraun]

Changed branch from u/jdemeyer/psi2 to bdce5dc