gh-37094: Fix bug with pari interface in `classical_modular_polynomial`
    
There is a bug in the new `classical_modular_polynomial` function
created in #36190.

A simple example to reproduce the bug is the following

```
p = 2^216 * 3^137 - 1
F.<i> = GF(p^2, modulus=[1,0,1])
E = EllipticCurve(F, [0, 6, 0, 1, 0])
classical_modular_polynomial(2, E.j_invariant())
```

This will fail with a `TypeError`.

The bug is due to the interface with the pari function `polmodular`. In
particular, contrary to the [documentation](http://pari.math.u-
bordeaux.fr/dochtml/ref/Polynomials_and_power_series.html#polmodular),
the `polmodular` function will only evaluate $\Phi_\ell(j, Y)$ for
$j$-invariants that are defined over $\mathbb F_p$, and not over any
generic finite field.
If however `j.parent()` is a generic finite field, but `j` itself is
defined over the prime field, pari will evaluate that and the current
implementation fails to convert back to a sage polynomial.

The proposed fix is to cast the result of the pari call to `ZZ['Y']`,
since the result of `polmodular` currently returns a polynomial in the
correct $\mathbb Z/n\mathbb Z[Y]$.

#sd123
    
URL: #37094
Reported by: Riccardo Zanotto
Reviewer(s): Frédéric Chapoton, Peter Bruin, Riccardo Zanotto

Author: Release Manager

src/sage/schemes/elliptic_curves/mod_poly.py

@@ -85,7 +85,24 @@ def classical_modular_polynomial(l, j=None):
         sage: j = GF(q).random_element()
         sage: l = random_prime(50)
         sage: Y = polygen(parent(j), 'Y')
-        sage: classical_modular_polynomial(l,j) == classical_modular_polynomial(l)(j,Y)
+        sage: classical_modular_polynomial(l, j) == classical_modular_polynomial(l)(j, Y)
+        True
+        sage: p = 2^216 * 3^137 - 1
+        sage: F.<i> = GF((p,2), modulus=[1,0,1])
+        sage: l = random_prime(50)
+        sage: j = F.random_element()
+        sage: Y = polygen(parent(j), 'Y')
+        sage: classical_modular_polynomial(l, j) == classical_modular_polynomial(l)(j, Y)
+        True
+        sage: E = EllipticCurve(F, [0, 6, 0, 1, 0])
+        sage: j = E.j_invariant()
+        sage: l = random_prime(50)
+        sage: classical_modular_polynomial(l, j) == classical_modular_polynomial(l)(j, Y)
+        True
+        sage: R.<Y> = QQ['Y']
+        sage: j = QQ(1/2)
+        sage: l = random_prime(50)
+        sage: classical_modular_polynomial(l, j) == classical_modular_polynomial(l)(j, Y)
         True
     """
     l = ZZ(l)
@@ -133,12 +150,14 @@ def classical_modular_polynomial(l, j=None):
     # Now try to get the instantiated modular polynomial directly from PARI.
     # This should be slightly more efficient (in particular regarding memory
     # usage) than computing and evaluating the generic modular polynomial.
+    # This currently only works if we are over Z/nZ.
     try:
         pari_Phi = pari.polmodular(l, 0, j)
+        return R(pari_Phi)
     except PariError:
         pass
-    else:
-        return R(pari_Phi)
+    except TypeError:
+        return R(ZZ['Y'](pari_Phi))
 
     # Nothing worked. Fall back to computing the generic modular polynomial
     # and simply evaluating it.