xref: /freebsd/lib/libmp/mpasbn.c (revision f0cfa1b168014f56c02b83e5f28412cc5f78d117)
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