gh-36269: `coefficients` for `QuasiModularFormsElement`
    
<!-- ^^^^^
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"
-->
<!-- Describe your changes here in detail -->

This patch implements `_compute` and `coefficients` methods to
`QuasiModularFormsElement`, which already exists in
`ModularFormsElement`. I made this for my own usage.
(This is my first PR for Sage, so comments are welcome!)

<!-- Why is this change required? What problem does it solve? -->
<!-- If this PR resolves an open issue, please link to it here. For
example "Fixes #12345". -->
<!-- If your change requires a documentation PR, please link it
appropriately. -->

### 📝 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! -->
<!-- Feel free to remove irrelevant items. -->

- [x] The title is concise, informative, and self-explanatory.
- [x] The description explains in detail what this PR is about.
- [x] I have linked a relevant issue or discussion.
- [x] I have created tests covering the changes.
- [x] I have updated the documentation accordingly.

### ⌛ Dependencies

<!-- List all open PRs that this PR logically depends on
- #12345: short description why this is a dependency
- #34567: ...
-->

<!-- If you're unsure about any of these, don't hesitate to ask. We're
here to help! -->
    
URL: #36269
Reported by: Seewoo Lee
Reviewer(s): David Ayotte, Seewoo Lee

src/sage/modular/quasimodform/element.py

@@ -5,9 +5,11 @@
 AUTHORS:
 
 - DAVID AYOTTE (2021-03-18): initial version
+- Seewoo Lee (2023-09): coefficients method
 """
 # ****************************************************************************
 #       Copyright (C) 2021 David Ayotte
+#                     2023 Seewoo Lee <[email protected]>
 #
 #
 # This program is free software: you can redistribute it and/or modify
@@ -772,3 +774,73 @@ def derivative(self):
         hom_comp = self.homogeneous_components()
 
         return sum(f.serre_derivative() + R(k) * u * f * E2 for k, f in hom_comp.items())
+
+    def _compute(self, X):
+        r"""
+        Compute the coefficients of `q^n` of the q-expansion of this,
+        graded quasimodular form for `n` in the list `X`.
+
+        The results are not cached.  (Use coefficients for cached results).
+
+        EXAMPLES::
+
+            sage: E2 = QuasiModularForms(1).0
+            sage: E2.q_expansion(10)
+            1 - 24*q - 72*q^2 - 96*q^3 - 168*q^4 - 144*q^5 - 288*q^6 - 192*q^7 - 360*q^8 - 312*q^9 + O(q^10)
+            sage: E2._compute([3, 6])
+            [-96, -288]
+            sage: E2._compute([])
+            []
+        """
+        if not isinstance(X, list) or not X:
+            return []
+        bound = max(X)
+        q_exp = self.q_expansion(bound + 1)
+        return [q_exp[i] for i in X]
+
+    def coefficients(self, X):
+        r"""
+        Return the coefficients of `q^n` of the q-expansion of this,
+        graded quasimodular form for `n` in the list `X`.
+
+        If X is an integer, return coefficients for indices from 1
+        to X. This method caches the result.
+
+        EXAMPLES::
+
+            sage: E2, E4 = QuasiModularForms(1).0, QuasiModularForms(1).1
+            sage: f = E2^2
+            sage: g = E2^3 * E4
+            sage: f.coefficients(10)
+            [-48, 432, 3264, 9456, 21600, 39744, 66432, 105840, 147984, 220320]
+            sage: f.coefficients([0,1])
+            [1, -48]
+            sage: f.coefficients([0,1,2,3])
+            [1, -48, 432, 3264]
+            sage: f.coefficients([2,3])
+            [432, 3264]
+            sage: g.coefficients(10)
+            [168,
+             -13608,
+             210336,
+             1805496,
+             -22562064,
+             -322437024,
+             -2063087808,
+             -9165872520,
+             -32250917496,
+             -96383477232]
+            sage: g.coefficients([3, 7])
+            [210336, -2063087808]
+        """
+        try:
+            self.__coefficients
+        except AttributeError:
+            self.__coefficients = {}
+        if isinstance(X, Integer):
+            X = list(range(1, X + 1))
+        Y = [n for n in X if n not in self.__coefficients]
+        v = self._compute(Y)
+        for i in range(len(v)):
+            self.__coefficients[Y[i]] = v[i]
+        return [self.__coefficients[x] for x in X]