crypto/ed25519/internal/edwards25519/edwards25519.go

   1 package edwards25519
   2
   3 // This code is a port of the public domain, “ref10” implementation of ed25519
   4 // from SUPERCOP.
   5
   6 // FieldElement represents an element of the field GF(2^255 - 19).  An element
   7 // t, entries t[0]...t[9], represents the integer t[0]+2^26 t[1]+2^51 t[2]+2^77
   8 // t[3]+2^102 t[4]+...+2^230 t[9].  Bounds on each t[i] vary depending on
   9 // context.
  10 type FieldElement [10]int32
  11
  12 var zero FieldElement
  13
  14 func FeZero(fe *FieldElement) {
  15         copy(fe[:], zero[:])
  16 }
  17
  18 func FeOne(fe *FieldElement) {
  19         FeZero(fe)
  20         fe[0] = 1
  21 }
  22
  23 func FeAdd(dst, a, b *FieldElement) {
  24         dst[0] = a[0] + b[0]
  25         dst[1] = a[1] + b[1]
  26         dst[2] = a[2] + b[2]
  27         dst[3] = a[3] + b[3]
  28         dst[4] = a[4] + b[4]
  29         dst[5] = a[5] + b[5]
  30         dst[6] = a[6] + b[6]
  31         dst[7] = a[7] + b[7]
  32         dst[8] = a[8] + b[8]
  33         dst[9] = a[9] + b[9]
  34 }
  35
  36 func FeSub(dst, a, b *FieldElement) {
  37         dst[0] = a[0] - b[0]
  38         dst[1] = a[1] - b[1]
  39         dst[2] = a[2] - b[2]
  40         dst[3] = a[3] - b[3]
  41         dst[4] = a[4] - b[4]
  42         dst[5] = a[5] - b[5]
  43         dst[6] = a[6] - b[6]
  44         dst[7] = a[7] - b[7]
  45         dst[8] = a[8] - b[8]
  46         dst[9] = a[9] - b[9]
  47 }
  48
  49 func FeCopy(dst, src *FieldElement) {
  50         copy(dst[:], src[:])
  51 }
  52
  53 // Replace (f,g) with (g,g) if b == 1;
  54 // replace (f,g) with (f,g) if b == 0.
  55 //
  56 // Preconditions: b in {0,1}.
  57 func FeCMove(f, g *FieldElement, b int32) {
  58         b = -b
  59         f[0] ^= b & (f[0] ^ g[0])
  60         f[1] ^= b & (f[1] ^ g[1])
  61         f[2] ^= b & (f[2] ^ g[2])
  62         f[3] ^= b & (f[3] ^ g[3])
  63         f[4] ^= b & (f[4] ^ g[4])
  64         f[5] ^= b & (f[5] ^ g[5])
  65         f[6] ^= b & (f[6] ^ g[6])
  66         f[7] ^= b & (f[7] ^ g[7])
  67         f[8] ^= b & (f[8] ^ g[8])
  68         f[9] ^= b & (f[9] ^ g[9])
  69 }
  70
  71 func load3(in []byte) int64 {
  72         var r int64
  73         r = int64(in[0])
  74         r |= int64(in[1]) << 8
  75         r |= int64(in[2]) << 16
  76         return r
  77 }
  78
  79 func load4(in []byte) int64 {
  80         var r int64
  81         r = int64(in[0])
  82         r |= int64(in[1]) << 8
  83         r |= int64(in[2]) << 16
  84         r |= int64(in[3]) << 24
  85         return r
  86 }
  87
  88 func FeFromBytes(dst *FieldElement, src *[32]byte) {
  89         h0 := load4(src[:])
  90         h1 := load3(src[4:]) << 6
  91         h2 := load3(src[7:]) << 5
  92         h3 := load3(src[10:]) << 3
  93         h4 := load3(src[13:]) << 2
  94         h5 := load4(src[16:])
  95         h6 := load3(src[20:]) << 7
  96         h7 := load3(src[23:]) << 5
  97         h8 := load3(src[26:]) << 4
  98         h9 := (load3(src[29:]) & 8388607) << 2
  99
 100         FeCombine(dst, h0, h1, h2, h3, h4, h5, h6, h7, h8, h9)
 101 }
 102
 103 // FeToBytes marshals h to s.
 104 // Preconditions:
 105 //   |h| bounded by 1.1*2^25,1.1*2^24,1.1*2^25,1.1*2^24,etc.
 106 //
 107 // Write p=2^255-19; q=floor(h/p).
 108 // Basic claim: q = floor(2^(-255)(h + 19 2^(-25)h9 + 2^(-1))).
 109 //
 110 // Proof:
 111 //   Have |h|<=p so |q|<=1 so |19^2 2^(-255) q|<1/4.
 112 //   Also have |h-2^230 h9|<2^230 so |19 2^(-255)(h-2^230 h9)|<1/4.
 113 //
 114 //   Write y=2^(-1)-19^2 2^(-255)q-19 2^(-255)(h-2^230 h9).
 115 //   Then 0<y<1.
 116 //
 117 //   Write r=h-pq.
 118 //   Have 0<=r<=p-1=2^255-20.
 119 //   Thus 0<=r+19(2^-255)r<r+19(2^-255)2^255<=2^255-1.
 120 //
 121 //   Write x=r+19(2^-255)r+y.
 122 //   Then 0<x<2^255 so floor(2^(-255)x) = 0 so floor(q+2^(-255)x) = q.
 123 //
 124 //   Have q+2^(-255)x = 2^(-255)(h + 19 2^(-25) h9 + 2^(-1))
 125 //   so floor(2^(-255)(h + 19 2^(-25) h9 + 2^(-1))) = q.
 126 func FeToBytes(s *[32]byte, h *FieldElement) {
 127         var carry [10]int32
 128
 129         q := (19*h[9] + (1 << 24)) >> 25
 130         q = (h[0] + q) >> 26
 131         q = (h[1] + q) >> 25
 132         q = (h[2] + q) >> 26
 133         q = (h[3] + q) >> 25
 134         q = (h[4] + q) >> 26
 135         q = (h[5] + q) >> 25
 136         q = (h[6] + q) >> 26
 137         q = (h[7] + q) >> 25
 138         q = (h[8] + q) >> 26
 139         q = (h[9] + q) >> 25
 140
 141         // Goal: Output h-(2^255-19)q, which is between 0 and 2^255-20.
 142         h[0] += 19 * q
 143         // Goal: Output h-2^255 q, which is between 0 and 2^255-20.
 144
 145         carry[0] = h[0] >> 26
 146         h[1] += carry[0]
 147         h[0] -= carry[0] << 26
 148         carry[1] = h[1] >> 25
 149         h[2] += carry[1]
 150         h[1] -= carry[1] << 25
 151         carry[2] = h[2] >> 26
 152         h[3] += carry[2]
 153         h[2] -= carry[2] << 26
 154         carry[3] = h[3] >> 25
 155         h[4] += carry[3]
 156         h[3] -= carry[3] << 25
 157         carry[4] = h[4] >> 26
 158         h[5] += carry[4]
 159         h[4] -= carry[4] << 26
 160         carry[5] = h[5] >> 25
 161         h[6] += carry[5]
 162         h[5] -= carry[5] << 25
 163         carry[6] = h[6] >> 26
 164         h[7] += carry[6]
 165         h[6] -= carry[6] << 26
 166         carry[7] = h[7] >> 25
 167         h[8] += carry[7]
 168         h[7] -= carry[7] << 25
 169         carry[8] = h[8] >> 26
 170         h[9] += carry[8]
 171         h[8] -= carry[8] << 26
 172         carry[9] = h[9] >> 25
 173         h[9] -= carry[9] << 25
 174         // h10 = carry9
 175
 176         // Goal: Output h[0]+...+2^255 h10-2^255 q, which is between 0 and 2^255-20.
 177         // Have h[0]+...+2^230 h[9] between 0 and 2^255-1;
 178         // evidently 2^255 h10-2^255 q = 0.
 179         // Goal: Output h[0]+...+2^230 h[9].
 180
 181         s[0] = byte(h[0] >> 0)
 182         s[1] = byte(h[0] >> 8)
 183         s[2] = byte(h[0] >> 16)
 184         s[3] = byte((h[0] >> 24) | (h[1] << 2))
 185         s[4] = byte(h[1] >> 6)
 186         s[5] = byte(h[1] >> 14)
 187         s[6] = byte((h[1] >> 22) | (h[2] << 3))
 188         s[7] = byte(h[2] >> 5)
 189         s[8] = byte(h[2] >> 13)
 190         s[9] = byte((h[2] >> 21) | (h[3] << 5))
 191         s[10] = byte(h[3] >> 3)
 192         s[11] = byte(h[3] >> 11)
 193         s[12] = byte((h[3] >> 19) | (h[4] << 6))
 194         s[13] = byte(h[4] >> 2)
 195         s[14] = byte(h[4] >> 10)
 196         s[15] = byte(h[4] >> 18)
 197         s[16] = byte(h[5] >> 0)
 198         s[17] = byte(h[5] >> 8)
 199         s[18] = byte(h[5] >> 16)
 200         s[19] = byte((h[5] >> 24) | (h[6] << 1))
 201         s[20] = byte(h[6] >> 7)
 202         s[21] = byte(h[6] >> 15)
 203         s[22] = byte((h[6] >> 23) | (h[7] << 3))
 204         s[23] = byte(h[7] >> 5)
 205         s[24] = byte(h[7] >> 13)
 206         s[25] = byte((h[7] >> 21) | (h[8] << 4))
 207         s[26] = byte(h[8] >> 4)
 208         s[27] = byte(h[8] >> 12)
 209         s[28] = byte((h[8] >> 20) | (h[9] << 6))
 210         s[29] = byte(h[9] >> 2)
 211         s[30] = byte(h[9] >> 10)
 212         s[31] = byte(h[9] >> 18)
 213 }
 214
 215 func FeIsNegative(f *FieldElement) byte {
 216         var s [32]byte
 217         FeToBytes(&s, f)
 218         return s[0] & 1
 219 }
 220
 221 func FeIsNonZero(f *FieldElement) int32 {
 222         var s [32]byte
 223         FeToBytes(&s, f)
 224         var x uint8
 225         for _, b := range s {
 226                 x |= b
 227         }
 228         x |= x >> 4
 229         x |= x >> 2
 230         x |= x >> 1
 231         return int32(x & 1)
 232 }
 233
 234 // FeNeg sets h = -f
 235 //
 236 // Preconditions:
 237 //    |f| bounded by 1.1*2^25,1.1*2^24,1.1*2^25,1.1*2^24,etc.
 238 //
 239 // Postconditions:
 240 //    |h| bounded by 1.1*2^25,1.1*2^24,1.1*2^25,1.1*2^24,etc.
 241 func FeNeg(h, f *FieldElement) {
 242         h[0] = -f[0]
 243         h[1] = -f[1]
 244         h[2] = -f[2]
 245         h[3] = -f[3]
 246         h[4] = -f[4]
 247         h[5] = -f[5]
 248         h[6] = -f[6]
 249         h[7] = -f[7]
 250         h[8] = -f[8]
 251         h[9] = -f[9]
 252 }
 253
 254 func FeCombine(h *FieldElement, h0, h1, h2, h3, h4, h5, h6, h7, h8, h9 int64) {
 255         var c0, c1, c2, c3, c4, c5, c6, c7, c8, c9 int64
 256
 257         /*
 258           |h0| <= (1.1*1.1*2^52*(1+19+19+19+19)+1.1*1.1*2^50*(38+38+38+38+38))
 259             i.e. |h0| <= 1.2*2^59; narrower ranges for h2, h4, h6, h8
 260           |h1| <= (1.1*1.1*2^51*(1+1+19+19+19+19+19+19+19+19))
 261             i.e. |h1| <= 1.5*2^58; narrower ranges for h3, h5, h7, h9
 262         */
 263
 264         c0 = (h0 + (1 << 25)) >> 26
 265         h1 += c0
 266         h0 -= c0 << 26
 267         c4 = (h4 + (1 << 25)) >> 26
 268         h5 += c4
 269         h4 -= c4 << 26
 270         /* |h0| <= 2^25 */
 271         /* |h4| <= 2^25 */
 272         /* |h1| <= 1.51*2^58 */
 273         /* |h5| <= 1.51*2^58 */
 274
 275         c1 = (h1 + (1 << 24)) >> 25
 276         h2 += c1
 277         h1 -= c1 << 25
 278         c5 = (h5 + (1 << 24)) >> 25
 279         h6 += c5
 280         h5 -= c5 << 25
 281         /* |h1| <= 2^24; from now on fits into int32 */
 282         /* |h5| <= 2^24; from now on fits into int32 */
 283         /* |h2| <= 1.21*2^59 */
 284         /* |h6| <= 1.21*2^59 */
 285
 286         c2 = (h2 + (1 << 25)) >> 26
 287         h3 += c2
 288         h2 -= c2 << 26
 289         c6 = (h6 + (1 << 25)) >> 26
 290         h7 += c6
 291         h6 -= c6 << 26
 292         /* |h2| <= 2^25; from now on fits into int32 unchanged */
 293         /* |h6| <= 2^25; from now on fits into int32 unchanged */
 294         /* |h3| <= 1.51*2^58 */
 295         /* |h7| <= 1.51*2^58 */
 296
 297         c3 = (h3 + (1 << 24)) >> 25
 298         h4 += c3
 299         h3 -= c3 << 25
 300         c7 = (h7 + (1 << 24)) >> 25
 301         h8 += c7
 302         h7 -= c7 << 25
 303         /* |h3| <= 2^24; from now on fits into int32 unchanged */
 304         /* |h7| <= 2^24; from now on fits into int32 unchanged */
 305         /* |h4| <= 1.52*2^33 */
 306         /* |h8| <= 1.52*2^33 */
 307
 308         c4 = (h4 + (1 << 25)) >> 26
 309         h5 += c4
 310         h4 -= c4 << 26
 311         c8 = (h8 + (1 << 25)) >> 26
 312         h9 += c8
 313         h8 -= c8 << 26
 314         /* |h4| <= 2^25; from now on fits into int32 unchanged */
 315         /* |h8| <= 2^25; from now on fits into int32 unchanged */
 316         /* |h5| <= 1.01*2^24 */
 317         /* |h9| <= 1.51*2^58 */
 318
 319         c9 = (h9 + (1 << 24)) >> 25
 320         h0 += c9 * 19
 321         h9 -= c9 << 25
 322         /* |h9| <= 2^24; from now on fits into int32 unchanged */
 323         /* |h0| <= 1.8*2^37 */
 324
 325         c0 = (h0 + (1 << 25)) >> 26
 326         h1 += c0
 327         h0 -= c0 << 26
 328         /* |h0| <= 2^25; from now on fits into int32 unchanged */
 329         /* |h1| <= 1.01*2^24 */
 330
 331         h[0] = int32(h0)
 332         h[1] = int32(h1)
 333         h[2] = int32(h2)
 334         h[3] = int32(h3)
 335         h[4] = int32(h4)
 336         h[5] = int32(h5)
 337         h[6] = int32(h6)
 338         h[7] = int32(h7)
 339         h[8] = int32(h8)
 340         h[9] = int32(h9)
 341 }
 342
 343 // FeMul calculates h = f * g
 344 // Can overlap h with f or g.
 345 //
 346 // Preconditions:
 347 //    |f| bounded by 1.1*2^26,1.1*2^25,1.1*2^26,1.1*2^25,etc.
 348 //    |g| bounded by 1.1*2^26,1.1*2^25,1.1*2^26,1.1*2^25,etc.
 349 //
 350 // Postconditions:
 351 //    |h| bounded by 1.1*2^25,1.1*2^24,1.1*2^25,1.1*2^24,etc.
 352 //
 353 // Notes on implementation strategy:
 354 //
 355 // Using schoolbook multiplication.
 356 // Karatsuba would save a little in some cost models.
 357 //
 358 // Most multiplications by 2 and 19 are 32-bit precomputations;
 359 // cheaper than 64-bit postcomputations.
 360 //
 361 // There is one remaining multiplication by 19 in the carry chain;
 362 // one *19 precomputation can be merged into this,
 363 // but the resulting data flow is considerably less clean.
 364 //
 365 // There are 12 carries below.
 366 // 10 of them are 2-way parallelizable and vectorizable.
 367 // Can get away with 11 carries, but then data flow is much deeper.
 368 //
 369 // With tighter constraints on inputs, can squeeze carries into int32.
 370 func FeMul(h, f, g *FieldElement) {
 371         f0 := int64(f[0])
 372         f1 := int64(f[1])
 373         f2 := int64(f[2])
 374         f3 := int64(f[3])
 375         f4 := int64(f[4])
 376         f5 := int64(f[5])
 377         f6 := int64(f[6])
 378         f7 := int64(f[7])
 379         f8 := int64(f[8])
 380         f9 := int64(f[9])
 381
 382         f1_2 := int64(2 * f[1])
 383         f3_2 := int64(2 * f[3])
 384         f5_2 := int64(2 * f[5])
 385         f7_2 := int64(2 * f[7])
 386         f9_2 := int64(2 * f[9])
 387
 388         g0 := int64(g[0])
 389         g1 := int64(g[1])
 390         g2 := int64(g[2])
 391         g3 := int64(g[3])
 392         g4 := int64(g[4])
 393         g5 := int64(g[5])
 394         g6 := int64(g[6])
 395         g7 := int64(g[7])
 396         g8 := int64(g[8])
 397         g9 := int64(g[9])
 398
 399         g1_19 := int64(19 * g[1]) /* 1.4*2^29 */
 400         g2_19 := int64(19 * g[2]) /* 1.4*2^30; still ok */
 401         g3_19 := int64(19 * g[3])
 402         g4_19 := int64(19 * g[4])
 403         g5_19 := int64(19 * g[5])
 404         g6_19 := int64(19 * g[6])
 405         g7_19 := int64(19 * g[7])
 406         g8_19 := int64(19 * g[8])
 407         g9_19 := int64(19 * g[9])
 408
 409         h0 := f0*g0 + f1_2*g9_19 + f2*g8_19 + f3_2*g7_19 + f4*g6_19 + f5_2*g5_19 + f6*g4_19 + f7_2*g3_19 + f8*g2_19 + f9_2*g1_19
 410         h1 := f0*g1 + f1*g0 + f2*g9_19 + f3*g8_19 + f4*g7_19 + f5*g6_19 + f6*g5_19 + f7*g4_19 + f8*g3_19 + f9*g2_19
 411         h2 := f0*g2 + f1_2*g1 + f2*g0 + f3_2*g9_19 + f4*g8_19 + f5_2*g7_19 + f6*g6_19 + f7_2*g5_19 + f8*g4_19 + f9_2*g3_19
 412         h3 := f0*g3 + f1*g2 + f2*g1 + f3*g0 + f4*g9_19 + f5*g8_19 + f6*g7_19 + f7*g6_19 + f8*g5_19 + f9*g4_19
 413         h4 := f0*g4 + f1_2*g3 + f2*g2 + f3_2*g1 + f4*g0 + f5_2*g9_19 + f6*g8_19 + f7_2*g7_19 + f8*g6_19 + f9_2*g5_19
 414         h5 := f0*g5 + f1*g4 + f2*g3 + f3*g2 + f4*g1 + f5*g0 + f6*g9_19 + f7*g8_19 + f8*g7_19 + f9*g6_19
 415         h6 := f0*g6 + f1_2*g5 + f2*g4 + f3_2*g3 + f4*g2 + f5_2*g1 + f6*g0 + f7_2*g9_19 + f8*g8_19 + f9_2*g7_19
 416         h7 := f0*g7 + f1*g6 + f2*g5 + f3*g4 + f4*g3 + f5*g2 + f6*g1 + f7*g0 + f8*g9_19 + f9*g8_19
 417         h8 := f0*g8 + f1_2*g7 + f2*g6 + f3_2*g5 + f4*g4 + f5_2*g3 + f6*g2 + f7_2*g1 + f8*g0 + f9_2*g9_19
 418         h9 := f0*g9 + f1*g8 + f2*g7 + f3*g6 + f4*g5 + f5*g4 + f6*g3 + f7*g2 + f8*g1 + f9*g0
 419
 420         FeCombine(h, h0, h1, h2, h3, h4, h5, h6, h7, h8, h9)
 421 }
 422
 423 func feSquare(f *FieldElement) (h0, h1, h2, h3, h4, h5, h6, h7, h8, h9 int64) {
 424         f0 := int64(f[0])
 425         f1 := int64(f[1])
 426         f2 := int64(f[2])
 427         f3 := int64(f[3])
 428         f4 := int64(f[4])
 429         f5 := int64(f[5])
 430         f6 := int64(f[6])
 431         f7 := int64(f[7])
 432         f8 := int64(f[8])
 433         f9 := int64(f[9])
 434         f0_2 := int64(2 * f[0])
 435         f1_2 := int64(2 * f[1])
 436         f2_2 := int64(2 * f[2])
 437         f3_2 := int64(2 * f[3])
 438         f4_2 := int64(2 * f[4])
 439         f5_2 := int64(2 * f[5])
 440         f6_2 := int64(2 * f[6])
 441         f7_2 := int64(2 * f[7])
 442         f5_38 := 38 * f5 // 1.31*2^30
 443         f6_19 := 19 * f6 // 1.31*2^30
 444         f7_38 := 38 * f7 // 1.31*2^30
 445         f8_19 := 19 * f8 // 1.31*2^30
 446         f9_38 := 38 * f9 // 1.31*2^30
 447
 448         h0 = f0*f0 + f1_2*f9_38 + f2_2*f8_19 + f3_2*f7_38 + f4_2*f6_19 + f5*f5_38
 449         h1 = f0_2*f1 + f2*f9_38 + f3_2*f8_19 + f4*f7_38 + f5_2*f6_19
 450         h2 = f0_2*f2 + f1_2*f1 + f3_2*f9_38 + f4_2*f8_19 + f5_2*f7_38 + f6*f6_19
 451         h3 = f0_2*f3 + f1_2*f2 + f4*f9_38 + f5_2*f8_19 + f6*f7_38
 452         h4 = f0_2*f4 + f1_2*f3_2 + f2*f2 + f5_2*f9_38 + f6_2*f8_19 + f7*f7_38
 453         h5 = f0_2*f5 + f1_2*f4 + f2_2*f3 + f6*f9_38 + f7_2*f8_19
 454         h6 = f0_2*f6 + f1_2*f5_2 + f2_2*f4 + f3_2*f3 + f7_2*f9_38 + f8*f8_19
 455         h7 = f0_2*f7 + f1_2*f6 + f2_2*f5 + f3_2*f4 + f8*f9_38
 456         h8 = f0_2*f8 + f1_2*f7_2 + f2_2*f6 + f3_2*f5_2 + f4*f4 + f9*f9_38
 457         h9 = f0_2*f9 + f1_2*f8 + f2_2*f7 + f3_2*f6 + f4_2*f5
 458
 459         return
 460 }
 461
 462 // FeSquare calculates h = f*f. Can overlap h with f.
 463 //
 464 // Preconditions:
 465 //    |f| bounded by 1.1*2^26,1.1*2^25,1.1*2^26,1.1*2^25,etc.
 466 //
 467 // Postconditions:
 468 //    |h| bounded by 1.1*2^25,1.1*2^24,1.1*2^25,1.1*2^24,etc.
 469 func FeSquare(h, f *FieldElement) {
 470         h0, h1, h2, h3, h4, h5, h6, h7, h8, h9 := feSquare(f)
 471         FeCombine(h, h0, h1, h2, h3, h4, h5, h6, h7, h8, h9)
 472 }
 473
 474 // FeSquare2 sets h = 2 * f * f
 475 //
 476 // Can overlap h with f.
 477 //
 478 // Preconditions:
 479 //    |f| bounded by 1.65*2^26,1.65*2^25,1.65*2^26,1.65*2^25,etc.
 480 //
 481 // Postconditions:
 482 //    |h| bounded by 1.01*2^25,1.01*2^24,1.01*2^25,1.01*2^24,etc.
 483 // See fe_mul.c for discussion of implementation strategy.
 484 func FeSquare2(h, f *FieldElement) {
 485         h0, h1, h2, h3, h4, h5, h6, h7, h8, h9 := feSquare(f)
 486
 487         h0 += h0
 488         h1 += h1
 489         h2 += h2
 490         h3 += h3
 491         h4 += h4
 492         h5 += h5
 493         h6 += h6
 494         h7 += h7
 495         h8 += h8
 496         h9 += h9
 497
 498         FeCombine(h, h0, h1, h2, h3, h4, h5, h6, h7, h8, h9)
 499 }
 500
 501 func FeInvert(out, z *FieldElement) {
 502         var t0, t1, t2, t3 FieldElement
 503         var i int
 504
 505         FeSquare(&t0, z)        // 2^1
 506         FeSquare(&t1, &t0)      // 2^2
 507         for i = 1; i < 2; i++ { // 2^3
 508                 FeSquare(&t1, &t1)
 509         }
 510         FeMul(&t1, z, &t1)      // 2^3 + 2^0
 511         FeMul(&t0, &t0, &t1)    // 2^3 + 2^1 + 2^0
 512         FeSquare(&t2, &t0)      // 2^4 + 2^2 + 2^1
 513         FeMul(&t1, &t1, &t2)    // 2^4 + 2^3 + 2^2 + 2^1 + 2^0
 514         FeSquare(&t2, &t1)      // 5,4,3,2,1
 515         for i = 1; i < 5; i++ { // 9,8,7,6,5
 516                 FeSquare(&t2, &t2)
 517         }
 518         FeMul(&t1, &t2, &t1)     // 9,8,7,6,5,4,3,2,1,0
 519         FeSquare(&t2, &t1)       // 10..1
 520         for i = 1; i < 10; i++ { // 19..10
 521                 FeSquare(&t2, &t2)
 522         }
 523         FeMul(&t2, &t2, &t1)     // 19..0
 524         FeSquare(&t3, &t2)       // 20..1
 525         for i = 1; i < 20; i++ { // 39..20
 526                 FeSquare(&t3, &t3)
 527         }
 528         FeMul(&t2, &t3, &t2)     // 39..0
 529         FeSquare(&t2, &t2)       // 40..1
 530         for i = 1; i < 10; i++ { // 49..10
 531                 FeSquare(&t2, &t2)
 532         }
 533         FeMul(&t1, &t2, &t1)     // 49..0
 534         FeSquare(&t2, &t1)       // 50..1
 535         for i = 1; i < 50; i++ { // 99..50
 536                 FeSquare(&t2, &t2)
 537         }
 538         FeMul(&t2, &t2, &t1)      // 99..0
 539         FeSquare(&t3, &t2)        // 100..1
 540         for i = 1; i < 100; i++ { // 199..100
 541                 FeSquare(&t3, &t3)
 542         }
 543         FeMul(&t2, &t3, &t2)     // 199..0
 544         FeSquare(&t2, &t2)       // 200..1
 545         for i = 1; i < 50; i++ { // 249..50
 546                 FeSquare(&t2, &t2)
 547         }
 548         FeMul(&t1, &t2, &t1)    // 249..0
 549         FeSquare(&t1, &t1)      // 250..1
 550         for i = 1; i < 5; i++ { // 254..5
 551                 FeSquare(&t1, &t1)
 552         }
 553         FeMul(out, &t1, &t0) // 254..5,3,1,0
 554 }
 555
 556 func fePow22523(out, z *FieldElement) {
 557         var t0, t1, t2 FieldElement
 558         var i int
 559
 560         FeSquare(&t0, z)
 561         for i = 1; i < 1; i++ {
 562                 FeSquare(&t0, &t0)
 563         }
 564         FeSquare(&t1, &t0)
 565         for i = 1; i < 2; i++ {
 566                 FeSquare(&t1, &t1)
 567         }
 568         FeMul(&t1, z, &t1)
 569         FeMul(&t0, &t0, &t1)
 570         FeSquare(&t0, &t0)
 571         for i = 1; i < 1; i++ {
 572                 FeSquare(&t0, &t0)
 573         }
 574         FeMul(&t0, &t1, &t0)
 575         FeSquare(&t1, &t0)
 576         for i = 1; i < 5; i++ {
 577                 FeSquare(&t1, &t1)
 578         }
 579         FeMul(&t0, &t1, &t0)
 580         FeSquare(&t1, &t0)
 581         for i = 1; i < 10; i++ {
 582                 FeSquare(&t1, &t1)
 583         }
 584         FeMul(&t1, &t1, &t0)
 585         FeSquare(&t2, &t1)
 586         for i = 1; i < 20; i++ {
 587                 FeSquare(&t2, &t2)
 588         }
 589         FeMul(&t1, &t2, &t1)
 590         FeSquare(&t1, &t1)
 591         for i = 1; i < 10; i++ {
 592                 FeSquare(&t1, &t1)
 593         }
 594         FeMul(&t0, &t1, &t0)
 595         FeSquare(&t1, &t0)
 596         for i = 1; i < 50; i++ {
 597                 FeSquare(&t1, &t1)
 598         }
 599         FeMul(&t1, &t1, &t0)
 600         FeSquare(&t2, &t1)
 601         for i = 1; i < 100; i++ {
 602                 FeSquare(&t2, &t2)
 603         }
 604         FeMul(&t1, &t2, &t1)
 605         FeSquare(&t1, &t1)
 606         for i = 1; i < 50; i++ {
 607                 FeSquare(&t1, &t1)
 608         }
 609         FeMul(&t0, &t1, &t0)
 610         FeSquare(&t0, &t0)
 611         for i = 1; i < 2; i++ {
 612                 FeSquare(&t0, &t0)
 613         }
 614         FeMul(out, &t0, z)
 615 }
 616
 617 // Group elements are members of the elliptic curve -x^2 + y^2 = 1 + d * x^2 *
 618 // y^2 where d = -121665/121666.
 619 //
 620 // Several representations are used:
 621 //   ProjectiveGroupElement: (X:Y:Z) satisfying x=X/Z, y=Y/Z
 622 //   ExtendedGroupElement: (X:Y:Z:T) satisfying x=X/Z, y=Y/Z, XY=ZT
 623 //   CompletedGroupElement: ((X:Z),(Y:T)) satisfying x=X/Z, y=Y/T
 624 //   PreComputedGroupElement: (y+x,y-x,2dxy)
 625
 626 type ProjectiveGroupElement struct {
 627         X, Y, Z FieldElement
 628 }
 629
 630 type ExtendedGroupElement struct {
 631         X, Y, Z, T FieldElement
 632 }
 633
 634 type CompletedGroupElement struct {
 635         X, Y, Z, T FieldElement
 636 }
 637
 638 type PreComputedGroupElement struct {
 639         yPlusX, yMinusX, xy2d FieldElement
 640 }
 641
 642 type CachedGroupElement struct {
 643         yPlusX, yMinusX, Z, T2d FieldElement
 644 }
 645
 646 func (p *ProjectiveGroupElement) Zero() {
 647         FeZero(&p.X)
 648         FeOne(&p.Y)
 649         FeOne(&p.Z)
 650 }
 651
 652 func (p *ProjectiveGroupElement) Double(r *CompletedGroupElement) {
 653         var t0 FieldElement
 654
 655         FeSquare(&r.X, &p.X)
 656         FeSquare(&r.Z, &p.Y)
 657         FeSquare2(&r.T, &p.Z)
 658         FeAdd(&r.Y, &p.X, &p.Y)
 659         FeSquare(&t0, &r.Y)
 660         FeAdd(&r.Y, &r.Z, &r.X)
 661         FeSub(&r.Z, &r.Z, &r.X)
 662         FeSub(&r.X, &t0, &r.Y)
 663         FeSub(&r.T, &r.T, &r.Z)
 664 }
 665
 666 func (p *ProjectiveGroupElement) ToBytes(s *[32]byte) {
 667         var recip, x, y FieldElement
 668
 669         FeInvert(&recip, &p.Z)
 670         FeMul(&x, &p.X, &recip)
 671         FeMul(&y, &p.Y, &recip)
 672         FeToBytes(s, &y)
 673         s[31] ^= FeIsNegative(&x) << 7
 674 }
 675
 676 func (p *ExtendedGroupElement) Zero() {
 677         FeZero(&p.X)
 678         FeOne(&p.Y)
 679         FeOne(&p.Z)
 680         FeZero(&p.T)
 681 }
 682
 683 func (p *ExtendedGroupElement) Double(r *CompletedGroupElement) {
 684         var q ProjectiveGroupElement
 685         p.ToProjective(&q)
 686         q.Double(r)
 687 }
 688
 689 func (p *ExtendedGroupElement) ToCached(r *CachedGroupElement) {
 690         FeAdd(&r.yPlusX, &p.Y, &p.X)
 691         FeSub(&r.yMinusX, &p.Y, &p.X)
 692         FeCopy(&r.Z, &p.Z)
 693         FeMul(&r.T2d, &p.T, &d2)
 694 }
 695
 696 func (p *ExtendedGroupElement) ToProjective(r *ProjectiveGroupElement) {
 697         FeCopy(&r.X, &p.X)
 698         FeCopy(&r.Y, &p.Y)
 699         FeCopy(&r.Z, &p.Z)
 700 }
 701
 702 func (p *ExtendedGroupElement) ToBytes(s *[32]byte) {
 703         var recip, x, y FieldElement
 704
 705         FeInvert(&recip, &p.Z)
 706         FeMul(&x, &p.X, &recip)
 707         FeMul(&y, &p.Y, &recip)
 708         FeToBytes(s, &y)
 709         s[31] ^= FeIsNegative(&x) << 7
 710 }
 711
 712 func (p *ExtendedGroupElement) FromBytes(s *[32]byte) bool {
 713         var u, v, v3, vxx, check FieldElement
 714
 715         FeFromBytes(&p.Y, s)
 716         FeOne(&p.Z)
 717         FeSquare(&u, &p.Y)
 718         FeMul(&v, &u, &d)
 719         FeSub(&u, &u, &p.Z) // y = y^2-1
 720         FeAdd(&v, &v, &p.Z) // v = dy^2+1
 721
 722         FeSquare(&v3, &v)
 723         FeMul(&v3, &v3, &v) // v3 = v^3
 724         FeSquare(&p.X, &v3)
 725         FeMul(&p.X, &p.X, &v)
 726         FeMul(&p.X, &p.X, &u) // x = uv^7
 727
 728         fePow22523(&p.X, &p.X) // x = (uv^7)^((q-5)/8)
 729         FeMul(&p.X, &p.X, &v3)
 730         FeMul(&p.X, &p.X, &u) // x = uv^3(uv^7)^((q-5)/8)
 731
 732         var tmpX, tmp2 [32]byte
 733
 734         FeSquare(&vxx, &p.X)
 735         FeMul(&vxx, &vxx, &v)
 736         FeSub(&check, &vxx, &u) // vx^2-u
 737         if FeIsNonZero(&check) == 1 {
 738                 FeAdd(&check, &vxx, &u) // vx^2+u
 739                 if FeIsNonZero(&check) == 1 {
 740                         return false
 741                 }
 742                 FeMul(&p.X, &p.X, &SqrtM1)
 743
 744                 FeToBytes(&tmpX, &p.X)
 745                 for i, v := range tmpX {
 746                         tmp2[31-i] = v
 747                 }
 748         }
 749
 750         if FeIsNegative(&p.X) != (s[31] >> 7) {
 751                 FeNeg(&p.X, &p.X)
 752         }
 753
 754         FeMul(&p.T, &p.X, &p.Y)
 755         return true
 756 }
 757
 758 func (p *CompletedGroupElement) ToProjective(r *ProjectiveGroupElement) {
 759         FeMul(&r.X, &p.X, &p.T)
 760         FeMul(&r.Y, &p.Y, &p.Z)
 761         FeMul(&r.Z, &p.Z, &p.T)
 762 }
 763
 764 func (p *CompletedGroupElement) ToExtended(r *ExtendedGroupElement) {
 765         FeMul(&r.X, &p.X, &p.T)
 766         FeMul(&r.Y, &p.Y, &p.Z)
 767         FeMul(&r.Z, &p.Z, &p.T)
 768         FeMul(&r.T, &p.X, &p.Y)
 769 }
 770
 771 func (p *PreComputedGroupElement) Zero() {
 772         FeOne(&p.yPlusX)
 773         FeOne(&p.yMinusX)
 774         FeZero(&p.xy2d)
 775 }
 776
 777 func geAdd(r *CompletedGroupElement, p *ExtendedGroupElement, q *CachedGroupElement) {
 778         var t0 FieldElement
 779
 780         FeAdd(&r.X, &p.Y, &p.X)
 781         FeSub(&r.Y, &p.Y, &p.X)
 782         FeMul(&r.Z, &r.X, &q.yPlusX)
 783         FeMul(&r.Y, &r.Y, &q.yMinusX)
 784         FeMul(&r.T, &q.T2d, &p.T)
 785         FeMul(&r.X, &p.Z, &q.Z)
 786         FeAdd(&t0, &r.X, &r.X)
 787         FeSub(&r.X, &r.Z, &r.Y)
 788         FeAdd(&r.Y, &r.Z, &r.Y)
 789         FeAdd(&r.Z, &t0, &r.T)
 790         FeSub(&r.T, &t0, &r.T)
 791 }
 792
 793 func geSub(r *CompletedGroupElement, p *ExtendedGroupElement, q *CachedGroupElement) {
 794         var t0 FieldElement
 795
 796         FeAdd(&r.X, &p.Y, &p.X)
 797         FeSub(&r.Y, &p.Y, &p.X)
 798         FeMul(&r.Z, &r.X, &q.yMinusX)
 799         FeMul(&r.Y, &r.Y, &q.yPlusX)
 800         FeMul(&r.T, &q.T2d, &p.T)
 801         FeMul(&r.X, &p.Z, &q.Z)
 802         FeAdd(&t0, &r.X, &r.X)
 803         FeSub(&r.X, &r.Z, &r.Y)
 804         FeAdd(&r.Y, &r.Z, &r.Y)
 805         FeSub(&r.Z, &t0, &r.T)
 806         FeAdd(&r.T, &t0, &r.T)
 807 }
 808
 809 func geMixedAdd(r *CompletedGroupElement, p *ExtendedGroupElement, q *PreComputedGroupElement) {
 810         var t0 FieldElement
 811
 812         FeAdd(&r.X, &p.Y, &p.X)
 813         FeSub(&r.Y, &p.Y, &p.X)
 814         FeMul(&r.Z, &r.X, &q.yPlusX)
 815         FeMul(&r.Y, &r.Y, &q.yMinusX)
 816         FeMul(&r.T, &q.xy2d, &p.T)
 817         FeAdd(&t0, &p.Z, &p.Z)
 818         FeSub(&r.X, &r.Z, &r.Y)
 819         FeAdd(&r.Y, &r.Z, &r.Y)
 820         FeAdd(&r.Z, &t0, &r.T)
 821         FeSub(&r.T, &t0, &r.T)
 822 }
 823
 824 func geMixedSub(r *CompletedGroupElement, p *ExtendedGroupElement, q *PreComputedGroupElement) {
 825         var t0 FieldElement
 826
 827         FeAdd(&r.X, &p.Y, &p.X)
 828         FeSub(&r.Y, &p.Y, &p.X)
 829         FeMul(&r.Z, &r.X, &q.yMinusX)
 830         FeMul(&r.Y, &r.Y, &q.yPlusX)
 831         FeMul(&r.T, &q.xy2d, &p.T)
 832         FeAdd(&t0, &p.Z, &p.Z)
 833         FeSub(&r.X, &r.Z, &r.Y)
 834         FeAdd(&r.Y, &r.Z, &r.Y)
 835         FeSub(&r.Z, &t0, &r.T)
 836         FeAdd(&r.T, &t0, &r.T)
 837 }
 838
 839 func slide(r *[256]int8, a *[32]byte) {
 840         for i := range r {
 841                 r[i] = int8(1 & (a[i>>3] >> uint(i&7)))
 842         }
 843
 844         for i := range r {
 845                 if r[i] != 0 {
 846                         for b := 1; b <= 6 && i+b < 256; b++ {
 847                                 if r[i+b] != 0 {
 848                                         if r[i]+(r[i+b]<<uint(b)) <= 15 {
 849                                                 r[i] += r[i+b] << uint(b)
 850                                                 r[i+b] = 0
 851                                         } else if r[i]-(r[i+b]<<uint(b)) >= -15 {
 852                                                 r[i] -= r[i+b] << uint(b)
 853                                                 for k := i + b; k < 256; k++ {
 854                                                         if r[k] == 0 {
 855                                                                 r[k] = 1
 856                                                                 break
 857                                                         }
 858                                                         r[k] = 0
 859                                                 }
 860                                         } else {
 861                                                 break
 862                                         }
 863                                 }
 864                         }
 865                 }
 866         }
 867 }
 868
 869 // GeDoubleScalarMultVartime sets r = a*A + b*B
 870 // where a = a[0]+256*a[1]+...+256^31 a[31].
 871 // and b = b[0]+256*b[1]+...+256^31 b[31].
 872 // B is the Ed25519 base point (x,4/5) with x positive.
 873 func GeDoubleScalarMultVartime(r *ProjectiveGroupElement, a *[32]byte, A *ExtendedGroupElement, b *[32]byte) {
 874         var aSlide, bSlide [256]int8
 875         var Ai [8]CachedGroupElement // A,3A,5A,7A,9A,11A,13A,15A
 876         var t CompletedGroupElement
 877         var u, A2 ExtendedGroupElement
 878         var i int
 879
 880         slide(&aSlide, a)
 881         slide(&bSlide, b)
 882
 883         A.ToCached(&Ai[0])
 884         A.Double(&t)
 885         t.ToExtended(&A2)
 886
 887         for i := 0; i < 7; i++ {
 888                 geAdd(&t, &A2, &Ai[i])
 889                 t.ToExtended(&u)
 890                 u.ToCached(&Ai[i+1])
 891         }
 892
 893         r.Zero()
 894
 895         for i = 255; i >= 0; i-- {
 896                 if aSlide[i] != 0 || bSlide[i] != 0 {
 897                         break
 898                 }
 899         }
 900
 901         for ; i >= 0; i-- {
 902                 r.Double(&t)
 903
 904                 if aSlide[i] > 0 {
 905                         t.ToExtended(&u)
 906                         geAdd(&t, &u, &Ai[aSlide[i]/2])
 907                 } else if aSlide[i] < 0 {
 908                         t.ToExtended(&u)
 909                         geSub(&t, &u, &Ai[(-aSlide[i])/2])
 910                 }
 911
 912                 if bSlide[i] > 0 {
 913                         t.ToExtended(&u)
 914                         geMixedAdd(&t, &u, &bi[bSlide[i]/2])
 915                 } else if bSlide[i] < 0 {
 916                         t.ToExtended(&u)
 917                         geMixedSub(&t, &u, &bi[(-bSlide[i])/2])
 918                 }
 919
 920                 t.ToProjective(r)
 921         }
 922 }
 923
 924 // equal returns 1 if b == c and 0 otherwise, assuming that b and c are
 925 // non-negative.
 926 func equal(b, c int32) int32 {
 927         x := uint32(b ^ c)
 928         x--
 929         return int32(x >> 31)
 930 }
 931
 932 // negative returns 1 if b < 0 and 0 otherwise.
 933 func negative(b int32) int32 {
 934         return (b >> 31) & 1
 935 }
 936
 937 func PreComputedGroupElementCMove(t, u *PreComputedGroupElement, b int32) {
 938         FeCMove(&t.yPlusX, &u.yPlusX, b)
 939         FeCMove(&t.yMinusX, &u.yMinusX, b)
 940         FeCMove(&t.xy2d, &u.xy2d, b)
 941 }
 942
 943 func selectPoint(t *PreComputedGroupElement, pos int32, b int32) {
 944         var minusT PreComputedGroupElement
 945         bNegative := negative(b)
 946         bAbs := b - (((-bNegative) & b) << 1)
 947
 948         t.Zero()
 949         for i := int32(0); i < 8; i++ {
 950                 PreComputedGroupElementCMove(t, &base[pos][i], equal(bAbs, i+1))
 951         }
 952         FeCopy(&minusT.yPlusX, &t.yMinusX)
 953         FeCopy(&minusT.yMinusX, &t.yPlusX)
 954         FeNeg(&minusT.xy2d, &t.xy2d)
 955         PreComputedGroupElementCMove(t, &minusT, bNegative)
 956 }
 957
 958 // GeScalarMultBase computes h = a*B, where
 959 //   a = a[0]+256*a[1]+...+256^31 a[31]
 960 //   B is the Ed25519 base point (x,4/5) with x positive.
 961 //
 962 // Preconditions:
 963 //   a[31] <= 127
 964 func GeScalarMultBase(h *ExtendedGroupElement, a *[32]byte) {
 965         var e [64]int8
 966
 967         for i, v := range a {
 968                 e[2*i] = int8(v & 15)
 969                 e[2*i+1] = int8((v >> 4) & 15)
 970         }
 971
 972         // each e[i] is between 0 and 15 and e[63] is between 0 and 7.
 973
 974         carry := int8(0)
 975         for i := 0; i < 63; i++ {
 976                 e[i] += carry
 977                 carry = (e[i] + 8) >> 4
 978                 e[i] -= carry << 4
 979         }
 980         e[63] += carry
 981         // each e[i] is between -8 and 8.
 982
 983         h.Zero()
 984         var t PreComputedGroupElement
 985         var r CompletedGroupElement
 986         for i := int32(1); i < 64; i += 2 {
 987                 selectPoint(&t, i/2, int32(e[i]))
 988                 geMixedAdd(&r, h, &t)
 989                 r.ToExtended(h)
 990         }
 991
 992         var s ProjectiveGroupElement
 993
 994         h.Double(&r)
 995         r.ToProjective(&s)
 996         s.Double(&r)
 997         r.ToProjective(&s)
 998         s.Double(&r)
 999         r.ToProjective(&s)
1000         s.Double(&r)
1001         r.ToExtended(h)
1002
1003         for i := int32(0); i < 64; i += 2 {
1004                 selectPoint(&t, i/2, int32(e[i]))
1005                 geMixedAdd(&r, h, &t)
1006                 r.ToExtended(h)
1007         }
1008 }
1009
1010 // The scalars are GF(2^252 + 27742317777372353535851937790883648493).
1011
1012 // Input:
1013 //   a[0]+256*a[1]+...+256^31*a[31] = a
1014 //   b[0]+256*b[1]+...+256^31*b[31] = b
1015 //   c[0]+256*c[1]+...+256^31*c[31] = c
1016 //
1017 // Output:
1018 //   s[0]+256*s[1]+...+256^31*s[31] = (ab+c) mod l
1019 //   where l = 2^252 + 27742317777372353535851937790883648493.
1020 func ScMulAdd(s, a, b, c *[32]byte) {
1021         a0 := 2097151 & load3(a[:])
1022         a1 := 2097151 & (load4(a[2:]) >> 5)
1023         a2 := 2097151 & (load3(a[5:]) >> 2)
1024         a3 := 2097151 & (load4(a[7:]) >> 7)
1025         a4 := 2097151 & (load4(a[10:]) >> 4)
1026         a5 := 2097151 & (load3(a[13:]) >> 1)
1027         a6 := 2097151 & (load4(a[15:]) >> 6)
1028         a7 := 2097151 & (load3(a[18:]) >> 3)
1029         a8 := 2097151 & load3(a[21:])
1030         a9 := 2097151 & (load4(a[23:]) >> 5)
1031         a10 := 2097151 & (load3(a[26:]) >> 2)
1032         a11 := (load4(a[28:]) >> 7)
1033         b0 := 2097151 & load3(b[:])
1034         b1 := 2097151 & (load4(b[2:]) >> 5)
1035         b2 := 2097151 & (load3(b[5:]) >> 2)
1036         b3 := 2097151 & (load4(b[7:]) >> 7)
1037         b4 := 2097151 & (load4(b[10:]) >> 4)
1038         b5 := 2097151 & (load3(b[13:]) >> 1)
1039         b6 := 2097151 & (load4(b[15:]) >> 6)
1040         b7 := 2097151 & (load3(b[18:]) >> 3)
1041         b8 := 2097151 & load3(b[21:])
1042         b9 := 2097151 & (load4(b[23:]) >> 5)
1043         b10 := 2097151 & (load3(b[26:]) >> 2)
1044         b11 := (load4(b[28:]) >> 7)
1045         c0 := 2097151 & load3(c[:])
1046         c1 := 2097151 & (load4(c[2:]) >> 5)
1047         c2 := 2097151 & (load3(c[5:]) >> 2)
1048         c3 := 2097151 & (load4(c[7:]) >> 7)
1049         c4 := 2097151 & (load4(c[10:]) >> 4)
1050         c5 := 2097151 & (load3(c[13:]) >> 1)
1051         c6 := 2097151 & (load4(c[15:]) >> 6)
1052         c7 := 2097151 & (load3(c[18:]) >> 3)
1053         c8 := 2097151 & load3(c[21:])
1054         c9 := 2097151 & (load4(c[23:]) >> 5)
1055         c10 := 2097151 & (load3(c[26:]) >> 2)
1056         c11 := (load4(c[28:]) >> 7)
1057         var carry [23]int64
1058
1059         s0 := c0 + a0*b0
1060         s1 := c1 + a0*b1 + a1*b0
1061         s2 := c2 + a0*b2 + a1*b1 + a2*b0
1062         s3 := c3 + a0*b3 + a1*b2 + a2*b1 + a3*b0
1063         s4 := c4 + a0*b4 + a1*b3 + a2*b2 + a3*b1 + a4*b0
1064         s5 := c5 + a0*b5 + a1*b4 + a2*b3 + a3*b2 + a4*b1 + a5*b0
1065         s6 := c6 + a0*b6 + a1*b5 + a2*b4 + a3*b3 + a4*b2 + a5*b1 + a6*b0
1066         s7 := c7 + a0*b7 + a1*b6 + a2*b5 + a3*b4 + a4*b3 + a5*b2 + a6*b1 + a7*b0
1067         s8 := c8 + a0*b8 + a1*b7 + a2*b6 + a3*b5 + a4*b4 + a5*b3 + a6*b2 + a7*b1 + a8*b0
1068         s9 := c9 + a0*b9 + a1*b8 + a2*b7 + a3*b6 + a4*b5 + a5*b4 + a6*b3 + a7*b2 + a8*b1 + a9*b0
1069         s10 := c10 + a0*b10 + a1*b9 + a2*b8 + a3*b7 + a4*b6 + a5*b5 + a6*b4 + a7*b3 + a8*b2 + a9*b1 + a10*b0
1070         s11 := c11 + a0*b11 + a1*b10 + a2*b9 + a3*b8 + a4*b7 + a5*b6 + a6*b5 + a7*b4 + a8*b3 + a9*b2 + a10*b1 + a11*b0
1071         s12 := a1*b11 + a2*b10 + a3*b9 + a4*b8 + a5*b7 + a6*b6 + a7*b5 + a8*b4 + a9*b3 + a10*b2 + a11*b1
1072         s13 := a2*b11 + a3*b10 + a4*b9 + a5*b8 + a6*b7 + a7*b6 + a8*b5 + a9*b4 + a10*b3 + a11*b2
1073         s14 := a3*b11 + a4*b10 + a5*b9 + a6*b8 + a7*b7 + a8*b6 + a9*b5 + a10*b4 + a11*b3
1074         s15 := a4*b11 + a5*b10 + a6*b9 + a7*b8 + a8*b7 + a9*b6 + a10*b5 + a11*b4
1075         s16 := a5*b11 + a6*b10 + a7*b9 + a8*b8 + a9*b7 + a10*b6 + a11*b5
1076         s17 := a6*b11 + a7*b10 + a8*b9 + a9*b8 + a10*b7 + a11*b6
1077         s18 := a7*b11 + a8*b10 + a9*b9 + a10*b8 + a11*b7
1078         s19 := a8*b11 + a9*b10 + a10*b9 + a11*b8
1079         s20 := a9*b11 + a10*b10 + a11*b9
1080         s21 := a10*b11 + a11*b10
1081         s22 := a11 * b11
1082         s23 := int64(0)
1083
1084         carry[0] = (s0 + (1 << 20)) >> 21
1085         s1 += carry[0]
1086         s0 -= carry[0] << 21
1087         carry[2] = (s2 + (1 << 20)) >> 21
1088         s3 += carry[2]
1089         s2 -= carry[2] << 21
1090         carry[4] = (s4 + (1 << 20)) >> 21
1091         s5 += carry[4]
1092         s4 -= carry[4] << 21
1093         carry[6] = (s6 + (1 << 20)) >> 21
1094         s7 += carry[6]
1095         s6 -= carry[6] << 21
1096         carry[8] = (s8 + (1 << 20)) >> 21
1097         s9 += carry[8]
1098         s8 -= carry[8] << 21
1099         carry[10] = (s10 + (1 << 20)) >> 21
1100         s11 += carry[10]
1101         s10 -= carry[10] << 21
1102         carry[12] = (s12 + (1 << 20)) >> 21
1103         s13 += carry[12]
1104         s12 -= carry[12] << 21
1105         carry[14] = (s14 + (1 << 20)) >> 21
1106         s15 += carry[14]
1107         s14 -= carry[14] << 21
1108         carry[16] = (s16 + (1 << 20)) >> 21
1109         s17 += carry[16]
1110         s16 -= carry[16] << 21
1111         carry[18] = (s18 + (1 << 20)) >> 21
1112         s19 += carry[18]
1113         s18 -= carry[18] << 21
1114         carry[20] = (s20 + (1 << 20)) >> 21
1115         s21 += carry[20]
1116         s20 -= carry[20] << 21
1117         carry[22] = (s22 + (1 << 20)) >> 21
1118         s23 += carry[22]
1119         s22 -= carry[22] << 21
1120
1121         carry[1] = (s1 + (1 << 20)) >> 21
1122         s2 += carry[1]
1123         s1 -= carry[1] << 21
1124         carry[3] = (s3 + (1 << 20)) >> 21
1125         s4 += carry[3]
1126         s3 -= carry[3] << 21
1127         carry[5] = (s5 + (1 << 20)) >> 21
1128         s6 += carry[5]
1129         s5 -= carry[5] << 21
1130         carry[7] = (s7 + (1 << 20)) >> 21
1131         s8 += carry[7]
1132         s7 -= carry[7] << 21
1133         carry[9] = (s9 + (1 << 20)) >> 21
1134         s10 += carry[9]
1135         s9 -= carry[9] << 21
1136         carry[11] = (s11 + (1 << 20)) >> 21
1137         s12 += carry[11]
1138         s11 -= carry[11] << 21
1139         carry[13] = (s13 + (1 << 20)) >> 21
1140         s14 += carry[13]
1141         s13 -= carry[13] << 21
1142         carry[15] = (s15 + (1 << 20)) >> 21
1143         s16 += carry[15]
1144         s15 -= carry[15] << 21
1145         carry[17] = (s17 + (1 << 20)) >> 21
1146         s18 += carry[17]
1147         s17 -= carry[17] << 21
1148         carry[19] = (s19 + (1 << 20)) >> 21
1149         s20 += carry[19]
1150         s19 -= carry[19] << 21
1151         carry[21] = (s21 + (1 << 20)) >> 21
1152         s22 += carry[21]
1153         s21 -= carry[21] << 21
1154
1155         s11 += s23 * 666643
1156         s12 += s23 * 470296
1157         s13 += s23 * 654183
1158         s14 -= s23 * 997805
1159         s15 += s23 * 136657
1160         s16 -= s23 * 683901
1161         s23 = 0
1162
1163         s10 += s22 * 666643
1164         s11 += s22 * 470296
1165         s12 += s22 * 654183
1166         s13 -= s22 * 997805
1167         s14 += s22 * 136657
1168         s15 -= s22 * 683901
1169         s22 = 0
1170
1171         s9 += s21 * 666643
1172         s10 += s21 * 470296
1173         s11 += s21 * 654183
1174         s12 -= s21 * 997805
1175         s13 += s21 * 136657
1176         s14 -= s21 * 683901
1177         s21 = 0
1178
1179         s8 += s20 * 666643
1180         s9 += s20 * 470296
1181         s10 += s20 * 654183
1182         s11 -= s20 * 997805
1183         s12 += s20 * 136657
1184         s13 -= s20 * 683901
1185         s20 = 0
1186
1187         s7 += s19 * 666643
1188         s8 += s19 * 470296
1189         s9 += s19 * 654183
1190         s10 -= s19 * 997805
1191         s11 += s19 * 136657
1192         s12 -= s19 * 683901
1193         s19 = 0
1194
1195         s6 += s18 * 666643
1196         s7 += s18 * 470296
1197         s8 += s18 * 654183
1198         s9 -= s18 * 997805
1199         s10 += s18 * 136657
1200         s11 -= s18 * 683901
1201         s18 = 0
1202
1203         carry[6] = (s6 + (1 << 20)) >> 21
1204         s7 += carry[6]
1205         s6 -= carry[6] << 21
1206         carry[8] = (s8 + (1 << 20)) >> 21
1207         s9 += carry[8]
1208         s8 -= carry[8] << 21
1209         carry[10] = (s10 + (1 << 20)) >> 21
1210         s11 += carry[10]
1211         s10 -= carry[10] << 21
1212         carry[12] = (s12 + (1 << 20)) >> 21
1213         s13 += carry[12]
1214         s12 -= carry[12] << 21
1215         carry[14] = (s14 + (1 << 20)) >> 21
1216         s15 += carry[14]
1217         s14 -= carry[14] << 21
1218         carry[16] = (s16 + (1 << 20)) >> 21
1219         s17 += carry[16]
1220         s16 -= carry[16] << 21
1221
1222         carry[7] = (s7 + (1 << 20)) >> 21
1223         s8 += carry[7]
1224         s7 -= carry[7] << 21
1225         carry[9] = (s9 + (1 << 20)) >> 21
1226         s10 += carry[9]
1227         s9 -= carry[9] << 21
1228         carry[11] = (s11 + (1 << 20)) >> 21
1229         s12 += carry[11]
1230         s11 -= carry[11] << 21
1231         carry[13] = (s13 + (1 << 20)) >> 21
1232         s14 += carry[13]
1233         s13 -= carry[13] << 21
1234         carry[15] = (s15 + (1 << 20)) >> 21
1235         s16 += carry[15]
1236         s15 -= carry[15] << 21
1237
1238         s5 += s17 * 666643
1239         s6 += s17 * 470296
1240         s7 += s17 * 654183
1241         s8 -= s17 * 997805
1242         s9 += s17 * 136657
1243         s10 -= s17 * 683901
1244         s17 = 0
1245
1246         s4 += s16 * 666643
1247         s5 += s16 * 470296
1248         s6 += s16 * 654183
1249         s7 -= s16 * 997805
1250         s8 += s16 * 136657
1251         s9 -= s16 * 683901
1252         s16 = 0
1253
1254         s3 += s15 * 666643
1255         s4 += s15 * 470296
1256         s5 += s15 * 654183
1257         s6 -= s15 * 997805
1258         s7 += s15 * 136657
1259         s8 -= s15 * 683901
1260         s15 = 0
1261
1262         s2 += s14 * 666643
1263         s3 += s14 * 470296
1264         s4 += s14 * 654183
1265         s5 -= s14 * 997805
1266         s6 += s14 * 136657
1267         s7 -= s14 * 683901
1268         s14 = 0
1269
1270         s1 += s13 * 666643
1271         s2 += s13 * 470296
1272         s3 += s13 * 654183
1273         s4 -= s13 * 997805
1274         s5 += s13 * 136657
1275         s6 -= s13 * 683901
1276         s13 = 0
1277
1278         s0 += s12 * 666643
1279         s1 += s12 * 470296
1280         s2 += s12 * 654183
1281         s3 -= s12 * 997805
1282         s4 += s12 * 136657
1283         s5 -= s12 * 683901
1284         s12 = 0
1285
1286         carry[0] = (s0 + (1 << 20)) >> 21
1287         s1 += carry[0]
1288         s0 -= carry[0] << 21
1289         carry[2] = (s2 + (1 << 20)) >> 21
1290         s3 += carry[2]
1291         s2 -= carry[2] << 21
1292         carry[4] = (s4 + (1 << 20)) >> 21
1293         s5 += carry[4]
1294         s4 -= carry[4] << 21
1295         carry[6] = (s6 + (1 << 20)) >> 21
1296         s7 += carry[6]
1297         s6 -= carry[6] << 21
1298         carry[8] = (s8 + (1 << 20)) >> 21
1299         s9 += carry[8]
1300         s8 -= carry[8] << 21
1301         carry[10] = (s10 + (1 << 20)) >> 21
1302         s11 += carry[10]
1303         s10 -= carry[10] << 21
1304
1305         carry[1] = (s1 + (1 << 20)) >> 21
1306         s2 += carry[1]
1307         s1 -= carry[1] << 21
1308         carry[3] = (s3 + (1 << 20)) >> 21
1309         s4 += carry[3]
1310         s3 -= carry[3] << 21
1311         carry[5] = (s5 + (1 << 20)) >> 21
1312         s6 += carry[5]
1313         s5 -= carry[5] << 21
1314         carry[7] = (s7 + (1 << 20)) >> 21
1315         s8 += carry[7]
1316         s7 -= carry[7] << 21
1317         carry[9] = (s9 + (1 << 20)) >> 21
1318         s10 += carry[9]
1319         s9 -= carry[9] << 21
1320         carry[11] = (s11 + (1 << 20)) >> 21
1321         s12 += carry[11]
1322         s11 -= carry[11] << 21
1323
1324         s0 += s12 * 666643
1325         s1 += s12 * 470296
1326         s2 += s12 * 654183
1327         s3 -= s12 * 997805
1328         s4 += s12 * 136657
1329         s5 -= s12 * 683901
1330         s12 = 0
1331
1332         carry[0] = s0 >> 21
1333         s1 += carry[0]
1334         s0 -= carry[0] << 21
1335         carry[1] = s1 >> 21
1336         s2 += carry[1]
1337         s1 -= carry[1] << 21
1338         carry[2] = s2 >> 21
1339         s3 += carry[2]
1340         s2 -= carry[2] << 21
1341         carry[3] = s3 >> 21
1342         s4 += carry[3]
1343         s3 -= carry[3] << 21
1344         carry[4] = s4 >> 21
1345         s5 += carry[4]
1346         s4 -= carry[4] << 21
1347         carry[5] = s5 >> 21
1348         s6 += carry[5]
1349         s5 -= carry[5] << 21
1350         carry[6] = s6 >> 21
1351         s7 += carry[6]
1352         s6 -= carry[6] << 21
1353         carry[7] = s7 >> 21
1354         s8 += carry[7]
1355         s7 -= carry[7] << 21
1356         carry[8] = s8 >> 21
1357         s9 += carry[8]
1358         s8 -= carry[8] << 21
1359         carry[9] = s9 >> 21
1360         s10 += carry[9]
1361         s9 -= carry[9] << 21
1362         carry[10] = s10 >> 21
1363         s11 += carry[10]
1364         s10 -= carry[10] << 21
1365         carry[11] = s11 >> 21
1366         s12 += carry[11]
1367         s11 -= carry[11] << 21
1368
1369         s0 += s12 * 666643
1370         s1 += s12 * 470296
1371         s2 += s12 * 654183
1372         s3 -= s12 * 997805
1373         s4 += s12 * 136657
1374         s5 -= s12 * 683901
1375         s12 = 0
1376
1377         carry[0] = s0 >> 21
1378         s1 += carry[0]
1379         s0 -= carry[0] << 21
1380         carry[1] = s1 >> 21
1381         s2 += carry[1]
1382         s1 -= carry[1] << 21
1383         carry[2] = s2 >> 21
1384         s3 += carry[2]
1385         s2 -= carry[2] << 21
1386         carry[3] = s3 >> 21
1387         s4 += carry[3]
1388         s3 -= carry[3] << 21
1389         carry[4] = s4 >> 21
1390         s5 += carry[4]
1391         s4 -= carry[4] << 21
1392         carry[5] = s5 >> 21
1393         s6 += carry[5]
1394         s5 -= carry[5] << 21
1395         carry[6] = s6 >> 21
1396         s7 += carry[6]
1397         s6 -= carry[6] << 21
1398         carry[7] = s7 >> 21
1399         s8 += carry[7]
1400         s7 -= carry[7] << 21
1401         carry[8] = s8 >> 21
1402         s9 += carry[8]
1403         s8 -= carry[8] << 21
1404         carry[9] = s9 >> 21
1405         s10 += carry[9]
1406         s9 -= carry[9] << 21
1407         carry[10] = s10 >> 21
1408         s11 += carry[10]
1409         s10 -= carry[10] << 21
1410
1411         s[0] = byte(s0 >> 0)
1412         s[1] = byte(s0 >> 8)
1413         s[2] = byte((s0 >> 16) | (s1 << 5))
1414         s[3] = byte(s1 >> 3)
1415         s[4] = byte(s1 >> 11)
1416         s[5] = byte((s1 >> 19) | (s2 << 2))
1417         s[6] = byte(s2 >> 6)
1418         s[7] = byte((s2 >> 14) | (s3 << 7))
1419         s[8] = byte(s3 >> 1)
1420         s[9] = byte(s3 >> 9)
1421         s[10] = byte((s3 >> 17) | (s4 << 4))
1422         s[11] = byte(s4 >> 4)
1423         s[12] = byte(s4 >> 12)
1424         s[13] = byte((s4 >> 20) | (s5 << 1))
1425         s[14] = byte(s5 >> 7)
1426         s[15] = byte((s5 >> 15) | (s6 << 6))
1427         s[16] = byte(s6 >> 2)
1428         s[17] = byte(s6 >> 10)
1429         s[18] = byte((s6 >> 18) | (s7 << 3))
1430         s[19] = byte(s7 >> 5)
1431         s[20] = byte(s7 >> 13)
1432         s[21] = byte(s8 >> 0)
1433         s[22] = byte(s8 >> 8)
1434         s[23] = byte((s8 >> 16) | (s9 << 5))
1435         s[24] = byte(s9 >> 3)
1436         s[25] = byte(s9 >> 11)
1437         s[26] = byte((s9 >> 19) | (s10 << 2))
1438         s[27] = byte(s10 >> 6)
1439         s[28] = byte((s10 >> 14) | (s11 << 7))
1440         s[29] = byte(s11 >> 1)
1441         s[30] = byte(s11 >> 9)
1442         s[31] = byte(s11 >> 17)
1443 }
1444
1445 // Input:
1446 //   s[0]+256*s[1]+...+256^63*s[63] = s
1447 //
1448 // Output:
1449 //   s[0]+256*s[1]+...+256^31*s[31] = s mod l
1450 //   where l = 2^252 + 27742317777372353535851937790883648493.
1451 func ScReduce(out *[32]byte, s *[64]byte) {
1452         s0 := 2097151 & load3(s[:])
1453         s1 := 2097151 & (load4(s[2:]) >> 5)
1454         s2 := 2097151 & (load3(s[5:]) >> 2)
1455         s3 := 2097151 & (load4(s[7:]) >> 7)
1456         s4 := 2097151 & (load4(s[10:]) >> 4)
1457         s5 := 2097151 & (load3(s[13:]) >> 1)
1458         s6 := 2097151 & (load4(s[15:]) >> 6)
1459         s7 := 2097151 & (load3(s[18:]) >> 3)
1460         s8 := 2097151 & load3(s[21:])
1461         s9 := 2097151 & (load4(s[23:]) >> 5)
1462         s10 := 2097151 & (load3(s[26:]) >> 2)
1463         s11 := 2097151 & (load4(s[28:]) >> 7)
1464         s12 := 2097151 & (load4(s[31:]) >> 4)
1465         s13 := 2097151 & (load3(s[34:]) >> 1)
1466         s14 := 2097151 & (load4(s[36:]) >> 6)
1467         s15 := 2097151 & (load3(s[39:]) >> 3)
1468         s16 := 2097151 & load3(s[42:])
1469         s17 := 2097151 & (load4(s[44:]) >> 5)
1470         s18 := 2097151 & (load3(s[47:]) >> 2)
1471         s19 := 2097151 & (load4(s[49:]) >> 7)
1472         s20 := 2097151 & (load4(s[52:]) >> 4)
1473         s21 := 2097151 & (load3(s[55:]) >> 1)
1474         s22 := 2097151 & (load4(s[57:]) >> 6)
1475         s23 := (load4(s[60:]) >> 3)
1476
1477         s11 += s23 * 666643
1478         s12 += s23 * 470296
1479         s13 += s23 * 654183
1480         s14 -= s23 * 997805
1481         s15 += s23 * 136657
1482         s16 -= s23 * 683901
1483         s23 = 0
1484
1485         s10 += s22 * 666643
1486         s11 += s22 * 470296
1487         s12 += s22 * 654183
1488         s13 -= s22 * 997805
1489         s14 += s22 * 136657
1490         s15 -= s22 * 683901
1491         s22 = 0
1492
1493         s9 += s21 * 666643
1494         s10 += s21 * 470296
1495         s11 += s21 * 654183
1496         s12 -= s21 * 997805
1497         s13 += s21 * 136657
1498         s14 -= s21 * 683901
1499         s21 = 0
1500
1501         s8 += s20 * 666643
1502         s9 += s20 * 470296
1503         s10 += s20 * 654183
1504         s11 -= s20 * 997805
1505         s12 += s20 * 136657
1506         s13 -= s20 * 683901
1507         s20 = 0
1508
1509         s7 += s19 * 666643
1510         s8 += s19 * 470296
1511         s9 += s19 * 654183
1512         s10 -= s19 * 997805
1513         s11 += s19 * 136657
1514         s12 -= s19 * 683901
1515         s19 = 0
1516
1517         s6 += s18 * 666643
1518         s7 += s18 * 470296
1519         s8 += s18 * 654183
1520         s9 -= s18 * 997805
1521         s10 += s18 * 136657
1522         s11 -= s18 * 683901
1523         s18 = 0
1524
1525         var carry [17]int64
1526
1527         carry[6] = (s6 + (1 << 20)) >> 21
1528         s7 += carry[6]
1529         s6 -= carry[6] << 21
1530         carry[8] = (s8 + (1 << 20)) >> 21
1531         s9 += carry[8]
1532         s8 -= carry[8] << 21
1533         carry[10] = (s10 + (1 << 20)) >> 21
1534         s11 += carry[10]
1535         s10 -= carry[10] << 21
1536         carry[12] = (s12 + (1 << 20)) >> 21
1537         s13 += carry[12]
1538         s12 -= carry[12] << 21
1539         carry[14] = (s14 + (1 << 20)) >> 21
1540         s15 += carry[14]
1541         s14 -= carry[14] << 21
1542         carry[16] = (s16 + (1 << 20)) >> 21
1543         s17 += carry[16]
1544         s16 -= carry[16] << 21
1545
1546         carry[7] = (s7 + (1 << 20)) >> 21
1547         s8 += carry[7]
1548         s7 -= carry[7] << 21
1549         carry[9] = (s9 + (1 << 20)) >> 21
1550         s10 += carry[9]
1551         s9 -= carry[9] << 21
1552         carry[11] = (s11 + (1 << 20)) >> 21
1553         s12 += carry[11]
1554         s11 -= carry[11] << 21
1555         carry[13] = (s13 + (1 << 20)) >> 21
1556         s14 += carry[13]
1557         s13 -= carry[13] << 21
1558         carry[15] = (s15 + (1 << 20)) >> 21
1559         s16 += carry[15]
1560         s15 -= carry[15] << 21
1561
1562         s5 += s17 * 666643
1563         s6 += s17 * 470296
1564         s7 += s17 * 654183
1565         s8 -= s17 * 997805
1566         s9 += s17 * 136657
1567         s10 -= s17 * 683901
1568         s17 = 0
1569
1570         s4 += s16 * 666643
1571         s5 += s16 * 470296
1572         s6 += s16 * 654183
1573         s7 -= s16 * 997805
1574         s8 += s16 * 136657
1575         s9 -= s16 * 683901
1576         s16 = 0
1577
1578         s3 += s15 * 666643
1579         s4 += s15 * 470296
1580         s5 += s15 * 654183
1581         s6 -= s15 * 997805
1582         s7 += s15 * 136657
1583         s8 -= s15 * 683901
1584         s15 = 0
1585
1586         s2 += s14 * 666643
1587         s3 += s14 * 470296
1588         s4 += s14 * 654183
1589         s5 -= s14 * 997805
1590         s6 += s14 * 136657
1591         s7 -= s14 * 683901
1592         s14 = 0
1593
1594         s1 += s13 * 666643
1595         s2 += s13 * 470296
1596         s3 += s13 * 654183
1597         s4 -= s13 * 997805
1598         s5 += s13 * 136657
1599         s6 -= s13 * 683901
1600         s13 = 0
1601
1602         s0 += s12 * 666643
1603         s1 += s12 * 470296
1604         s2 += s12 * 654183
1605         s3 -= s12 * 997805
1606         s4 += s12 * 136657
1607         s5 -= s12 * 683901
1608         s12 = 0
1609
1610         carry[0] = (s0 + (1 << 20)) >> 21
1611         s1 += carry[0]
1612         s0 -= carry[0] << 21
1613         carry[2] = (s2 + (1 << 20)) >> 21
1614         s3 += carry[2]
1615         s2 -= carry[2] << 21
1616         carry[4] = (s4 + (1 << 20)) >> 21
1617         s5 += carry[4]
1618         s4 -= carry[4] << 21
1619         carry[6] = (s6 + (1 << 20)) >> 21
1620         s7 += carry[6]
1621         s6 -= carry[6] << 21
1622         carry[8] = (s8 + (1 << 20)) >> 21
1623         s9 += carry[8]
1624         s8 -= carry[8] << 21
1625         carry[10] = (s10 + (1 << 20)) >> 21
1626         s11 += carry[10]
1627         s10 -= carry[10] << 21
1628
1629         carry[1] = (s1 + (1 << 20)) >> 21
1630         s2 += carry[1]
1631         s1 -= carry[1] << 21
1632         carry[3] = (s3 + (1 << 20)) >> 21
1633         s4 += carry[3]
1634         s3 -= carry[3] << 21
1635         carry[5] = (s5 + (1 << 20)) >> 21
1636         s6 += carry[5]
1637         s5 -= carry[5] << 21
1638         carry[7] = (s7 + (1 << 20)) >> 21
1639         s8 += carry[7]
1640         s7 -= carry[7] << 21
1641         carry[9] = (s9 + (1 << 20)) >> 21
1642         s10 += carry[9]
1643         s9 -= carry[9] << 21
1644         carry[11] = (s11 + (1 << 20)) >> 21
1645         s12 += carry[11]
1646         s11 -= carry[11] << 21
1647
1648         s0 += s12 * 666643
1649         s1 += s12 * 470296
1650         s2 += s12 * 654183
1651         s3 -= s12 * 997805
1652         s4 += s12 * 136657
1653         s5 -= s12 * 683901
1654         s12 = 0
1655
1656         carry[0] = s0 >> 21
1657         s1 += carry[0]
1658         s0 -= carry[0] << 21
1659         carry[1] = s1 >> 21
1660         s2 += carry[1]
1661         s1 -= carry[1] << 21
1662         carry[2] = s2 >> 21
1663         s3 += carry[2]
1664         s2 -= carry[2] << 21
1665         carry[3] = s3 >> 21
1666         s4 += carry[3]
1667         s3 -= carry[3] << 21
1668         carry[4] = s4 >> 21
1669         s5 += carry[4]
1670         s4 -= carry[4] << 21
1671         carry[5] = s5 >> 21
1672         s6 += carry[5]
1673         s5 -= carry[5] << 21
1674         carry[6] = s6 >> 21
1675         s7 += carry[6]
1676         s6 -= carry[6] << 21
1677         carry[7] = s7 >> 21
1678         s8 += carry[7]
1679         s7 -= carry[7] << 21
1680         carry[8] = s8 >> 21
1681         s9 += carry[8]
1682         s8 -= carry[8] << 21
1683         carry[9] = s9 >> 21
1684         s10 += carry[9]
1685         s9 -= carry[9] << 21
1686         carry[10] = s10 >> 21
1687         s11 += carry[10]
1688         s10 -= carry[10] << 21
1689         carry[11] = s11 >> 21
1690         s12 += carry[11]
1691         s11 -= carry[11] << 21
1692
1693         s0 += s12 * 666643
1694         s1 += s12 * 470296
1695         s2 += s12 * 654183
1696         s3 -= s12 * 997805
1697         s4 += s12 * 136657
1698         s5 -= s12 * 683901
1699         s12 = 0
1700
1701         carry[0] = s0 >> 21
1702         s1 += carry[0]
1703         s0 -= carry[0] << 21
1704         carry[1] = s1 >> 21
1705         s2 += carry[1]
1706         s1 -= carry[1] << 21
1707         carry[2] = s2 >> 21
1708         s3 += carry[2]
1709         s2 -= carry[2] << 21
1710         carry[3] = s3 >> 21
1711         s4 += carry[3]
1712         s3 -= carry[3] << 21
1713         carry[4] = s4 >> 21
1714         s5 += carry[4]
1715         s4 -= carry[4] << 21
1716         carry[5] = s5 >> 21
1717         s6 += carry[5]
1718         s5 -= carry[5] << 21
1719         carry[6] = s6 >> 21
1720         s7 += carry[6]
1721         s6 -= carry[6] << 21
1722         carry[7] = s7 >> 21
1723         s8 += carry[7]
1724         s7 -= carry[7] << 21
1725         carry[8] = s8 >> 21
1726         s9 += carry[8]
1727         s8 -= carry[8] << 21
1728         carry[9] = s9 >> 21
1729         s10 += carry[9]
1730         s9 -= carry[9] << 21
1731         carry[10] = s10 >> 21
1732         s11 += carry[10]
1733         s10 -= carry[10] << 21
1734
1735         out[0] = byte(s0 >> 0)
1736         out[1] = byte(s0 >> 8)
1737         out[2] = byte((s0 >> 16) | (s1 << 5))
1738         out[3] = byte(s1 >> 3)
1739         out[4] = byte(s1 >> 11)
1740         out[5] = byte((s1 >> 19) | (s2 << 2))
1741         out[6] = byte(s2 >> 6)
1742         out[7] = byte((s2 >> 14) | (s3 << 7))
1743         out[8] = byte(s3 >> 1)
1744         out[9] = byte(s3 >> 9)
1745         out[10] = byte((s3 >> 17) | (s4 << 4))
1746         out[11] = byte(s4 >> 4)
1747         out[12] = byte(s4 >> 12)
1748         out[13] = byte((s4 >> 20) | (s5 << 1))
1749         out[14] = byte(s5 >> 7)
1750         out[15] = byte((s5 >> 15) | (s6 << 6))
1751         out[16] = byte(s6 >> 2)
1752         out[17] = byte(s6 >> 10)
1753         out[18] = byte((s6 >> 18) | (s7 << 3))
1754         out[19] = byte(s7 >> 5)
1755         out[20] = byte(s7 >> 13)
1756         out[21] = byte(s8 >> 0)
1757         out[22] = byte(s8 >> 8)
1758         out[23] = byte((s8 >> 16) | (s9 << 5))
1759         out[24] = byte(s9 >> 3)
1760         out[25] = byte(s9 >> 11)
1761         out[26] = byte((s9 >> 19) | (s10 << 2))
1762         out[27] = byte(s10 >> 6)
1763         out[28] = byte((s10 >> 14) | (s11 << 7))
1764         out[29] = byte(s11 >> 1)
1765         out[30] = byte(s11 >> 9)
1766         out[31] = byte(s11 >> 17)
1767 }