lib/libmp/mpasbn.c

/*
 * Copyright (c) 2001 Dima Dorfman.
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 */

/*
 * This is the traditional Berkeley MP library implemented in terms of
 * the OpenSSL BIGNUM library.  It was written to replace libgmp, and
 * is meant to be as compatible with the latter as feasible.
 *
 * There seems to be a lack of documentation for the Berkeley MP
 * interface.  All I could find was libgmp documentation (which didn't
 * talk about the semantics of the functions) and an old SunOS 4.1
 * manual page from 1989.  The latter wasn't very detailed, either,
 * but at least described what the function's arguments were.  In
 * general the interface seems to be archaic, somewhat poorly
 * designed, and poorly, if at all, documented.  It is considered
 * harmful.
 *
 * Miscellaneous notes on this implementation:
 *
 *  - The SunOS manual page mentioned above indicates that if an error
 *  occurs, the library should "produce messages and core images."
 *  Given that most of the functions don't have return values (and
 *  thus no sane way of alerting the caller to an error), this seems
 *  reasonable.  The MPERR and MPERRX macros call warn and warnx,
 *  respectively, then abort().
 *
 *  - All the functions which take an argument to be "filled in"
 *  assume that the argument has been initialized by one of the *tom()
 *  routines before being passed to it.  I never saw this documented
 *  anywhere, but this seems to be consistent with the way this
 *  library is used.
 *
 *  - msqrt() is the only routine which had to be implemented which
 *  doesn't have a close counterpart in the OpenSSL BIGNUM library.
 *  It was implemented by hand using Newton's recursive formula.
 *  Doing it this way, although more error-prone, has the positive
 *  sideaffect of testing a lot of other functions; if msqrt()
 *  produces the correct results, most of the other routines will as
 *  well.
 *
 *  - Internal-use-only routines (i.e., those defined here statically
 *  and not in mp.h) have an underscore prepended to their name (this
 *  is more for aesthetical reasons than technical).  All such
 *  routines take an extra argument, 'msg', that denotes what they
 *  should call themselves in an error message.  This is so a user
 *  doesn't get an error message from a function they didn't call.
 */

#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");

#include <ctype.h>
#include <err.h>
#include <errno.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <openssl/crypto.h>
#include <openssl/err.h>

#include "mp.h"

#define MPERR(s)        do { warn s; abort(); } while (0)
#define MPERRX(s)       do { warnx s; abort(); } while (0)
#define BN_ERRCHECK(msg, expr) do {             \
        if (!(expr)) _bnerr(msg);               \
} while (0)

static void _bnerr(const char *);
static MINT *_dtom(const char *, const char *);
static MINT *_itom(const char *, short);
static void _madd(const char *, const MINT *, const MINT *, MINT *);
static int _mcmpa(const char *, const MINT *, const MINT *);
static void _mdiv(const char *, const MINT *, const MINT *, MINT *, MINT *,
                BN_CTX *);
static void _mfree(const char *, MINT *);
static void _moveb(const char *, const BIGNUM *, MINT *);
static void _movem(const char *, const MINT *, MINT *);
static void _msub(const char *, const MINT *, const MINT *, MINT *);
static char *_mtod(const char *, const MINT *);
static char *_mtox(const char *, const MINT *);
static void _mult(const char *, const MINT *, const MINT *, MINT *, BN_CTX *);
static void _sdiv(const char *, const MINT *, short, MINT *, short *, BN_CTX *);
static MINT *_xtom(const char *, const char *);

/*
 * Report an error from one of the BN_* functions using MPERRX.
 */
static void
_bnerr(const char *msg)
{

        ERR_load_crypto_strings();
        MPERRX(("%s: %s", msg, ERR_reason_error_string(ERR_get_error())));
}

/*
 * Convert a decimal string to an MINT.
 */
static MINT *
_dtom(const char *msg, const char *s)
{
        MINT *mp;

        mp = malloc(sizeof(*mp));
        if (mp == NULL)
                MPERR(("%s", msg));
        mp->bn = BN_new();
        if (mp->bn == NULL)
                _bnerr(msg);
        BN_ERRCHECK(msg, BN_dec2bn(&mp->bn, s));
        return (mp);
}

/*
 * Compute the greatest common divisor of mp1 and mp2; result goes in rmp.
 */
void
mp_gcd(const MINT *mp1, const MINT *mp2, MINT *rmp)
{
        BIGNUM b;
        BN_CTX *c;

        c = BN_CTX_new();
        if (c == NULL)
                _bnerr("gcd");
        BN_init(&b);
        BN_ERRCHECK("gcd", BN_gcd(&b, mp1->bn, mp2->bn, c));
        _moveb("gcd", &b, rmp);
        BN_free(&b);
        BN_CTX_free(c);
}

/*
 * Make an MINT out of a short integer.  Return value must be mfree()'d.
 */
static MINT *
_itom(const char *msg, short n)
{
        MINT *mp;
        char *s;

        asprintf(&s, "%x", n);
        if (s == NULL)
                MPERR(("%s", msg));
        mp = _xtom(msg, s);
        free(s);
        return (mp);
}

MINT *
mp_itom(short n)
{

        return (_itom("itom", n));
}

/*
 * Compute rmp=mp1+mp2.
 */
static void
_madd(const char *msg, const MINT *mp1, const MINT *mp2, MINT *rmp)
{
        BIGNUM b;

        BN_init(&b);
        BN_ERRCHECK(msg, BN_add(&b, mp1->bn, mp2->bn));
        _moveb(msg, &b, rmp);
        BN_free(&b);
}

void
mp_madd(const MINT *mp1, const MINT *mp2, MINT *rmp)
{

        _madd("madd", mp1, mp2, rmp);
}

/*
 * Return -1, 0, or 1 if mp1<mp2, mp1==mp2, or mp1>mp2, respectivley.
 */
int
mp_mcmp(const MINT *mp1, const MINT *mp2)
{

        return (BN_cmp(mp1->bn, mp2->bn));
}

/*
 * Same as mcmp but compares absolute values.
 */
static int
_mcmpa(const char *msg __unused, const MINT *mp1, const MINT *mp2)
{

        return (BN_ucmp(mp1->bn, mp2->bn));
}

/*
 * Compute qmp=nmp/dmp and rmp=nmp%dmp.
 */
static void
_mdiv(const char *msg, const MINT *nmp, const MINT *dmp, MINT *qmp, MINT *rmp,
    BN_CTX *c)
{
        BIGNUM q, r;

        BN_init(&r);
        BN_init(&q);
        BN_ERRCHECK(msg, BN_div(&q, &r, nmp->bn, dmp->bn, c));
        _moveb(msg, &q, qmp);
        _moveb(msg, &r, rmp);
        BN_free(&q);
        BN_free(&r);
}

void
mp_mdiv(const MINT *nmp, const MINT *dmp, MINT *qmp, MINT *rmp)
{
        BN_CTX *c;

        c = BN_CTX_new();
        if (c == NULL)
                _bnerr("mdiv");
        _mdiv("mdiv", nmp, dmp, qmp, rmp, c);
        BN_CTX_free(c);
}

/*
 * Free memory associated with an MINT.
 */
static void
_mfree(const char *msg __unused, MINT *mp)
{

        BN_clear(mp->bn);
        BN_free(mp->bn);
        free(mp);
}

void
mp_mfree(MINT *mp)
{

        _mfree("mfree", mp);
}

/*
 * Read an integer from standard input and stick the result in mp.
 * The input is treated to be in base 10.  This must be the silliest
 * API in existence; why can't the program read in a string and call
 * xtom()?  (Or if base 10 is desires, perhaps dtom() could be
 * exported.)
 */
void
mp_min(MINT *mp)
{
        MINT *rmp;
        char *line, *nline;
        size_t linelen;

        line = fgetln(stdin, &linelen);
        if (line == NULL)
                MPERR(("min"));
        nline = malloc(linelen);
        if (nline == NULL)
                MPERR(("min"));
        strncpy(nline, line, linelen);
        nline[linelen] = '\0';
        rmp = _dtom("min", nline);
        _movem("min", rmp, mp);
        _mfree("min", rmp);
        free(nline);
}

/*
 * Print the value of mp to standard output in base 10.  See blurb
 * above min() for why this is so useless.
 */
void
mp_mout(const MINT *mp)
{
        char *s;

        s = _mtod("mout", mp);
        printf("%s", s);
        free(s);
}

/*
 * Set the value of tmp to the value of smp (i.e., tmp=smp).
 */
void
mp_move(const MINT *smp, MINT *tmp)
{

        _movem("move", smp, tmp);
}


/*
 * Internal routine to set the value of tmp to that of sbp.
 */
static void
_moveb(const char *msg, const BIGNUM *sbp, MINT *tmp)
{

        BN_ERRCHECK(msg, BN_copy(tmp->bn, sbp));
}

/*
 * Internal routine to set the value of tmp to that of smp.
 */
static void
_movem(const char *msg, const MINT *smp, MINT *tmp)
{

        BN_ERRCHECK(msg, BN_copy(tmp->bn, smp->bn));
}

/*
 * Compute the square root of nmp and put the result in xmp.  The
 * remainder goes in rmp.  Should satisfy: rmp=nmp-(xmp*xmp).
 *
 * Note that the OpenSSL BIGNUM library does not have a square root
 * function, so this had to be implemented by hand using Newton's
 * recursive formula:
 *
 *              x = (x + (n / x)) / 2
 *
 * where x is the square root of the positive number n.  In the
 * beginning, x should be a reasonable guess, but the value 1,
 * although suboptimal, works, too; this is that is used below.
 */
void
mp_msqrt(const MINT *nmp, MINT *xmp, MINT *rmp)
{
        BN_CTX *c;
        MINT *tolerance;
        MINT *ox, *x;
        MINT *z1, *z2, *z3;
        short i;

        c = BN_CTX_new();
        if (c == NULL)
                _bnerr("msqrt");
        tolerance = _itom("msqrt", 1);
        x = _itom("msqrt", 1);
        ox = _itom("msqrt", 0);
        z1 = _itom("msqrt", 0);
        z2 = _itom("msqrt", 0);
        z3 = _itom("msqrt", 0);
        do {
                _movem("msqrt", x, ox);
                _mdiv("msqrt", nmp, x, z1, z2, c);
                _madd("msqrt", x, z1, z2);
                _sdiv("msqrt", z2, 2, x, &i, c);
                _msub("msqrt", ox, x, z3);
        } while (_mcmpa("msqrt", z3, tolerance) == 1);
        _movem("msqrt", x, xmp);
        _mult("msqrt", x, x, z1, c);
        _msub("msqrt", nmp, z1, z2);
        _movem("msqrt", z2, rmp);
        _mfree("msqrt", tolerance);
        _mfree("msqrt", ox);
        _mfree("msqrt", x);
        _mfree("msqrt", z1);
        _mfree("msqrt", z2);
        _mfree("msqrt", z3);
        BN_CTX_free(c);
}

/*
 * Compute rmp=mp1-mp2.
 */
static void
_msub(const char *msg, const MINT *mp1, const MINT *mp2, MINT *rmp)
{
        BIGNUM b;

        BN_init(&b);
        BN_ERRCHECK(msg, BN_sub(&b, mp1->bn, mp2->bn));
        _moveb(msg, &b, rmp);
        BN_free(&b);
}

void
mp_msub(const MINT *mp1, const MINT *mp2, MINT *rmp)
{

        _msub("msub", mp1, mp2, rmp);
}

/*
 * Return a decimal representation of mp.  Return value must be
 * free()'d.
 */
static char *
_mtod(const char *msg, const MINT *mp)
{
        char *s, *s2;

        s = BN_bn2dec(mp->bn);
        if (s == NULL)
                _bnerr(msg);
        asprintf(&s2, "%s", s);
        if (s2 == NULL)
                MPERR(("%s", msg));
        OPENSSL_free(s);
        return (s2);
}

/*
 * Return a hexadecimal representation of mp.  Return value must be
 * free()'d.
 */
static char *
_mtox(const char *msg, const MINT *mp)
{
        char *p, *s, *s2;
        int len;

        s = BN_bn2hex(mp->bn);
        if (s == NULL)
                _bnerr(msg);
        asprintf(&s2, "%s", s);
        if (s2 == NULL)
                MPERR(("%s", msg));
        OPENSSL_free(s);

        /*
         * This is a kludge for libgmp compatibility.  The latter's
         * implementation of this function returns lower-case letters,
         * but BN_bn2hex returns upper-case.  Some programs (e.g.,
         * newkey(1)) are sensitive to this.  Although it's probably
         * their fault, it's nice to be compatible.
         */
        len = strlen(s2);
        for (p = s2; p < s2 + len; p++)
                *p = tolower(*p);

        return (s2);
}

char *
mp_mtox(const MINT *mp)
{

        return (_mtox("mtox", mp));
}

/*
 * Compute rmp=mp1*mp2.
 */
static void
_mult(const char *msg, const MINT *mp1, const MINT *mp2, MINT *rmp, BN_CTX *c)
{
        BIGNUM b;

        BN_init(&b);
        BN_ERRCHECK(msg, BN_mul(&b, mp1->bn, mp2->bn, c));
        _moveb(msg, &b, rmp);
        BN_free(&b);
}

void
mp_mult(const MINT *mp1, const MINT *mp2, MINT *rmp)
{
        BN_CTX *c;

        c = BN_CTX_new();
        if (c == NULL)
                _bnerr("mult");
        _mult("mult", mp1, mp2, rmp, c);
        BN_CTX_free(c);
}

/*
 * Compute rmp=(bmp^emp)mod mmp.  (Note that here and above rpow() '^'
 * means 'raise to power', not 'bitwise XOR'.)
 */
void
mp_pow(const MINT *bmp, const MINT *emp, const MINT *mmp, MINT *rmp)
{
        BIGNUM b;
        BN_CTX *c;

        c = BN_CTX_new();
        if (c == NULL)
                _bnerr("pow");
        BN_init(&b);
        BN_ERRCHECK("pow", BN_mod_exp(&b, bmp->bn, emp->bn, mmp->bn, c));
        _moveb("pow", &b, rmp);
        BN_free(&b);
        BN_CTX_free(c);
}

/*
 * Compute rmp=bmp^e.  (See note above pow().)
 */
void
mp_rpow(const MINT *bmp, short e, MINT *rmp)
{
        MINT *emp;
        BIGNUM b;
        BN_CTX *c;

        c = BN_CTX_new();
        if (c == NULL)
                _bnerr("rpow");
        BN_init(&b);
        emp = _itom("rpow", e);
        BN_ERRCHECK("rpow", BN_exp(&b, bmp->bn, emp->bn, c));
        _moveb("rpow", &b, rmp);
        _mfree("rpow", emp);
        BN_free(&b);
        BN_CTX_free(c);
}

/*
 * Compute qmp=nmp/d and ro=nmp%d.
 */
static void
_sdiv(const char *msg, const MINT *nmp, short d, MINT *qmp, short *ro,
    BN_CTX *c)
{
        MINT *dmp, *rmp;
        BIGNUM q, r;
        char *s;

        BN_init(&q);
        BN_init(&r);
        dmp = _itom(msg, d);
        rmp = _itom(msg, 0);
        BN_ERRCHECK(msg, BN_div(&q, &r, nmp->bn, dmp->bn, c));
        _moveb(msg, &q, qmp);
        _moveb(msg, &r, rmp);
        s = _mtox(msg, rmp);
        errno = 0;
        *ro = strtol(s, NULL, 16);
        if (errno != 0)
                MPERR(("%s underflow or overflow", msg));
        free(s);
        _mfree(msg, dmp);
        _mfree(msg, rmp);
        BN_free(&r);
        BN_free(&q);
}

void
mp_sdiv(const MINT *nmp, short d, MINT *qmp, short *ro)
{
        BN_CTX *c;

        c = BN_CTX_new();
        if (c == NULL)
                _bnerr("sdiv");
        _sdiv("sdiv", nmp, d, qmp, ro, c);
        BN_CTX_free(c);
}

/*
 * Convert a hexadecimal string to an MINT.
 */
static MINT *
_xtom(const char *msg, const char *s)
{
        MINT *mp;

        mp = malloc(sizeof(*mp));
        if (mp == NULL)
                MPERR(("%s", msg));
        mp->bn = BN_new();
        if (mp->bn == NULL)
                _bnerr(msg);
        BN_ERRCHECK(msg, BN_hex2bn(&mp->bn, s));
        return (mp);
}

MINT *
mp_xtom(const char *s)
{

        return (_xtom("xtom", s));
}
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	}