Use extended Euclidean algorithm for modular inverse

ed25519.py

@@ -39,22 +39,20 @@ def pow2(x, p):
         p -= 1
     return x
 
-
 def inv(z):
-    """$= z^{-1} \mod q$, for z != 0"""
-    # Adapted from curve25519_athlon.c in djb's Curve25519.
-    z2 = z * z % q                                # 2
-    z9 = pow2(z2, 2) * z % q                      # 9
-    z11 = z9 * z2 % q                             # 11
-    z2_5_0 = (z11 * z11) % q * z9 % q             # 31 == 2^5 - 2^0
-    z2_10_0 = pow2(z2_5_0, 5) * z2_5_0 % q        # 2^10 - 2^0
-    z2_20_0 = pow2(z2_10_0, 10) * z2_10_0 % q     # ...
-    z2_40_0 = pow2(z2_20_0, 20) * z2_20_0 % q
-    z2_50_0 = pow2(z2_40_0, 10) * z2_10_0 % q
-    z2_100_0 = pow2(z2_50_0, 50) * z2_50_0 % q
-    z2_200_0 = pow2(z2_100_0, 100) * z2_100_0 % q
-    z2_250_0 = pow2(z2_200_0, 50) * z2_50_0 % q   # 2^250 - 2^0
-    return pow2(z2_250_0, 5) * z11 % q            # 2^255 - 2^5 + 11 = q - 2
+    """
+    Use the extended Euclidean algorithm to find the inverse mod q.
+
+    See, for example, H. Cohen, A Course in Computational Algebraic Number
+    Theory, Algorithm 1.3.6
+    """
+    d, div = z, q
+    u = 1
+    r = 0
+    while div != 0:
+        (qq, div), d = divmod(d, div), div
+        r, u = u - qq * r, r
+    return u
 
 
 d = -121665 * inv(121666)
@@ -92,7 +90,7 @@ def scalarmult(P, e):
     if e == 0:
         return (0, 1)
 
-    Q = scalarmult(P, e // 2)
+    Q = scalarmult(P, e / 2)
     Q = edwards(Q, Q)
 
     if e & 1: