--- /dev/null
+/* -*-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 <assert.h>
+#include <stdio.h>
+
+#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 <mLib/alloc.h>
+#include <mLib/testrig.h>
+
+#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 -------------------------------------------------*/
--- /dev/null
+/* -*-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 <stdio.h>
+#include <stdlib.h>
+
+#include <mLib/bits.h>
+
+/*----- 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 -------------------------------------------------*/
--- /dev/null
+/* -*-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 <mLib/alloc.h>
+#include <mLib/dstr.h>
+#include <mLib/quis.h>
+#include <mLib/testrig.h>
+
+#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 -------------------------------------------------*/
--- /dev/null
+/* -*-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 <assert.h>
+
+#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 <mLib/alloc.h>
+#include <mLib/dstr.h>
+#include <mLib/quis.h>
+#include <mLib/testrig.h>
+
+#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 -------------------------------------------------*/
--- /dev/null
+/* -*-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
projects / u / mdw / catacomb / commitdiff