From ae747c9bab929bc794ddad9f59c381657d347d1f Mon Sep 17 00:00:00 2001 From: mdw Date: Sun, 8 Oct 2000 15:49:37 +0000 Subject: [PATCH] First glimmerings of binary polynomial arithmetic. --- gfx-kmul.c | 260 +++++++++++++++++++++++++++++++++++++ gfx-sqr-mktab.c | 99 ++++++++++++++ gfx-sqr.c | 222 +++++++++++++++++++++++++++++++ gfx.c | 397 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ gfx.h | 170 ++++++++++++++++++++++++ 5 files changed, 1148 insertions(+) create mode 100644 gfx-kmul.c create mode 100644 gfx-sqr-mktab.c create mode 100644 gfx-sqr.c create mode 100644 gfx.c create mode 100644 gfx.h diff --git a/gfx-kmul.c b/gfx-kmul.c new file mode 100644 index 0000000..dc2e524 --- /dev/null +++ b/gfx-kmul.c @@ -0,0 +1,260 @@ +/* -*-c-*- + * + * $Id: gfx-kmul.c,v 1.1 2000/10/08 15:49:37 mdw Exp $ + * + * Karatsuba's multiplication algorithm on binary polynomials + * + * (c) 2000 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: gfx-kmul.c,v $ + * Revision 1.1 2000/10/08 15:49:37 mdw + * First glimmerings of binary polynomial arithmetic. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include +#include + +#include "gfx.h" +#include "karatsuba.h" + +/*----- Tweakables --------------------------------------------------------*/ + +#ifdef TEST_RIG +# undef GFK_THRESH +# define GFK_THRESH 1 +#endif + +/*----- Main code ---------------------------------------------------------*/ + +/* --- @gfx_kmul@ --- * + * + * Arguments: @mpw *dv, *dvl@ = pointer to destination buffer + * @const mpw *av, *avl@ = pointer to first argument + * @const mpw *bv, *bvl@ = pointer to second argument + * @mpw *sv, *svl@ = pointer to scratch workspace + * + * Returns: --- + * + * Use: Multiplies two binary polynomials using Karatsuba's + * algorithm. This is rather faster than traditional long + * multiplication (e.g., @gfx_umul@) on polynomials with large + * degree, although more expensive on small ones. + * + * The destination must be twice as large as the larger + * argument. The scratch space must be twice as large as the + * larger argument. + */ + +void gfx_kmul(mpw *dv, mpw *dvl, + const mpw *av, const mpw *avl, + const mpw *bv, const mpw *bvl, + mpw *sv, mpw *svl) +{ + const mpw *avm, *bvm; + size_t m; + + /* --- Dispose of easy cases to @mpx_umul@ --- * + * + * Karatsuba is only a win on large numbers, because of all the + * recursiveness and bookkeeping. The recursive calls make a quick check + * to see whether to bottom out to @gfx_umul@ which should help quite a + * lot, but sometimes the only way to know is to make sure... + */ + + MPX_SHRINK(av, avl); + MPX_SHRINK(bv, bvl); + + if (avl - av <= GFK_THRESH || bvl - bv <= GFK_THRESH) { + gfx_mul(dv, dvl, av, avl, bv, bvl); + return; + } + + /* --- How the algorithm works --- * + * + * Let %$A = xb + y$% and %$B = ub + v$%. Then, simply by expanding, + * %$AB = x u b^2 + b(x v + y u) + y v$%. That's not helped any, because + * I've got four multiplications, each four times easier than the one I + * started with. However, note that I can rewrite the coefficient of %$b$% + * as %$xv + yu = (x + y)(u + v) - xu - yv$%. The terms %$xu$% and %$yv$% + * I've already calculated, and that leaves only one more multiplication to + * do. So now I have three multiplications, each four times easier, and + * that's a win. + */ + + /* --- First things --- * + * + * Sort out where to break the factors in half. I'll choose the midpoint + * of the larger one, since this minimizes the amount of work I have to do + * most effectively. + */ + + if (avl - av > bvl - bv) { + m = (avl - av + 1) >> 1; + avm = av + m; + if (bvl - bv > m) + bvm = bv + m; + else + bvm = bvl; + } else { + m = (bvl - bv + 1) >> 1; + bvm = bv + m; + if (avl - av > m) + avm = av + m; + else + avm = avl; + } + + assert(((void)"Destination too small for Karatsuba gf-multiply", + dvl - dv >= 4 * m)); + assert(((void)"Not enough workspace for Karatsuba gf-multiply", + svl - sv >= 4 * m)); + + /* --- Sort out the middle term --- */ + + { + mpw *bsv = sv + m, *ssv = bsv + m; + mpw *rdv = dv + m, *rdvl = rdv + 2 * m; + + UXOR2(sv, bsv, av, avm, avm, avl); + UXOR2(bsv, ssv, bv, bvm, bvm, bvl); + if (m > GFK_THRESH) + gfx_kmul(rdv, rdvl, sv, bsv, bsv, ssv, ssv, svl); + else + gfx_mul(rdv, rdvl, sv, bsv, bsv, ssv); + } + + /* --- Sort out the other two terms --- */ + + { + mpw *svm = sv + m, *ssv = svm + m; + mpw *tdv = dv + m; + mpw *rdv = tdv + m; + + if (avl == avm || bvl == bvm) + MPX_ZERO(rdv + m, dvl); + else { + if (m > GFK_THRESH) + gfx_kmul(sv, ssv, avm, avl, bvm, bvl, ssv, svl); + else + gfx_mul(sv, ssv, avm, avl, bvm, bvl); + MPX_COPY(rdv + m, dvl, svm, ssv); + UXOR(rdv, sv, svm); + UXOR(tdv, sv, ssv); + } + + if (m > GFK_THRESH) + gfx_kmul(sv, ssv, av, avm, bv, bvm, ssv, svl); + else + gfx_mul(sv, ssv, av, avm, bv, bvm); + MPX_COPY(dv, tdv, sv, svm); + UXOR(tdv, sv, ssv); + UXOR(tdv, svm, ssv); + } +} + +/*----- Test rig ----------------------------------------------------------*/ + +#ifdef TEST_RIG + +#include +#include + +#define ALLOC(v, vl, sz) do { \ + size_t _sz = (sz); \ + mpw *_vv = xmalloc(MPWS(_sz)); \ + mpw *_vvl = _vv + _sz; \ + (v) = _vv; \ + (vl) = _vvl; \ +} while (0) + +#define LOAD(v, vl, d) do { \ + const dstr *_d = (d); \ + mpw *_v, *_vl; \ + ALLOC(_v, _vl, MPW_RQ(_d->len)); \ + mpx_loadb(_v, _vl, _d->buf, _d->len); \ + (v) = _v; \ + (vl) = _vl; \ +} while (0) + +#define MAX(x, y) ((x) > (y) ? (x) : (y)) + +static void dumpmp(const char *msg, const mpw *v, const mpw *vl) +{ + fputs(msg, stderr); + MPX_SHRINK(v, vl); + while (v < vl) + fprintf(stderr, " %08lx", (unsigned long)*--vl); + fputc('\n', stderr); +} + +static int mul(dstr *v) +{ + mpw *a, *al; + mpw *b, *bl; + mpw *c, *cl; + mpw *d, *dl; + mpw *s, *sl; + size_t m; + int ok = 1; + + LOAD(a, al, &v[0]); + LOAD(b, bl, &v[1]); + LOAD(c, cl, &v[2]); + m = MAX(al - a, bl - b) + 1; + ALLOC(d, dl, 2 * m); + ALLOC(s, sl, 2 * m); + + gfx_kmul(d, dl, a, al, b, bl, s, sl); + if (!mpx_ueq(d, dl, c, cl)) { + fprintf(stderr, "\n*** mul failed\n"); + dumpmp(" a", a, al); + dumpmp(" b", b, bl); + dumpmp("expected", c, cl); + dumpmp(" result", d, dl); + ok = 0; + } + + free(a); free(b); free(c); free(d); free(s); + return (ok); +} + +static test_chunk defs[] = { + { "mul", mul, { &type_hex, &type_hex, &type_hex, 0 } }, + { 0, 0, { 0 } } +}; + +int main(int argc, char *argv[]) +{ + test_run(argc, argv, defs, SRCDIR"/tests/gfx"); + return (0); +} + +#endif + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/gfx-sqr-mktab.c b/gfx-sqr-mktab.c new file mode 100644 index 0000000..70801c1 --- /dev/null +++ b/gfx-sqr-mktab.c @@ -0,0 +1,99 @@ +/* -*-c-*- + * + * $Id: gfx-sqr-mktab.c,v 1.1 2000/10/08 15:49:37 mdw Exp $ + * + * Build table for squaring of binary polynomials + * + * (c) 2000 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: gfx-sqr-mktab.c,v $ + * Revision 1.1 2000/10/08 15:49:37 mdw + * First glimmerings of binary polynomial arithmetic. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include +#include + +#include + +/*----- Main code ---------------------------------------------------------*/ + +static void mktab(uint16 *t) +{ + unsigned i, j, x; + + for (i = 0; i < 256; i++) { + x = 0; + for (j = 0; j < 8; j++) { + if (i & (1 << j)) + x |= 1 << (2 * j); + } + t[i] = x; + } +} + +int main(void) +{ + uint16 t[256]; + unsigned i; + + mktab(t); +fputs("\ +/* -*-c-*-\n\ + *\n\ + * Bit spacing table for binary polynomial squaring\n\ + */\n\ +\n\ +#ifndef GFX_SQR_TAB_H\n\ +#define GFX_SQR_TAB_H\n\ +\n\ +#define GFX_SQRTAB { \\\n\ + ", stdout); + + for (i = 0; i < 256; i++) { + printf("0x%04x", t[i]); + if (i == 255) + puts(" \\\n}\n"); + else if (i % 8 == 7) + fputs(", \\\n ", stdout); + else + fputs(", ", stdout); + } + + fputs("#endif\n", stdout); + + if (fclose(stdout)) { + fprintf(stderr, "error writing data\n"); + exit(EXIT_FAILURE); + } + + return (0); +} + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/gfx-sqr.c b/gfx-sqr.c new file mode 100644 index 0000000..b253a32 --- /dev/null +++ b/gfx-sqr.c @@ -0,0 +1,222 @@ +/* -*-c-*- + * + * $Id: gfx-sqr.c,v 1.1 2000/10/08 15:49:37 mdw Exp $ + * + * Sqaring binary polynomials + * + * (c) 2000 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: gfx-sqr.c,v $ + * Revision 1.1 2000/10/08 15:49:37 mdw + * First glimmerings of binary polynomial arithmetic. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include "mpx.h" +/* #include "gfx.h" */ +#include "gfx-sqr-tab.h" + +/*----- Static variables --------------------------------------------------*/ + +static uint16 tab[256] = GFX_SQRTAB; + +/*----- Main code ---------------------------------------------------------*/ + +/* --- @gfx_sqr@ --- * + * + * Arguments: @mpw *dv, *dvl@ = destination vector base and limit + * @const mpw *av, *avl@ = argument vector base and limit + * + * Returns: --- + * + * Use: Performs squaring of binary polynomials. + */ + +void gfx_sqr(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl) +{ + mpd a = 0, aa = 0; + unsigned b = 0, bb = 0; + + /* --- Simple stuff --- */ + + if (dv >= dvl) + return; + MPX_SHRINK(av, avl); + + /* --- The main algorithm --- * + * + * Our method depends on the fact that, in a field of characteristic 2, we + * have that %$(a + b)^2 = a^2 + b^2$%. Thus, to square a polynomial, it's + * sufficient just to put a zero bit between each of the bits of the + * original argument. We use a precomputed table for this, and work on + * entire octets at a time. Life is more complicated because we've got to + * be careful of bizarre architectures which don't have words with a + * multiple of 8 bits in them. + */ + + for (;;) { + + /* --- Input buffering --- */ + + if (b < 8) { + if (av >= avl) + break; + a |= *av++ << b; + b += MPW_BITS; + } + + /* --- Do the work in the middle --- */ + + aa |= (mpd)(tab[U8(a)]) << bb; + bb += 16; + a >>= 8; + b -= 8; + + /* --- Output buffering --- */ + + if (bb > MPW_BITS) { + *dv++ = MPW(aa); + if (dv >= dvl) + return; + aa >>= MPW_BITS; + bb -= MPW_BITS; + } + } + + /* --- Flush the input buffer --- */ + + if (b) for (;;) { + aa |= (mpd)(tab[U8(a)]) << bb; + bb += 16; + if (bb > MPW_BITS) { + *dv++ = MPW(aa); + if (dv >= dvl) + return; + aa >>= MPW_BITS; + bb -= MPW_BITS; + } + a >>= 8; + if (b <= 8) + break; + else + b -= 8; + } + + /* --- Flush the output buffer --- */ + + if (bb) for (;;) { + *dv++ = MPW(aa); + if (dv >= dvl) + return; + aa >>= MPW_BITS; + if (bb <= MPW_BITS) + break; + else + bb -= MPW_BITS; + } + + /* --- Zero the rest of everything --- */ + + MPX_ZERO(dv, dvl); +} + +/*----- Test rig ----------------------------------------------------------*/ + +#ifdef TEST_RIG + +#include +#include +#include +#include + +#define ALLOC(v, vl, sz) do { \ + size_t _sz = (sz); \ + mpw *_vv = xmalloc(MPWS(_sz)); \ + mpw *_vvl = _vv + _sz; \ + (v) = _vv; \ + (vl) = _vvl; \ +} while (0) + +#define LOAD(v, vl, d) do { \ + const dstr *_d = (d); \ + mpw *_v, *_vl; \ + ALLOC(_v, _vl, MPW_RQ(_d->len)); \ + mpx_loadb(_v, _vl, _d->buf, _d->len); \ + (v) = _v; \ + (vl) = _vl; \ +} while (0) + +#define MAX(x, y) ((x) > (y) ? (x) : (y)) + +static void dumpmp(const char *msg, const mpw *v, const mpw *vl) +{ + fputs(msg, stderr); + MPX_SHRINK(v, vl); + while (v < vl) + fprintf(stderr, " %08lx", (unsigned long)*--vl); + fputc('\n', stderr); +} + +static int vsqr(dstr *v) +{ + mpw *a, *al; + mpw *b, *bl; + mpw *d, *dl; + int ok = 1; + + LOAD(a, al, &v[0]); + LOAD(b, bl, &v[1]); + ALLOC(d, dl, 2 * (al - a)); + + gfx_sqr(d, dl, a, al); + if (!mpx_ueq(d, dl, b, bl)) { + fprintf(stderr, "\n*** vsqr failed\n"); + dumpmp(" a", a, al); + dumpmp("expected", b, bl); + dumpmp(" result", d, dl); + ok = 0; + } + + free(a); free(b); free(d); + return (ok); +} + +static test_chunk defs[] = { + { "sqr", vsqr, { &type_hex, &type_hex, 0 } }, + { 0, 0, { 0 } } +}; + +int main(int argc, char *argv[]) +{ + test_run(argc, argv, defs, SRCDIR"/tests/gfx"); + return (0); +} + +#endif + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/gfx.c b/gfx.c new file mode 100644 index 0000000..97320f9 --- /dev/null +++ b/gfx.c @@ -0,0 +1,397 @@ +/* -*-c-*- + * + * $Id: gfx.c,v 1.1 2000/10/08 15:49:37 mdw Exp $ + * + * Low-level arithmetic on binary polynomials + * + * (c) 2000 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: gfx.c,v $ + * Revision 1.1 2000/10/08 15:49:37 mdw + * First glimmerings of binary polynomial arithmetic. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include + +#include "mpx.h" +#include "mpscan.h" + +/*----- Main code ---------------------------------------------------------*/ + +/* --- @gfx_add@ --- * + * + * Arguments: @mpw *dv, *dvl@ = destination vector base and limit + * @const mpw *av, *avl@ = first addend vector base and limit + * @const mpw *bv, *bvl@ = second addend vector base and limit + * + * Returns: --- + * + * Use: Adds two %$\gf{2}$% polynomials. This is the same as + * subtraction. + */ + +void gfx_add(mpw *dv, mpw *dvl, + const mpw *av, const mpw *avl, + const mpw *bv, const mpw *bvl) +{ + MPX_SHRINK(av, avl); + MPX_SHRINK(bv, bvl); + + while (av < avl || bv < bvl) { + mpw a, b; + if (dv >= dvl) + return; + a = (av < avl) ? *av++ : 0; + b = (bv < bvl) ? *bv++ : 0; + *dv++ = a ^ b; + } + if (dv < dvl) + MPX_ZERO(dv, dvl); +} + +/* --- @gfx_acc@ --- * + * + * Arguments: @mpw *dv, *dvl@ = destination vector base and limit + * @const mpw *av, *avl@ = addend vector base and limit + * + * Returns: --- + * + * Use: Adds the addend into the destination. This is considerably + * faster than the three-address add call. + */ + +void gfx_acc(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl) +{ + size_t dlen, alen; + + MPX_SHRINK(av, avl); + dlen = dvl - dv; + alen = avl - av; + if (dlen < alen) + avl = av + dlen; + while (av < avl) + *dv++ ^= *av++; +} + +/* --- @gfx_accshift@ --- * + * + * Arguments: @mpw *dv, *dvl@ = destination vector base and limit + * @const mpw *av, *avl@ = addend vector base and limit + * @size_t n@ = number of bits to shift + * + * Returns: --- + * + * Use: Shifts the argument left by %$n$% places and adds it to the + * destination. This is a primitive used by multiplication and + * division. + */ + +void gfx_accshift(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, size_t n) +{ + size_t c = n / MPW_BITS; + mpw x = 0, y; + size_t dlen, alen; + + /* --- Sort out the shift amounts --- */ + + if (dvl - dv < c) + return; + dv += c; + n %= MPW_BITS; + if (!n) { + gfx_acc(dv, dvl, av, avl); + return; + } + c = MPW_BITS - n; + + /* --- Sort out vector lengths --- */ + + MPX_SHRINK(av, avl); + dlen = dvl - dv; + alen = avl - av; + if (dlen < alen) + avl = av + dlen; + + /* --- Now do the hard work --- */ + + while (av < avl) { + y = *av++; + *dv++ ^= MPW((y << n) | (x >> c)); + x = y; + } + if (dv < dvl) + *dv++ ^= x >> c; +} + +/* --- @gfx_mul@ --- * + * + * Arguments: @mpw *dv, *dvl@ = destination vector base and limit + * @const mpw *av, *avl@ = first argument vector base and limit + * @const mpw *bv, *bvl@ = second argument vector base and limit + * + * Returns: --- + * + * Use: Does multiplication of polynomials over %$\gf{2}$%. + */ + +void gfx_mul(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, + const mpw *bv, const mpw *bvl) +{ + mpscan sc; + const mpw *v; + mpw *vv; + mpw z; + mpd x, y; + + MPX_SHRINK(av, avl); + MPX_SHRINK(bv, bvl); + MPSCAN_INITX(&sc, av, avl); + MPX_ZERO(dv, dvl); + + while (bv < bvl && dv < dvl) { + x = 0; + for (v = av, vv = dv++; v < avl && vv < dvl; v++) { + z = *bv; y = *v; + while (z) { + if (z & 1u) x ^= y; + z >>= 1; y <<= 1; + } + *vv++ ^= MPW(x); + x >>= MPW_BITS; + } + if (vv < dvl) + *vv++ = MPW(x); + bv++; + } +} + +/* --- @gfx_div@ --- * + * + * Arguments: @mpw *qv, *qvl@ = quotient vector base and limit + * @mpw *rv, *rvl@ = dividend/remainder vector base and limit + * @const mpw *dv, *dvl@ = divisor vector base and limit + * + * Returns: --- + * + * Use: Performs division on polynomials over %$\gf{2}$%. + */ + +void gfx_div(mpw *qv, mpw *qvl, mpw *rv, mpw *rvl, + const mpw *dv, const mpw *dvl) +{ + size_t dlen, rlen, qlen; + size_t dbits; + mpw *rvv, *rvd; + unsigned rvm, n, qi; + mpw q; + + MPX_SHRINK(rv, rvl); + MPX_SHRINK(dv, dvl); + assert(((void)"division by zero in gfx_div", dv < dvl)); + MPX_BITS(dbits, dv, dvl); + dlen = dvl - dv; + rlen = rvl - rv; + qlen = qvl - qv; + + MPX_ZERO(qv, qvl); + if (dlen > rlen) + return; + rvd = rvl - dlen; + rvv = rvl - 1; + rvm = 1 << (MPW_BITS - 1); + n = MPW_BITS - (dbits % MPW_BITS); + if (n == MPW_BITS) + n = 0; + q = 0; + qi = rvd - rv; + + for (;;) { + q <<= 1; + if (*rvv & rvm) { + q |= 1; + gfx_accshift(rvd, rvl, dv, dvl, n); + } + rvm >>= 1; + if (!rvm) { + rvm = 1 << (MPW_BITS - 1); + rvv--; + } + if (n) + n--; + else { + if (qi < qlen) + qv[qi] = q; + q = 0; + qi--; + if (rvd == rv) + break; + n = MPW_BITS - 1; + rvd--; + } + } +} + +/*----- Test rig ----------------------------------------------------------*/ + +#ifdef TEST_RIG + +#include +#include +#include +#include + +#define ALLOC(v, vl, sz) do { \ + size_t _sz = (sz); \ + mpw *_vv = xmalloc(MPWS(_sz)); \ + mpw *_vvl = _vv + _sz; \ + (v) = _vv; \ + (vl) = _vvl; \ +} while (0) + +#define LOAD(v, vl, d) do { \ + const dstr *_d = (d); \ + mpw *_v, *_vl; \ + ALLOC(_v, _vl, MPW_RQ(_d->len)); \ + mpx_loadb(_v, _vl, _d->buf, _d->len); \ + (v) = _v; \ + (vl) = _vl; \ +} while (0) + +#define MAX(x, y) ((x) > (y) ? (x) : (y)) + +static void dumpmp(const char *msg, const mpw *v, const mpw *vl) +{ + fputs(msg, stderr); + MPX_SHRINK(v, vl); + while (v < vl) + fprintf(stderr, " %08lx", (unsigned long)*--vl); + fputc('\n', stderr); +} + +static int vadd(dstr *v) +{ + mpw *a, *al; + mpw *b, *bl; + mpw *c, *cl; + mpw *d, *dl; + int ok = 1; + + LOAD(a, al, &v[0]); + LOAD(b, bl, &v[1]); + LOAD(c, cl, &v[2]); + ALLOC(d, dl, MAX(al - a, bl - b) + 1); + + gfx_add(d, dl, a, al, b, bl); + if (!mpx_ueq(d, dl, c, cl)) { + fprintf(stderr, "\n*** vadd failed\n"); + dumpmp(" a", a, al); + dumpmp(" b", b, bl); + dumpmp("expected", c, cl); + dumpmp(" result", d, dl); + ok = 0; + } + + free(a); free(b); free(c); free(d); + return (ok); +} + +static int vmul(dstr *v) +{ + mpw *a, *al; + mpw *b, *bl; + mpw *c, *cl; + mpw *d, *dl; + int ok = 1; + + LOAD(a, al, &v[0]); + LOAD(b, bl, &v[1]); + LOAD(c, cl, &v[2]); + ALLOC(d, dl, (al - a) + (bl - b)); + + gfx_mul(d, dl, a, al, b, bl); + if (!mpx_ueq(d, dl, c, cl)) { + fprintf(stderr, "\n*** vmul failed\n"); + dumpmp(" a", a, al); + dumpmp(" b", b, bl); + dumpmp("expected", c, cl); + dumpmp(" result", d, dl); + ok = 0; + } + + free(a); free(b); free(c); free(d); + return (ok); +} + +static int vdiv(dstr *v) +{ + mpw *a, *al; + mpw *b, *bl; + mpw *q, *ql; + mpw *r, *rl; + mpw *qq, *qql; + int ok = 1; + + ALLOC(a, al, MPW_RQ(v[0].len) + 2); mpx_loadb(a, al, v[0].buf, v[0].len); + LOAD(b, bl, &v[1]); + LOAD(q, ql, &v[2]); + LOAD(r, rl, &v[3]); + ALLOC(qq, qql, al - a); + + gfx_div(qq, qql, a, al, b, bl); + if (!mpx_ueq(qq, qql, q, ql) || + !mpx_ueq(a, al, r, rl)) { + fprintf(stderr, "\n*** vdiv failed\n"); + dumpmp(" divisor", b, bl); + dumpmp("expect r", r, rl); + dumpmp("result r", a, al); + dumpmp("expect q", q, ql); + dumpmp("result q", qq, qql); + ok = 0; + } + + free(a); free(b); free(r); free(q); free(qq); + return (ok); +} + +static test_chunk defs[] = { + { "add", vadd, { &type_hex, &type_hex, &type_hex, 0 } }, + { "mul", vmul, { &type_hex, &type_hex, &type_hex, 0 } }, + { "div", vdiv, { &type_hex, &type_hex, &type_hex, &type_hex, 0 } }, + { 0, 0, { 0 } } +}; + +int main(int argc, char *argv[]) +{ + test_run(argc, argv, defs, SRCDIR"/tests/gfx"); + return (0); +} + +#endif + +/*----- That's all, folks -------------------------------------------------*/ diff --git a/gfx.h b/gfx.h new file mode 100644 index 0000000..778baee --- /dev/null +++ b/gfx.h @@ -0,0 +1,170 @@ +/* -*-c-*- + * + * $Id: gfx.h,v 1.1 2000/10/08 15:49:37 mdw Exp $ + * + * Low-level arithmetic on binary polynomials + * + * (c) 2000 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: gfx.h,v $ + * Revision 1.1 2000/10/08 15:49:37 mdw + * First glimmerings of binary polynomial arithmetic. + * + */ + +#ifndef CATACOMB_GFX_H +#define CATACOMB_GFX_H + +#ifdef __cplusplus + extern "C" { +#endif + +/*----- Header files ------------------------------------------------------*/ + +#ifndef CATACOMB_MPX_H +# include "mpx.h" +#endif + +/*----- Functions provided ------------------------------------------------*/ + +/* --- @gfx_add@ --- * + * + * Arguments: @mpw *dv, *dvl@ = destination vector base and limit + * @const mpw *av, *avl@ = first addend vector base and limit + * @const mpw *bv, *bvl@ = second addend vector base and limit + * + * Returns: --- + * + * Use: Adds two %$\gf{2}$% polynomials. This is the same as + * subtraction. + */ + +extern void gfx_add(mpw */*dv*/, mpw */*dvl*/, + const mpw */*av*/, const mpw */*avl*/, + const mpw */*bv*/, const mpw */*bvl*/); + +/* --- @gfx_acc@ --- * + * + * Arguments: @mpw *dv, *dvl@ = destination vector base and limit + * @const mpw *av, *avl@ = addend vector base and limit + * + * Returns: --- + * + * Use: Adds the addend into the destination. This is considerably + * faster than the three-address add call. + */ + +extern void gfx_acc(mpw */*dv*/, mpw */*dvl*/, + const mpw */*av*/, const mpw */*avl*/); + +/* --- @gfx_accshift@ --- * + * + * Arguments: @mpw *dv, *dvl@ = destination vector base and limit + * @const mpw *av, *avl@ = addend vector base and limit + * @size_t n@ = number of bits to shift + * + * Returns: --- + * + * Use: Shifts the argument left by %$n$% places and adds it to the + * destination. This is a primitive used by multiplication and + * division. + */ + +extern void gfx_accshift(mpw */*dv*/, mpw */*dvl*/, + const mpw */*av*/, const mpw */*avl*/, + size_t /*n*/); + +/* --- @gfx_mul@ --- * + * + * Arguments: @mpw *dv, *dvl@ = destination vector base and limit + * @const mpw *av, *avl@ = first argument vector base and limit + * @const mpw *bv, *bvl@ = second argument vector base and limit + * + * Returns: --- + * + * Use: Does multiplication of polynomials over %$\gf{2}$%. + */ + +extern void gfx_mul(mpw */*dv*/, mpw */*dvl*/, + const mpw */*av*/, const mpw */*avl*/, + const mpw */*bv*/, const mpw */*bvl*/); + +/* --- @gfx_div@ --- * + * + * Arguments: @mpw *qv, *qvl@ = quotient vector base and limit + * @mpw *rv, *rvl@ = dividend/remainder vector base and limit + * @const mpw *dv, *dvl@ = divisor vector base and limit + * + * Returns: --- + * + * Use: Performs division on polynomials over %$\gf{2}$%. + */ + +extern void gfx_div(mpw */*qv*/, mpw */*qvl*/, mpw */*rv*/, mpw */*rvl*/, + const mpw */*dv*/, const mpw */*dvl*/); + +/*----- Karatsuba multiplication algorithms -------------------------------*/ + +/* --- @GFK_THRESH@ --- * + * + * This is the limiting length for using Karatsuba algorithms. It's best to + * use the simpler classical multiplication method on numbers smaller than + * this. + */ + +#define GFK_THRESH 2 + +/* --- @gfx_kmul@ --- * + * + * Arguments: @mpw *dv, *dvl@ = pointer to destination buffer + * @const mpw *av, *avl@ = pointer to first argument + * @const mpw *bv, *bvl@ = pointer to second argument + * @mpw *sv, *svl@ = pointer to scratch workspace + * + * Returns: --- + * + * Use: Multiplies two binary polynomials using Karatsuba's + * algorithm. This is rather faster than traditional long + * multiplication (e.g., @gfx_umul@) on polynomials with large + * degree, although more expensive on small ones. + * + * The destination must be twice as large as the larger + * argument. The scratch space must be twice as large as the + * larger argument. + */ + +extern void gfx_kmul(mpw */*dv*/, mpw */*dvl*/, + const mpw */*av*/, const mpw */*avl*/, + const mpw */*bv*/, const mpw */*bvl*/, + mpw */*sv*/, mpw */*svl*/); + +/*----- That's all, folks -------------------------------------------------*/ + +#ifdef __cplusplus + } +#endif + +#endif -- 2.11.0