1 /*- 2 * SPDX-License-Identifier: BSD-2-Clause-FreeBSD 3 * 4 * Copyright (c) 2001 Dima Dorfman. 5 * All rights reserved. 6 * 7 * Redistribution and use in source and binary forms, with or without 8 * modification, are permitted provided that the following conditions 9 * are met: 10 * 1. Redistributions of source code must retain the above copyright 11 * notice, this list of conditions and the following disclaimer. 12 * 2. Redistributions in binary form must reproduce the above copyright 13 * notice, this list of conditions and the following disclaimer in the 14 * documentation and/or other materials provided with the distribution. 15 * 16 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 17 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 18 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 19 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 20 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 21 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 22 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 23 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 24 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 25 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 26 * SUCH DAMAGE. 27 */ 28 29 /* 30 * This is the traditional Berkeley MP library implemented in terms of 31 * the OpenSSL BIGNUM library. It was written to replace libgmp, and 32 * is meant to be as compatible with the latter as feasible. 33 * 34 * There seems to be a lack of documentation for the Berkeley MP 35 * interface. All I could find was libgmp documentation (which didn't 36 * talk about the semantics of the functions) and an old SunOS 4.1 37 * manual page from 1989. The latter wasn't very detailed, either, 38 * but at least described what the function's arguments were. In 39 * general the interface seems to be archaic, somewhat poorly 40 * designed, and poorly, if at all, documented. It is considered 41 * harmful. 42 * 43 * Miscellaneous notes on this implementation: 44 * 45 * - The SunOS manual page mentioned above indicates that if an error 46 * occurs, the library should "produce messages and core images." 47 * Given that most of the functions don't have return values (and 48 * thus no sane way of alerting the caller to an error), this seems 49 * reasonable. The MPERR and MPERRX macros call warn and warnx, 50 * respectively, then abort(). 51 * 52 * - All the functions which take an argument to be "filled in" 53 * assume that the argument has been initialized by one of the *tom() 54 * routines before being passed to it. I never saw this documented 55 * anywhere, but this seems to be consistent with the way this 56 * library is used. 57 * 58 * - msqrt() is the only routine which had to be implemented which 59 * doesn't have a close counterpart in the OpenSSL BIGNUM library. 60 * It was implemented by hand using Newton's recursive formula. 61 * Doing it this way, although more error-prone, has the positive 62 * sideaffect of testing a lot of other functions; if msqrt() 63 * produces the correct results, most of the other routines will as 64 * well. 65 * 66 * - Internal-use-only routines (i.e., those defined here statically 67 * and not in mp.h) have an underscore prepended to their name (this 68 * is more for aesthetical reasons than technical). All such 69 * routines take an extra argument, 'msg', that denotes what they 70 * should call themselves in an error message. This is so a user 71 * doesn't get an error message from a function they didn't call. 72 */ 73 74 #include <sys/cdefs.h> 75 __FBSDID("$FreeBSD$"); 76 77 #include <ctype.h> 78 #include <err.h> 79 #include <errno.h> 80 #include <stdio.h> 81 #include <stdlib.h> 82 #include <string.h> 83 84 #include <openssl/crypto.h> 85 #include <openssl/err.h> 86 87 #include "mp.h" 88 89 #define MPERR(s) do { warn s; abort(); } while (0) 90 #define MPERRX(s) do { warnx s; abort(); } while (0) 91 #define BN_ERRCHECK(msg, expr) do { \ 92 if (!(expr)) _bnerr(msg); \ 93 } while (0) 94 95 static void _bnerr(const char *); 96 static MINT *_dtom(const char *, const char *); 97 static MINT *_itom(const char *, short); 98 static void _madd(const char *, const MINT *, const MINT *, MINT *); 99 static int _mcmpa(const char *, const MINT *, const MINT *); 100 static void _mdiv(const char *, const MINT *, const MINT *, MINT *, MINT *, 101 BN_CTX *); 102 static void _mfree(const char *, MINT *); 103 static void _moveb(const char *, const BIGNUM *, MINT *); 104 static void _movem(const char *, const MINT *, MINT *); 105 static void _msub(const char *, const MINT *, const MINT *, MINT *); 106 static char *_mtod(const char *, const MINT *); 107 static char *_mtox(const char *, const MINT *); 108 static void _mult(const char *, const MINT *, const MINT *, MINT *, BN_CTX *); 109 static void _sdiv(const char *, const MINT *, short, MINT *, short *, BN_CTX *); 110 static MINT *_xtom(const char *, const char *); 111 112 /* 113 * Report an error from one of the BN_* functions using MPERRX. 114 */ 115 static void 116 _bnerr(const char *msg) 117 { 118 119 ERR_load_crypto_strings(); 120 MPERRX(("%s: %s", msg, ERR_reason_error_string(ERR_get_error()))); 121 } 122 123 /* 124 * Convert a decimal string to an MINT. 125 */ 126 static MINT * 127 _dtom(const char *msg, const char *s) 128 { 129 MINT *mp; 130 131 mp = malloc(sizeof(*mp)); 132 if (mp == NULL) 133 MPERR(("%s", msg)); 134 mp->bn = BN_new(); 135 if (mp->bn == NULL) 136 _bnerr(msg); 137 BN_ERRCHECK(msg, BN_dec2bn(&mp->bn, s)); 138 return (mp); 139 } 140 141 /* 142 * Compute the greatest common divisor of mp1 and mp2; result goes in rmp. 143 */ 144 void 145 mp_gcd(const MINT *mp1, const MINT *mp2, MINT *rmp) 146 { 147 BIGNUM *b; 148 BN_CTX *c; 149 150 b = NULL; 151 c = BN_CTX_new(); 152 if (c != NULL) 153 b = BN_new(); 154 if (c == NULL || b == NULL) 155 _bnerr("gcd"); 156 BN_ERRCHECK("gcd", BN_gcd(b, mp1->bn, mp2->bn, c)); 157 _moveb("gcd", b, rmp); 158 BN_free(b); 159 BN_CTX_free(c); 160 } 161 162 /* 163 * Make an MINT out of a short integer. Return value must be mfree()'d. 164 */ 165 static MINT * 166 _itom(const char *msg, short n) 167 { 168 MINT *mp; 169 char *s; 170 171 asprintf(&s, "%x", n); 172 if (s == NULL) 173 MPERR(("%s", msg)); 174 mp = _xtom(msg, s); 175 free(s); 176 return (mp); 177 } 178 179 MINT * 180 mp_itom(short n) 181 { 182 183 return (_itom("itom", n)); 184 } 185 186 /* 187 * Compute rmp=mp1+mp2. 188 */ 189 static void 190 _madd(const char *msg, const MINT *mp1, const MINT *mp2, MINT *rmp) 191 { 192 BIGNUM *b; 193 194 b = BN_new(); 195 if (b == NULL) 196 _bnerr(msg); 197 BN_ERRCHECK(msg, BN_add(b, mp1->bn, mp2->bn)); 198 _moveb(msg, b, rmp); 199 BN_free(b); 200 } 201 202 void 203 mp_madd(const MINT *mp1, const MINT *mp2, MINT *rmp) 204 { 205 206 _madd("madd", mp1, mp2, rmp); 207 } 208 209 /* 210 * Return -1, 0, or 1 if mp1<mp2, mp1==mp2, or mp1>mp2, respectivley. 211 */ 212 int 213 mp_mcmp(const MINT *mp1, const MINT *mp2) 214 { 215 216 return (BN_cmp(mp1->bn, mp2->bn)); 217 } 218 219 /* 220 * Same as mcmp but compares absolute values. 221 */ 222 static int 223 _mcmpa(const char *msg __unused, const MINT *mp1, const MINT *mp2) 224 { 225 226 return (BN_ucmp(mp1->bn, mp2->bn)); 227 } 228 229 /* 230 * Compute qmp=nmp/dmp and rmp=nmp%dmp. 231 */ 232 static void 233 _mdiv(const char *msg, const MINT *nmp, const MINT *dmp, MINT *qmp, MINT *rmp, 234 BN_CTX *c) 235 { 236 BIGNUM *q, *r; 237 238 q = NULL; 239 r = BN_new(); 240 if (r != NULL) 241 q = BN_new(); 242 if (r == NULL || q == NULL) 243 _bnerr(msg); 244 BN_ERRCHECK(msg, BN_div(q, r, nmp->bn, dmp->bn, c)); 245 _moveb(msg, q, qmp); 246 _moveb(msg, r, rmp); 247 BN_free(q); 248 BN_free(r); 249 } 250 251 void 252 mp_mdiv(const MINT *nmp, const MINT *dmp, MINT *qmp, MINT *rmp) 253 { 254 BN_CTX *c; 255 256 c = BN_CTX_new(); 257 if (c == NULL) 258 _bnerr("mdiv"); 259 _mdiv("mdiv", nmp, dmp, qmp, rmp, c); 260 BN_CTX_free(c); 261 } 262 263 /* 264 * Free memory associated with an MINT. 265 */ 266 static void 267 _mfree(const char *msg __unused, MINT *mp) 268 { 269 270 BN_clear(mp->bn); 271 BN_free(mp->bn); 272 free(mp); 273 } 274 275 void 276 mp_mfree(MINT *mp) 277 { 278 279 _mfree("mfree", mp); 280 } 281 282 /* 283 * Read an integer from standard input and stick the result in mp. 284 * The input is treated to be in base 10. This must be the silliest 285 * API in existence; why can't the program read in a string and call 286 * xtom()? (Or if base 10 is desires, perhaps dtom() could be 287 * exported.) 288 */ 289 void 290 mp_min(MINT *mp) 291 { 292 MINT *rmp; 293 char *line, *nline; 294 size_t linelen; 295 296 line = fgetln(stdin, &linelen); 297 if (line == NULL) 298 MPERR(("min")); 299 nline = malloc(linelen + 1); 300 if (nline == NULL) 301 MPERR(("min")); 302 memcpy(nline, line, linelen); 303 nline[linelen] = '\0'; 304 rmp = _dtom("min", nline); 305 _movem("min", rmp, mp); 306 _mfree("min", rmp); 307 free(nline); 308 } 309 310 /* 311 * Print the value of mp to standard output in base 10. See blurb 312 * above min() for why this is so useless. 313 */ 314 void 315 mp_mout(const MINT *mp) 316 { 317 char *s; 318 319 s = _mtod("mout", mp); 320 printf("%s", s); 321 free(s); 322 } 323 324 /* 325 * Set the value of tmp to the value of smp (i.e., tmp=smp). 326 */ 327 void 328 mp_move(const MINT *smp, MINT *tmp) 329 { 330 331 _movem("move", smp, tmp); 332 } 333 334 335 /* 336 * Internal routine to set the value of tmp to that of sbp. 337 */ 338 static void 339 _moveb(const char *msg, const BIGNUM *sbp, MINT *tmp) 340 { 341 342 BN_ERRCHECK(msg, BN_copy(tmp->bn, sbp)); 343 } 344 345 /* 346 * Internal routine to set the value of tmp to that of smp. 347 */ 348 static void 349 _movem(const char *msg, const MINT *smp, MINT *tmp) 350 { 351 352 BN_ERRCHECK(msg, BN_copy(tmp->bn, smp->bn)); 353 } 354 355 /* 356 * Compute the square root of nmp and put the result in xmp. The 357 * remainder goes in rmp. Should satisfy: rmp=nmp-(xmp*xmp). 358 * 359 * Note that the OpenSSL BIGNUM library does not have a square root 360 * function, so this had to be implemented by hand using Newton's 361 * recursive formula: 362 * 363 * x = (x + (n / x)) / 2 364 * 365 * where x is the square root of the positive number n. In the 366 * beginning, x should be a reasonable guess, but the value 1, 367 * although suboptimal, works, too; this is that is used below. 368 */ 369 void 370 mp_msqrt(const MINT *nmp, MINT *xmp, MINT *rmp) 371 { 372 BN_CTX *c; 373 MINT *tolerance; 374 MINT *ox, *x; 375 MINT *z1, *z2, *z3; 376 short i; 377 378 c = BN_CTX_new(); 379 if (c == NULL) 380 _bnerr("msqrt"); 381 tolerance = _itom("msqrt", 1); 382 x = _itom("msqrt", 1); 383 ox = _itom("msqrt", 0); 384 z1 = _itom("msqrt", 0); 385 z2 = _itom("msqrt", 0); 386 z3 = _itom("msqrt", 0); 387 do { 388 _movem("msqrt", x, ox); 389 _mdiv("msqrt", nmp, x, z1, z2, c); 390 _madd("msqrt", x, z1, z2); 391 _sdiv("msqrt", z2, 2, x, &i, c); 392 _msub("msqrt", ox, x, z3); 393 } while (_mcmpa("msqrt", z3, tolerance) == 1); 394 _movem("msqrt", x, xmp); 395 _mult("msqrt", x, x, z1, c); 396 _msub("msqrt", nmp, z1, z2); 397 _movem("msqrt", z2, rmp); 398 _mfree("msqrt", tolerance); 399 _mfree("msqrt", ox); 400 _mfree("msqrt", x); 401 _mfree("msqrt", z1); 402 _mfree("msqrt", z2); 403 _mfree("msqrt", z3); 404 BN_CTX_free(c); 405 } 406 407 /* 408 * Compute rmp=mp1-mp2. 409 */ 410 static void 411 _msub(const char *msg, const MINT *mp1, const MINT *mp2, MINT *rmp) 412 { 413 BIGNUM *b; 414 415 b = BN_new(); 416 if (b == NULL) 417 _bnerr(msg); 418 BN_ERRCHECK(msg, BN_sub(b, mp1->bn, mp2->bn)); 419 _moveb(msg, b, rmp); 420 BN_free(b); 421 } 422 423 void 424 mp_msub(const MINT *mp1, const MINT *mp2, MINT *rmp) 425 { 426 427 _msub("msub", mp1, mp2, rmp); 428 } 429 430 /* 431 * Return a decimal representation of mp. Return value must be 432 * free()'d. 433 */ 434 static char * 435 _mtod(const char *msg, const MINT *mp) 436 { 437 char *s, *s2; 438 439 s = BN_bn2dec(mp->bn); 440 if (s == NULL) 441 _bnerr(msg); 442 asprintf(&s2, "%s", s); 443 if (s2 == NULL) 444 MPERR(("%s", msg)); 445 OPENSSL_free(s); 446 return (s2); 447 } 448 449 /* 450 * Return a hexadecimal representation of mp. Return value must be 451 * free()'d. 452 */ 453 static char * 454 _mtox(const char *msg, const MINT *mp) 455 { 456 char *p, *s, *s2; 457 int len; 458 459 s = BN_bn2hex(mp->bn); 460 if (s == NULL) 461 _bnerr(msg); 462 asprintf(&s2, "%s", s); 463 if (s2 == NULL) 464 MPERR(("%s", msg)); 465 OPENSSL_free(s); 466 467 /* 468 * This is a kludge for libgmp compatibility. The latter's 469 * implementation of this function returns lower-case letters, 470 * but BN_bn2hex returns upper-case. Some programs (e.g., 471 * newkey(1)) are sensitive to this. Although it's probably 472 * their fault, it's nice to be compatible. 473 */ 474 len = strlen(s2); 475 for (p = s2; p < s2 + len; p++) 476 *p = tolower(*p); 477 478 return (s2); 479 } 480 481 char * 482 mp_mtox(const MINT *mp) 483 { 484 485 return (_mtox("mtox", mp)); 486 } 487 488 /* 489 * Compute rmp=mp1*mp2. 490 */ 491 static void 492 _mult(const char *msg, const MINT *mp1, const MINT *mp2, MINT *rmp, BN_CTX *c) 493 { 494 BIGNUM *b; 495 496 b = BN_new(); 497 if (b == NULL) 498 _bnerr(msg); 499 BN_ERRCHECK(msg, BN_mul(b, mp1->bn, mp2->bn, c)); 500 _moveb(msg, b, rmp); 501 BN_free(b); 502 } 503 504 void 505 mp_mult(const MINT *mp1, const MINT *mp2, MINT *rmp) 506 { 507 BN_CTX *c; 508 509 c = BN_CTX_new(); 510 if (c == NULL) 511 _bnerr("mult"); 512 _mult("mult", mp1, mp2, rmp, c); 513 BN_CTX_free(c); 514 } 515 516 /* 517 * Compute rmp=(bmp^emp)mod mmp. (Note that here and above rpow() '^' 518 * means 'raise to power', not 'bitwise XOR'.) 519 */ 520 void 521 mp_pow(const MINT *bmp, const MINT *emp, const MINT *mmp, MINT *rmp) 522 { 523 BIGNUM *b; 524 BN_CTX *c; 525 526 b = NULL; 527 c = BN_CTX_new(); 528 if (c != NULL) 529 b = BN_new(); 530 if (c == NULL || b == NULL) 531 _bnerr("pow"); 532 BN_ERRCHECK("pow", BN_mod_exp(b, bmp->bn, emp->bn, mmp->bn, c)); 533 _moveb("pow", b, rmp); 534 BN_free(b); 535 BN_CTX_free(c); 536 } 537 538 /* 539 * Compute rmp=bmp^e. (See note above pow().) 540 */ 541 void 542 mp_rpow(const MINT *bmp, short e, MINT *rmp) 543 { 544 MINT *emp; 545 BIGNUM *b; 546 BN_CTX *c; 547 548 b = NULL; 549 c = BN_CTX_new(); 550 if (c != NULL) 551 b = BN_new(); 552 if (c == NULL || b == NULL) 553 _bnerr("rpow"); 554 emp = _itom("rpow", e); 555 BN_ERRCHECK("rpow", BN_exp(b, bmp->bn, emp->bn, c)); 556 _moveb("rpow", b, rmp); 557 _mfree("rpow", emp); 558 BN_free(b); 559 BN_CTX_free(c); 560 } 561 562 /* 563 * Compute qmp=nmp/d and ro=nmp%d. 564 */ 565 static void 566 _sdiv(const char *msg, const MINT *nmp, short d, MINT *qmp, short *ro, 567 BN_CTX *c) 568 { 569 MINT *dmp, *rmp; 570 BIGNUM *q, *r; 571 char *s; 572 573 r = NULL; 574 q = BN_new(); 575 if (q != NULL) 576 r = BN_new(); 577 if (q == NULL || r == NULL) 578 _bnerr(msg); 579 dmp = _itom(msg, d); 580 rmp = _itom(msg, 0); 581 BN_ERRCHECK(msg, BN_div(q, r, nmp->bn, dmp->bn, c)); 582 _moveb(msg, q, qmp); 583 _moveb(msg, r, rmp); 584 s = _mtox(msg, rmp); 585 errno = 0; 586 *ro = strtol(s, NULL, 16); 587 if (errno != 0) 588 MPERR(("%s underflow or overflow", msg)); 589 free(s); 590 _mfree(msg, dmp); 591 _mfree(msg, rmp); 592 BN_free(r); 593 BN_free(q); 594 } 595 596 void 597 mp_sdiv(const MINT *nmp, short d, MINT *qmp, short *ro) 598 { 599 BN_CTX *c; 600 601 c = BN_CTX_new(); 602 if (c == NULL) 603 _bnerr("sdiv"); 604 _sdiv("sdiv", nmp, d, qmp, ro, c); 605 BN_CTX_free(c); 606 } 607 608 /* 609 * Convert a hexadecimal string to an MINT. 610 */ 611 static MINT * 612 _xtom(const char *msg, const char *s) 613 { 614 MINT *mp; 615 616 mp = malloc(sizeof(*mp)); 617 if (mp == NULL) 618 MPERR(("%s", msg)); 619 mp->bn = BN_new(); 620 if (mp->bn == NULL) 621 _bnerr(msg); 622 BN_ERRCHECK(msg, BN_hex2bn(&mp->bn, s)); 623 return (mp); 624 } 625 626 MINT * 627 mp_xtom(const char *s) 628 { 629 630 return (_xtom("xtom", s)); 631 } 632