lib/libmp/mpasbn.c

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