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