diff --git a/src/doc/en/reference/references/index.rst b/src/doc/en/reference/references/index.rst
index fd3963ab275..d810bed4869 100644
--- a/src/doc/en/reference/references/index.rst
+++ b/src/doc/en/reference/references/index.rst
@@ -1044,6 +1044,11 @@ REFERENCES:
             *Guide to using plantri*, version 5.0, 2016.
             http://cs.anu.edu.au/~bdm/plantri/plantri-guide.txt
 
+.. [BM2021] Alin Bostan and Ryuhei Mori,
+            *A simple and fast algorithm for computing the N-th term of a linearly recurrent sequence*,
+            Proceedings of Symposium on Simplicity in Algorithms (SOSA), pp. 118--132, 2021.
+            :doi:`10.1137/1.9781611976496.14`
+
 .. [BMBFLR2008] A. Blondin-Massé, S. Brlek, A. Frosini, S. Labbé,
                 S. Rinaldi, *Reconstructing words from a fixed
                 palindromic length sequence*, Proc. TCS 2008, 5th IFIP
diff --git a/src/sage/rings/cfinite_sequence.py b/src/sage/rings/cfinite_sequence.py
index 26cce303581..4d9353e57b8 100644
--- a/src/sage/rings/cfinite_sequence.py
+++ b/src/sage/rings/cfinite_sequence.py
@@ -391,7 +391,6 @@ def __init__(self, parent, ogf):
             rem = num % den
             if den != 1:
                 self._a = R(num / den).list()
-                self._aa = (rem.valuation() * [0] + R(rem / den).list())[:self._deg]  # needed for _get_item_
             else:
                 self._a = num.list()
             if len(self._a) < alen:
@@ -645,33 +644,41 @@ def __getitem__(self, key):
             [0, 0, 1, 2, 3, 4, 5, 6, 7, 8]
             sage: s = C(x^3 * (1 - x)^-2); s[0:10]
             [0, 0, 0, 1, 2, 3, 4, 5, 6, 7]
+            sage: s = C(1/(1-x^1000)); s[10^18]
+            1
+            sage: s = C(1/(1-x^1000)); s[10^20]
+            1
+
+        REFERENCES:
+
+        - [BM2021]_
         """
         if isinstance(key, slice):
             m = max(key.start, key.stop)
             return [self[ii] for ii in range(*key.indices(m + 1))]
         elif isinstance(key, Integral):
-            from sage.matrix.constructor import Matrix
-            d = self._deg
-            if (self._off <= key and key < self._off + len(self._a)):
-                return self._a[key - self._off]
-            elif d == 0:
+            n = key - self._off
+            if n < 0:
                 return 0
-            (quo, rem) = self.numerator().quo_rem(self.denominator())
-            wp = quo[key - self._off]
-            if key < self._off:
-                return wp
-            A = Matrix(QQ, 1, d, self._c)
-            B = Matrix.identity(QQ, d - 1)
-            C = Matrix(QQ, d - 1, 1, 0)
-            if quo == 0:
-                off = self._off
-                V = Matrix(QQ, d, 1, self._a[:d][::-1])
+            den = self.denominator()
+            num = self.numerator()
+            if self._off >= 0:
+                num = num.shift(-self._off)
             else:
-                off = 0
-                V = Matrix(QQ, d, 1, self._aa[:d][::-1])
-            M = Matrix.block([[A], [B, C]], subdivide=False)
-
-            return wp + list(M ** (key - off) * V)[d - 1][0]
+                den = den.shift(self._off)
+            (quo, num) = num.quo_rem(den)
+            if quo.degree() < n:
+                wp = 0
+            else:
+                wp = quo[n]
+            P = self.parent().polynomial_ring()
+            x = P.gen()
+            while n:
+                nden = den(-x)
+                num = P((num * nden).list()[n % 2::2])
+                den = P((den * nden).list()[::2])
+                n //= 2
+            return wp + num[0] / den[0]
         else:
             raise TypeError("invalid argument type")