contrib/bearssl/test/test_math.c

   1 /*
   2  * Copyright (c) 2016 Thomas Pornin <pornin@bolet.org>
   3  *
   4  * Permission is hereby granted, free of charge, to any person obtaining
   5  * a copy of this software and associated documentation files (the
   6  * "Software"), to deal in the Software without restriction, including
   7  * without limitation the rights to use, copy, modify, merge, publish,
   8  * distribute, sublicense, and/or sell copies of the Software, and to
   9  * permit persons to whom the Software is furnished to do so, subject to
  10  * the following conditions:
  11  *
  12  * The above copyright notice and this permission notice shall be
  13  * included in all copies or substantial portions of the Software.
  14  *
  15  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  16  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  17  * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  18  * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
  19  * BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
  20  * ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
  21  * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  22  * SOFTWARE.
  23  */
  24
  25 #include <stdio.h>
  26 #include <stdlib.h>
  27 #include <string.h>
  28 #include <stdarg.h>
  29 #include <time.h>
  30
  31 #include <gmp.h>
  32
  33 #include "bearssl.h"
  34 #include "inner.h"
  35
  36 /*
  37  * Pointers to implementations.
  38  */
  39 typedef struct {
  40         uint32_t word_size;
  41         void (*zero)(uint32_t *x, uint32_t bit_len);
  42         void (*decode)(uint32_t *x, const void *src, size_t len);
  43         uint32_t (*decode_mod)(uint32_t *x,
  44                 const void *src, size_t len, const uint32_t *m);
  45         void (*reduce)(uint32_t *x, const uint32_t *a, const uint32_t *m);
  46         void (*decode_reduce)(uint32_t *x,
  47                 const void *src, size_t len, const uint32_t *m);
  48         void (*encode)(void *dst, size_t len, const uint32_t *x);
  49         uint32_t (*add)(uint32_t *a, const uint32_t *b, uint32_t ctl);
  50         uint32_t (*sub)(uint32_t *a, const uint32_t *b, uint32_t ctl);
  51         uint32_t (*ninv)(uint32_t x);
  52         void (*montymul)(uint32_t *d, const uint32_t *x, const uint32_t *y,
  53                 const uint32_t *m, uint32_t m0i);
  54         void (*to_monty)(uint32_t *x, const uint32_t *m);
  55         void (*from_monty)(uint32_t *x, const uint32_t *m, uint32_t m0i);
  56         void (*modpow)(uint32_t *x, const unsigned char *e, size_t elen,
  57                 const uint32_t *m, uint32_t m0i, uint32_t *t1, uint32_t *t2);
  58 } int_impl;
  59
  60 static const int_impl i31_impl = {
  61         31,
  62         &br_i31_zero,
  63         &br_i31_decode,
  64         &br_i31_decode_mod,
  65         &br_i31_reduce,
  66         &br_i31_decode_reduce,
  67         &br_i31_encode,
  68         &br_i31_add,
  69         &br_i31_sub,
  70         &br_i31_ninv31,
  71         &br_i31_montymul,
  72         &br_i31_to_monty,
  73         &br_i31_from_monty,
  74         &br_i31_modpow
  75 };
  76 static const int_impl i32_impl = {
  77         32,
  78         &br_i32_zero,
  79         &br_i32_decode,
  80         &br_i32_decode_mod,
  81         &br_i32_reduce,
  82         &br_i32_decode_reduce,
  83         &br_i32_encode,
  84         &br_i32_add,
  85         &br_i32_sub,
  86         &br_i32_ninv32,
  87         &br_i32_montymul,
  88         &br_i32_to_monty,
  89         &br_i32_from_monty,
  90         &br_i32_modpow
  91 };
  92
  93 static const int_impl *impl;
  94
  95 static gmp_randstate_t RNG;
  96
  97 /*
  98  * Get a random prime of length 'size' bits. This function also guarantees
  99  * that x-1 is not a multiple of 65537.
 100  */
 101 static void
 102 rand_prime(mpz_t x, int size)
 103 {
 104         for (;;) {
 105                 mpz_urandomb(x, RNG, size - 1);
 106                 mpz_setbit(x, 0);
 107                 mpz_setbit(x, size - 1);
 108                 if (mpz_probab_prime_p(x, 50)) {
 109                         mpz_sub_ui(x, x, 1);
 110                         if (mpz_divisible_ui_p(x, 65537)) {
 111                                 continue;
 112                         }
 113                         mpz_add_ui(x, x, 1);
 114                         return;
 115                 }
 116         }
 117 }
 118
 119 /*
 120  * Print out a GMP integer (for debug).
 121  */
 122 static void
 123 print_z(mpz_t z)
 124 {
 125         unsigned char zb[1000];
 126         size_t zlen, k;
 127
 128         mpz_export(zb, &zlen, 1, 1, 0, 0, z);
 129         if (zlen == 0) {
 130                 printf(" 00");
 131                 return;
 132         }
 133         if ((zlen & 3) != 0) {
 134                 k = 4 - (zlen & 3);
 135                 memmove(zb + k, zb, zlen);
 136                 memset(zb, 0, k);
 137                 zlen += k;
 138         }
 139         for (k = 0; k < zlen; k += 4) {
 140                 printf(" %02X%02X%02X%02X",
 141                         zb[k], zb[k + 1], zb[k + 2], zb[k + 3]);
 142         }
 143 }
 144
 145 /*
 146  * Print out an i31 or i32 integer (for debug).
 147  */
 148 static void
 149 print_u(uint32_t *x)
 150 {
 151         size_t k;
 152
 153         if (x[0] == 0) {
 154                 printf(" 00000000 (0, 0)");
 155                 return;
 156         }
 157         for (k = (x[0] + 31) >> 5; k > 0; k --) {
 158                 printf(" %08lX", (unsigned long)x[k]);
 159         }
 160         printf(" (%u, %u)", (unsigned)(x[0] >> 5), (unsigned)(x[0] & 31));
 161 }
 162
 163 /*
 164  * Check that an i31/i32 number and a GMP number are equal.
 165  */
 166 static void
 167 check_eqz(uint32_t *x, mpz_t z)
 168 {
 169         unsigned char xb[1000];
 170         unsigned char zb[1000];
 171         size_t xlen, zlen;
 172         int good;
 173
 174         xlen = ((x[0] + 31) & ~(uint32_t)31) >> 3;
 175         impl->encode(xb, xlen, x);
 176         mpz_export(zb, &zlen, 1, 1, 0, 0, z);
 177         good = 1;
 178         if (xlen < zlen) {
 179                 good = 0;
 180         } else if (xlen > zlen) {
 181                 size_t u;
 182
 183                 for (u = xlen; u > zlen; u --) {
 184                         if (xb[xlen - u] != 0) {
 185                                 good = 0;
 186                                 break;
 187                         }
 188                 }
 189         }
 190         good = good && memcmp(xb + xlen - zlen, zb, zlen) == 0;
 191         if (!good) {
 192                 size_t u;
 193
 194                 printf("Mismatch:\n");
 195                 printf("  x = ");
 196                 print_u(x);
 197                 printf("\n");
 198                 printf("  ex = ");
 199                 for (u = 0; u < xlen; u ++) {
 200                         printf("%02X", xb[u]);
 201                 }
 202                 printf("\n");
 203                 printf("  z = ");
 204                 print_z(z);
 205                 printf("\n");
 206                 exit(EXIT_FAILURE);
 207         }
 208 }
 209
 210 /* obsolete
 211 static void
 212 mp_to_br(uint32_t *mx, uint32_t x_bitlen, mpz_t x)
 213 {
 214         uint32_t x_ebitlen;
 215         size_t xlen;
 216
 217         if (mpz_sizeinbase(x, 2) > x_bitlen) {
 218                 abort();
 219         }
 220         x_ebitlen = ((x_bitlen / 31) << 5) + (x_bitlen % 31);
 221         br_i31_zero(mx, x_ebitlen);
 222         mpz_export(mx + 1, &xlen, -1, sizeof *mx, 0, 1, x);
 223 }
 224 */
 225
 226 static void
 227 test_modint(void)
 228 {
 229         int i, j, k;
 230         mpz_t p, a, b, v, t1;
 231
 232         printf("Test modular integers: ");
 233         fflush(stdout);
 234
 235         gmp_randinit_mt(RNG);
 236         mpz_init(p);
 237         mpz_init(a);
 238         mpz_init(b);
 239         mpz_init(v);
 240         mpz_init(t1);
 241         mpz_set_ui(t1, (unsigned long)time(NULL));
 242         gmp_randseed(RNG, t1);
 243         for (k = 2; k <= 128; k ++) {
 244                 for (i = 0; i < 10; i ++) {
 245                         unsigned char ep[100], ea[100], eb[100], ev[100];
 246                         size_t plen, alen, blen, vlen;
 247                         uint32_t mp[40], ma[40], mb[40], mv[60], mx[100];
 248                         uint32_t mt1[40], mt2[40], mt3[40];
 249                         uint32_t ctl;
 250                         uint32_t mp0i;
 251
 252                         rand_prime(p, k);
 253                         mpz_urandomm(a, RNG, p);
 254                         mpz_urandomm(b, RNG, p);
 255                         mpz_urandomb(v, RNG, k + 60);
 256                         if (mpz_sgn(b) == 0) {
 257                                 mpz_set_ui(b, 1);
 258                         }
 259                         mpz_export(ep, &plen, 1, 1, 0, 0, p);
 260                         mpz_export(ea, &alen, 1, 1, 0, 0, a);
 261                         mpz_export(eb, &blen, 1, 1, 0, 0, b);
 262                         mpz_export(ev, &vlen, 1, 1, 0, 0, v);
 263
 264                         impl->decode(mp, ep, plen);
 265                         if (impl->decode_mod(ma, ea, alen, mp) != 1) {
 266                                 printf("Decode error\n");
 267                                 printf("  ea = ");
 268                                 print_z(a);
 269                                 printf("\n");
 270                                 printf("  p = ");
 271                                 print_u(mp);
 272                                 printf("\n");
 273                                 exit(EXIT_FAILURE);
 274                         }
 275                         mp0i = impl->ninv(mp[1]);
 276                         if (impl->decode_mod(mb, eb, blen, mp) != 1) {
 277                                 printf("Decode error\n");
 278                                 printf("  eb = ");
 279                                 print_z(b);
 280                                 printf("\n");
 281                                 printf("  p = ");
 282                                 print_u(mp);
 283                                 printf("\n");
 284                                 exit(EXIT_FAILURE);
 285                         }
 286                         impl->decode(mv, ev, vlen);
 287                         check_eqz(mp, p);
 288                         check_eqz(ma, a);
 289                         check_eqz(mb, b);
 290                         check_eqz(mv, v);
 291
 292                         impl->decode_mod(ma, ea, alen, mp);
 293                         impl->decode_mod(mb, eb, blen, mp);
 294                         ctl = impl->add(ma, mb, 1);
 295                         ctl |= impl->sub(ma, mp, 0) ^ (uint32_t)1;
 296                         impl->sub(ma, mp, ctl);
 297                         mpz_add(t1, a, b);
 298                         mpz_mod(t1, t1, p);
 299                         check_eqz(ma, t1);
 300
 301                         impl->decode_mod(ma, ea, alen, mp);
 302                         impl->decode_mod(mb, eb, blen, mp);
 303                         impl->add(ma, mp, impl->sub(ma, mb, 1));
 304                         mpz_sub(t1, a, b);
 305                         mpz_mod(t1, t1, p);
 306                         check_eqz(ma, t1);
 307
 308                         impl->decode_reduce(ma, ev, vlen, mp);
 309                         mpz_mod(t1, v, p);
 310                         check_eqz(ma, t1);
 311
 312                         impl->decode(mv, ev, vlen);
 313                         impl->reduce(ma, mv, mp);
 314                         mpz_mod(t1, v, p);
 315                         check_eqz(ma, t1);
 316
 317                         impl->decode_mod(ma, ea, alen, mp);
 318                         impl->to_monty(ma, mp);
 319                         mpz_mul_2exp(t1, a, ((k + impl->word_size - 1)
 320                                 / impl->word_size) * impl->word_size);
 321                         mpz_mod(t1, t1, p);
 322                         check_eqz(ma, t1);
 323                         impl->from_monty(ma, mp, mp0i);
 324                         check_eqz(ma, a);
 325
 326                         impl->decode_mod(ma, ea, alen, mp);
 327                         impl->decode_mod(mb, eb, blen, mp);
 328                         impl->to_monty(ma, mp);
 329                         impl->montymul(mt1, ma, mb, mp, mp0i);
 330                         mpz_mul(t1, a, b);
 331                         mpz_mod(t1, t1, p);
 332                         check_eqz(mt1, t1);
 333
 334                         impl->decode_mod(ma, ea, alen, mp);
 335                         impl->modpow(ma, ev, vlen, mp, mp0i, mt1, mt2);
 336                         mpz_powm(t1, a, v, p);
 337                         check_eqz(ma, t1);
 338
 339                         /*
 340                         br_modint_decode(ma, mp, ea, alen);
 341                         br_modint_decode(mb, mp, eb, blen);
 342                         if (!br_modint_div(ma, mb, mp, mt1, mt2, mt3)) {
 343                                 fprintf(stderr, "division failed\n");
 344                                 exit(EXIT_FAILURE);
 345                         }
 346                         mpz_sub_ui(t1, p, 2);
 347                         mpz_powm(t1, b, t1, p);
 348                         mpz_mul(t1, a, t1);
 349                         mpz_mod(t1, t1, p);
 350                         check_eqz(ma, t1);
 351
 352                         br_modint_decode(ma, mp, ea, alen);
 353                         br_modint_decode(mb, mp, eb, blen);
 354                         for (j = 0; j <= (2 * k + 5); j ++) {
 355                                 br_int_add(mx, j, ma, mb);
 356                                 mpz_add(t1, a, b);
 357                                 mpz_tdiv_r_2exp(t1, t1, j);
 358                                 check_eqz(mx, t1);
 359
 360                                 br_int_mul(mx, j, ma, mb);
 361                                 mpz_mul(t1, a, b);
 362                                 mpz_tdiv_r_2exp(t1, t1, j);
 363                                 check_eqz(mx, t1);
 364                         }
 365                         */
 366                 }
 367                 printf(".");
 368                 fflush(stdout);
 369         }
 370         mpz_clear(p);
 371         mpz_clear(a);
 372         mpz_clear(b);
 373         mpz_clear(v);
 374         mpz_clear(t1);
 375
 376         printf(" done.\n");
 377         fflush(stdout);
 378 }
 379
 380 #if 0
 381 static void
 382 test_RSA_core(void)
 383 {
 384         int i, j, k;
 385         mpz_t n, e, d, p, q, dp, dq, iq, t1, t2, phi;
 386
 387         printf("Test RSA core: ");
 388         fflush(stdout);
 389
 390         gmp_randinit_mt(RNG);
 391         mpz_init(n);
 392         mpz_init(e);
 393         mpz_init(d);
 394         mpz_init(p);
 395         mpz_init(q);
 396         mpz_init(dp);
 397         mpz_init(dq);
 398         mpz_init(iq);
 399         mpz_init(t1);
 400         mpz_init(t2);
 401         mpz_init(phi);
 402         mpz_set_ui(t1, (unsigned long)time(NULL));
 403         gmp_randseed(RNG, t1);
 404
 405         /*
 406          * To test corner cases, we want to try RSA keys such that the
 407          * lengths of both factors can be arbitrary modulo 2^32. Factors
 408          * p and q need not be of the same length; p can be greater than
 409          * q and q can be greater than p.
 410          *
 411          * To keep computation time reasonable, we use p and q factors of
 412          * less than 128 bits; this is way too small for secure RSA,
 413          * but enough to exercise all code paths (since we work only with
 414          * 32-bit words).
 415          */
 416         for (i = 64; i <= 96; i ++) {
 417                 rand_prime(p, i);
 418                 for (j = i - 33; j <= i + 33; j ++) {
 419                         uint32_t mp[40], mq[40], mdp[40], mdq[40], miq[40];
 420
 421                         /*
 422                          * Generate a RSA key pair, with p of length i bits,
 423                          * and q of length j bits.
 424                          */
 425                         do {
 426                                 rand_prime(q, j);
 427                         } while (mpz_cmp(p, q) == 0);
 428                         mpz_mul(n, p, q);
 429                         mpz_set_ui(e, 65537);
 430                         mpz_sub_ui(t1, p, 1);
 431                         mpz_sub_ui(t2, q, 1);
 432                         mpz_mul(phi, t1, t2);
 433                         mpz_invert(d, e, phi);
 434                         mpz_mod(dp, d, t1);
 435                         mpz_mod(dq, d, t2);
 436                         mpz_invert(iq, q, p);
 437
 438                         /*
 439                          * Convert the key pair elements to BearSSL arrays.
 440                          */
 441                         mp_to_br(mp, mpz_sizeinbase(p, 2), p);
 442                         mp_to_br(mq, mpz_sizeinbase(q, 2), q);
 443                         mp_to_br(mdp, mpz_sizeinbase(dp, 2), dp);
 444                         mp_to_br(mdq, mpz_sizeinbase(dq, 2), dq);
 445                         mp_to_br(miq, mp[0], iq);
 446
 447                         /*
 448                          * Compute and check ten public/private operations.
 449                          */
 450                         for (k = 0; k < 10; k ++) {
 451                                 uint32_t mx[40];
 452
 453                                 mpz_urandomm(t1, RNG, n);
 454                                 mpz_powm(t2, t1, e, n);
 455                                 mp_to_br(mx, mpz_sizeinbase(n, 2), t2);
 456                                 br_rsa_private_core(mx, mp, mq, mdp, mdq, miq);
 457                                 check_eqz(mx, t1);
 458                         }
 459                 }
 460                 printf(".");
 461                 fflush(stdout);
 462         }
 463
 464         printf(" done.\n");
 465         fflush(stdout);
 466 }
 467 #endif
 468
 469 int
 470 main(void)
 471 {
 472         printf("===== i32 ======\n");
 473         impl = &i32_impl;
 474         test_modint();
 475         printf("===== i31 ======\n");
 476         impl = &i31_impl;
 477         test_modint();
 478         /*
 479         test_RSA_core();
 480         */
 481         return 0;
 482 }