crypto/ecc.c

   1 /*
   2  * Copyright (c) 2013, Kenneth MacKay
   3  * All rights reserved.
   4  *
   5  * Redistribution and use in source and binary forms, with or without
   6  * modification, are permitted provided that the following conditions are
   7  * met:
   8  *  * Redistributions of source code must retain the above copyright
   9  *   notice, this list of conditions and the following disclaimer.
  10  *  * Redistributions in binary form must reproduce the above copyright
  11  *    notice, this list of conditions and the following disclaimer in the
  12  *    documentation and/or other materials provided with the distribution.
  13  *
  14  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  15  * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  16  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  17  * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
  18  * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
  19  * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
  20  * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  21  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  22  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  23  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  24  * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  25  */
  26
  27 #include <linux/random.h>
  28 #include <linux/slab.h>
  29 #include <linux/swab.h>
  30 #include <linux/fips.h>
  31 #include <crypto/ecdh.h>
  32 #include <crypto/rng.h>
  33
  34 #include "ecc.h"
  35 #include "ecc_curve_defs.h"
  36
  37 typedef struct {
  38         u64 m_low;
  39         u64 m_high;
  40 } uint128_t;
  41
  42 static inline const struct ecc_curve *ecc_get_curve(unsigned int curve_id)
  43 {
  44         switch (curve_id) {
  45         /* In FIPS mode only allow P256 and higher */
  46         case ECC_CURVE_NIST_P192:
  47                 return fips_enabled ? NULL : &nist_p192;
  48         case ECC_CURVE_NIST_P256:
  49                 return &nist_p256;
  50         default:
  51                 return NULL;
  52         }
  53 }
  54
  55 static u64 *ecc_alloc_digits_space(unsigned int ndigits)
  56 {
  57         size_t len = ndigits * sizeof(u64);
  58
  59         if (!len)
  60                 return NULL;
  61
  62         return kmalloc(len, GFP_KERNEL);
  63 }
  64
  65 static void ecc_free_digits_space(u64 *space)
  66 {
  67         kzfree(space);
  68 }
  69
  70 static struct ecc_point *ecc_alloc_point(unsigned int ndigits)
  71 {
  72         struct ecc_point *p = kmalloc(sizeof(*p), GFP_KERNEL);
  73
  74         if (!p)
  75                 return NULL;
  76
  77         p->x = ecc_alloc_digits_space(ndigits);
  78         if (!p->x)
  79                 goto err_alloc_x;
  80
  81         p->y = ecc_alloc_digits_space(ndigits);
  82         if (!p->y)
  83                 goto err_alloc_y;
  84
  85         p->ndigits = ndigits;
  86
  87         return p;
  88
  89 err_alloc_y:
  90         ecc_free_digits_space(p->x);
  91 err_alloc_x:
  92         kfree(p);
  93         return NULL;
  94 }
  95
  96 static void ecc_free_point(struct ecc_point *p)
  97 {
  98         if (!p)
  99                 return;
 100
 101         kzfree(p->x);
 102         kzfree(p->y);
 103         kzfree(p);
 104 }
 105
 106 static void vli_clear(u64 *vli, unsigned int ndigits)
 107 {
 108         int i;
 109
 110         for (i = 0; i < ndigits; i++)
 111                 vli[i] = 0;
 112 }
 113
 114 /* Returns true if vli == 0, false otherwise. */
 115 static bool vli_is_zero(const u64 *vli, unsigned int ndigits)
 116 {
 117         int i;
 118
 119         for (i = 0; i < ndigits; i++) {
 120                 if (vli[i])
 121                         return false;
 122         }
 123
 124         return true;
 125 }
 126
 127 /* Returns nonzero if bit bit of vli is set. */
 128 static u64 vli_test_bit(const u64 *vli, unsigned int bit)
 129 {
 130         return (vli[bit / 64] & ((u64)1 << (bit % 64)));
 131 }
 132
 133 /* Counts the number of 64-bit "digits" in vli. */
 134 static unsigned int vli_num_digits(const u64 *vli, unsigned int ndigits)
 135 {
 136         int i;
 137
 138         /* Search from the end until we find a non-zero digit.
 139          * We do it in reverse because we expect that most digits will
 140          * be nonzero.
 141          */
 142         for (i = ndigits - 1; i >= 0 && vli[i] == 0; i--);
 143
 144         return (i + 1);
 145 }
 146
 147 /* Counts the number of bits required for vli. */
 148 static unsigned int vli_num_bits(const u64 *vli, unsigned int ndigits)
 149 {
 150         unsigned int i, num_digits;
 151         u64 digit;
 152
 153         num_digits = vli_num_digits(vli, ndigits);
 154         if (num_digits == 0)
 155                 return 0;
 156
 157         digit = vli[num_digits - 1];
 158         for (i = 0; digit; i++)
 159                 digit >>= 1;
 160
 161         return ((num_digits - 1) * 64 + i);
 162 }
 163
 164 /* Sets dest = src. */
 165 static void vli_set(u64 *dest, const u64 *src, unsigned int ndigits)
 166 {
 167         int i;
 168
 169         for (i = 0; i < ndigits; i++)
 170                 dest[i] = src[i];
 171 }
 172
 173 /* Returns sign of left - right. */
 174 static int vli_cmp(const u64 *left, const u64 *right, unsigned int ndigits)
 175 {
 176         int i;
 177
 178         for (i = ndigits - 1; i >= 0; i--) {
 179                 if (left[i] > right[i])
 180                         return 1;
 181                 else if (left[i] < right[i])
 182                         return -1;
 183         }
 184
 185         return 0;
 186 }
 187
 188 /* Computes result = in << c, returning carry. Can modify in place
 189  * (if result == in). 0 < shift < 64.
 190  */
 191 static u64 vli_lshift(u64 *result, const u64 *in, unsigned int shift,
 192                       unsigned int ndigits)
 193 {
 194         u64 carry = 0;
 195         int i;
 196
 197         for (i = 0; i < ndigits; i++) {
 198                 u64 temp = in[i];
 199
 200                 result[i] = (temp << shift) | carry;
 201                 carry = temp >> (64 - shift);
 202         }
 203
 204         return carry;
 205 }
 206
 207 /* Computes vli = vli >> 1. */
 208 static void vli_rshift1(u64 *vli, unsigned int ndigits)
 209 {
 210         u64 *end = vli;
 211         u64 carry = 0;
 212
 213         vli += ndigits;
 214
 215         while (vli-- > end) {
 216                 u64 temp = *vli;
 217                 *vli = (temp >> 1) | carry;
 218                 carry = temp << 63;
 219         }
 220 }
 221
 222 /* Computes result = left + right, returning carry. Can modify in place. */
 223 static u64 vli_add(u64 *result, const u64 *left, const u64 *right,
 224                    unsigned int ndigits)
 225 {
 226         u64 carry = 0;
 227         int i;
 228
 229         for (i = 0; i < ndigits; i++) {
 230                 u64 sum;
 231
 232                 sum = left[i] + right[i] + carry;
 233                 if (sum != left[i])
 234                         carry = (sum < left[i]);
 235
 236                 result[i] = sum;
 237         }
 238
 239         return carry;
 240 }
 241
 242 /* Computes result = left - right, returning borrow. Can modify in place. */
 243 static u64 vli_sub(u64 *result, const u64 *left, const u64 *right,
 244                    unsigned int ndigits)
 245 {
 246         u64 borrow = 0;
 247         int i;
 248
 249         for (i = 0; i < ndigits; i++) {
 250                 u64 diff;
 251
 252                 diff = left[i] - right[i] - borrow;
 253                 if (diff != left[i])
 254                         borrow = (diff > left[i]);
 255
 256                 result[i] = diff;
 257         }
 258
 259         return borrow;
 260 }
 261
 262 static uint128_t mul_64_64(u64 left, u64 right)
 263 {
 264         u64 a0 = left & 0xffffffffull;
 265         u64 a1 = left >> 32;
 266         u64 b0 = right & 0xffffffffull;
 267         u64 b1 = right >> 32;
 268         u64 m0 = a0 * b0;
 269         u64 m1 = a0 * b1;
 270         u64 m2 = a1 * b0;
 271         u64 m3 = a1 * b1;
 272         uint128_t result;
 273
 274         m2 += (m0 >> 32);
 275         m2 += m1;
 276
 277         /* Overflow */
 278         if (m2 < m1)
 279                 m3 += 0x100000000ull;
 280
 281         result.m_low = (m0 & 0xffffffffull) | (m2 << 32);
 282         result.m_high = m3 + (m2 >> 32);
 283
 284         return result;
 285 }
 286
 287 static uint128_t add_128_128(uint128_t a, uint128_t b)
 288 {
 289         uint128_t result;
 290
 291         result.m_low = a.m_low + b.m_low;
 292         result.m_high = a.m_high + b.m_high + (result.m_low < a.m_low);
 293
 294         return result;
 295 }
 296
 297 static void vli_mult(u64 *result, const u64 *left, const u64 *right,
 298                      unsigned int ndigits)
 299 {
 300         uint128_t r01 = { 0, 0 };
 301         u64 r2 = 0;
 302         unsigned int i, k;
 303
 304         /* Compute each digit of result in sequence, maintaining the
 305          * carries.
 306          */
 307         for (k = 0; k < ndigits * 2 - 1; k++) {
 308                 unsigned int min;
 309
 310                 if (k < ndigits)
 311                         min = 0;
 312                 else
 313                         min = (k + 1) - ndigits;
 314
 315                 for (i = min; i <= k && i < ndigits; i++) {
 316                         uint128_t product;
 317
 318                         product = mul_64_64(left[i], right[k - i]);
 319
 320                         r01 = add_128_128(r01, product);
 321                         r2 += (r01.m_high < product.m_high);
 322                 }
 323
 324                 result[k] = r01.m_low;
 325                 r01.m_low = r01.m_high;
 326                 r01.m_high = r2;
 327                 r2 = 0;
 328         }
 329
 330         result[ndigits * 2 - 1] = r01.m_low;
 331 }
 332
 333 static void vli_square(u64 *result, const u64 *left, unsigned int ndigits)
 334 {
 335         uint128_t r01 = { 0, 0 };
 336         u64 r2 = 0;
 337         int i, k;
 338
 339         for (k = 0; k < ndigits * 2 - 1; k++) {
 340                 unsigned int min;
 341
 342                 if (k < ndigits)
 343                         min = 0;
 344                 else
 345                         min = (k + 1) - ndigits;
 346
 347                 for (i = min; i <= k && i <= k - i; i++) {
 348                         uint128_t product;
 349
 350                         product = mul_64_64(left[i], left[k - i]);
 351
 352                         if (i < k - i) {
 353                                 r2 += product.m_high >> 63;
 354                                 product.m_high = (product.m_high << 1) |
 355                                                  (product.m_low >> 63);
 356                                 product.m_low <<= 1;
 357                         }
 358
 359                         r01 = add_128_128(r01, product);
 360                         r2 += (r01.m_high < product.m_high);
 361                 }
 362
 363                 result[k] = r01.m_low;
 364                 r01.m_low = r01.m_high;
 365                 r01.m_high = r2;
 366                 r2 = 0;
 367         }
 368
 369         result[ndigits * 2 - 1] = r01.m_low;
 370 }
 371
 372 /* Computes result = (left + right) % mod.
 373  * Assumes that left < mod and right < mod, result != mod.
 374  */
 375 static void vli_mod_add(u64 *result, const u64 *left, const u64 *right,
 376                         const u64 *mod, unsigned int ndigits)
 377 {
 378         u64 carry;
 379
 380         carry = vli_add(result, left, right, ndigits);
 381
 382         /* result > mod (result = mod + remainder), so subtract mod to
 383          * get remainder.
 384          */
 385         if (carry || vli_cmp(result, mod, ndigits) >= 0)
 386                 vli_sub(result, result, mod, ndigits);
 387 }
 388
 389 /* Computes result = (left - right) % mod.
 390  * Assumes that left < mod and right < mod, result != mod.
 391  */
 392 static void vli_mod_sub(u64 *result, const u64 *left, const u64 *right,
 393                         const u64 *mod, unsigned int ndigits)
 394 {
 395         u64 borrow = vli_sub(result, left, right, ndigits);
 396
 397         /* In this case, p_result == -diff == (max int) - diff.
 398          * Since -x % d == d - x, we can get the correct result from
 399          * result + mod (with overflow).
 400          */
 401         if (borrow)
 402                 vli_add(result, result, mod, ndigits);
 403 }
 404
 405 /* Computes p_result = p_product % curve_p.
 406  * See algorithm 5 and 6 from
 407  * http://www.isys.uni-klu.ac.at/PDF/2001-0126-MT.pdf
 408  */
 409 static void vli_mmod_fast_192(u64 *result, const u64 *product,
 410                               const u64 *curve_prime, u64 *tmp)
 411 {
 412         const unsigned int ndigits = 3;
 413         int carry;
 414
 415         vli_set(result, product, ndigits);
 416
 417         vli_set(tmp, &product[3], ndigits);
 418         carry = vli_add(result, result, tmp, ndigits);
 419
 420         tmp[0] = 0;
 421         tmp[1] = product[3];
 422         tmp[2] = product[4];
 423         carry += vli_add(result, result, tmp, ndigits);
 424
 425         tmp[0] = tmp[1] = product[5];
 426         tmp[2] = 0;
 427         carry += vli_add(result, result, tmp, ndigits);
 428
 429         while (carry || vli_cmp(curve_prime, result, ndigits) != 1)
 430                 carry -= vli_sub(result, result, curve_prime, ndigits);
 431 }
 432
 433 /* Computes result = product % curve_prime
 434  * from http://www.nsa.gov/ia/_files/nist-routines.pdf
 435  */
 436 static void vli_mmod_fast_256(u64 *result, const u64 *product,
 437                               const u64 *curve_prime, u64 *tmp)
 438 {
 439         int carry;
 440         const unsigned int ndigits = 4;
 441
 442         /* t */
 443         vli_set(result, product, ndigits);
 444
 445         /* s1 */
 446         tmp[0] = 0;
 447         tmp[1] = product[5] & 0xffffffff00000000ull;
 448         tmp[2] = product[6];
 449         tmp[3] = product[7];
 450         carry = vli_lshift(tmp, tmp, 1, ndigits);
 451         carry += vli_add(result, result, tmp, ndigits);
 452
 453         /* s2 */
 454         tmp[1] = product[6] << 32;
 455         tmp[2] = (product[6] >> 32) | (product[7] << 32);
 456         tmp[3] = product[7] >> 32;
 457         carry += vli_lshift(tmp, tmp, 1, ndigits);
 458         carry += vli_add(result, result, tmp, ndigits);
 459
 460         /* s3 */
 461         tmp[0] = product[4];
 462         tmp[1] = product[5] & 0xffffffff;
 463         tmp[2] = 0;
 464         tmp[3] = product[7];
 465         carry += vli_add(result, result, tmp, ndigits);
 466
 467         /* s4 */
 468         tmp[0] = (product[4] >> 32) | (product[5] << 32);
 469         tmp[1] = (product[5] >> 32) | (product[6] & 0xffffffff00000000ull);
 470         tmp[2] = product[7];
 471         tmp[3] = (product[6] >> 32) | (product[4] << 32);
 472         carry += vli_add(result, result, tmp, ndigits);
 473
 474         /* d1 */
 475         tmp[0] = (product[5] >> 32) | (product[6] << 32);
 476         tmp[1] = (product[6] >> 32);
 477         tmp[2] = 0;
 478         tmp[3] = (product[4] & 0xffffffff) | (product[5] << 32);
 479         carry -= vli_sub(result, result, tmp, ndigits);
 480
 481         /* d2 */
 482         tmp[0] = product[6];
 483         tmp[1] = product[7];
 484         tmp[2] = 0;
 485         tmp[3] = (product[4] >> 32) | (product[5] & 0xffffffff00000000ull);
 486         carry -= vli_sub(result, result, tmp, ndigits);
 487
 488         /* d3 */
 489         tmp[0] = (product[6] >> 32) | (product[7] << 32);
 490         tmp[1] = (product[7] >> 32) | (product[4] << 32);
 491         tmp[2] = (product[4] >> 32) | (product[5] << 32);
 492         tmp[3] = (product[6] << 32);
 493         carry -= vli_sub(result, result, tmp, ndigits);
 494
 495         /* d4 */
 496         tmp[0] = product[7];
 497         tmp[1] = product[4] & 0xffffffff00000000ull;
 498         tmp[2] = product[5];
 499         tmp[3] = product[6] & 0xffffffff00000000ull;
 500         carry -= vli_sub(result, result, tmp, ndigits);
 501
 502         if (carry < 0) {
 503                 do {
 504                         carry += vli_add(result, result, curve_prime, ndigits);
 505                 } while (carry < 0);
 506         } else {
 507                 while (carry || vli_cmp(curve_prime, result, ndigits) != 1)
 508                         carry -= vli_sub(result, result, curve_prime, ndigits);
 509         }
 510 }
 511
 512 /* Computes result = product % curve_prime
 513  *  from http://www.nsa.gov/ia/_files/nist-routines.pdf
 514 */
 515 static bool vli_mmod_fast(u64 *result, u64 *product,
 516                           const u64 *curve_prime, unsigned int ndigits)
 517 {
 518         u64 tmp[2 * ECC_MAX_DIGITS];
 519
 520         switch (ndigits) {
 521         case 3:
 522                 vli_mmod_fast_192(result, product, curve_prime, tmp);
 523                 break;
 524         case 4:
 525                 vli_mmod_fast_256(result, product, curve_prime, tmp);
 526                 break;
 527         default:
 528                 pr_err("unsupports digits size!\n");
 529                 return false;
 530         }
 531
 532         return true;
 533 }
 534
 535 /* Computes result = (left * right) % curve_prime. */
 536 static void vli_mod_mult_fast(u64 *result, const u64 *left, const u64 *right,
 537                               const u64 *curve_prime, unsigned int ndigits)
 538 {
 539         u64 product[2 * ECC_MAX_DIGITS];
 540
 541         vli_mult(product, left, right, ndigits);
 542         vli_mmod_fast(result, product, curve_prime, ndigits);
 543 }
 544
 545 /* Computes result = left^2 % curve_prime. */
 546 static void vli_mod_square_fast(u64 *result, const u64 *left,
 547                                 const u64 *curve_prime, unsigned int ndigits)
 548 {
 549         u64 product[2 * ECC_MAX_DIGITS];
 550
 551         vli_square(product, left, ndigits);
 552         vli_mmod_fast(result, product, curve_prime, ndigits);
 553 }
 554
 555 #define EVEN(vli) (!(vli[0] & 1))
 556 /* Computes result = (1 / p_input) % mod. All VLIs are the same size.
 557  * See "From Euclid's GCD to Montgomery Multiplication to the Great Divide"
 558  * https://labs.oracle.com/techrep/2001/smli_tr-2001-95.pdf
 559  */
 560 static void vli_mod_inv(u64 *result, const u64 *input, const u64 *mod,
 561                         unsigned int ndigits)
 562 {
 563         u64 a[ECC_MAX_DIGITS], b[ECC_MAX_DIGITS];
 564         u64 u[ECC_MAX_DIGITS], v[ECC_MAX_DIGITS];
 565         u64 carry;
 566         int cmp_result;
 567
 568         if (vli_is_zero(input, ndigits)) {
 569                 vli_clear(result, ndigits);
 570                 return;
 571         }
 572
 573         vli_set(a, input, ndigits);
 574         vli_set(b, mod, ndigits);
 575         vli_clear(u, ndigits);
 576         u[0] = 1;
 577         vli_clear(v, ndigits);
 578
 579         while ((cmp_result = vli_cmp(a, b, ndigits)) != 0) {
 580                 carry = 0;
 581
 582                 if (EVEN(a)) {
 583                         vli_rshift1(a, ndigits);
 584
 585                         if (!EVEN(u))
 586                                 carry = vli_add(u, u, mod, ndigits);
 587
 588                         vli_rshift1(u, ndigits);
 589                         if (carry)
 590                                 u[ndigits - 1] |= 0x8000000000000000ull;
 591                 } else if (EVEN(b)) {
 592                         vli_rshift1(b, ndigits);
 593
 594                         if (!EVEN(v))
 595                                 carry = vli_add(v, v, mod, ndigits);
 596
 597                         vli_rshift1(v, ndigits);
 598                         if (carry)
 599                                 v[ndigits - 1] |= 0x8000000000000000ull;
 600                 } else if (cmp_result > 0) {
 601                         vli_sub(a, a, b, ndigits);
 602                         vli_rshift1(a, ndigits);
 603
 604                         if (vli_cmp(u, v, ndigits) < 0)
 605                                 vli_add(u, u, mod, ndigits);
 606
 607                         vli_sub(u, u, v, ndigits);
 608                         if (!EVEN(u))
 609                                 carry = vli_add(u, u, mod, ndigits);
 610
 611                         vli_rshift1(u, ndigits);
 612                         if (carry)
 613                                 u[ndigits - 1] |= 0x8000000000000000ull;
 614                 } else {
 615                         vli_sub(b, b, a, ndigits);
 616                         vli_rshift1(b, ndigits);
 617
 618                         if (vli_cmp(v, u, ndigits) < 0)
 619                                 vli_add(v, v, mod, ndigits);
 620
 621                         vli_sub(v, v, u, ndigits);
 622                         if (!EVEN(v))
 623                                 carry = vli_add(v, v, mod, ndigits);
 624
 625                         vli_rshift1(v, ndigits);
 626                         if (carry)
 627                                 v[ndigits - 1] |= 0x8000000000000000ull;
 628                 }
 629         }
 630
 631         vli_set(result, u, ndigits);
 632 }
 633
 634 /* ------ Point operations ------ */
 635
 636 /* Returns true if p_point is the point at infinity, false otherwise. */
 637 static bool ecc_point_is_zero(const struct ecc_point *point)
 638 {
 639         return (vli_is_zero(point->x, point->ndigits) &&
 640                 vli_is_zero(point->y, point->ndigits));
 641 }
 642
 643 /* Point multiplication algorithm using Montgomery's ladder with co-Z
 644  * coordinates. From http://eprint.iacr.org/2011/338.pdf
 645  */
 646
 647 /* Double in place */
 648 static void ecc_point_double_jacobian(u64 *x1, u64 *y1, u64 *z1,
 649                                       u64 *curve_prime, unsigned int ndigits)
 650 {
 651         /* t1 = x, t2 = y, t3 = z */
 652         u64 t4[ECC_MAX_DIGITS];
 653         u64 t5[ECC_MAX_DIGITS];
 654
 655         if (vli_is_zero(z1, ndigits))
 656                 return;
 657
 658         /* t4 = y1^2 */
 659         vli_mod_square_fast(t4, y1, curve_prime, ndigits);
 660         /* t5 = x1*y1^2 = A */
 661         vli_mod_mult_fast(t5, x1, t4, curve_prime, ndigits);
 662         /* t4 = y1^4 */
 663         vli_mod_square_fast(t4, t4, curve_prime, ndigits);
 664         /* t2 = y1*z1 = z3 */
 665         vli_mod_mult_fast(y1, y1, z1, curve_prime, ndigits);
 666         /* t3 = z1^2 */
 667         vli_mod_square_fast(z1, z1, curve_prime, ndigits);
 668
 669         /* t1 = x1 + z1^2 */
 670         vli_mod_add(x1, x1, z1, curve_prime, ndigits);
 671         /* t3 = 2*z1^2 */
 672         vli_mod_add(z1, z1, z1, curve_prime, ndigits);
 673         /* t3 = x1 - z1^2 */
 674         vli_mod_sub(z1, x1, z1, curve_prime, ndigits);
 675         /* t1 = x1^2 - z1^4 */
 676         vli_mod_mult_fast(x1, x1, z1, curve_prime, ndigits);
 677
 678         /* t3 = 2*(x1^2 - z1^4) */
 679         vli_mod_add(z1, x1, x1, curve_prime, ndigits);
 680         /* t1 = 3*(x1^2 - z1^4) */
 681         vli_mod_add(x1, x1, z1, curve_prime, ndigits);
 682         if (vli_test_bit(x1, 0)) {
 683                 u64 carry = vli_add(x1, x1, curve_prime, ndigits);
 684
 685                 vli_rshift1(x1, ndigits);
 686                 x1[ndigits - 1] |= carry << 63;
 687         } else {
 688                 vli_rshift1(x1, ndigits);
 689         }
 690         /* t1 = 3/2*(x1^2 - z1^4) = B */
 691
 692         /* t3 = B^2 */
 693         vli_mod_square_fast(z1, x1, curve_prime, ndigits);
 694         /* t3 = B^2 - A */
 695         vli_mod_sub(z1, z1, t5, curve_prime, ndigits);
 696         /* t3 = B^2 - 2A = x3 */
 697         vli_mod_sub(z1, z1, t5, curve_prime, ndigits);
 698         /* t5 = A - x3 */
 699         vli_mod_sub(t5, t5, z1, curve_prime, ndigits);
 700         /* t1 = B * (A - x3) */
 701         vli_mod_mult_fast(x1, x1, t5, curve_prime, ndigits);
 702         /* t4 = B * (A - x3) - y1^4 = y3 */
 703         vli_mod_sub(t4, x1, t4, curve_prime, ndigits);
 704
 705         vli_set(x1, z1, ndigits);
 706         vli_set(z1, y1, ndigits);
 707         vli_set(y1, t4, ndigits);
 708 }
 709
 710 /* Modify (x1, y1) => (x1 * z^2, y1 * z^3) */
 711 static void apply_z(u64 *x1, u64 *y1, u64 *z, u64 *curve_prime,
 712                     unsigned int ndigits)
 713 {
 714         u64 t1[ECC_MAX_DIGITS];
 715
 716         vli_mod_square_fast(t1, z, curve_prime, ndigits);    /* z^2 */
 717         vli_mod_mult_fast(x1, x1, t1, curve_prime, ndigits); /* x1 * z^2 */
 718         vli_mod_mult_fast(t1, t1, z, curve_prime, ndigits);  /* z^3 */
 719         vli_mod_mult_fast(y1, y1, t1, curve_prime, ndigits); /* y1 * z^3 */
 720 }
 721
 722 /* P = (x1, y1) => 2P, (x2, y2) => P' */
 723 static void xycz_initial_double(u64 *x1, u64 *y1, u64 *x2, u64 *y2,
 724                                 u64 *p_initial_z, u64 *curve_prime,
 725                                 unsigned int ndigits)
 726 {
 727         u64 z[ECC_MAX_DIGITS];
 728
 729         vli_set(x2, x1, ndigits);
 730         vli_set(y2, y1, ndigits);
 731
 732         vli_clear(z, ndigits);
 733         z[0] = 1;
 734
 735         if (p_initial_z)
 736                 vli_set(z, p_initial_z, ndigits);
 737
 738         apply_z(x1, y1, z, curve_prime, ndigits);
 739
 740         ecc_point_double_jacobian(x1, y1, z, curve_prime, ndigits);
 741
 742         apply_z(x2, y2, z, curve_prime, ndigits);
 743 }
 744
 745 /* Input P = (x1, y1, Z), Q = (x2, y2, Z)
 746  * Output P' = (x1', y1', Z3), P + Q = (x3, y3, Z3)
 747  * or P => P', Q => P + Q
 748  */
 749 static void xycz_add(u64 *x1, u64 *y1, u64 *x2, u64 *y2, u64 *curve_prime,
 750                      unsigned int ndigits)
 751 {
 752         /* t1 = X1, t2 = Y1, t3 = X2, t4 = Y2 */
 753         u64 t5[ECC_MAX_DIGITS];
 754
 755         /* t5 = x2 - x1 */
 756         vli_mod_sub(t5, x2, x1, curve_prime, ndigits);
 757         /* t5 = (x2 - x1)^2 = A */
 758         vli_mod_square_fast(t5, t5, curve_prime, ndigits);
 759         /* t1 = x1*A = B */
 760         vli_mod_mult_fast(x1, x1, t5, curve_prime, ndigits);
 761         /* t3 = x2*A = C */
 762         vli_mod_mult_fast(x2, x2, t5, curve_prime, ndigits);
 763         /* t4 = y2 - y1 */
 764         vli_mod_sub(y2, y2, y1, curve_prime, ndigits);
 765         /* t5 = (y2 - y1)^2 = D */
 766         vli_mod_square_fast(t5, y2, curve_prime, ndigits);
 767
 768         /* t5 = D - B */
 769         vli_mod_sub(t5, t5, x1, curve_prime, ndigits);
 770         /* t5 = D - B - C = x3 */
 771         vli_mod_sub(t5, t5, x2, curve_prime, ndigits);
 772         /* t3 = C - B */
 773         vli_mod_sub(x2, x2, x1, curve_prime, ndigits);
 774         /* t2 = y1*(C - B) */
 775         vli_mod_mult_fast(y1, y1, x2, curve_prime, ndigits);
 776         /* t3 = B - x3 */
 777         vli_mod_sub(x2, x1, t5, curve_prime, ndigits);
 778         /* t4 = (y2 - y1)*(B - x3) */
 779         vli_mod_mult_fast(y2, y2, x2, curve_prime, ndigits);
 780         /* t4 = y3 */
 781         vli_mod_sub(y2, y2, y1, curve_prime, ndigits);
 782
 783         vli_set(x2, t5, ndigits);
 784 }
 785
 786 /* Input P = (x1, y1, Z), Q = (x2, y2, Z)
 787  * Output P + Q = (x3, y3, Z3), P - Q = (x3', y3', Z3)
 788  * or P => P - Q, Q => P + Q
 789  */
 790 static void xycz_add_c(u64 *x1, u64 *y1, u64 *x2, u64 *y2, u64 *curve_prime,
 791                        unsigned int ndigits)
 792 {
 793         /* t1 = X1, t2 = Y1, t3 = X2, t4 = Y2 */
 794         u64 t5[ECC_MAX_DIGITS];
 795         u64 t6[ECC_MAX_DIGITS];
 796         u64 t7[ECC_MAX_DIGITS];
 797
 798         /* t5 = x2 - x1 */
 799         vli_mod_sub(t5, x2, x1, curve_prime, ndigits);
 800         /* t5 = (x2 - x1)^2 = A */
 801         vli_mod_square_fast(t5, t5, curve_prime, ndigits);
 802         /* t1 = x1*A = B */
 803         vli_mod_mult_fast(x1, x1, t5, curve_prime, ndigits);
 804         /* t3 = x2*A = C */
 805         vli_mod_mult_fast(x2, x2, t5, curve_prime, ndigits);
 806         /* t4 = y2 + y1 */
 807         vli_mod_add(t5, y2, y1, curve_prime, ndigits);
 808         /* t4 = y2 - y1 */
 809         vli_mod_sub(y2, y2, y1, curve_prime, ndigits);
 810
 811         /* t6 = C - B */
 812         vli_mod_sub(t6, x2, x1, curve_prime, ndigits);
 813         /* t2 = y1 * (C - B) */
 814         vli_mod_mult_fast(y1, y1, t6, curve_prime, ndigits);
 815         /* t6 = B + C */
 816         vli_mod_add(t6, x1, x2, curve_prime, ndigits);
 817         /* t3 = (y2 - y1)^2 */
 818         vli_mod_square_fast(x2, y2, curve_prime, ndigits);
 819         /* t3 = x3 */
 820         vli_mod_sub(x2, x2, t6, curve_prime, ndigits);
 821
 822         /* t7 = B - x3 */
 823         vli_mod_sub(t7, x1, x2, curve_prime, ndigits);
 824         /* t4 = (y2 - y1)*(B - x3) */
 825         vli_mod_mult_fast(y2, y2, t7, curve_prime, ndigits);
 826         /* t4 = y3 */
 827         vli_mod_sub(y2, y2, y1, curve_prime, ndigits);
 828
 829         /* t7 = (y2 + y1)^2 = F */
 830         vli_mod_square_fast(t7, t5, curve_prime, ndigits);
 831         /* t7 = x3' */
 832         vli_mod_sub(t7, t7, t6, curve_prime, ndigits);
 833         /* t6 = x3' - B */
 834         vli_mod_sub(t6, t7, x1, curve_prime, ndigits);
 835         /* t6 = (y2 + y1)*(x3' - B) */
 836         vli_mod_mult_fast(t6, t6, t5, curve_prime, ndigits);
 837         /* t2 = y3' */
 838         vli_mod_sub(y1, t6, y1, curve_prime, ndigits);
 839
 840         vli_set(x1, t7, ndigits);
 841 }
 842
 843 static void ecc_point_mult(struct ecc_point *result,
 844                            const struct ecc_point *point, const u64 *scalar,
 845                            u64 *initial_z, const struct ecc_curve *curve,
 846                            unsigned int ndigits)
 847 {
 848         /* R0 and R1 */
 849         u64 rx[2][ECC_MAX_DIGITS];
 850         u64 ry[2][ECC_MAX_DIGITS];
 851         u64 z[ECC_MAX_DIGITS];
 852         u64 sk[2][ECC_MAX_DIGITS];
 853         u64 *curve_prime = curve->p;
 854         int i, nb;
 855         int num_bits;
 856         int carry;
 857
 858         carry = vli_add(sk[0], scalar, curve->n, ndigits);
 859         vli_add(sk[1], sk[0], curve->n, ndigits);
 860         scalar = sk[!carry];
 861         num_bits = sizeof(u64) * ndigits * 8 + 1;
 862
 863         vli_set(rx[1], point->x, ndigits);
 864         vli_set(ry[1], point->y, ndigits);
 865
 866         xycz_initial_double(rx[1], ry[1], rx[0], ry[0], initial_z, curve_prime,
 867                             ndigits);
 868
 869         for (i = num_bits - 2; i > 0; i--) {
 870                 nb = !vli_test_bit(scalar, i);
 871                 xycz_add_c(rx[1 - nb], ry[1 - nb], rx[nb], ry[nb], curve_prime,
 872                            ndigits);
 873                 xycz_add(rx[nb], ry[nb], rx[1 - nb], ry[1 - nb], curve_prime,
 874                          ndigits);
 875         }
 876
 877         nb = !vli_test_bit(scalar, 0);
 878         xycz_add_c(rx[1 - nb], ry[1 - nb], rx[nb], ry[nb], curve_prime,
 879                    ndigits);
 880
 881         /* Find final 1/Z value. */
 882         /* X1 - X0 */
 883         vli_mod_sub(z, rx[1], rx[0], curve_prime, ndigits);
 884         /* Yb * (X1 - X0) */
 885         vli_mod_mult_fast(z, z, ry[1 - nb], curve_prime, ndigits);
 886         /* xP * Yb * (X1 - X0) */
 887         vli_mod_mult_fast(z, z, point->x, curve_prime, ndigits);
 888
 889         /* 1 / (xP * Yb * (X1 - X0)) */
 890         vli_mod_inv(z, z, curve_prime, point->ndigits);
 891
 892         /* yP / (xP * Yb * (X1 - X0)) */
 893         vli_mod_mult_fast(z, z, point->y, curve_prime, ndigits);
 894         /* Xb * yP / (xP * Yb * (X1 - X0)) */
 895         vli_mod_mult_fast(z, z, rx[1 - nb], curve_prime, ndigits);
 896         /* End 1/Z calculation */
 897
 898         xycz_add(rx[nb], ry[nb], rx[1 - nb], ry[1 - nb], curve_prime, ndigits);
 899
 900         apply_z(rx[0], ry[0], z, curve_prime, ndigits);
 901
 902         vli_set(result->x, rx[0], ndigits);
 903         vli_set(result->y, ry[0], ndigits);
 904 }
 905
 906 static inline void ecc_swap_digits(const u64 *in, u64 *out,
 907                                    unsigned int ndigits)
 908 {
 909         int i;
 910
 911         for (i = 0; i < ndigits; i++)
 912                 out[i] = __swab64(in[ndigits - 1 - i]);
 913 }
 914
 915 static int __ecc_is_key_valid(const struct ecc_curve *curve,
 916                               const u64 *private_key, unsigned int ndigits)
 917 {
 918         u64 one[ECC_MAX_DIGITS] = { 1, };
 919         u64 res[ECC_MAX_DIGITS];
 920
 921         if (!private_key)
 922                 return -EINVAL;
 923
 924         if (curve->g.ndigits != ndigits)
 925                 return -EINVAL;
 926
 927         /* Make sure the private key is in the range [2, n-3]. */
 928         if (vli_cmp(one, private_key, ndigits) != -1)
 929                 return -EINVAL;
 930         vli_sub(res, curve->n, one, ndigits);
 931         vli_sub(res, res, one, ndigits);
 932         if (vli_cmp(res, private_key, ndigits) != 1)
 933                 return -EINVAL;
 934
 935         return 0;
 936 }
 937
 938 int ecc_is_key_valid(unsigned int curve_id, unsigned int ndigits,
 939                      const u64 *private_key, unsigned int private_key_len)
 940 {
 941         int nbytes;
 942         const struct ecc_curve *curve = ecc_get_curve(curve_id);
 943
 944         nbytes = ndigits << ECC_DIGITS_TO_BYTES_SHIFT;
 945
 946         if (private_key_len != nbytes)
 947                 return -EINVAL;
 948
 949         return __ecc_is_key_valid(curve, private_key, ndigits);
 950 }
 951
 952 /*
 953  * ECC private keys are generated using the method of extra random bits,
 954  * equivalent to that described in FIPS 186-4, Appendix B.4.1.
 955  *
 956  * d = (c mod(n–1)) + 1    where c is a string of random bits, 64 bits longer
 957  *                         than requested
 958  * 0 <= c mod(n-1) <= n-2  and implies that
 959  * 1 <= d <= n-1
 960  *
 961  * This method generates a private key uniformly distributed in the range
 962  * [1, n-1].
 963  */
 964 int ecc_gen_privkey(unsigned int curve_id, unsigned int ndigits, u64 *privkey)
 965 {
 966         const struct ecc_curve *curve = ecc_get_curve(curve_id);
 967         u64 priv[ECC_MAX_DIGITS];
 968         unsigned int nbytes = ndigits << ECC_DIGITS_TO_BYTES_SHIFT;
 969         unsigned int nbits = vli_num_bits(curve->n, ndigits);
 970         int err;
 971
 972         /* Check that N is included in Table 1 of FIPS 186-4, section 6.1.1 */
 973         if (nbits < 160 || ndigits > ARRAY_SIZE(priv))
 974                 return -EINVAL;
 975
 976         /*
 977          * FIPS 186-4 recommends that the private key should be obtained from a
 978          * RBG with a security strength equal to or greater than the security
 979          * strength associated with N.
 980          *
 981          * The maximum security strength identified by NIST SP800-57pt1r4 for
 982          * ECC is 256 (N >= 512).
 983          *
 984          * This condition is met by the default RNG because it selects a favored
 985          * DRBG with a security strength of 256.
 986          */
 987         if (crypto_get_default_rng())
 988                 return -EFAULT;
 989
 990         err = crypto_rng_get_bytes(crypto_default_rng, (u8 *)priv, nbytes);
 991         crypto_put_default_rng();
 992         if (err)
 993                 return err;
 994
 995         /* Make sure the private key is in the valid range. */
 996         if (__ecc_is_key_valid(curve, priv, ndigits))
 997                 return -EINVAL;
 998
 999         ecc_swap_digits(priv, privkey, ndigits);
1000
1001         return 0;
1002 }
1003
1004 int ecc_make_pub_key(unsigned int curve_id, unsigned int ndigits,
1005                      const u64 *private_key, u64 *public_key)
1006 {
1007         int ret = 0;
1008         struct ecc_point *pk;
1009         u64 priv[ECC_MAX_DIGITS];
1010         const struct ecc_curve *curve = ecc_get_curve(curve_id);
1011
1012         if (!private_key || !curve || ndigits > ARRAY_SIZE(priv)) {
1013                 ret = -EINVAL;
1014                 goto out;
1015         }
1016
1017         ecc_swap_digits(private_key, priv, ndigits);
1018
1019         pk = ecc_alloc_point(ndigits);
1020         if (!pk) {
1021                 ret = -ENOMEM;
1022                 goto out;
1023         }
1024
1025         ecc_point_mult(pk, &curve->g, priv, NULL, curve, ndigits);
1026         if (ecc_point_is_zero(pk)) {
1027                 ret = -EAGAIN;
1028                 goto err_free_point;
1029         }
1030
1031         ecc_swap_digits(pk->x, public_key, ndigits);
1032         ecc_swap_digits(pk->y, &public_key[ndigits], ndigits);
1033
1034 err_free_point:
1035         ecc_free_point(pk);
1036 out:
1037         return ret;
1038 }
1039
1040 /* SP800-56A section 5.6.2.3.4 partial verification: ephemeral keys only */
1041 static int ecc_is_pubkey_valid_partial(const struct ecc_curve *curve,
1042                                        struct ecc_point *pk)
1043 {
1044         u64 yy[ECC_MAX_DIGITS], xxx[ECC_MAX_DIGITS], w[ECC_MAX_DIGITS];
1045
1046         /* Check 1: Verify key is not the zero point. */
1047         if (ecc_point_is_zero(pk))
1048                 return -EINVAL;
1049
1050         /* Check 2: Verify key is in the range [1, p-1]. */
1051         if (vli_cmp(curve->p, pk->x, pk->ndigits) != 1)
1052                 return -EINVAL;
1053         if (vli_cmp(curve->p, pk->y, pk->ndigits) != 1)
1054                 return -EINVAL;
1055
1056         /* Check 3: Verify that y^2 == (x^3 + a·x + b) mod p */
1057         vli_mod_square_fast(yy, pk->y, curve->p, pk->ndigits); /* y^2 */
1058         vli_mod_square_fast(xxx, pk->x, curve->p, pk->ndigits); /* x^2 */
1059         vli_mod_mult_fast(xxx, xxx, pk->x, curve->p, pk->ndigits); /* x^3 */
1060         vli_mod_mult_fast(w, curve->a, pk->x, curve->p, pk->ndigits); /* a·x */
1061         vli_mod_add(w, w, curve->b, curve->p, pk->ndigits); /* a·x + b */
1062         vli_mod_add(w, w, xxx, curve->p, pk->ndigits); /* x^3 + a·x + b */
1063         if (vli_cmp(yy, w, pk->ndigits) != 0) /* Equation */
1064                 return -EINVAL;
1065
1066         return 0;
1067
1068 }
1069
1070 int crypto_ecdh_shared_secret(unsigned int curve_id, unsigned int ndigits,
1071                               const u64 *private_key, const u64 *public_key,
1072                               u64 *secret)
1073 {
1074         int ret = 0;
1075         struct ecc_point *product, *pk;
1076         u64 priv[ECC_MAX_DIGITS];
1077         u64 rand_z[ECC_MAX_DIGITS];
1078         unsigned int nbytes;
1079         const struct ecc_curve *curve = ecc_get_curve(curve_id);
1080
1081         if (!private_key || !public_key || !curve ||
1082             ndigits > ARRAY_SIZE(priv) || ndigits > ARRAY_SIZE(rand_z)) {
1083                 ret = -EINVAL;
1084                 goto out;
1085         }
1086
1087         nbytes = ndigits << ECC_DIGITS_TO_BYTES_SHIFT;
1088
1089         get_random_bytes(rand_z, nbytes);
1090
1091         pk = ecc_alloc_point(ndigits);
1092         if (!pk) {
1093                 ret = -ENOMEM;
1094                 goto out;
1095         }
1096
1097         ecc_swap_digits(public_key, pk->x, ndigits);
1098         ecc_swap_digits(&public_key[ndigits], pk->y, ndigits);
1099         ret = ecc_is_pubkey_valid_partial(curve, pk);
1100         if (ret)
1101                 goto err_alloc_product;
1102
1103         ecc_swap_digits(private_key, priv, ndigits);
1104
1105         product = ecc_alloc_point(ndigits);
1106         if (!product) {
1107                 ret = -ENOMEM;
1108                 goto err_alloc_product;
1109         }
1110
1111         ecc_point_mult(product, pk, priv, rand_z, curve, ndigits);
1112
1113         ecc_swap_digits(product->x, secret, ndigits);
1114
1115         if (ecc_point_is_zero(product))
1116                 ret = -EFAULT;
1117
1118         ecc_free_point(product);
1119 err_alloc_product:
1120         ecc_free_point(pk);
1121 out:
1122         return ret;
1123 }