gh-39028: Numerical evaluation of modular form & implement conversion of more modular form spaces to Pari
    
As in the title.

Currently modular form can be evaluated at symbolic `q` that is a
generator of a power series, I think we can also allow evaluating it at
a numerical `q`.

### 📝 Checklist

<!-- Put an `x` in all the boxes that apply. -->

- [x] The title is concise and informative.
- [x] The description explains in detail what this PR is about.
- [x] I have linked a relevant issue or discussion. (not aware of one)
- [x] I have created tests covering the changes.
- [ ] I have updated the documentation and checked the documentation
preview.

### ⌛ Dependencies

<!-- List all open PRs that this PR logically depends on. For example,
-->
<!-- - #12345: short description why this is a dependency -->
<!-- - #34567: ... -->
    
URL: #39028
Reported by: user202729
Reviewer(s): Travis Scrimshaw, user202729

Author: Release Manager

src/sage/modular/modform/ambient_eps.py

@@ -288,3 +288,18 @@ def hecke_module_of_level(self, N):
             return constructor.ModularForms(self.character().restrict(N), self.weight(), self.base_ring(), prec=self.prec())
         else:
             raise ValueError("N (=%s) must be a divisor or a multiple of the level of self (=%s)" % (N, self.level()))
+
+    def _pari_init_(self):
+        """
+        Conversion to Pari.
+
+        EXAMPLES::
+
+            sage: m = ModularForms(DirichletGroup(17).0^2, 2)
+            sage: pari.mfdim(m)
+            3
+            sage: pari.mfparams(m)
+            [17, 2, Mod(9, 17), 4, t^4 + 1]
+        """
+        from sage.libs.pari import pari
+        return pari.mfinit([self.level(), self.weight(), self.character()], 4)

src/sage/modular/modform/cuspidal_submodule.py

@@ -377,6 +377,21 @@ def _compute_q_expansion_basis(self, prec=None):
         return [weight1.modular_ratio_to_prec(chi, f, prec) for f in
             weight1.hecke_stable_subspace(chi)]
 
+    def _pari_init_(self):
+        """
+        Conversion to Pari.
+
+        EXAMPLES::
+
+            sage: A = CuspForms(DirichletGroup(23, QQ).0, 1)
+            sage: pari.mfparams(A)
+            [23, 1, -23, 1, t + 1]
+            sage: pari.mfdim(A)
+            1
+        """
+        from sage.libs.pari import pari
+        return pari.mfinit([self.level(), self.weight(), self.character()], 1)
+
 
 class CuspidalSubmodule_wt1_gH(CuspidalSubmodule):
     r"""

src/sage/modular/modform/eisenstein_submodule.py

@@ -586,6 +586,22 @@ class EisensteinSubmodule_eps(EisensteinSubmodule_params):
         q^5 + (zeta3 + 1)*q^8 + O(q^10)
         ]
     """
+    def _pari_init_(self):
+        """
+        Conversion to Pari.
+
+        EXAMPLES::
+
+            sage: e = DirichletGroup(27,CyclotomicField(3)).0**2
+            sage: M = ModularForms(e,2,prec=10).eisenstein_subspace()
+            sage: pari.mfdim(M)
+            6
+            sage: pari.mfparams(M)
+            [27, 2, Mod(10, 27), 3, t^2 + t + 1]
+        """
+        from sage.libs.pari import pari
+        return pari.mfinit([self.level(), self.weight(), self.character()], 3)
+
     # TODO
     # def _compute_q_expansion_basis(self, prec):
     #     B = EisensteinSubmodule_params._compute_q_expansion_basis(self, prec)

src/sage/modular/modform/element.py

@@ -219,6 +219,38 @@ def _repr_(self):
         """
         return str(self.q_expansion())
 
+    def _pari_init_(self):
+        """
+        Conversion to Pari.
+
+        TESTS::
+
+            sage: M = EisensteinForms(96, 2)
+            sage: M.6
+            O(q^6)
+            sage: M.7
+            O(q^6)
+            sage: pari(M.6) == pari(M.7)
+            False
+            sage: pari(M.6).mfcoefs(10)
+            [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0]
+
+            sage: M = ModularForms(DirichletGroup(17).0^2, 2)
+            sage: pari(M.0).mfcoefs(5)
+            [0, 1, Mod(-t^3 + t^2 - 1, t^4 + 1), Mod(t^3 - t^2 - t - 1, t^4 + 1), Mod(2*t^3 - t^2 + 2*t, t^4 + 1), Mod(-t^3 - t^2, t^4 + 1)]
+            sage: M.0.qexp(5)
+            q + (-zeta8^3 + zeta8^2 - 1)*q^2 + (zeta8^3 - zeta8^2 - zeta8 - 1)*q^3 + (2*zeta8^3 - zeta8^2 + 2*zeta8)*q^4 + O(q^5)
+        """
+        from sage.libs.pari import pari
+        from sage.rings.number_field.number_field_element import NumberFieldElement
+        M = pari(self.parent())
+        f = self.qexp(self.parent().sturm_bound())
+        coefficients = [
+            x.__pari__('t') if isinstance(x, NumberFieldElement) else x
+            for x in f]
+        # we cannot compute pari(f) directly because we need to set the variable name as t
+        return M.mflinear(M.mftobasis(coefficients + [0] * (f.prec() - len(coefficients))))
+
     def __call__(self, x, prec=None):
         """
         Evaluate the `q`-expansion of this modular form at x.
@@ -233,9 +265,150 @@ def __call__(self, x, prec=None):
 
             sage: f(0)
             0
-        """
+
+        Evaluate numerically::
+
+            sage: f = ModularForms(1, 12).0
+            sage: f(0.3)  # rel tol 1e-12
+            2.34524576548591e-6
+            sage: f = EisensteinForms(1, 4).0
+            sage: f(0.9)  # rel tol 1e-12
+            1.26475942209241e7
+
+        TESTS::
+
+            sage: f = ModularForms(96, 2).0
+            sage: f(0.3)  # rel tol 1e-12
+            0.299999997396191
+            sage: f(0.0+0.0*I)
+            0
+
+        For simplicity, ``float`` or ``complex`` input are converted to ``CC``, except for
+        input ``0`` where exact result is returned::
+
+            sage: result = f(0.3r); result  # rel tol 1e-12
+            0.299999997396191
+            sage: result.parent()
+            Complex Field with 53 bits of precision
+            sage: result = f(0.3r + 0.3jr); result  # rel tol 1e-12
+            0.299999359878484 + 0.299999359878484*I
+            sage: result.parent()
+            Complex Field with 53 bits of precision
+
+        Symbolic numerical values use precision of ``CC`` by default::
+
+            sage: f(sqrt(1/2))  # rel tol 1e-12
+            0.700041406692037
+            sage: f(sqrt(1/2)*QQbar.zeta(8))  # rel tol 1e-12
+            0.496956554651376 + 0.496956554651376*I
+
+        Higher precision::
+
+            sage: f(ComplexField(128)(0.3))  # rel tol 1e-36
+            0.29999999739619131029285166058750164058
+            sage: f(ComplexField(128)(1+2*I)/3)  # rel tol 1e-36
+            0.32165384572356882556790532669389900691 + 0.67061244638367586302820790711257777390*I
+
+        Confirm numerical evaluation matches the q-expansion::
+
+            sage: f = EisensteinForms(1, 4).0
+            sage: f(0.3)  # rel tol 1e-12
+            741.741819297986
+            sage: f.qexp(50).polynomial()(0.3)  # rel tol 1e-12
+            741.741819297986
+
+        With a nontrivial character::
+
+            sage: M = ModularForms(DirichletGroup(17).0^2, 2)
+            sage: M.0(0.5)  # rel tol 1e-12
+            0.166916655031616 + 0.0111529051752428*I
+            sage: M.0.qexp(60).polynomial()(0.5)  # rel tol 1e-12
+            0.166916655031616 + 0.0111529051752428*I
+
+        Higher precision::
+
+            sage: f(ComplexField(128)(1+2*I)/3)  # rel tol 1e-36
+            429.19994832206294278688085399056359632 - 786.15736284188243351153830824852974995*I
+            sage: f.qexp(400).polynomial()(ComplexField(128)(1+2*I)/3)  # rel tol 1e-36
+            429.19994832206294278688085399056359631 - 786.15736284188243351153830824852974999*I
+
+        Check ``SR`` does not make the result lose precision::
+
+            sage: f(ComplexField(128)(1+2*I)/3 + x - x)  # rel tol 1e-36
+            429.19994832206294278688085399056359632 - 786.15736284188243351153830824852974995*I
+        """
+        from sage.rings.integer import Integer
+        from sage.misc.functional import log
+        from sage.structure.element import parent
+        from sage.rings.complex_mpfr import ComplexNumber
+        from sage.rings.cc import CC
+        from sage.rings.real_mpfr import RealNumber
+        from sage.symbolic.constants import pi
+        from sage.rings.imaginary_unit import I  # import from here instead of sage.symbolic.constants to avoid cast to SR
+        from sage.symbolic.expression import Expression
+        if isinstance(x, Expression):
+            try:
+                x = x.pyobject()
+            except TypeError:
+                pass
+        if x in CC:
+            if x == 0:
+                return self.qexp(1)[0]
+            if not isinstance(x, (RealNumber, ComplexNumber)):
+                x = CC(x)  # might lose precision if this is done unconditionally (TODO what about interval and ball types?)
+        if isinstance(x, (RealNumber, ComplexNumber)):
+            return self.eval_at_tau(log(x)/(2*parent(x)(pi)*I))  # cast to parent(x) to force numerical evaluation of pi
         return self.q_expansion(prec)(x)
 
+    def eval_at_tau(self, tau):
+        r"""
+        Evaluate this modular form at the half-period ratio `\tau`.
+        This is related to `q` by `q = e^{2\pi i \tau}`.
+
+        EXAMPLES::
+
+            sage: f = ModularForms(1, 12).0
+            sage: f.eval_at_tau(0.3 * I)  # rel tol 1e-12
+            0.00150904633897550
+
+        TESTS:
+
+        Symbolic numerical values use precision of ``CC`` by default::
+
+            sage: f.eval_at_tau(sqrt(1/5)*I)  # rel tol 1e-12
+            0.0123633234207127
+            sage: f.eval_at_tau(sqrt(1/2)*QQbar.zeta(8))  # rel tol 1e-12
+            -0.114263670441098
+
+        For simplicity, ``complex`` input are converted to ``CC``::
+
+            sage: result = f.eval_at_tau(0.3jr); result  # rel tol 1e-12
+            0.00150904633897550
+            sage: result.parent()
+            Complex Field with 53 bits of precision
+
+        Check ``SR`` does not make the result lose precision::
+
+            sage: f = EisensteinForms(1, 4).0
+            sage: f.eval_at_tau(ComplexField(128)(1+2*I)/3 + x - x)  # rel tol 1e-36
+            -1.0451570582202060056197878314286036966 + 2.7225112098519803098203933583286590274*I
+        """
+        from sage.libs.pari.convert_sage import gen_to_sage
+        from sage.libs.pari import pari
+        from sage.rings.cc import CC
+        from sage.rings.complex_mpfr import ComplexNumber, ComplexField
+        from sage.rings.real_mpfr import RealNumber
+        from sage.symbolic.expression import Expression
+        if isinstance(tau, Expression):
+            try:
+                tau = tau.pyobject()
+            except TypeError:
+                pass
+        if not isinstance(tau, (RealNumber, ComplexNumber)):
+            tau = CC(tau)
+        precision = tau.prec()
+        return ComplexField(precision)(pari.mfeval(self.parent(), self, tau, precision=precision))
+
     @cached_method
     def valuation(self):
         """