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 c = BN_CTX_new(); 151 if (c == NULL) 152 _bnerr("gcd"); 153 BN_init(&b); 154 BN_ERRCHECK("gcd", BN_gcd(&b, mp1->bn, mp2->bn, c)); 155 _moveb("gcd", &b, rmp); 156 BN_free(&b); 157 BN_CTX_free(c); 158 } 159 160 /* 161 * Make an MINT out of a short integer. Return value must be mfree()'d. 162 */ 163 static MINT * 164 _itom(const char *msg, short n) 165 { 166 MINT *mp; 167 char *s; 168 169 asprintf(&s, "%x", n); 170 if (s == NULL) 171 MPERR(("%s", msg)); 172 mp = _xtom(msg, s); 173 free(s); 174 return (mp); 175 } 176 177 MINT * 178 mp_itom(short n) 179 { 180 181 return (_itom("itom", n)); 182 } 183 184 /* 185 * Compute rmp=mp1+mp2. 186 */ 187 static void 188 _madd(const char *msg, const MINT *mp1, const MINT *mp2, MINT *rmp) 189 { 190 BIGNUM b; 191 192 BN_init(&b); 193 BN_ERRCHECK(msg, BN_add(&b, mp1->bn, mp2->bn)); 194 _moveb(msg, &b, rmp); 195 BN_free(&b); 196 } 197 198 void 199 mp_madd(const MINT *mp1, const MINT *mp2, MINT *rmp) 200 { 201 202 _madd("madd", mp1, mp2, rmp); 203 } 204 205 /* 206 * Return -1, 0, or 1 if mp1<mp2, mp1==mp2, or mp1>mp2, respectivley. 207 */ 208 int 209 mp_mcmp(const MINT *mp1, const MINT *mp2) 210 { 211 212 return (BN_cmp(mp1->bn, mp2->bn)); 213 } 214 215 /* 216 * Same as mcmp but compares absolute values. 217 */ 218 static int 219 _mcmpa(const char *msg __unused, const MINT *mp1, const MINT *mp2) 220 { 221 222 return (BN_ucmp(mp1->bn, mp2->bn)); 223 } 224 225 /* 226 * Compute qmp=nmp/dmp and rmp=nmp%dmp. 227 */ 228 static void 229 _mdiv(const char *msg, const MINT *nmp, const MINT *dmp, MINT *qmp, MINT *rmp, 230 BN_CTX *c) 231 { 232 BIGNUM q, r; 233 234 BN_init(&r); 235 BN_init(&q); 236 BN_ERRCHECK(msg, BN_div(&q, &r, nmp->bn, dmp->bn, c)); 237 _moveb(msg, &q, qmp); 238 _moveb(msg, &r, rmp); 239 BN_free(&q); 240 BN_free(&r); 241 } 242 243 void 244 mp_mdiv(const MINT *nmp, const MINT *dmp, MINT *qmp, MINT *rmp) 245 { 246 BN_CTX *c; 247 248 c = BN_CTX_new(); 249 if (c == NULL) 250 _bnerr("mdiv"); 251 _mdiv("mdiv", nmp, dmp, qmp, rmp, c); 252 BN_CTX_free(c); 253 } 254 255 /* 256 * Free memory associated with an MINT. 257 */ 258 static void 259 _mfree(const char *msg __unused, MINT *mp) 260 { 261 262 BN_clear(mp->bn); 263 BN_free(mp->bn); 264 free(mp); 265 } 266 267 void 268 mp_mfree(MINT *mp) 269 { 270 271 _mfree("mfree", mp); 272 } 273 274 /* 275 * Read an integer from standard input and stick the result in mp. 276 * The input is treated to be in base 10. This must be the silliest 277 * API in existence; why can't the program read in a string and call 278 * xtom()? (Or if base 10 is desires, perhaps dtom() could be 279 * exported.) 280 */ 281 void 282 mp_min(MINT *mp) 283 { 284 MINT *rmp; 285 char *line, *nline; 286 size_t linelen; 287 288 line = fgetln(stdin, &linelen); 289 if (line == NULL) 290 MPERR(("min")); 291 nline = malloc(linelen + 1); 292 if (nline == NULL) 293 MPERR(("min")); 294 memcpy(nline, line, linelen); 295 nline[linelen] = '\0'; 296 rmp = _dtom("min", nline); 297 _movem("min", rmp, mp); 298 _mfree("min", rmp); 299 free(nline); 300 } 301 302 /* 303 * Print the value of mp to standard output in base 10. See blurb 304 * above min() for why this is so useless. 305 */ 306 void 307 mp_mout(const MINT *mp) 308 { 309 char *s; 310 311 s = _mtod("mout", mp); 312 printf("%s", s); 313 free(s); 314 } 315 316 /* 317 * Set the value of tmp to the value of smp (i.e., tmp=smp). 318 */ 319 void 320 mp_move(const MINT *smp, MINT *tmp) 321 { 322 323 _movem("move", smp, tmp); 324 } 325 326 327 /* 328 * Internal routine to set the value of tmp to that of sbp. 329 */ 330 static void 331 _moveb(const char *msg, const BIGNUM *sbp, MINT *tmp) 332 { 333 334 BN_ERRCHECK(msg, BN_copy(tmp->bn, sbp)); 335 } 336 337 /* 338 * Internal routine to set the value of tmp to that of smp. 339 */ 340 static void 341 _movem(const char *msg, const MINT *smp, MINT *tmp) 342 { 343 344 BN_ERRCHECK(msg, BN_copy(tmp->bn, smp->bn)); 345 } 346 347 /* 348 * Compute the square root of nmp and put the result in xmp. The 349 * remainder goes in rmp. Should satisfy: rmp=nmp-(xmp*xmp). 350 * 351 * Note that the OpenSSL BIGNUM library does not have a square root 352 * function, so this had to be implemented by hand using Newton's 353 * recursive formula: 354 * 355 * x = (x + (n / x)) / 2 356 * 357 * where x is the square root of the positive number n. In the 358 * beginning, x should be a reasonable guess, but the value 1, 359 * although suboptimal, works, too; this is that is used below. 360 */ 361 void 362 mp_msqrt(const MINT *nmp, MINT *xmp, MINT *rmp) 363 { 364 BN_CTX *c; 365 MINT *tolerance; 366 MINT *ox, *x; 367 MINT *z1, *z2, *z3; 368 short i; 369 370 c = BN_CTX_new(); 371 if (c == NULL) 372 _bnerr("msqrt"); 373 tolerance = _itom("msqrt", 1); 374 x = _itom("msqrt", 1); 375 ox = _itom("msqrt", 0); 376 z1 = _itom("msqrt", 0); 377 z2 = _itom("msqrt", 0); 378 z3 = _itom("msqrt", 0); 379 do { 380 _movem("msqrt", x, ox); 381 _mdiv("msqrt", nmp, x, z1, z2, c); 382 _madd("msqrt", x, z1, z2); 383 _sdiv("msqrt", z2, 2, x, &i, c); 384 _msub("msqrt", ox, x, z3); 385 } while (_mcmpa("msqrt", z3, tolerance) == 1); 386 _movem("msqrt", x, xmp); 387 _mult("msqrt", x, x, z1, c); 388 _msub("msqrt", nmp, z1, z2); 389 _movem("msqrt", z2, rmp); 390 _mfree("msqrt", tolerance); 391 _mfree("msqrt", ox); 392 _mfree("msqrt", x); 393 _mfree("msqrt", z1); 394 _mfree("msqrt", z2); 395 _mfree("msqrt", z3); 396 BN_CTX_free(c); 397 } 398 399 /* 400 * Compute rmp=mp1-mp2. 401 */ 402 static void 403 _msub(const char *msg, const MINT *mp1, const MINT *mp2, MINT *rmp) 404 { 405 BIGNUM b; 406 407 BN_init(&b); 408 BN_ERRCHECK(msg, BN_sub(&b, mp1->bn, mp2->bn)); 409 _moveb(msg, &b, rmp); 410 BN_free(&b); 411 } 412 413 void 414 mp_msub(const MINT *mp1, const MINT *mp2, MINT *rmp) 415 { 416 417 _msub("msub", mp1, mp2, rmp); 418 } 419 420 /* 421 * Return a decimal representation of mp. Return value must be 422 * free()'d. 423 */ 424 static char * 425 _mtod(const char *msg, const MINT *mp) 426 { 427 char *s, *s2; 428 429 s = BN_bn2dec(mp->bn); 430 if (s == NULL) 431 _bnerr(msg); 432 asprintf(&s2, "%s", s); 433 if (s2 == NULL) 434 MPERR(("%s", msg)); 435 OPENSSL_free(s); 436 return (s2); 437 } 438 439 /* 440 * Return a hexadecimal representation of mp. Return value must be 441 * free()'d. 442 */ 443 static char * 444 _mtox(const char *msg, const MINT *mp) 445 { 446 char *p, *s, *s2; 447 int len; 448 449 s = BN_bn2hex(mp->bn); 450 if (s == NULL) 451 _bnerr(msg); 452 asprintf(&s2, "%s", s); 453 if (s2 == NULL) 454 MPERR(("%s", msg)); 455 OPENSSL_free(s); 456 457 /* 458 * This is a kludge for libgmp compatibility. The latter's 459 * implementation of this function returns lower-case letters, 460 * but BN_bn2hex returns upper-case. Some programs (e.g., 461 * newkey(1)) are sensitive to this. Although it's probably 462 * their fault, it's nice to be compatible. 463 */ 464 len = strlen(s2); 465 for (p = s2; p < s2 + len; p++) 466 *p = tolower(*p); 467 468 return (s2); 469 } 470 471 char * 472 mp_mtox(const MINT *mp) 473 { 474 475 return (_mtox("mtox", mp)); 476 } 477 478 /* 479 * Compute rmp=mp1*mp2. 480 */ 481 static void 482 _mult(const char *msg, const MINT *mp1, const MINT *mp2, MINT *rmp, BN_CTX *c) 483 { 484 BIGNUM b; 485 486 BN_init(&b); 487 BN_ERRCHECK(msg, BN_mul(&b, mp1->bn, mp2->bn, c)); 488 _moveb(msg, &b, rmp); 489 BN_free(&b); 490 } 491 492 void 493 mp_mult(const MINT *mp1, const MINT *mp2, MINT *rmp) 494 { 495 BN_CTX *c; 496 497 c = BN_CTX_new(); 498 if (c == NULL) 499 _bnerr("mult"); 500 _mult("mult", mp1, mp2, rmp, c); 501 BN_CTX_free(c); 502 } 503 504 /* 505 * Compute rmp=(bmp^emp)mod mmp. (Note that here and above rpow() '^' 506 * means 'raise to power', not 'bitwise XOR'.) 507 */ 508 void 509 mp_pow(const MINT *bmp, const MINT *emp, const MINT *mmp, MINT *rmp) 510 { 511 BIGNUM b; 512 BN_CTX *c; 513 514 c = BN_CTX_new(); 515 if (c == NULL) 516 _bnerr("pow"); 517 BN_init(&b); 518 BN_ERRCHECK("pow", BN_mod_exp(&b, bmp->bn, emp->bn, mmp->bn, c)); 519 _moveb("pow", &b, rmp); 520 BN_free(&b); 521 BN_CTX_free(c); 522 } 523 524 /* 525 * Compute rmp=bmp^e. (See note above pow().) 526 */ 527 void 528 mp_rpow(const MINT *bmp, short e, MINT *rmp) 529 { 530 MINT *emp; 531 BIGNUM b; 532 BN_CTX *c; 533 534 c = BN_CTX_new(); 535 if (c == NULL) 536 _bnerr("rpow"); 537 BN_init(&b); 538 emp = _itom("rpow", e); 539 BN_ERRCHECK("rpow", BN_exp(&b, bmp->bn, emp->bn, c)); 540 _moveb("rpow", &b, rmp); 541 _mfree("rpow", emp); 542 BN_free(&b); 543 BN_CTX_free(c); 544 } 545 546 /* 547 * Compute qmp=nmp/d and ro=nmp%d. 548 */ 549 static void 550 _sdiv(const char *msg, const MINT *nmp, short d, MINT *qmp, short *ro, 551 BN_CTX *c) 552 { 553 MINT *dmp, *rmp; 554 BIGNUM q, r; 555 char *s; 556 557 BN_init(&q); 558 BN_init(&r); 559 dmp = _itom(msg, d); 560 rmp = _itom(msg, 0); 561 BN_ERRCHECK(msg, BN_div(&q, &r, nmp->bn, dmp->bn, c)); 562 _moveb(msg, &q, qmp); 563 _moveb(msg, &r, rmp); 564 s = _mtox(msg, rmp); 565 errno = 0; 566 *ro = strtol(s, NULL, 16); 567 if (errno != 0) 568 MPERR(("%s underflow or overflow", msg)); 569 free(s); 570 _mfree(msg, dmp); 571 _mfree(msg, rmp); 572 BN_free(&r); 573 BN_free(&q); 574 } 575 576 void 577 mp_sdiv(const MINT *nmp, short d, MINT *qmp, short *ro) 578 { 579 BN_CTX *c; 580 581 c = BN_CTX_new(); 582 if (c == NULL) 583 _bnerr("sdiv"); 584 _sdiv("sdiv", nmp, d, qmp, ro, c); 585 BN_CTX_free(c); 586 } 587 588 /* 589 * Convert a hexadecimal string to an MINT. 590 */ 591 static MINT * 592 _xtom(const char *msg, const char *s) 593 { 594 MINT *mp; 595 596 mp = malloc(sizeof(*mp)); 597 if (mp == NULL) 598 MPERR(("%s", msg)); 599 mp->bn = BN_new(); 600 if (mp->bn == NULL) 601 _bnerr(msg); 602 BN_ERRCHECK(msg, BN_hex2bn(&mp->bn, s)); 603 return (mp); 604 } 605 606 MINT * 607 mp_xtom(const char *s) 608 { 609 610 return (_xtom("xtom", s)); 611 } 612