mdw@git.distorted.org.uk Git - catacomb/blob - utils/curve25519.sage

   1 #! /usr/local/bin/sage
   2 ### -*- mode: python; coding: utf-8 -*-
   3
   4 ###--------------------------------------------------------------------------
   5 ### Some general utilities.
   6
   7 def ld(v):
   8   return 0 + sum(ord(v[i]) << 8*i for i in xrange(len(v)))
   9
  10 def st(x, n):
  11   return ''.join(chr((x >> 8*i)&0xff) for i in xrange(n))
  12
  13 ###--------------------------------------------------------------------------
  14 ### Define the curve.
  15
  16 p = 2^255 - 19; k = GF(p)
  17 A = k(486662); A0 = (A - 2)/4
  18 E = EllipticCurve(k, [0, A, 0, 1, 0]); P = E.lift_x(9)
  19 l = 2^252 + 27742317777372353535851937790883648493
  20
  21 assert is_prime(l)
  22 assert (l*P).is_zero()
  23 assert (p + 1 - 8*l)^2 <= 4*p
  24
  25 ###--------------------------------------------------------------------------
  26 ### Example points from `Cryptography in NaCl'.
  27
  28 x = ld(map(chr, [0x70,0x07,0x6d,0x0a,0x73,0x18,0xa5,0x7d
  29 ,0x3c,0x16,0xc1,0x72,0x51,0xb2,0x66,0x45
  30 ,0xdf,0x4c,0x2f,0x87,0xeb,0xc0,0x99,0x2a
  31 ,0xb1,0x77,0xfb,0xa5,0x1d,0xb9,0x2c,0x6a]))
  32 y = ld(map(chr, [0x58,0xab,0x08,0x7e,0x62,0x4a,0x8a,0x4b
  33 ,0x79,0xe1,0x7f,0x8b,0x83,0x80,0x0e,0xe6
  34 ,0x6f,0x3b,0xb1,0x29,0x26,0x18,0xb6,0xfd
  35 ,0x1c,0x2f,0x8b,0x27,0xff,0x88,0xe0,0x6b]))
  36 X = x*P
  37 Y = y*P
  38 Z = x*Y
  39 assert Z == y*X
  40
  41 ###--------------------------------------------------------------------------
  42 ### Arithmetic implementation.
  43
  44 def sqrn(x, n):
  45   for i in xrange(n): x = x*x
  46   return x
  47
  48 def inv(x):
  49   t2 = sqrn(x, 1)        #   1 | 2
  50   u = sqrn(t2, 2)        #   3 | 8
  51   t = u*x                #   4 | 9
  52   t11 = t*t2             #   5 | 11
  53   u = sqrn(t11, 1)       #   6 | 22
  54   t = u*t                #   7 | 2^5 - 1 = 31
  55   u = sqrn(t, 5)         #  12 | 2^10 - 2^5
  56   t2p10m1 = u*t          #  13 | 2^10 - 1
  57   u = sqrn(t2p10m1, 10)  #  23 | 2^20 - 2^10
  58   t = u*t2p10m1          #  24 | 2^20 - 1
  59   u = sqrn(t, 20)        #  44 | 2^40 - 2^20
  60   t = u*t                #  45 | 2^40 - 1
  61   u = sqrn(t, 10)        #  55 | 2^50 - 2^10
  62   t2p50m1 = u*t2p10m1    #  56 | 2^50 - 1
  63   u = sqrn(t2p50m1, 50)  # 106 | 2^100 - 2^50
  64   t = u*t2p50m1          # 107 | 2^100 - 1
  65   u = sqrn(t, 100)       # 207 | 2^200 - 2^100
  66   t = u*t                # 208 | 2^200 - 1
  67   u = sqrn(t, 50)        # 258 | 2^250 - 2^50
  68   t = u*t2p50m1          # 259 | 2^250 - 1
  69   u = sqrn(t, 5)         # 264 | 2^255 - 2^5
  70   t = u*t11              # 265 | 2^255 - 21
  71   return t
  72
  73 assert inv(k(9))*9 == 1
  74
  75 ###--------------------------------------------------------------------------
  76 ### The Montgomery ladder.
  77
  78 A0 = (A - 2)/4
  79
  80 def x25519(n, x1):
  81
  82   ## Let Q = (x_1 : y_1 : 1) be an input point.  We calculate
  83   ## n Q = (x_n : y_n : z_n), returning x_n/z_n (unless z_n = 0,
  84   ## in which case we return zero).
  85   ##
  86   ## We're given that n = 2^254 + n'_254, where 0 <= n'_254 < 2^254.
  87   bb = n.bits()
  88   x, z = 1, 0
  89   u, w = x1, 1
  90
  91   ## Initially, let i = 255.
  92   for i in xrange(len(bb) - 1, -1, -1):
  93
  94     ## Split n = n_i 2^i + n'_i, where 0 <= n'_i < 2^i, so n_0 = n.
  95     ## We have x, z = x_{n_{i+1}}, z_{n_{i+1}}, and
  96     ## u, w = x_{n_{i+1}+1}, z_{n_{i+1}+1}.
  97     ## Now either n_i = 2 n_{i+1} or n_i = 2 n_{i+1} + 1, depending
  98     ## on bit i of n.
  99
 100     ## Swap (x : z) and (u : w) if bit i of n is set.
 101     if bb[i]: x, z, u, w = u, w, x, z
 102
 103     ## Do the ladder step.
 104     xmz, xpz = x - z, x + z
 105     umw, upw = u - w, u + w
 106     xmz2, xpz2 = xmz^2, xpz^2
 107     xpz2mxmz2 = xpz2 - xmz2
 108     xmzupw, xpzumw = xmz*upw, xpz*umw
 109     x, z = xmz2*xpz2, xpz2mxmz2*(xpz2 + A0*xpz2mxmz2)
 110     u, w = (xmzupw + xpzumw)^2, x1*(xmzupw - xpzumw)^2
 111
 112     ## Finally, unswap.
 113     if bb[i]: x, z, u, w = u, w, x, z
 114
 115   ## Almost done.
 116   return x*inv(z)
 117
 118 assert x25519(y, k(9)) == Y[0]
 119 assert x25519(x, Y[0]) == x25519(y, X[0]) == Z[0]
 120
 121 ###----- That's all, folks --------------------------------------------------