lib/libmp/mpasbn.c

   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);
 290         if (nline == NULL)
 291                 MPERR(("min"));
 292         strncpy(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 }