gh-35260: Implement the logarithm and the exponential of a Drinfeld module
    
<!-- ^^^^^
Please provide a concise, informative and self-explanatory title.
Don't put issue numbers in there, do this in the PR body below.
For example, instead of "Fixes #1234" use "Introduce new method to
calculate 1+1"
-->
### 📚 Description

The goal of this PR is to implement the logarithm and the exponential of
a Drinfeld module.

### Background material

Let $\phi:T \mapsto \gamma(T) + g_1 \tau + \cdots + g_r \tau^r$ with
$g_r\neq 0$ be a rank $r$ Drinfeld module and let $\Lambda\subset
\mathbb{C}\_{\infty}$ be the associated rank $r$ lattice. Then recall
that the exponential $e\_{\phi}:\mathbb{C}\_{\infty} \rightarrow
\mathbb{C}\_{\infty}$ of the Drinfeld module $\phi$ is defined by:

$$
e_{\phi}(z) := z\prod_{\lambda\in \Lambda\setminus \{0\}} \left( 1 -
z/\lambda \right).
$$

It can be viewed as a power series in $z$:

$$
e_{\phi}(z) = z + \alpha_1 z^{q} + \alpha_2 z^{q^2} + \alpha_3 z^{q^3} +
\cdots
$$

Moreover, it satisfies the following functional equation:

$$
e_{\phi}(az) = \phi_a(e_{\phi}(z)).
$$

The logarithm of $\phi$ is defined to be the compositional inverse of
$e_{\phi}$, denoted $\mathrm{log}\_{\phi}$. By applying
$\mathrm{log}\_{\phi}$ on both side of the functional equation, we
obtain the following recurrence relation:

$$
\beta_k = \frac{1}{\gamma(T) - \gamma(T)^{q^i - 1}} \sum_{i = 0}^{k-1}
\beta_i g_{k-i}^{q^i}
$$

where $\beta_k$ is the $q^k$-th coefficient of $\mathrm{log}_{\phi}$ and
$\beta_0 = 1$.

### Examples

```
sage: A = GF(2)['T']
sage: K.<T> = Frac(A)
sage: phi = DrinfeldModule(A, [T, 1])
sage: q = A.base_ring().cardinality()
sage: log = phi.logarithm(); log
z + ((1/(T^2+T))*z^2) + ((1/(T^6+T^5+T^3+T^2))*z^4) + O(z^7)
sage: log[q] == -1/((T**q - T))
True
sage: log[q**2] == 1/((T**q - T)*(T**(q**2) - T))
True
sage: log[q**3] == -1/((T**q - T)*(T**(q**2) - T)*(T**(q**3) - T))
True
```

### Ideas

The idea here is to first compute the logarithm of a Drinfeld module
using lazy power series and then compute its compositional inverse to
find its exponential.


<!-- Describe your changes here in detail -->
<!-- Why is this change required? What problem does it solve? -->
<!-- If it resolves an open issue, please link to the issue here. For
example "Closes #1337" -->

### 📝 Checklist

<!-- Put an `x` in all the boxes that apply. -->
<!-- If your change requires a documentation PR, please link it
appropriately -->
<!-- If you're unsure about any of these, don't hesitate to ask. We're
here to help! -->

- [x] I have made sure that the title is self-explanatory and the
description concisely explains the PR.
- [ ] I have linked an issue or discussion.
- [x] I have created tests covering the changes.
- [x] I have updated the documentation accordingly.

### ⌛ Dependencies
<!-- List all open pull requests that this PR logically depends on -->
<!--
- #xyz: short description why this is a dependency
- #abc: ...
-->

Depends on #35026 because this PR implements Drinfeld modules.

CC: @xcaruso @ymusleh @kryzar
    
URL: #35260
Reported by: David Ayotte
Reviewer(s): Antoine Leudière, David Ayotte, Xavier Caruso

Author: Release Manager

src/sage/rings/function_field/drinfeld_modules/drinfeld_module.py

@@ -26,10 +26,13 @@
 
 from sage.categories.drinfeld_modules import DrinfeldModules
 from sage.categories.homset import Hom
+from sage.misc.cachefunc import cached_method
 from sage.misc.latex import latex
 from sage.misc.latex import latex_variable_name
+from sage.misc.lazy_import import lazy_import
 from sage.misc.lazy_string import _LazyString
 from sage.rings.integer import Integer
+from sage.rings.integer_ring import ZZ
 from sage.rings.polynomial.ore_polynomial_element import OrePolynomial
 from sage.rings.polynomial.polynomial_ring import PolynomialRing_general
 from sage.rings.ring_extension import RingExtension_generic
@@ -38,6 +41,8 @@
 from sage.structure.sequence import Sequence
 from sage.structure.unique_representation import UniqueRepresentation
 
+lazy_import('sage.rings.lazy_series_ring', 'LazyPowerSeriesRing')
+
 
 class DrinfeldModule(Parent, UniqueRepresentation):
     r"""
@@ -958,6 +963,134 @@ def coefficients(self, sparse=True):
         """
         return self._gen.coefficients(sparse=sparse)
 
+    @cached_method
+    def _compute_coefficient_exp(self, k):
+        r"""
+        Return the `q^k`-th coefficient of the exponential of this Drinfeld module.
+
+        INPUT:
+
+        - ``k`` (integer) -- the index of the coefficient
+
+        TESTS::
+
+            sage: A = GF(2)['T']
+            sage: K.<T> = Frac(A)
+            sage: phi = DrinfeldModule(A, [T, 1])
+            sage: q = A.base_ring().cardinality()
+            sage: phi._compute_coefficient_exp(0)
+            1
+            sage: phi._compute_coefficient_exp(1)
+            1/(T^2 + T)
+            sage: phi._compute_coefficient_exp(2)
+            1/(T^8 + T^6 + T^5 + T^3)
+            sage: phi._compute_coefficient_exp(3)
+            1/(T^24 + T^20 + T^18 + T^17 + T^14 + T^13 + T^11 + T^7)
+        """
+        k = ZZ(k)
+        if k.is_zero():
+            return self._base.one()
+        q = self._Fq.cardinality()
+        c = self._base.zero()
+        for i in range(k):
+            j = k - i
+            c += self._compute_coefficient_exp(i)*self._compute_coefficient_log(j)**(q**i)
+        return -c
+
+    def exponential(self, name='z'):
+        r"""
+        Return the exponential of this Drinfeld module.
+
+        Note that the exponential is only defined when the
+        `\mathbb{F}_q[T]`-characteristic is zero.
+
+        INPUT:
+
+        - ``name`` (string, default: ``'z'``) -- the name of the
+          generator of the lazy power series ring.
+
+        OUTPUT:
+
+        A lazy power series over the base field.
+
+        EXAMPLES::
+
+            sage: A = GF(2)['T']
+            sage: K.<T> = Frac(A)
+            sage: phi = DrinfeldModule(A, [T, 1])
+            sage: q = A.base_ring().cardinality()
+            sage: exp = phi.exponential(); exp
+            z + ((1/(T^2+T))*z^2) + ((1/(T^8+T^6+T^5+T^3))*z^4) + O(z^8)
+
+        The exponential is returned as a lazy power series, meaning that
+        any of its coefficients can be computed on demands::
+
+            sage: exp[2^4]
+            1/(T^64 + T^56 + T^52 + ... + T^27 + T^23 + T^15)
+            sage: exp[2^5]
+            1/(T^160 + T^144 + T^136 + ... + T^55 + T^47 + T^31)
+
+        Example in higher rank::
+
+            sage: A = GF(5)['T']
+            sage: K.<T> = Frac(A)
+            sage: phi = DrinfeldModule(A, [T, T^2, T + T^2 + T^4, 1])
+            sage: exp = phi.exponential(); exp
+            z + ((T/(T^4+4))*z^5) + O(z^8)
+
+        The exponential is the compositional inverse of the logarithm
+        (see :meth:`logarithm`)::
+
+            sage: log = phi.logarithm(); log
+            z + ((4*T/(T^4+4))*z^5) + O(z^8)
+            sage: exp.compose(log)
+            z + O(z^8)
+            sage: log.compose(exp)
+            z + O(z^8)
+
+        ::
+
+            sage: Fq.<w> = GF(3)
+            sage: A = Fq['T']
+            sage: phi = DrinfeldModule(A, [w, 1])
+            sage: phi.exponential()
+            Traceback (most recent call last):
+            ...
+            ValueError: characteristic must be zero (=T + 2)
+
+        TESTS::
+
+            sage: A = GF(2)['T']
+            sage: K.<T> = Frac(A)
+            sage: phi = DrinfeldModule(A, [T, 1])
+            sage: exp = phi.exponential()
+            sage: exp[2] == 1/(T**q - T)  # expected value
+            True
+            sage: exp[2^2] == 1/((T**(q**2) - T)*(T**q - T)**q)  # expected value
+            True
+            sage: exp[2^3] == 1/((T**(q**3) - T)*(T**(q**2) - T)**q*(T**q - T)**(q**2))  # expected value
+            True
+
+        REFERENCE:
+
+        See section 4.6 of [Gos1998]_ for the definition of the
+        exponential.
+        """
+        if self.category()._characteristic:
+            raise ValueError(f"characteristic must be zero (={self.characteristic()})")
+        L = LazyPowerSeriesRing(self._base, name)
+        zero = self._base.zero()
+        q = self._Fq.cardinality()
+
+        def coeff_exp(k):
+            # Return the k-th coefficient of the exponential.
+            k = ZZ(k)
+            if k.is_power_of(q):
+                return self._compute_coefficient_exp(k.log(q))
+            else:
+                return zero
+        return L(coeff_exp, valuation=1)
+
     def gen(self):
         r"""
         Return the generator of the Drinfeld module.
@@ -1103,6 +1236,120 @@ def j_invariant(self):
         q = self._Fq.order()
         return (g**(q+1)) / delta
 
+    @cached_method
+    def _compute_coefficient_log(self, k):
+        r"""
+        Return the `q^k`-th coefficient of the logarithm of this Drinfeld module.
+
+        TESTS::
+
+            sage: A = GF(2)['T']
+            sage: K.<T> = Frac(A)
+            sage: phi = DrinfeldModule(A, [T, 1])
+            sage: q = A.base_ring().cardinality()
+            sage: phi._compute_coefficient_log(0)
+            1
+            sage: phi._compute_coefficient_log(1)
+            1/(T^2 + T)
+            sage: phi._compute_coefficient_log(2)
+            1/(T^6 + T^5 + T^3 + T^2)
+            sage: phi._compute_coefficient_log(3)
+            1/(T^14 + T^13 + T^11 + T^10 + T^7 + T^6 + T^4 + T^3)
+        """
+        k = ZZ(k)
+        if k.is_zero():
+            return self._base.one()
+        r = self._gen.degree()
+        T = self._gen[0]
+        q = self._Fq.cardinality()
+        c = self._base.zero()
+        for i in range(k):
+            j = k - i
+            if j < r + 1:
+                c += self._compute_coefficient_log(i)*self._gen[j]**(q**i)
+        return c/(T - T**(q**k))
+
+    def logarithm(self, name='z'):
+        r"""
+        Return the logarithm of the given Drinfeld module.
+
+        By definition, the logarithm is the compositional inverse of the
+        exponential (see :meth:`exponential`). Note that the logarithm
+        is only defined when the `\mathbb{F}_q[T]`-characteristic is
+        zero.
+
+        INPUT:
+
+        - ``name`` (string, default: ``'z'``) -- the name of the
+          generator of the lazy power series ring.
+
+        OUTPUT:
+
+        A lazy power series over the base field.
+
+        EXAMPLES::
+
+            sage: A = GF(2)['T']
+            sage: K.<T> = Frac(A)
+            sage: phi = DrinfeldModule(A, [T, 1])
+            sage: log = phi.logarithm(); log
+            z + ((1/(T^2+T))*z^2) + ((1/(T^6+T^5+T^3+T^2))*z^4) + O(z^8)
+
+        The logarithm is returned as a lazy power series, meaning that
+        any of its coefficients can be computed on demands::
+
+            sage: log[2^4]
+            1/(T^30 + T^29 + T^27 + ... + T^7 + T^5 + T^4)
+            sage: log[2^5]
+            1/(T^62 + T^61 + T^59 + ... + T^8 + T^6 + T^5)
+
+        Example in higher rank::
+
+            sage: A = GF(5)['T']
+            sage: K.<T> = Frac(A)
+            sage: phi = DrinfeldModule(A, [T, T^2, T + T^2 + T^4, 1])
+            sage: phi.logarithm()
+            z + ((4*T/(T^4+4))*z^5) + O(z^8)
+
+        TESTS::
+
+            sage: A = GF(2)['T']
+            sage: K.<T> = Frac(A)
+            sage: phi = DrinfeldModule(A, [T, 1])
+            sage: q = 2
+            sage: log[2] == -1/((T**q - T))  # expected value
+            True
+            sage: log[2**2] == 1/((T**q - T)*(T**(q**2) - T))  # expected value
+            True
+            sage: log[2**3] == -1/((T**q - T)*(T**(q**2) - T)*(T**(q**3) - T))  # expected value
+            True
+
+        ::
+
+            sage: Fq.<w> = GF(3)
+            sage: A = Fq['T']
+            sage: phi = DrinfeldModule(A, [w, 1])
+            sage: phi.logarithm()
+            Traceback (most recent call last):
+            ...
+            ValueError: characteristic must be zero (=T + 2)
+        """
+        if self.category()._characteristic:
+            raise ValueError(f"characteristic must be zero (={self.characteristic()})")
+        L = LazyPowerSeriesRing(self._base, name)
+        zero = self._base.zero()
+        q = self._Fq.cardinality()
+
+        def coeff_log(k):
+            # Return the k-th coefficient of the logarithm
+            k = ZZ(k)
+            if k.is_power_of(q):
+                return self._compute_coefficient_log(k.log(q))
+            else:
+                return self._base.zero()
+        return L(coeff_log, valuation=1)
+
+
     def morphism(self):
         r"""
         Return the morphism object that defines the Drinfeld module.