tinc-vpn.org Git - tinc/blob - src/ed25519/fe.c

   1 #include "fixedint.h"
   2 #include "fe.h"
   3
   4
   5 /*
   6     helper functions
   7 */
   8 static uint64_t load_3(const unsigned char *in) {
   9         uint64_t result;
  10
  11         result = in[0];
  12         result |= shlu64(in[1], 8);
  13         result |= shlu64(in[2], 16);
  14
  15         return result;
  16 }
  17
  18 static uint64_t load_4(const unsigned char *in) {
  19         uint64_t result;
  20
  21         result = in[0];
  22         result |= shlu64(in[1], 8);
  23         result |= shlu64(in[2], 16);
  24         result |= shlu64(in[3], 24);
  25
  26         return result;
  27 }
  28
  29
  30
  31 /*
  32     h = 0
  33 */
  34
  35 void fe_0(fe h) {
  36         h[0] = 0;
  37         h[1] = 0;
  38         h[2] = 0;
  39         h[3] = 0;
  40         h[4] = 0;
  41         h[5] = 0;
  42         h[6] = 0;
  43         h[7] = 0;
  44         h[8] = 0;
  45         h[9] = 0;
  46 }
  47
  48
  49
  50 /*
  51     h = 1
  52 */
  53
  54 void fe_1(fe h) {
  55         h[0] = 1;
  56         h[1] = 0;
  57         h[2] = 0;
  58         h[3] = 0;
  59         h[4] = 0;
  60         h[5] = 0;
  61         h[6] = 0;
  62         h[7] = 0;
  63         h[8] = 0;
  64         h[9] = 0;
  65 }
  66
  67
  68
  69 /*
  70     h = f + g
  71     Can overlap h with f or g.
  72
  73     Preconditions:
  74        |f| bounded by 1.1*2^25,1.1*2^24,1.1*2^25,1.1*2^24,etc.
  75        |g| bounded by 1.1*2^25,1.1*2^24,1.1*2^25,1.1*2^24,etc.
  76
  77     Postconditions:
  78        |h| bounded by 1.1*2^26,1.1*2^25,1.1*2^26,1.1*2^25,etc.
  79 */
  80
  81 void fe_add(fe h, const fe f, const fe g) {
  82         int32_t f0 = f[0];
  83         int32_t f1 = f[1];
  84         int32_t f2 = f[2];
  85         int32_t f3 = f[3];
  86         int32_t f4 = f[4];
  87         int32_t f5 = f[5];
  88         int32_t f6 = f[6];
  89         int32_t f7 = f[7];
  90         int32_t f8 = f[8];
  91         int32_t f9 = f[9];
  92         int32_t g0 = g[0];
  93         int32_t g1 = g[1];
  94         int32_t g2 = g[2];
  95         int32_t g3 = g[3];
  96         int32_t g4 = g[4];
  97         int32_t g5 = g[5];
  98         int32_t g6 = g[6];
  99         int32_t g7 = g[7];
 100         int32_t g8 = g[8];
 101         int32_t g9 = g[9];
 102         int32_t h0 = f0 + g0;
 103         int32_t h1 = f1 + g1;
 104         int32_t h2 = f2 + g2;
 105         int32_t h3 = f3 + g3;
 106         int32_t h4 = f4 + g4;
 107         int32_t h5 = f5 + g5;
 108         int32_t h6 = f6 + g6;
 109         int32_t h7 = f7 + g7;
 110         int32_t h8 = f8 + g8;
 111         int32_t h9 = f9 + g9;
 112
 113         h[0] = h0;
 114         h[1] = h1;
 115         h[2] = h2;
 116         h[3] = h3;
 117         h[4] = h4;
 118         h[5] = h5;
 119         h[6] = h6;
 120         h[7] = h7;
 121         h[8] = h8;
 122         h[9] = h9;
 123 }
 124
 125
 126
 127 /*
 128     Replace (f,g) with (g,g) if b == 1;
 129     replace (f,g) with (f,g) if b == 0.
 130
 131     Preconditions: b in {0,1}.
 132 */
 133
 134 void fe_cmov(fe f, const fe g, unsigned int b) {
 135         int32_t f0 = f[0];
 136         int32_t f1 = f[1];
 137         int32_t f2 = f[2];
 138         int32_t f3 = f[3];
 139         int32_t f4 = f[4];
 140         int32_t f5 = f[5];
 141         int32_t f6 = f[6];
 142         int32_t f7 = f[7];
 143         int32_t f8 = f[8];
 144         int32_t f9 = f[9];
 145         int32_t g0 = g[0];
 146         int32_t g1 = g[1];
 147         int32_t g2 = g[2];
 148         int32_t g3 = g[3];
 149         int32_t g4 = g[4];
 150         int32_t g5 = g[5];
 151         int32_t g6 = g[6];
 152         int32_t g7 = g[7];
 153         int32_t g8 = g[8];
 154         int32_t g9 = g[9];
 155         int32_t x0 = f0 ^ g0;
 156         int32_t x1 = f1 ^ g1;
 157         int32_t x2 = f2 ^ g2;
 158         int32_t x3 = f3 ^ g3;
 159         int32_t x4 = f4 ^ g4;
 160         int32_t x5 = f5 ^ g5;
 161         int32_t x6 = f6 ^ g6;
 162         int32_t x7 = f7 ^ g7;
 163         int32_t x8 = f8 ^ g8;
 164         int32_t x9 = f9 ^ g9;
 165
 166         b = (unsigned int)(- (int) b);  /* silence warning */
 167         x0 &= b;
 168         x1 &= b;
 169         x2 &= b;
 170         x3 &= b;
 171         x4 &= b;
 172         x5 &= b;
 173         x6 &= b;
 174         x7 &= b;
 175         x8 &= b;
 176         x9 &= b;
 177
 178         f[0] = f0 ^ x0;
 179         f[1] = f1 ^ x1;
 180         f[2] = f2 ^ x2;
 181         f[3] = f3 ^ x3;
 182         f[4] = f4 ^ x4;
 183         f[5] = f5 ^ x5;
 184         f[6] = f6 ^ x6;
 185         f[7] = f7 ^ x7;
 186         f[8] = f8 ^ x8;
 187         f[9] = f9 ^ x9;
 188 }
 189
 190 /*
 191     Replace (f,g) with (g,f) if b == 1;
 192     replace (f,g) with (f,g) if b == 0.
 193
 194     Preconditions: b in {0,1}.
 195 */
 196
 197 void fe_cswap(fe f, fe g, unsigned int b) {
 198         int32_t f0 = f[0];
 199         int32_t f1 = f[1];
 200         int32_t f2 = f[2];
 201         int32_t f3 = f[3];
 202         int32_t f4 = f[4];
 203         int32_t f5 = f[5];
 204         int32_t f6 = f[6];
 205         int32_t f7 = f[7];
 206         int32_t f8 = f[8];
 207         int32_t f9 = f[9];
 208         int32_t g0 = g[0];
 209         int32_t g1 = g[1];
 210         int32_t g2 = g[2];
 211         int32_t g3 = g[3];
 212         int32_t g4 = g[4];
 213         int32_t g5 = g[5];
 214         int32_t g6 = g[6];
 215         int32_t g7 = g[7];
 216         int32_t g8 = g[8];
 217         int32_t g9 = g[9];
 218         int32_t x0 = f0 ^ g0;
 219         int32_t x1 = f1 ^ g1;
 220         int32_t x2 = f2 ^ g2;
 221         int32_t x3 = f3 ^ g3;
 222         int32_t x4 = f4 ^ g4;
 223         int32_t x5 = f5 ^ g5;
 224         int32_t x6 = f6 ^ g6;
 225         int32_t x7 = f7 ^ g7;
 226         int32_t x8 = f8 ^ g8;
 227         int32_t x9 = f9 ^ g9;
 228         b = -b;
 229         x0 &= b;
 230         x1 &= b;
 231         x2 &= b;
 232         x3 &= b;
 233         x4 &= b;
 234         x5 &= b;
 235         x6 &= b;
 236         x7 &= b;
 237         x8 &= b;
 238         x9 &= b;
 239         f[0] = f0 ^ x0;
 240         f[1] = f1 ^ x1;
 241         f[2] = f2 ^ x2;
 242         f[3] = f3 ^ x3;
 243         f[4] = f4 ^ x4;
 244         f[5] = f5 ^ x5;
 245         f[6] = f6 ^ x6;
 246         f[7] = f7 ^ x7;
 247         f[8] = f8 ^ x8;
 248         f[9] = f9 ^ x9;
 249         g[0] = g0 ^ x0;
 250         g[1] = g1 ^ x1;
 251         g[2] = g2 ^ x2;
 252         g[3] = g3 ^ x3;
 253         g[4] = g4 ^ x4;
 254         g[5] = g5 ^ x5;
 255         g[6] = g6 ^ x6;
 256         g[7] = g7 ^ x7;
 257         g[8] = g8 ^ x8;
 258         g[9] = g9 ^ x9;
 259 }
 260
 261
 262
 263 /*
 264     h = f
 265 */
 266
 267 void fe_copy(fe h, const fe f) {
 268         int32_t f0 = f[0];
 269         int32_t f1 = f[1];
 270         int32_t f2 = f[2];
 271         int32_t f3 = f[3];
 272         int32_t f4 = f[4];
 273         int32_t f5 = f[5];
 274         int32_t f6 = f[6];
 275         int32_t f7 = f[7];
 276         int32_t f8 = f[8];
 277         int32_t f9 = f[9];
 278
 279         h[0] = f0;
 280         h[1] = f1;
 281         h[2] = f2;
 282         h[3] = f3;
 283         h[4] = f4;
 284         h[5] = f5;
 285         h[6] = f6;
 286         h[7] = f7;
 287         h[8] = f8;
 288         h[9] = f9;
 289 }
 290
 291
 292
 293 /*
 294     Ignores top bit of h.
 295 */
 296
 297 void fe_frombytes(fe h, const unsigned char *s) {
 298         int64_t h0 = load_4(s);
 299         int64_t h1 = load_3(s + 4) << 6;
 300         int64_t h2 = load_3(s + 7) << 5;
 301         int64_t h3 = load_3(s + 10) << 3;
 302         int64_t h4 = load_3(s + 13) << 2;
 303         int64_t h5 = load_4(s + 16);
 304         int64_t h6 = load_3(s + 20) << 7;
 305         int64_t h7 = load_3(s + 23) << 5;
 306         int64_t h8 = load_3(s + 26) << 4;
 307         int64_t h9 = (load_3(s + 29) & 8388607) << 2;
 308         int64_t carry0;
 309         int64_t carry1;
 310         int64_t carry2;
 311         int64_t carry3;
 312         int64_t carry4;
 313         int64_t carry5;
 314         int64_t carry6;
 315         int64_t carry7;
 316         int64_t carry8;
 317         int64_t carry9;
 318
 319         carry9 = (h9 + (1L << 24)) >> 25;
 320         h0 += carry9 * 19;
 321         h9 -= shl64(carry9, 25);
 322         carry1 = (h1 + (1L << 24)) >> 25;
 323         h2 += carry1;
 324         h1 -= shl64(carry1, 25);
 325         carry3 = (h3 + (1L << 24)) >> 25;
 326         h4 += carry3;
 327         h3 -= shl64(carry3, 25);
 328         carry5 = (h5 + (1L << 24)) >> 25;
 329         h6 += carry5;
 330         h5 -= shl64(carry5, 25);
 331         carry7 = (h7 + (1L << 24)) >> 25;
 332         h8 += carry7;
 333         h7 -= shl64(carry7, 25);
 334         carry0 = (h0 + (1L << 25)) >> 26;
 335         h1 += carry0;
 336         h0 -= shl64(carry0, 26);
 337         carry2 = (h2 + (1L << 25)) >> 26;
 338         h3 += carry2;
 339         h2 -= shl64(carry2, 26);
 340         carry4 = (h4 + (1L << 25)) >> 26;
 341         h5 += carry4;
 342         h4 -= shl64(carry4, 26);
 343         carry6 = (h6 + (1L << 25)) >> 26;
 344         h7 += carry6;
 345         h6 -= shl64(carry6, 26);
 346         carry8 = (h8 + (1L << 25)) >> 26;
 347         h9 += carry8;
 348         h8 -= shl64(carry8, 26);
 349
 350         h[0] = h0;
 351         h[1] = h1;
 352         h[2] = h2;
 353         h[3] = h3;
 354         h[4] = h4;
 355         h[5] = h5;
 356         h[6] = h6;
 357         h[7] = h7;
 358         h[8] = h8;
 359         h[9] = h9;
 360 }
 361
 362
 363
 364 void fe_invert(fe out, const fe z) {
 365         fe t0;
 366         fe t1;
 367         fe t2;
 368         fe t3;
 369         int i;
 370
 371         fe_sq(t0, z);
 372
 373         for(i = 1; i < 1; ++i) {
 374                 fe_sq(t0, t0);
 375         }
 376
 377         fe_sq(t1, t0);
 378
 379         for(i = 1; i < 2; ++i) {
 380                 fe_sq(t1, t1);
 381         }
 382
 383         fe_mul(t1, z, t1);
 384         fe_mul(t0, t0, t1);
 385         fe_sq(t2, t0);
 386
 387         for(i = 1; i < 1; ++i) {
 388                 fe_sq(t2, t2);
 389         }
 390
 391         fe_mul(t1, t1, t2);
 392         fe_sq(t2, t1);
 393
 394         for(i = 1; i < 5; ++i) {
 395                 fe_sq(t2, t2);
 396         }
 397
 398         fe_mul(t1, t2, t1);
 399         fe_sq(t2, t1);
 400
 401         for(i = 1; i < 10; ++i) {
 402                 fe_sq(t2, t2);
 403         }
 404
 405         fe_mul(t2, t2, t1);
 406         fe_sq(t3, t2);
 407
 408         for(i = 1; i < 20; ++i) {
 409                 fe_sq(t3, t3);
 410         }
 411
 412         fe_mul(t2, t3, t2);
 413         fe_sq(t2, t2);
 414
 415         for(i = 1; i < 10; ++i) {
 416                 fe_sq(t2, t2);
 417         }
 418
 419         fe_mul(t1, t2, t1);
 420         fe_sq(t2, t1);
 421
 422         for(i = 1; i < 50; ++i) {
 423                 fe_sq(t2, t2);
 424         }
 425
 426         fe_mul(t2, t2, t1);
 427         fe_sq(t3, t2);
 428
 429         for(i = 1; i < 100; ++i) {
 430                 fe_sq(t3, t3);
 431         }
 432
 433         fe_mul(t2, t3, t2);
 434         fe_sq(t2, t2);
 435
 436         for(i = 1; i < 50; ++i) {
 437                 fe_sq(t2, t2);
 438         }
 439
 440         fe_mul(t1, t2, t1);
 441         fe_sq(t1, t1);
 442
 443         for(i = 1; i < 5; ++i) {
 444                 fe_sq(t1, t1);
 445         }
 446
 447         fe_mul(out, t1, t0);
 448 }
 449
 450
 451
 452 /*
 453     return 1 if f is in {1,3,5,...,q-2}
 454     return 0 if f is in {0,2,4,...,q-1}
 455
 456     Preconditions:
 457        |f| bounded by 1.1*2^26,1.1*2^25,1.1*2^26,1.1*2^25,etc.
 458 */
 459
 460 int fe_isnegative(const fe f) {
 461         unsigned char s[32];
 462
 463         fe_tobytes(s, f);
 464
 465         return s[0] & 1;
 466 }
 467
 468
 469
 470 /*
 471     return 1 if f == 0
 472     return 0 if f != 0
 473
 474     Preconditions:
 475        |f| bounded by 1.1*2^26,1.1*2^25,1.1*2^26,1.1*2^25,etc.
 476 */
 477
 478 int fe_isnonzero(const fe f) {
 479         unsigned char s[32];
 480         unsigned char r;
 481
 482         fe_tobytes(s, f);
 483
 484         r = s[0];
 485 #define F(i) r |= s[i]
 486         F(1);
 487         F(2);
 488         F(3);
 489         F(4);
 490         F(5);
 491         F(6);
 492         F(7);
 493         F(8);
 494         F(9);
 495         F(10);
 496         F(11);
 497         F(12);
 498         F(13);
 499         F(14);
 500         F(15);
 501         F(16);
 502         F(17);
 503         F(18);
 504         F(19);
 505         F(20);
 506         F(21);
 507         F(22);
 508         F(23);
 509         F(24);
 510         F(25);
 511         F(26);
 512         F(27);
 513         F(28);
 514         F(29);
 515         F(30);
 516         F(31);
 517 #undef F
 518
 519         return r != 0;
 520 }
 521
 522
 523
 524 /*
 525     h = f * g
 526     Can overlap h with f or g.
 527
 528     Preconditions:
 529        |f| bounded by 1.65*2^26,1.65*2^25,1.65*2^26,1.65*2^25,etc.
 530        |g| bounded by 1.65*2^26,1.65*2^25,1.65*2^26,1.65*2^25,etc.
 531
 532     Postconditions:
 533        |h| bounded by 1.01*2^25,1.01*2^24,1.01*2^25,1.01*2^24,etc.
 534     */
 535
 536 /*
 537 Notes on implementation strategy:
 538
 539 Using schoolbook multiplication.
 540 Karatsuba would save a little in some cost models.
 541
 542 Most multiplications by 2 and 19 are 32-bit precomputations;
 543 cheaper than 64-bit postcomputations.
 544
 545 There is one remaining multiplication by 19 in the carry chain;
 546 one *19 precomputation can be merged into this,
 547 but the resulting data flow is considerably less clean.
 548
 549 There are 12 carries below.
 550 10 of them are 2-way parallelizable and vectorizable.
 551 Can get away with 11 carries, but then data flow is much deeper.
 552
 553 With tighter constraints on inputs can squeeze carries into int32.
 554 */
 555
 556 void fe_mul(fe h, const fe f, const fe g) {
 557         int32_t f0 = f[0];
 558         int32_t f1 = f[1];
 559         int32_t f2 = f[2];
 560         int32_t f3 = f[3];
 561         int32_t f4 = f[4];
 562         int32_t f5 = f[5];
 563         int32_t f6 = f[6];
 564         int32_t f7 = f[7];
 565         int32_t f8 = f[8];
 566         int32_t f9 = f[9];
 567         int32_t g0 = g[0];
 568         int32_t g1 = g[1];
 569         int32_t g2 = g[2];
 570         int32_t g3 = g[3];
 571         int32_t g4 = g[4];
 572         int32_t g5 = g[5];
 573         int32_t g6 = g[6];
 574         int32_t g7 = g[7];
 575         int32_t g8 = g[8];
 576         int32_t g9 = g[9];
 577         int32_t g1_19 = 19 * g1; /* 1.959375*2^29 */
 578         int32_t g2_19 = 19 * g2; /* 1.959375*2^30; still ok */
 579         int32_t g3_19 = 19 * g3;
 580         int32_t g4_19 = 19 * g4;
 581         int32_t g5_19 = 19 * g5;
 582         int32_t g6_19 = 19 * g6;
 583         int32_t g7_19 = 19 * g7;
 584         int32_t g8_19 = 19 * g8;
 585         int32_t g9_19 = 19 * g9;
 586         int32_t f1_2 = 2 * f1;
 587         int32_t f3_2 = 2 * f3;
 588         int32_t f5_2 = 2 * f5;
 589         int32_t f7_2 = 2 * f7;
 590         int32_t f9_2 = 2 * f9;
 591         int64_t f0g0    = f0   * (int64_t) g0;
 592         int64_t f0g1    = f0   * (int64_t) g1;
 593         int64_t f0g2    = f0   * (int64_t) g2;
 594         int64_t f0g3    = f0   * (int64_t) g3;
 595         int64_t f0g4    = f0   * (int64_t) g4;
 596         int64_t f0g5    = f0   * (int64_t) g5;
 597         int64_t f0g6    = f0   * (int64_t) g6;
 598         int64_t f0g7    = f0   * (int64_t) g7;
 599         int64_t f0g8    = f0   * (int64_t) g8;
 600         int64_t f0g9    = f0   * (int64_t) g9;
 601         int64_t f1g0    = f1   * (int64_t) g0;
 602         int64_t f1g1_2  = f1_2 * (int64_t) g1;
 603         int64_t f1g2    = f1   * (int64_t) g2;
 604         int64_t f1g3_2  = f1_2 * (int64_t) g3;
 605         int64_t f1g4    = f1   * (int64_t) g4;
 606         int64_t f1g5_2  = f1_2 * (int64_t) g5;
 607         int64_t f1g6    = f1   * (int64_t) g6;
 608         int64_t f1g7_2  = f1_2 * (int64_t) g7;
 609         int64_t f1g8    = f1   * (int64_t) g8;
 610         int64_t f1g9_38 = f1_2 * (int64_t) g9_19;
 611         int64_t f2g0    = f2   * (int64_t) g0;
 612         int64_t f2g1    = f2   * (int64_t) g1;
 613         int64_t f2g2    = f2   * (int64_t) g2;
 614         int64_t f2g3    = f2   * (int64_t) g3;
 615         int64_t f2g4    = f2   * (int64_t) g4;
 616         int64_t f2g5    = f2   * (int64_t) g5;
 617         int64_t f2g6    = f2   * (int64_t) g6;
 618         int64_t f2g7    = f2   * (int64_t) g7;
 619         int64_t f2g8_19 = f2   * (int64_t) g8_19;
 620         int64_t f2g9_19 = f2   * (int64_t) g9_19;
 621         int64_t f3g0    = f3   * (int64_t) g0;
 622         int64_t f3g1_2  = f3_2 * (int64_t) g1;
 623         int64_t f3g2    = f3   * (int64_t) g2;
 624         int64_t f3g3_2  = f3_2 * (int64_t) g3;
 625         int64_t f3g4    = f3   * (int64_t) g4;
 626         int64_t f3g5_2  = f3_2 * (int64_t) g5;
 627         int64_t f3g6    = f3   * (int64_t) g6;
 628         int64_t f3g7_38 = f3_2 * (int64_t) g7_19;
 629         int64_t f3g8_19 = f3   * (int64_t) g8_19;
 630         int64_t f3g9_38 = f3_2 * (int64_t) g9_19;
 631         int64_t f4g0    = f4   * (int64_t) g0;
 632         int64_t f4g1    = f4   * (int64_t) g1;
 633         int64_t f4g2    = f4   * (int64_t) g2;
 634         int64_t f4g3    = f4   * (int64_t) g3;
 635         int64_t f4g4    = f4   * (int64_t) g4;
 636         int64_t f4g5    = f4   * (int64_t) g5;
 637         int64_t f4g6_19 = f4   * (int64_t) g6_19;
 638         int64_t f4g7_19 = f4   * (int64_t) g7_19;
 639         int64_t f4g8_19 = f4   * (int64_t) g8_19;
 640         int64_t f4g9_19 = f4   * (int64_t) g9_19;
 641         int64_t f5g0    = f5   * (int64_t) g0;
 642         int64_t f5g1_2  = f5_2 * (int64_t) g1;
 643         int64_t f5g2    = f5   * (int64_t) g2;
 644         int64_t f5g3_2  = f5_2 * (int64_t) g3;
 645         int64_t f5g4    = f5   * (int64_t) g4;
 646         int64_t f5g5_38 = f5_2 * (int64_t) g5_19;
 647         int64_t f5g6_19 = f5   * (int64_t) g6_19;
 648         int64_t f5g7_38 = f5_2 * (int64_t) g7_19;
 649         int64_t f5g8_19 = f5   * (int64_t) g8_19;
 650         int64_t f5g9_38 = f5_2 * (int64_t) g9_19;
 651         int64_t f6g0    = f6   * (int64_t) g0;
 652         int64_t f6g1    = f6   * (int64_t) g1;
 653         int64_t f6g2    = f6   * (int64_t) g2;
 654         int64_t f6g3    = f6   * (int64_t) g3;
 655         int64_t f6g4_19 = f6   * (int64_t) g4_19;
 656         int64_t f6g5_19 = f6   * (int64_t) g5_19;
 657         int64_t f6g6_19 = f6   * (int64_t) g6_19;
 658         int64_t f6g7_19 = f6   * (int64_t) g7_19;
 659         int64_t f6g8_19 = f6   * (int64_t) g8_19;
 660         int64_t f6g9_19 = f6   * (int64_t) g9_19;
 661         int64_t f7g0    = f7   * (int64_t) g0;
 662         int64_t f7g1_2  = f7_2 * (int64_t) g1;
 663         int64_t f7g2    = f7   * (int64_t) g2;
 664         int64_t f7g3_38 = f7_2 * (int64_t) g3_19;
 665         int64_t f7g4_19 = f7   * (int64_t) g4_19;
 666         int64_t f7g5_38 = f7_2 * (int64_t) g5_19;
 667         int64_t f7g6_19 = f7   * (int64_t) g6_19;
 668         int64_t f7g7_38 = f7_2 * (int64_t) g7_19;
 669         int64_t f7g8_19 = f7   * (int64_t) g8_19;
 670         int64_t f7g9_38 = f7_2 * (int64_t) g9_19;
 671         int64_t f8g0    = f8   * (int64_t) g0;
 672         int64_t f8g1    = f8   * (int64_t) g1;
 673         int64_t f8g2_19 = f8   * (int64_t) g2_19;
 674         int64_t f8g3_19 = f8   * (int64_t) g3_19;
 675         int64_t f8g4_19 = f8   * (int64_t) g4_19;
 676         int64_t f8g5_19 = f8   * (int64_t) g5_19;
 677         int64_t f8g6_19 = f8   * (int64_t) g6_19;
 678         int64_t f8g7_19 = f8   * (int64_t) g7_19;
 679         int64_t f8g8_19 = f8   * (int64_t) g8_19;
 680         int64_t f8g9_19 = f8   * (int64_t) g9_19;
 681         int64_t f9g0    = f9   * (int64_t) g0;
 682         int64_t f9g1_38 = f9_2 * (int64_t) g1_19;
 683         int64_t f9g2_19 = f9   * (int64_t) g2_19;
 684         int64_t f9g3_38 = f9_2 * (int64_t) g3_19;
 685         int64_t f9g4_19 = f9   * (int64_t) g4_19;
 686         int64_t f9g5_38 = f9_2 * (int64_t) g5_19;
 687         int64_t f9g6_19 = f9   * (int64_t) g6_19;
 688         int64_t f9g7_38 = f9_2 * (int64_t) g7_19;
 689         int64_t f9g8_19 = f9   * (int64_t) g8_19;
 690         int64_t f9g9_38 = f9_2 * (int64_t) g9_19;
 691         int64_t h0 = f0g0 + f1g9_38 + f2g8_19 + f3g7_38 + f4g6_19 + f5g5_38 + f6g4_19 + f7g3_38 + f8g2_19 + f9g1_38;
 692         int64_t h1 = f0g1 + f1g0   + f2g9_19 + f3g8_19 + f4g7_19 + f5g6_19 + f6g5_19 + f7g4_19 + f8g3_19 + f9g2_19;
 693         int64_t h2 = f0g2 + f1g1_2 + f2g0   + f3g9_38 + f4g8_19 + f5g7_38 + f6g6_19 + f7g5_38 + f8g4_19 + f9g3_38;
 694         int64_t h3 = f0g3 + f1g2   + f2g1   + f3g0   + f4g9_19 + f5g8_19 + f6g7_19 + f7g6_19 + f8g5_19 + f9g4_19;
 695         int64_t h4 = f0g4 + f1g3_2 + f2g2   + f3g1_2 + f4g0   + f5g9_38 + f6g8_19 + f7g7_38 + f8g6_19 + f9g5_38;
 696         int64_t h5 = f0g5 + f1g4   + f2g3   + f3g2   + f4g1   + f5g0   + f6g9_19 + f7g8_19 + f8g7_19 + f9g6_19;
 697         int64_t h6 = f0g6 + f1g5_2 + f2g4   + f3g3_2 + f4g2   + f5g1_2 + f6g0   + f7g9_38 + f8g8_19 + f9g7_38;
 698         int64_t h7 = f0g7 + f1g6   + f2g5   + f3g4   + f4g3   + f5g2   + f6g1   + f7g0   + f8g9_19 + f9g8_19;
 699         int64_t h8 = f0g8 + f1g7_2 + f2g6   + f3g5_2 + f4g4   + f5g3_2 + f6g2   + f7g1_2 + f8g0   + f9g9_38;
 700         int64_t h9 = f0g9 + f1g8   + f2g7   + f3g6   + f4g5   + f5g4   + f6g3   + f7g2   + f8g1   + f9g0   ;
 701         int64_t carry0;
 702         int64_t carry1;
 703         int64_t carry2;
 704         int64_t carry3;
 705         int64_t carry4;
 706         int64_t carry5;
 707         int64_t carry6;
 708         int64_t carry7;
 709         int64_t carry8;
 710         int64_t carry9;
 711
 712         carry0 = (h0 + (1L << 25)) >> 26;
 713         h1 += carry0;
 714         h0 -= shl64(carry0, 26);
 715         carry4 = (h4 + (1L << 25)) >> 26;
 716         h5 += carry4;
 717         h4 -= shl64(carry4, 26);
 718
 719         carry1 = (h1 + (1L << 24)) >> 25;
 720         h2 += carry1;
 721         h1 -= shl64(carry1, 25);
 722         carry5 = (h5 + (1L << 24)) >> 25;
 723         h6 += carry5;
 724         h5 -= shl64(carry5, 25);
 725
 726         carry2 = (h2 + (1L << 25)) >> 26;
 727         h3 += carry2;
 728         h2 -= shl64(carry2, 26);
 729         carry6 = (h6 + (1L << 25)) >> 26;
 730         h7 += carry6;
 731         h6 -= shl64(carry6, 26);
 732
 733         carry3 = (h3 + (1L << 24)) >> 25;
 734         h4 += carry3;
 735         h3 -= shl64(carry3, 25);
 736         carry7 = (h7 + (1L << 24)) >> 25;
 737         h8 += carry7;
 738         h7 -= shl64(carry7, 25);
 739
 740         carry4 = (h4 + (1L << 25)) >> 26;
 741         h5 += carry4;
 742         h4 -= shl64(carry4, 26);
 743         carry8 = (h8 + (1L << 25)) >> 26;
 744         h9 += carry8;
 745         h8 -= shl64(carry8, 26);
 746
 747         carry9 = (h9 + (1L << 24)) >> 25;
 748         h0 += carry9 * 19;
 749         h9 -= shl64(carry9, 25);
 750
 751         carry0 = (h0 + (1L << 25)) >> 26;
 752         h1 += carry0;
 753         h0 -= shl64(carry0, 26);
 754
 755         h[0] = (int32_t) h0;
 756         h[1] = (int32_t) h1;
 757         h[2] = (int32_t) h2;
 758         h[3] = (int32_t) h3;
 759         h[4] = (int32_t) h4;
 760         h[5] = (int32_t) h5;
 761         h[6] = (int32_t) h6;
 762         h[7] = (int32_t) h7;
 763         h[8] = (int32_t) h8;
 764         h[9] = (int32_t) h9;
 765 }
 766
 767
 768 /*
 769 h = f * 121666
 770 Can overlap h with f.
 771
 772 Preconditions:
 773    |f| bounded by 1.1*2^26,1.1*2^25,1.1*2^26,1.1*2^25,etc.
 774
 775 Postconditions:
 776    |h| bounded by 1.1*2^25,1.1*2^24,1.1*2^25,1.1*2^24,etc.
 777 */
 778
 779 void fe_mul121666(fe h, fe f) {
 780         int32_t f0 = f[0];
 781         int32_t f1 = f[1];
 782         int32_t f2 = f[2];
 783         int32_t f3 = f[3];
 784         int32_t f4 = f[4];
 785         int32_t f5 = f[5];
 786         int32_t f6 = f[6];
 787         int32_t f7 = f[7];
 788         int32_t f8 = f[8];
 789         int32_t f9 = f[9];
 790         int64_t h0 = f0 * (int64_t) 121666;
 791         int64_t h1 = f1 * (int64_t) 121666;
 792         int64_t h2 = f2 * (int64_t) 121666;
 793         int64_t h3 = f3 * (int64_t) 121666;
 794         int64_t h4 = f4 * (int64_t) 121666;
 795         int64_t h5 = f5 * (int64_t) 121666;
 796         int64_t h6 = f6 * (int64_t) 121666;
 797         int64_t h7 = f7 * (int64_t) 121666;
 798         int64_t h8 = f8 * (int64_t) 121666;
 799         int64_t h9 = f9 * (int64_t) 121666;
 800         int64_t carry0;
 801         int64_t carry1;
 802         int64_t carry2;
 803         int64_t carry3;
 804         int64_t carry4;
 805         int64_t carry5;
 806         int64_t carry6;
 807         int64_t carry7;
 808         int64_t carry8;
 809         int64_t carry9;
 810
 811         carry9 = (h9 + (int64_t)(1 << 24)) >> 25;
 812         h0 += carry9 * 19;
 813         h9 -= shl64(carry9, 25);
 814         carry1 = (h1 + (int64_t)(1 << 24)) >> 25;
 815         h2 += carry1;
 816         h1 -= shl64(carry1, 25);
 817         carry3 = (h3 + (int64_t)(1 << 24)) >> 25;
 818         h4 += carry3;
 819         h3 -= shl64(carry3, 25);
 820         carry5 = (h5 + (int64_t)(1 << 24)) >> 25;
 821         h6 += carry5;
 822         h5 -= shl64(carry5, 25);
 823         carry7 = (h7 + (int64_t)(1 << 24)) >> 25;
 824         h8 += carry7;
 825         h7 -= shl64(carry7, 25);
 826
 827         carry0 = (h0 + (int64_t)(1 << 25)) >> 26;
 828         h1 += carry0;
 829         h0 -= shl64(carry0, 26);
 830         carry2 = (h2 + (int64_t)(1 << 25)) >> 26;
 831         h3 += carry2;
 832         h2 -= shl64(carry2, 26);
 833         carry4 = (h4 + (int64_t)(1 << 25)) >> 26;
 834         h5 += carry4;
 835         h4 -= shl64(carry4, 26);
 836         carry6 = (h6 + (int64_t)(1 << 25)) >> 26;
 837         h7 += carry6;
 838         h6 -= shl64(carry6, 26);
 839         carry8 = (h8 + (int64_t)(1 << 25)) >> 26;
 840         h9 += carry8;
 841         h8 -= shl64(carry8, 26);
 842
 843         h[0] = h0;
 844         h[1] = h1;
 845         h[2] = h2;
 846         h[3] = h3;
 847         h[4] = h4;
 848         h[5] = h5;
 849         h[6] = h6;
 850         h[7] = h7;
 851         h[8] = h8;
 852         h[9] = h9;
 853 }
 854
 855
 856 /*
 857 h = -f
 858
 859 Preconditions:
 860    |f| bounded by 1.1*2^25,1.1*2^24,1.1*2^25,1.1*2^24,etc.
 861
 862 Postconditions:
 863    |h| bounded by 1.1*2^25,1.1*2^24,1.1*2^25,1.1*2^24,etc.
 864 */
 865
 866 void fe_neg(fe h, const fe f) {
 867         int32_t f0 = f[0];
 868         int32_t f1 = f[1];
 869         int32_t f2 = f[2];
 870         int32_t f3 = f[3];
 871         int32_t f4 = f[4];
 872         int32_t f5 = f[5];
 873         int32_t f6 = f[6];
 874         int32_t f7 = f[7];
 875         int32_t f8 = f[8];
 876         int32_t f9 = f[9];
 877         int32_t h0 = -f0;
 878         int32_t h1 = -f1;
 879         int32_t h2 = -f2;
 880         int32_t h3 = -f3;
 881         int32_t h4 = -f4;
 882         int32_t h5 = -f5;
 883         int32_t h6 = -f6;
 884         int32_t h7 = -f7;
 885         int32_t h8 = -f8;
 886         int32_t h9 = -f9;
 887
 888         h[0] = h0;
 889         h[1] = h1;
 890         h[2] = h2;
 891         h[3] = h3;
 892         h[4] = h4;
 893         h[5] = h5;
 894         h[6] = h6;
 895         h[7] = h7;
 896         h[8] = h8;
 897         h[9] = h9;
 898 }
 899
 900
 901 void fe_pow22523(fe out, const fe z) {
 902         fe t0;
 903         fe t1;
 904         fe t2;
 905         int i;
 906         fe_sq(t0, z);
 907
 908         for(i = 1; i < 1; ++i) {
 909                 fe_sq(t0, t0);
 910         }
 911
 912         fe_sq(t1, t0);
 913
 914         for(i = 1; i < 2; ++i) {
 915                 fe_sq(t1, t1);
 916         }
 917
 918         fe_mul(t1, z, t1);
 919         fe_mul(t0, t0, t1);
 920         fe_sq(t0, t0);
 921
 922         for(i = 1; i < 1; ++i) {
 923                 fe_sq(t0, t0);
 924         }
 925
 926         fe_mul(t0, t1, t0);
 927         fe_sq(t1, t0);
 928
 929         for(i = 1; i < 5; ++i) {
 930                 fe_sq(t1, t1);
 931         }
 932
 933         fe_mul(t0, t1, t0);
 934         fe_sq(t1, t0);
 935
 936         for(i = 1; i < 10; ++i) {
 937                 fe_sq(t1, t1);
 938         }
 939
 940         fe_mul(t1, t1, t0);
 941         fe_sq(t2, t1);
 942
 943         for(i = 1; i < 20; ++i) {
 944                 fe_sq(t2, t2);
 945         }
 946
 947         fe_mul(t1, t2, t1);
 948         fe_sq(t1, t1);
 949
 950         for(i = 1; i < 10; ++i) {
 951                 fe_sq(t1, t1);
 952         }
 953
 954         fe_mul(t0, t1, t0);
 955         fe_sq(t1, t0);
 956
 957         for(i = 1; i < 50; ++i) {
 958                 fe_sq(t1, t1);
 959         }
 960
 961         fe_mul(t1, t1, t0);
 962         fe_sq(t2, t1);
 963
 964         for(i = 1; i < 100; ++i) {
 965                 fe_sq(t2, t2);
 966         }
 967
 968         fe_mul(t1, t2, t1);
 969         fe_sq(t1, t1);
 970
 971         for(i = 1; i < 50; ++i) {
 972                 fe_sq(t1, t1);
 973         }
 974
 975         fe_mul(t0, t1, t0);
 976         fe_sq(t0, t0);
 977
 978         for(i = 1; i < 2; ++i) {
 979                 fe_sq(t0, t0);
 980         }
 981
 982         fe_mul(out, t0, z);
 983         return;
 984 }
 985
 986
 987 /*
 988 h = f * f
 989 Can overlap h with f.
 990
 991 Preconditions:
 992    |f| bounded by 1.65*2^26,1.65*2^25,1.65*2^26,1.65*2^25,etc.
 993
 994 Postconditions:
 995    |h| bounded by 1.01*2^25,1.01*2^24,1.01*2^25,1.01*2^24,etc.
 996 */
 997
 998 /*
 999 See fe_mul.c for discussion of implementation strategy.
1000 */
1001
1002 void fe_sq(fe h, const fe f) {
1003         int32_t f0 = f[0];
1004         int32_t f1 = f[1];
1005         int32_t f2 = f[2];
1006         int32_t f3 = f[3];
1007         int32_t f4 = f[4];
1008         int32_t f5 = f[5];
1009         int32_t f6 = f[6];
1010         int32_t f7 = f[7];
1011         int32_t f8 = f[8];
1012         int32_t f9 = f[9];
1013         int32_t f0_2 = 2 * f0;
1014         int32_t f1_2 = 2 * f1;
1015         int32_t f2_2 = 2 * f2;
1016         int32_t f3_2 = 2 * f3;
1017         int32_t f4_2 = 2 * f4;
1018         int32_t f5_2 = 2 * f5;
1019         int32_t f6_2 = 2 * f6;
1020         int32_t f7_2 = 2 * f7;
1021         int32_t f5_38 = 38 * f5; /* 1.959375*2^30 */
1022         int32_t f6_19 = 19 * f6; /* 1.959375*2^30 */
1023         int32_t f7_38 = 38 * f7; /* 1.959375*2^30 */
1024         int32_t f8_19 = 19 * f8; /* 1.959375*2^30 */
1025         int32_t f9_38 = 38 * f9; /* 1.959375*2^30 */
1026         int64_t f0f0    = f0   * (int64_t) f0;
1027         int64_t f0f1_2  = f0_2 * (int64_t) f1;
1028         int64_t f0f2_2  = f0_2 * (int64_t) f2;
1029         int64_t f0f3_2  = f0_2 * (int64_t) f3;
1030         int64_t f0f4_2  = f0_2 * (int64_t) f4;
1031         int64_t f0f5_2  = f0_2 * (int64_t) f5;
1032         int64_t f0f6_2  = f0_2 * (int64_t) f6;
1033         int64_t f0f7_2  = f0_2 * (int64_t) f7;
1034         int64_t f0f8_2  = f0_2 * (int64_t) f8;
1035         int64_t f0f9_2  = f0_2 * (int64_t) f9;
1036         int64_t f1f1_2  = f1_2 * (int64_t) f1;
1037         int64_t f1f2_2  = f1_2 * (int64_t) f2;
1038         int64_t f1f3_4  = f1_2 * (int64_t) f3_2;
1039         int64_t f1f4_2  = f1_2 * (int64_t) f4;
1040         int64_t f1f5_4  = f1_2 * (int64_t) f5_2;
1041         int64_t f1f6_2  = f1_2 * (int64_t) f6;
1042         int64_t f1f7_4  = f1_2 * (int64_t) f7_2;
1043         int64_t f1f8_2  = f1_2 * (int64_t) f8;
1044         int64_t f1f9_76 = f1_2 * (int64_t) f9_38;
1045         int64_t f2f2    = f2   * (int64_t) f2;
1046         int64_t f2f3_2  = f2_2 * (int64_t) f3;
1047         int64_t f2f4_2  = f2_2 * (int64_t) f4;
1048         int64_t f2f5_2  = f2_2 * (int64_t) f5;
1049         int64_t f2f6_2  = f2_2 * (int64_t) f6;
1050         int64_t f2f7_2  = f2_2 * (int64_t) f7;
1051         int64_t f2f8_38 = f2_2 * (int64_t) f8_19;
1052         int64_t f2f9_38 = f2   * (int64_t) f9_38;
1053         int64_t f3f3_2  = f3_2 * (int64_t) f3;
1054         int64_t f3f4_2  = f3_2 * (int64_t) f4;
1055         int64_t f3f5_4  = f3_2 * (int64_t) f5_2;
1056         int64_t f3f6_2  = f3_2 * (int64_t) f6;
1057         int64_t f3f7_76 = f3_2 * (int64_t) f7_38;
1058         int64_t f3f8_38 = f3_2 * (int64_t) f8_19;
1059         int64_t f3f9_76 = f3_2 * (int64_t) f9_38;
1060         int64_t f4f4    = f4   * (int64_t) f4;
1061         int64_t f4f5_2  = f4_2 * (int64_t) f5;
1062         int64_t f4f6_38 = f4_2 * (int64_t) f6_19;
1063         int64_t f4f7_38 = f4   * (int64_t) f7_38;
1064         int64_t f4f8_38 = f4_2 * (int64_t) f8_19;
1065         int64_t f4f9_38 = f4   * (int64_t) f9_38;
1066         int64_t f5f5_38 = f5   * (int64_t) f5_38;
1067         int64_t f5f6_38 = f5_2 * (int64_t) f6_19;
1068         int64_t f5f7_76 = f5_2 * (int64_t) f7_38;
1069         int64_t f5f8_38 = f5_2 * (int64_t) f8_19;
1070         int64_t f5f9_76 = f5_2 * (int64_t) f9_38;
1071         int64_t f6f6_19 = f6   * (int64_t) f6_19;
1072         int64_t f6f7_38 = f6   * (int64_t) f7_38;
1073         int64_t f6f8_38 = f6_2 * (int64_t) f8_19;
1074         int64_t f6f9_38 = f6   * (int64_t) f9_38;
1075         int64_t f7f7_38 = f7   * (int64_t) f7_38;
1076         int64_t f7f8_38 = f7_2 * (int64_t) f8_19;
1077         int64_t f7f9_76 = f7_2 * (int64_t) f9_38;
1078         int64_t f8f8_19 = f8   * (int64_t) f8_19;
1079         int64_t f8f9_38 = f8   * (int64_t) f9_38;
1080         int64_t f9f9_38 = f9   * (int64_t) f9_38;
1081         int64_t h0 = f0f0  + f1f9_76 + f2f8_38 + f3f7_76 + f4f6_38 + f5f5_38;
1082         int64_t h1 = f0f1_2 + f2f9_38 + f3f8_38 + f4f7_38 + f5f6_38;
1083         int64_t h2 = f0f2_2 + f1f1_2 + f3f9_76 + f4f8_38 + f5f7_76 + f6f6_19;
1084         int64_t h3 = f0f3_2 + f1f2_2 + f4f9_38 + f5f8_38 + f6f7_38;
1085         int64_t h4 = f0f4_2 + f1f3_4 + f2f2   + f5f9_76 + f6f8_38 + f7f7_38;
1086         int64_t h5 = f0f5_2 + f1f4_2 + f2f3_2 + f6f9_38 + f7f8_38;
1087         int64_t h6 = f0f6_2 + f1f5_4 + f2f4_2 + f3f3_2 + f7f9_76 + f8f8_19;
1088         int64_t h7 = f0f7_2 + f1f6_2 + f2f5_2 + f3f4_2 + f8f9_38;
1089         int64_t h8 = f0f8_2 + f1f7_4 + f2f6_2 + f3f5_4 + f4f4   + f9f9_38;
1090         int64_t h9 = f0f9_2 + f1f8_2 + f2f7_2 + f3f6_2 + f4f5_2;
1091         int64_t carry0;
1092         int64_t carry1;
1093         int64_t carry2;
1094         int64_t carry3;
1095         int64_t carry4;
1096         int64_t carry5;
1097         int64_t carry6;
1098         int64_t carry7;
1099         int64_t carry8;
1100         int64_t carry9;
1101         carry0 = (h0 + (1L << 25)) >> 26;
1102         h1 += carry0;
1103         h0 -= shl64(carry0, 26);
1104         carry4 = (h4 + (1L << 25)) >> 26;
1105         h5 += carry4;
1106         h4 -= shl64(carry4, 26);
1107         carry1 = (h1 + (1L << 24)) >> 25;
1108         h2 += carry1;
1109         h1 -= shl64(carry1, 25);
1110         carry5 = (h5 + (1L << 24)) >> 25;
1111         h6 += carry5;
1112         h5 -= shl64(carry5, 25);
1113         carry2 = (h2 + (1L << 25)) >> 26;
1114         h3 += carry2;
1115         h2 -= shl64(carry2, 26);
1116         carry6 = (h6 + (1L << 25)) >> 26;
1117         h7 += carry6;
1118         h6 -= shl64(carry6, 26);
1119         carry3 = (h3 + (1L << 24)) >> 25;
1120         h4 += carry3;
1121         h3 -= shl64(carry3, 25);
1122         carry7 = (h7 + (1L << 24)) >> 25;
1123         h8 += carry7;
1124         h7 -= shl64(carry7, 25);
1125         carry4 = (h4 + (1L << 25)) >> 26;
1126         h5 += carry4;
1127         h4 -= shl64(carry4, 26);
1128         carry8 = (h8 + (1L << 25)) >> 26;
1129         h9 += carry8;
1130         h8 -= shl64(carry8, 26);
1131         carry9 = (h9 + (1L << 24)) >> 25;
1132         h0 += carry9 * 19;
1133         h9 -= shl64(carry9, 25);
1134         carry0 = (h0 + (1L << 25)) >> 26;
1135         h1 += carry0;
1136         h0 -= shl64(carry0, 26);
1137         h[0] = (int32_t) h0;
1138         h[1] = (int32_t) h1;
1139         h[2] = (int32_t) h2;
1140         h[3] = (int32_t) h3;
1141         h[4] = (int32_t) h4;
1142         h[5] = (int32_t) h5;
1143         h[6] = (int32_t) h6;
1144         h[7] = (int32_t) h7;
1145         h[8] = (int32_t) h8;
1146         h[9] = (int32_t) h9;
1147 }
1148
1149
1150 /*
1151 h = 2 * f * f
1152 Can overlap h with f.
1153
1154 Preconditions:
1155    |f| bounded by 1.65*2^26,1.65*2^25,1.65*2^26,1.65*2^25,etc.
1156
1157 Postconditions:
1158    |h| bounded by 1.01*2^25,1.01*2^24,1.01*2^25,1.01*2^24,etc.
1159 */
1160
1161 /*
1162 See fe_mul.c for discussion of implementation strategy.
1163 */
1164
1165 void fe_sq2(fe h, const fe f) {
1166         int32_t f0 = f[0];
1167         int32_t f1 = f[1];
1168         int32_t f2 = f[2];
1169         int32_t f3 = f[3];
1170         int32_t f4 = f[4];
1171         int32_t f5 = f[5];
1172         int32_t f6 = f[6];
1173         int32_t f7 = f[7];
1174         int32_t f8 = f[8];
1175         int32_t f9 = f[9];
1176         int32_t f0_2 = 2 * f0;
1177         int32_t f1_2 = 2 * f1;
1178         int32_t f2_2 = 2 * f2;
1179         int32_t f3_2 = 2 * f3;
1180         int32_t f4_2 = 2 * f4;
1181         int32_t f5_2 = 2 * f5;
1182         int32_t f6_2 = 2 * f6;
1183         int32_t f7_2 = 2 * f7;
1184         int32_t f5_38 = 38 * f5; /* 1.959375*2^30 */
1185         int32_t f6_19 = 19 * f6; /* 1.959375*2^30 */
1186         int32_t f7_38 = 38 * f7; /* 1.959375*2^30 */
1187         int32_t f8_19 = 19 * f8; /* 1.959375*2^30 */
1188         int32_t f9_38 = 38 * f9; /* 1.959375*2^30 */
1189         int64_t f0f0    = f0   * (int64_t) f0;
1190         int64_t f0f1_2  = f0_2 * (int64_t) f1;
1191         int64_t f0f2_2  = f0_2 * (int64_t) f2;
1192         int64_t f0f3_2  = f0_2 * (int64_t) f3;
1193         int64_t f0f4_2  = f0_2 * (int64_t) f4;
1194         int64_t f0f5_2  = f0_2 * (int64_t) f5;
1195         int64_t f0f6_2  = f0_2 * (int64_t) f6;
1196         int64_t f0f7_2  = f0_2 * (int64_t) f7;
1197         int64_t f0f8_2  = f0_2 * (int64_t) f8;
1198         int64_t f0f9_2  = f0_2 * (int64_t) f9;
1199         int64_t f1f1_2  = f1_2 * (int64_t) f1;
1200         int64_t f1f2_2  = f1_2 * (int64_t) f2;
1201         int64_t f1f3_4  = f1_2 * (int64_t) f3_2;
1202         int64_t f1f4_2  = f1_2 * (int64_t) f4;
1203         int64_t f1f5_4  = f1_2 * (int64_t) f5_2;
1204         int64_t f1f6_2  = f1_2 * (int64_t) f6;
1205         int64_t f1f7_4  = f1_2 * (int64_t) f7_2;
1206         int64_t f1f8_2  = f1_2 * (int64_t) f8;
1207         int64_t f1f9_76 = f1_2 * (int64_t) f9_38;
1208         int64_t f2f2    = f2   * (int64_t) f2;
1209         int64_t f2f3_2  = f2_2 * (int64_t) f3;
1210         int64_t f2f4_2  = f2_2 * (int64_t) f4;
1211         int64_t f2f5_2  = f2_2 * (int64_t) f5;
1212         int64_t f2f6_2  = f2_2 * (int64_t) f6;
1213         int64_t f2f7_2  = f2_2 * (int64_t) f7;
1214         int64_t f2f8_38 = f2_2 * (int64_t) f8_19;
1215         int64_t f2f9_38 = f2   * (int64_t) f9_38;
1216         int64_t f3f3_2  = f3_2 * (int64_t) f3;
1217         int64_t f3f4_2  = f3_2 * (int64_t) f4;
1218         int64_t f3f5_4  = f3_2 * (int64_t) f5_2;
1219         int64_t f3f6_2  = f3_2 * (int64_t) f6;
1220         int64_t f3f7_76 = f3_2 * (int64_t) f7_38;
1221         int64_t f3f8_38 = f3_2 * (int64_t) f8_19;
1222         int64_t f3f9_76 = f3_2 * (int64_t) f9_38;
1223         int64_t f4f4    = f4   * (int64_t) f4;
1224         int64_t f4f5_2  = f4_2 * (int64_t) f5;
1225         int64_t f4f6_38 = f4_2 * (int64_t) f6_19;
1226         int64_t f4f7_38 = f4   * (int64_t) f7_38;
1227         int64_t f4f8_38 = f4_2 * (int64_t) f8_19;
1228         int64_t f4f9_38 = f4   * (int64_t) f9_38;
1229         int64_t f5f5_38 = f5   * (int64_t) f5_38;
1230         int64_t f5f6_38 = f5_2 * (int64_t) f6_19;
1231         int64_t f5f7_76 = f5_2 * (int64_t) f7_38;
1232         int64_t f5f8_38 = f5_2 * (int64_t) f8_19;
1233         int64_t f5f9_76 = f5_2 * (int64_t) f9_38;
1234         int64_t f6f6_19 = f6   * (int64_t) f6_19;
1235         int64_t f6f7_38 = f6   * (int64_t) f7_38;
1236         int64_t f6f8_38 = f6_2 * (int64_t) f8_19;
1237         int64_t f6f9_38 = f6   * (int64_t) f9_38;
1238         int64_t f7f7_38 = f7   * (int64_t) f7_38;
1239         int64_t f7f8_38 = f7_2 * (int64_t) f8_19;
1240         int64_t f7f9_76 = f7_2 * (int64_t) f9_38;
1241         int64_t f8f8_19 = f8   * (int64_t) f8_19;
1242         int64_t f8f9_38 = f8   * (int64_t) f9_38;
1243         int64_t f9f9_38 = f9   * (int64_t) f9_38;
1244         int64_t h0 = f0f0  + f1f9_76 + f2f8_38 + f3f7_76 + f4f6_38 + f5f5_38;
1245         int64_t h1 = f0f1_2 + f2f9_38 + f3f8_38 + f4f7_38 + f5f6_38;
1246         int64_t h2 = f0f2_2 + f1f1_2 + f3f9_76 + f4f8_38 + f5f7_76 + f6f6_19;
1247         int64_t h3 = f0f3_2 + f1f2_2 + f4f9_38 + f5f8_38 + f6f7_38;
1248         int64_t h4 = f0f4_2 + f1f3_4 + f2f2   + f5f9_76 + f6f8_38 + f7f7_38;
1249         int64_t h5 = f0f5_2 + f1f4_2 + f2f3_2 + f6f9_38 + f7f8_38;
1250         int64_t h6 = f0f6_2 + f1f5_4 + f2f4_2 + f3f3_2 + f7f9_76 + f8f8_19;
1251         int64_t h7 = f0f7_2 + f1f6_2 + f2f5_2 + f3f4_2 + f8f9_38;
1252         int64_t h8 = f0f8_2 + f1f7_4 + f2f6_2 + f3f5_4 + f4f4   + f9f9_38;
1253         int64_t h9 = f0f9_2 + f1f8_2 + f2f7_2 + f3f6_2 + f4f5_2;
1254         int64_t carry0;
1255         int64_t carry1;
1256         int64_t carry2;
1257         int64_t carry3;
1258         int64_t carry4;
1259         int64_t carry5;
1260         int64_t carry6;
1261         int64_t carry7;
1262         int64_t carry8;
1263         int64_t carry9;
1264         h0 += h0;
1265         h1 += h1;
1266         h2 += h2;
1267         h3 += h3;
1268         h4 += h4;
1269         h5 += h5;
1270         h6 += h6;
1271         h7 += h7;
1272         h8 += h8;
1273         h9 += h9;
1274         carry0 = (h0 + (1L << 25)) >> 26;
1275         h1 += carry0;
1276         h0 -= shl64(carry0, 26);
1277         carry4 = (h4 + (1L << 25)) >> 26;
1278         h5 += carry4;
1279         h4 -= shl64(carry4, 26);
1280         carry1 = (h1 + (1L << 24)) >> 25;
1281         h2 += carry1;
1282         h1 -= shl64(carry1, 25);
1283         carry5 = (h5 + (1L << 24)) >> 25;
1284         h6 += carry5;
1285         h5 -= shl64(carry5, 25);
1286         carry2 = (h2 + (1L << 25)) >> 26;
1287         h3 += carry2;
1288         h2 -= shl64(carry2, 26);
1289         carry6 = (h6 + (1L << 25)) >> 26;
1290         h7 += carry6;
1291         h6 -= shl64(carry6, 26);
1292         carry3 = (h3 + (1L << 24)) >> 25;
1293         h4 += carry3;
1294         h3 -= shl64(carry3, 25);
1295         carry7 = (h7 + (1L << 24)) >> 25;
1296         h8 += carry7;
1297         h7 -= shl64(carry7, 25);
1298         carry4 = (h4 + (1L << 25)) >> 26;
1299         h5 += carry4;
1300         h4 -= shl64(carry4, 26);
1301         carry8 = (h8 + (1L << 25)) >> 26;
1302         h9 += carry8;
1303         h8 -= shl64(carry8, 26);
1304         carry9 = (h9 + (1L << 24)) >> 25;
1305         h0 += carry9 * 19;
1306         h9 -= shl64(carry9, 25);
1307         carry0 = (h0 + (1L << 25)) >> 26;
1308         h1 += carry0;
1309         h0 -= shl64(carry0, 26);
1310         h[0] = (int32_t) h0;
1311         h[1] = (int32_t) h1;
1312         h[2] = (int32_t) h2;
1313         h[3] = (int32_t) h3;
1314         h[4] = (int32_t) h4;
1315         h[5] = (int32_t) h5;
1316         h[6] = (int32_t) h6;
1317         h[7] = (int32_t) h7;
1318         h[8] = (int32_t) h8;
1319         h[9] = (int32_t) h9;
1320 }
1321
1322
1323 /*
1324 h = f - g
1325 Can overlap h with f or g.
1326
1327 Preconditions:
1328    |f| bounded by 1.1*2^25,1.1*2^24,1.1*2^25,1.1*2^24,etc.
1329    |g| bounded by 1.1*2^25,1.1*2^24,1.1*2^25,1.1*2^24,etc.
1330
1331 Postconditions:
1332    |h| bounded by 1.1*2^26,1.1*2^25,1.1*2^26,1.1*2^25,etc.
1333 */
1334
1335 void fe_sub(fe h, const fe f, const fe g) {
1336         int32_t f0 = f[0];
1337         int32_t f1 = f[1];
1338         int32_t f2 = f[2];
1339         int32_t f3 = f[3];
1340         int32_t f4 = f[4];
1341         int32_t f5 = f[5];
1342         int32_t f6 = f[6];
1343         int32_t f7 = f[7];
1344         int32_t f8 = f[8];
1345         int32_t f9 = f[9];
1346         int32_t g0 = g[0];
1347         int32_t g1 = g[1];
1348         int32_t g2 = g[2];
1349         int32_t g3 = g[3];
1350         int32_t g4 = g[4];
1351         int32_t g5 = g[5];
1352         int32_t g6 = g[6];
1353         int32_t g7 = g[7];
1354         int32_t g8 = g[8];
1355         int32_t g9 = g[9];
1356         int32_t h0 = f0 - g0;
1357         int32_t h1 = f1 - g1;
1358         int32_t h2 = f2 - g2;
1359         int32_t h3 = f3 - g3;
1360         int32_t h4 = f4 - g4;
1361         int32_t h5 = f5 - g5;
1362         int32_t h6 = f6 - g6;
1363         int32_t h7 = f7 - g7;
1364         int32_t h8 = f8 - g8;
1365         int32_t h9 = f9 - g9;
1366
1367         h[0] = h0;
1368         h[1] = h1;
1369         h[2] = h2;
1370         h[3] = h3;
1371         h[4] = h4;
1372         h[5] = h5;
1373         h[6] = h6;
1374         h[7] = h7;
1375         h[8] = h8;
1376         h[9] = h9;
1377 }
1378
1379
1380
1381 /*
1382 Preconditions:
1383   |h| bounded by 1.1*2^26,1.1*2^25,1.1*2^26,1.1*2^25,etc.
1384
1385 Write p=2^255-19; q=floor(h/p).
1386 Basic claim: q = floor(2^(-255)(h + 19 2^(-25)h9 + 2^(-1))).
1387
1388 Proof:
1389   Have |h|<=p so |q|<=1 so |19^2 2^(-255) q|<1/4.
1390   Also have |h-2^230 h9|<2^231 so |19 2^(-255)(h-2^230 h9)|<1/4.
1391
1392   Write y=2^(-1)-19^2 2^(-255)q-19 2^(-255)(h-2^230 h9).
1393   Then 0<y<1.
1394
1395   Write r=h-pq.
1396   Have 0<=r<=p-1=2^255-20.
1397   Thus 0<=r+19(2^-255)r<r+19(2^-255)2^255<=2^255-1.
1398
1399   Write x=r+19(2^-255)r+y.
1400   Then 0<x<2^255 so floor(2^(-255)x) = 0 so floor(q+2^(-255)x) = q.
1401
1402   Have q+2^(-255)x = 2^(-255)(h + 19 2^(-25) h9 + 2^(-1))
1403   so floor(2^(-255)(h + 19 2^(-25) h9 + 2^(-1))) = q.
1404 */
1405
1406 void fe_tobytes(unsigned char *s, const fe h) {
1407         int32_t h0 = h[0];
1408         int32_t h1 = h[1];
1409         int32_t h2 = h[2];
1410         int32_t h3 = h[3];
1411         int32_t h4 = h[4];
1412         int32_t h5 = h[5];
1413         int32_t h6 = h[6];
1414         int32_t h7 = h[7];
1415         int32_t h8 = h[8];
1416         int32_t h9 = h[9];
1417         int32_t q;
1418         int32_t carry0;
1419         int32_t carry1;
1420         int32_t carry2;
1421         int32_t carry3;
1422         int32_t carry4;
1423         int32_t carry5;
1424         int32_t carry6;
1425         int32_t carry7;
1426         int32_t carry8;
1427         int32_t carry9;
1428         q = (19 * h9 + (((int32_t) 1) << 24)) >> 25;
1429         q = (h0 + q) >> 26;
1430         q = (h1 + q) >> 25;
1431         q = (h2 + q) >> 26;
1432         q = (h3 + q) >> 25;
1433         q = (h4 + q) >> 26;
1434         q = (h5 + q) >> 25;
1435         q = (h6 + q) >> 26;
1436         q = (h7 + q) >> 25;
1437         q = (h8 + q) >> 26;
1438         q = (h9 + q) >> 25;
1439         /* Goal: Output h-(2^255-19)q, which is between 0 and 2^255-20. */
1440         h0 += 19 * q;
1441         /* Goal: Output h-2^255 q, which is between 0 and 2^255-20. */
1442         carry0 = h0 >> 26;
1443         h1 += carry0;
1444         h0 -= shl32(carry0, 26);
1445         carry1 = h1 >> 25;
1446         h2 += carry1;
1447         h1 -= shl32(carry1, 25);
1448         carry2 = h2 >> 26;
1449         h3 += carry2;
1450         h2 -= shl32(carry2, 26);
1451         carry3 = h3 >> 25;
1452         h4 += carry3;
1453         h3 -= shl32(carry3, 25);
1454         carry4 = h4 >> 26;
1455         h5 += carry4;
1456         h4 -= shl32(carry4, 26);
1457         carry5 = h5 >> 25;
1458         h6 += carry5;
1459         h5 -= shl32(carry5, 25);
1460         carry6 = h6 >> 26;
1461         h7 += carry6;
1462         h6 -= shl32(carry6, 26);
1463         carry7 = h7 >> 25;
1464         h8 += carry7;
1465         h7 -= shl32(carry7, 25);
1466         carry8 = h8 >> 26;
1467         h9 += carry8;
1468         h8 -= shl32(carry8, 26);
1469         carry9 = h9 >> 25;
1470         h9 -= shl32(carry9, 25);
1471
1472         /* h10 = carry9 */
1473         /*
1474         Goal: Output h0+...+2^255 h10-2^255 q, which is between 0 and 2^255-20.
1475         Have h0+...+2^230 h9 between 0 and 2^255-1;
1476         evidently 2^255 h10-2^255 q = 0.
1477         Goal: Output h0+...+2^230 h9.
1478         */
1479         s[0] = (unsigned char)(h0 >> 0);
1480         s[1] = (unsigned char)(h0 >> 8);
1481         s[2] = (unsigned char)(h0 >> 16);
1482         s[3] = (unsigned char)((h0 >> 24) | shl32(h1, 2));
1483         s[4] = (unsigned char)(h1 >> 6);
1484         s[5] = (unsigned char)(h1 >> 14);
1485         s[6] = (unsigned char)((h1 >> 22) | shl32(h2, 3));
1486         s[7] = (unsigned char)(h2 >> 5);
1487         s[8] = (unsigned char)(h2 >> 13);
1488         s[9] = (unsigned char)((h2 >> 21) | shl32(h3, 5));
1489         s[10] = (unsigned char)(h3 >> 3);
1490         s[11] = (unsigned char)(h3 >> 11);
1491         s[12] = (unsigned char)((h3 >> 19) | shl32(h4, 6));
1492         s[13] = (unsigned char)(h4 >> 2);
1493         s[14] = (unsigned char)(h4 >> 10);
1494         s[15] = (unsigned char)(h4 >> 18);
1495         s[16] = (unsigned char)(h5 >> 0);
1496         s[17] = (unsigned char)(h5 >> 8);
1497         s[18] = (unsigned char)(h5 >> 16);
1498         s[19] = (unsigned char)((h5 >> 24) | shl32(h6, 1));
1499         s[20] = (unsigned char)(h6 >> 7);
1500         s[21] = (unsigned char)(h6 >> 15);
1501         s[22] = (unsigned char)((h6 >> 23) | shl32(h7, 3));
1502         s[23] = (unsigned char)(h7 >> 5);
1503         s[24] = (unsigned char)(h7 >> 13);
1504         s[25] = (unsigned char)((h7 >> 21) | shl32(h8, 4));
1505         s[26] = (unsigned char)(h8 >> 4);
1506         s[27] = (unsigned char)(h8 >> 12);
1507         s[28] = (unsigned char)((h8 >> 20) | shl32(h9, 6));
1508         s[29] = (unsigned char)(h9 >> 2);
1509         s[30] = (unsigned char)(h9 >> 10);
1510         s[31] = (unsigned char)(h9 >> 18);
1511 }