diff --git a/DoubleAndAdd.py b/DoubleAndAdd.py
index be5ade5..180bd89 100644
--- a/DoubleAndAdd.py
+++ b/DoubleAndAdd.py
@@ -9,6 +9,7 @@ def bits(n):
         yield n & 1
         n >>= 1
 
+
 def double_and_add(n, x):
     """
     Returns the result of n * x, computed using
diff --git a/FormalDerivative.py b/FormalDerivative.py
index 00f2288..e2440f2 100644
--- a/FormalDerivative.py
+++ b/FormalDerivative.py
@@ -1,12 +1,13 @@
 def FormalDerivative(Fx):
-  return [(x[0] * x[1], x[1] - 1) for x in Fx if x[1] > 0]
+    return [(x[0] * x[1], x[1] - 1) for x in Fx if x[1] > 0]
 
 
 def polyrep(Fx):
-  s = ["%dx ^ %d" % x for x in Fx]
-  return " + ".join(s)
+    s = ["%dx ^ %d" % x for x in Fx]
+    return " + ".join(s)
 
-fx = [(-1,6),(1,0)]
-print(polyrep(fx))
+
+fx = [(-1, 6), (1, 0)]
+print((polyrep(fx)))
 fx = FormalDerivative(fx)
-print(polyrep(fx))
+print((polyrep(fx)))
diff --git a/GaussGF2.py b/GaussGF2.py
index 79c7b02..0daef0f 100644
--- a/GaussGF2.py
+++ b/GaussGF2.py
@@ -2,24 +2,25 @@
 # Author Dario Clavijo 2021
 # based on https://www.cs.umd.edu/~gasarch/TOPICS/factoring/fastgauss.pdf, page2
 
+
 def Gauss_GF2(A):
-  h = len(A)
-  m = len(A[0])
-  marks = [False] * h
-  for j in range(0,m):
-    for i in range(0,h):
-      if A[i][j] == 1:
-          marks[i] = True
-          for k in range(j+1,m):
-            if A[i][k] == 1:
-              A[i][k] = (A[i][j] ^ A[i][k]) 
-          break
-  return marks, A
-              
-if __name__ == "__main__":
-  A  = [[1,1,0,0],[1,1,0,1],[0,1,1,1],[0,0,1,0],[0,0,0,1]]
-  marks,A = Gauss_GF2(A)
-  print(marks)
-  for row in A:
-    print(row)
+    h = len(A)
+    m = len(A[0])
+    marks = [False] * h
+    for j in range(0, m):
+        for i in range(0, h):
+            if A[i][j] == 1:
+                marks[i] = True
+                for k in range(j + 1, m):
+                    if A[i][k] == 1:
+                        A[i][k] = A[i][j] ^ A[i][k]
+                break
+    return marks, A
 
+
+if __name__ == "__main__":
+    A = [[1, 1, 0, 0], [1, 1, 0, 1], [0, 1, 1, 1], [0, 0, 1, 0], [0, 0, 0, 1]]
+    marks, A = Gauss_GF2(A)
+    print(marks)
+    for row in A:
+        print(row)
diff --git a/anomaly_detection.py b/anomaly_detection.py
index 7d08078..debc9ac 100644
--- a/anomaly_detection.py
+++ b/anomaly_detection.py
@@ -107,14 +107,14 @@ def test():
     MAD anomaly detection: [(5, 12, 3.8333333333333344)]
     """
     S = [2, 3, 5, 2, 3, 12, 5, 3, 4]
-    print("Series:", S)
-    print("Mean:", mean(S))
-    print("StdDev:", stddev(S))
-    print("Bounds:", bounds(S))
-    print("simple anomaly detection:", list(simple_anonaly_detection(S)))
-    print("zvalues:", list(zvalue(S)))
-    print("zvalue anomaly detection:", list(zvalue_anomaly_detection(S)))
-    print("MAD anomaly detection:", list(MAD_anomaly_detection(S)))
+    print(("Series:", S))
+    print(("Mean:", mean(S)))
+    print(("StdDev:", stddev(S)))
+    print(("Bounds:", bounds(S)))
+    print(("simple anomaly detection:", list(simple_anonaly_detection(S))))
+    print(("zvalues:", list(zvalue(S))))
+    print(("zvalue anomaly detection:", list(zvalue_anomaly_detection(S))))
+    print(("MAD anomaly detection:", list(MAD_anomaly_detection(S))))
 
 
 if __name__ == "__main__":
diff --git a/archimedes_pi.py b/archimedes_pi.py
index 2313e34..894814e 100644
--- a/archimedes_pi.py
+++ b/archimedes_pi.py
@@ -14,13 +14,13 @@ def approximation(sides):
     pi_guess = (inside + outside) / 2
     accuracy = (1 - (pi_guess - pi) / pi) * 100
 
-    print "sides:", sides
-    print "angle:", angle
-    print "inside:", inside
-    print "outside:", outside,
-    print "pi_guess:", pi_guess
-    print "accuracy:", accuracy
+    print("sides:", sides)
+    print("angle:", angle)
+    print("inside:", inside)
+    print("outside:", outside, end=" ")
+    print("pi_guess:", pi_guess)
+    print("accuracy:", accuracy)
 
 
-for i in xrange(10):
-    approximation(10 ** i)
+for i in range(10):
+    approximation(10**i)
diff --git a/barret.py b/barret.py
index dd093eb..c806126 100644
--- a/barret.py
+++ b/barret.py
@@ -1,26 +1,27 @@
-
-class BarretReduccer():
+class BarretReduccer:
     """
     https://en.wikipedia.org/wiki/Barrett_reduction
     """
+
     def __init__(self, m):
-        if m <=0: raise ValueError("Modulus must be positive")
-        if m & (m - 1) == 0: raise ValueError("Modulus must not be a power of 2")
+        if m <= 0:
+            raise ValueError("Modulus must be positive")
+        if m & (m - 1) == 0:
+            raise ValueError("Modulus must not be a power of 2")
         self.m = m
         self.shift = m.bit_length() << 1
         self.q = (1 << self.shift) // m
 
     def reduce(self, a):
-        #if 0 <= a <= self.m**2: raise ValueError("a must be >0 and < n^2")
+        # if 0 <= a <= self.m**2: raise ValueError("a must be >0 and < n^2")
         a -= ((a * self.q) >> self.shift) * self.m
-        if a >= self.m: a -= self.m
+        if a >= self.m:
+            a -= self.m
         return a
 
-def test():
-  br = BarretReduccer(12)
-  print(br.reduce(6))
-  print(br.reduce(13))
-  print(br.reduce(18))
 
-
-  
+def test():
+    br = BarretReduccer(12)
+    print((br.reduce(6)))
+    print((br.reduce(13)))
+    print((br.reduce(18)))
diff --git a/binGCD.py b/binGCD.py
index d6048c7..4d5856b 100644
--- a/binGCD.py
+++ b/binGCD.py
@@ -18,9 +18,9 @@ def gcd(a, b):
 
 
 def test():
-    print(gcd(115, 5))
-    print(gcd(5, 7))
-    print(gcd(10, 4))
+    print((gcd(115, 5)))
+    print((gcd(5, 7)))
+    print((gcd(10, 4)))
 
     a = (
         37975227936943673922808872755445627854565536638199
@@ -30,7 +30,7 @@ def test():
         40094690950920881030683735292761468389214899724061
         * 5846418214406154678836553182979162384198610505601062333
     )
-    print(gcd(a, b))
+    print((gcd(a, b)))
 
 
 if __name__ == "__main__":
diff --git a/binary_polinomial_factoring.sage b/binary_polinomial_factoring.sage
index 8bf21d5..b9be04f 100644
--- a/binary_polinomial_factoring.sage
+++ b/binary_polinomial_factoring.sage
@@ -64,7 +64,7 @@ def test():
                 ff += 1
             else:
                 nf += 1
-            print(n, i, ff, nf, f)
+            print((n, i, ff, nf, f))
         n += 1
 
 
diff --git a/birthdayparadox.py b/birthdayparadox.py
index ad380ff..3efd816 100644
--- a/birthdayparadox.py
+++ b/birthdayparadox.py
@@ -1,11 +1,11 @@
-def bpp(g,w):
-  """
-  Compute the probability of all elements in group g in the whole set w to have unique birthdays.
-  """
-  proba = 1
-  for i in range(w, w - g, -1):
-    proba *= (i/w)
-  return proba
+def bpp(g, w):
+    """
+    Compute the probability of all elements in group g in the whole set w to have unique birthdays.
+    """
+    proba = 1
+    for i in range(w, w - g, -1):
+        proba *= i / w
+    return proba
 
 
-print(bpp(23,365))
+print((bpp(23, 365)))
diff --git a/chirp_zeta_transform.py b/chirp_zeta_transform.py
index 05e95d0..3e9cf96 100644
--- a/chirp_zeta_transform.py
+++ b/chirp_zeta_transform.py
@@ -1,6 +1,6 @@
 # Given the following paper: https://www.nature.com/articles/s41598-019-50234-9
 
-#TODO: Implimentand test functions
+# TODO: Implimentand test functions
 # Found function implementation: ToeplitzMultiplyE
 # Suplemental paper (at bottom of previous link):
 # https://static-content.springer.com/esm/art%3A10.1038%2Fs41598-019-50234-9/MediaObjects/41598_2019_50234_MOESM1_ESM.pdf
@@ -13,36 +13,42 @@
 #
 # Depends on FFT and IFFT (import fftw)
 
-def CZT(x,M,W,A):
+
+def CZT(x, M, W, A):
     N = len(x)
-    X, r , c = [] * N, [] * N, [] * M
+    X, r, c = [] * N, [] * N, [] * M
     for k in range(0, N - 1):
         k2 = k * k
         X[k] = W[((k2) >> 1)] * A[-k] * x[k]
         r[k] = W[-((k2) >> 1)]
-    for k in range(0,M-1):
+    for k in range(0, M - 1):
         c[k] = W[-((k * k) >> 1)]
-    X = ToeplitzMultiplyE(r,c,X)
+    X = ToeplitzMultiplyE(r, c, X)
     for k in range(0, M - 1):
         X[k] = W[((k * k) >> 1)] * X[k]
     return X
 
-def ICZT(X,N,W,A):
+
+def ICZT(X, N, W, A):
     M = len(X)
     if M != N:
-        raise("M must == to N")
+        raise ("M must == to N")
     n = N
     x, p, u, z = [] * n, [] * n, [] * n, [] * n
     for k in range(0, n - 1):
         x[k] = W[((k * k) >> 1)] * X[k]
     p[0] = 1
-    for k in range(0,n-1):
+    for k in range(0, n - 1):
         p[k] = (p[k - 1]) * (W[k] - 1)
-    for k in range(0, n-1, 2):
-        u[k] = ((W[((k * k) << 1) - ((n << 1) - 1) + n * (n - 1)]) / (p[n - k - 1] * p[k]))
-    for k in range(1, n-1, 2):
-        u[k] = ((W[((k * k) << 1) - ((n << 1) - 1) + n * (n - 1)]) / (p[n - k - 1] * p[k])) * -1
-    for j = range(n-1,0,-1):
+    for k in range(0, n - 1, 2):
+        u[k] = (W[((k * k) << 1) - ((n << 1) - 1) + n * (n - 1)]) / (
+            p[n - k - 1] * p[k]
+        )
+    for k in range(1, n - 1, 2):
+        u[k] = (
+            (W[((k * k) << 1) - ((n << 1) - 1) + n * (n - 1)]) / (p[n - k - 1] * p[k])
+        ) * -1
+    for j in range(n - 1, 0, -1):
         u1[k] = u[k - j]
     u2 = u[0] + [0] * (n - 1)
     x1 = ToeplitzMultiplyE(u1, z, x)
@@ -54,5 +60,3 @@ def ICZT(X,N,W,A):
     for k in range(0, n - 1):
         x[k] = A[k] * W[-((k * k) >> 1)] * x[k]
     return x
-
-
diff --git a/contfrac.py b/contfrac.py
index fcd2b9c..f1c0c11 100644
--- a/contfrac.py
+++ b/contfrac.py
@@ -17,5 +17,5 @@ def cont_frac(a, b):
     return coef
 
 
-print(cont_frac(649, 200))
-print(cont_frac(17993, 90581))
+print((cont_frac(649, 200)))
+print((cont_frac(17993, 90581)))
diff --git a/dft.py b/dft.py
index 0ea7422..a245857 100644
--- a/dft.py
+++ b/dft.py
@@ -14,13 +14,13 @@ def dft_slow(signal):
 
     def Fk(signal, k):
         Freq = 0.0
-        for n in xrange(0, N):
+        for n in range(0, N):
             coef = math.exp((-2 * math.pi * k * n) / N)
             Freq += signal[n] * coef
         return Freq
 
     histogram = []
-    for k in xrange(0, N):
+    for k in range(0, N):
         histogram.append(Fk(signal, k))
     return histogram
 
@@ -30,4 +30,4 @@ def nyquist_norm(histogram):
     return [histogram[n] * 2 for n in range(0, N / 2)]
 
 
-print nyquist_norm(dft_slow(signal))
+print((nyquist_norm(dft_slow(signal))))
diff --git a/dx.py b/dx.py
index 6d8cd4d..89e344f 100644
--- a/dx.py
+++ b/dx.py
@@ -20,7 +20,7 @@ def dx(x):
 
 
 def g(x):
-    return x ** 2
+    return x**2
 
 
 dg = derivative(g)
@@ -42,10 +42,10 @@ def intf(b, a=0):
     return intf
 
 
-print "g(x)=", [g(x) for x in range(11)]
+print(("g(x)=", [g(x) for x in range(11)]))
 
-print "df g(x)=", [dg(x) for x in range(11)]
+print(("df g(x)=", [dg(x) for x in range(11)]))
 
 inv1 = integral(derivative(g))
 
-print "integral(df g(x))=", [inv1(x) for x in range(11)]
+print(("integral(df g(x))=", [inv1(x) for x in range(11)]))
diff --git a/e.py b/e.py
index 15c83a5..0a0abea 100644
--- a/e.py
+++ b/e.py
@@ -18,7 +18,7 @@ def direct_e():
 
 
 def taylor_e():
-    return sum(1.0 / (math.factorial(n)) for n in xrange(0, 15))
+    return sum(1.0 / (math.factorial(n)) for n in range(0, 15))
 
 
 # by limits definition:
@@ -30,29 +30,31 @@ def taylor_e():
 def lim_ddx_e(x):
     h = 1.0 / (534645555)
     return (
-        (math.e ** x * math.e ** h) - math.e ** x
+        (math.e**x * math.e**h) - math.e**x
     ) / h  # = math.e * ((math.e ** h - 1.0)/h)
 
 
 def test():
-  print(math.e)
-  print(direct_e())
-  print(taylor_e())
-  print(lim_ddx_e(1))
+    print((math.e))
+    print((direct_e()))
+    print((taylor_e()))
+    print((lim_ddx_e(1)))
 
 
-i = complex(0,1)
+i = complex(0, 1)
 pi = math.pi
 
+
 def exp(n, precision=100):
-    return sum((n ** x) / math.factorial(x) for x in range(0, precision))
+    return sum((n**x) / math.factorial(x) for x in range(0, precision))
+
 
 def test_exp():
-  print(exp(0))
-  print(exp(0.5))
-  print(exp(1))
-  print(exp(pi))
-  print(exp(i*pi))
+    print((exp(0)))
+    print((exp(0.5)))
+    print((exp(1)))
+    print((exp(pi)))
+    print((exp(i * pi)))
 
 
 test_exp()
diff --git a/euclid.py b/euclid.py
index 983def3..3c8138d 100644
--- a/euclid.py
+++ b/euclid.py
@@ -1,13 +1,12 @@
-
-
 def gcd1(a, b):
     while b != 0:
-        t = b 
-        b = a % b 
-        a = t 
+        t = b
+        b = a % b
+        a = t
     return a
 
-def gcd2(a,b):
+
+def gcd2(a, b):
     while a != b:
         if a > b:
             a = a - b
@@ -15,8 +14,9 @@ def gcd2(a,b):
             b = b - a
     return a
 
+
 def gcd3(a, b):
-    if b = 0:
-        return a; 
-    else
+    if b == 0:
+        return a
+    else:
         return gcd(b, a % b)
diff --git a/euler_integer_factorization.py b/euler_integer_factorization.py
index 4f2ac08..1b92e73 100644
--- a/euler_integer_factorization.py
+++ b/euler_integer_factorization.py
@@ -11,14 +11,14 @@ def euler(n):
     sol = []
 
     while a < end and len(sol) < 2:
-        b2 = n - a ** 2
+        b2 = n - a**2
         b = gmpy2.isqrt(b2)
-        b3 = b2 - b ** 2
+        b3 = b2 - b**2
         # print(b2,b,b3)
         if b3 == 0 and a != firstb and b != firstb:
             firstb = b
             sol.append([b, a])
-            print(b, a)
+            print((b, a))
         a += 1
 
     if len(sol) < 2:
@@ -38,11 +38,11 @@ def euler(n):
 
     # n = (k**2 + h**2)*(m**2 + l**2)
 
-    return (k ** 2 + h ** 2) // 4, (l ** 2 + m ** 2) // 4
+    return (k**2 + h**2) // 4, (l**2 + m**2) // 4
 
 
 import sys
 
 r = euler(int(sys.argv[1]))
 print(r)
-print(r[0] * r[1])
+print((r[0] * r[1]))
diff --git a/fermat_factor.py b/fermat_factor.py
index cd89709..cd8334e 100644
--- a/fermat_factor.py
+++ b/fermat_factor.py
@@ -8,13 +8,13 @@
 
 def fermat(n):
     a = gmpy2.isqrt(n)
-    if a ** 2 == n:
+    if a**2 == n:
         return a, a
-    b2 = (a ** 2) - n
+    b2 = (a**2) - n
     step = 0
     while not gmpy2.is_square(b2):
         a += 1
-        b2 = (a ** 2) - n
+        b2 = (a**2) - n
         # step += 1
         # print(n,step,a,b2)
         ib2 = gmpy2.isqrt(b2)
@@ -23,10 +23,10 @@ def fermat(n):
 
 def test():
     n = 5959
-    print("-" * 40)
+    print(("-" * 40))
     a, b = fermat(n)
-    print("-" * 40)
-    print(n, a, b, n == a * b)
+    print(("-" * 40))
+    print((n, a, b, n == a * b))
 
 
 test()
diff --git a/fft.py b/fft.py
index ea51c0c..0f59917 100644
--- a/fft.py
+++ b/fft.py
@@ -1,16 +1,18 @@
 # Author Dario Clavijo 2023
 # based on https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm
 from gmpy2 import *
+
 pi_i = (-2j) * complex(gmpy2.const_pi())
 
+
 def fft(X):
-  if (N:=len(X)) == 1:
+    if (N := len(X)) == 1:
+        return X
+    N2 = N >> 1
+    E = fft(X[::2])
+    O = fft(X[1::2])
+    for k in range(0, N2):
+        q = exp((pi_i * k) / N) * O[k]
+        X[k] = E[k] + q
+        X[k + N2] = E[k] - q
     return X
-  N2 = N >> 1
-  E = fft(X[::2])
-  O = fft(X[1::2])
-  for k in range(0, N2):
-    q = exp((pi_i * k) / N ) * O[k]
-    X[k] = E[k] + q
-    X[k + N2] = E[k] - q
-  return X
diff --git a/fibonacci_eigen.py b/fibonacci_eigen.py
index 3b7e1d5..2cd4ef9 100644
--- a/fibonacci_eigen.py
+++ b/fibonacci_eigen.py
@@ -34,36 +34,40 @@
 λ2 = (1 - sqrt(5)) / 2
 """
 
+
 def find_fibonacci_eigenvalues():
-  "It computes lambda(λ) from: (A-λI)v = 0 " 
-  _lambda = symbols('λ') # Lambda escalar placeholder
-  A = Matrix(([1, 1], [1, 0])) # Fibonacci general equation matrix
-  I = Matrix(([1, 0], [0, 1])) # Identity matrix
-  L = _lambda * I # Lambda diagonal 
-  E = A-L # Resulting matrix = A-λI
-  print("A:",A,"\nI:",I,"\nL:",L,"\nE:",E)
-  D = E.det() # compute the matrix determinant
-  print("D:",D)
-  lambdas = solve(D) # Get roots, solve equation by bhaskara
-  la0 = lambdas[0] # root 0
-  la1 = lambdas[1] # root 1
-  print("λ0:",la0,"\nλ1:",la1)
-  return (la0,la1)
-
-def binet_fib_eig(la0,la1,k):
-  ik = pow(la0,k) - pow(la1,k) / (la0 -la1) # Binet's formula: (a^n-b^n)/(a-b)
-  return round(ik)
+    "It computes lambda(λ) from: (A-λI)v = 0"
+    _lambda = symbols("λ")  # Lambda escalar placeholder
+    A = Matrix(([1, 1], [1, 0]))  # Fibonacci general equation matrix
+    I = Matrix(([1, 0], [0, 1]))  # Identity matrix
+    L = _lambda * I  # Lambda diagonal
+    E = A - L  # Resulting matrix = A-λI
+    print(("A:", A, "\nI:", I, "\nL:", L, "\nE:", E))
+    D = E.det()  # compute the matrix determinant
+    print(("D:", D))
+    lambdas = solve(D)  # Get roots, solve equation by bhaskara
+    la0 = lambdas[0]  # root 0
+    la1 = lambdas[1]  # root 1
+    print(("λ0:", la0, "\nλ1:", la1))
+    return (la0, la1)
+
+
+def binet_fib_eig(la0, la1, k):
+    ik = pow(la0, k) - pow(la1, k) / (la0 - la1)  # Binet's formula: (a^n-b^n)/(a-b)
+    return round(ik)
+
 
 def test():
-  la0,la1 = find_fibonacci_eigenvalues()
-  for i in range(3,30):
-    a = fib(i)
-    b = binet_fib_eig(la0,la1,i)
-    print(i,a,b,a==b)
+    la0, la1 = find_fibonacci_eigenvalues()
+    for i in range(3, 30):
+        a = fib(i)
+        b = binet_fib_eig(la0, la1, i)
+        print((i, a, b, a == b))
 
 
 def test2():
-  return
+    return
+
 
-#if __name__ == "__main__":
+# if __name__ == "__main__":
 #  test()
diff --git a/int2limbs.py b/int2limbs.py
index 2a09d99..6e39fec 100644
--- a/int2limbs.py
+++ b/int2limbs.py
@@ -1,11 +1,11 @@
 def int_to_limbs(n, base):
     """
     Converts an integer to its limbs (digits in a given base).
-    
+
     Args:
     - n: the integer to convert
     - base: the base in which to represent the integer
-    
+
     Returns:
     A list of limbs (digits) representing the integer in the given base.
     """
@@ -14,16 +14,16 @@ def int_to_limbs(n, base):
         n, remainder = divmod(n, base)
         limbs.append(remainder)
     return limbs[::-1]
-    
+
 
 def limbs_to_int(limbs, base):
     """
     Converts a list of limbs (digits in a given base) to an integer.
-    
+
     Args:
     - limbs: a list of digits representing the integer in the given base
     - base: the base in which the integer is represented
-    
+
     Returns:
     An integer representing the value of the given limbs in the given base.
     """
diff --git a/iscongruent.py b/iscongruent.py
index 9a289b3..daae28c 100644
--- a/iscongruent.py
+++ b/iscongruent.py
@@ -1,9 +1,9 @@
 def iscongruent(a, b, n):
     """a and b are congruent modulo n only if the absolute difference of both divides n"""
-    #return (abs(a - b) // n) == 1
+    # return (abs(a - b) // n) == 1
     return (abs(a - b) % n) == 0
 
 
 if __name__ == "__main__":
-    print(iscongruent(4, 6, 12))
-    print(iscongruent(5, 12, 60))
+    print((iscongruent(4, 6, 12)))
+    print((iscongruent(5, 12, 60)))
diff --git a/karatsuba.py b/karatsuba.py
index fd22653..f6a91da 100644
--- a/karatsuba.py
+++ b/karatsuba.py
@@ -1,17 +1,21 @@
 import math
 from functools import cache
 
-digits = lambda n: int(math.log10(n))+1
+digits = lambda n: int(math.log10(n)) + 1
+
 
 def pow10(x, p):
-  # Computes x*(10^p) efficiently
-  for _ in range(1, p): x = (x << 3) + (x << 1)
-  return x
+    # Computes x*(10^p) efficiently
+    for _ in range(1, p):
+        x = (x << 3) + (x << 1)
+    return x
+
 
 @cache
 def karatsuba(x, y):
-    if (x < 10 or y < 10): return x * y 
-    m = max(digits(x), digits(y)) 
+    if x < 10 or y < 10:
+        return x * y
+    m = max(digits(x), digits(y))
     m2 = m // 2
     m3 = pow10(10, m2)
     a, b = divmod(x, m3)
@@ -21,6 +25,10 @@ def karatsuba(x, y):
     z2 = karatsuba(a, c)
     return (z2 * (pow10(10, (m2 << 1)))) + ((z1 - z2 - z0) * m3) + z0
 
-a,b = 132621000145968689486197181919225338865594781061448283090178076569, 6476706838600125160840627922710127689
-print(karatsuba(a,b))
-print(a*b)
+
+a, b = (
+    132621000145968689486197181919225338865594781061448283090178076569,
+    6476706838600125160840627922710127689,
+)
+print((karatsuba(a, b)))
+print((a * b))
diff --git a/lazy_erathostenes.py b/lazy_erathostenes.py
index a65d4b7..222f9a6 100644
--- a/lazy_erathostenes.py
+++ b/lazy_erathostenes.py
@@ -1,13 +1,16 @@
 from itertools import islice
 
+
 def nats(n):
-  yield n
-  yield from nats(n+1)
+    yield n
+    yield from nats(n + 1)
+
 
 def sieve(s):
-  #print("sieve",s)
-  n = next(s)
-  yield n
-  yield from sieve(i for i in s if i % n != 0)
+    # print("sieve",s)
+    n = next(s)
+    yield n
+    yield from sieve(i for i in s if i % n != 0)
+
 
-print(list(islice(sieve(nats(2)), 30)))
+print((list(islice(sieve(nats(2)), 30))))
diff --git a/legendre_symbol.py b/legendre_symbol.py
index a816bd3..e672f07 100644
--- a/legendre_symbol.py
+++ b/legendre_symbol.py
@@ -19,7 +19,7 @@ def factor(n0):
                 factors.append(i)
     if not factors:
         factors = [n]
-    print("factor(%d)=%s" % (n0, str(factors)))
+    print(("factor(%d)=%s" % (n0, str(factors))))
     return factors
 
 
@@ -29,33 +29,33 @@ def legendre_naive(p, q):
     else:
         r = -1
         for x in range(0, q):
-            y = (x ** 2) % q
+            y = (x**2) % q
             z = p % q
             # print x,y,z
             if y == z:
                 r = 1
                 break
-    print("legendre_naive(%d|%d)=%d" % (p, q, r))
+    print(("legendre_naive(%d|%d)=%d" % (p, q, r)))
     return r
 
 
 def legendre_prop0(p, q):
-    print("legendre_prop0(%d|%d)" % (p, q))
+    print(("legendre_prop0(%d|%d)" % (p, q)))
     phi = (p - 1) * (q - 1)
     _pow = pow(-1, int(phi // 4))
     return _pow * legendre_naive(q, p)
 
 
 def legendre_prop1(p, q):
-    print("legendre_prop1 (2| %d)" % q)
+    print(("legendre_prop1 (2| %d)" % q))
     if p == 2:
-        return pow((-1), (((q ** 2) - 1) // 8))
+        return pow((-1), (((q**2) - 1) // 8))
     else:
         raise Exception("p != 2")
 
 
 def legendre_prop2(ab, q):
-    print("legendre_prop2(%d,%d)" % (ab, q))
+    print(("legendre_prop2(%d,%d)" % (ab, q)))
     tmp = 1
     for f in factor(ab):
         if f == 2:
@@ -67,7 +67,7 @@ def legendre_prop2(ab, q):
 
 def _legendre(p, q):
     r = p % q
-    print(">legendre(%d|%d) r=%d" % (p, q, r))
+    print((">legendre(%d|%d) r=%d" % (p, q, r)))
     if r == 0:
         ret = 0
     else:
@@ -90,7 +90,7 @@ def _legendre(p, q):
                     tmp *= _legendre(q, f)
             # tmp = legendre_prop2(p1,q)
             ret = tmp
-    print("<legendre(%d|%d) r=%d -> %d" % (p, q, r, ret))
+    print(("<legendre(%d|%d) r=%d -> %d" % (p, q, r, ret)))
     return ret
 
 
diff --git a/legendre_symbol_small.py b/legendre_symbol_small.py
index 0e19fd3..8b45132 100644
--- a/legendre_symbol_small.py
+++ b/legendre_symbol_small.py
@@ -1,21 +1,23 @@
 #!/usr/bin/env python3
 # Author Dario Clavijo 2021
 
-def legendre_symbol(a,p)
-  """ legendre symbol computation"""
-  
-  ls = pow(a,(p-1)//2,p)
-  
-  if ls == p - 1:
-     return -1
-  else:
-     return ls
+
+def legendre_symbol(a, p):
+    """legendre symbol computation"""
+
+    ls = pow(a, (p - 1) // 2, p)
+
+    if ls == p - 1:
+        return -1
+    else:
+        return ls
+
 
 def test():
-  for a in range(1,10):
-    for p in range(1,10):
-      print(a,p,legendre_symbol(a,p))
+    for a in range(1, 10):
+        for p in range(1, 10):
+            print((a, p, legendre_symbol(a, p)))
+
 
 if __name__ == "__main__":
-  test()
-  
+    test()
diff --git a/leyland_primes.py b/leyland_primes.py
index 79a0cc3..db04709 100644
--- a/leyland_primes.py
+++ b/leyland_primes.py
@@ -11,7 +11,7 @@ def leyland_primes_naive(m):
     tmp = []
     for i in range(2, m):
         for j in range(2, m):
-            x = i ** j + j ** i
+            x = i**j + j**i
             if gmpy2.is_prime(x):
                 if x not in tmp:
                     tmp.append(x)
@@ -24,11 +24,11 @@ def leyland_primes_oeis(N):
     n = 1
     n1 = 1
     while n < N:
-        n = 2 * n1 ** n1
+        n = 2 * n1**n1
         k = n1 + 1
         while k < N:
             if gmpy2.gcd(n, k) == 1:
-                a = n ** k + k ** n
+                a = n**k + k**n
                 if a > N:
                     break
                 if gmpy2.is_prime(a):
@@ -41,4 +41,4 @@ def leyland_primes_oeis(N):
 
 
 if __name__ == "__main__":
-    print(leyland_primes_oeis(10 ** 10))
+    print((leyland_primes_oeis(10**10)))
diff --git a/ln2.py b/ln2.py
index 08bd60a..03c9b15 100644
--- a/ln2.py
+++ b/ln2.py
@@ -9,4 +9,4 @@ def taylor_ln2():
     return sum
 
 
-print taylor_ln2()
+print((taylor_ln2()))
diff --git a/math_rels.py b/math_rels.py
index b91ec1a..19dc618 100644
--- a/math_rels.py
+++ b/math_rels.py
@@ -27,7 +27,7 @@ def close(a, b, e):
 t = -min(-a, a)
 x = max(-a, a)
 y = abs(a)
-z = math.sqrt(a ** 2)
+z = math.sqrt(a**2)
 w = a * sgn(a)
 
-print t == x == y == z == w
+print((t == x == y == z == w))
diff --git a/mean.py b/mean.py
index a7405d3..f688d87 100644
--- a/mean.py
+++ b/mean.py
@@ -20,4 +20,4 @@ def GM(data):
 # harmonic mean
 def HM(data):
     accum = sum(1 / data[i] for i in range(0, len(data) - 1))
-    return len(data) * (accum ** -1)
+    return len(data) * (accum**-1)
diff --git a/minmax_math.py b/minmax_math.py
index 616df3f..f780bb8 100644
--- a/minmax_math.py
+++ b/minmax_math.py
@@ -11,5 +11,5 @@ def _max(a, b):
     return 0.5 * (a + b + math.sqrt((a - b) ** 2))
 
 
-print _min(5, 7)
-print _max(5, 7)
+print((_min(5, 7)))
+print((_max(5, 7)))
diff --git a/modulo.py b/modulo.py
index 47dc2ba..2f98240 100644
--- a/modulo.py
+++ b/modulo.py
@@ -23,7 +23,7 @@ def test():
         a = subMod(i, x)
         r, m = divMod(i, x)
         b = x % i
-        print(x, i, a, b, a == b, r, m, x == (r * i) + m)
+        print((x, i, a, b, a == b, r, m, x == (r * i) + m))
 
 
 if __name__ == "__main__":
diff --git a/multiplicativeOrder.py b/multiplicativeOrder.py
index b9695d7..2d7c0c1 100644
--- a/multiplicativeOrder.py
+++ b/multiplicativeOrder.py
@@ -43,11 +43,11 @@ def CyclicGroup(g, n):
 
 
 if __name__ == "__main__":
-    print("ZnMultGroup(10)", ZnMultGroup(10))
-    print("ZnMultGroup(11)", ZnMultGroup(11))
-    print("CyclicGroup(2,11)", CyclicGroup(2, 11))
+    print(("ZnMultGroup(10)", ZnMultGroup(10)))
+    print(("ZnMultGroup(11)", ZnMultGroup(11)))
+    print(("CyclicGroup(2,11)", CyclicGroup(2, 11)))
     for i in range(2, 10):
         for j in range(i, 10):
             if gcd(i, j) == 1:
                 # print("Zn:",j,ZnMultGroup(j))
-                print("multiplicativeOrder", i, j, multiplicativeOrder(i, j))
+                print(("multiplicativeOrder", i, j, multiplicativeOrder(i, j)))
diff --git a/newton_method.py b/newton_method.py
index edee424..88abcf0 100644
--- a/newton_method.py
+++ b/newton_method.py
@@ -1,27 +1,35 @@
 """
 https://en.wikipedia.org/wiki/Newton%27s_method
 """
-def f(x):             
-	return x**2 - 2   # f(x) = x^2 - 2
+
+
+def f(x):
+    return x**2 - 2  # f(x) = x^2 - 2
+
 
 def f_prime(x):
-	return 2*x        # f'(x) = 2x
+    return 2 * x  # f'(x) = 2x
+
 
 def newtons_method(
-    x0,               # The initial guess
-    f,                # The function whose root we are trying to find
-    f_prime,          # The derivative of the function
-    tolerance,        # 7-digit accuracy is desired
-    epsilon,          # Do not divide by a number smaller than this
-    max_iterations,   # The maximum number of iterations to execute
-    ):
-	for _ in range(max_iterations):
-		y = f(x0)
-		yprime = f_prime(x0)
-		if abs(yprime) < epsilon:       # Stop if the denominator is too small
-		    break
-		x1 = x0 - y / yprime            # Do Newton's computation
-		if abs(x1 - x0) <= tolerance:   # Stop when the result is within the desired tolerance
-		    return x1                   # x1 is a solution within tolerance and maximum number of iterations
-		x0 = x1                         # Update x0 to start the process again
-	return None
+    x0,  # The initial guess
+    f,  # The function whose root we are trying to find
+    f_prime,  # The derivative of the function
+    tolerance,  # 7-digit accuracy is desired
+    epsilon,  # Do not divide by a number smaller than this
+    max_iterations,  # The maximum number of iterations to execute
+):
+    for _ in range(max_iterations):
+        y = f(x0)
+        yprime = f_prime(x0)
+        if abs(yprime) < epsilon:  # Stop if the denominator is too small
+            break
+        x1 = x0 - y / yprime  # Do Newton's computation
+        if (
+            abs(x1 - x0) <= tolerance
+        ):  # Stop when the result is within the desired tolerance
+            return (
+                x1  # x1 is a solution within tolerance and maximum number of iterations
+            )
+        x0 = x1  # Update x0 to start the process again
+    return None
diff --git a/oeis_closed_form.py b/oeis_closed_form.py
index f2507cc..a4337c6 100644
--- a/oeis_closed_form.py
+++ b/oeis_closed_form.py
@@ -4,44 +4,47 @@
 import requests
 import sqlite3
 
+
 def getSequence(id):
-  f = requests.get(f"https://oeis.org/search?fmt=json&q=id:{id}")
-  doc = json.loads(f.content)
-  return [int(x) for x in doc['results'][0]['data'].split(",")]
+    f = requests.get(f"https://oeis.org/search?fmt=json&q=id:{id}")
+    doc = json.loads(f.content)
+    return [int(x) for x in doc["results"][0]["data"].split(",")]
+
 
 def guessSequence(lst):
-  C = CFiniteSequences(QQ)
-  if (s:= C.guess(lst)) == 0:
-    return
-  else:
-    return s.closed_form()
+    C = CFiniteSequences(QQ)
+    if (s := C.guess(lst)) == 0:
+        return
+    else:
+        return s.closed_form()
 
 
-def checkSequence(id,items=10):
-  data = getSequence(id)
-  seq = data[:items]
-  print(id,seq)
-  if len(seq) > 7 and (r:=guessSequence(seq)) is not None:
-    seq = data[:items * 5]
-    return guessSequence(seq)
+def checkSequence(id, items=10):
+    data = getSequence(id)
+    seq = data[:items]
+    print((id, seq))
+    if len(seq) > 7 and (r := guessSequence(seq)) is not None:
+        seq = data[: items * 5]
+        return guessSequence(seq)
 
 
 def procfile():
-  IDS=[line.rstrip() for line in open(sys.argv[1],'r').readlines()]
-  for id in IDS:
-    print(id, getSequence(id))
+    IDS = [line.rstrip() for line in open(sys.argv[1], "r").readlines()]
+    for id in IDS:
+        print((id, getSequence(id)))
 
 
 def proc2():
-  #conn = sqlite3.connect('oeis.db')
-  #cur = conn.cursor()
-  #cur.execute("CREATE TABLE sequence(id, name,data,formula)")
-  #for row in cur.execute("SELECT id,name,data,formula FROM sequence ORDER BY id"):
-  n = 1
-  while True:
-    id = "A%06d" % n
-    print(id, checkSequence(id))
-    n+=1
+    # conn = sqlite3.connect('oeis.db')
+    # cur = conn.cursor()
+    # cur.execute("CREATE TABLE sequence(id, name,data,formula)")
+    # for row in cur.execute("SELECT id,name,data,formula FROM sequence ORDER BY id"):
+    n = 1
+    while True:
+        id = "A%06d" % n
+        print((id, checkSequence(id)))
+        n += 1
+
 
 proc2()
-#procfile()
+# procfile()
diff --git a/pascal_triangle.py b/pascal_triangle.py
index d83f536..442cc8b 100644
--- a/pascal_triangle.py
+++ b/pascal_triangle.py
@@ -10,4 +10,4 @@ def hexify(i):
 
 
 for i in range(0, 257):
-    print bin(nk(256, i))
+    print((bin(nk(256, i))))
diff --git a/pisano_period.py b/pisano_period.py
index 24056ab..5c1c8e2 100644
--- a/pisano_period.py
+++ b/pisano_period.py
@@ -7,10 +7,11 @@ def pisano_period(m):
     """
     A001175
     """
-    if m == 1: return 1
+    if m == 1:
+        return 1
     a = b = 1
     n = i = 0
-    while i <= (m ** 2):
+    while i <= (m**2):
         b, a = a, (a + b) % m
         if a == 1 and b == 0:
             return n + 2
@@ -19,4 +20,4 @@ def pisano_period(m):
 
 if __name__ == "__main__":
     for i in range(1, 100):
-        print(pisano_period(i))
+        print((pisano_period(i)))
diff --git a/polysolv.py b/polysolv.py
index e9b8f1e..53750d7 100644
--- a/polysolv.py
+++ b/polysolv.py
@@ -72,29 +72,29 @@ def poly_to_text(poly, grade=0):
 
 
 def poly_synthetic_div_complete_step(poly, grade):
-    print "-" * 60
-    print "Grade:", grade
-    print "P = ", poly_to_text(poly, grade)
-    print "Coefs: P =", poly
+    print(("-" * 60))
+    print(("Grade:", grade))
+    print(("P = ", poly_to_text(poly, grade)))
+    print(("Coefs: P =", poly))
     rationals = get_rationals(poly, grade)
-    print "rationals:", rationals
+    print(("rationals:", rationals))
     term = []
     if len(rationals) > 0:
         for Q in rationals:
             R = poly_synthetic_div(poly, Q)
             if R[-1] == 0:
-                print "Q =", Q
-                print ("%d is divisor" % Q)
-                print "Coefs: R =", R
-                print "R = ", poly_to_text(R, grade - 1)
+                print(("Q =", Q))
+                print(("%d is divisor" % Q))
+                print(("Coefs: R =", R))
+                print(("R = ", poly_to_text(R, grade - 1)))
             if R[-1] == 0:
                 term = ["(x + %d)" % -Q]
                 # break
                 return R[:-1], term
-        print "No divisor found"
+        print("No divisor found")
         return None
     else:
-        print "No rationals"
+        print("No rationals")
     return None
 
 
@@ -109,9 +109,9 @@ def factor_poly(Poly, Grade):
         else:
             break
     terms += ["(%s)" % poly_to_text(P, grade)]
-    print "=" * 60
+    print(("=" * 60))
     s = sanitize("".join(terms))
-    print "Result:", s
+    print(("Result:", s))
 
 
 P = [1, 1, -11, -5, 30]
diff --git a/powmod.py b/powmod.py
index 13e2a35..daeefa0 100644
--- a/powmod.py
+++ b/powmod.py
@@ -16,8 +16,8 @@ def test(A):
     print(A)
     a = powmod(A[0], A[1], A[2])
     b = pow(A[0], A[1], A[2])
-    print(a, b)
+    print((a, b))
     return a == b
 
 
-print(test((13, 484, 497)))
+print((test((13, 484, 497))))
diff --git a/prosthaphaeresis.py b/prosthaphaeresis.py
index d936a56..b5cc2a9 100644
--- a/prosthaphaeresis.py
+++ b/prosthaphaeresis.py
@@ -41,8 +41,8 @@ def pmultcossin(a, b):
 a = 0.17
 b = 0.37
 
-print a * b
-print pmultcos(a, b)
-print pmultsin(a, b)
-print pmultsincos(a, b)
-print pmultcossin(a, b)
+print((a * b))
+print((pmultcos(a, b)))
+print((pmultsin(a, b)))
+print((pmultsincos(a, b)))
+print((pmultcossin(a, b)))
diff --git a/quartersquare.py b/quartersquare.py
index cffe48e..44205e6 100644
--- a/quartersquare.py
+++ b/quartersquare.py
@@ -16,6 +16,6 @@ def logmult(a, b):
 
 a = 255
 b = 255
-print logmult(a, b)
-print qsqr(a, b)
-print a * b
+print((logmult(a, b)))
+print((qsqr(a, b)))
+print((a * b))
diff --git a/ramanujan_sum.py b/ramanujan_sum.py
index 4cdf0a0..4d02ec6 100644
--- a/ramanujan_sum.py
+++ b/ramanujan_sum.py
@@ -5,8 +5,10 @@
 
 ipi2 = (2j) * complex(gmpy2.const_pi())
 
-def c(q,n):
-  nipi2 = ipi2 * n
-  return sum(exp(nipi2 * (a/q)) for a in range(1, q+1) if gcd(a,q) == 1)
 
-print([c(1,n).real for n in range(1, 31)])
+def c(q, n):
+    nipi2 = ipi2 * n
+    return sum(exp(nipi2 * (a / q)) for a in range(1, q + 1) if gcd(a, q) == 1)
+
+
+print([c(1, n).real for n in range(1, 31)])
diff --git a/redc.py b/redc.py
index 8377bd6..0e9a6a0 100644
--- a/redc.py
+++ b/redc.py
@@ -1,13 +1,15 @@
 import gmpy2
+from functools import reduce
+
 
 class MontgomeryReducer:
-    
-    def __init__(self, m, gmpy = None):
-        if m is None or m < 3 or m & 1 == 0: raise ValueError("Modulus must be >= 3 and odd.")
+    def __init__(self, m, gmpy=None):
+        if m is None or m < 3 or m & 1 == 0:
+            raise ValueError("Modulus must be >= 3 and odd.")
         self.mod = m
         self.gmpy = gmpy
         self.bits = ((m.bit_length() >> 3) + 1) << 3
-        self.reducer = 1 << self.bits 
+        self.reducer = 1 << self.bits
         self.mask = self.reducer - 1
         if self.gmpy is None:
             self.one = self.reducer % m
@@ -16,13 +18,14 @@ def __init__(self, m, gmpy = None):
         else:
             self.one = self.gmpy.f_mod(self.reducer, m)
             self.reciprocal = self.gmpy.powmod(self.reducer, -1, m)
-            self.factor = self.gmpy.f_div(gmpy2.mul(self.reducer, self.reciprocal) - 1, m)
+            self.factor = self.gmpy.f_div(
+                gmpy2.mul(self.reducer, self.reciprocal) - 1, m
+            )
 
-    
     def encode(self, x):
         if self.gmpy is None:
             return (x << self.bits) % self.mod
-        else:  
+        else:
             return self.gmpy.f_mod(x << self.bits, self.mod)
 
     def decode(self, x):
@@ -38,14 +41,15 @@ def reduce(self, op):
         else:
             tmp = self.gmpy.mul((op & self.mask), self.factor) & self.mask
             red = op + self.gmpy.mul(tmp, self.mod) >> self.bits
-        if (red > self.mod): red -= self.mod
+        if red > self.mod:
+            red -= self.mod
         return red
 
     def mul(self, x, y):
         if self.gmpy is None:
-          return self.reduce(x * y)
+            return self.reduce(x * y)
         else:
-          return self.reduce(self.gmpy.mul(x,y))
+            return self.reduce(self.gmpy.mul(x, y))
 
     def sum(self, x, y):
         return self.reduce(x + y)
@@ -58,50 +62,55 @@ def __mul__(self, other):
             return MontgomeryReducer(self.mod * other.mod)
         else:
             return MontgomeryReducer(self.gmpy.mul(self.mod, other.mod))
-    
+
     def pow(self, x, y):
-        if y < 0: return gmpy2.powmod(x, y, self.mod)
+        if y < 0:
+            return gmpy2.powmod(x, y, self.mod)
         z = self.one
         while y == 1:
-            if y & 1 == 1: z = self.mul(z, x)
+            if y & 1 == 1:
+                z = self.mul(z, x)
             x = self.mul(x, x)
             y >>= 1
         return z
 
+
 def test():
     import gmpy2
     from operator import mul
     from functools import reduce
+
     p = 2
-    TMP=[]
-    for _ in range(1,210):
+    TMP = []
+    for _ in range(1, 210):
         p = int(gmpy2.next_prime(p))
-        TMP.append(MontgomeryReducer(p,gmpy=gmpy2))
-    C = reduce(mul,TMP)
-    print(C.mod)
-    print(C.decode(C.encode(1)))
+        TMP.append(MontgomeryReducer(p, gmpy=gmpy2))
+    C = reduce(mul, TMP)
+    print((C.mod))
+    print((C.decode(C.encode(1))))
     c2 = C.encode(2)
     c3 = C.encode(13)
-    print(C.decode(C.mul(c2, c3)))
-    print(C.decode(C.decode(c2 * c3)))
+    print((C.decode(C.mul(c2, c3))))
+    print((C.decode(C.decode(c2 * c3))))
 
     L = 4
     x = 65539
     import time
+
     t0 = time.time()
-    for i in range(1, 10 ** L):
+    for i in range(1, 10**L):
         C.decode(C.pow(C.encode(x), i))
     t1 = time.time()
-    print("Montgomery powmod",t1 - t0)
+    print(("Montgomery powmod", t1 - t0))
     m = C.mod
-    for i in range(1, 10 ** L):
+    for i in range(1, 10**L):
         pow(x, i, m)
     t2 = time.time()
-    print("python powmod", t2 - t1)
-    for i in range(1, 10 ** L):
+    print(("python powmod", t2 - t1))
+    for i in range(1, 10**L):
         gmpy2.powmod(x, i, m)
     t3 = time.time()
-    print("gmpy2 powmod",t3 - t2)
+    print(("gmpy2 powmod", t3 - t2))
 
 
 test()
diff --git a/roots.py b/roots.py
index 4b425ce..ade6de0 100644
--- a/roots.py
+++ b/roots.py
@@ -9,5 +9,5 @@ def nth_root(n, x):
     return math.e ** ((1.0 / n) * math.log(x))
 
 
-print nth_root(2, 2)
-print nth_root(3, 8)
+print((nth_root(2, 2)))
+print((nth_root(3, 8)))
diff --git a/ruffini.py b/ruffini.py
index 553ef22..66a91e8 100644
--- a/ruffini.py
+++ b/ruffini.py
@@ -7,7 +7,7 @@
 
 # tells if a n is prime or not
 def isPrime(n):
-    return all(n % i != 0 for i in range(2, int(n ** 0.5) + 1))
+    return all(n % i != 0 for i in range(2, int(n**0.5) + 1))
 
 
 # finds the prime factors of n
@@ -51,18 +51,18 @@ def ruffini_step(P, D):
 
 
 def ruffini_test(P):
-    print "P=", P
+    print(("P=", P))
     TI = P[len(P) - 1]
     D = [1] + factors(abs(TI)) + [abs(TI)]
     D = addnegatives(D)
     D.sort()
-    print "D=", D
+    print(("D=", D))
 
     tmp = P
     for k in range(len(P) - 3):
         ret = ruffini_step(tmp, D)
         tmp = ret[1]
-        print "Grade=", len(P) - k, "root=", ret[0], "coefs=", ret[1]
+        print(("Grade=", len(P) - k, "root=", ret[0], "coefs=", ret[1]))
 
 
 # P = [1,-7,13,23,-78]
diff --git a/seqmult.py b/seqmult.py
index 94496f3..17098b3 100644
--- a/seqmult.py
+++ b/seqmult.py
@@ -9,9 +9,9 @@ def SeqMult(s):
         if l % 2 != 0:
             s += [1]
             l += 1
-        s = list(map(mul, s[:l // 2], s[l // 2 :]))
+        s = list(map(mul, s[: l // 2], s[l // 2 :]))
         l = len(s)
     return s[0]
 
 
-print(SeqMult([1, 2, 3, 4, 5]) == (1 * 2 * 3 * 4 * 5))
+print((SeqMult([1, 2, 3, 4, 5]) == (1 * 2 * 3 * 4 * 5)))
diff --git a/sieve.py b/sieve.py
index db31a26..9eda6b8 100644
--- a/sieve.py
+++ b/sieve.py
@@ -7,10 +7,10 @@ def sieve_Erathostenes(n):
     for i in range(2, int(math.sqrt(n))):
         if sieves[i] == True:
             for j in range(0, n):
-                h = (i ** 2) + (j * i)
+                h = (i**2) + (j * i)
                 if h < n:
                     sieves[h] = False
     return [k for k in range(2, n) if sieves[k] == True]
 
 
-print sieve_Erathostenes(121)
+print((sieve_Erathostenes(121)))
diff --git a/sign.py b/sign.py
index 1527146..67580b7 100644
--- a/sign.py
+++ b/sign.py
@@ -17,6 +17,6 @@ def almostEqual(a, b, epsilon):
     return _abs(a - b) <= epsilon
 
 
-print sign(1)
-print sign(0)
-print sign(-1)
+print((sign(1)))
+print((sign(0)))
+print((sign(-1)))
diff --git a/sqrttest.py b/sqrttest.py
index fed36e7..0d8bf00 100644
--- a/sqrttest.py
+++ b/sqrttest.py
@@ -16,14 +16,15 @@ def sqrt(x):
 
 
 def sqrt(n):
-  x = n
-  y = 1
-  e = 0.00000001
-  while (x - y) > e:
-    x = (x + y) / 2
-    y = n / x
-  return x
-  
+    x = n
+    y = 1
+    e = 0.00000001
+    while (x - y) > e:
+        x = (x + y) / 2
+        y = n / x
+    return x
+
+
 # given the x**2-2=0 condition we bruteforce an aproximation to sqrt(2)
 
 
@@ -38,17 +39,22 @@ def sqrttest():
     last_aprox = None
     while error != 0:
         i = i + 1
-        error = (x ** target) - target
+        error = (x**target) - target
         x = x - (error * step)
 
         if error < lower_error:
             lower_error = abs(error)
 
-        print "Count: %d, Aprox: %s, Error: %s, LowerError: %s" % (
-            i,
-            x,
-            error,
-            lower_error,
+        print(
+            (
+                "Count: %d, Aprox: %s, Error: %s, LowerError: %s"
+                % (
+                    i,
+                    x,
+                    error,
+                    lower_error,
+                )
+            )
         )
 
 
diff --git a/subGCD.py b/subGCD.py
index b1d0ddd..3ec6f8e 100644
--- a/subGCD.py
+++ b/subGCD.py
@@ -20,9 +20,9 @@ def gcd(a, b):
 
 
 def test():
-    print(gcd(115, 5))
-    print(gcd(5, 7))
-    print(gcd(10, 4))
+    print((gcd(115, 5)))
+    print((gcd(5, 7)))
+    print((gcd(10, 4)))
 
     a = (
         37975227936943673922808872755445627854565536638199
@@ -32,7 +32,7 @@ def test():
         40094690950920881030683735292761468389214899724061
         * 5846418214406154678836553182979162384198610505601062333
     )
-    print(gcd(a, b))
+    print((gcd(a, b)))
 
 
 if __name__ == "__main__":
diff --git a/test_inv.py b/test_inv.py
index 52424f1..fff132c 100644
--- a/test_inv.py
+++ b/test_inv.py
@@ -1,26 +1,26 @@
 from gmpy2 import invert
 from timing import timing
 
+
 @timing
 def compute_modinv_1_n(n, p):
-  """
-  https://codeforces.com/blog/entry/83075
-  """
-  inv = [0, 1]
-  inv.extend((p - p // i) * inv[p % i] % p for i in range(2, n))
-  return inv
+    """
+    https://codeforces.com/blog/entry/83075
+    """
+    inv = [0, 1]
+    inv.extend((p - p // i) * inv[p % i] % p for i in range(2, n))
+    return inv
+
 
 @timing
 def compute_modinv_gmpy_1_n(n, p):
-  inv = [0]
-  inv.extend(int(invert(i, p)) for i in range(1, n))
-  return inv
+    inv = [0]
+    inv.extend(int(invert(i, p)) for i in range(1, n))
+    return inv
 
 
-n, p = 10**4, 10**9+7
-#n, p = 10, 65537
+n, p = 10**4, 10**9 + 7
+# n, p = 10, 65537
 
 a = compute_modinv_1_n(n, p)
 b = compute_modinv_gmpy_1_n(n, p)
-
-
diff --git a/timing.py b/timing.py
index 751cf1d..18fc78d 100644
--- a/timing.py
+++ b/timing.py
@@ -1,12 +1,16 @@
 from functools import wraps
 from time import time
+
+
 def timing(f):
     @wraps(f)
     def wrap(*args, **kw):
         ts = time()
         result = f(*args, **kw)
         te = time()
-        print ('func:%r args:[%r, %r] took: %2.4f sec' % \
-          (f.__name__, args, kw, te-ts))
+        print(
+            ("func:%r args:[%r, %r] took: %2.4f sec" % (f.__name__, args, kw, te - ts))
+        )
         return result
+
     return wrap
diff --git a/vectors.py b/vectors.py
index 3f07592..4742707 100644
--- a/vectors.py
+++ b/vectors.py
@@ -60,22 +60,22 @@ def magnitude(vec):
 
 def tests():
     vec = [0, 4, -3]
-    print length(vec)
-    print normalize(vec)
+    print((length(vec)))
+    print((normalize(vec)))
 
     s = 3
     k = [1, 2]
     j = [2, 3]
 
     tmp = multiply(k, j)
-    print tmp
-    print multiplyScalar(tmp, s)
+    print(tmp)
+    print((multiplyScalar(tmp, s)))
 
     k = [0, 1, 0]
     j = [1, 0, 0]
 
-    print dot(k, j)
-    print cross(k, j)
+    print((dot(k, j)))
+    print((cross(k, j)))
 
-    print magnitude(j)
-    print magnitude(k)
+    print((magnitude(j)))
+    print((magnitude(k)))
diff --git a/wilsons_theorem.py b/wilsons_theorem.py
index 58c6499..7d3a88e 100644
--- a/wilsons_theorem.py
+++ b/wilsons_theorem.py
@@ -6,4 +6,4 @@ def wilson_is_prime(n):
 
 
 for i in range(1, 100):
-    print(i, wilson_is_prime(i), is_prime(i))
+    print((i, wilson_is_prime(i), is_prime(i)))