X-Git-Url: https://git.distorted.org.uk/~mdw/catacomb/blobdiff_plain/9c1437f372e62f0b3b3a7162aabee73bdc96ce4b..a90d420cbe87490c844ae422c966e746d3134b07:/utils/curve25519.sage

diff --git a/utils/curve25519.sage b/utils/curve25519.sage
index 416f059f..6fa3cd34 100644
--- a/utils/curve25519.sage
+++ b/utils/curve25519.sage
@@ -115,7 +115,7 @@ def inv(x):
   t = u*t11              # 265 | 2^255 - 21
   return t
 
-def quosqrt(x, y):
+def quosqrt_djb(x, y):
 
   ## First, some preliminary values.
   y2 = sqrn(y, 1)        #   1 | 0, 2
@@ -162,14 +162,66 @@ def quosqrt(x, y):
   elif t == -x: return beta*sqrtm1
   else: raise ValueError, 'not a square'
 
+def quosqrt_mdw(x, y):
+  v = x*y
+
+  ## Now we calculate w = v^{3*2^252 - 8}.  This will be explained later.
+  u = sqrn(v, 1)         #   1 | 2
+  t = u*v                #   2 | 3
+  u = sqrn(t, 2)         #   4 | 12
+  t15 = u*t              #   5 | 15
+  u = sqrn(t15, 1)       #   6 | 30
+  t = u*v                #   7 | 31 = 2^5 - 1
+  u = sqrn(t, 5)         #  12 | 2^10 - 2^5
+  t = u*t                #  13 | 2^10 - 1
+  u = sqrn(t, 10)        #  23 | 2^20 - 2^10
+  u = u*t                #  24 | 2^20 - 1
+  u = sqrn(u, 10)        #  34 | 2^30 - 2^10
+  t = u*t                #  35 | 2^30 - 1
+  u = sqrn(t, 1)         #  36 | 2^31 - 2
+  t = u*v                #  37 | 2^31 - 1
+  u = sqrn(t, 31)        #  68 | 2^62 - 2^31
+  t = u*t                #  69 | 2^62 - 1
+  u = sqrn(t, 62)        # 131 | 2^124 - 2^62
+  t = u*t                # 132 | 2^124 - 1
+  u = sqrn(t, 124)       # 256 | 2^248 - 2^124
+  t = u*t                # 257 | 2^248 - 1
+  u = sqrn(t, 1)         # 258 | 2^249 - 2
+  t = u*v                # 259 | 2^249 - 1
+  t = sqrn(t, 3)         # 262 | 2^252 - 8
+  u = sqrn(t, 1)         # 263 | 2^253 - 16
+  t = u*t                # 264 | 3*2^252 - 24
+  t = t*t15              # 265 | 3*2^252 - 9
+  w = t*v                # 266 | 3*2^252 - 8
+
+  ## Awesome.  Now let me explain.  Let v be a square in GF(p), and let w =
+  ## v^(3*2^252 - 8).  In particular, let's consider
+  ##
+  ##    v^2 w^4 = v^2 v^{3*2^254 - 32} = (v^{2^254 - 10})^3
+  ##
+  ## But 2^254 - 10 = ((2^255 - 19) - 1)/2 = (p - 1)/2.  Since v is a square,
+  ## it has order dividing (p - 1)/2, and therefore v^2 w^4 = 1 and
+  ##
+  ##    w^4 = 1/v^2
+  ##
+  ## That in turn implies that w^2 = ±1/v.  Now, recall that v = x y, and let
+  ## w' = w x.  Then w'^2 = ±x^2/v = ±x/y.  If y w'^2 = x then we set
+  ## z = w', since we have z^2 = x/y; otherwise let z = i w', where i^2 = -1,
+  ## so z^2 = -w^2 = x/y, and we're done.
+  t = w*x
+  u = y*t^2
+  if u == x: return t
+  elif u == -x: return t*sqrtm1
+  else: raise ValueError, 'not a square'
+
+quosqrt = quosqrt_mdw
+
 assert inv(k(9))*9 == 1
 assert 5*quosqrt(k(4), k(5))^2 == 4
 
 ###--------------------------------------------------------------------------
 ### The Montgomery ladder.
 
-A0 = (A - 2)/4
-
 def x25519(n, x1):
 
   ## Let Q = (x_1 : y_1 : 1) be an input point.  We calculate
@@ -215,7 +267,7 @@ assert x25519(x, Y[0]) == x25519(y, X[0]) == Z[0]
 ### Edwards curve parameters and conversion.
 
 a = k(-1)
-d = k(-A0/(A0 + 1))
+d = -A0/(A0 + 1)
 
 def mont_to_ed(u, v):
   return sqrt(-A - 2)*u/v, (u - 1)/(u + 1)