2 * SPDX-License-Identifier: BSD-2-Clause-FreeBSD
4 * Copyright (c) 2001 Dima Dorfman.
7 * Redistribution and use in source and binary forms, with or without
8 * modification, are permitted provided that the following conditions
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.
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
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.
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
43 * Miscellaneous notes on this implementation:
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().
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
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
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.
74 #include <sys/cdefs.h>
75 __FBSDID("$FreeBSD$");
84 #include <openssl/crypto.h>
85 #include <openssl/err.h>
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); \
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 *,
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 *);
113 * Report an error from one of the BN_* functions using MPERRX.
116 _bnerr(const char *msg)
119 ERR_load_crypto_strings();
120 MPERRX(("%s: %s", msg, ERR_reason_error_string(ERR_get_error())));
124 * Convert a decimal string to an MINT.
127 _dtom(const char *msg, const char *s)
131 mp = malloc(sizeof(*mp));
137 BN_ERRCHECK(msg, BN_dec2bn(&mp->bn, s));
142 * Compute the greatest common divisor of mp1 and mp2; result goes in rmp.
145 mp_gcd(const MINT *mp1, const MINT *mp2, MINT *rmp)
154 BN_ERRCHECK("gcd", BN_gcd(&b, mp1->bn, mp2->bn, c));
155 _moveb("gcd", &b, rmp);
161 * Make an MINT out of a short integer. Return value must be mfree()'d.
164 _itom(const char *msg, short n)
169 asprintf(&s, "%x", n);
181 return (_itom("itom", n));
185 * Compute rmp=mp1+mp2.
188 _madd(const char *msg, const MINT *mp1, const MINT *mp2, MINT *rmp)
193 BN_ERRCHECK(msg, BN_add(&b, mp1->bn, mp2->bn));
194 _moveb(msg, &b, rmp);
199 mp_madd(const MINT *mp1, const MINT *mp2, MINT *rmp)
202 _madd("madd", mp1, mp2, rmp);
206 * Return -1, 0, or 1 if mp1<mp2, mp1==mp2, or mp1>mp2, respectivley.
209 mp_mcmp(const MINT *mp1, const MINT *mp2)
212 return (BN_cmp(mp1->bn, mp2->bn));
216 * Same as mcmp but compares absolute values.
219 _mcmpa(const char *msg __unused, const MINT *mp1, const MINT *mp2)
222 return (BN_ucmp(mp1->bn, mp2->bn));
226 * Compute qmp=nmp/dmp and rmp=nmp%dmp.
229 _mdiv(const char *msg, const MINT *nmp, const MINT *dmp, MINT *qmp, MINT *rmp,
236 BN_ERRCHECK(msg, BN_div(&q, &r, nmp->bn, dmp->bn, c));
237 _moveb(msg, &q, qmp);
238 _moveb(msg, &r, rmp);
244 mp_mdiv(const MINT *nmp, const MINT *dmp, MINT *qmp, MINT *rmp)
251 _mdiv("mdiv", nmp, dmp, qmp, rmp, c);
256 * Free memory associated with an MINT.
259 _mfree(const char *msg __unused, MINT *mp)
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
288 line = fgetln(stdin, &linelen);
291 nline = malloc(linelen + 1);
294 memcpy(nline, line, linelen);
295 nline[linelen] = '\0';
296 rmp = _dtom("min", nline);
297 _movem("min", rmp, mp);
303 * Print the value of mp to standard output in base 10. See blurb
304 * above min() for why this is so useless.
307 mp_mout(const MINT *mp)
311 s = _mtod("mout", mp);
317 * Set the value of tmp to the value of smp (i.e., tmp=smp).
320 mp_move(const MINT *smp, MINT *tmp)
323 _movem("move", smp, tmp);
328 * Internal routine to set the value of tmp to that of sbp.
331 _moveb(const char *msg, const BIGNUM *sbp, MINT *tmp)
334 BN_ERRCHECK(msg, BN_copy(tmp->bn, sbp));
338 * Internal routine to set the value of tmp to that of smp.
341 _movem(const char *msg, const MINT *smp, MINT *tmp)
344 BN_ERRCHECK(msg, BN_copy(tmp->bn, smp->bn));
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).
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
355 * x = (x + (n / x)) / 2
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.
362 mp_msqrt(const MINT *nmp, MINT *xmp, MINT *rmp)
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);
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);
400 * Compute rmp=mp1-mp2.
403 _msub(const char *msg, const MINT *mp1, const MINT *mp2, MINT *rmp)
408 BN_ERRCHECK(msg, BN_sub(&b, mp1->bn, mp2->bn));
409 _moveb(msg, &b, rmp);
414 mp_msub(const MINT *mp1, const MINT *mp2, MINT *rmp)
417 _msub("msub", mp1, mp2, rmp);
421 * Return a decimal representation of mp. Return value must be
425 _mtod(const char *msg, const MINT *mp)
429 s = BN_bn2dec(mp->bn);
432 asprintf(&s2, "%s", s);
440 * Return a hexadecimal representation of mp. Return value must be
444 _mtox(const char *msg, const MINT *mp)
449 s = BN_bn2hex(mp->bn);
452 asprintf(&s2, "%s", s);
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.
465 for (p = s2; p < s2 + len; p++)
472 mp_mtox(const MINT *mp)
475 return (_mtox("mtox", mp));
479 * Compute rmp=mp1*mp2.
482 _mult(const char *msg, const MINT *mp1, const MINT *mp2, MINT *rmp, BN_CTX *c)
487 BN_ERRCHECK(msg, BN_mul(&b, mp1->bn, mp2->bn, c));
488 _moveb(msg, &b, rmp);
493 mp_mult(const MINT *mp1, const MINT *mp2, MINT *rmp)
500 _mult("mult", mp1, mp2, rmp, c);
505 * Compute rmp=(bmp^emp)mod mmp. (Note that here and above rpow() '^'
506 * means 'raise to power', not 'bitwise XOR'.)
509 mp_pow(const MINT *bmp, const MINT *emp, const MINT *mmp, MINT *rmp)
518 BN_ERRCHECK("pow", BN_mod_exp(&b, bmp->bn, emp->bn, mmp->bn, c));
519 _moveb("pow", &b, rmp);
525 * Compute rmp=bmp^e. (See note above pow().)
528 mp_rpow(const MINT *bmp, short e, MINT *rmp)
538 emp = _itom("rpow", e);
539 BN_ERRCHECK("rpow", BN_exp(&b, bmp->bn, emp->bn, c));
540 _moveb("rpow", &b, rmp);
547 * Compute qmp=nmp/d and ro=nmp%d.
550 _sdiv(const char *msg, const MINT *nmp, short d, MINT *qmp, short *ro,
561 BN_ERRCHECK(msg, BN_div(&q, &r, nmp->bn, dmp->bn, c));
562 _moveb(msg, &q, qmp);
563 _moveb(msg, &r, rmp);
566 *ro = strtol(s, NULL, 16);
568 MPERR(("%s underflow or overflow", msg));
577 mp_sdiv(const MINT *nmp, short d, MINT *qmp, short *ro)
584 _sdiv("sdiv", nmp, d, qmp, ro, c);
589 * Convert a hexadecimal string to an MINT.
592 _xtom(const char *msg, const char *s)
596 mp = malloc(sizeof(*mp));
602 BN_ERRCHECK(msg, BN_hex2bn(&mp->bn, s));
607 mp_xtom(const char *s)
610 return (_xtom("xtom", s));