Opened 3 years ago

Closed 3 years ago

# Optimize Psi2()

Reported by: Owned by: embray minor sage-7.4 elliptic curves Erik Bray, Jeroen Demeyer Jeroen Demeyer, Erik Bray N/A bdce5dc (Commits) bdce5dced1e9f811616836253c4e7cd45b71fcf3

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.

### comment:1 Changed 3 years ago by embray

• Status changed from new to needs_review

### comment:2 Changed 3 years ago by jdemeyer

• Description modified (diff)

### comment:3 Changed 3 years ago by jdemeyer

• Description modified (diff)

### comment:4 Changed 3 years ago by jdemeyer

• Branch changed from u/embray/psi2 to u/jdemeyer/psi2

### comment:5 Changed 3 years ago by jdemeyer

• Commit changed from 9b22d75ac27fcae0d07065daccab96ba14562f79 to 10ff2ec9ca61fb39187b9ab2fe0fb1965251ca24

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`

### comment:6 Changed 3 years ago by embray

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

### comment:7 Changed 3 years ago by git

• Commit changed from 10ff2ec9ca61fb39187b9ab2fe0fb1965251ca24 to bdce5dced1e9f811616836253c4e7cd45b71fcf3

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

 ​bdce5dc `Optimize Psi2()`

### comment:8 Changed 3 years ago by jdemeyer

• Authors changed from Erik Bray to Erik Bray, Jeroen Demeyer
• Description modified (diff)
• Summary changed from Make Psi2 function less likely to break the stack to Optimize Psi2()

### comment:9 Changed 3 years ago by jdemeyer

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

### comment:10 follow-up: ↓ 11 Changed 3 years ago by 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 :)

Version 0, edited 3 years ago by embray (next)

### comment:11 in reply to: ↑ 10 Changed 3 years ago by jdemeyer

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 :)

1. 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)`.
1. I remove `raise ValueError("%s must be one of %s."%(l,hyperelliptic_primes))` since that exception is already raised by `_hyperelliptic_isogeny_data`.
1. I rename `y` -> `v` to avoid confusion (the first `v` and the second `v` are really the same thing conceptually).
1. I remove the variable `x` from the first `R` since it's never used. This is probably where the main speed-up comes from.
1. I replace `s = [1]` by `s = [K(1)]` to ensure that all elements of `s` have `K` as parent.
1. 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`.
1. 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`).

### comment:12 Changed 3 years ago by embray

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

### comment:13 Changed 3 years ago by embray

• Status changed from needs_review to positive_review

### comment:14 Changed 3 years ago by vbraun

• Status changed from positive_review to needs_work

Reviewer name..

### comment:15 Changed 3 years ago by jdemeyer

• Reviewers set to Jeroen Demeyer, Erik Bray
• Status changed from needs_work to positive_review

### comment:16 Changed 3 years ago by vbraun

• Branch changed from u/jdemeyer/psi2 to bdce5dced1e9f811616836253c4e7cd45b71fcf3
• Resolution set to fixed
• Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.