diff --git a/darmonpoints/Fdoms/fdom.gp b/darmonpoints/Fdoms/fdom.gp
index 65f461c..ce92180 100644
--- a/darmonpoints/Fdoms/fdom.gp
+++ b/darmonpoints/Fdoms/fdom.gp
@@ -115,7 +115,7 @@ addhelp(fdom, "This package can be used to compute fundamental domains for Shimu
 	install("enum_successrate","GGGLD0,G,p","enum_successrate","libfdom.so");
 	addhelp(enum_successrate,"Inputs A, p, C, Ntests, {R=0}: quaternion algebra A corresponding to a Shimura curve, upper half plane point p, positive real number C, positive integer Ntests, positive real R.\n Computes the small norm 1 elements of A (<=C) Ntests times, where we pick z_1=0 and z_2 a random point in the hyperbolic disc of radius R. If R=0, we auto-set R to be the same R as the algfdom method. We output the pair [obs, exp], of the number of found norm 1 elements, and the expected number.");
 	install("enum_successrate_range","GGGGLLD0,G,DsD1,L,D1,L,p","enum_successrate_range","libfdom.so");
-	addhelp(enum_successrate_range,"Inputs A, p, Cmin, Cmax, ntrials, Ntests, {R=0}, {fname=NULL}, {compile=1}, {WSL=1}: q-alg A corresp. to a Shimura curve, upper half plane point p, 0<Cmin<Cmax, ntrials>1 and Ntests>0 integers, R>=0, fname=file name, compile and WSL=0, 1.\n Runs enum_successrate on ntrials trials of C between Cmin and Cmax. This prints the results to the file pltos/build/fname.dat, and retuns [A, B, R^2], where the expected trend line is A+B*C, and R^2 is the R^2 value of this trendline with the data. If compile=1, we create and compile a LaTeX (pgfplots) of the curve. If WSL=1, we also display said plot, assuming we are using Windows Subsystem for Linux.");
+	addhelp(enum_successrate_range,"Inputs A, p, Cmin, Cmax, ntrials, Ntests, {R=0}, {fname=NULL}, {compile=1}, {WSL=1}: q-alg A corresp. to a Shimura curve, upper half plane point p, 0<Cmin<Cmax, ntrials>1 and Ntests>0 integers, R>=0, fname=file name, compile and WSL=0, 1.\n Runs enum_successrate on ntrials trials of C between Cmin and Cmax. This prints the results to the file pltos/build/fname.dat, and returns [A, B, R^2], where the expected trend line is A+B*C, and R^2 is the R^2 value of this trendline with the data. If compile=1, we create and compile a LaTeX (pgfplots) of the curve. If WSL=1, we also display said plot, assuming we are using Windows Subsystem for Linux.");
 	install("enum_time","GGGD300,L,p","enum_time","libfdom.so");
 	addhelp(enum_time,"Inputs A, p, Cset, {mintesttime=300}: quaternion algebra A corresponding to a Shimura curve, upper half plane point p, vector of positive real numbers, mintesttime positive integer.\n This computes how long the call to algsmallnorm1elts(A, p, C, z1, z2) takes for all C in Cset, and returns a column vector of the timings. This does NOT take into account time spent initializing things related to the algebra (e.g. cholesky of the norm form), since this can be computed once and reused many times. If the time taken is <mintesttime (in milliseconds), we repeat the test K times until we have taken at least mintesttime, and divide the final result by K. A larger value of mintesttime will produce more accurate results, but will take longer.");
 	install("enum_time_range","GGGGLD300,L,DsD1,L,D1,L,p","enum_time_range","libfdom.so");
@@ -131,7 +131,7 @@ addhelp(fdom, "This package can be used to compute fundamental domains for Shimu
 
 	\\REGRESSIONS & PLOTS
 	install("OLS","GGD1,L,","OLS","libfdom.so");
-	addhelp(OLS,"Inputs X, y, {retrsqr=1}:  m*n matrix X with top row being all 1's, length n column vector y, retrsqr=0, 1.\n Performs ordinary least squares regression on the data, where the n inputs are the columns of X, and the outputs are the entries of y. We must include a constant term, hence why the first row of X must be all 1's. If retrsqr=1, returns [pararms, R^2], and otherwise returns params, where params is the length m column vector of best fit parameters.");
+	addhelp(OLS,"Inputs X, y, {retrsqr=1}:  m*n matrix X with top row being all 1's, length n column vector y, retrsqr=0, 1.\n Performs ordinary least squares regression on the data, where the n inputs are the columns of X, and the outputs are the entries of y. We must include a constant term, hence why the first row of X must be all 1's. If retrsqr=1, returns [params, R^2], and otherwise returns params, where params is the length m column vector of best fit parameters.");
 	install("OLS_nointercept","GGD1,L,","OLS_nointercept","libfdom.so");
 	addhelp(OLS_nointercept,"Inputs X, y, {retrsqr=1}: vector X, column vector y (of same length), retrsqr=0, 1.\n Performs ordinary least squares regression on the data assuming that y[i]=c*X[i], i.e. the y-intercept is 0. Returns c if retrsqr=0, or [c, R^2] otherwise.");
 	install("OLS_single","GGD1,L,","OLS_single","libfdom.so");
diff --git a/darmonpoints/arithgroup.py b/darmonpoints/arithgroup.py
index 19cf802..d9efbe3 100644
--- a/darmonpoints/arithgroup.py
+++ b/darmonpoints/arithgroup.py
@@ -121,15 +121,15 @@ def intersect_geodesic_arcs(x1, x2, y1, y2):
     r"""
     TESTS::
 
-    sage: from darmonpoints.arithgroup import intersect_geodesic_arcs
-    sage: intersect_geodesic_arcs(1,3,2,4)
-    1/2*I*sqrt(3) + 5/2
-    sage: print(intersect_geodesic_arcs(-1, 1, 0, AA(-1).sqrt()))
-    None
-    sage: intersect_geodesic_arcs(-1, 1, 0, 2*AA(-1).sqrt())
-    I
-    sage: intersect_geodesic_arcs(-3, 3, 2*AA(-1).sqrt(), Infinity)
-    3*I
+        sage: from darmonpoints.arithgroup import intersect_geodesic_arcs
+        sage: intersect_geodesic_arcs(1,3,2,4)
+        1/2*I*sqrt(3) + 5/2
+        sage: print(intersect_geodesic_arcs(-1, 1, 0, AA(-1).sqrt()))
+        None
+        sage: intersect_geodesic_arcs(-1, 1, 0, 2*AA(-1).sqrt())
+        I
+        sage: intersect_geodesic_arcs(-3, 3, 2*AA(-1).sqrt(), Infinity)
+        3*I
     """
     # verbose('Entering intersect_geodesic_arcs')
     e1 = geodesic_circle(x1, x2)
@@ -1810,7 +1810,7 @@ def compute_quadratic_embedding(self, D, return_generator=False, **kwargs):
         try:
             _, iota = self.magma.Embed(OK_magma, O_magma, nvals=2)
         except RuntimeError:
-            print("An error ocurred!")
+            print("An error occurred!")
             print("OK_magma = %s" % OK_magma)
             print("O_magma =" % O_magma)
             raise RuntimeError("Error while computing quadratic embedding")
@@ -2538,7 +2538,7 @@ def __init__(
         }
         verbose("Done computing coset reps.")
         ## compute the generators of H
-        verbose("Computing the auxiliar data...")
+        verbose("Computing the auxiliary data...")
         (
             self._gens_dict_auxiliary,
             self._gens_matrices_auxiliary,
@@ -2744,14 +2744,14 @@ def _get_word_rep_auxiliary(self, h, check=False):
 
         Firstly, we write this as h = 1.h. Then we write h = gh', where g in Gens(G) (so we must be
         able to solve the word problem for G). Then write 1.g = zp', so that
-            h = z * p' * h'. Now iterate. We will end up with z_1 z_2 ... z_t p_0, where p_0 = id rep.
+        h = z * p' * h'. Now iterate. We will end up with z_1 z_2 ... z_t p_0, where p_0 = id rep.
 
         OUTPUT:
 
         - a list of integers in {-t,-t+1,...,t-1,t}, where the output of _generators_auxiliary is
           [a_1,...,a_t].
 
-        For example,
+        For example::
 
            h = abc in H, a,b,c in Gens(G)
            h = 1.abc
diff --git a/darmonpoints/arithgroup_element.py b/darmonpoints/arithgroup_element.py
index d97ba5b..d58136a 100644
--- a/darmonpoints/arithgroup_element.py
+++ b/darmonpoints/arithgroup_element.py
@@ -46,7 +46,7 @@ def __init__(self, parent, word_rep=None, quaternion_rep=None, check=False):
         INPUT:
 
         - a list of the form [(g1,a1),(g2,a2),...,(gn,an)] where the gi are indices
-          refering to fixed generators, and the ai are integers, or
+          referring to fixed generators, and the ai are integers, or
           an element of the quaternion algebra ``self.parent().quaternion_algebra()``.
 
         """
@@ -281,7 +281,7 @@ def find_bounding_cycle(self, G, npow=1):
             # If a < 0 use the relation g^a = -g^(-a) + del(g^a|g^(-a))
             ans.extend([(npow, [gaq**-1], ga)] if a < 0 else [])
             # By the above line we have to deal with g^a with -g^abs(a) if a <0
-            # We add the corresponding boundaries, which we will substract if a > 0 and add if a < 0
+            # We add the corresponding boundaries, which we will subtract if a > 0 and add if a < 0
             ans.extend(
                 [(-sgn(a) * npow, [gq**j for j in range(1, abs(a))], g)]
                 if abs(a) > 1
diff --git a/darmonpoints/arithgroup_generic.py b/darmonpoints/arithgroup_generic.py
index 1d7efbe..4816f1e 100644
--- a/darmonpoints/arithgroup_generic.py
+++ b/darmonpoints/arithgroup_generic.py
@@ -526,7 +526,7 @@ def decompose_into_commutators(self, x):
         commutator_list = []
         for i in range(len(self.gens())):
             while True:
-                # Find the first occurence of generator i
+                # Find the first occurrence of generator i
                 try:
                     idx = [x[0] for x in oldword[1:]].index(i) + 1
                 except ValueError:
@@ -638,7 +638,7 @@ def compute_cusp_stabiliser(self, cusp_matrix):
         We know the stabiliser of infinity is given by matrices of form
         (u, a; 0, u^-1), so a finite index subgroup is generated by (1, alpha; 0, 1)
         and (1, 1; 0, 1) for K = Q(alpha). Given the cusp, we use a matrix
-        sending infinty to that cusp, and the conjugate by it, before taking powers
+        sending infinity to that cusp, and the conjugate by it, before taking powers
         to ensure the result is integral and lies in Gamma_0(N).
 
         Input:
@@ -662,7 +662,7 @@ def compute_cusp_stabiliser(self, cusp_matrix):
             infinity_gens = [matrix([[1, 1], [0, 1]])]
             N_ideal = ZZ.ideal(P.N())
 
-        ## Initilise (empty) list of generators of Stab_Gamma(cusp)
+        ## Initialise (empty) list of generators of Stab_Gamma(cusp)
         cusp_gens = []
 
         ## Loop over all the generators of stab at infinity, conjugate into stab at cusp
@@ -689,12 +689,13 @@ def cusp_reduction_table(self):
 
         Takes as input the object representing P^1(O_F/N), where F is a number field
         (that is possibly Q), and N is some ideal in the field.  Runs the following algorithm:
-                - take a remaining element C = (c:d) of P^1(O_F/N);
-                - add this to the set of cusps, declaring it to be our chosen rep;
-                - run through every translate C' = (c':d') of C under the stabiliser of infinity, and
-                        remove this translate from the set of remaining elements;
-                - store the matrix T in the stabiliser such that C' * T = C (as elements in P^1)
-                        in the dictionary, with key C'.
+
+        - take a remaining element C = (c:d) of P^1(O_F/N);
+        - add this to the set of cusps, declaring it to be our chosen rep;
+        - run through every translate C' = (c':d') of C under the stabiliser of infinity, and
+          remove this translate from the set of remaining elements;
+        - store the matrix T in the stabiliser such that C' * T = C (as elements in P^1)
+          in the dictionary, with key C'.
         """
         P = self.get_P1List()
         if hasattr(P.N(), "number_field"):
diff --git a/darmonpoints/cohomology_arithmetic.py b/darmonpoints/cohomology_arithmetic.py
index 26d4d6c..30e738a 100644
--- a/darmonpoints/cohomology_arithmetic.py
+++ b/darmonpoints/cohomology_arithmetic.py
@@ -468,18 +468,19 @@ def BI(self, h, j=None):  # Returns \int_{h \Z_p} z^j \Phi\{0 \to \infy\}
         Input a 2x2 matrix h in SL2(OK) (which embeds as an element of Sigma_0(p)), and a value j.
 
         Options for j:
-            - classical case: specify a non-negative integer j. Then returns the value
-                    BI_{h,j} := Int_{h.Zp} z^j . d Phi{0 --> infty},
-                that is, the value of the distribution Phi{0 --> infty} at the function z^j x
-                the indicator function of the open set h.Zp.
-
-            - Bianchi case: specify a tuple (k,l). Then returns the value
-                    BI_{h,j} := Int_{h.(Zp x Zp)} x^k y^l . d Phi{0 --> infty},
-                that is, the value of the distribution Phi{0 --> infty} at the function x^k y^l x
-                the indicator function of the open set h.(Zp x Zp).
-
-            - do not specify j. Then returns the the distribution mu whose moments are
-                BI_{h,j}.
+
+        - classical case: specify a non-negative integer j. Then returns the value
+                BI_{h,j} := Int_{h.Zp} z^j . d Phi{0 --> infty},
+            that is, the value of the distribution Phi{0 --> infty} at the function z^j x
+            the indicator function of the open set h.Zp.
+
+        - Bianchi case: specify a tuple (k,l). Then returns the value
+                BI_{h,j} := Int_{h.(Zp x Zp)} x^k y^l . d Phi{0 --> infty},
+            that is, the value of the distribution Phi{0 --> infty} at the function x^k y^l x
+            the indicator function of the open set h.(Zp x Zp).
+
+        - do not specify j. Then returns the the distribution mu whose moments are
+            BI_{h,j}.
 
         """
         V = self.parent().coefficient_module()  ## Module V in H^1(G,V)
diff --git a/darmonpoints/darmonpoints.py b/darmonpoints/darmonpoints.py
index e90366f..e406598 100644
--- a/darmonpoints/darmonpoints.py
+++ b/darmonpoints/darmonpoints.py
@@ -178,28 +178,28 @@ def darmon_point(
 
     We first need to import the module::
 
-    sage: from darmonpoints.darmonpoints import darmon_point
+        sage: from darmonpoints.darmonpoints import darmon_point
 
     A first example (Stark--Heegner point)::
 
-    sage: from darmonpoints.darmonpoints import darmon_point
-    sage: darmon_point(7,EllipticCurve('35a1'),41,20, cohomological=False, use_magma=False, use_ps_dists = True)[0]
-    Starting computation of the Darmon point
-    ...
-    -70*alpha + 449
+        sage: from darmonpoints.darmonpoints import darmon_point
+        sage: darmon_point(7,EllipticCurve('35a1'),41,20, cohomological=False, use_magma=False, use_ps_dists = True)[0]
+        Starting computation of the Darmon point
+        ...
+        -70*alpha + 449
 
     A quaternionic (Greenberg) point::
 
-    sage: darmon_point(13,EllipticCurve('78a1'),5,20) # long time # optional - magma
+        sage: darmon_point(13,EllipticCurve('78a1'),5,20) # long time # optional - magma
 
     A Darmon point over a cubic (1,1) field::
 
-    sage: F.<r> = NumberField(x^3 - x^2 - x + 2)
-    sage: E = EllipticCurve([-r -1, -r, -r - 1,-r - 1, 0])
-    sage: N = E.conductor()
-    sage: P = F.ideal(r^2 - 2*r - 1)
-    sage: beta = -3*r^2 + 9*r - 6
-    sage: darmon_point(P,E,beta,20) # long time # optional - magma
+        sage: F.<r> = NumberField(x^3 - x^2 - x + 2)
+        sage: E = EllipticCurve([-r -1, -r, -r - 1,-r - 1, 0])
+        sage: N = E.conductor()
+        sage: P = F.ideal(r^2 - 2*r - 1)
+        sage: beta = -3*r^2 + 9*r - 6
+        sage: darmon_point(P,E,beta,20) # long time # optional - magma
 
     """
     # global G, Coh, phiE, Phi, dK, J, J1, cycleGn, nn, Jlist
diff --git a/darmonpoints/darmonvonk.py b/darmonpoints/darmonvonk.py
index 6a66f50..34ef65e 100644
--- a/darmonpoints/darmonvonk.py
+++ b/darmonpoints/darmonvonk.py
@@ -292,11 +292,11 @@ def darmon_vonk_point(
     r"""
     EXAMPLES ::
 
-    sage: from darmonpoints.darmonvonk import darmon_vonk_point
-    sage: J = darmon_vonk_point(5, 1, 3, 13, 60, parity='+', recognize_point='algdep',magma=magma) # optional - magma
-    #### Starting computation of the Darmon-Vonk point ####
-    ...
-    f = 7*x^2 + 11*x + 7
+        sage: from darmonpoints.darmonvonk import darmon_vonk_point
+        sage: J = darmon_vonk_point(5, 1, 3, 13, 60, parity='+', recognize_point='algdep',magma=magma) # optional - magma
+        #### Starting computation of the Darmon-Vonk point ####
+        ...
+        f = 7*x^2 + 11*x + 7
     """
     if magma is None:
         from sage.interfaces.magma import Magma
diff --git a/darmonpoints/findcurve.py b/darmonpoints/findcurve.py
index 05beb1d..f0ef7f9 100644
--- a/darmonpoints/findcurve.py
+++ b/darmonpoints/findcurve.py
@@ -41,29 +41,29 @@ def find_curve(
 
     First example::
 
-    sage: from darmonpoints.findcurve import find_curve
-    sage: find_curve(5,6,30,20) # long time # optional - magma
-    # B = F<i,j,k>, with i^2 = -1 and j^2 = 3
-    ...
-    '(1, 0, 1, -289, 1862)'
+        sage: from darmonpoints.findcurve import find_curve
+        sage: find_curve(5,6,30,20) # long time # optional - magma
+        # B = F<i,j,k>, with i^2 = -1 and j^2 = 3
+        ...
+        '(1, 0, 1, -289, 1862)'
 
     A second example, now over a real quadratic::
 
-    sage: from darmonpoints.findcurve import find_curve
-    sage: F.<r> = QuadraticField(5)
-    sage: P = F.ideal(3/2*r + 1/2)
-    sage: D = F.ideal(3)
-    sage: find_curve(P,D,P*D,30,ramification_at_infinity = F.real_places()[:1]) # long time # optional - magma
-    ...
+        sage: from darmonpoints.findcurve import find_curve
+        sage: F.<r> = QuadraticField(5)
+        sage: P = F.ideal(3/2*r + 1/2)
+        sage: D = F.ideal(3)
+        sage: find_curve(P,D,P*D,30,ramification_at_infinity = F.real_places()[:1]) # long time # optional - magma
+        ...
 
     Now over a cubic of mixed signature::
 
-    sage: from darmonpoints.findcurve import find_curve
-    sage: F.<r> = NumberField(x^3 -3)
-    sage: P = F.ideal(r-2)
-    sage: D = F.ideal(r-1)
-    sage: find_curve(P,D,P*D,30) # long time # optional - magma
-    ...
+        sage: from darmonpoints.findcurve import find_curve
+        sage: F.<r> = NumberField(x^3 -3)
+        sage: P = F.ideal(r-2)
+        sage: D = F.ideal(r-1)
+        sage: find_curve(P,D,P*D,30) # long time # optional - magma
+        ...
 
     """
     config = configparser.ConfigParser()
diff --git a/darmonpoints/integrals.py b/darmonpoints/integrals.py
index f365ff2..9d5a91e 100644
--- a/darmonpoints/integrals.py
+++ b/darmonpoints/integrals.py
@@ -436,11 +436,11 @@ def riemann_sum(G, phi, hc, depth=1, mult=False, progress_bar=False, K=None):
         if hce == 0:
             continue
         # verbose('hc = %s'%hce)
-        te = sample_point(G, e, prec)
-        if te == Infinity:
+        t_e = sample_point(G, e, prec)
+        if t_e == Infinity:
             continue
         if mult:
-            res *= phi(K(te)) ** hce
+            res *= phi(K(t_e)) ** hce
         else:
-            res += phi(K(te)) * hce
+            res += phi(K(t_e)) * hce
     return res
diff --git a/darmonpoints/limits.py b/darmonpoints/limits.py
index 7a19599..26eb1b7 100644
--- a/darmonpoints/limits.py
+++ b/darmonpoints/limits.py
@@ -19,7 +19,7 @@
 
 
 ##########################################################################################
-# now the new functions that we need...they follow pretty close the article we're writting
+# now the new functions that we need...they follow pretty close the article we're writing
 ##########################################################################################
 def factorize_matrix(m, M):
     # assert is_in_Gamma_1(m,M,determinant_condition = False)
diff --git a/darmonpoints/my_p1list_nf.py b/darmonpoints/my_p1list_nf.py
index cf754c5..35ee6f5 100644
--- a/darmonpoints/my_p1list_nf.py
+++ b/darmonpoints/my_p1list_nf.py
@@ -190,7 +190,7 @@ def list(self):
 
     def normalize(self, c, d, with_scalar=False):
         r"""
-        [We have gutted this function from its orginal state!]
+        [We have gutted this function from its original state!]
         [Note: there is now no check!]
 
         Return a normalized  element of (a canonical representative of an element
diff --git a/darmonpoints/ocbianchi.py b/darmonpoints/ocbianchi.py
index 5772ad8..82ada47 100644
--- a/darmonpoints/ocbianchi.py
+++ b/darmonpoints/ocbianchi.py
@@ -94,7 +94,7 @@ def _an_element_(self):
 
     def _element_constructor_(self, first, second):
         r"""
-        Element constructor. Takes as input tuple (g,h) of elments
+        Element constructor. Takes as input tuple (g,h) of elements
         of Sigma_0(p) and represents the element (g,h) in the product.
         """
         return Sigma0SquaredElement(self, (self._Sigma0(first), self._Sigma0(second)))
@@ -528,7 +528,7 @@ def evaluate_at_poly(self, P, R=None, depth=None):
         ## For each monomial x^iy^j in the polynomial, multip] ly the coefficient of X^iY^j (in mu) by the
         ## coefficient of x^iy^j (in f) and take the sum. This is our final value
         ##          --> P.coefficients is a dictionary which has monomials as keys; we generate monomials using exponents.
-        ##         --> self._moments takes as input an index and spits out the cofficient. So generate the index from the exponent.
+        ##         --> self._moments takes as input an index and spits out the coefficient. So generate the index from the exponent.
         coefficient_list = []
         for polx in P.padded_list(self._depth):
             coefficient_list.extend(polx.padded_list(self._depth))
@@ -624,7 +624,7 @@ class BianchiDistributions(Module, UniqueRepresentation):
 
     def __init__(self, p, depth, act_on_left=False, adjuster=None):
         self._dimension = (
-            0  ## Hack!! Dimension was being called before it was intialised
+            0  ## Hack!! Dimension was being called before it was initialised
         )
         self._Rmod = ZpCA(p, depth - 1)  ## create Zp
         Module.__init__(self, base=self._Rmod)
diff --git a/darmonpoints/padicperiods.py b/darmonpoints/padicperiods.py
index c18420c..b6c4b4d 100644
--- a/darmonpoints/padicperiods.py
+++ b/darmonpoints/padicperiods.py
@@ -981,7 +981,7 @@ def find_igusa_invariants(
     parallelize=True,
 ):
     fwrite("# Trying to recognize invariants...", outfile)
-    Pring = embedding.domain()
+    P_ring = embedding.domain()
     if prec is None:
         prec = a.parent().precision_cap()
     Tlist = []
@@ -1001,7 +1001,7 @@ def find_igusa_invariants(
         fwrite("# Doing matrix %s / %s ( = %s)" % (ii, len(Tlist), tt.list()), outfile)
         Lp = a + b * tt
         inp_vec = [
-            (Lp, ordmat, prec, Pring, cheatjs, embedding, 3, list_I10, Pgen, outfile)
+            (Lp, ordmat, prec, P_ring, cheatjs, embedding, 3, list_I10, Pgen, outfile)
             for ordmat in all_possible_ordmats(Lp, 20)
         ]
 
@@ -1279,7 +1279,7 @@ def guess_equation(
         F = QQ
         P = Pgen
         Pnrm = Pgen
-        Pring = QQ
+        P_ring = QQ
         D = Dgen
         Np = Npgen
         Sinv_places = []
@@ -1295,7 +1295,7 @@ def guess_equation(
         r = F.gen()
         P = F.ideal(Pgen)
         Pnrm = P.norm()
-        Pring = P.ring()
+        P_ring = P.ring()
         D = F.ideal(Dgen)
         Np = F.ideal(Npgen)
         if Sinf is None:
@@ -1798,13 +1798,16 @@ def find_kadziela_matrices(M, T):
     """
     The matrix M describes the relation between periods (A,B,D)^t
     and the periods (A0,B0)^t, where (A,B,D) are the periods of
-    the Teitelbaum periods, and (A0,B0) are the Darmon ones.
+    the Teitelbaum periods, and (A0,B0) are the Darmon ones. ::
+
            (A,B,D)^t = M * (A0,B0)^t
+
     The matrix T describes the action of Hecke on homology.
     That is, the first column of T describes the image of T
     on the first basis vector.
 
-    The output are matrices X and Y such that
+    The output are matrices X and Y such that::
+
          X * matrix(2,2,[A,B,B,D]) = matrix(2,2,[A0,B0,C0,D0]) * Y
 
     """
diff --git a/darmonpoints/schottky.py b/darmonpoints/schottky.py
index bc75bd1..58c077c 100644
--- a/darmonpoints/schottky.py
+++ b/darmonpoints/schottky.py
@@ -735,23 +735,22 @@ def theta(self, prec, a, b=None, **kwargs):
 
         EXAMPLES ::
 
-        sage: from darmonpoints.schottky import *
-        sage: p = 3
-        sage: prec = 20
-        sage: working_prec = 200
-        sage: K = Qp(p,working_prec)
-        sage: g1 = matrix(K, 2, 2, [-5,32,-8,35])
-        sage: g2 = matrix(K, 2, 2, [-13,80,-8,43])
-        sage: G = SchottkyGroup(K, (g1, g2))
-        sage: a = 23
-        sage: b = 14
-        sage: z0 = K(8)
-        sage: m = 10
-        sage: Tg = G.theta_naive(m, z=z0, a=a,b=b)
-        sage: T = G.theta(prec, a, b).improve(m)
-        sage: (T(z0) / Tg - 1).valuation() > m
-        True
-
+            sage: from darmonpoints.schottky import *
+            sage: p = 3
+            sage: prec = 20
+            sage: working_prec = 200
+            sage: K = Qp(p,working_prec)
+            sage: g1 = matrix(K, 2, 2, [-5,32,-8,35])
+            sage: g2 = matrix(K, 2, 2, [-13,80,-8,43])
+            sage: G = SchottkyGroup(K, (g1, g2))
+            sage: a = 23
+            sage: b = 14
+            sage: z0 = K(8)
+            sage: m = 10
+            sage: Tg = G.theta_naive(m, z=z0, a=a,b=b)
+            sage: T = G.theta(prec, a, b).improve(m)
+            sage: (T(z0) / Tg - 1).valuation() > m
+            True
         """
         K = self.base_ring()
         if b is not None:
@@ -841,26 +840,26 @@ def period(self, i, j, prec, **kwargs):
 
         EXAMPLES ::
 
-        sage: from darmonpoints.schottky import *
-        sage: p = 3
-        sage: prec = 10
-        sage: working_prec = 200
-        sage: K = Qp(p,working_prec)
-        sage: h1 = matrix(K, 2, 2, [-5,32,-8,35])
-        sage: h2 = matrix(K, 2, 2, [-13,80,-8,43])
-        sage: G = SchottkyGroup(K, (h1,h2))
-        sage: q00g = G.period_naive(0, 0, prec)
-        sage: q01g = G.period_naive(0, 1, prec)
-        sage: q11g = G.period_naive(1, 1, prec)
-        sage: q00 = G.period(0,0, prec)
-        sage: q01 = G.period(0,1, prec)
-        sage: q11 = G.period(1,1, prec)
-        sage: (q00g-q00).valuation() >= prec
-        True
-        sage: (q01g-q01).valuation() >= prec
-        True
-        sage: (q11g-q11).valuation() >= prec
-        True
+            sage: from darmonpoints.schottky import *
+            sage: p = 3
+            sage: prec = 10
+            sage: working_prec = 200
+            sage: K = Qp(p,working_prec)
+            sage: h1 = matrix(K, 2, 2, [-5,32,-8,35])
+            sage: h2 = matrix(K, 2, 2, [-13,80,-8,43])
+            sage: G = SchottkyGroup(K, (h1,h2))
+            sage: q00g = G.period_naive(0, 0, prec)
+            sage: q01g = G.period_naive(0, 1, prec)
+            sage: q11g = G.period_naive(1, 1, prec)
+            sage: q00 = G.period(0,0, prec)
+            sage: q01 = G.period(0,1, prec)
+            sage: q11 = G.period(1,1, prec)
+            sage: (q00g-q00).valuation() >= prec
+            True
+            sage: (q01g-q01).valuation() >= prec
+            True
+            sage: (q11g-q11).valuation() >= prec
+            True
         """
         if i in ZZ:
             assert j in ZZ
diff --git a/darmonpoints/util.py b/darmonpoints/util.py
index effc0ce..10fb26d 100644
--- a/darmonpoints/util.py
+++ b/darmonpoints/util.py
@@ -684,7 +684,7 @@ def polynomial_roots_old(f, K):
             rt = hensel_lift(f, r)
             if f(rt).valuation() < K.precision_cap() // 2:  # DEBUG!!
                 raise RuntimeError(
-                    "Couln't reach target valuation (%s), got only %s"
+                    "Could not reach target valuation (%s), got only %s"
                     % (K.precision_cap(), f(rt).valuation())
                 )
             ans.append(rt)
@@ -2166,7 +2166,7 @@ def simplification_isomorphism(G, return_inverse=False):
           To:   Finitely presented group <  |  >
           Defn: x |--> 1
 
-    ALGORITM:
+    ALGORITHM:
 
     Uses GAP.
     """