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
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
### 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)