pd.if.org Git - zpackage/blob - crypto/libeddsa/lib/fld.c

   1 /*
   2  * code for the field GF(2^255-19).
   3  *
   4  * This code is in public domain.
   5  *
   6  * Philipp Lay <philipp.lay@illunis.net>
   7  */
   8
   9 #include "bitness.h"
  10 #include "fld.h"
  11
  12
  13 #ifdef USE_64BIT
  14
  15 /*
  16  * 64bit implementation
  17  */
  18
  19 /*
  20  * exported field constants
  21  */
  22
  23 /* con_d = - 121665/121666 (mod q) */
  24 const fld_t con_d = {
  25         929955233495203, 466365720129213, 1662059464998953, 2033849074728123,
  26         1442794654840575 };
  27
  28 /* con_2d = 2 * con_d (mod q) */
  29 const fld_t con_2d = {
  30         1859910466990425, 932731440258426, 1072319116312658, 1815898335770999,
  31         633789495995903 };
  32
  33 /* con_m2d = -2 * con_d (mod q) */
  34 const fld_t con_m2d = {
  35         391889346694804, 1319068373426821, 1179480697372589, 435901477914248,
  36         1618010317689344 };
  37
  38 /* con_j^2 = 1 (mod q) */
  39 const fld_t con_j = {
  40         1718705420411056, 234908883556509, 2233514472574048, 2117202627021982,
  41         765476049583133 };
  42
  43
  44
  45 /*
  46  * fld_reduce - returns the smallest non-negative representation of x modulo q
  47  *              with 0 <= x[i] <= 2^51 - 1 for 0 <= i <= 4.
  48  *
  49  * This function assumes
  50  *   abs(x[i]) <= 2^63 - 2^12 = 2^12 * (2^51 - 1),
  51  * for 0 <= i <= 4 to work properly.
  52  */
  53 void
  54 fld_reduce(fld_t res, const fld_t x)
  55 {
  56         /* first carry round with offset 19 */
  57         res[0] = x[0] + 19;
  58         res[1] = x[1] + (res[0] >> FLD_LIMB_BITS);
  59         res[2] = x[2] + (res[1] >> FLD_LIMB_BITS);
  60         res[3] = x[3] + (res[2] >> FLD_LIMB_BITS);
  61         res[4] = x[4] + (res[3] >> FLD_LIMB_BITS);
  62
  63         /* -2^12 <= (res[4] >> FLD_LIMB_BITS) <= 2^12-1 */
  64         res[0] = (res[0] & FLD_LIMB_MASK) + 19*(res[4] >> FLD_LIMB_BITS);
  65         res[1] &= FLD_LIMB_MASK;
  66         res[2] &= FLD_LIMB_MASK;
  67         res[3] &= FLD_LIMB_MASK;
  68         res[4] &= FLD_LIMB_MASK;
  69
  70         /* now we have
  71          *   -19*2^12 <= res[0] <= 2^51-1 + 19*(2^12-1),
  72          * and
  73          *   0 <= res[i] <= 2^51-1 for 1 <= i <= 4.
  74          */
  75
  76         /* second round */
  77         res[1] += (res[0] >> FLD_LIMB_BITS);
  78         res[2] += (res[1] >> FLD_LIMB_BITS);
  79         res[3] += (res[2] >> FLD_LIMB_BITS);
  80         res[4] += (res[3] >> FLD_LIMB_BITS);
  81
  82         /* -1 <= (res[4] >> FLD_LIMB_BITS) <= 1 */
  83         res[0] = (res[0] & FLD_LIMB_MASK) + 19*(res[4] >> FLD_LIMB_BITS);
  84
  85         res[1] &= FLD_LIMB_MASK;
  86         res[2] &= FLD_LIMB_MASK;
  87         res[3] &= FLD_LIMB_MASK;
  88         res[4] &= FLD_LIMB_MASK;
  89
  90         /* the second round yields
  91          *   -19 <= res[0] <= 2^51-1 + 19
  92          * and
  93          *   0 <= res[i] <= 2^51-1 for 1 <= i <= 4.
  94          */
  95         res[0] -= 19;
  96
  97         /* with the offset removed we have
  98          *   -38 <= res[0] <= 2^51-1
  99          * and need a third round to assure that res[0] is non-negative.
 100          *
 101          * please note, that no positive carry is possible at this point.
 102          */
 103
 104         res[1] += (res[0] >> FLD_LIMB_BITS);
 105         res[2] += (res[1] >> FLD_LIMB_BITS);
 106         res[3] += (res[2] >> FLD_LIMB_BITS);
 107         res[4] += (res[3] >> FLD_LIMB_BITS);
 108
 109         /* -1 <= (res[4] >> FLD_LIMB_BITS) <= 0, because of non-positive carry */
 110         res[0] = (res[0] & FLD_LIMB_MASK) + 19*(res[4] >> FLD_LIMB_BITS);
 111
 112         res[1] &= FLD_LIMB_MASK;
 113         res[2] &= FLD_LIMB_MASK;
 114         res[3] &= FLD_LIMB_MASK;
 115         res[4] &= FLD_LIMB_MASK;
 116
 117
 118         /* if res[0] was negative before this round we had an negative carry and
 119          * have now
 120          *   2^51 - 38 - 19 <= res[0] <= 2^51 - 1.
 121          *
 122          * so in any case it is
 123          *   0 <= res[0] <= 2^51 - 1
 124          * and
 125          *   0 <= res[i] <= 2^51 - 1 for 1 <= i <= 4
 126          * as wished.
 127          *
 128          * moreover res is the smallest non-negative representant of x modulo q.
 129          */
 130 }
 131
 132 /*
 133  * fld_import - import an 256bit, unsigned, little-endian integer into
 134  * our signed 51-bit limb format and reduce modulo 2^255-19.
 135  */
 136 void
 137 fld_import(fld_t dst, const uint8_t src[32])
 138 {
 139         int i;
 140         uint64_t tmp = 0;
 141         int fill = 0;
 142
 143         for (i = 0; i < FLD_LIMB_NUM; i++) {
 144                 for (;fill < FLD_LIMB_BITS; fill += 8)
 145                         tmp |= (uint64_t)*src++ << fill;
 146
 147                 dst[i] = tmp & FLD_LIMB_MASK;
 148                 tmp >>= FLD_LIMB_BITS;
 149                 fill -= FLD_LIMB_BITS;
 150         }
 151         dst[0] += 19*tmp;
 152
 153         /* dst is now reduced and partially carried (first limb may
 154          * use 52bits instead of 51).
 155          */
 156 }
 157
 158 /*
 159  * fld_export - export our internal format to a 256bit, unsigned,
 160  * little-endian packed format.
 161  */
 162 void
 163 fld_export(uint8_t dst[32], const fld_t src)
 164 {
 165         int i;
 166         fld_t tmp;
 167         uint64_t foo = 0;
 168         int fill = 0;
 169
 170         fld_reduce(tmp, src);
 171
 172         for (i = 0; i < FLD_LIMB_NUM; i++) {
 173                 foo |= (uint64_t)tmp[i] << fill;
 174                 for (fill += FLD_LIMB_BITS; fill >= 8; fill -= 8, foo >>= 8)
 175                         *dst++ = foo & 0xff;
 176         }
 177         *dst++ = foo & 0xff;
 178 }
 179
 180 /*
 181  * fld_scale - multiply e by scalar s and reduce modulo q.
 182  */
 183 void
 184 fld_scale(fld_t res, const fld_t e, limb_t s)
 185 {
 186         llimb_t carry;
 187
 188         carry = (llimb_t)s*e[0];
 189         res[0] = carry & FLD_LIMB_MASK;
 190
 191         carry = (carry >> FLD_LIMB_BITS) + (llimb_t)s*e[1];
 192         res[1] = carry & FLD_LIMB_MASK;
 193
 194         carry = (carry >> FLD_LIMB_BITS) + (llimb_t)s*e[2];
 195         res[2] = carry & FLD_LIMB_MASK;
 196
 197         carry = (carry >> FLD_LIMB_BITS) + (llimb_t)s*e[3];
 198         res[3] = carry & FLD_LIMB_MASK;
 199
 200         carry = (carry >> FLD_LIMB_BITS) + (llimb_t)s*e[4];
 201         res[4] = carry & FLD_LIMB_MASK;
 202
 203         res[0] += 19*(carry >> FLD_LIMB_BITS);
 204 }
 205
 206 /*
 207  * fld_mul - multiply a with b and reduce modulo q.
 208  */
 209 void
 210 fld_mul(fld_t res, const fld_t a, const fld_t b)
 211 {
 212         limb_t a19_1, a19_2, a19_3, a19_4, tmp;
 213         llimb_t c[5];
 214
 215         a19_1 = 19*a[1];
 216         a19_2 = 19*a[2];
 217         a19_3 = 19*a[3];
 218         a19_4 = 19*a[4];
 219
 220         c[0] = (llimb_t)a[0]*b[0] + (llimb_t)a19_1*b[4] + (llimb_t)a19_2*b[3]
 221                 + (llimb_t)a19_3*b[2] + (llimb_t)a19_4*b[1];
 222         c[1] = (llimb_t)a[0]*b[1] + (llimb_t)a[1]*b[0] + (llimb_t)a19_2*b[4]
 223                 + (llimb_t)a19_3*b[3] + (llimb_t)a19_4*b[2];
 224         c[2] = (llimb_t)a[0]*b[2] + (llimb_t)a[1]*b[1] + (llimb_t)a[2]*b[0]
 225                 + (llimb_t)a19_3*b[4] + (llimb_t)a19_4*b[3];
 226         c[3] = (llimb_t)a[0]*b[3] + (llimb_t)a[1]*b[2] + (llimb_t)a[2]*b[1]
 227                 + (llimb_t)a[3]*b[0] + (llimb_t)a19_4*b[4];
 228         c[4] = (llimb_t)a[0]*b[4] + (llimb_t)a[1]*b[3] + (llimb_t)a[2]*b[2]
 229                 + (llimb_t)a[3]*b[1] + (llimb_t)a[4]*b[0];
 230
 231
 232         c[1] += c[0] >> FLD_LIMB_BITS;
 233         c[2] += c[1] >> FLD_LIMB_BITS;
 234         c[3] += c[2] >> FLD_LIMB_BITS;
 235         c[4] += c[3] >> FLD_LIMB_BITS;
 236
 237         tmp = (c[0] & FLD_LIMB_MASK) + 19*(c[4] >> FLD_LIMB_BITS);
 238         res[1] = (c[1] & FLD_LIMB_MASK) + (tmp >> FLD_LIMB_BITS);
 239
 240         res[0] = tmp & FLD_LIMB_MASK;
 241         res[2] = c[2] & FLD_LIMB_MASK;
 242         res[3] = c[3] & FLD_LIMB_MASK;
 243         res[4] = c[4] & FLD_LIMB_MASK;
 244 }
 245
 246 /*
 247  * fld_sq - square x and reduce modulo q.
 248  */
 249 void
 250 fld_sq(fld_t res, const fld_t x)
 251 {
 252         limb_t x2_1, x2_2, x2_3, x2_4, x19_3, x19_4, tmp;
 253         llimb_t c[5];
 254
 255         x2_1 = 2*x[1];
 256         x2_2 = 2*x[2];
 257         x2_3 = 2*x[3];
 258         x2_4 = 2*x[4];
 259         x19_3 = 19*x[3];
 260         x19_4 = 19*x[4];
 261
 262         c[0] = (llimb_t)x[0]*x[0] + (llimb_t)x2_1*x19_4 + (llimb_t)x2_2*x19_3;
 263         c[1] = (llimb_t)x[0]*x2_1 + (llimb_t)x2_2*x19_4 + (llimb_t)x19_3*x[3];
 264         c[2] = (llimb_t)x[0]*x2_2 + (llimb_t)x[1]*x[1] + (llimb_t)x2_3*x19_4;
 265         c[3] = (llimb_t)x[0]*x2_3 + (llimb_t)x2_1*x[2] + (llimb_t)x19_4*x[4];
 266         c[4] = (llimb_t)x[0]*x2_4 + (llimb_t)x2_1*x[3] + (llimb_t)x[2]*x[2];
 267
 268         c[1] += c[0] >> FLD_LIMB_BITS;
 269         c[2] += c[1] >> FLD_LIMB_BITS;
 270         c[3] += c[2] >> FLD_LIMB_BITS;
 271         c[4] += c[3] >> FLD_LIMB_BITS;
 272
 273         tmp = (c[0] & FLD_LIMB_MASK) + 19*(c[4] >> FLD_LIMB_BITS);
 274         res[1] = (c[1] & FLD_LIMB_MASK) + (tmp >> FLD_LIMB_BITS);
 275
 276         res[0] = tmp & FLD_LIMB_MASK;
 277         res[2] = c[2] & FLD_LIMB_MASK;
 278         res[3] = c[3] & FLD_LIMB_MASK;
 279         res[4] = c[4] & FLD_LIMB_MASK;
 280 }
 281
 282 #else           /* USE_64BIT */
 283
 284
 285 /*
 286  * 32bit implementation
 287  */
 288
 289
 290 /*
 291  * exported field constants
 292  */
 293
 294 /* con_d = - 121665/121666 (mod q) */
 295 const fld_t con_d = { 56195235, 13857412, 51736253, 6949390, 114729, 24766616,
 296                       60832955, 30306712, 48412415, 21499315 };
 297
 298 /* con_2d = 2 * con_d (mod q) */
 299 const fld_t con_2d = { 45281625, 27714825, 36363642, 13898781, 229458,
 300                        15978800, 54557047, 27058993, 29715967, 9444199 };
 301 /* con_m2d = -2 * con_d (mod q) */
 302 const fld_t con_m2d = { 21827220, 5839606, 30745221, 19655650, 66879405,
 303                         17575631, 12551816, 6495438, 37392896, 24110232 };
 304 /* j^2 = 1 (mod q) */
 305 const fld_t con_j = { 34513072, 25610706, 9377949, 3500415, 12389472, 33281959,
 306                       41962654, 31548777, 326685, 11406482 };
 307
 308
 309 /*
 310  * macro for doing one carry-reduce round
 311  *
 312  * tmp is being used as carry-register and could be of type
 313  * limb_t or llimb_t.
 314  *
 315  * off is normally zero, but could be used to wrap-around values x
 316  * with q <= x < 2^255 (see fld_reduce for an example).
 317  */
 318 #define CARRY(dst, src, tmp, off)                                       \
 319 do {                                                                    \
 320         int _ii;                                                        \
 321         (tmp) = (off);                                                  \
 322         (tmp) <<= FLD_LIMB_BITS(0);                                     \
 323         for (_ii = 0; _ii < FLD_LIMB_NUM; _ii += 2) {                   \
 324                 (tmp) = ((tmp) >> FLD_LIMB_BITS(0)) + (src)[_ii];       \
 325                 (dst)[_ii] = (tmp) & FLD_LIMB_MASK(1);                  \
 326                 (tmp) = ((tmp) >> FLD_LIMB_BITS(1)) + (src)[_ii+1];     \
 327                 (dst)[_ii+1] = (tmp) & FLD_LIMB_MASK(0);                \
 328         }                                                               \
 329         (dst)[0] += 19*((tmp) >> FLD_LIMB_BITS(0));                     \
 330 } while(0)
 331
 332
 333 /*
 334  * fld_reduce - returns the smallest non-negative representation of x modulo q
 335  *              with 0 <= x[i] <= 2^26 - 1 for i % 2 == 0
 336  *              and  0 <= x[i] <= 2^25 - 1 for i % 2 == 1.
 337  *
 338  * assumes:
 339  *   abs(x[i]) <= 2^31 - 2^5 = 32 * (2^26 - 1)
 340  */
 341 void
 342 fld_reduce(fld_t res, const fld_t x)
 343 {
 344         limb_t tmp;
 345
 346         CARRY(res, x, tmp, 19);
 347         /* we have
 348          *   -2^26 - 19*2^5 <= res[0] <= 2^26 - 1 + 19*(2^5 - 1)
 349          * and for i >= 1:
 350          *   0 <= res[i] <= 2^26-1, i % 2 == 0,
 351          *   0 <= res[i] <= 2^25-1, i % 2 == 1.
 352          */
 353
 354         CARRY(res, res, tmp, 0);
 355         /* now
 356          *   -2^26 - 19 <= res[0] <= 2^26 - 1 + 19
 357          * holds and, as before,
 358          *   0 <= res[i] <= 2^26-1, i % 2 == 0
 359          *   0 <= res[i] <= 2^25-1, i % 2 == 1,
 360          * for i >= 1.
 361          */
 362
 363         /* next round we will first remove our offset resulting in
 364          *   -2^26 - 38 <= res[0] <= 2^26 - 1,
 365          * therefor only a negative carry could appear.
 366          */
 367         CARRY(res, res, tmp, -19);
 368         /* now we have
 369          *   0 <= res[i] <= 2^26-1, i % 2 == 0
 370          *   0 <= res[i] <= 2^25-1, i % 2 == 1
 371          * for all limbs as wished.
 372          *
 373          * if a carry had happend, we even know
 374          *   2^26 - 38 - 19 <= res[0] <= 2^26 - 1.
 375          */
 376 }
 377
 378 /*
 379  * fld_import - import an 256bit, unsigned, little-endian integer into
 380  * our internal fld_t format.
 381  */
 382 void
 383 fld_import(fld_t dst, const uint8_t src[32])
 384 {
 385         int i;
 386         uint32_t foo = 0;
 387         int fill = 0;
 388         int d = 1;
 389
 390         for (i = 0; i < FLD_LIMB_NUM; i++) {
 391                 for (; fill < FLD_LIMB_BITS(d); fill += 8)
 392                         foo |= (uint32_t)*src++ << fill;
 393                 dst[i] = foo & FLD_LIMB_MASK(d);
 394
 395                 foo >>= FLD_LIMB_BITS(d);
 396                 fill -= FLD_LIMB_BITS(d);
 397                 d = 1-d;
 398         }
 399         dst[0] += 19*foo;
 400 }
 401
 402 /*
 403  * fld_export - export a field element into 256bit little-endian encoded form.
 404  */
 405 void
 406 fld_export(uint8_t dst[32], const fld_t src)
 407 {
 408         uint32_t foo;
 409         fld_t tmp;
 410         int fill, i;
 411
 412         fld_reduce(tmp, src);
 413
 414         for (i = 0, fill = 0, foo = 0; i < FLD_LIMB_NUM; i += 2) {
 415                 foo |= (tmp[i] & FLD_LIMB_MASK(1)) << fill;
 416                 for (fill += FLD_LIMB_BITS(1); fill >= 8; fill -= 8, foo >>= 8)
 417                         *dst++ = foo & 0xff;
 418
 419                 foo |= (tmp[i+1] & FLD_LIMB_MASK(0)) << fill;
 420                 for (fill += FLD_LIMB_BITS(0); fill >= 8; fill -= 8, foo >>= 8)
 421                         *dst++ = foo & 0xff;
 422         }
 423         *dst++ = foo & 0xff;
 424 }
 425
 426 /*
 427  * fld_scale - multiply e by scalar s and reduce modulo q.
 428  */
 429 void
 430 fld_scale(fld_t dst, const fld_t e, limb_t x)
 431 {
 432         llimb_t tmp;
 433         int i;
 434
 435         for (tmp = 0, i = 0; i < FLD_LIMB_NUM; i += 2) {
 436                 tmp = (tmp >> FLD_LIMB_BITS(0)) + (llimb_t) x*e[i];
 437                 dst[i] = tmp & FLD_LIMB_MASK(1);
 438                 tmp = (tmp >> FLD_LIMB_BITS(1)) + ((llimb_t) x*e[i+1]);
 439                 dst[i+1] = tmp & FLD_LIMB_MASK(0);
 440         }
 441         dst[0] += 19*(tmp >> FLD_LIMB_BITS(0));
 442 }
 443
 444 /*
 445  * fld_mul - multiply a with b and reduce modulo q.
 446  */
 447 void
 448 fld_mul(fld_t dst, const fld_t a, const fld_t b)
 449 {
 450         llimb_t tmp;
 451         llimb_t c[10];
 452
 453         c[0] = (llimb_t)a[0]*b[0];
 454         c[1] = (llimb_t)a[0]*b[1] + (llimb_t)a[1]*b[0];
 455         c[2] = (llimb_t)a[0]*b[2] + (llimb_t)2*a[1]*b[1] + (llimb_t)a[2]*b[0];
 456         c[3] = (llimb_t)a[0]*b[3] + (llimb_t)a[1]*b[2] + (llimb_t)a[2]*b[1]
 457                 + (llimb_t)a[3]*b[0];
 458         c[4] = (llimb_t)a[0]*b[4] + (llimb_t)2*a[1]*b[3] + (llimb_t)a[2]*b[2]
 459                 + (llimb_t)2*a[3]*b[1] + (llimb_t)a[4]*b[0];
 460         c[5] = (llimb_t)a[0]*b[5] + (llimb_t)a[1]*b[4] + (llimb_t)a[2]*b[3]
 461                 + (llimb_t)a[3]*b[2] + (llimb_t)a[4]*b[1] + (llimb_t)a[5]*b[0];
 462         c[6] = (llimb_t)a[0]*b[6] + (llimb_t)2*a[1]*b[5] + (llimb_t)a[2]*b[4]
 463                 + (llimb_t)2*a[3]*b[3] + (llimb_t)a[4]*b[2] + (llimb_t)2*a[5]*b[1]
 464                 + (llimb_t)a[6]*b[0];
 465         c[7] = (llimb_t)a[0]*b[7] + (llimb_t)a[1]*b[6] + (llimb_t)a[2]*b[5]
 466                 + (llimb_t)a[3]*b[4] + (llimb_t)a[4]*b[3] + (llimb_t)a[5]*b[2]
 467                 + (llimb_t)a[6]*b[1] + (llimb_t)a[7]*b[0];
 468         c[8] = (llimb_t)a[0]*b[8] + (llimb_t)2*a[1]*b[7] + (llimb_t)a[2]*b[6]
 469                 + (llimb_t)2*a[3]*b[5] + (llimb_t)a[4]*b[4] + (llimb_t)2*a[5]*b[3]
 470                 + (llimb_t)a[6]*b[2] + (llimb_t)2*a[7]*b[1] + (llimb_t)a[8]*b[0];
 471         c[9] = (llimb_t)a[0]*b[9] + (llimb_t)a[1]*b[8] + (llimb_t)a[2]*b[7]
 472                 + (llimb_t)a[3]*b[6] + (llimb_t)a[4]*b[5] + (llimb_t)a[5]*b[4]
 473                 + (llimb_t)a[6]*b[3] + (llimb_t)a[7]*b[2] + (llimb_t)a[8]*b[1]
 474                 + (llimb_t)a[9]*b[0];
 475
 476         c[0] += 19 * ((llimb_t)2*a[1]*b[9] + (llimb_t)a[2]*b[8] + (llimb_t)2*a[3]*b[7]
 477                       + (llimb_t)a[4]*b[6] + (llimb_t)2*a[5]*b[5] + (llimb_t)a[6]*b[4]
 478                       + (llimb_t)2*a[7]*b[3] + (llimb_t)a[8]*b[2] + (llimb_t)2*a[9]*b[1]);
 479         c[1] += 19 * ((llimb_t)a[2]*b[9] + (llimb_t)a[3]*b[8] + (llimb_t)a[4]*b[7]
 480                       + (llimb_t)a[5]*b[6] + (llimb_t)a[6]*b[5] + (llimb_t)a[7]*b[4]
 481                       + (llimb_t)a[8]*b[3] + (llimb_t)a[9]*b[2]);
 482         c[2] += 19 * ((llimb_t)2*a[3]*b[9] + (llimb_t)a[4]*b[8] + (llimb_t)2*a[5]*b[7]
 483                       + (llimb_t)a[6]*b[6] + (llimb_t)2*a[7]*b[5] + (llimb_t)a[8]*b[4]
 484                       + (llimb_t)2*a[9]*b[3]);
 485         c[3] += 19 * ((llimb_t)a[4]*b[9] + (llimb_t)a[5]*b[8] + (llimb_t)a[6]*b[7]
 486                       + (llimb_t)a[7]*b[6] + (llimb_t)a[8]*b[5] + (llimb_t)a[9]*b[4]);
 487         c[4] += 19 * ((llimb_t)2*a[5]*b[9] + (llimb_t)a[6]*b[8] + (llimb_t)2*a[7]*b[7]
 488                       + (llimb_t)a[8]*b[6] + (llimb_t)2*a[9]*b[5]);
 489         c[5] += 19 * ((llimb_t)a[6]*b[9] + (llimb_t)a[7]*b[8] + (llimb_t)a[8]*b[7]
 490                       + (llimb_t)a[9]*b[6]);
 491         c[6] += 19 * ((llimb_t)2*a[7]*b[9] + (llimb_t)a[8]*b[8] + (llimb_t)2*a[9]*b[7]);
 492         c[7] += 19 * ((llimb_t)a[8]*b[9] + (llimb_t)a[9]*b[8]);
 493         c[8] += 19 * ((llimb_t)2*a[9]*b[9]);
 494
 495         CARRY(c, c, tmp, 0);
 496         CARRY(dst, c, tmp, 0);
 497 }
 498
 499 /*
 500  * fld_sq - square x and reduce modulo q.
 501  */
 502 void
 503 fld_sq(fld_t dst, const fld_t a)
 504 {
 505         llimb_t tmp;
 506         llimb_t c[10];
 507
 508         c[0] = (llimb_t)a[0]*a[0];
 509         c[1] = (llimb_t)2*a[0]*a[1];
 510         c[2] = 2*((llimb_t)a[0]*a[2] + (llimb_t)a[1]*a[1]);
 511         c[3] = 2*((llimb_t)a[0]*a[3] + (llimb_t)a[1]*a[2]);
 512         c[4] = 2*((llimb_t)a[0]*a[4] + (llimb_t)2*a[1]*a[3]) + (llimb_t)a[2]*a[2];
 513         c[5] = 2*((llimb_t)a[0]*a[5] + (llimb_t)a[1]*a[4] + (llimb_t)a[2]*a[3]);
 514         c[6] = 2*((llimb_t)a[0]*a[6] + (llimb_t)2*a[1]*a[5] + (llimb_t)a[2]*a[4] + (llimb_t)a[3]*a[3]);
 515         c[7] = 2*((llimb_t)a[0]*a[7] + (llimb_t)a[1]*a[6] + (llimb_t)a[2]*a[5] + (llimb_t)a[3]*a[4]);
 516         c[8] = 2*((llimb_t)a[0]*a[8] + (llimb_t)2*a[1]*a[7] + (llimb_t)a[2]*a[6] + (llimb_t)2*a[3]*a[5]) + (llimb_t)a[4]*a[4];
 517         c[9] = 2*((llimb_t)a[0]*a[9] + (llimb_t)a[1]*a[8] + (llimb_t)a[2]*a[7] + (llimb_t)a[3]*a[6] + (llimb_t)a[4]*a[5]);
 518
 519         c[0] += 19*2*((llimb_t)2*a[1]*a[9] + (llimb_t)a[2]*a[8] + (llimb_t)2*a[3]*a[7] + (llimb_t)a[4]*a[6] + (llimb_t)a[5]*a[5]);
 520         c[1] += 19*2*((llimb_t)a[2]*a[9] + (llimb_t)a[3]*a[8] + (llimb_t)a[4]*a[7] + (llimb_t)a[5]*a[6]);
 521         c[2] += 19*(2*((llimb_t)2*a[3]*a[9] + (llimb_t)a[4]*a[8] + (llimb_t)2*a[5]*a[7]) + (llimb_t)a[6]*a[6]);
 522         c[3] += 19*2*((llimb_t)a[4]*a[9] + (llimb_t)a[5]*a[8] + (llimb_t)a[6]*a[7]);
 523         c[4] += 19*(2*((llimb_t)2*a[5]*a[9] + (llimb_t)a[6]*a[8] + (llimb_t)a[7]*a[7]));
 524         c[5] += 19*2*((llimb_t)a[6]*a[9] + (llimb_t)a[7]*a[8]);
 525         c[6] += 19*((llimb_t)2*2*a[7]*a[9] + (llimb_t)a[8]*a[8]);
 526         c[7] += 19*2*(llimb_t)a[8]*a[9];
 527         c[8] += 19*2*(llimb_t)a[9]*a[9];
 528
 529         CARRY(c, c, tmp, 0);
 530         CARRY(dst, c, tmp, 0);
 531 }
 532
 533 #endif          /* USE_64BIT */
 534
 535
 536 /*
 537  * common code
 538  */
 539
 540
 541 /*
 542  * fld_eq - compares two field elements in a time-constant manner.
 543  *
 544  * returns 1 if a == b and 0 otherwise.
 545  */
 546 int
 547 fld_eq(const fld_t a, const fld_t b)
 548 {
 549         fld_t tmp;
 550         limb_t res;
 551         int i;
 552
 553         /* tmp <- a - b */
 554         fld_sub(tmp, a, b);
 555         fld_reduce(tmp, tmp);
 556
 557         /* and check tmp for zero */
 558         res = tmp[0];
 559         for (i = 1; i < FLD_LIMB_NUM; i++)
 560                 res |= tmp[i];
 561         for (i = 4*sizeof(limb_t); i > 0; i >>= 1)
 562                 res |= res >> i;
 563
 564         /* now (res & 1) is zero iff tmp is zero */
 565         res = ~res & 1;
 566
 567         return res;
 568 }
 569
 570
 571 /*
 572  * fld_inv - inverts z modulo q.
 573  *
 574  * this code is taken from nacl. it works by taking z to the q-2
 575  * power. by lagrange's theorem (aka 'fermat's little theorem' in this
 576  * special case) this gives us z^-1 modulo q.
 577  */
 578 void
 579 fld_inv(fld_t res, const fld_t z)
 580 {
 581         fld_t z2;
 582         fld_t z9;
 583         fld_t z11;
 584         fld_t z2_5_0;
 585         fld_t z2_10_0;
 586         fld_t z2_20_0;
 587         fld_t z2_50_0;
 588         fld_t z2_100_0;
 589         fld_t t0;
 590         fld_t t1;
 591         int i;
 592
 593         /* 2 */ fld_sq(z2, z);
 594         /* 4 */ fld_sq(t1, z2);
 595         /* 8 */ fld_sq(t0, t1);
 596         /* 9 */ fld_mul(z9,t0, z);
 597         /* 11 */ fld_mul(z11, z9, z2);
 598         /* 22 */ fld_sq(t0, z11);
 599         /* 2^5 - 2^0 = 31 */ fld_mul(z2_5_0, t0, z9);
 600
 601         /* 2^6 - 2^1 */ fld_sq(t0, z2_5_0);
 602         /* 2^7 - 2^2 */ fld_sq(t1, t0);
 603         /* 2^8 - 2^3 */ fld_sq(t0, t1);
 604         /* 2^9 - 2^4 */ fld_sq(t1, t0);
 605         /* 2^10 - 2^5 */ fld_sq(t0, t1);
 606         /* 2^10 - 2^0 */ fld_mul(z2_10_0, t0, z2_5_0);
 607
 608         /* 2^11 - 2^1 */ fld_sq(t0, z2_10_0);
 609         /* 2^12 - 2^2 */ fld_sq(t1, t0);
 610         /* 2^20 - 2^10 */ for (i = 2; i < 10; i += 2) { fld_sq(t0, t1); fld_sq(t1, t0); }
 611         /* 2^20 - 2^0 */ fld_mul(z2_20_0, t1, z2_10_0);
 612
 613         /* 2^21 - 2^1 */ fld_sq(t0, z2_20_0);
 614         /* 2^22 - 2^2 */ fld_sq(t1, t0);
 615         /* 2^40 - 2^20 */ for (i = 2;i < 20;i += 2) { fld_sq(t0, t1); fld_sq(t1, t0); }
 616         /* 2^40 - 2^0 */ fld_mul(t0, t1, z2_20_0);
 617
 618         /* 2^41 - 2^1 */ fld_sq(t1, t0);
 619         /* 2^42 - 2^2 */ fld_sq(t0, t1);
 620         /* 2^50 - 2^10 */ for (i = 2; i < 10; i += 2) { fld_sq(t1, t0); fld_sq(t0, t1); }
 621         /* 2^50 - 2^0 */ fld_mul(z2_50_0, t0, z2_10_0);
 622
 623         /* 2^51 - 2^1 */ fld_sq(t0, z2_50_0);
 624         /* 2^52 - 2^2 */ fld_sq(t1, t0);
 625         /* 2^100 - 2^50 */ for (i = 2; i < 50; i += 2) { fld_sq(t0, t1); fld_sq(t1, t0); }
 626         /* 2^100 - 2^0 */ fld_mul(z2_100_0, t1, z2_50_0);
 627
 628         /* 2^101 - 2^1 */ fld_sq(t1, z2_100_0);
 629         /* 2^102 - 2^2 */ fld_sq(t0, t1);
 630         /* 2^200 - 2^100 */ for (i = 2; i < 100; i += 2) { fld_sq(t1, t0); fld_sq(t0, t1); }
 631         /* 2^200 - 2^0 */ fld_mul(t1, t0, z2_100_0);
 632
 633         /* 2^201 - 2^1 */ fld_sq(t0, t1);
 634         /* 2^202 - 2^2 */ fld_sq(t1, t0);
 635         /* 2^250 - 2^50 */ for (i = 2; i < 50; i += 2) { fld_sq(t0, t1); fld_sq(t1, t0); }
 636         /* 2^250 - 2^0 */ fld_mul(t0, t1, z2_50_0);
 637
 638         /* 2^251 - 2^1 */ fld_sq(t1, t0);
 639         /* 2^252 - 2^2 */ fld_sq(t0, t1);
 640
 641         /* 2^253 - 2^3 */ fld_sq(t1, t0);
 642         /* 2^254 - 2^4 */ fld_sq(t0, t1);
 643         /* 2^255 - 2^5 */ fld_sq(t1, t0);
 644         /* 2^255 - 21 */ fld_mul(res, t1, z11);
 645 }
 646
 647
 648 /*
 649  * fld_pow2523 - compute z^((q-5)/8) modulo q, ie (z*res)^2 is either z
 650  * or -z modulo q.
 651  *
 652  * this function is used to mix a square-root modulo q with an invertation in
 653  * ed_import. see the ed25519 paper for an explanation.
 654  *
 655  * this code is, like fld_inv, taken from nacl.
 656  */
 657 void
 658 fld_pow2523(fld_t res, const fld_t z)
 659 {
 660         fld_t z2;
 661         fld_t z9;
 662         fld_t z2_5_0;
 663         fld_t z2_10_0;
 664         fld_t z2_20_0;
 665         fld_t z2_50_0;
 666         fld_t z2_100_0;
 667         fld_t t;
 668         int i;
 669
 670         /* 2 */ fld_sq(z2, z);
 671         /* 4 */ fld_sq(t, z2);
 672         /* 8 */ fld_sq(t, t);
 673         /* 9 */ fld_mul(z9, t, z);
 674         /* 11 */ fld_mul(t, z9, z2);
 675         /* 22 */ fld_sq(t, t);
 676         /* 2^5 - 2^0 = 31 */ fld_mul(z2_5_0, t, z9);
 677
 678         /* 2^6 - 2^1 */ fld_sq(t, z2_5_0);
 679         /* 2^10 - 2^5 */ for (i = 1;i < 5;i++) { fld_sq(t, t); }
 680         /* 2^10 - 2^0 */ fld_mul(z2_10_0, t, z2_5_0);
 681
 682         /* 2^11 - 2^1 */ fld_sq(t, z2_10_0);
 683         /* 2^20 - 2^10 */ for (i = 1;i < 10;i++) { fld_sq(t, t); }
 684         /* 2^20 - 2^0 */ fld_mul(z2_20_0, t, z2_10_0);
 685
 686         /* 2^21 - 2^1 */ fld_sq(t, z2_20_0);
 687         /* 2^40 - 2^20 */ for (i = 1;i < 20;i++) { fld_sq(t, t); }
 688         /* 2^40 - 2^0 */ fld_mul(t, t, z2_20_0);
 689
 690         /* 2^41 - 2^1 */ fld_sq(t, t);
 691         /* 2^50 - 2^10 */ for (i = 1;i < 10;i++) { fld_sq(t, t); }
 692         /* 2^50 - 2^0 */ fld_mul(z2_50_0, t, z2_10_0);
 693
 694         /* 2^51 - 2^1 */ fld_sq(t, z2_50_0);
 695         /* 2^100 - 2^50 */ for (i = 1;i < 50;i++) { fld_sq(t, t); }
 696         /* 2^100 - 2^0 */ fld_mul(z2_100_0, t, z2_50_0);
 697
 698         /* 2^101 - 2^1 */ fld_sq(t, z2_100_0);
 699         /* 2^200 - 2^100 */ for (i = 1;i < 100;i++) { fld_sq(t, t); }
 700         /* 2^200 - 2^0 */ fld_mul(t, t, z2_100_0);
 701
 702         /* 2^201 - 2^1 */ fld_sq(t, t);
 703         /* 2^250 - 2^50 */ for (i = 1;i < 50;i++) { fld_sq(t, t); }
 704         /* 2^250 - 2^0 */ fld_mul(t, t, z2_50_0);
 705
 706         /* 2^251 - 2^1 */ fld_sq(t, t);
 707         /* 2^252 - 2^2 */ fld_sq(t, t);
 708         /* 2^252 - 3 */ fld_mul(res, t, z);
 709 }