--- /dev/null
+/* -*-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 -------------------------------------------------*/
--- /dev/null
+/* -*-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 -------------------------------------------------*/
--- /dev/null
+/* -*-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 -------------------------------------------------*/
--- /dev/null
+/* -*-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 -------------------------------------------------*/
--- /dev/null
+/* -*-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 -------------------------------------------------*/
--- /dev/null
+/* -*-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 -------------------------------------------------*/
--- /dev/null
+/* -*-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 <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <mLib/dstr.h>
+#include <mLib/report.h>
+#include <mLib/sub.h>
+#include <mLib/testrig.h>
+
+#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 -------------------------------------------------*/
/* -*-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
*/
/*----- 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.
*
*/
/*----- Header files ------------------------------------------------------*/
-#include <stdlib.h>
+#include <assert.h>
#include <string.h>
-#ifndef MPTYPES_H
-# include "mptypes.h"
+#include <mLib/sub.h>
+
+#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 <mLib/testrig.h>
+
+#ifndef MPTEXT_H
+# include "mptext.h"
+#endif
+
+extern const test_type type_mp;
/*----- That's all, folks -------------------------------------------------*/
--- /dev/null
+/* -*-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
--- /dev/null
+/* -*-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 <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <mLib/alloc.h>
+#include <mLib/sub.h>
+
+#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 -------------------------------------------------*/
--- /dev/null
+/* -*-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
--- /dev/null
+/* -*-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 -------------------------------------------------*/
--- /dev/null
+/* -*-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
--- /dev/null
+/* -*-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 <mLib/dstr.h>
+
+#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 -------------------------------------------------*/
--- /dev/null
+/* -*-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 <stdio.h>
+
+#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 -------------------------------------------------*/
--- /dev/null
+/* -*-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 <string.h>
+
+#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 -------------------------------------------------*/
--- /dev/null
+/* -*-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 <ctype.h>
+#include <stdio.h>
+
+#include <mLib/darray.h>
+
+#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 <mLib/testrig.h>
+
+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 -------------------------------------------------*/
--- /dev/null
+/* -*-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 <stdio.h>
+
+/* --- 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 <mLib/dstr.h>
+
+/* --- 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
--- /dev/null
+# 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;
+}
--- /dev/null
+# 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;
+}
--- /dev/null
+# 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
+}
projects / u / mdw / catacomb / commitdiff