From: mdw Date: Wed, 17 Nov 1999 18:02:17 +0000 (+0000) Subject: New multiprecision integer arithmetic suite. X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/commitdiff_plain/d3409d5ecf2492cff862616de72a580d1a8e8dc0 New multiprecision integer arithmetic suite. --- diff --git a/mp-arith.c b/mp-arith.c new file mode 100644 index 0000000..d8381ee --- /dev/null +++ b/mp-arith.c @@ -0,0 +1,422 @@ +/* -*-c-*- + * + * $Id: mp-arith.c,v 1.1 1999/11/17 18:02:16 mdw Exp $ + * + * Basic arithmetic on multiprecision integers + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: mp-arith.c,v $ + * Revision 1.1 1999/11/17 18:02:16 mdw + * New multiprecision integer arithmetic suite. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include "mp.h" + +/*----- Main code ---------------------------------------------------------*/ + +/* --- @mp_2c@ --- * + * + * Arguments: @mp *a@ = source + * + * Returns: Result, @a@ converted to two's complement notation. + */ + +mp *mp_2c(mp *d, mp *a) +{ + if (!(a->f & MP_NEG)) + return (MP_COPY(a)); + + MP_MODIFY(d, MP_LEN(a)); + mpx_2c(d->v, d->vl, a->v, a->vl); + d->f = a->f & MP_BURN; + MP_SHRINK(d); + return (d); +} + +/* --- @mp_sm@ --- * + * + * Arguments: @mp *d@ = destination + * @mp *a@ = source + * + * Returns: Result, @a@ converted to the native signed-magnitude + * notation. + */ + +mp *mp_sm(mp *d, mp *a) +{ + if (!MP_LEN(a) || a->vl[-1] < MPW_MAX / 2) + return (MP_COPY(a)); + + MP_MODIFY(d, MP_LEN(a)); + mpx_2c(d->v, d->vl, a->v, a->vl); + d->f = (a->f & (MP_BURN | MP_NEG)) ^ MP_NEG; + MP_SHRINK(d); + return (d); +} + +/* --- @mp_lsl@ --- * + * + * Arguments: @mp *d@ = destination + * @const mp *a@ = source + * @size_t n@ = number of bits to move + * + * Returns: Result, @a@ shifted left by @n@. + */ + +mp *mp_lsl(mp *d, const mp *a, size_t n) +{ + MP_MODIFY(d, MP_LEN(a) + (n + MPW_BITS - 1) / MPW_BITS); + mpx_lsl(d->v, d->vl, a->v, a->vl, n); + d->f = a->f & (MP_NEG | MP_BURN); + MP_SHRINK(d); + return (d); +} + +/* --- @mp_lsr@ --- * + * + * Arguments: @mp *d@ = destination + * @const mp *a@ = source + * @size_t n@ = number of bits to move + * + * Returns: Result, @a@ shifted left by @n@. + */ + +mp *mp_lsr(mp *d, const mp *a, size_t n) +{ + MP_MODIFY(d, MP_LEN(a)); + mpx_lsr(d->v, d->vl, a->v, a->vl, n); + d->f = a->f & (MP_NEG | MP_BURN); + MP_SHRINK(d); + return (d); +} + +/* --- @mp_cmp@ --- * + * + * Arguments: @const mp *a, *b@ = two numbers + * + * Returns: Less than, equal to or greater than zero, according to + * whether @a@ is less than, equal to or greater than @b@. + */ + +int mp_cmp(const mp *a, const mp *b) +{ + if (!((a->f ^ b->f) & MP_NEG)) + return (mpx_ucmp(a->v, a->vl, b->v, b->vl)); + else if (a->f & MP_NEG) + return (-1); + else + return (+1); +} + +/* --- @mp_add@ --- * + * + * Arguments: @mp *d@ = destination + * @const mp *a, *b@ = sources + * + * Returns: Result, @a@ added to @b@. + */ + +mp *mp_add(mp *d, const mp *a, const mp *b) +{ + MP_MODIFY(d, (MP_LEN(a) > MP_LEN(b) ? MP_LEN(a) : MP_LEN(b)) + 1); + if (!((a->f ^ b->f) & MP_NEG)) + mpx_uadd(d->v, d->vl, a->v, a->vl, b->v, b->vl); + else { + if (MPX_UCMP(a->v, a->vl, <, b->v, b->vl)) { + const mp *t = a; a = b; b = t; + } + mpx_usub(d->v, d->vl, a->v, a->vl, b->v, b->vl); + } + d->f = ((a->f | b->f) & MP_BURN) | (a->f & MP_NEG); + MP_SHRINK(d); + return (d); +} + +/* --- @mp_sub@ --- * + * + * Arguments: @mp *d@ = destination + * @const mp *a, *b@ = sources + * + * Returns: Result, @b@ subtracted from @a@. + */ + +mp *mp_sub(mp *d, const mp *a, const mp *b) +{ + unsigned sgn = 0; + MP_MODIFY(d, (MP_LEN(a) > MP_LEN(b) ? MP_LEN(a) : MP_LEN(b)) + 1); + if ((a->f ^ b->f) & MP_NEG) + mpx_uadd(d->v, d->vl, a->v, a->vl, b->v, b->vl); + else { + if (MPX_UCMP(a->v, a->vl, <, b->v, b->vl)) { + const mp *t = a; a = b; b = t; + sgn = MP_NEG; + } + mpx_usub(d->v, d->vl, a->v, a->vl, b->v, b->vl); + } + d->f = ((a->f | b->f) & MP_BURN) | ((a->f ^ sgn) & MP_NEG); + MP_SHRINK(d); + return (d); +} + +/* --- @mp_mul@ --- * + * + * Arguments: @mp *d@ = destination + * @const mp *a, *b@ = sources + * + * Returns: Result, @a@ multiplied by @b@. + */ + +mp *mp_mul(mp *d, const mp *a, const mp *b) +{ + if (d == a || d == b) + d = MP_NEW; + MP_MODIFY(d, MP_LEN(a) + MP_LEN(b)); + mpx_umul(d->v, d->vl, a->v, a->vl, b->v, b->vl); + d->f = ((a->f | b->f) & MP_BURN) | ((a->f ^ b->f) & MP_NEG); + MP_SHRINK(d); + return (d); +} + +/* --- @mp_sqr@ --- * + * + * Arguments: @mp *d@ = destination + * @const mp *a@ = source + * + * Returns: Result, @a@ squared. + */ + +mp *mp_sqr(mp *d, const mp *a) +{ + if (d == a) + d = MP_NEW; + MP_MODIFY(d, 2 * MP_LEN(a)); + mpx_usqr(d->v, d->vl, a->v, a->vl); + d->f = a->f & MP_BURN; + MP_SHRINK(d); + return (d); +} + +/* --- @mp_div@ --- * + * + * Arguments: @mp **qq, **rr@ = destination, quotient and remainder + * @const mp *a, *b@ = sources + * + * Use: Calculates the quotient and remainder when @a@ is divided by + * @b@. The destinations @*qq@ and @*rr@ must be distinct. + * Either of @qq@ or @rr@ may be null to indicate that the + * result is irrelevant. (Discarding both results is silly.) + * There is a performance advantage if @a == *rr@. + * + * The behaviour when @a@ and @b@ have the same sign is + * straightforward. When the signs differ, this implementation + * chooses @r@ to have the same sign as @b@, rather than the + * more normal choice that the remainder has the same sign as + * the dividend. This makes modular arithmetic a little more + * straightforward. + */ + +void mp_div(mp **qq, mp **rr, const mp *a, const mp *b) + { + mp *r = rr ? *rr : MP_NEW; + mp *q = qq ? *qq : MP_NEW; + mpw *sv, *svl; + + /* --- Set up some temporary workspace --- */ + + { + size_t rq = MP_LEN(b) + 1; + sv = MP_ALLOC(rq); + svl = sv + rq; + } + + /* --- Set the remainder up right --- * + * + * Just in case the divisor is larger, be able to cope with this. It's not + * important in @mpx_udiv@, but it is here because of the sign correction. + */ + + { + size_t rq = MP_LEN(a) + 2; + if (MP_LEN(b) > rq) + rq = MP_LEN(b); + + if (r == a) { + MP_SPLIT(r); + MP_ENSURE(r, MP_LEN(r) + 2); + } else { + if (r == b) + r = MP_NEW; + MP_MODIFY(r, MP_LEN(a) + 2); + memcpy(r->v, a->v, MPWS(MP_LEN(a))); + memset(r->v + MP_LEN(a), 0, MPWS(2)); + } + } + + /* --- Fix up the quotient too --- */ + + if (q == a || q == b) + q = MP_NEW; + MP_MODIFY(q, MP_LEN(a)); + + /* --- Perform the calculation --- */ + + mpx_udiv(q->v, q->vl, r->v, r->vl, b->v, b->vl, sv, svl); + + /* --- Sort out the sign of the results --- * + * + * If the signs of the arguments differ, and the remainder is nonzero, I + * must add one to the absolute value of the quotient and subtract the + * remainder from @b@. + */ + + q->f = ((a->f | b->f) & MP_BURN) | ((a->f ^ b->f) & MP_NEG); + if (q->f & MP_NEG) { + mpw *v = r->v; + while (v < r->vl) { + if (*v) { + MPX_UADDN(q->v, q->vl, 1); + mpx_usub(r->v, r->vl, b->v, b->vl, r->v, r->vl); + break; + } + } + } + + r->f = ((a->f | b->f) & MP_BURN) | (b->f & MP_NEG); + + /* --- Store the return values --- */ + + if (!qq) + MP_DROP(q); + else { + MP_SHRINK(q); + *qq = q; + } + + if (!rr) + MP_DROP(r); + else { + MP_SHRINK(r); + *rr = r; + } + + MP_FREE(sv); +} + +/*----- Test rig ----------------------------------------------------------*/ + +#ifdef TEST_RIG + +static int verify(const char *op, mp *expect, mp *result, mp *a, mp *b) +{ + if (MP_CMP(expect, !=, result)) { + fprintf(stderr, "\n*** %s failed", op); + fputs("\n*** a = ", stderr); mp_writefile(a, stderr, 10); + fputs("\n*** b = ", stderr); mp_writefile(b, stderr, 10); + fputs("\n*** result = ", stderr); mp_writefile(result, stderr, 10); + fputs("\n*** expect = ", stderr); mp_writefile(expect, stderr, 10); + fputc('\n', stderr); + return (0); + } + return (1); +} + +#define RIG(name, op) \ + static int t ## name(dstr *v) \ + { \ + mp *a = *(mp **)v[0].buf; \ + mpw n = *(int *)v[1].buf; \ + mp b; \ + mp *r = *(mp **)v[2].buf; \ + mp *c = op(MP_NEW, a, n); \ + int ok; \ + mp_build(&b, &n, &n + 1); \ + ok = verify(#name, r, c, a, &b); \ + mp_drop(a); mp_drop(c); mp_drop(r); \ + return (ok); \ + } + +RIG(lsl, mp_lsl) +RIG(lsr, mp_lsr) + +#undef RIG + +#define RIG(name, op) \ + static int t ## name(dstr *v) \ + { \ + mp *a = *(mp **)v[0].buf; \ + mp *b = *(mp **)v[1].buf; \ + mp *r = *(mp **)v[2].buf; \ + mp *c = op(MP_NEW, a, b); \ + int ok = verify(#name, r, c, a, b); \ + mp_drop(a); mp_drop(b); mp_drop(c); mp_drop(r); \ + return (ok); \ + } + +RIG(add, mp_add) +RIG(sub, mp_sub) +RIG(mul, mp_mul) + +#undef RIG + +static int tdiv(dstr *v) +{ + mp *a = *(mp **)v[0].buf; + mp *b = *(mp **)v[1].buf; + mp *q = *(mp **)v[2].buf; + mp *r = *(mp **)v[3].buf; + mp *c = MP_NEW, *d = MP_NEW; + int ok = 1; + mp_div(&c, &d, a, b); + ok &= verify("div(quotient)", q, c, a, b); + ok &= verify("div(remainder)", r, d, a, b); + mp_drop(a); mp_drop(b); mp_drop(c); mp_drop(d); mp_drop(r); mp_drop(q); + return (ok); +} + +static test_chunk tests[] = { + { "lsl", tlsl, { &type_mp, &type_mp, &type_mp, 0 } }, + { "lsr", tlsr, { &type_mp, &type_mp, &type_mp, 0 } }, + { "add", tadd, { &type_mp, &type_mp, &type_mp, 0 } }, + { "sub", tsub, { &type_mp, &type_mp, &type_mp, 0 } }, + { "mul", tmul, { &type_mp, &type_mp, &type_mp, 0 } }, + { "div", tdiv, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } }, + { 0, 0, { 0 } }, +}; + +int main(int argc, char *argv[]) +{ + sub_init(); + test_run(argc, argv, tests, SRCDIR "/tests/mp"); + return (0); +} + +#endif + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/mp-const.c b/mp-const.c new file mode 100644 index 0000000..2119470 --- /dev/null +++ b/mp-const.c @@ -0,0 +1,57 @@ +/* -*-c-*- + * + * $Id: mp-const.c,v 1.1 1999/11/17 18:02:16 mdw Exp $ + * + * Useful multiprecision constants + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: mp-const.c,v $ + * Revision 1.1 1999/11/17 18:02:16 mdw + * New multiprecision integer arithmetic suite. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include "mp.h" + +/*----- Global variables --------------------------------------------------*/ + +static mpw mpw_const[] = { 1, 2, 3, 4, 5, 10 }; + +mp mp_const[] = { + { &mpw_const[0], &mpw_const[0], 0, MP_CONST, 0 }, + { &mpw_const[0], &mpw_const[0] + 1, 1, MP_CONST, 0 }, + { &mpw_const[1], &mpw_const[1] + 1, 1, MP_CONST, 0 }, + { &mpw_const[2], &mpw_const[2] + 1, 1, MP_CONST, 0 }, + { &mpw_const[3], &mpw_const[3] + 1, 1, MP_CONST, 0 }, + { &mpw_const[4], &mpw_const[4] + 1, 1, MP_CONST, 0 }, + { &mpw_const[5], &mpw_const[5] + 1, 1, MP_CONST, 0 }, + { &mpw_const[0], &mpw_const[0] + 1, 1, MP_CONST | MP_NEG, 0 }, +}; + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/mp-gcd.c b/mp-gcd.c new file mode 100644 index 0000000..8298698 --- /dev/null +++ b/mp-gcd.c @@ -0,0 +1,316 @@ +/* -*-c-*- + * + * $Id: mp-gcd.c,v 1.1 1999/11/17 18:02:16 mdw Exp $ + * + * Extended GCD calculation + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: mp-gcd.c,v $ + * Revision 1.1 1999/11/17 18:02:16 mdw + * New multiprecision integer arithmetic suite. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include "mp.h" + +/*----- Main code ---------------------------------------------------------*/ + +/* --- @mp_gcd@ --- * + * + * Arguments: @mp **gcd, **xx, **yy@ = where to write the results + * @mp *a, *b@ = sources (must be nonzero) + * + * Returns: --- + * + * Use: Calculates @gcd(a, b)@, and two numbers @x@ and @y@ such that + * @ax + by = gcd(a, b)@. This is useful for computing modular + * inverses. Neither @a@ nor @b@ may be zero. Note that, + * unlike @mp_div@ for example, it is not possible to specify + * explicit destinations -- new MPs are always allocated. + */ + +void mp_gcd(mp **gcd, mp **xx, mp **yy, mp *a, mp *b) +{ + mp *X = MP_ONE, *Y = MP_ZERO; + mp *x = MP_ZERO, *y = MP_ONE; + mp *u, *v; + size_t shift = 0; + int ext = xx || yy; + int swap = 0; + + /* --- Ensure that @a@ is larger than @b@ --- */ + + if (MP_CMP(a, <, b)) { + { mp *t = a; a = b; b = t; } + swap = 1; + } + + /* --- Take a reference to the arguments --- */ + + a = MP_COPY(a); + b = MP_COPY(b); + + /* --- Make sure @a@ and @b@ are not both even --- */ + + if (((a->v[0] | b->v[0]) & 1) == 0) { + mpscan asc, bsc; + + /* --- Break off my copies --- */ + + MP_SPLIT(a); + MP_SPLIT(b); + MP_SCAN(&asc, a); + MP_SCAN(&bsc, b); + + /* --- Start scanning --- */ + + for (;;) { + if (!MP_STEP(&asc) || !MP_STEP(&bsc)) + assert(((void)"zero argument passed to mp_gcd", 0)); + if (MP_BIT(&asc) || MP_BIT(&bsc)) + break; + shift++; + } + + /* --- Shift @a@ and @b@ down --- */ + + a = mp_lsr(a, a, shift); + b = mp_lsr(b, b, shift); + } + + /* --- Set up @u@ and @v@ --- */ + + u = MP_COPY(a); + v = MP_COPY(b); + + /* --- Start the main loop --- */ + + for (;;) { + + /* --- While @u@ is even --- */ + + { + mpscan sc, xsc, ysc; + size_t n = 0, nn = 0; + + MP_SCAN(&sc, u); + MP_SCAN(&xsc, X); MP_SCAN(&ysc, Y); + for (;;) { + MP_STEP(&sc); + MP_STEP(&xsc); MP_STEP(&ysc); + if (MP_BIT(&sc)) + break; + if (ext && (MP_BIT(&xsc) | MP_BIT(&ysc))) { + if (n) { + X = mp_lsr(X, X, n); + Y = mp_lsr(Y, Y, n); + n = 0; + } + X = mp_add(X, X, b); + Y = mp_sub(Y, Y, a); + MP_SCAN(&xsc, X); + MP_SCAN(&ysc, Y); + MP_STEP(&xsc); MP_STEP(&ysc); + } + n++; nn++; + } + + if (nn) { + u = mp_lsr(u, u, nn); + if (ext && n) { + X = mp_lsr(X, X, n); + Y = mp_lsr(Y, Y, n); + } + } + } + + /* --- While @v@ is even --- */ + + { + mpscan sc, xsc, ysc; + size_t n = 0, nn = 0; + + MP_SCAN(&sc, v); + MP_SCAN(&xsc, x); MP_SCAN(&ysc, y); + for (;;) { + MP_STEP(&sc); + MP_STEP(&xsc); MP_STEP(&ysc); + if (MP_BIT(&sc)) + break; + if (ext && (MP_BIT(&xsc) | MP_BIT(&ysc))) { + if (n) { + x = mp_lsr(x, x, n); + y = mp_lsr(y, y, n); + n = 0; + } + x = mp_add(x, x, b); + y = mp_sub(y, y, a); + MP_SCAN(&xsc, x); + MP_SCAN(&ysc, y); + MP_STEP(&xsc); MP_STEP(&ysc); + } + n++; nn++; + } + + if (nn) { + v = mp_lsr(v, v, nn); + if (ext && n) { + x = mp_lsr(x, x, n); + y = mp_lsr(y, y, n); + } + } + } + + /* --- End-of-loop fiddling --- */ + + if (MP_CMP(u, >=, v)) { + u = mp_sub(u, u, v); + if (ext) { + X = mp_sub(X, X, x); + Y = mp_sub(Y, Y, y); + } + } else { + v = mp_sub(v, v, u); + if (ext) { + x = mp_sub(x, x, X); + y = mp_sub(y, y, Y); + } + } + + if (MP_CMP(u, ==, MP_ZERO)) + break; + } + + /* --- Write the results out --- */ + + if (gcd) + *gcd = mp_lsl(v, v, shift); + else + MP_DROP(v); + + /* --- Perform a little normalization --- * + * + * Ensure that the coefficient returned is positive, if there is only one. + * If there are two, favour @y@. + */ + + if (ext) { + if (swap) { + mp *t = x; x = y; y = t; + } + if (yy) { + if (y->f & MP_NEG) { + y = mp_add(y, y, a); + x = mp_sub(x, x, b); + } + } else if (x->f & MP_NEG) + x = mp_add(x, x, b); + + if (xx) *xx = x; else MP_DROP(x); + if (yy) *yy = y; else MP_DROP(y); + } + + MP_DROP(u); + MP_DROP(X); MP_DROP(Y); + MP_DROP(a); MP_DROP(b); +} + +/*----- Test rig ----------------------------------------------------------*/ + +#ifdef TEST_RIG + +static int gcd(dstr *v) +{ + int ok = 1; + mp *a = *(mp **)v[0].buf; + mp *b = *(mp **)v[1].buf; + mp *g = *(mp **)v[2].buf; + mp *x = *(mp **)v[3].buf; + mp *y = *(mp **)v[4].buf; + + mp *gg, *xx, *yy; + mp_gcd(&gg, &xx, &yy, a, b); + if (MP_CMP(x, !=, xx)) { + fputs("\n*** mp_gcd(x) failed", stderr); + fputs("\na = ", stderr); mp_writefile(a, stderr, 10); + fputs("\nb = ", stderr); mp_writefile(b, stderr, 10); + fputs("\nexpect = ", stderr); mp_writefile(x, stderr, 10); + fputs("\nresult = ", stderr); mp_writefile(xx, stderr, 10); + fputc('\n', stderr); + ok = 0; + } + if (MP_CMP(y, !=, yy)) { + fputs("\n*** mp_gcd(y) failed", stderr); + fputs("\na = ", stderr); mp_writefile(a, stderr, 10); + fputs("\nb = ", stderr); mp_writefile(b, stderr, 10); + fputs("\nexpect = ", stderr); mp_writefile(y, stderr, 10); + fputs("\nresult = ", stderr); mp_writefile(yy, stderr, 10); + fputc('\n', stderr); + ok = 0; + } + + if (!ok) { + mp *ax = mp_mul(MP_NEW, a, xx); + mp *by = mp_mul(MP_NEW, b, yy); + ax = mp_add(ax, ax, by); + if (MP_CMP(ax, ==, gg)) + fputs("\n*** (Alternative result found.)\n", stderr); + MP_DROP(ax); + MP_DROP(by); + } + + if (MP_CMP(g, !=, gg)) { + fputs("\n*** mp_gcd(gcd) failed", stderr); + fputs("\na = ", stderr); mp_writefile(a, stderr, 10); + fputs("\nb = ", stderr); mp_writefile(b, stderr, 10); + fputs("\nexpect = ", stderr); mp_writefile(g, stderr, 10); + fputs("\nresult = ", stderr); mp_writefile(gg, stderr, 10); + fputc('\n', stderr); + ok = 0; + } + MP_DROP(a); MP_DROP(b); MP_DROP(g); MP_DROP(x); MP_DROP(y); + MP_DROP(gg); MP_DROP(xx); MP_DROP(yy); + return (ok); +} + +static test_chunk tests[] = { + { "gcd", gcd, { &type_mp, &type_mp, &type_mp, &type_mp, &type_mp, 0 } }, + { 0, 0, { 0 } } +}; + +int main(int argc, char *argv[]) +{ + sub_init(); + test_run(argc, argv, tests, SRCDIR "/tests/mp"); + return (0); +} + +#endif + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/mp-io.c b/mp-io.c new file mode 100644 index 0000000..0aa2bca --- /dev/null +++ b/mp-io.c @@ -0,0 +1,149 @@ +/* -*-c-*- + * + * $Id: mp-io.c,v 1.1 1999/11/17 18:02:16 mdw Exp $ + * + * Loading and storing of multiprecision integers + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: mp-io.c,v $ + * Revision 1.1 1999/11/17 18:02:16 mdw + * New multiprecision integer arithmetic suite. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include "mp.h" + +/*----- Main code ---------------------------------------------------------*/ + +/* --- @mp_octets@ --- * + * + * Arguments: @const mp *m@ = a multiprecision integer + * + * Returns: The number of octets required to represent @m@. + * + * Use: Calculates the external storage required for a multiprecision + * integer. + */ + +size_t mp_octets(const mp *m) +{ + size_t sz; + MPX_OCTETS(sz, m->v, m->vl); + return (sz); +} + +/* --- @mp_loadl@ --- * + * + * Arguments: @mp *d@ = destination + * @const void *pv@ = pointer to source data + * @size_t sz@ = size of the source data + * + * Returns: Resulting multiprecision number. + * + * Use: Loads a multiprecision number from an array of octets. The + * first byte in the array is the least significant. More + * formally, if the bytes are %$b_0, b_1, \ldots, b_{n-1}$% + * then the result is %$N = \sum_{0 \le i < n} b_i 2^{8i}$%. + */ + +mp *mp_loadl(mp *d, const void *pv, size_t sz) +{ + MP_MODIFY(d, MPW_RQ(sz)); + mpx_loadl(d->v, d->vl, pv, sz); + mp_shrink(d); + return (d); +} + +/* --- @mp_storel@ --- * + * + * Arguments: @const mp *m@ = source + * @void *pv@ = pointer to output array + * @size_t sz@ = size of the output array + * + * Returns: --- + * + * Use: Stores a multiprecision number in an array of octets. The + * first byte in the array is the least significant. If the + * array is too small to represent the number, high-order bits + * are truncated; if the array is too large, high order bytes + * are filled with zeros. More formally, if the number is + * %$N = \sum{0 \le i} b_i 2^{8i}$% where %$0 \le b_i < 256$%, + * then the array is %$b_0, b_1, \ldots, b_{n-1}$%. + */ + +void mp_storel(const mp *m, void *pv, size_t sz) +{ + mpx_storel(m->v, m->vl, pv, sz); +} + +/* --- @mp_loadb@ --- * + * + * Arguments: @mp *d@ = destination + * @const void *pv@ = pointer to source data + * @size_t sz@ = size of the source data + * + * Returns: Resulting multiprecision number. + * + * Use: Loads a multiprecision number from an array of octets. The + * last byte in the array is the least significant. More + * formally, if the bytes are %$b_{n-1}, b_{n-2}, \ldots, b_0$% + * then the result is %$N = \sum_{0 \le i < n} b_i 2^{8i}$%. + */ + +mp *mp_loadb(mp *d, const void *pv, size_t sz) +{ + MP_MODIFY(d, MPW_RQ(sz)); + mpx_loadb(d->v, d->vl, pv, sz); + mp_shrink(d); + return (d); +} + +/* --- @mp_storeb@ --- * + * + * Arguments: @const mp *m@ = source + * @void *pv@ = pointer to output array + * @size_t sz@ = size of the output array + * + * Returns: --- + * + * Use: Stores a multiprecision number in an array of octets. The + * last byte in the array is the least significant. If the + * array is too small to represent the number, high-order bits + * are truncated; if the array is too large, high order bytes + * are filled with zeros. More formally, if the number is + * %$N = \sum{0 \le i} b_i 2^{8i}$% where %$0 \le b_i < 256$%, + * then the array is %$b_{n-1}, b_{n-2}, \ldots, b_0$%. + */ + +void mp_storeb(const mp *m, void *pv, size_t sz) +{ + mpx_storeb(m->v, m->vl, pv, sz); +} + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/mp-mem.c b/mp-mem.c new file mode 100644 index 0000000..2f75365 --- /dev/null +++ b/mp-mem.c @@ -0,0 +1,187 @@ +/* -*-c-*- + * + * $Id: mp-mem.c,v 1.1 1999/11/17 18:02:16 mdw Exp $ + * + * Memory management for multiprecision numbers + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: mp-mem.c,v $ + * Revision 1.1 1999/11/17 18:02:16 mdw + * New multiprecision integer arithmetic suite. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include "mp.h" + +/*----- Main code ---------------------------------------------------------*/ + +/* --- @mp_create@ --- * + * + * Arguments: @size_t sz@ = size of vector required + * + * Returns: Pointer to pristine new MP structure with enough memory + * bolted onto it. + * + * Use: Creates a new multiprecision integer, initially zero. The + * integer has a single reference. + */ + +mp *mp_create(size_t sz) +{ + mp *m = CREATE(mp); + m->v = MP_ALLOC(sz); + m->vl = m->v + sz; + m->sz = sz; + m->f = MP_UNDEF; + m->ref = 1; + return (m); +} + +/* --- @mp_build@ --- * + * + * Arguments: @mp *m@ = pointer to an MP block to fill in + * @mpw *v@ = pointer to a word array + * @mpw *vl@ = pointer just past end of array + * + * Returns: --- + * + * Use: Creates a multiprecision integer representing some smallish + * number. You must provide storage for the number and dispose + * of it when you've finished with it. The number is marked as + * constant while it exists. + */ + +void mp_build(mp *m, mpw *v, mpw *vl) +{ + m->v = v; + m->vl = vl; + m->sz = vl - v; + m->f = MP_CONST; + m->ref = 1; +} + +/* --- @mp_destroy@ --- * + * + * Arguments: @mp *m@ = pointer to a multiprecision integer + * + * Returns: --- + * + * Use: Destroys a multiprecision integer. The reference count isn't + * checked. Don't use this function if you don't know what + * you're doing: use @mp_drop@ instead. + */ + +void mp_destroy(mp *m) +{ + if (m->f & MP_CONST) + return; + if (m->f & MP_BURN) + memset(m->v, 0, MPWS(m->sz)); + MP_FREE(m->v); + DESTROY(m); +} + +/* --- @mp_copy@ --- * + * + * Arguments: @mp *m@ = pointer to a multiprecision integer + * + * Returns: A copy of the given multiprecision integer. + * + * Use: Copies the given integer. In fact you just get another + * reference to the same old one again. + */ + +mp *mp_copy(mp *m) { return MP_COPY(m); } + +/* --- @mp_drop@ --- * + * + * Arguments: @mp *m@ = pointer to a multiprecision integer + * + * Returns: --- + * + * Use: Drops a reference to an integer which isn't wanted any more. + * If there are no more references, the integer is destroyed. + */ + +void mp_drop(mp *m) { MP_DROP(m); } + +/* --- @mp_split@ --- * + * + * Arguments: @mp *m@ = pointer to a multiprecision integer + * + * Returns: A reference to the same integer, possibly with a different + * address. + * + * Use: Splits off a modifiable version of the integer referred to. + */ + +mp *mp_split(mp *m) { MP_SPLIT(m); return (m); } + +/* --- @mp_resize@ --- * + * + * Arguments: @mp *m@ = pointer to a multiprecision integer + * @size_t sz@ = new size + * + * Returns: --- + * + * Use: Resizes the vector containing the integer's digits. The new + * size must be at least as large as the current integer's + * length. This isn't really intended for client use. + */ + +void mp_resize(mp *m, size_t sz) { MP_RESIZE(m, sz); } + +/* --- @mp_ensure@ --- * + * + * Arguments: @mp *m@ = pointer to a multiprecision integer + * @size_t sz@ = required size + * + * Returns: --- + * + * Use: Ensures that the integer has enough space for @sz@ digits. + * The value is not changed. + */ + +void mp_ensure(mp *m, size_t sz) { MP_ENSURE(m, sz); } + +/* --- @mp_modify@ --- * + * + * Arguments: @mp *m@ = pointer to a multiprecision integer + * @size_t sz@ = size required + * + * Returns: Pointer to the integer (possibly different). + * + * Use: Prepares an integer to be overwritten. It's split off from + * other references to the same integer, and sufficient space is + * allocated. + */ + +mp *mp_modify(mp *m, size_t sz) { MP_MODIFY(m, sz); return (m); } + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/mp-misc.c b/mp-misc.c new file mode 100644 index 0000000..7f5513f --- /dev/null +++ b/mp-misc.c @@ -0,0 +1,109 @@ +/* -*-c-*- + * + * $Id: mp-misc.c,v 1.1 1999/11/17 18:02:16 mdw Exp $ + * + * Miscellaneous multiprecision support functions + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: mp-misc.c,v $ + * Revision 1.1 1999/11/17 18:02:16 mdw + * New multiprecision integer arithmetic suite. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include "mp.h" + +/*----- Paranoia management -----------------------------------------------*/ + +/* --- @mp_burn@ --- * + * + * Arguments: @mp *m@ = pointer to a multiprecision integer + * + * Returns: --- + * + * Use: Marks the integer as `burn-after-use'. When the integer's + * memory is deallocated, it is deleted so that traces can't + * remain in the swap file. In theory. + */ + +void mp_burn(mp *m) +{ + m->f |= MP_BURN; +} + +/*----- Basic manipulation ------------------------------------------------*/ + +/* --- @mp_shrink@ --- * + * + * Arguments: @mp *m@ = pointer to a multiprecision integer + * + * Returns: --- + * + * Use: Reduces the recorded length of an integer. This doesn't + * reduce the amount of memory used, although it can improve + * performance a bit. To reduce memory, use @mp_minimize@ + * instead. This can't change the value of an integer, and is + * therefore safe to use even when there are multiple + * references. + */ + +void mp_shrink(mp *m) { MP_SHRINK(m); } + +/* --- @mp_minimize@ --- * + * + * Arguments: @mp *m@ = pointer to a multiprecision integer + * + * Returns: --- + * + * Use: Reduces the amount of memory an integer uses. It's best to + * do this to numbers which aren't going to change in the + * future. + */ + +void mp_minimize(mp *m) +{ + MP_SHRINK(m); + MP_RESIZE(m, MP_LEN(m)); +} + +/*----- Bit scanning ------------------------------------------------------*/ + +/* --- @mp_scan@ --- * + * + * Arguments: @mpscan *sc@ = pointer to bitscanner block + * @const mp *m@ = pointer to a multiprecision integer + * + * Returns: --- + * + * Use: Initializes a bitscanner on a multiprecision integer. + */ + +void mp_scan(mpscan *sc, const mp *m) { MP_SCAN(sc, m); } + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/mp-test.c b/mp-test.c new file mode 100644 index 0000000..018b537 --- /dev/null +++ b/mp-test.c @@ -0,0 +1,72 @@ +/* -*-c-*- + * + * $Id: mp-test.c,v 1.1 1999/11/17 18:02:16 mdw Exp $ + * + * Testing functionality for multiprecision integers + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: mp-test.c,v $ + * Revision 1.1 1999/11/17 18:02:16 mdw + * New multiprecision integer arithmetic suite. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include +#include +#include + +#include +#include +#include +#include + +#include "mp.h" +#include "mptext.h" + +/*----- The `MP' testrig data type ----------------------------------------*/ + +static void mp_cvt(const char *buf, dstr *d) +{ + mp *m; + DENSURE(d, sizeof(mp *)); + m = mp_readstring(MP_NEW, (char *)buf, 0, 0); + if (!m) + die(1, "bad integer `%s'", buf); + *(mp **)d->buf = m; +} + +static void mp_dump(dstr *d, FILE *fp) +{ + mp *m = *(mp **)d->buf; + mp_writefile(m, fp, 10); +} + +const test_type type_mp = { mp_cvt, mp_dump }; + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/mp.h b/mp.h index fcccee2..72aa675 100644 --- a/mp.h +++ b/mp.h @@ -1,8 +1,8 @@ /* -*-c-*- * - * $Id: mp.h,v 1.1 1999/09/03 08:41:12 mdw Exp $ + * $Id: mp.h,v 1.2 1999/11/17 18:02:16 mdw Exp $ * - * Multiprecision arithmetic + * Simple multiprecision arithmetic * * (c) 1999 Straylight/Edgeware */ @@ -30,8 +30,8 @@ /*----- Revision history --------------------------------------------------* * * $Log: mp.h,v $ - * Revision 1.1 1999/09/03 08:41:12 mdw - * Initial import. + * Revision 1.2 1999/11/17 18:02:16 mdw + * New multiprecision integer arithmetic suite. * */ @@ -44,401 +44,592 @@ /*----- Header files ------------------------------------------------------*/ -#include +#include #include -#ifndef MPTYPES_H -# include "mptypes.h" +#include + +#ifndef MPW_H +# include "mpw.h" #endif #ifndef MPX_H # include "mpx.h" #endif -#ifndef MPSCAN_H -# include "mpscan.h" -#endif - /*----- Data structures ---------------------------------------------------*/ typedef struct mp { - mpw *v; /* Vector of words */ - unsigned f; /* Various flags */ - size_t sz; /* Size allocated to word vector */ - size_t len; /* Length of current word vector */ + mpw *v, *vl; + size_t sz; + unsigned f; + unsigned ref; } mp; -enum { - MPF_SIGN = 1, /* Sign bit */ - MPF_BURN = 2 /* Burn word vector after use */ -}; +#define MP_NEG 1u +#define MP_BURN 2u +#define MP_CONST 4u +#define MP_UNDEF 8u -typedef struct mp_bitscan { - const mp *x; /* Pointer to target MP */ - size_t i; /* Index into MP vector */ - mpw w; /* Current value being examined */ - unsigned bits; /* Number of bits left in @w@ */ -} mp_bitscan; +/*----- Useful constants --------------------------------------------------*/ -/*----- External variables ------------------------------------------------*/ +extern mp mp_const[]; -extern mp mp_std; +#define MP_ZERO (&mp_const[0]) +#define MP_ONE (&mp_const[1]) +#define MP_TWO (&mp_const[2]) +#define MP_THREE (&mp_const[3]) +#define MP_FOUR (&mp_const[4]) +#define MP_FIVE (&mp_const[5]) +#define MP_TEN (&mp_const[6]) +#define MP_MONE (&mp_const[7]) -#define MP_ZERO (mp_std + 0) -#define MP_ONE (mp_std + 1) -#define MP_TWO (mp_std + 2) -#define MP_THREE (mp_std + 3) -#define MP_FOUR (mp_std + 4) -#define MP_FIVE (mp_std + 5) -#define MP_TEN (mp_std + 6) -#define MP_MONE (mp_std + 7) +#define MP_NEW ((mp *)0) -/*----- Memory allocation and low-level fiddling --------------------------*/ +/*----- Memory allocation hooks -------------------------------------------*/ -/* --- @mp_create@ --- * +#ifndef MPARENA_H +# include "mparena.h" +#endif + +/* --- @MP_ARENA@ --- * * - * Arguments @mp *x@ = pointer to MP head + * This selects where memory is allocated from. Tweak to use more fancy + * things like custom arenas. + */ + +#ifndef MP_ARENA +# define MP_ARENA MPARENA_GLOBAL +#endif + +/* --- @MP_ALLOC@ --- * * - * Returns: --- + * Arguments: @size_t sz@ = size required + * + * Returns: Pointer to an allocated vector of the requested size. * - * Use: Initializes an MP ready for use. The initial value is zero. + * Use: Hook for vector allocation. */ -#define MP_INIT { 0, 0, 0, 0 } - -extern void mp_create(mp */*x*/); +#ifndef MP_ALLOC +# define MP_ALLOC(sz) mpalloc(MP_ARENA, (sz)) +#endif -/* --- @MP_BURN@ --- * +/* --- @MP_FREE@ --- * + * + * Arguments: @mpw *v@ = pointer to vector * - * Arguments: @x@ = pointer to MP head + * Returns: --- * - * Use: Burns the contents of the MP, if it contains sensitive data. + * Use: Hook for vector deallocation. */ -#define MP_BURN(x) do { \ - mp *_y = (x) \ - if (_y->v && _y->f & mpf_burn) { \ - memset(_y->v, 0, _y->sz * sizeof(mpw)); \ - _y->f &= ~MPF_BURN; \ - } \ -} while (0) +#ifndef MP_FREE +# define MP_FREE(v) mpfree(MP_ARENA, (v)) +#endif -/* --- @mp_free@ --- * +/*----- Paranoia management -----------------------------------------------*/ + +/* --- @mp_burn@ --- * * - * Arguments: @mp *x@ = pointer to MP head + * Arguments: @mp *m@ = pointer to a multiprecision integer * * Returns: --- * - * Use: Releases the memory used by an MP. + * Use: Marks the integer as `burn-after-use'. When the integer's + * memory is deallocated, it is deleted so that traces can't + * remain in the swap file. In theory. */ -#define MP_DESTROY(x) do { \ - mp *_x = (x); \ - MP_BURN(_x); \ - if (_x->v) \ - free(_x->v); \ - _x->v = 0; \ - _x->f = 0; \ - _x->sz = 0; \ - _x->len = 0; \ -} while (0) - -extern void mp_free(mp */*x*/); +extern void mp_burn(mp */*m*/); -/* --- @MP_ENSURE@ --- * +/*----- Trivial macros ----------------------------------------------------*/ + +/* --- @MP_LEN@ --- * * - * Arguments: @x@ = pointer to MP head - * @len@ = length required (in words) + * Arguments: @mp *m@ = pointer to a multiprecision integer * - * Use: Ensures that the MP has enough memory to store a @len@-word - * value. + * Returns: Length of the integer, in words. */ -#define MP_ENSURE(x, len) do { \ - mp *_x = (x); \ - size_t _len = (len); \ - if (_x->sz < _len) \ - mp_resize(_x, _len); \ -} while (0) +#define MP_LEN(m) ((m)->vl - ((m)->v)) -/* --- @mp_resize@ --- * +/*----- Memory management and reference counting --------------------------*/ + +/* --- @mp_create@ --- * * - * Arguments: @mp *x@ = pointer to MP head - * @size_t sz@ = size required (in words) + * Arguments: @size_t sz@ = size of vector required * - * Returns: --- + * Returns: Pointer to pristine new MP structure with enough memory + * bolted onto it. * - * Use: Resizes the MP so that its word vector has space for - * exactly @sz@ words. + * Use: Creates a new multiprecision integer with indeterminate + * contents. The integer has a single reference. */ -extern void mp_resize(mp */*x*/, size_t /*sz*/); +extern mp *mp_create(size_t /*sz*/); -/* --- @mp_norm@ --- * +/* --- @mp_build@ --- * * - * Arguments: @mp *x@ = pointer to MP head + * Arguments: @mp *m@ = pointer to an MP block to fill in + * @mpw *v@ = pointer to a word array + * @mpw *vl@ = pointer just past end of array * * Returns: --- * - * Use: `Normalizes' an MP. Fixes the @len@ field so that it's - * correct. Assumes that @len@ is either already correct or - * too large. + * Use: Creates a multiprecision integer representing some smallish + * number. You must provide storage for the number and dispose + * of it when you've finished with it. The number is marked as + * constant while it exists. */ -#define MP_NORM(x) do { \ - mp *_y = (x); \ - MPX_LEN(_y->len, _y->v, _y->len); \ -} while (0) +extern void mp_build(mp */*m*/, mpw */*v*/, mpw */*vl*/); -extern void mp_norm(mp */*x*/); - -/* --- @mp_dump@ --- * +/* --- @mp_destroy@ --- * * - * Arguments: @mp *x@ = pointer to MP head - * @FILE *fp@ = pointer to stream to write on + * Arguments: @mp *m@ = pointer to a multiprecision integer * * Returns: --- * - * Use: Dumps an MP to a stream. + * Use: Destroys a multiprecision integer. The reference count isn't + * checked. Don't use this function if you don't know what + * you're doing: use @mp_drop@ instead. */ -extern void mp_dump(mp */*x*/, FILE */*fp*/); +extern void mp_destroy(mp */*m*/); -/* --- @MP_PARANOID@ --- * +/* --- @mp_copy@ --- * * - * Arguments: @x@ = pointer to MP head + * Arguments: @mp *m@ = pointer to a multiprecision integer * - * Use: Marks the MP as containing sensitive data which should be - * burnt when no longer required. + * Returns: A copy of the given multiprecision integer. + * + * Use: Copies the given integer. In fact you just get another + * reference to the same old one again. */ -#define MP_PARANOID(x) ((x)->f |= MPF_BURN) +extern mp *mp_copy(mp */*m*/); -/* --- @mp_copy@ --- * +#define MP_COPY(m) ((m)->ref++, (m)) + +/* --- @mp_drop@ --- * * - * Arguments: @mp *d@ = pointer to MP head for destination - * @const mp *s@ = pointer to MP head for source + * Arguments: @mp *m@ = pointer to a multiprecision integer * * Returns: --- * - * Use: Copies an MP. + * Use: Drops a reference to an integer which isn't wanted any more. + * If there are no more references, the integer is destroyed. */ -extern void mp_copy(mp */*d*/, const mp */*s*/); +extern void mp_drop(mp */*m*/); -/* --- @mp_bits@ --- * +#define MP_DROP(m) do { \ + mp *_mm = (m); \ + if (_mm->ref > 1) \ + _mm->ref--; \ + else if (!(_mm->f & MP_CONST)) \ + mp_destroy(_mm); \ +} while (0) + +/* --- @mp_split@ --- * * - * Arguments: @mp *x@ = pointer to MP head + * Arguments: @mp *m@ = pointer to a multiprecision integer * - * Returns: Length of the number in bits. + * Returns: A reference to the same integer, possibly with a different + * address. * - * Use: Calculates the number of bits required to represent a number. - * The number must be normalized. + * Use: Splits off a modifiable version of the integer referred to. */ -unsigned long mp_bits(mp */*x*/); +extern mp *mp_split(mp */*m*/); + +#define MP_SPLIT(m) do { \ + mp *_mm = (m); \ + if ((_mm->f & MP_CONST) || _mm->ref != 1) { \ + mp *_dd = mp_create(_mm->sz); \ + _dd->vl = _dd->v + MP_LEN(_mm); \ + _dd->f = _mm->f & (MP_NEG | MP_BURN); \ + memcpy(_dd->v, _mm->v, MPWS(MP_LEN(_mm))); \ + _dd->ref = 1; \ + _mm->ref--; \ + (m) = _dd; \ + } \ +} while (0) -/* --- @mp_octets@ --- * +/* --- @mp_resize@ --- * * - * Arguments: @mp *x@ = pointer to MP head + * Arguments: @mp *m@ = pointer to a multiprecision integer + * @size_t sz@ = new size * - * Returns: Length of number in octets. + * Returns: --- * - * Use: Calculates the number of octets required to represent a - * number. The number must be normalized. + * Use: Resizes the vector containing the integer's digits. The new + * size must be at least as large as the current integer's + * length. The integer's length is increased and new digits are + * filled with zeroes. This isn't really intended for client + * use. */ -extern size_t mp_octets(mp */*x*/); - -/*----- Loading and storing as binary data --------------------------------*/ +extern void mp_resize(mp */*m*/, size_t /*sz*/); + +#define MP_RESIZE(m, ssz) do { \ + mp *_m = (m); \ + size_t _sz = (ssz); \ + size_t _len = MP_LEN(_m); \ + mpw *_v = MP_ALLOC(_sz); \ + memcpy(_v, _m->v, MPWS(_len)); \ + if (_m->f & MP_BURN) \ + memset(_m->v, 0, MPWS(_m->sz)); \ + MP_FREE(_m->v); \ + _m->v = _v; \ + _m->vl = _v + _len; \ + _m->sz = _sz; \ +} while (0) -/* --- @mp_storel@ --- * +/* --- @mp_ensure@ --- * * - * Arguments: @mp *x@ = pointer to MP head - * @octet *p@ = pointer to octet array - * @size_t sz@ = size of octet array + * Arguments: @mp *m@ = pointer to a multiprecision integer + * @size_t sz@ = required size * * Returns: --- * - * Use: Stores an MP in an octet array, least significant octet - * first. High-end octets are silently discarded if there - * isn't enough space for them. + * Use: Ensures that the integer has enough space for @sz@ digits. + * The value is not changed. */ -extern void mp_storel(mp */*x*/, octet */*p*/, size_t /*sz*/); +extern void mp_ensure(mp */*m*/, size_t /*sz*/); + +#define MP_ENSURE(m, ssz) do { \ + mp *_mm = (m); \ + size_t _ssz = (ssz); \ + size_t _len = MP_LEN(_mm); \ + if (_ssz > _mm->sz) \ + MP_RESIZE(_mm, _ssz); \ + if (!(_mm->f & MP_UNDEF) && _ssz > _len) { \ + memset(_mm->vl, 0, MPWS(_ssz - _len)); \ + _mm->vl = _mm->v + _ssz; \ + } \ +} while (0) -/* --- @mp_loadl@ --- * +/* --- @mp_modify@ --- * * - * Arguments: @mp *x@ = pointer to MP head - * @const octet *p@ = pointer to octet array - * @size_t sz@ = size of octet array + * Arguments: @mp *m@ = pointer to a multiprecision integer + * @size_t sz@ = size required * - * Returns: --- + * Returns: Pointer to the integer (possibly different). * - * Use: Loads an MP in an octet array, least significant octet - * first. + * Use: Prepares an integer to be overwritten. It's split off from + * other references to the same integer, and sufficient space is + * allocated. */ -extern void mp_loadl(mp */*x*/, const octet */*p*/, size_t /*sz*/); +extern mp *mp_modify(mp */*m*/, size_t /*sz*/); -/* --- @mp_storeb@ --- * +#define MP_MODIFY(m, sz) do { \ + size_t _rq = (sz); \ + mp *_m = (m); \ + if (_m == MP_NEW) \ + _m = mp_create(_rq); \ + else { \ + MP_SPLIT(_m); \ + MP_ENSURE(_m, _rq); \ + } \ + _m->vl = _m->v + _rq; \ + (m) = _m; \ +} while (0) + +/*----- Size manipulation -------------------------------------------------*/ + +/* --- @mp_shrink@ --- * * - * Arguments: @mp *x@ = pointer to MP head - * @octet *p@ = pointer to octet array - * @size_t sz@ = size of octet array + * Arguments: @mp *m@ = pointer to a multiprecision integer * * Returns: --- * - * Use: Stores an MP in an octet array, most significant octet - * first. High-end octets are silently discarded if there - * isn't enough space for them. + * Use: Reduces the recorded length of an integer. This doesn't + * reduce the amount of memory used, although it can improve + * performance a bit. To reduce memory, use @mp_minimize@ + * instead. This can't change the value of an integer, and is + * therefore safe to use even when there are multiple + * references. */ -extern void mp_storeb(mp */*x*/, octet */*p*/, size_t /*sz*/); +extern void mp_shrink(mp */*m*/); -/* --- @mp_loadb@ --- * +#define MP_SHRINK(m) do { \ + mp *_mm = (m); \ + MPX_SHRINK(_mm->v, _mm->vl); \ + if (!MP_LEN(_mm)) \ + _mm->f &= ~MP_NEG; \ +} while (0) + +/* --- @mp_minimize@ --- * * - * Arguments: @mp *x@ = pointer to MP head - * @const octet *p@ = pointer to octet array - * @size_t sz@ = size of octet array + * Arguments: @mp *m@ = pointer to a multiprecision integer * * Returns: --- * - * Use: Loads an MP in an octet array, most significant octet - * first. + * Use: Reduces the amount of memory an integer uses. It's best to + * do this to numbers which aren't going to change in the + * future. */ -extern void mp_loadb(mp */*x*/, const octet */*p*/, size_t /*sz*/); +extern void mp_minimize(mp */*m*/); -/*----- Iterating through bits --------------------------------------------*/ +/*----- Bit scanning ------------------------------------------------------*/ -/* --- @mp_mkbitscan@ --- * +#ifndef MPSCAN_H +# include "mpscan.h" +#endif + +/* --- @mp_scan@ --- * * - * Arguments: @mp_bitscan *sc@ = pointer to bitscan object - * @const mp *x@ = pointer to MP head + * Arguments: @mpscan *sc@ = pointer to bitscanner block + * @const mp *m@ = pointer to a multiprecision integer * * Returns: --- * - * Use: Initializes a bitscan object. + * Use: Initializes a bitscanner on a multiprecision integer. */ -extern void mp_mkbitscan(mp_bitscan */*sc*/, const mp */*x*/); +extern void mp_scan(mpscan */*sc*/, const mp */*m*/); + +#define MP_SCAN(sc, m) do { \ + mp *_mm = (m); \ + mpscan *_sc = (sc); \ + MPSCAN_INITX(_sc, _mm->v, _mm->vl); \ +} while (0) + +/* --- Other bitscanning aliases --- */ + +#define mp_step mpscan_step +#define mp_bit mpscan_bit + +#define MP_STEP MPSCAN_STEP +#define MP_BIT MPSCAN_BIT + +/*----- Loading and storing -----------------------------------------------*/ -/* --- @mp_bstep@ --- * +/* --- @mp_octets@ --- * * - * Arguments: @mp_bitscan *sc@ = pointer to bitscanner object + * Arguments: @const mp *m@ = a multiprecision integer * - * Returns: Nonzero if there is another bit to read. + * Returns: The number of octets required to represent @m@. * - * Use: Steps on to the next bit, and tells the caller whether one - * exists. + * Use: Calculates the external storage required for a multiprecision + * integer. */ -extern int mp_bstep(mp_bitscan */*sc*/); +extern size_t mp_octets(const mp *m); -/* --- @mp_bit@ --- * +/* --- @mp_loadl@ --- * * - * Arguments: @const mp_bitscan *sc@ = pointer to bitscanner + * Arguments: @mp *d@ = destination + * @const void *pv@ = pointer to source data + * @size_t sz@ = size of the source data * - * Returns: Current bit value. + * Returns: Resulting multiprecision number. * - * Use: Returns the value of the current bit. + * Use: Loads a multiprecision number from an array of octets. The + * first byte in the array is the least significant. More + * formally, if the bytes are %$b_0, b_1, \ldots, b_{n-1}$% + * then the result is %$N = \sum_{0 \le i < n} b_i 2^{8i}$%. */ -#define MP_BIT(sc) ((sc)->w & 1) +extern mp *mp_loadl(mp */*d*/, const void */*pv*/, size_t /*sz*/); -extern int mp_bit(const mp_bitscan */*sc*/); +/* --- @mp_storel@ --- * + * + * Arguments: @const mp *m@ = source + * @void *pv@ = pointer to output array + * @size_t sz@ = size of the output array + * + * Returns: --- + * + * Use: Stores a multiprecision number in an array of octets. The + * first byte in the array is the least significant. If the + * array is too small to represent the number, high-order bits + * are truncated; if the array is too large, high order bytes + * are filled with zeros. More formally, if the number is + * %$N = \sum{0 \le i} b_i 2^{8i}$% where %$0 \le b_i < 256$%, + * then the array is %$b_0, b_1, \ldots, b_{n-1}$%. + */ -/*----- Shifting ----------------------------------------------------------*/ +extern void mp_storel(const mp */*m*/, void */*pv*/, size_t /*sz*/); -/* --- @mp_lsl@ --- * +/* --- @mp_loadb@ --- * * - * Arguments: @mp *d@ = pointer to MP head of destination - * @const mp *x@ = pointer to MP head of source - * @size_t n@ = number of bits to shift + * Arguments: @mp *d@ = destination + * @const void *pv@ = pointer to source data + * @size_t sz@ = size of the source data * - * Returns: --- + * Returns: Resulting multiprecision number. * - * Use: Shifts a number left by a given number of bit positions. + * Use: Loads a multiprecision number from an array of octets. The + * last byte in the array is the least significant. More + * formally, if the bytes are %$b_{n-1}, b_{n-2}, \ldots, b_0$% + * then the result is %$N = \sum_{0 \le i < n} b_i 2^{8i}$%. */ -extern void mp_lsl(mp */*d*/, const mp */*x*/, size_t /*n*/); +extern mp *mp_loadb(mp */*d*/, const void */*pv*/, size_t /*sz*/); -/* --- @mp_lsr@ --- * +/* --- @mp_storeb@ --- * * - * Arguments: @mp *d@ = pointer to MP head of destination - * @const mp *x@ = pointer to MP head of source - * @size_t n@ = number of bits to shift + * Arguments: @const mp *m@ = source + * @void *pv@ = pointer to output array + * @size_t sz@ = size of the output array * * Returns: --- * - * Use: Shifts a number right by a given number of bit positions. + * Use: Stores a multiprecision number in an array of octets. The + * last byte in the array is the least significant. If the + * array is too small to represent the number, high-order bits + * are truncated; if the array is too large, high order bytes + * are filled with zeros. More formally, if the number is + * %$N = \sum{0 \le i} b_i 2^{8i}$% where %$0 \le b_i < 256$%, + * then the array is %$b_{n-1}, b_{n-2}, \ldots, b_0$%. */ -extern void mp_lsr(mp */*d*/, const mp */*x*/, size_t /*n*/); +extern void mp_storeb(const mp */*m*/, void */*pv*/, size_t /*sz*/); -/*----- Unsigned arithmetic -----------------------------------------------*/ +/*----- Simple arithmetic -------------------------------------------------*/ -/* --- @mp_uadd@ --- * +/* --- @mp_2c@ --- * * - * Arguments: @const mp *d@ = pointers to MP head of destination - * @const mp *x, *y@ = pointers to MP heads of operands + * Arguments: @mp *d@ = destination + * @mp *a@ = source * - * Returns: --- + * Returns: Result, @a@ converted to two's complement notation. + */ + +extern mp *mp_2c(mp */*d*/, mp */*a*/); + +/* --- @mp_sm@ --- * + * + * Arguments: @mp *d@ = destination + * @mp *a@ = source * - * Use: Performs unsigned MP addition. + * Returns: Result, @a@ converted to the native signed-magnitude + * notation. */ -extern void mp_uadd(mp */*d*/, const mp */*x*/, const mp */*y*/); +extern mp *mp_sm(mp */*d*/, mp */*a*/); -/* --- @mp_usub@ --- * +/* --- @mp_lsl@ --- * * - * Arguments: @const mp *d@ = pointers to MP head of destination - * @const mp *x, *y@ = pointers to MP heads of operands + * Arguments: @mp *d@ = destination + * @const mp *a@ = source + * @size_t n@ = number of bits to move * - * Returns: --- + * Returns: Result, @a@ shifted left by @n@. + */ + +extern mp *mp_lsl(mp */*d*/, const mp */*a*/, size_t /*n*/); + +/* --- @mp_lsr@ --- * + * + * Arguments: @mp *d@ = destination + * @const mp *a@ = source + * @size_t n@ = number of bits to move * - * Use: Performs unsigned MP subtraction. + * Returns: Result, @a@ shifted left by @n@. */ -extern void mp_usub(mp */*d*/, const mp */*x*/, const mp */*y*/); +extern mp *mp_lsr(mp */*d*/, const mp */*a*/, size_t /*n*/); -/* --- @mp_ucmp@ --- * +/* --- @mp_cmp@ --- * * - * Arguments: @const mp *x, *y@ = pointers to MP heads of operands + * Arguments: @const mp *a, *b@ = two numbers * - * Returns: Less than, equal to, or greater than zero. + * Returns: Less than, equal to or greater than zero, according to + * whether @a@ is less than, equal to or greater than @b@. + */ + +extern int mp_cmp(const mp */*a*/, const mp */*b*/); + +#define MP_CMP(a, op, b) (mp_cmp((a), (b)) op 0) + +/* --- @mp_add@ --- * * - * Use: Performs unsigned MP comparison. + * Arguments: @mp *d@ = destination + * @const mp *a, *b@ = sources + * + * Returns: Result, @a@ added to @b@. */ -#define MP_UCMP(x, op, y) (mp_ucmp((x), (y)) op 0) +extern mp *mp_add(mp */*d*/, const mp */*a*/, const mp */*b*/); -extern int mp_ucmp(const mp */*x*/, const mp */*y*/); +/* --- @mp_sub@ --- * + * + * Arguments: @mp *d@ = destination + * @const mp *a, *b@ = sources + * + * Returns: Result, @b@ subtracted from @a@. + */ -/* --- @mp_umul@ --- * +extern mp *mp_sub(mp */*d*/, const mp */*a*/, const mp */*b*/); + +/* --- @mp_mul@ --- * * - * Arguments: @mp *d@ = pointer to MP head of destination - * @const mp *x, *y@ = pointes to MP heads of operands + * Arguments: @mp *d@ = destination + * @const mp *a, *b@ = sources * - * Returns: --- + * Returns: Result, @a@ multiplied by @b@. + */ + +extern mp *mp_mul(mp */*d*/, const mp */*a*/, const mp */*b*/); + +/* --- @mp_sqr@ --- * + * + * Arguments: @mp *d@ = destination + * @const mp *a@ = source + * + * Returns: Result, @a@ squared. + */ + +extern mp *mp_sqr(mp */*d*/, const mp */*a*/); + +/* --- @mp_div@ --- * * - * Use: Performs unsigned MP multiplication. + * Arguments: @mp **qq, **rr@ = destination, quotient and remainder + * @const mp *a, *b@ = sources + * + * Use: Calculates the quotient and remainder when @a@ is divided by + * @b@. */ -extern void mp_umul(mp */*d*/, const mp */*x*/, const mp */*y*/); +extern void mp_div(mp **/*qq*/, mp **/*rr*/, + const mp */*a*/, const mp */*b*/); + +/*----- More advanced algorithms ------------------------------------------*/ -/* --- @mp_udiv@ --- * +/* --- @mp_gcd@ --- * * - * Arguments: @mp *q, *r@ = pointers to MP heads for quotient, remainder - * @const mp *x, *y@ = pointers to MP heads for operands + * Arguments: @mp **gcd, **xx, **yy@ = where to write the results + * @mp *a, *b@ = sources (must be nonzero) * * Returns: --- * - * Use: Performs unsigned MP division. + * Use: Calculates @gcd(a, b)@, and two numbers @x@ and @y@ such that + * @ax + by = gcd(a, b)@. This is useful for computing modular + * inverses. Neither @a@ nor @b@ may be zero. Note that, + * unlike @mp_div@ for example, it is not possible to specify + * explicit destinations -- new MPs are always allocated. */ -extern void mp_udiv(mp */*q*/, mp */*r*/, const mp */*x*/, const mp */*y*/); +extern void mp_gcd(mp **/*gcd*/, mp **/*xx*/, mp **/*yy*/, + mp */*a*/, mp */*b*/); + +/*----- Test harness support ----------------------------------------------*/ + +#include + +#ifndef MPTEXT_H +# include "mptext.h" +#endif + +extern const test_type type_mp; /*----- That's all, folks -------------------------------------------------*/ diff --git a/mpalloc.h b/mpalloc.h new file mode 100644 index 0000000..a56ed00 --- /dev/null +++ b/mpalloc.h @@ -0,0 +1,127 @@ +/* -*-c-*- + * + * $Id: mpalloc.h,v 1.1 1999/11/17 18:02:16 mdw Exp $ + * + * Allocation and freeing of MP buffers + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: mpalloc.h,v $ + * Revision 1.1 1999/11/17 18:02:16 mdw + * New multiprecision integer arithmetic suite. + * + */ + +#ifndef MPARENA_H +#define MPARENA_H + +#ifdef __cplusplus + extern "C" { +#endif + +/*----- Header files ------------------------------------------------------*/ + +#ifndef MPW_H +# include "mpw.h" +#endif + +/*----- Data structures ---------------------------------------------------*/ + +typedef struct mparena_node { + struct mparena_node *left, *right; + mpw *v; +} mparena_node; + +typedef struct mparena { + mparena_node *root; +} mparena_arena; + +/*----- Magical constants -------------------------------------------------*/ + +#define MPARENA_GLOBAL ((mparena *)0) + +/*----- Functions provided ------------------------------------------------*/ + +/* --- @mparena_create@ --- * + * + * Arguments: @mparena *a@ = pointer to arena block + * + * Returns: --- + * + * Use: Initializes an MP arena so that blocks can be allocated from + * it. + */ + +extern void mparena_create(mparena */*a*/); + +#define MPARENA_INIT { 0 } + +/* --- @mparena_destroy@ --- * + * + * Arguments: @mparena *a@ = pointer to arena block + * + * Returns: --- + * + * Use: Frees an MP arena, and all the vectors held within it. The + * blocks which are currently allocated can be freed into some + * other arena. + */ + +extern void mparena_destroy(mparena */*a*/); + +/* --- @mp_alloc@ --- * + * + * Arguments: @mparena *a@ = pointer to arena block + * @size_t n@ = number of digits required + * + * Returns: Pointer to a suitably sized block. + * + * Use: Allocates a lump of data suitable for use as an array of MP + * digits. + */ + +extern mpw *mp_alloc(mparena */*a*/, size_t /*n*/); + +/* --- @mp_free@ --- * + * + * Arguments: @mparena *a@ = pointer to arena block + * @mpw *v@ = pointer to allocated vector + * + * Returns: --- + * + * Use: Returns an MP vector to an arena. It doesn't have to be + * returned to the arena from which it was allocated. + */ + +extern mpw *mp_free(mparena */*a*/, mpw */*v*/); + +/*----- That's all, folks -------------------------------------------------*/ + +#ifdef __cplusplus + } +#endif + +#endif diff --git a/mparena.c b/mparena.c new file mode 100644 index 0000000..4e02454 --- /dev/null +++ b/mparena.c @@ -0,0 +1,266 @@ +/* -*-c-*- + * + * $Id: mparena.c,v 1.1 1999/11/17 18:02:16 mdw Exp $ + * + * Allocation and freeing of MP buffers + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: mparena.c,v $ + * Revision 1.1 1999/11/17 18:02:16 mdw + * New multiprecision integer arithmetic suite. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include +#include +#include + +#include +#include + +#include "mparena.h" + +/*----- Default allocator -------------------------------------------------*/ + +static void *defalloc(mparena *a, size_t sz) { return xmalloc(sz); } +static void deffree(mparena *a, void *p) { free(p); } + +mparena_ops mparena_defops = { defalloc, deffree }; + +/*----- Static variables --------------------------------------------------*/ + +static mparena arena = { 0, &mparena_defops }; + +#define MPARENA_RESOLVE(a) do { \ + if ((a) == MPARENA_GLOBAL) \ + (a) = &arena; \ +} while (0) + +/*----- Main code ---------------------------------------------------------*/ + +/* --- @tdump@ --- * + * + * Arguments: @mparena_node *n@ = pointer to tree node to dump + * + * Returns: --- + * + * Use: Recursively dumps out the allocation tree. + */ + +static void tdump(mparena_node *n) +{ + if (!n) + putchar('*'); + else { + putchar('('); + tdump(n->left); + printf(", %u, ", n->v[0]); + tdump(n->right); + putchar(')'); + } +} + +/* --- @mparena_create@ --- * + * + * Arguments: @mparena *a@ = pointer to arena block + * + * Returns: --- + * + * Use: Initializes an MP arena so that blocks can be allocated from + * it. + */ + +void mparena_create(mparena *a) +{ + a->root = 0; + a->ops = &mparena_defops; +} + +/* --- @mparena_setops@ --- * + * + * Arguments: @mparena *a@ = pointer to arena block + * @mparena_ops *ops@ = pointer to operations block or null + * + * Returns: The previous operations block. + * + * Use: Sets or queries the operations attached to an arena. + */ + +mparena_ops *mparena_setops(mparena *a, mparena_ops *ops) +{ + mparena_ops *o; + MPARENA_RESOLVE(a); + o = a->ops; + if (ops) + a->ops = ops; + return (0); +} + +/* --- @mparena_destroy@ --- * + * + * Arguments: @mparena *a@ = pointer to arena block + * + * Returns: --- + * + * Use: Frees an MP arena, and all the vectors held within it. The + * blocks which are currently allocated can be freed into some + * other arena. + */ + +static void tfree(mparena *a, mparena_node *n) +{ + a->ops->free(a, n->v); + if (n->left) + tfree(a, n->left); + if (n->right) + tfree(a, n->right); + DESTROY(n); +} + +void mparena_destroy(mparena *a) +{ + tfree(a, a->root); + a->root = 0; +} + +/* --- @mpalloc@ --- * + * + * Arguments: @mparena *a@ = pointer to arena block + * @size_t sz@ = number of digits required + * + * Returns: Pointer to a suitably sized block. + * + * Use: Allocates a lump of data suitable for use as an array of MP + * digits. + */ + +mpw *mpalloc(mparena *a, size_t sz) +{ + mparena_node **nn, *n; + mpw *v; + + MPARENA_RESOLVE(a); + nn = &a->root; + +#ifdef notdef + printf("*** alloc %u\n", sz); + tdump(a->root); putchar('\n'); +#endif + + /* --- First, find a block which is big enough --- */ + +again: + n = *nn; + if (!n) { + v = a->ops->alloc(a, MPWS(sz + 1)); + v[0] = sz; + return (v + 1); + } + if (n->v[0] < sz) { + nn = &n->right; + goto again; + } + + /* --- Now try to find a smaller block which is suitable --- */ + + while (n->left && n->left->v[0] >= sz) { + nn = &n->left; + n = *nn; + } + + /* --- If the block we've got is still too large, start digging --- */ + + if (n->v[0] >= sz * 2) { + nn = &n->left; + goto again; + } + + /* --- I've now found a suitable block --- */ + + v = n->v; + + /* --- Remove this node from the tree --- */ + + if (!n->left) + *nn = n->right; + else if (!n->right) + *nn = n->left; + else { + mparena_node *left = n->left; + mparena_node *p = *nn = n->right; + while (p->left) + p = p->left; + p->left = left; + } + + /* --- Get rid of this node now --- */ + + DESTROY(n); + return (v + 1); +} + +/* --- @mpfree@ --- * + * + * Arguments: @mparena *a@ = pointer to arena block + * @mpw *v@ = pointer to allocated vector + * + * Returns: --- + * + * Use: Returns an MP vector to an arena. It doesn't have to be + * returned to the arena from which it was allocated. + */ + +void mpfree(mparena *a, mpw *v) +{ + mparena_node **nn, *n; + size_t sz = *--v; + + MPARENA_RESOLVE(a); + nn = &a->root; + + while (*nn) { + n = *nn; + if (n->v[0] > sz) + nn = &n->left; + else + nn = &n->right; + } + + n = CREATE(mparena_node); + n->left = n->right = 0; + n->v = v; + *nn = n; + +#ifdef notdef + printf("*** free %u\n", sz); + tdump(a->root); putchar('\n'); +#endif +} + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/mparena.h b/mparena.h new file mode 100644 index 0000000..61c735e --- /dev/null +++ b/mparena.h @@ -0,0 +1,167 @@ +/* -*-c-*- + * + * $Id: mparena.h,v 1.1 1999/11/17 18:02:16 mdw Exp $ + * + * Allocation and freeing of MP buffers + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: mparena.h,v $ + * Revision 1.1 1999/11/17 18:02:16 mdw + * New multiprecision integer arithmetic suite. + * + */ + +#ifndef MPARENA_H +#define MPARENA_H + +#ifdef __cplusplus + extern "C" { +#endif + +/*----- Header files ------------------------------------------------------*/ + +#ifndef MPW_H +# include "mpw.h" +#endif + +/*----- Data structures ---------------------------------------------------*/ + +/* --- @mparena_node@ --- * + * + * For internal use by the MP arena manager. The free blocks are held in a + * binary tree by size, held in the first digit of each vector. + */ + +typedef struct mparena_node { + struct mparena_node *left, *right; + mpw *v; +} mparena_node; + +/* --- @mparena@ --- * + * + * The actual arena. + */ + +typedef struct mparena { + mparena_node *root; + struct mparena_ops *ops; +} mparena; + +/* --- @mparena_ops@ --- * + * + * Operations required for an arena memory manager. The default manager just + * calls @xmalloc@ and @free@, although it's possible to envisage a more + * paranoid implementation which allocates locked memory pages. Switch them + * over with @mparena_setops@. It's usual to only do this when you've + * attached your extra state to the end of the @mparena@ structure. + */ + +typedef struct mparena_ops { + void *(*alloc)(mparena */*a*/, size_t /*sz*/); + void (*free)(mparena */*a*/, void */*p*/); +} mparena_ops; + +/*----- Magical constants -------------------------------------------------*/ + +#define MPARENA_GLOBAL ((mparena *)0) + +extern mparena_ops mparena_defaultops; + +/*----- Functions provided ------------------------------------------------*/ + +/* --- @mparena_create@ --- * + * + * Arguments: @mparena *a@ = pointer to arena block + * + * Returns: --- + * + * Use: Initializes an MP arena so that blocks can be allocated from + * it. + */ + +extern void mparena_create(mparena */*a*/); + +#define MPARENA_INIT { 0, &mparena_defaultops } + +/* --- @mparena_setops@ --- * + * + * Arguments: @mparena *a@ = pointer to arena block + * @mparena_ops *ops@ = pointer to operations block or null + * + * Returns: The previous operations block. + * + * Use: Sets or queries the operations attached to an arena. + */ + +extern mparena_ops *mparena_setops(mparena */*a*/, mparena_ops */*ops*/); + +/* --- @mparena_destroy@ --- * + * + * Arguments: @mparena *a@ = pointer to arena block + * + * Returns: --- + * + * Use: Frees an MP arena, and all the vectors held within it. The + * blocks which are currently allocated can be freed into some + * other arena. + */ + +extern void mparena_destroy(mparena */*a*/); + +/* --- @mpalloc@ --- * + * + * Arguments: @mparena *a@ = pointer to arena block + * @size_t sz@ = number of digits required + * + * Returns: Pointer to a suitably sized block. + * + * Use: Allocates a lump of data suitable for use as an array of MP + * digits. + */ + +extern mpw *mpalloc(mparena */*a*/, size_t /*sz*/); + +/* --- @mpfree@ --- * + * + * Arguments: @mparena *a@ = pointer to arena block + * @mpw *v@ = pointer to allocated vector + * + * Returns: --- + * + * Use: Returns an MP vector to an arena. It doesn't have to be + * returned to the arena from which it was allocated. + */ + +extern void mpfree(mparena */*a*/, mpw */*v*/); + +/*----- That's all, folks -------------------------------------------------*/ + +#ifdef __cplusplus + } +#endif + +#endif diff --git a/mpmont.c b/mpmont.c new file mode 100644 index 0000000..d924009 --- /dev/null +++ b/mpmont.c @@ -0,0 +1,419 @@ +/* -*-c-*- + * + * $Id: mpmont.c,v 1.1 1999/11/17 18:02:16 mdw Exp $ + * + * Montgomery reduction + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: mpmont.c,v $ + * Revision 1.1 1999/11/17 18:02:16 mdw + * New multiprecision integer arithmetic suite. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include "mp.h" +#include "mpmont.h" + +/*----- Main code ---------------------------------------------------------*/ + +/* --- @mpmont_create@ --- * + * + * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context + * @mp *m@ = modulus to use + * + * Returns: --- + * + * Use: Initializes a Montgomery reduction context ready for use. + */ + +void mpmont_create(mpmont *mm, mp *m) +{ + /* --- Take a copy of the modulus --- */ + + mp_shrink(m); + mm->m = MP_COPY(m); + + /* --- Find the magic value @mi@ --- * + * + * This is a slightly grungy way of solving the problem, but it does work. + */ + + { + mpw av[2] = { 0, 1 }; + mp a, b; + mp *i; + mpw mi; + + mp_build(&a, av, av + 2); + mp_build(&b, m->v, m->v + 1); + mp_gcd(0, 0, &i, &a, &b); + mi = i->v[0]; + if (!(i->f & MP_NEG)) + mi = MPW(-mi); + mm->mi = mi; + MP_DROP(i); + } + + /* --- Discover the values %$R \bmod m$% and %$R^2 \bmod m$% --- */ + + { + size_t l = MP_LEN(m); + mp *r = mp_create(l + 1); + + mm->shift = l * MPW_BITS; + MPX_ZERO(r->v, r->vl - 1); + r->vl[-1] = 1; + mm->r = mm->r2 = MP_NEW; + mp_div(0, &mm->r, r, m); + r = mp_sqr(r, mm->r); + mp_div(0, &mm->r2, r, m); + MP_DROP(r); + } +} + +/* --- @mpmont_destroy@ --- * + * + * Arguments: @mpmont *mm@ = pointer to a Montgomery reduction context + * + * Returns: --- + * + * Use: Disposes of a context when it's no longer of any use to + * anyone. + */ + +void mpmont_destroy(mpmont *mm) +{ + MP_DROP(mm->m); + MP_DROP(mm->r); + MP_DROP(mm->r2); +} + +/* --- @mpmont_reduce@ --- * + * + * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context + * @mp *d@ = destination + * @const mp *a@ = source, assumed positive + * + * Returns: Result, %$a R^{-1} \bmod m$%. + */ + +mp *mpmont_reduce(mpmont *mm, mp *d, const mp *a) +{ + mpw *dv, *dvl; + const mpw *mv, *mvl; + size_t n; + + /* --- Initial conditioning of the arguments --- */ + + n = MP_LEN(mm->m); + + if (d == a) + MP_MODIFY(d, 2 * n); + else { + MP_MODIFY(d, 2 * n); + memcpy(d->v, a->v, MPWS(MP_LEN(a))); + memset(d->v + MP_LEN(a), 0, MPWS(MP_LEN(d) - MP_LEN(a))); + } + + dv = d->v; dvl = d->vl; + mv = mm->m->v; mvl = mm->m->vl; + + /* --- Let's go to work --- */ + + while (n--) { + mpw u = MPW(*dv * mm->mi); + MPX_UMLAN(dv, dvl, mv, mvl, u); + dv++; + } + + /* --- Done --- */ + + memmove(d->v, dv, MPWS(dvl - dv)); + d->vl -= dv - d->v; + MP_SHRINK(d); + d->f = a->f & MP_BURN; + return (d); +} + +/* --- @mpmont_mul@ --- * + * + * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context + * @mp *d@ = destination + * @const mp *a, *b@ = sources, assumed positive + * + * Returns: Result, %$a b R^{-1} \bmod m$%. + */ + +mp *mpmont_mul(mpmont *mm, mp *d, const mp *a, const mp *b) +{ + mpw *dv, *dvl; + const mpw *av, *avl; + const mpw *bv, *bvl; + const mpw *mv, *mvl; + mpw y; + size_t n, i; + + /* --- Initial conditioning of the arguments --- */ + + if (MP_LEN(a) > MP_LEN(b)) { + const mp *t = a; a = b; b = t; + } + n = MP_LEN(mm->m); + + MP_MODIFY(d, 2 * n + 1); + dv = d->v; dvl = d->vl; + MPX_ZERO(dv, dvl); + av = a->v; avl = a->vl; + bv = b->v; bvl = b->vl; + mv = mm->m->v; mvl = mm->m->vl; + y = *bv; + + /* --- Montgomery multiplication phase --- */ + + i = 0; + while (i < n && av < avl) { + mpw x = *av++; + mpw u = MPW((*dv + x * y) * mm->mi); + MPX_UMLAN(dv, dvl, bv, bvl, x); + MPX_UMLAN(dv, dvl, mv, mvl, u); + dv++; + i++; + } + + /* --- Simpler Montgomery reduction phase --- */ + + while (i < n) { + mpw u = MPW(*dv * mm->mi); + MPX_UMLAN(dv, dvl, mv, mvl, u); + dv++; + i++; + } + + /* --- Done --- */ + + memmove(d->v, dv, MPWS(dvl - dv)); + d->vl -= dv - d->v; + MP_SHRINK(d); + d->f = (a->f | b->f) & MP_BURN; + return (d); +} + +/* --- @mpmont_exp@ --- * + * + * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context + * @const mp *a@ = base + * @const mp *e@ = exponent + * + * Returns: Result, %$a^e \bmod m$%. + */ + +mp *mpmont_exp(mpmont *mm, const mp *a, const mp *e) +{ + mpscan sc; + mp *ar = mpmont_mul(mm, MP_NEW, a, mm->r2); + mp *d = MP_COPY(mm->r); + + mp_scan(&sc, e); + + if (MP_STEP(&sc)) { + for (;;) { + mp *dd; + if (MP_BIT(&sc)) { + dd = mpmont_mul(mm, MP_NEW, d, ar); + MP_DROP(d); + d = dd; + } + if (!MP_STEP(&sc)) + break; + dd = mpmont_mul(mm, MP_NEW, ar, ar); + MP_DROP(ar); + ar = dd; + } + } + MP_DROP(ar); + return (mpmont_reduce(mm, d, d)); +} + +/*----- Test rig ----------------------------------------------------------*/ + +#ifdef TEST_RIG + +static int tcreate(dstr *v) +{ + mp *m = *(mp **)v[0].buf; + mp *mi = *(mp **)v[1].buf; + mp *r = *(mp **)v[2].buf; + mp *r2 = *(mp **)v[3].buf; + + mpmont mm; + int ok = 1; + + mpmont_create(&mm, m); + + if (mm.mi != mi->v[0]) { + fprintf(stderr, "\n*** bad mi: found %lu, expected %lu", + (unsigned long)mm.mi, (unsigned long)mi->v[0]); + fputs("\nm = ", stderr); mp_writefile(m, stderr, 10); + fputc('\n', stderr); + ok = 0; + } + + if (MP_CMP(mm.r, !=, r)) { + fputs("\n*** bad r", stderr); + fputs("\nm = ", stderr); mp_writefile(m, stderr, 10); + fputs("\nexpected ", stderr); mp_writefile(r, stderr, 10); + fputs("\n found ", stderr); mp_writefile(r, stderr, 10); + fputc('\n', stderr); + ok = 0; + } + + if (MP_CMP(mm.r2, !=, r2)) { + fputs("\n*** bad r2", stderr); + fputs("\nm = ", stderr); mp_writefile(m, stderr, 10); + fputs("\nexpected ", stderr); mp_writefile(r2, stderr, 10); + fputs("\n found ", stderr); mp_writefile(r2, stderr, 10); + fputc('\n', stderr); + ok = 0; + } + + MP_DROP(m); + MP_DROP(mi); + MP_DROP(r); + MP_DROP(r2); + mpmont_destroy(&mm); + return (ok); +} + +static int tmul(dstr *v) +{ + mp *m = *(mp **)v[0].buf; + mp *a = *(mp **)v[1].buf; + mp *b = *(mp **)v[2].buf; + mp *r = *(mp **)v[3].buf; + mp *mr, *qr; + int ok = 1; + + mpmont mm; + mpmont_create(&mm, m); + + { + mp *ar = mpmont_mul(&mm, MP_NEW, a, mm.r2); + mp *br = mpmont_mul(&mm, MP_NEW, b, mm.r2); + mr = mpmont_mul(&mm, MP_NEW, ar, br); + mr = mpmont_reduce(&mm, mr, mr); + MP_DROP(ar); MP_DROP(br); + } + + qr = mp_mul(MP_NEW, a, b); + mp_div(0, &qr, qr, m); + + if (MP_CMP(qr, !=, r)) { + fputs("\n*** classical modmul failed", stderr); + fputs("\n m = ", stderr); mp_writefile(m, stderr, 10); + fputs("\n a = ", stderr); mp_writefile(a, stderr, 10); + fputs("\n b = ", stderr); mp_writefile(b, stderr, 10); + fputs("\n r = ", stderr); mp_writefile(r, stderr, 10); + fputs("\nqr = ", stderr); mp_writefile(qr, stderr, 10); + fputc('\n', stderr); + ok = 0; + } + + if (MP_CMP(mr, !=, r)) { + fputs("\n*** montgomery modmul failed", stderr); + fputs("\n m = ", stderr); mp_writefile(m, stderr, 10); + fputs("\n a = ", stderr); mp_writefile(a, stderr, 10); + fputs("\n b = ", stderr); mp_writefile(b, stderr, 10); + fputs("\n r = ", stderr); mp_writefile(r, stderr, 10); + fputs("\nmr = ", stderr); mp_writefile(mr, stderr, 10); + fputc('\n', stderr); + ok = 0; + } + + MP_DROP(m); + MP_DROP(a); + MP_DROP(b); + MP_DROP(r); + MP_DROP(mr); + MP_DROP(qr); + mpmont_destroy(&mm); + return ok; +} + +static int texp(dstr *v) +{ + mp *m = *(mp **)v[0].buf; + mp *a = *(mp **)v[1].buf; + mp *b = *(mp **)v[2].buf; + mp *r = *(mp **)v[3].buf; + mp *mr; + int ok = 1; + + mpmont mm; + mpmont_create(&mm, m); + + mr = mpmont_exp(&mm, a, b); + + if (MP_CMP(mr, !=, r)) { + fputs("\n*** montgomery modexp failed", stderr); + fputs("\n m = ", stderr); mp_writefile(m, stderr, 10); + fputs("\n a = ", stderr); mp_writefile(a, stderr, 10); + fputs("\n e = ", stderr); mp_writefile(b, stderr, 10); + fputs("\n r = ", stderr); mp_writefile(r, stderr, 10); + fputs("\nmr = ", stderr); mp_writefile(mr, stderr, 10); + fputc('\n', stderr); + ok = 0; + } + + MP_DROP(m); + MP_DROP(a); + MP_DROP(b); + MP_DROP(r); + MP_DROP(mr); + mpmont_destroy(&mm); + return ok; +} + + +static test_chunk tests[] = { + { "create", tcreate, { &type_mp, &type_mp, &type_mp, &type_mp } }, + { "mul", tmul, { &type_mp, &type_mp, &type_mp, &type_mp } }, + { "exp", texp, { &type_mp, &type_mp, &type_mp, &type_mp } }, + { 0, 0, { 0 } }, +}; + +int main(int argc, char *argv[]) +{ + sub_init(); + test_run(argc, argv, tests, SRCDIR "/tests/mpmont"); + return (0); +} + +#endif + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/mpmont.h b/mpmont.h new file mode 100644 index 0000000..42aa17c --- /dev/null +++ b/mpmont.h @@ -0,0 +1,114 @@ +/* -*-c-*- + * + * $Id: mpmont.h,v 1.1 1999/11/17 18:02:16 mdw Exp $ + * + * Montgomery reduction + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: mpmont.h,v $ + * Revision 1.1 1999/11/17 18:02:16 mdw + * New multiprecision integer arithmetic suite. + * + */ + +#ifndef MPMONT_H +#define MPMONT_H + +#ifdef __cplusplus + extern "C" { +#endif + +/*----- Header files ------------------------------------------------------*/ + +#ifndef MP_H +# include "mp.h" +#endif + +/*----- What's going on here? ---------------------------------------------* + * + * Given a little bit of precomputation, Montgomery reduction enables modular + * reductions of products to be calculated rather rapidly, without recourse + * to annoying things like division. + * + * Before starting, you need to do a little work. In particular, the + * following things need to be worked out: + * + * * %$m$%, which is the modulus you'll be working with. + * + * * %$b$%, the radix of the number system you're in (here, it's + * @MPW_MAX + 1@). + * + * * %$-m^{-1} \bmod b$%, a useful number for the reduction step. (This + * means that the modulus mustn't be even. This shouldn't be a problem.) + * + * * %$R = b^n > m > b^{n - 1}$%, or at least %$\log_2 R$%. + * + * * %$R \bmod m$% and %$R^2 \bmod m$%, which are useful when doing + * calculations such as exponentiation. + * + * The result of a Montgomery reduction of %$x$% is %$x R^{-1} \bmod m$%, + * which doesn't look ever-so useful. The trick is to initially apply a + * factor of %$R$% to all of your numbers so that when you multiply and + * perform a Montgomery reduction you get %$(xR \cdot yR)R^{-1} \bmod m$%, + * which is just %$xyR \bmod m$%. Thanks to distributivity, even additions + * and subtractions can be performed on numbers in this form -- the extra + * factor of %$R$% just runs through all the calculations until it's finally + * stripped out by a final reduction operation. + */ + +/*----- Data structures ---------------------------------------------------*/ + +/* --- A Montgomery reduction context --- */ + +typedef struct mpmont { + mp *m; /* Modulus */ + mpw mi; /* %$-m^{-1} \bmod b$% */ + size_t shift; /* %$\log_2 R$% */ + mp *r, *r2; /* %$R \bmod m$%, %$R^2 \bmod m$% */ +} mpmont; + +/*----- Functions provided ------------------------------------------------*/ + +/* --- @mpmont_create@ --- * + * + * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context + * @mp *m@ = modulus to use + * + * Returns: --- + * + * Use: Initializes a Montgomery reduction context ready for use. + */ + +extern void mpmont_create(mpmont */*mm*/, mp */*m*/); + +/*----- That's all, folks -------------------------------------------------*/ + +#ifdef __cplusplus + } +#endif + +#endif diff --git a/mptext-dstr.c b/mptext-dstr.c new file mode 100644 index 0000000..df473f9 --- /dev/null +++ b/mptext-dstr.c @@ -0,0 +1,96 @@ +/* -*-c-*- + * + * $Id: mptext-dstr.c,v 1.1 1999/11/17 18:02:16 mdw Exp $ + * + * Reading and writing large integers on strings + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: mptext-dstr.c,v $ + * Revision 1.1 1999/11/17 18:02:16 mdw + * New multiprecision integer arithmetic suite. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include + +#include "mptext.h" + +/*----- Main code ---------------------------------------------------------*/ + +/* --- Operations table --- */ + +static int get(void *p) +{ + mptext_dstrctx *c = p; + if (c->i >= c->d->len) + return (EOF); + return (c->d->buf[c->i++]); +} + +static void unget(int ch, void *p) +{ + mptext_dstrctx *c = p; + if (ch == EOF || c->i == 0) + return; + c->i--; +} + +static int put(char *s, size_t sz, void *p) +{ + mptext_dstrctx *c = p; + DPUTM(c->d, s, sz); + return (0); +} + +const mptext_ops mptext_dstrops = { get, unget, put }; + +/* --- Convenience functions --- */ + +mp *mp_readdstr(mp *m, dstr *d, size_t *off, int radix) +{ + mptext_dstrctx c; + c.d = d; + c.i = off ? *off : 0; + m = mp_read(m, radix, &mptext_dstrops, &c); + if (off) + *off = c.i; + return (m); +} + +int mp_writedstr(mp *m, dstr *d, int radix) +{ + mptext_dstrctx c; + int rc; + c.d = d; + rc = mp_write(m, radix, &mptext_dstrops, &c); + DPUTZ(d); + return (rc); +} + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/mptext-file.c b/mptext-file.c new file mode 100644 index 0000000..ed6f700 --- /dev/null +++ b/mptext-file.c @@ -0,0 +1,72 @@ +/* -*-c-*- + * + * $Id: mptext-file.c,v 1.1 1999/11/17 18:02:16 mdw Exp $ + * + * Reading and writing large integers on files + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: mptext-file.c,v $ + * Revision 1.1 1999/11/17 18:02:16 mdw + * New multiprecision integer arithmetic suite. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include + +#include "mptext.h" + +/*----- Main code ---------------------------------------------------------*/ + +/* --- Operations table --- */ + +static int get(void *p) { FILE *fp = p; return getc(fp); } + +static void unget(int ch, void *p) { FILE *fp = p; ungetc(ch, fp); } + +static int put(char *s, size_t sz, void *p) +{ + FILE *fp = p; + return (fwrite(s, 1, sz, fp) != sz); +} + +const mptext_ops mptext_fileops = { get, unget, put }; + +/* --- Convenience functions --- */ + +mp *mp_readfile(mp *m, FILE *fp, int radix) +{ + return mp_read(m, radix, &mptext_fileops, fp); +} + +int mp_writefile(mp *m, FILE *fp, int radix) +{ + return mp_write(m, radix, &mptext_fileops, fp); +} + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/mptext-string.c b/mptext-string.c new file mode 100644 index 0000000..962fc96 --- /dev/null +++ b/mptext-string.c @@ -0,0 +1,106 @@ +/* -*-c-*- + * + * $Id: mptext-string.c,v 1.1 1999/11/17 18:02:16 mdw Exp $ + * + * Reading and writing large integers on strings + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: mptext-string.c,v $ + * Revision 1.1 1999/11/17 18:02:16 mdw + * New multiprecision integer arithmetic suite. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include + +#include "mptext.h" + +/*----- Main code ---------------------------------------------------------*/ + +/* --- Operations table --- */ + +static int get(void *p) +{ + mptext_stringctx *c = p; + if (c->buf >= c->lim) + return (EOF); + return (*c->buf++); +} + +static void unget(int ch, void *p) +{ + mptext_stringctx *c = p; + if (ch != EOF) + c->buf--; +} + +static int put(char *s, size_t sz, void *p) +{ + mptext_stringctx *c = p; + int rc = 0; + if (sz > c->lim - c->buf) { + sz = c->lim - c->buf; + rc = EOF; + } + if (sz) { + memcpy(c->buf, s, sz); + c->buf += sz; + } + return (rc); +} + +const mptext_ops mptext_stringops = { get, unget, put }; + +/* --- Convenience functions --- */ + +mp *mp_readstring(mp *m, const char *p, char **end, int radix) +{ + mptext_stringctx c; + c.buf = (char *)p; + c.lim = c.buf + strlen(p); + m = mp_read(m, radix, &mptext_stringops, &c); + if (end) + *end = c.buf; + return (m); +} + +int mp_writestring(mp *m, char *p, size_t sz, int radix) +{ + mptext_stringctx c; + int rc; + if (!sz) + return (0); + c.buf = p; + c.lim = p + sz - 1; + rc = mp_write(m, radix, &mptext_stringops, &c); + *c.buf = 0; + return (rc); +} + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/mptext.c b/mptext.c new file mode 100644 index 0000000..dd78bde --- /dev/null +++ b/mptext.c @@ -0,0 +1,314 @@ +/* -*-c-*- + * + * $Id: mptext.c,v 1.1 1999/11/17 18:02:16 mdw Exp $ + * + * Textual representation of multiprecision numbers + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: mptext.c,v $ + * Revision 1.1 1999/11/17 18:02:16 mdw + * New multiprecision integer arithmetic suite. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include +#include + +#include + +#include "mp.h" +#include "mptext.h" + +/*----- Data structures ---------------------------------------------------*/ + +#ifndef CHARV + DA_DECL(charv, char); +# define CHARV +#endif + +/*----- Main code ---------------------------------------------------------*/ + +/* --- @mp_read@ --- * + * + * Arguments: @mp *m@ = destination multiprecision number + * @int radix@ = base to assume for data (or zero to guess) + * @const mptext_ops *ops@ = pointer to operations block + * @void *p@ = data for the operations block + * + * Returns: The integer read, or zero if it didn't work. + * + * Use: Reads an integer from some source. If the @radix@ is + * specified, the number is assumed to be given in that radix, + * with the letters `a' (either upper- or lower-case) upwards + * standing for digits greater than 9. Otherwise, base 10 is + * assumed unless the number starts with `0' (octal), `0x' (hex) + * or `nnn_' (base `nnn'). An arbitrary amount of whitespace + * before the number is ignored. + */ + +mp *mp_read(mp *m, int radix, const mptext_ops *ops, void *p) +{ + int r; + int ch; + unsigned f = 0; + + enum { + f_neg = 1u, + f_ok = 2u + }; + + /* --- Initialize the destination number --- */ + + MP_MODIFY(m, 16); + m->vl = m->v; + m->f &= ~MP_UNDEF; + + /* --- Read an initial character --- */ + + ch = ops->get(p); + while (isspace(ch)) + ch = ops->get(p); + + /* --- Handle an initial sign --- */ + + if (ch == '-') { + f |= f_neg; + ch = ops->get(p); + while (isspace(ch)) + ch = ops->get(p); + } + + /* --- If the radix is zero, look for leading zeros --- */ + + if (radix) + r = -1; + else if (ch != '0') { + radix = 10; + r = 0; + } else { + ch = ops->get(p); + if (ch == 'x') { + ch = ops->get(p); + radix = 16; + } else { + radix = 8; + f |= f_ok; + } + r = -1; + } + + /* --- Time to start --- */ + + for (;; ch = ops->get(p)) { + int x; + + /* --- An underscore indicates a numbered base --- */ + + if (ch == '_' && r > 0 && r <= 36) { + radix = r; + m->vl = m->v; + r = -1; + f &= ~f_ok; + continue; + } + + /* --- Check that the character is a digit and in range --- */ + + if (!isalnum(ch)) + break; + if (ch >= '0' && ch <= '9') + x = ch - '0'; + else { + ch = tolower(ch); + if (ch >= 'a' && ch <= 'z') /* ASCII dependent! */ + x = ch - 'a' + 10; + else + break; + } + + /* --- Sort out what to do with the character --- */ + + if (x >= 10 && r >= 0) + r = -1; + if (x >= radix) + break; + + if (r >= 0) + r = r * 10 + x; + + /* --- Stick the character on the end of my integer --- */ + + MP_ENSURE(m, MP_LEN(m) + 1); + MPX_UMULN(m->v, m->vl, m->v, m->vl - 1, radix); + MPX_UADDN(m->v, m->vl, x); + MP_SHRINK(m); + f |= f_ok; + } + + ops->unget(ch, p); + + /* --- Bail out if the number was bad --- */ + + if (!(f & f_ok)) { + MP_DROP(m); + return (0); + } + + /* --- Set the sign and return --- */ + + m->f = 0; + if (f & f_neg) + m->f |= MP_NEG; + return (m); +} + +/* --- @mp_write@ --- * + * + * Arguments: @mp *m@ = pointer to a multi-precision integer + * @int radix@ = radix to use when writing the number out + * @const mptext_ops *ops@ = pointer to an operations block + * @void *p@ = data for the operations block + * + * Returns: Zero if it worked, nonzero otherwise. + * + * Use: Writes a large integer in textual form. + */ + +int mp_write(mp *m, int radix, const mptext_ops *ops, void *p) +{ + charv v = DA_INIT; + mpw b = radix; + mp bb; + + /* --- Set various things up --- */ + + m = MP_COPY(m); + mp_build(&bb, &b, &b + 1); + + /* --- If the number is negative, sort that out --- */ + + if (m->f & MP_NEG) { + if (ops->put("-", 1, p)) + return (EOF); + MP_SPLIT(m); + m->f &= ~MP_NEG; + } + + /* --- Write digits to a temporary array --- */ + + do { + mp *q = MP_NEW, *r = MP_NEW; + int ch; + mpw x; + + mp_div(&q, &r, m, &bb); + x = r->v[0]; + if (x < 10) + ch = '0' + x; + else + ch = 'a' + x - 10; + DA_UNSHIFT(&v, ch); + MP_DROP(m); + MP_DROP(r); + m = q; + } while (MP_CMP(m, !=, MP_ZERO)); + + /* --- Finished that --- */ + + MP_DROP(m); + + { + int rc = ops->put(DA(&v), DA_LEN(&v), p); + DA_DESTROY(&v); + return (rc ? EOF : 0); + } +} + +/*----- Test rig ----------------------------------------------------------*/ + +#ifdef TEST_RIG + +#include + +static int verify(dstr *v) +{ + int ok = 1; + int ib = *(int *)v[0].buf, ob = *(int *)v[2].buf; + dstr d = DSTR_INIT; + mp *m = mp_readdstr(MP_NEW, &v[1], 0, ib); + if (m) { + if (!ob) { + fprintf(stderr, "*** unexpected successful parse\n" + "*** input [%i] = %s\n", + ib, v[1].buf); + mp_writedstr(m, &d, 10); + fprintf(stderr, "*** (value = %s)\n", d.buf); + ok = 0; + } else { + mp_writedstr(m, &d, ob); + if (d.len != v[3].len || memcmp(d.buf, v[3].buf, d.len) != 0) { + fprintf(stderr, "*** failed read or write\n" + "*** input [%i] = %s\n" + "*** output [%i] = %s\n" + "*** expected [%i] = %s\n", + ib, v[1].buf, ob, d.buf, ob, v[3].buf); + ok = 0; + } + } + mp_drop(m); + } else { + if (ob) { + fprintf(stderr, "*** unexpected parse failure\n" + "*** input [%i] = %s\n" + "*** expected [%i] = %s\n", + ib, v[1].buf, ob, v[3].buf); + ok = 0; + } + } + + dstr_destroy(&d); + return (ok); +} + +static test_chunk tests[] = { + { "mptext", verify, + { &type_int, &type_string, &type_int, &type_string, 0 } }, + { 0, 0, { 0 } } +}; + +int main(int argc, char *argv[]) +{ + sub_init(); + test_run(argc, argv, tests, SRCDIR "/tests/mptext"); + return (0); +} + +#endif + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/mptext.h b/mptext.h new file mode 100644 index 0000000..667c098 --- /dev/null +++ b/mptext.h @@ -0,0 +1,160 @@ +/* -*-c-*- + * + * $Id: mptext.h,v 1.1 1999/11/17 18:02:16 mdw Exp $ + * + * Textual representation of multiprecision numbers + * + * (c) 1999 Straylight/Edgeware + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This file is part of Catacomb. + * + * Catacomb is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as + * published by the Free Software Foundation; either version 2 of the + * License, or (at your option) any later version. + * + * Catacomb is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with Catacomb; if not, write to the Free + * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, + * MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: mptext.h,v $ + * Revision 1.1 1999/11/17 18:02:16 mdw + * New multiprecision integer arithmetic suite. + * + */ + +#ifndef MPTEXT_H +#define MPTEXT_H + +#ifdef __cplusplus + extern "C" { +#endif + +/*----- Header files ------------------------------------------------------*/ + +#ifndef MP_H +# include "mp.h" +#endif + +/*----- Data structures ---------------------------------------------------*/ + +typedef struct mptext_ops { + int (*get)(void */*p*/); + void (*unget)(int /*ch*/, void */*p*/); + int (*put)(char */*s*/, size_t /*len*/, void */*p*/); +} mptext_ops; + +/*----- Functions provided ------------------------------------------------*/ + +/* --- @mp_read@ --- * + * + * Arguments: @mp *m@ = destination multiprecision number + * @int radix@ = base to assume for data (or zero to guess) + * @const mptext_ops *ops@ = pointer to operations block + * @void *p@ = data for the operations block + * + * Returns: The integer read, or zero if it didn't work. + * + * Use: Reads an integer from some source. If the @radix@ is + * specified, the number is assumed to be given in that radix, + * with the letters `a' (either upper- or lower-case) upwards + * standing for digits greater than 9. Otherwise, base 10 is + * assumed unless the number starts with `0' (octal), `0x' (hex) + * or `nnn_' (base `nnn'). An arbitrary amount of whitespace + * before the number is ignored. + */ + +extern mp *mp_read(mp */*m*/, int /*radix*/, + const mptext_ops */*ops*/, void */*p*/); + +/* --- @mp_write@ --- * + * + * Arguments: @mp *m@ = pointer to a multi-precision integer + * @int radix@ = radix to use when writing the number out + * @const mptext_ops *ops@ = pointer to an operations block + * @void *p@ = data for the operations block + * + * Returns: Zero if it worked, nonzero otherwise. + * + * Use: Writes a large integer in textual form. + */ + +extern int mp_write(mp */*m*/, int /*radix*/, + const mptext_ops */*ops*/, void */*p*/); + +/*----- File I/O ----------------------------------------------------------*/ + +#include + +/* --- Operations table --- * + * + * The @mptext_fileops@ expect the pointer argument to be a @FILE *@. + */ + +extern const mptext_ops mptext_fileops; + +/* --- Convenience functions --- */ + +extern mp *mp_readfile(mp */*m*/, FILE */*fp*/, int /*radix*/); +extern int mp_writefile(mp */*m*/, FILE */*fp*/, int /*radix*/); + +/*----- String I/O --------------------------------------------------------*/ + +/* --- Context format --- */ + +typedef struct mptext_stringctx { + char *buf; + char *lim; +} mptext_stringctx; + +/* --- Operations table --- */ + +extern const mptext_ops mptext_stringops; + +/* --- Convenience functions --- */ + +extern mp *mp_readstring(mp */*m*/, const char */*p*/, char **/*end*/, + int /*radix*/); +extern int mp_writestring(mp */*m*/, char */*p*/, size_t /*sz*/, + int /*radix*/); + +/*----- Dynamic string I/O ------------------------------------------------*/ + +#include + +/* --- Context format --- */ + +typedef struct mptext_dstrctx { + dstr *d; + size_t i; +} mptext_dstrctx; + +/* --- Operations table --- */ + +extern const mptext_ops mptext_dstrops; + +/* --- Convenience functions --- */ + +extern mp *mp_readdstr(mp */*m*/, dstr */*d*/, size_t */*off*/, + int /*radix*/); +extern int mp_writedstr(mp */*m*/, dstr */*d*/, int /*radix*/); + +/*----- That's all, folks -------------------------------------------------*/ + +#ifdef __cplusplus + } +#endif + +#endif diff --git a/tests/mp b/tests/mp new file mode 100644 index 0000000..27854a2 --- /dev/null +++ b/tests/mp @@ -0,0 +1,57 @@ +# Test vectors for MP functions +# +# $Id: mp,v 1.1 1999/11/17 18:02:17 mdw Exp $ + +add { + 5 4 9; 5 -4 1; -5 4 -1; -5 -4 -9; + 0xffffffff 1 0x100000000; +} + +sub { + 5 4 1; 5 -4 9; -5 4 -9; -5 -4 -1; + 4 5 -1; 4 -5 9; -4 5 -9; -4 -5 1; +} + +mul { + 5 4 20; -5 4 -20; 5 -4 -20; -5 -4 20; + 0x10000 0x10000 0x100000000; +} + +div { + 9 4 2 1; -9 4 -3 3; 9 -4 -3 -3; -9 -4 2 -1; +} + +gcd { + 16 12 4 -2 3; + 12 16 4 -1 1; + 693 609 21 -181 206; + 4398082908043 90980984098081324 1 -32483863573352089 1570292150447; + + 829561629303257626084392170900075 32498098450983560651904114638965 + 5 -22841190347053190672253237276815 583054885752979049202923618992482; + + 5509672937670943767152343650729669537671508 + 398326674296699796695672966992514673531 + 17 + -191606556147997561067126486929677861359 + 2650310725368604614586643627755316700713319; + + 324098408098290809832490802984098208098324 + 23430980840982340982098409823089098443 + 1 + -4158709420138833210339208344965073815 + 57523460582278135926717203882531035926727; + + # --- RSA test --- + # + # The first number is (p - 1)(q - 1) from `mpmont'. The second is a + # random number (it's actually prime, but that doesn't matter) which I + # can use as an RSA encryption exponent. The last is the partner + # decryption exponent, produced using the extended GCD algorithm. + + 665251164384574309450646977867043764321191240895546832784045453360 + 5945908509680983480596809586040589085680968709809890671 + 1 + -4601007896041464028712478963832994007038251361995647370 + 514778499400157641662814932021958856708417966520837469125919104431; +} diff --git a/tests/mpmont b/tests/mpmont new file mode 100644 index 0000000..1cb2a70 --- /dev/null +++ b/tests/mpmont @@ -0,0 +1,53 @@ +# Test vectors for Montgomery reduction +# +# $Id: mpmont,v 1.1 1999/11/17 18:02:17 mdw Exp $ + +create { + 340809809850981098423498794792349 # m + 266454859 # -m^{-1} mod b + 130655606683780235388773757767708 # R mod m + 237786678640282040194246459306177; # R^2 mod m +} + +mul { + 43289823545 + 234324324 + 6456542564 + 10807149256; +} + +exp { + 4325987397987458979875737589783 + 435365332435654643667 + 8745435676786567758678547 + 2439674515119108242643169132064; + + # --- Quick RSA test --- + + 905609324890967090294090970600361 # This is p + 3 + 905609324890967090294090970600360 # This is (p - 1) + 1; # Fermat test: p is prime + + 734589569806680985408670989082927 # This is q + 5 + 734589569806680985408670989082926 # And this is (q - 1) + 1; # Fermat again: q is prime + + # --- Encrypt a message --- + # + # The public and private exponents are from the GCD test. The message + # is just obvious. The modulus is the product of the two primes above. + + 665251164384574309450646977867045404520085938543622535546005136647 + 123456789012345678901234567890123456789012345678901234567890 + 5945908509680983480596809586040589085680968709809890671 + 25906467774034212974484417859588980567136610347807401817990462701; + + # --- And decrypt it again --- + + 665251164384574309450646977867045404520085938543622535546005136647 + 25906467774034212974484417859588980567136610347807401817990462701 + 514778499400157641662814932021958856708417966520837469125919104431 + 123456789012345678901234567890123456789012345678901234567890; +} diff --git a/tests/mptext b/tests/mptext new file mode 100644 index 0000000..d77582b --- /dev/null +++ b/tests/mptext @@ -0,0 +1,30 @@ +# Test vectors for MP textual I/O +# +# $Id: mptext,v 1.1 1999/11/17 18:02:17 mdw Exp $ + +mptext { + + # --- Perfectly valid things --- + + 10 0 10 0; + 0 0 10 0; + 10 52 10 52; + 10 654365464655464577673765769678 10 654365464655464577673765769678; + 10 654365464655464577673765769678 16 8425e6d06f272b9a2d73ed1ce; + 16 8425E6D06F272B9A2D73ED1CE 10 654365464655464577673765769678; + 0 654365464655464577673765769678 16 8425e6d06f272b9a2d73ed1ce; + 0 16_8425E6D06F272B9A2D73ED1CE 10 654365464655464577673765769678; + 0 "-0x8425E6D06F272B9A2D73ED1CE" 10 "-654365464655464577673765769678"; + 8 "-366570443501403714657464766613" 10 "-596569802840985608098409867"; + 0 0366570443501403714657464766613 10 596569802840985608098409867; + + + # --- Bogus things --- + + 10 "" 0 0; # Empty string fails + 10 foo 0 0; # Non-numeric character + 10 134f 10 134; # Stop parsing when reaching `f' + 4 12345 10 27; # Stop parsing when reaching `4' + 0 37_ 10 37; # 37 is an invalid base, so stop at `_' + 0 36_ 0 0; # 36 is a valid base, so restart and fail +}