diff --git a/src/doc/en/reference/references/index.rst b/src/doc/en/reference/references/index.rst
index 52f520395d6..76d56d4cb65 100644
--- a/src/doc/en/reference/references/index.rst
+++ b/src/doc/en/reference/references/index.rst
@@ -3038,6 +3038,9 @@ REFERENCES:
 
 .. [Har1994] Frank Harary. *Graph Theory*. Reading, MA: Addison-Wesley, 1994.
 
+.. [Har2009] Harvey, David. *Efficient computation of p-adic heights*.
+             LMS J. Comput. Math. **11** (2008), 40–59.
+
 .. [Harako2020] Shuichi Harako. *The second homology group of the commutative
                 case of Kontsevich's symplectic derivation Lie algebra*.
                 Preprint, 2020, :arxiv:`2006.06064`.
@@ -4737,17 +4740,21 @@ REFERENCES:
               sharpening of the Parikh mapping*. Theoret. Informatics Appl. 35
               (2001) 551-564.
 
+.. [MST2006] Barry Mazur, William Stein, John Tate.
+             *Computation of p-adic heights and log convergence*.
+             Doc. Math. 2006, Extra Vol., 577-614.
+
 .. [MSZ2013] Michael Maschler, Solan Eilon, and Zamir Shmuel. *Game
              Theory*. Cambridge: Cambridge University Press,
              (2013). ISBN 9781107005488.
 
-.. [MT1991] Mazur, B., & Tate, J. (1991). The `p`-adic sigma
-            function. Duke Mathematical Journal, 62(3), 663-688.
+.. [MT1991] Mazur, B., & Tate, J. (1991). *The `p`-adic sigma
+            function*. Duke Mathematical Journal, **62** (3), 663-688.
 
-.. [MTT1986] \B. Mazur, J. Tate, and J. Teitelbaum, On `p`-adic
-             analogues of the conjectures of Birch and
-             Swinnerton-Dyer, Inventiones Mathematicae 84, (1986),
-             1-48.
+.. [MTT1986] \B. Mazur, J. Tate, and J. Teitelbaum,
+             *On `p`-adic analogues of the conjectures of Birch and
+             Swinnerton-Dyer*,
+             Inventiones Mathematicae **84**, (1986), 1-48.
 
 .. [Mu1997] Murty, M. Ram. *Congruences between modular forms*. In "Analytic
             Number Theory" (ed. Y. Motohashi), London Math. Soc. Lecture Notes
diff --git a/src/sage/schemes/elliptic_curves/BSD.py b/src/sage/schemes/elliptic_curves/BSD.py
index f5522de2e00..85130f0a5f9 100644
--- a/src/sage/schemes/elliptic_curves/BSD.py
+++ b/src/sage/schemes/elliptic_curves/BSD.py
@@ -3,6 +3,7 @@
 
 from sage.arith.misc import prime_divisors
 from sage.rings.integer_ring import ZZ
+from sage.rings.rational_field import QQ
 from sage.rings.infinity import Infinity
 from sage.rings.number_field.number_field import QuadraticField
 from sage.functions.other import ceil
@@ -89,11 +90,17 @@ def simon_two_descent_work(E, two_tor_rk):
         sage: from sage.schemes.elliptic_curves.BSD import simon_two_descent_work
         sage: E = EllipticCurve('14a')
         sage: simon_two_descent_work(E, E.two_torsion_rank())
+        doctest:warning
+        ...
+        DeprecationWarning: Use E.rank(algorithm="pari") instead, as this script has been ported over to pari.
+        See https://github.com/sagemath/sage/issues/35621 for details.
         (0, 0, 0, 0, [])
         sage: E = EllipticCurve('37a')
         sage: simon_two_descent_work(E, E.two_torsion_rank())
         (1, 1, 0, 0, [(0 : 0 : 1)])
     """
+    from sage.misc.superseded import deprecation
+    deprecation(35621, 'Use the two-descent in pari instead, as this script has been ported over to pari.')
     rank_lower_bd, two_sel_rk, gens = E.simon_two_descent()
     rank_upper_bd = two_sel_rk - two_tor_rk
     gens = [P for P in gens if P.additive_order() == Infinity]
@@ -140,6 +147,55 @@ def mwrank_two_descent_work(E, two_tor_rk):
     sha2_upper_bd = MWRC.selmer_rank() - two_tor_rk - rank_lower_bd
     return rank_lower_bd, rank_upper_bd, sha2_lower_bd, sha2_upper_bd, gens
 
+def pari_two_descent_work(E):
+    r"""
+    Prepare the output from pari by two-isogeny.
+
+    INPUT:
+
+    - ``E`` -- an elliptic curve
+
+    OUTPUT: A tuple of 5 elements with the first 4 being integers.
+
+    - a lower bound on the rank
+
+    - an upper bound on the rank
+
+    - a lower bound on the rank of Sha[2]
+
+    - an upper bound on the rank of Sha[2]
+
+    - a list of the generators found
+
+    EXAMPLES::
+
+        sage: from sage.schemes.elliptic_curves.BSD import pari_two_descent_work
+        sage: E = EllipticCurve('14a')
+        sage: pari_two_descent_work(E)
+        (0, 0, 0, 0, [])
+        sage: E = EllipticCurve('37a')
+        sage: pari_two_descent_work(E) # random, up to sign
+        (1, 1, 0, 0, [(0 : -1 : 1)])
+        sage: E = EllipticCurve('210e7')
+        sage: pari_two_descent_work(E)
+        (0, 2, 0, 2, [])
+        sage: E = EllipticCurve('66b3')
+        sage: pari_two_descent_work(E)
+        (0, 0, 2, 2, [])
+
+    """
+    ep = E.pari_curve()
+    lower, rank_upper_bd, s, pts = ep.ellrank()
+    gens = sorted([E.point([QQ(x[0]),QQ(x[1])], check=True) for x in pts])
+    gens = E.saturation(gens)[0]
+    # this is explained in the pari-gp documentation:
+    # s is the dimension of Sha[2]/2Sha[4],
+    # which is a lower bound for dim Sha[2]
+    # dim Sha[2] = dim Sel2 - rank E(Q) - dim tors
+    # rank_upper_bd = dim Sel_2 - dim tors - s
+    sha_upper_bd = rank_upper_bd - len(gens) + s
+    return len(gens), rank_upper_bd, s, sha_upper_bd, gens
+
 
 def native_two_isogeny_descent_work(E, two_tor_rk):
     """
@@ -254,7 +310,7 @@ def prove_BSD(E, verbosity=0, two_desc='mwrank', proof=None, secs_hi=5,
       - 2: print information about remaining primes
 
     - ``two_desc`` -- string (default ``'mwrank'``), what to use for the
-      two-descent. Options are ``'mwrank', 'simon', 'sage'``
+      two-descent. Options are ``'mwrank', 'pari', 'sage'``
 
     - ``proof`` -- bool or None (default: None, see
       proof.elliptic_curve or sage.structure.proof). If False, this
@@ -317,7 +373,7 @@ def prove_BSD(E, verbosity=0, two_desc='mwrank', proof=None, secs_hi=5,
         True for p = 3 by Kolyvagin bound
         True for p = 5 by Kolyvagin bound
         []
-        sage: E.prove_BSD(two_desc='simon')
+        sage: E.prove_BSD(two_desc='pari')
         []
 
     A rank two curve::
@@ -433,6 +489,15 @@ def prove_BSD(E, verbosity=0, two_desc='mwrank', proof=None, secs_hi=5,
         p = 2: True by 2-descent
         True for p not in {2} by Kolyvagin.
         []
+
+    ::
+
+        sage: E = EllipticCurve('66b3')
+        sage: E.prove_BSD(two_desc="pari",verbosity=1)
+        p = 2: True by 2-descent
+        True for p not in {2} by Kolyvagin.
+        []
+
     """
     if proof is None:
         from sage.structure.proof.proof import get_flag
@@ -461,8 +526,8 @@ def prove_BSD(E, verbosity=0, two_desc='mwrank', proof=None, secs_hi=5,
 
     if two_desc == 'mwrank':
         M = mwrank_two_descent_work(BSD.curve, BSD.two_tor_rk)
-    elif two_desc == 'simon':
-        M = simon_two_descent_work(BSD.curve, BSD.two_tor_rk)
+    elif two_desc == 'pari':
+        M = pari_two_descent_work(BSD.curve)
     elif two_desc == 'sage':
         M = native_two_isogeny_descent_work(BSD.curve, BSD.two_tor_rk)
     else:
diff --git a/src/sage/schemes/elliptic_curves/constructor.py b/src/sage/schemes/elliptic_curves/constructor.py
index eae54a05fcd..770295962dd 100644
--- a/src/sage/schemes/elliptic_curves/constructor.py
+++ b/src/sage/schemes/elliptic_curves/constructor.py
@@ -308,7 +308,7 @@ class EllipticCurveFactory(UniqueFactory):
         TypeError: invalid input to EllipticCurve constructor
     """
     def create_key_and_extra_args(self, x=None, y=None, j=None, minimal_twist=True, **kwds):
-        """
+        r"""
         Return a ``UniqueFactory`` key and possibly extra parameters.
 
         INPUT: See the documentation for :class:`EllipticCurveFactory`.
@@ -329,7 +329,7 @@ def create_key_and_extra_args(self, x=None, y=None, j=None, minimal_twist=True,
             sage: EllipticCurve.create_key_and_extra_args(j=8000)
             ((Rational Field, (0, 1, 0, -3, 1)), {})
 
-        When constructing a curve over `\\QQ` from a Cremona or LMFDB
+        When constructing a curve over `\QQ` from a Cremona or LMFDB
         label, the invariants from the database are returned as
         ``extra_args``::
 
@@ -454,7 +454,7 @@ def create_key_and_extra_args(self, x=None, y=None, j=None, minimal_twist=True,
         return (R, tuple(R(a) for a in x)), kwds
 
     def create_object(self, version, key, **kwds):
-        """
+        r"""
         Create an object from a ``UniqueFactory`` key.
 
         EXAMPLES::
@@ -466,7 +466,7 @@ def create_object(self, version, key, **kwds):
         .. NOTE::
 
             Keyword arguments are currently only passed to the
-            constructor for elliptic curves over `\\QQ`; elliptic
+            constructor for elliptic curves over `\QQ`; elliptic
             curves over other fields do not support them.
         """
         R, x = key
@@ -494,7 +494,7 @@ def create_object(self, version, key, **kwds):
 
 
 def EllipticCurve_from_Weierstrass_polynomial(f):
-    """
+    r"""
     Return the elliptic curve defined by a cubic in (long) Weierstrass
     form.
 
@@ -534,7 +534,7 @@ def EllipticCurve_from_Weierstrass_polynomial(f):
     return EllipticCurve(coefficients_from_Weierstrass_polynomial(f))
 
 def coefficients_from_Weierstrass_polynomial(f):
-    """
+    r"""
     Return the coefficients `[a_1, a_2, a_3, a_4, a_6]` of a cubic in
     Weierstrass form.
 
@@ -582,7 +582,7 @@ def coefficients_from_Weierstrass_polynomial(f):
 
 
 def EllipticCurve_from_c4c6(c4, c6):
-    """
+    r"""
     Return an elliptic curve with given `c_4` and
     `c_6` invariants.
 
@@ -664,7 +664,7 @@ def EllipticCurve_from_j(j, minimal_twist=True):
 
 
 def coefficients_from_j(j, minimal_twist=True):
-    """
+    r"""
     Return Weierstrass coefficients `(a_1, a_2, a_3, a_4, a_6)` for an
     elliptic curve with given `j`-invariant.
 
@@ -680,7 +680,7 @@ def coefficients_from_j(j, minimal_twist=True):
         sage: coefficients_from_j(1)
         [1, 0, 0, 36, 3455]
 
-    The ``minimal_twist`` parameter (ignored except over `\\QQ` and
+    The ``minimal_twist`` parameter (ignored except over `\QQ` and
     True by default) controls whether or not a minimal twist is
     computed::
 
@@ -1238,7 +1238,7 @@ def EllipticCurve_from_cubic(F, P=None, morphism=True):
 
 
 def tangent_at_smooth_point(C,P):
-    """Return the tangent at the smooth point `P` of projective curve `C`.
+    r"""Return the tangent at the smooth point `P` of projective curve `C`.
 
     INPUT:
 
@@ -1278,7 +1278,7 @@ def tangent_at_smooth_point(C,P):
         return C.tangents(P,factor=False)[0]
 
 def chord_and_tangent(F, P):
-    """Return the third point of intersection of a cubic with the tangent at one point.
+    r"""Return the third point of intersection of a cubic with the tangent at one point.
 
     INPUT:
 
@@ -1352,7 +1352,7 @@ def chord_and_tangent(F, P):
 
 
 def projective_point(p):
-    """
+    r"""
     Return equivalent point with denominators removed
 
     INPUT:
@@ -1384,7 +1384,7 @@ def projective_point(p):
 
 
 def are_projectively_equivalent(P, Q, base_ring):
-    """
+    r"""
     Test whether ``P`` and ``Q`` are projectively equivalent.
 
     INPUT:
@@ -1409,7 +1409,7 @@ def are_projectively_equivalent(P, Q, base_ring):
 
 def EllipticCurves_with_good_reduction_outside_S(S=[], proof=None, verbose=False):
     r"""
-    Return a sorted list of all elliptic curves defined over `Q`
+    Return a sorted list of all elliptic curves defined over `\QQ`
     with good reduction outside the set `S` of primes.
 
     INPUT:
diff --git a/src/sage/schemes/elliptic_curves/ell_rational_field.py b/src/sage/schemes/elliptic_curves/ell_rational_field.py
index 7aa20a1c2c1..b5f5de2bd3f 100644
--- a/src/sage/schemes/elliptic_curves/ell_rational_field.py
+++ b/src/sage/schemes/elliptic_curves/ell_rational_field.py
@@ -1,5 +1,5 @@
 # -*- coding: utf-8 -*-
-"""
+r"""
 Elliptic curves over the rational numbers
 
 AUTHORS:
@@ -200,7 +200,7 @@ def __init__(self, ainvs, **kwds):
             self._set_torsion_order(kwds['torsion_order'])
 
     def _set_rank(self, r):
-        """
+        r"""
         Internal function to set the cached rank of this elliptic curve to
         ``r``.
 
@@ -221,7 +221,7 @@ def _set_rank(self, r):
         self.__rank = (Integer(r), True)
 
     def _set_torsion_order(self, t):
-        """
+        r"""
         Internal function to set the cached torsion order of this elliptic
         curve to ``t``.
 
@@ -242,7 +242,7 @@ def _set_torsion_order(self, t):
         self.__torsion_order = Integer(t)
 
     def _set_cremona_label(self, L):
-        """
+        r"""
         Internal function to set the cached label of this elliptic curve to
         ``L``.
 
@@ -267,7 +267,7 @@ def _set_cremona_label(self, L):
         self.__cremona_label = L
 
     def _set_conductor(self, N):
-        """
+        r"""
         Internal function to set the cached conductor of this elliptic
         curve to ``N.``
 
@@ -288,7 +288,7 @@ def _set_conductor(self, N):
         self.__conductor_pari = Integer(N)
 
     def _set_modular_degree(self, deg):
-        """
+        r"""
         Internal function to set the cached modular degree of this elliptic
         curve to ``deg``.
 
@@ -309,7 +309,7 @@ def _set_modular_degree(self, deg):
         self.__modular_degree = Integer(deg)
 
     def _set_gens(self, gens):
-        """
+        r"""
         Internal function to set the cached generators of this elliptic
         curve to ``gens``.
 
@@ -382,7 +382,7 @@ def is_p_integral(self, p):
         return bool(misc.mul([x.valuation(p) >= 0 for x in self.ainvs()]))
 
     def is_integral(self):
-        """
+        r"""
         Return ``True`` if this elliptic curve has integral coefficients (in
         Z).
 
@@ -476,7 +476,7 @@ def mwrank(self, options=''):
         return mwrank(list(self.a_invariants()))
 
     def conductor(self, algorithm="pari"):
-        """
+        r"""
         Return the conductor of the elliptic curve.
 
         INPUT:
@@ -576,7 +576,7 @@ def conductor(self, algorithm="pari"):
     ####################################################################
 
     def pari_curve(self):
-        """
+        r"""
         Return the PARI curve corresponding to this elliptic curve.
 
         EXAMPLES::
@@ -629,7 +629,7 @@ def pari_curve(self):
             return self._pari_curve
 
     def pari_mincurve(self):
-        """
+        r"""
         Return the PARI curve corresponding to a minimal model for this
         elliptic curve.
 
@@ -653,7 +653,7 @@ def pari_mincurve(self):
 
     @cached_method
     def database_attributes(self):
-        """
+        r"""
         Return a dictionary containing information about ``self`` in
         the elliptic curve database.
 
@@ -688,7 +688,7 @@ def database_attributes(self):
             raise LookupError("Cremona database does not contain entry for " + repr(self))
 
     def database_curve(self):
-        """
+        r"""
         Return the curve in the elliptic curve database isomorphic to this
         curve, if possible. Otherwise raise a ``LookupError`` exception.
 
@@ -760,7 +760,7 @@ def Np(self, p):
     #  Access to mwrank
     ####################################################################
     def mwrank_curve(self, verbose=False):
-        """
+        r"""
         Construct an mwrank_EllipticCurve from this elliptic curve
 
         The resulting mwrank_EllipticCurve has available methods from John
@@ -792,7 +792,7 @@ def two_descent(self, verbose=True,
                     second_limit=8,
                     n_aux=-1,
                     second_descent=1):
-        """
+        r"""
         Compute 2-descent data for this curve.
 
         INPUT:
@@ -1362,7 +1362,7 @@ def f(r):
         return f
 
     def pollack_stevens_modular_symbol(self, sign=0, implementation='eclib'):
-        """
+        r"""
         Create the modular symbol attached to the elliptic curve,
         suitable for overconvergent calculations.
 
@@ -1777,6 +1777,13 @@ def simon_two_descent(self, verbose=0, lim1=5, lim3=50, limtriv=3,
         Return lower and upper bounds on the rank of the Mordell-Weil
         group `E(\QQ)` and a list of points of infinite order.
 
+        .. WARNING::
+
+            This function is deprecated as the functionality of
+            Simon's script for elliptic curves over the rationals
+            has been ported over to pari.
+            Use :meth:`.rank` with the keyword ``algorithm='pari'`` instead.
+
         INPUT:
 
         - ``verbose`` -- 0, 1, 2, or 3 (default: 0), the verbosity level
@@ -1829,6 +1836,10 @@ def simon_two_descent(self, verbose=0, lim1=5, lim3=50, limtriv=3,
 
             sage: E = EllipticCurve('11a1')
             sage: E.simon_two_descent()
+            doctest:warning
+            ...
+            DeprecationWarning: Use E.rank(algorithm="pari") instead, as this script has been ported over to pari.
+            See https://github.com/sagemath/sage/issues/35621 for details.
             (0, 0, [])
             sage: E = EllipticCurve('37a1')
             sage: E.simon_two_descent()
@@ -1909,6 +1920,9 @@ def simon_two_descent(self, verbose=0, lim1=5, lim3=50, limtriv=3,
             sage: E.selmer_rank()  # uses mwrank
             1
         """
+        from sage.misc.superseded import deprecation
+        deprecation(35621, 'Use E.rank(algorithm="pari") instead, as this script has been ported over to pari.')
+
         t = EllipticCurve_number_field.simon_two_descent(self, verbose=verbose,
                                                          lim1=lim1, lim3=lim3, limtriv=limtriv,
                                                          maxprob=maxprob, limbigprime=limbigprime,
@@ -1978,8 +1992,9 @@ def three_selmer_rank(self, algorithm='UseSUnits'):
     def rank(self, use_database=True, verbose=False,
              only_use_mwrank=True,
              algorithm='mwrank_lib',
-             proof=None):
-        """
+             proof=None,
+             pari_effort=0):
+        r"""
         Return the rank of this elliptic curve, assuming no conjectures.
 
         If we fail to provably compute the rank, raises a RuntimeError
@@ -1999,6 +2014,8 @@ def rank(self, use_database=True, verbose=False,
 
           - ``'mwrank_lib'`` -- call mwrank c library
 
+          - ``'pari'`` -- call ellrank in pari
+
         - ``only_use_mwrank`` -- (default: ``True``) if ``False`` try
           using analytic rank methods first
 
@@ -2006,9 +2023,15 @@ def rank(self, use_database=True, verbose=False,
           ``proof.elliptic_curve`` or ``sage.structure.proof``); note that
           results obtained from databases are considered ``proof=True``
 
+        - ``pari_effort`` -- (default: 0) parameter used in when
+          the algorithm ``pari`` is chosen. It measure of the effort
+          done to find rational points. Values up to 10 can be chosen;
+          the running times increase roughly like the cube of the
+          effort value.
+
         OUTPUT: the rank of the elliptic curve as :class:`Integer`
 
-        IMPLEMENTATION: Uses L-functions, mwrank, and databases.
+        IMPLEMENTATION: Uses L-functions, mwrank, pari, and databases.
 
         EXAMPLES::
 
@@ -2024,7 +2047,7 @@ def rank(self, use_database=True, verbose=False,
             4
             sage: EllipticCurve([0, 0, 1, -79, 342]).rank(proof=False)
             5
-            sage: EllipticCurve([0, 0, 1, -79, 342]).simon_two_descent()[0]  # long time (7s on sage.math, 2012)
+            sage: EllipticCurve([0, 0, 1, -79, 342]).rank(algorithm="pari")
             5
 
         Examples with denominators in defining equations::
@@ -2038,7 +2061,7 @@ def rank(self, use_database=True, verbose=False,
             sage: E.minimal_model().rank()
             1
 
-        A large example where mwrank doesn't determine the result with certainty::
+        A large example where mwrank doesn't determine the result with certainty, but pari does::
 
             sage: EllipticCurve([1,0,0,0,37455]).rank(proof=False)
             0
@@ -2046,6 +2069,8 @@ def rank(self, use_database=True, verbose=False,
             Traceback (most recent call last):
             ...
             RuntimeError: rank not provably correct (lower bound: 0)
+            sage: EllipticCurve([1,0,0,0,37455]).rank(algorithm="pari")
+            0
 
         TESTS::
 
@@ -2054,6 +2079,14 @@ def rank(self, use_database=True, verbose=False,
             ...
             ValueError: unknown algorithm 'garbage'
 
+        An example to check if the points are saturated::
+
+            sage: E = EllipticCurve([0,0, 1, -7, 6])
+            sage: E.gens(use_database=False, algorithm="pari") # random
+            [(2 : 0 : 1), (-1 : 3 : 1), (11 : 35 : 1)]
+            sage: E.saturation(_)[1]
+            1
+
         Since :trac:`23962`, the default is to use the Cremona
         database. We also check that the result is cached correctly::
 
@@ -2068,6 +2101,15 @@ def rank(self, use_database=True, verbose=False,
             0
             sage: E._EllipticCurve_rational_field__rank
             (0, True)
+
+        This example has Sha = Z/4 x Z/4 and the rank cannot be
+        determined using pari only::
+
+            sage: E =EllipticCurve([-113^2,0])
+            sage: E.rank(use_database=False, verbose=False, algorithm="pari")
+            Traceback (most recent call last):
+            ...
+            RuntimeError: rank not provably correct (lower bound: 0, upper bound:2). Hint: increase pari_effort.
         """
         if proof is None:
             from sage.structure.proof.proof import get_flag
@@ -2168,6 +2210,27 @@ def rank(self, use_database=True, verbose=False,
             self.__rank = (rank, proof)
             return rank
 
+        if algorithm == 'pari':
+            ep = self.pari_curve()
+            # if we know already some points in _known_points
+            # we can give them to pari to speed it up
+            kpts = [ [x[0],x[1]] for x in self._known_points ]
+            lower, upper, s, pts = ep.ellrank(pari_effort, kpts)
+            ge = sorted([self.point([QQ(x[0]),QQ(x[1])], check=True) for x in pts])
+            ge = self.saturation(ge)[0]
+            self._known_points = ge
+            # note that lower is only a conjectural
+            # lower bound for the rank, the only
+            # proven lower bound is #ge.
+            if len(ge) == upper:
+                verbose_verbose(f"rank {upper} unconditionally determined by pari")
+                rank = Integer(upper)
+                self.__rank = (rank, True)
+                self.__gens = (ge, True)
+                return rank
+            else:
+                verbose_verbose(f"Warning -- rank could not be determined by pari; ellrank returned {lower=}, {upper=}, {s=}, {pts=}", level=1)
+                raise RuntimeError(f"rank not provably correct (lower bound: {len(ge)}, upper bound:{upper}). Hint: increase pari_effort.")
         raise ValueError("unknown algorithm {!r}".format(algorithm))
 
     def gens(self, proof=None, **kwds):
@@ -2194,6 +2257,8 @@ def gens(self, proof=None, **kwds):
 
           - ``'mwrank_lib'`` -- call mwrank C library
 
+          - ``'pari'`` -- use ellrank in pari
+
         - ``only_use_mwrank`` -- bool (default True) if False, first
           attempts to use more naive, natively implemented methods
 
@@ -2207,6 +2272,12 @@ def gens(self, proof=None, **kwds):
           points found by two-descent in the Mordell-Weil group is
           greater than this, a warning message will be displayed.
 
+        - ``pari_effort`` -- (default: 0) parameter used in when
+          the algorithm ``pari`` is chosen. It measure of the effort
+          done to find rational points. Values up to 10 can be chosen,
+          the running times increase roughly like the cube of the
+          effort value.
+
         OUTPUT:
 
         - ``generators`` -- list of generators for the Mordell-Weil
@@ -2218,13 +2289,18 @@ def gens(self, proof=None, **kwds):
             :meth:`~gens_certain` method to find out afterwards
             whether the generators were proved.
 
-        IMPLEMENTATION: Uses Cremona's mwrank C library.
+        IMPLEMENTATION: Uses Cremona's mwrank C++ library or ellrank in pari.
 
         EXAMPLES::
 
             sage: E = EllipticCurve('389a')
             sage: E.gens()                 # random output
             [(-1 : 1 : 1), (0 : 0 : 1)]
+            sage: E.gens(algorithm="pari")    # random output
+            [(5/4 : 5/8 : 1), (0 : 0 : 1)]
+            sage: E = EllipticCurve([0,2429469980725060,0,275130703388172136833647756388,0])
+            sage: len(E.gens(algorithm="pari"))
+            14
 
         A non-integral example::
 
@@ -2240,6 +2316,9 @@ def gens(self, proof=None, **kwds):
              over Rational Field
             sage: E1.gens() # random (if database not used)
             [(-400 : 8000 : 1), (0 : -8000 : 1)]
+            sage: E1.gens(algorithm="pari")   #random
+            [(-400 : 8000 : 1), (0 : -8000 : 1)]
+
         """
         if proof is None:
             from sage.structure.proof.proof import get_flag
@@ -2266,7 +2345,8 @@ def _compute_gens(self, proof,
                       only_use_mwrank=True,
                       use_database=True,
                       descent_second_limit=12,
-                      sat_bound=1000):
+                      sat_bound=1000,
+                      pari_effort=0):
         r"""
         Return generators for the Mordell-Weil group `E(Q)` *modulo*
         torsion.
@@ -2289,6 +2369,25 @@ def _compute_gens(self, proof,
             sage: gens, proved = E._compute_gens(proof=False)
             sage: proved
             True
+
+        TESTS::
+
+            sage: E = EllipticCurve([-127^2,0])
+            sage: E.gens(use_database=False, algorithm="pari")
+            Traceback (most recent call last):
+            ...
+            RuntimeError: generators could not be determined. So far we found []. Hint: increase pari_effort.
+            sage: E.gens(use_database=False, algorithm="pari",pari_effort=4)
+            [(611429153205013185025/9492121848205441 : 15118836457596902442737698070880/924793900700594415341761 : 1)]
+
+            sage: E = EllipticCurve([-157^2,0])
+            sage: E.gens(use_database=False, algorithm="pari")
+            Traceback (most recent call last):
+            ...
+            RuntimeError: generators could not be determined. So far we found []. Hint: increase pari_effort.
+            sage: E.gens(use_database=False, algorithm="pari",pari_effort=10) # long time
+            [(-166136231668185267540804/2825630694251145858025 : 167661624456834335404812111469782006/150201095200135518108761470235125 : 1)]
+
         """
         # If the optional extended database is installed and an
         # isomorphic curve is in the database then its gens will be
@@ -2333,7 +2432,29 @@ def _compute_gens(self, proof,
             except RuntimeError:
                 pass
         # end if (not_use_mwrank)
-        if algorithm == "mwrank_lib":
+        if algorithm == "pari":
+            ep = self.pari_curve()
+            # if we know already some points in _known_points
+            # we can give them to pari to speed it up
+            kpts = [ [x[0],x[1]] for x in self._known_points ]
+            lower, upper, s, pts = ep.ellrank(pari_effort, kpts)
+            ge = sorted([self.point([QQ(x[0]),QQ(x[1])], check=True) for x in pts])
+            ge = self.saturation(ge)[0]
+            self._known_points = ge
+            # note that lower is only a conjectural
+            # lower bound for the rank, the only
+            # proven lower bound is #ge.
+            if len(ge) == upper:
+                verbose_verbose(f"rank {upper} unconditionally determined by pari")
+                rank = Integer(upper)
+                self.__rank = (rank, True)
+                if len(ge) == rank:
+                    self.__gens = (ge, True)
+                    return ge, True
+            # cases when we did not find all points
+            verbose_verbose(f"Warning -- generators could not be determined by pari; ellrank returned {lower=}, {upper=}, {s=}, {pts=}", level=1)
+            raise RuntimeError(f"generators could not be determined. So far we found {ge}. Hint: increase pari_effort.")
+        elif algorithm == "mwrank_lib":
             verbose_verbose("Calling mwrank C++ library.")
             if not self.is_integral():
                 xterm = 1
@@ -2398,7 +2519,7 @@ def _compute_gens(self, proof,
         return G, proved
 
     def gens_certain(self):
-        """
+        r"""
         Return ``True`` if the generators have been proven correct.
 
         EXAMPLES::
@@ -2853,86 +2974,118 @@ def point_search(self, height_limit, verbose=False, rank_bound=None):
             points = self.saturation(points, verbose=verbose)[0]
         return points
 
-    def selmer_rank(self):
+    def selmer_rank(self, algorithm="pari"):
         r"""
-        The rank of the 2-Selmer group of the curve.
+        Return the rank of the 2-Selmer group of the curve.
 
-        EXAMPLES: The following is the curve 960D1, which has rank 0, but
-        Sha of order 4.
+        INPUT:
 
-        ::
+        - ``algorithm`` -- (default:``'pari'``)
+          either ``'pari'`` or ``'mwrank'``
 
-            sage: E = EllipticCurve([0, -1, 0, -900, -10098])
-            sage: E.selmer_rank()
-            3
-
-        Here the Selmer rank is equal to the 2-torsion rank (=1) plus
-        the 2-rank of Sha (=2), and the rank itself is zero::
+        EXAMPLES:
+        This example has rank 1, Sha[2] of order 4 and
+        a single rational 2-torsion point::
 
-            sage: E.rank()
-            0
+            sage: E = EllipticCurve([1, 1, 1, 508, -2551])
+            sage: E.selmer_rank()
+            4
+            sage: E.selmer_rank(algorithm="mwrank")
+            4
 
-        In contrast, for the curve 571A, also with rank 0 and Sha of
-        order 4, we get a worse bound::
+        The following is the curve 960d1, which has rank 0, but
+        Sha of order 4::
 
-            sage: E = EllipticCurve([0, -1, 1, -929, -10595])
+            sage: E = EllipticCurve([0, -1, 0, -900, -10098])
             sage: E.selmer_rank()
-            2
-            sage: E.rank_bound()
-            2
+            3
+            sage: E.selmer_rank(algorithm="mwrank")
+            3
 
-        To establish that the rank is in fact 0 in this case, we would
-        need to carry out a higher descent::
+        This curve has rank 1, and 4 elements in Sha[2].
+        Yet the order of Sha is 16, so that group is the product
+        of two cyclic groups of order 4::
 
-            sage: E.three_selmer_rank() # optional - magma
-            0
+            sage: E = EllipticCurve([1, 0, 0, -150752, -22541610])
+            sage: E.selmer_rank()
+            4
 
-        Or use the L-function to compute the analytic rank::
+        Instead in this last example of rank 0, Sha is a product of four cyclic groups of order 2::
 
-            sage: E.rank(only_use_mwrank=False)
+            sage: E = EllipticCurve([1, 0, 0, -49280, -4214808])
+            sage: E.selmer_rank()
+            5
+            sage: E.rank()
             0
         """
         try:
             return self.__selmer_rank
         except AttributeError:
-            C = self.mwrank_curve()
-            self.__selmer_rank = C.selmer_rank()
-            return self.__selmer_rank
+            if algorithm=="pari":
+                ep = self.pari_curve()
+                lower, upper, s, pts = ep.ellrank()
+                tor = self.two_torsion_rank()
+                return upper + tor + s
+            elif algorithm=="mwrank":
+                C = self.mwrank_curve()
+                self.__selmer_rank = C.selmer_rank()
+                return self.__selmer_rank
+            else:
+                raise ValueError(f"unknown {algorithm=}")
 
-    def rank_bound(self):
+    def rank_bound(self, algorithm="pari"):
         r"""
-        Upper bound on the rank of the curve, computed using
-        2-descent.
+        Return the upper bound on the rank of the curve,
+        computed using a 2-descent.
+
+        INPUT:
+
+        - ``algorithm`` -- (default:``'pari'``)
+          either ``'pari'`` or ``'mwrank'``
 
         In many cases, this is the actual rank of the
-        curve.  If the curve has no 2-torsion it is the same as the
-        2-selmer rank.
+        curve.
 
-        EXAMPLES: The following is the curve 960D1, which has rank 0, but
-        Sha of order 4.
+        EXAMPLES::
 
-        ::
+            sage: E = EllipticCurve("389a1")
+            sage: E.rank_bound()
+            2
 
-            sage: E = EllipticCurve([0, -1, 0, -900, -10098])
+        The following is the curve 571a1, which has
+        rank 0, but Sha of order 4, yet pari, using the Cassels
+        pairing is able to show that the rank is 0.
+        The 2-descent in mwrank only determines a weaker upper bound::
+
+            sage: E = EllipticCurve([0, -1, 1, -929, -10595])
             sage: E.rank_bound()
             0
+            sage: E.rank_bound(algorithm="mwrank")
+            2
 
-        It gives 0 instead of 2, because it knows Sha is nontrivial. In
-        contrast, for the curve 571A, also with rank 0 and Sha of order 4,
-        we get a worse bound::
+        In the following last example, both algorithm only determine a rank bound larger than the actual rank::
 
-            sage: E = EllipticCurve([0, -1, 1, -929, -10595])
+            sage: E = EllipticCurve([1, 1, 1, -896670, -327184905])
             sage: E.rank_bound()
             2
-            sage: E.rank(only_use_mwrank=False)   # uses L-function
+            sage: E.rank_bound(algorithm="mwrank")
+            2
+            sage: E.rank(only_use_mwrank=False) # uses L-function
             0
         """
         try:
             return self.__rank_bound
         except AttributeError:
-            C = self.mwrank_curve()
-            self.__rank_bound = C.rank_bound()
-            return self.__rank_bound
+            if algorithm=="pari":
+                ep = self.pari_curve()
+                lower, upper, s, pts = ep.ellrank()
+                return upper
+            elif algorithm=="mwrank":
+                C = self.mwrank_curve()
+                self.__rank_bound = C.rank_bound()
+                return self.__rank_bound
+            else:
+                raise ValueError(f"unknown {algorithm=}")
 
     def an(self, n):
         r"""
@@ -2948,7 +3101,7 @@ def an(self, n):
         return Integer(self.pari_mincurve().ellak(n))
 
     def ap(self, p):
-        """
+        r"""
         The ``p``-th Fourier coefficient of the modular form corresponding to
         this elliptic curve, where ``p`` is prime.
 
@@ -3009,7 +3162,7 @@ def is_minimal(self):
         return self.ainvs() == self.minimal_model().ainvs()
 
     def is_p_minimal(self, p):
-        """
+        r"""
         Tests if curve is ``p``-minimal at a given prime ``p``.
 
         INPUT:
@@ -3188,7 +3341,7 @@ def tamagawa_exponent(self, p):
         return 4
 
     def tamagawa_product(self):
-        """
+        r"""
         Return the product of the Tamagawa numbers.
 
         EXAMPLES::
@@ -3204,7 +3357,7 @@ def tamagawa_product(self):
             return self.__tamagawa_product
 
     def real_components(self):
-        """
+        r"""
         Return the number of real components.
 
         EXAMPLES::
@@ -3361,7 +3514,7 @@ def elliptic_exponential(self, z, embedding=None):
         return self.period_lattice().elliptic_exponential(z)
 
     def lseries(self):
-        """
+        r"""
         Return the L-series of this elliptic curve.
 
         Further documentation is available for the functions which apply to
@@ -3911,7 +4064,7 @@ def congruence_number(self, M=1):
                 return self._generalized_congmod_numbers(M)["congnum"]
 
     def cremona_label(self, space=False):
-        """
+        r"""
         Return the Cremona label associated to (the minimal model) of this
         curve, if it is known. If not, raise a ``LookupError`` exception.
 
@@ -4003,7 +4156,7 @@ def reduction(self,p):
         return self.change_ring(rings.GF(p))
 
     def torsion_order(self):
-        """
+        r"""
         Return the order of the torsion subgroup.
 
         EXAMPLES::
@@ -4107,7 +4260,7 @@ def torsion_subgroup(self):
         return self.__torsion_subgroup
 
     def torsion_points(self):
-        """
+        r"""
         Return the torsion points of this elliptic curve as a sorted
         list.
 
@@ -4375,7 +4528,7 @@ def has_rational_cm(self, field=None):
             raise ValueError("Error in has_rational_cm: %s is not an extension field of QQ" % field)
 
     def quadratic_twist(self, D):
-        """
+        r"""
         Return the global minimal model of the quadratic twist of this
         curve by ``D``.
 
@@ -4746,7 +4899,7 @@ def isogenies_prime_degree(self, l=None):
             return isogs
 
     def is_isogenous(self, other, proof=True, maxp=200):
-        """
+        r"""
         Return whether or not self is isogenous to other.
 
         INPUT:
@@ -4824,7 +4977,7 @@ def is_isogenous(self, other, proof=True, maxp=200):
             return E2 in E1.isogeny_class().curves
 
     def isogeny_degree(self, other):
-        """
+        r"""
         Return the minimal degree of an isogeny between ``self`` and
         ``other``.
 
@@ -4949,7 +5102,7 @@ def isogeny_degree(self, other):
 #         """
 
     def optimal_curve(self):
-        """
+        r"""
         Given an elliptic curve that is in the installed Cremona
         database, return the optimal curve isogenous to it.
 
@@ -5289,7 +5442,7 @@ def is_ordinary(self, p, ell=None):
         return self.ap(ell) % p != 0
 
     def is_good(self, p, check=True):
-        """
+        r"""
         Return ``True`` if ``p`` is a prime of good reduction for `E`.
 
         INPUT:
@@ -5314,7 +5467,7 @@ def is_good(self, p, check=True):
         return self.conductor() % p != 0
 
     def is_supersingular(self, p, ell=None):
-        """
+        r"""
         Return ``True`` precisely when p is a prime of good reduction and the
         mod-``p`` representation attached to this elliptic curve is
         supersingular at ell.
@@ -5344,7 +5497,7 @@ def is_supersingular(self, p, ell=None):
         return self.is_good(p) and not self.is_ordinary(p, ell)
 
     def supersingular_primes(self, B):
-        """
+        r"""
         Return a list of all supersingular primes for this elliptic curve
         up to and possibly including B.
 
@@ -5379,7 +5532,7 @@ def supersingular_primes(self, B):
                       [P[i] for i in range(2,len(v)) if v[i] == 0 and N % P[i] != 0]
 
     def ordinary_primes(self, B):
-        """
+        r"""
         Return a list of all ordinary primes for this elliptic curve up to
         and possibly including B.
 
@@ -5463,7 +5616,7 @@ def eval_modular_form(self, points, order):
     ########################################################################
 
     def sha(self):
-        """
+        r"""
         Return an object of class
         'sage.schemes.elliptic_curves.sha_tate.Sha' attached to this
         elliptic curve.
@@ -5528,7 +5681,7 @@ def sha(self):
     matrix_of_frobenius = padics.matrix_of_frobenius
 
     def mod5family(self):
-        """
+        r"""
         Return the family of all elliptic curves with the same mod-5
         representation as ``self``.
 
@@ -5836,7 +5989,7 @@ def integral_x_coords_in_interval(self,xmin,xmax):
     prove_BSD = BSD.prove_BSD
 
     def integral_points(self, mw_base='auto', both_signs=False, verbose=False):
-        """
+        r"""
         Compute all integral points (up to sign) on this elliptic curve.
 
         INPUT:
@@ -6233,7 +6386,7 @@ def point_preprocessing(free,tor):
         return int_points
 
     def S_integral_points(self, S, mw_base='auto', both_signs=False, verbose=False, proof=None):
-        """
+        r"""
         Compute all S-integral points (up to sign) on this elliptic curve.
 
         INPUT:
@@ -6910,7 +7063,7 @@ def S_integral_x_coords_with_abs_bounded_by(abs_bound):
 
 
 def cremona_curves(conductors):
-    """
+    r"""
     Return iterator over all known curves (in database) with conductor
     in the list of conductors.
 
@@ -6943,7 +7096,7 @@ def cremona_curves(conductors):
     return sage.databases.cremona.CremonaDatabase().iter(conductors)
 
 def cremona_optimal_curves(conductors):
-    """
+    r"""
     Return iterator over all known optimal curves (in database) with
     conductor in the list of conductors.
 
diff --git a/src/sage/schemes/elliptic_curves/ell_torsion.py b/src/sage/schemes/elliptic_curves/ell_torsion.py
index 5cd89e9c2b5..9bd2b1e8f96 100644
--- a/src/sage/schemes/elliptic_curves/ell_torsion.py
+++ b/src/sage/schemes/elliptic_curves/ell_torsion.py
@@ -140,7 +140,7 @@ def __init__(self, E):
 
         INPUT:
 
-        - ``E`` -- An elliptic curve defined over a number field (including `\Q`)
+        - ``E`` -- An elliptic curve defined over a number field (including `\QQ`)
 
         EXAMPLES::
 
@@ -212,7 +212,7 @@ def __init__(self, E):
                                                     [T1, T2], structure)
 
     def _repr_(self):
-        """
+        r"""
         String representation of an instance of the EllipticCurveTorsionSubgroup class.
 
         EXAMPLES::
@@ -242,7 +242,7 @@ def __richcmp__(self, other, op):
         return richcmp(self.__E, other.__E, op)
 
     def curve(self):
-        """
+        r"""
         Return the curve of this torsion subgroup.
 
         EXAMPLES::
@@ -259,7 +259,7 @@ def curve(self):
 
     @cached_method
     def points(self):
-        """
+        r"""
         Return a list of all the points in this torsion subgroup.
 
         The list is cached.
diff --git a/src/sage/schemes/elliptic_curves/hom.py b/src/sage/schemes/elliptic_curves/hom.py
index 2ec9d2fd747..a425b8abd00 100644
--- a/src/sage/schemes/elliptic_curves/hom.py
+++ b/src/sage/schemes/elliptic_curves/hom.py
@@ -56,7 +56,7 @@ def __init__(self, *args, **kwds):
               From: Elliptic Curve defined by y^2 = x^3 + 5*x + 5 over Finite Field in z2 of size 257^2
               To:   Elliptic Curve defined by y^2 = x^3 + 151*x + 22 over Finite Field in z2 of size 257^2
             sage: E.isogeny(P, algorithm='velusqrt')  # indirect doctest
-            Elliptic-curve isogeny (using Îlu) of degree 127:
+            Elliptic-curve isogeny (using square-root Vélu) of degree 127:
               From: Elliptic Curve defined by y^2 = x^3 + 5*x + 5 over Finite Field in z2 of size 257^2
               To:   Elliptic Curve defined by y^2 = x^3 + 119*x + 231 over Finite Field in z2 of size 257^2
             sage: E.montgomery_model(morphism=True)   # indirect doctest
@@ -75,7 +75,7 @@ def __init__(self, *args, **kwds):
             self._domain._fetch_cached_order(self._codomain)
 
     def _repr_type(self):
-        """
+        r"""
         Return a textual representation of what kind of morphism
         this is. Used by :meth:`Morphism._repr_`.
 
@@ -89,7 +89,7 @@ def _repr_type(self):
 
     @staticmethod
     def _composition_impl(left, right):
-        """
+        r"""
         Called by :meth:`_composition_`.
 
         TESTS::
@@ -101,7 +101,7 @@ def _composition_impl(left, right):
         return NotImplemented
 
     def _composition_(self, other, homset):
-        """
+        r"""
         Return the composition of this elliptic-curve morphism
         with another elliptic-curve morphism.
 
@@ -133,7 +133,7 @@ def _composition_(self, other, homset):
 
     @staticmethod
     def _comparison_impl(left, right, op):
-        """
+        r"""
         Called by :meth:`_richcmp_`.
 
         TESTS::
@@ -626,7 +626,7 @@ def is_injective(self):
         return self.degree() == 1
 
     def is_zero(self):
-        """
+        r"""
         Check whether this elliptic-curve morphism is the zero map.
 
         .. NOTE::
@@ -664,7 +664,7 @@ def __neg__(self):
 
     @cached_method
     def __hash__(self):
-        """
+        r"""
         Return a hash value for this elliptic-curve morphism.
 
         ALGORITHM:
diff --git a/src/sage/schemes/elliptic_curves/hom_velusqrt.py b/src/sage/schemes/elliptic_curves/hom_velusqrt.py
index b60e551d876..4bc84d944f8 100644
--- a/src/sage/schemes/elliptic_curves/hom_velusqrt.py
+++ b/src/sage/schemes/elliptic_curves/hom_velusqrt.py
@@ -1,7 +1,8 @@
 r"""
-Îlu algorithm for elliptic-curve isogenies
+Square‑root Vélu algorithm for elliptic-curve isogenies
 
-The Îlu algorithm computes isogenies of elliptic curves in time `\tilde
+The square-root Vélu algorithm, also called the √élu algorithm,
+computes isogenies of elliptic curves in time `\tilde
 O(\sqrt\ell)` rather than naïvely `O(\ell)`, where `\ell` is the degree.
 
 The core idea is to reindex the points in the kernel subgroup in a
@@ -26,7 +27,7 @@
     10009
     sage: phi = EllipticCurveHom_velusqrt(E, K)
     sage: phi
-    Elliptic-curve isogeny (using Îlu) of degree 10009:
+    Elliptic-curve isogeny (using square-root Vélu) of degree 10009:
       From: Elliptic Curve defined by y^2 = x^3 + 5*x + 5 over Finite Field of size 6666679
       To:   Elliptic Curve defined by y^2 = x^3 + 227975*x + 3596133 over Finite Field of size 6666679
     sage: phi.codomain()
@@ -56,7 +57,7 @@
     sage: K = E(9091, 517864)
     sage: phi = EllipticCurveHom_velusqrt(E, K, model='montgomery')
     sage: phi
-    Elliptic-curve isogeny (using Îlu) of degree 2999:
+    Elliptic-curve isogeny (using square-root Vélu) of degree 2999:
       From: Elliptic Curve defined by y^2 = x^3 + x + 1 over Finite Field of size 6666679
       To:   Elliptic Curve defined by y^2 = x^3 + 1559358*x^2 + x over Finite Field of size 6666679
 
@@ -68,7 +69,7 @@
     sage: K.order()
     37
     sage: EllipticCurveHom_velusqrt(E, K)
-    Elliptic-curve isogeny (using Îlu) of degree 37:
+    Elliptic-curve isogeny (using square-root Vélu) of degree 37:
       From: Elliptic Curve defined by y^2 + x*y + 3*y = x^3 + 2*x^2 + 4*x + 5 over Finite Field of size 101
       To:   Elliptic Curve defined by y^2 = x^3 + 66*x + 86 over Finite Field of size 101
 
@@ -88,7 +89,7 @@
 
 Furthermore, the implementation is restricted to finite fields,
 since this appears to be the most relevant application for the
-Îlu algorithm::
+square-root Vélu algorithm::
 
     sage: E = EllipticCurve('26b1')
     sage: P = E(1,0)
@@ -534,7 +535,7 @@ def _point_outside_subgroup(P):
         Frobenius). Thus, once `\pi-1` can be represented in Sage,
         we may just return that in
         :meth:`~sage.schemes.elliptic_curves.ell_field.EllipticCurve_field.isogeny`
-        rather than insisting on using Îlu.
+        rather than insisting on using square-root Vélu.
     """
     E = P.curve()
     n = P.order()
@@ -555,7 +556,7 @@ def _point_outside_subgroup(P):
 class EllipticCurveHom_velusqrt(EllipticCurveHom):
     r"""
     This class implements separable odd-degree isogenies of elliptic
-    curves over finite fields using the Îlu algorithm.
+    curves over finite fields using the square-root Vélu algorithm.
 
     The complexity is `\tilde O(\sqrt{\ell})` base-field operations,
     where `\ell` is the degree.
@@ -578,7 +579,7 @@ class EllipticCurveHom_velusqrt(EllipticCurveHom):
         sage: E = EllipticCurve(F, [t,t])
         sage: K = E(2154*t^2 + 5711*t + 2899, 7340*t^2 + 4653*t + 6935)
         sage: phi = EllipticCurveHom_velusqrt(E, K); phi
-        Elliptic-curve isogeny (using Îlu) of degree 601:
+        Elliptic-curve isogeny (using square-root Vélu) of degree 601:
           From: Elliptic Curve defined by y^2 = x^3 + t*x + t over Finite Field in t of size 10009^3
           To:   Elliptic Curve defined by y^2 = x^3 + (263*t^2+3173*t+4759)*x + (3898*t^2+6111*t+9443) over Finite Field in t of size 10009^3
         sage: phi(K)
@@ -644,7 +645,7 @@ class EllipticCurveHom_velusqrt(EllipticCurveHom):
     """
     def __init__(self, E, P, *, codomain=None, model=None, Q=None):
         r"""
-        Initialize this Îlu isogeny from a kernel point of odd order.
+        Initialize this square-root Vélu isogeny from a kernel point of odd order.
 
         EXAMPLES::
 
@@ -652,7 +653,7 @@ def __init__(self, E, P, *, codomain=None, model=None, Q=None):
             sage: E = EllipticCurve(GF(71), [5,5])
             sage: P = E(-2, 22)
             sage: EllipticCurveHom_velusqrt(E, P)
-            Elliptic-curve isogeny (using Îlu) of degree 19:
+            Elliptic-curve isogeny (using square-root Vélu) of degree 19:
               From: Elliptic Curve defined by y^2 = x^3 + 5*x + 5 over Finite Field of size 71
               To:   Elliptic Curve defined by y^2 = x^3 + 13*x + 11 over Finite Field of size 71
 
@@ -661,16 +662,16 @@ def __init__(self, E, P, *, codomain=None, model=None, Q=None):
             sage: E.<P> = EllipticCurve(GF(419), [1,0])
             sage: K = 4*P
             sage: EllipticCurveHom_velusqrt(E, K)
-            Elliptic-curve isogeny (using Îlu) of degree 105:
+            Elliptic-curve isogeny (using square-root Vélu) of degree 105:
               From: Elliptic Curve defined by y^2 = x^3 + x over Finite Field of size 419
               To:   Elliptic Curve defined by y^2 = x^3 + 301*x + 86 over Finite Field of size 419
             sage: E2 = EllipticCurve(GF(419), [0,6,0,385,42])
             sage: EllipticCurveHom_velusqrt(E, K, codomain=E2)
-            Elliptic-curve isogeny (using Îlu) of degree 105:
+            Elliptic-curve isogeny (using square-root Vélu) of degree 105:
               From: Elliptic Curve defined by y^2 = x^3 + x over Finite Field of size 419
               To:   Elliptic Curve defined by y^2 = x^3 + 6*x^2 + 385*x + 42 over Finite Field of size 419
             sage: EllipticCurveHom_velusqrt(E, K, model="montgomery")
-            Elliptic-curve isogeny (using Îlu) of degree 105:
+            Elliptic-curve isogeny (using square-root Vélu) of degree 105:
               From: Elliptic Curve defined by y^2 = x^3 + x over Finite Field of size 419
               To:   Elliptic Curve defined by y^2 = x^3 + 6*x^2 + x over Finite Field of size 419
 
@@ -739,9 +740,10 @@ def __init__(self, E, P, *, codomain=None, model=None, Q=None):
 
     def _raw_eval(self, x, y=None):
         r"""
-        Evaluate the "inner" Îlu isogeny (i.e., without applying
-        pre- and post-isomorphism) at either just an `x`-coordinate
-        or a pair `(x,y)` of coordinates.
+        Evaluate the "inner" square-root Vélu isogeny
+        (i.e., without applying pre- and post-isomorphism)
+        at either just an `x`-coordinate or a pair
+        `(x,y)` of coordinates.
 
         If the given point lies in the kernel, the empty tuple
         ``()`` is returned.
@@ -812,7 +814,8 @@ def _raw_eval(self, x, y=None):
 
     def _compute_codomain(self, model=None):
         r"""
-        Helper method to compute the codomain of this Îlu isogeny
+        Helper method to compute the codomain of this
+        square-root Vélu isogeny
         once the data for :meth:`_raw_eval` has been initialized.
 
         Called by the constructor.
@@ -839,7 +842,7 @@ def _compute_codomain(self, model=None):
 
             sage: phi._compute_codomain('montgomery')
             sage: phi
-            Elliptic-curve isogeny (using Îlu) of degree 19:
+            Elliptic-curve isogeny (using square-root Vélu) of degree 19:
               From: Elliptic Curve defined by y^2 = x^3 + 5*x^2 + x over Finite Field of size 71
               To:   Elliptic Curve defined by y^2 = x^3 + 40*x^2 + x over Finite Field of size 71
 
@@ -849,7 +852,7 @@ def _compute_codomain(self, model=None):
             sage: E = EllipticCurve([3*t, 2*t+4, 3*t+2, t+4, 3*t])
             sage: K = E(3*t, 2)
             sage: EllipticCurveHom_velusqrt(E, K)   # indirect doctest
-            Elliptic-curve isogeny (using Îlu) of degree 19:
+            Elliptic-curve isogeny (using square-root Vélu) of degree 19:
               From: Elliptic Curve defined by y^2 + 3*t*x*y + (3*t+2)*y = x^3 + (2*t+4)*x^2 + (t+4)*x + 3*t over Finite Field in t of size 5^2
               To:   Elliptic Curve defined by y^2 = x^3 + (4*t+3)*x + 2 over Finite Field in t of size 5^2
         """
@@ -883,7 +886,7 @@ def _compute_codomain(self, model=None):
 
     def _eval(self, P):
         r"""
-        Evaluate this Îlu isogeny at a point.
+        Evaluate this square-root Vélu isogeny at a point.
 
         INPUT:
 
@@ -896,7 +899,7 @@ def _eval(self, P):
             sage: K = E(4, 19)
             sage: phi = EllipticCurveHom_velusqrt(E, K, model='montgomery')
             sage: phi
-            Elliptic-curve isogeny (using Îlu) of degree 19:
+            Elliptic-curve isogeny (using square-root Vélu) of degree 19:
               From: Elliptic Curve defined by y^2 = x^3 + 5*x^2 + x over Finite Field of size 71
               To:   Elliptic Curve defined by y^2 = x^3 + 40*x^2 + x over Finite Field of size 71
             sage: phi(K)
@@ -955,7 +958,7 @@ def _eval(self, P):
 
     def _repr_(self):
         r"""
-        Return basic information about this Îlu isogeny as a string.
+        Return basic information about this square-root Vélu isogeny as a string.
 
         EXAMPLES::
 
@@ -963,18 +966,18 @@ def _repr_(self):
             sage: E.<P> = EllipticCurve(GF(71), [5,5])
             sage: phi = EllipticCurveHom_velusqrt(E, P)
             sage: phi   # indirect doctest
-            Elliptic-curve isogeny (using Îlu) of degree 57:
+            Elliptic-curve isogeny (using square-root Vélu) of degree 57:
               From: Elliptic Curve defined by y^2 = x^3 + 5*x + 5 over Finite Field of size 71
               To:   Elliptic Curve defined by y^2 = x^3 + 19*x + 45 over Finite Field of size 71
         """
-        return f'Elliptic-curve isogeny (using Îlu) of degree {self._degree}:' \
+        return f'Elliptic-curve isogeny (using square-root Vélu) of degree {self._degree}:' \
                 f'\n  From: {self._domain}' \
                 f'\n  To:   {self._codomain}'
 
     @staticmethod
     def _comparison_impl(left, right, op):
         r"""
-        Compare a Îlu isogeny to another elliptic-curve morphism.
+        Compare a square-root Vélu isogeny to another elliptic-curve morphism.
 
         Called by :meth:`EllipticCurveHom._richcmp_`.
 
@@ -991,11 +994,11 @@ def _comparison_impl(left, right, op):
             sage: from sage.schemes.elliptic_curves.hom_velusqrt import EllipticCurveHom_velusqrt
             sage: E = EllipticCurve(GF(101), [5,5,5,5,5])
             sage: phi = EllipticCurveHom_velusqrt(E, E.lift_x(11)); phi
-            Elliptic-curve isogeny (using Îlu) of degree 59:
+            Elliptic-curve isogeny (using square-root Vélu) of degree 59:
               From: Elliptic Curve defined by y^2 + 5*x*y + 5*y = x^3 + 5*x^2 + 5*x + 5 over Finite Field of size 101
               To:   Elliptic Curve defined by y^2 = x^3 + 15*x + 25 over Finite Field of size 101
             sage: psi = EllipticCurveHom_velusqrt(E, E.lift_x(-1)); psi
-            Elliptic-curve isogeny (using Îlu) of degree 59:
+            Elliptic-curve isogeny (using square-root Vélu) of degree 59:
               From: Elliptic Curve defined by y^2 + 5*x*y + 5*y = x^3 + 5*x^2 + 5*x + 5 over Finite Field of size 101
               To:   Elliptic Curve defined by y^2 = x^3 + 15*x + 25 over Finite Field of size 101
             sage: phi == psi
@@ -1008,7 +1011,7 @@ def _comparison_impl(left, right, op):
     @cached_method
     def kernel_polynomial(self):
         r"""
-        Return the kernel polynomial of this Îlu isogeny.
+        Return the kernel polynomial of this square-root Vélu isogeny.
 
         .. NOTE::
 
@@ -1039,19 +1042,20 @@ def kernel_polynomial(self):
     @cached_method
     def dual(self):
         r"""
-        Return the dual of this Îlu isogeny as an :class:`EllipticCurveHom`.
+        Return the dual of this square-root Vélu
+        isogeny as an :class:`EllipticCurveHom`.
 
         .. NOTE::
 
             The dual is computed by :class:`EllipticCurveIsogeny`,
-            hence it does not benefit from the Îlu speedup.
+            hence it does not benefit from the square-root Vélu speedup.
 
         EXAMPLES::
 
             sage: E = EllipticCurve(GF(101^2), [1, 1, 1, 1, 1])
             sage: K = E.cardinality() // 11 * E.gens()[0]
             sage: phi = E.isogeny(K, algorithm='velusqrt'); phi
-            Elliptic-curve isogeny (using Îlu) of degree 11:
+            Elliptic-curve isogeny (using square-root Vélu) of degree 11:
               From: Elliptic Curve defined by y^2 + x*y + y = x^3 + x^2 + x + 1 over Finite Field in z2 of size 101^2
               To:   Elliptic Curve defined by y^2 = x^3 + 39*x + 40 over Finite Field in z2 of size 101^2
             sage: phi.dual()
@@ -1072,7 +1076,7 @@ def dual(self):
     @cached_method
     def rational_maps(self):
         r"""
-        Return the pair of explicit rational maps of this Îlu isogeny
+        Return the pair of explicit rational maps of this square-root Vélu isogeny
         as fractions of bivariate polynomials in `x` and `y`.
 
         .. NOTE::
@@ -1084,7 +1088,7 @@ def rational_maps(self):
             sage: E = EllipticCurve(GF(101^2), [1, 1, 1, 1, 1])
             sage: K = (E.cardinality() // 11) * E.gens()[0]
             sage: phi = E.isogeny(K, algorithm='velusqrt'); phi
-            Elliptic-curve isogeny (using Îlu) of degree 11:
+            Elliptic-curve isogeny (using square-root Vélu) of degree 11:
               From: Elliptic Curve defined by y^2 + x*y + y = x^3 + x^2 + x + 1 over Finite Field in z2 of size 101^2
               To:   Elliptic Curve defined by y^2 = x^3 + 39*x + 40 over Finite Field in z2 of size 101^2
             sage: phi.rational_maps()
@@ -1109,7 +1113,8 @@ def rational_maps(self):
     @cached_method
     def x_rational_map(self):
         r"""
-        Return the `x`-coordinate rational map of this Îlu isogeny
+        Return the `x`-coordinate rational map of
+        this square-root Vélu isogeny
         as a univariate rational function in `x`.
 
         .. NOTE::
@@ -1121,7 +1126,7 @@ def x_rational_map(self):
             sage: E = EllipticCurve(GF(101^2), [1, 1, 1, 1, 1])
             sage: K = (E.cardinality() // 11) * E.gens()[0]
             sage: phi = E.isogeny(K, algorithm='velusqrt'); phi
-            Elliptic-curve isogeny (using Îlu) of degree 11:
+            Elliptic-curve isogeny (using square-root Vélu) of degree 11:
               From: Elliptic Curve defined by y^2 + x*y + y = x^3 + x^2 + x + 1 over Finite Field in z2 of size 101^2
               To:   Elliptic Curve defined by y^2 = x^3 + 39*x + 40 over Finite Field in z2 of size 101^2
             sage: phi.x_rational_map()
@@ -1145,7 +1150,7 @@ def x_rational_map(self):
     def scaling_factor(self):
         r"""
         Return the Weierstrass scaling factor associated to this
-        Îlu isogeny.
+        square-root Vélu isogeny.
 
         The scaling factor is the constant `u` (in the base field)
         such that `\varphi^* \omega_2 = u \omega_1`, where
@@ -1158,7 +1163,7 @@ def scaling_factor(self):
             sage: E = EllipticCurve(GF(101^2), [1, 1, 1, 1, 1])
             sage: K = (E.cardinality() // 11) * E.gens()[0]
             sage: phi = E.isogeny(K, algorithm='velusqrt', model='montgomery'); phi
-            Elliptic-curve isogeny (using Îlu) of degree 11:
+            Elliptic-curve isogeny (using square-root Vélu) of degree 11:
               From: Elliptic Curve defined by y^2 + x*y + y = x^3 + x^2 + x + 1 over Finite Field in z2 of size 101^2
               To:   Elliptic Curve defined by y^2 = x^3 + 61*x^2 + x over Finite Field in z2 of size 101^2
             sage: phi.scaling_factor()
@@ -1187,7 +1192,7 @@ def is_separable(self):
 
 def _random_example_for_testing():
     r"""
-    Function to generate somewhat random valid Îlu inputs
+    Function to generate somewhat random valid Vélu inputs
     for testing purposes.
 
     EXAMPLES::
diff --git a/src/sage/schemes/elliptic_curves/padic_lseries.py b/src/sage/schemes/elliptic_curves/padic_lseries.py
index f264ca20434..7f32eb10a21 100644
--- a/src/sage/schemes/elliptic_curves/padic_lseries.py
+++ b/src/sage/schemes/elliptic_curves/padic_lseries.py
@@ -1375,7 +1375,7 @@ def _prec_bounds(self, n, prec):
         return [infinity] + [2 * e[j] - c0 for j in range(1, len(e))]
 
     def _poly(self, a):
-        """
+        r"""
         Given an element a in Qp[alpha] this returns the list
         containing the two coordinates in Qp.
 
@@ -1672,7 +1672,9 @@ def Dp_valued_height(self,prec=20):
         elog = Ehat.log(prec + Integer(3))
 
         # we will have to do it properly with David Harvey's _multiply_point()
-        n = LCM(E.tamagawa_numbers())
+        # import here to avoid circular import
+        from sage.schemes.elliptic_curves.padics import _multiple_to_make_good_reduction
+        n = _multiple_to_make_good_reduction(E)
         n = LCM(n, E.Np(p)) # allowed here because E has good reduction at p
 
         def height(P,check=True):
diff --git a/src/sage/schemes/elliptic_curves/padics.py b/src/sage/schemes/elliptic_curves/padics.py
index 30587d1cfa5..f59e4b823fc 100644
--- a/src/sage/schemes/elliptic_curves/padics.py
+++ b/src/sage/schemes/elliptic_curves/padics.py
@@ -2,7 +2,7 @@
 #
 # All these methods are imported in EllipticCurve_rational_field,
 # so there is no reason to add this module to the documentation.
-"""
+r"""
 Miscellaneous `p`-adic methods
 """
 
@@ -232,6 +232,8 @@ def padic_lseries(self, p, normalize=None, implementation='eclib',
 def padic_regulator(self, p, prec=20, height=None, check_hypotheses=True):
     r"""
     Compute the cyclotomic `p`-adic regulator of this curve.
+    The model of the curve needs to be integral and minimal at `p`.
+    Moreover the reduction at `p` should not be additive.
 
     INPUT:
 
@@ -246,16 +248,11 @@ def padic_regulator(self, p, prec=20, height=None, check_hypotheses=True):
     - ``check_hypotheses`` -- boolean, whether to check
       that this is a curve for which the p-adic height makes sense
 
-    OUTPUT: The p-adic cyclotomic regulator of this curve, to the
+    OUTPUT: The `p`-adic cyclotomic regulator of this curve, to the
     requested precision.
 
     If the rank is 0, we output 1.
 
-    .. TODO::
-
-        Remove restriction that curve must be in minimal
-        Weierstrass form. This is currently required for E.gens().
-
     AUTHORS:
 
     - Liang Xiao: original implementation at the 2006 MSRI
@@ -308,7 +305,7 @@ def padic_regulator(self, p, prec=20, height=None, check_hypotheses=True):
         sage: E.padic_regulator(7) == Em.padic_regulator(7)
         True
 
-    Allow a Python int as input::
+    Allow a python int as input::
 
         sage: E = EllipticCurve('37a')
         sage: E.padic_regulator(int(5))
@@ -343,9 +340,10 @@ def padic_regulator(self, p, prec=20, height=None, check_hypotheses=True):
 def padic_height_pairing_matrix(self, p, prec=20, height=None, check_hypotheses=True):
     r"""
     Computes the cyclotomic `p`-adic height pairing matrix of
-    this curve with respect to the basis self.gens() for the
-    Mordell-Weil group for a given odd prime p of good ordinary
+    this curve with respect to the basis ``self.gens()`` for the
+    Mordell-Weil group for a given odd prime `p` of good ordinary
     reduction.
+    The model needs to be integral and minimal at `p`.
 
     INPUT:
 
@@ -360,14 +358,9 @@ def padic_height_pairing_matrix(self, p, prec=20, height=None, check_hypotheses=
     - ``check_hypotheses`` -- boolean, whether to check
       that this is a curve for which the p-adic height makes sense
 
-    OUTPUT: The p-adic cyclotomic height pairing matrix of this curve
+    OUTPUT: The `p`-adic cyclotomic height pairing matrix of this curve
     to the given precision.
 
-    .. TODO::
-
-        remove restriction that curve must be in minimal
-        Weierstrass form. This is currently required for E.gens().
-
     AUTHORS:
 
     - David Harvey, Liang Xiao, Robert Bradshaw, Jennifer
@@ -428,14 +421,14 @@ def padic_height_pairing_matrix(self, p, prec=20, height=None, check_hypotheses=
 
 def _multiply_point(E, R, P, m):
     r"""
-    Computes coordinates of a multiple of P with entries in a ring.
+    Computes coordinates of a multiple of `P` with entries in a ring.
 
     INPUT:
 
-    - ``E`` -- elliptic curve over Q with integer
+    - ``E`` -- elliptic curve over `\QQ` with integer
       coefficients
 
-    - ``P`` -- a rational point on P that reduces to a
+    - ``P`` -- a rational point on `P` that reduces to a
       non-singular point at all primes
 
     - ``R`` -- a ring in which 2 is invertible (typically
@@ -455,8 +448,7 @@ def _multiply_point(E, R, P, m):
     that the sign for `b'` will match the sign for
     `d'`.
 
-    ALGORITHM: Proposition 9 of "Efficient Computation of p-adic
-    Heights" (David Harvey, to appear in LMS JCM).
+    ALGORITHM: Proposition 9 in [Har2009]_.
 
     Complexity is soft-`O(\log L \log m + \log^2 m)`.
 
@@ -584,12 +576,89 @@ def _multiply_point(E, R, P, m):
 
     return theta, omega, psi_m * d
 
+def _multiple_to_make_good_reduction(E):
+    r"""
+    Return the integer `n_2` such that for all points `P` in `E(\QQ)`
+    `n_2*P` has good reduction at all primes.
+    If the model is globally minimal the lcm of the
+    Tamagawa numbers will do, otherwise we have to take into
+    account the change of the model.
+
+    INPUT:
+
+    - ``E`` -- an elliptic curve over `\QQ`
+
+    OUTPUT:
+
+    - a positive integer ``n2``
+
+    EXAMPLE (:trac:`34790`)::
+
+        sage: from sage.schemes.elliptic_curves.padics import _multiple_to_make_good_reduction
+        sage: E = EllipticCurve([-1728,-100656])
+        sage: _multiple_to_make_good_reduction(E)
+        30
+
+    The number ``n_2`` is not always optimal but it is in this example.
+    The first multiple of the generator `P` with good reduction in this
+    non-minimal model is `30 P`.
+
+    TESTS::
+
+        sage: from sage.schemes.elliptic_curves.padics import _multiple_to_make_good_reduction
+        sage: E = EllipticCurve([1/2,1/3])
+        sage: _multiple_to_make_good_reduction(E)
+        Traceback (most recent call last):
+        ...
+        NotImplementedError: This only implemented for integral models. Please change the model first.
+        sage: E = EllipticCurve([0,3])
+        sage: _multiple_to_make_good_reduction(E)
+        1
+        sage: E = EllipticCurve([0,5^7]) # min eq is additive
+        sage: _multiple_to_make_good_reduction(E)
+        5
+        sage: E = EllipticCurve([7,0,0,0,7^7]) # min eq is split mult
+        sage: _multiple_to_make_good_reduction(E)
+        6
+        sage: E = EllipticCurve([0,-3^2,0,0,3^7]) # min eq is non-split mult
+        sage: _multiple_to_make_good_reduction(E)
+        4
+
+    """
+    if not E.is_integral():
+        st = ("This only implemented for integral models. "
+              "Please change the model first.")
+        raise NotImplementedError(st)
+    if E.is_minimal():
+        n2 = arith.LCM(E.tamagawa_numbers())
+    else:
+        # generalising to number fields one can get the u from local_data
+        Emin = E.global_minimal_model()
+        iota = E.isomorphism_to(Emin)
+        u = Integer(iota.u)
+        ps = u.prime_divisors()
+        li = []
+        for p in ps:
+            np = u.valuation(p)
+            if Emin.discriminant() %p != 0:
+                li.append(Emin.Np(p) * p**(np-1))
+            elif Emin.has_additive_reduction(p):
+                li.append(E.tamagawa_number(p) * p**np)
+            elif E.has_split_multiplicative_reduction(p):
+                li.append(E.tamagawa_number(p) * (p-1) * p**(np-1))
+            else: # non split
+                li.append(E.tamagawa_number(p) * (p+1) * p**(np-1))
+        otherbad = Integer(Emin.discriminant()).prime_divisors()
+        otherbad = [p for p in otherbad if u%p != 0 ]
+        li += [E.tamagawa_number(p) for p in otherbad]
+        n2 = arith.LCM(li)
+    return n2
 
 def padic_height(self, p, prec=20, sigma=None, check_hypotheses=True):
     r"""
-    Compute the cyclotomic p-adic height.
+    Compute the cyclotomic `p`-adic height.
 
-    The equation of the curve must be minimal at `p`.
+    The equation of the curve must be integral and minimal at `p`.
 
     INPUT:
 
@@ -606,7 +675,7 @@ def padic_height(self, p, prec=20, sigma=None, check_hypotheses=True):
 
     OUTPUT: A function that accepts two parameters:
 
-    - a Q-rational point on the curve whose height should be computed
+    - a `\QQ`-rational point on the curve whose height should be computed
 
     - optional boolean flag 'check': if False, it skips some input
       checking, and returns the p-adic height of that point to the
@@ -614,8 +683,8 @@ def padic_height(self, p, prec=20, sigma=None, check_hypotheses=True):
 
     - The normalization (sign and a factor 1/2 with respect to some other
       normalizations that appear in the literature) is chosen in such a way
-      as to make the p-adic Birch Swinnerton-Dyer conjecture hold as stated
-      in [Mazur-Tate-Teitelbaum].
+      as to make the `p`-adic Birch Swinnerton-Dyer conjecture hold as stated
+      in [MTT1986]_.
 
     AUTHORS:
 
@@ -734,11 +803,10 @@ def padic_height(self, p, prec=20, sigma=None, check_hypotheses=True):
 
     # else good ordinary case
 
-    # For notation and definitions, see "Efficient Computation of
-    # p-adic Heights", David Harvey (unpublished)
+    # For notation and definitions, see [Har2009]_.
 
     n1 = self.change_ring(rings.GF(p)).cardinality()
-    n2 = arith.LCM(self.tamagawa_numbers())
+    n2 = _multiple_to_make_good_reduction(self)
     n = arith.LCM(n1, n2)
     m = int(n / n2)
 
@@ -792,7 +860,7 @@ def height(P, check=True):
 
 def padic_height_via_multiply(self, p, prec=20, E2=None, check_hypotheses=True):
     r"""
-    Computes the cyclotomic p-adic height.
+    Computes the cyclotomic `p`-adic height.
 
     The equation of the curve must be minimal at `p`.
 
@@ -805,7 +873,7 @@ def padic_height_via_multiply(self, p, prec=20, E2=None, check_hypotheses=True):
 
     - ``E2`` -- precomputed value of E2. If not supplied,
       this function will call padic_E2 to compute it. The value supplied
-      must be correct mod `p^(prec-2)` (or slightly higher in the
+      must be correct mod `p^{prec-2}` (or slightly higher in the
       anomalous case; see the code for details).
 
     - ``check_hypotheses`` -- boolean, whether to check
@@ -813,22 +881,21 @@ def padic_height_via_multiply(self, p, prec=20, E2=None, check_hypotheses=True):
 
     OUTPUT: A function that accepts two parameters:
 
-    - a Q-rational point on the curve whose height should be computed
+    - a `\QQ`-rational point on the curve whose height should be computed
 
     - optional boolean flag 'check': if False, it skips some input
-      checking, and returns the p-adic height of that point to the
+      checking, and returns the `p`-adic height of that point to the
       desired precision.
 
     - The normalization (sign and a factor 1/2 with respect to some other
       normalizations that appear in the literature) is chosen in such a way
       as to make the p-adic Birch Swinnerton-Dyer conjecture hold as stated
-      in [Mazur-Tate-Teitelbaum].
+      in [MTT1986]_.
 
     AUTHORS:
 
     - David Harvey (2008-01): based on the padic_height() function,
-      using the algorithm of"Computing p-adic heights via
-      point multiplication"
+      using the algorithm of [Har2009]_.
 
     EXAMPLES::
 
@@ -879,11 +946,10 @@ def padic_height_via_multiply(self, p, prec=20, E2=None, check_hypotheses=True):
     if prec < 1:
         raise ValueError("prec (=%s) must be at least 1" % prec)
 
-    # For notation and definitions, see ``Computing p-adic heights via point
-    # multiplication'' (David Harvey, still in draft form)
+    # For notation and definitions, [Har2009]_
 
     n1 = self.change_ring(rings.GF(p)).cardinality()
-    n2 = arith.LCM(self.tamagawa_numbers())
+    n2 = _multiple_to_make_good_reduction(self)
     n = arith.LCM(n1, n2)
     m = int(n / n2)
 
@@ -938,9 +1004,9 @@ def height(P, check=True):
 
 def padic_sigma(self, p, N=20, E2=None, check=False, check_hypotheses=True):
     r"""
-    Computes the p-adic sigma function with respect to the standard
+    Computes the `p`-adic sigma function with respect to the standard
     invariant differential `dx/(2y + a_1 x + a_3)`, as
-    defined by Mazur and Tate, as a power series in the usual
+    defined by Mazur and Tate in [MT1991]_, as a power series in the usual
     uniformiser `t` at the origin.
 
     The equation of the curve must be minimal at `p`.
@@ -963,7 +1029,7 @@ def padic_sigma(self, p, N=20, E2=None, check=False, check_hypotheses=True):
 
     - ``differential equation`` -- note that this does NOT
       guarantee correctness of all the returned digits, but it comes
-      pretty close :-))
+      pretty close.
 
     - ``check_hypotheses`` -- boolean, whether to check
       that this is a curve for which the p-adic sigma function makes
@@ -984,9 +1050,9 @@ def padic_sigma(self, p, N=20, E2=None, check=False, check_hypotheses=True):
     `p`-adic digits).
 
     ALGORITHM: Described in "Efficient Computation of p-adic Heights"
-    (David Harvey), which is basically an optimised version of the
+    (David Harvey) [Har2009]_ which is basically an optimised version of the
     algorithm from "p-adic Heights and Log Convergence" (Mazur, Stein,
-    Tate).
+    Tate) [MST2006]_.
 
     Running time is soft-`O(N^2 \log p)`, plus whatever time is
     necessary to compute `E_2`.
@@ -1159,7 +1225,7 @@ def padic_sigma_truncated(self, p, N=20, lamb=0, E2=None, check_hypotheses=True)
     r"""
     Compute the p-adic sigma function with respect to the standard
     invariant differential `dx/(2y + a_1 x + a_3)`, as
-    defined by Mazur and Tate, as a power series in the usual
+    defined by Mazur and Tate in [MT1991]_, as a power series in the usual
     uniformiser `t` at the origin.
 
     The equation of the curve must be minimal at `p`.
@@ -1193,10 +1259,9 @@ def padic_sigma_truncated(self, p, N=20, lamb=0, E2=None, check_hypotheses=True)
     correct to precision `O(p^{N - 2 + (3 - j)(lamb + 1)})`.
 
     ALGORITHM: Described in "Efficient Computation of p-adic Heights"
-    (David Harvey, to appear in LMS JCM), which is basically an
+    [Har2009]_, which is basically an
     optimised version of the algorithm from "p-adic Heights and Log
-    Convergence" (Mazur, Stein, Tate), and "Computing p-adic heights
-    via point multiplication" (David Harvey, still draft form).
+    Convergence" (Mazur, Stein, Tate) [MST2006]_.
 
     Running time is soft-`O(N^2 \lambda^{-1} \log p)`, plus
     whatever time is necessary to compute `E_2`.
@@ -1378,7 +1443,7 @@ def padic_E2(self, p, prec=20, check=False, check_hypotheses=True, algorithm="au
     trick.
 
     EXAMPLES: Here is the example discussed in the paper "Computation
-    of p-adic Heights and Log Convergence" (Mazur, Stein, Tate)::
+    of p-adic Heights and Log Convergence" (Mazur, Stein, Tate) [MST2006]_::
 
         sage: EllipticCurve([-1, 1/4]).padic_E2(5)
         2 + 4*5 + 2*5^3 + 5^4 + 3*5^5 + 2*5^6 + 5^8 + 3*5^9 + 4*5^10 + 2*5^11 + 2*5^12 + 2*5^14 + 3*5^15 + 3*5^16 + 3*5^17 + 4*5^18 + 2*5^19 + O(5^20)
@@ -1673,7 +1738,7 @@ def _brent(F, p, N):
     some log-log factors.
 
     For more information, and a proof of the precision guarantees, see
-    Lemma 4 in "Efficient Computation of p-adic Heights" (David Harvey).
+    Lemma 4 in [Har2009]_.
 
     AUTHORS: