git.distorted.org.uk Git - u/mdw/catacomb/blob - gfx-kmul.c

   1 /* -*-c-*-
   2  *
   3  * $Id: gfx-kmul.c,v 1.1 2000/10/08 15:49:37 mdw Exp $
   4  *
   5  * Karatsuba's multiplication algorithm on binary polynomials
   6  *
   7  * (c) 2000 Straylight/Edgeware
   8  */
   9
  10 /*----- Licensing notice --------------------------------------------------*
  11  *
  12  * This file is part of Catacomb.
  13  *
  14  * Catacomb is free software; you can redistribute it and/or modify
  15  * it under the terms of the GNU Library General Public License as
  16  * published by the Free Software Foundation; either version 2 of the
  17  * License, or (at your option) any later version.
  18  *
  19  * Catacomb is distributed in the hope that it will be useful,
  20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  22  * GNU Library General Public License for more details.
  23  *
  24  * You should have received a copy of the GNU Library General Public
  25  * License along with Catacomb; if not, write to the Free
  26  * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  27  * MA 02111-1307, USA.
  28  */
  29
  30 /*----- Revision history --------------------------------------------------*
  31  *
  32  * $Log: gfx-kmul.c,v $
  33  * Revision 1.1  2000/10/08 15:49:37  mdw
  34  * First glimmerings of binary polynomial arithmetic.
  35  *
  36  */
  37
  38 /*----- Header files ------------------------------------------------------*/
  39
  40 #include <assert.h>
  41 #include <stdio.h>
  42
  43 #include "gfx.h"
  44 #include "karatsuba.h"
  45
  46 /*----- Tweakables --------------------------------------------------------*/
  47
  48 #ifdef TEST_RIG
  49 #  undef GFK_THRESH
  50 #  define GFK_THRESH 1
  51 #endif
  52
  53 /*----- Main code ---------------------------------------------------------*/
  54
  55 /* --- @gfx_kmul@ --- *
  56  *
  57  * Arguments:   @mpw *dv, *dvl@ = pointer to destination buffer
  58  *              @const mpw *av, *avl@ = pointer to first argument
  59  *              @const mpw *bv, *bvl@ = pointer to second argument
  60  *              @mpw *sv, *svl@ = pointer to scratch workspace
  61  *
  62  * Returns:     ---
  63  *
  64  * Use:         Multiplies two binary polynomials using Karatsuba's
  65  *              algorithm.  This is rather faster than traditional long
  66  *              multiplication (e.g., @gfx_umul@) on polynomials with large
  67  *              degree, although more expensive on small ones.
  68  *
  69  *              The destination must be twice as large as the larger
  70  *              argument.  The scratch space must be twice as large as the
  71  *              larger argument.
  72  */
  73
  74 void gfx_kmul(mpw *dv, mpw *dvl,
  75               const mpw *av, const mpw *avl,
  76               const mpw *bv, const mpw *bvl,
  77               mpw *sv, mpw *svl)
  78 {
  79   const mpw *avm, *bvm;
  80   size_t m;
  81
  82   /* --- Dispose of easy cases to @mpx_umul@ --- *
  83    *
  84    * Karatsuba is only a win on large numbers, because of all the
  85    * recursiveness and bookkeeping.  The recursive calls make a quick check
  86    * to see whether to bottom out to @gfx_umul@ which should help quite a
  87    * lot, but sometimes the only way to know is to make sure...
  88    */
  89
  90   MPX_SHRINK(av, avl);
  91   MPX_SHRINK(bv, bvl);
  92
  93   if (avl - av <= GFK_THRESH || bvl - bv <= GFK_THRESH) {
  94     gfx_mul(dv, dvl, av, avl, bv, bvl);
  95     return;
  96   }
  97
  98   /* --- How the algorithm works --- *
  99    *
 100    * Let %$A = xb + y$% and %$B = ub + v$%.  Then, simply by expanding,
 101    * %$AB = x u b^2 + b(x v + y u) + y v$%.  That's not helped any, because
 102    * I've got four multiplications, each four times easier than the one I
 103    * started with.  However, note that I can rewrite the coefficient of %$b$%
 104    * as %$xv + yu = (x + y)(u + v) - xu - yv$%.  The terms %$xu$% and %$yv$%
 105    * I've already calculated, and that leaves only one more multiplication to
 106    * do.  So now I have three multiplications, each four times easier, and
 107    * that's a win.
 108    */
 109
 110   /* --- First things --- *
 111    *
 112    * Sort out where to break the factors in half.  I'll choose the midpoint
 113    * of the larger one, since this minimizes the amount of work I have to do
 114    * most effectively.
 115    */
 116
 117   if (avl - av > bvl - bv) {
 118     m = (avl - av + 1) >> 1;
 119     avm = av + m;
 120     if (bvl - bv > m)
 121       bvm = bv + m;
 122     else
 123       bvm = bvl;
 124   } else {
 125     m = (bvl - bv + 1) >> 1;
 126     bvm = bv + m;
 127     if (avl - av > m)
 128       avm = av + m;
 129     else
 130       avm = avl;
 131   }
 132
 133   assert(((void)"Destination too small for Karatsuba gf-multiply",
 134           dvl - dv >= 4 * m));
 135   assert(((void)"Not enough workspace for Karatsuba gf-multiply",
 136           svl - sv >= 4 * m));
 137
 138   /* --- Sort out the middle term --- */
 139
 140   {
 141     mpw *bsv = sv + m, *ssv = bsv + m;
 142     mpw *rdv = dv + m, *rdvl = rdv + 2 * m;
 143
 144     UXOR2(sv, bsv, av, avm, avm, avl);
 145     UXOR2(bsv, ssv, bv, bvm, bvm, bvl);
 146     if (m > GFK_THRESH)
 147       gfx_kmul(rdv, rdvl, sv, bsv, bsv, ssv, ssv, svl);
 148     else
 149       gfx_mul(rdv, rdvl, sv, bsv, bsv, ssv);
 150   }
 151
 152   /* --- Sort out the other two terms --- */
 153
 154   {
 155     mpw *svm = sv + m, *ssv = svm + m;
 156     mpw *tdv = dv + m;
 157     mpw *rdv = tdv + m;
 158
 159     if (avl == avm || bvl == bvm)
 160       MPX_ZERO(rdv + m, dvl);
 161     else {
 162       if (m > GFK_THRESH)
 163         gfx_kmul(sv, ssv, avm, avl, bvm, bvl, ssv, svl);
 164       else
 165         gfx_mul(sv, ssv, avm, avl, bvm, bvl);
 166       MPX_COPY(rdv + m, dvl, svm, ssv);
 167       UXOR(rdv, sv, svm);
 168       UXOR(tdv, sv, ssv);
 169     }
 170
 171     if (m > GFK_THRESH)
 172       gfx_kmul(sv, ssv, av, avm, bv, bvm, ssv, svl);
 173     else
 174       gfx_mul(sv, ssv, av, avm, bv, bvm);
 175     MPX_COPY(dv, tdv, sv, svm);
 176     UXOR(tdv, sv, ssv);
 177     UXOR(tdv, svm, ssv);
 178   }
 179 }
 180
 181 /*----- Test rig ----------------------------------------------------------*/
 182
 183 #ifdef TEST_RIG
 184
 185 #include <mLib/alloc.h>
 186 #include <mLib/testrig.h>
 187
 188 #define ALLOC(v, vl, sz) do {                                           \
 189   size_t _sz = (sz);                                                    \
 190   mpw *_vv = xmalloc(MPWS(_sz));                                        \
 191   mpw *_vvl = _vv + _sz;                                                \
 192   (v) = _vv;                                                            \
 193   (vl) = _vvl;                                                          \
 194 } while (0)
 195
 196 #define LOAD(v, vl, d) do {                                             \
 197   const dstr *_d = (d);                                                 \
 198   mpw *_v, *_vl;                                                        \
 199   ALLOC(_v, _vl, MPW_RQ(_d->len));                                      \
 200   mpx_loadb(_v, _vl, _d->buf, _d->len);                                 \
 201   (v) = _v;                                                             \
 202   (vl) = _vl;                                                           \
 203 } while (0)
 204
 205 #define MAX(x, y) ((x) > (y) ? (x) : (y))
 206
 207 static void dumpmp(const char *msg, const mpw *v, const mpw *vl)
 208 {
 209   fputs(msg, stderr);
 210   MPX_SHRINK(v, vl);
 211   while (v < vl)
 212     fprintf(stderr, " %08lx", (unsigned long)*--vl);
 213   fputc('\n', stderr);
 214 }
 215
 216 static int mul(dstr *v)
 217 {
 218   mpw *a, *al;
 219   mpw *b, *bl;
 220   mpw *c, *cl;
 221   mpw *d, *dl;
 222   mpw *s, *sl;
 223   size_t m;
 224   int ok = 1;
 225
 226   LOAD(a, al, &v[0]);
 227   LOAD(b, bl, &v[1]);
 228   LOAD(c, cl, &v[2]);
 229   m = MAX(al - a, bl - b) + 1;
 230   ALLOC(d, dl, 2 * m);
 231   ALLOC(s, sl, 2 * m);
 232
 233   gfx_kmul(d, dl, a, al, b, bl, s, sl);
 234   if (!mpx_ueq(d, dl, c, cl)) {
 235     fprintf(stderr, "\n*** mul failed\n");
 236     dumpmp("       a", a, al);
 237     dumpmp("       b", b, bl);
 238     dumpmp("expected", c, cl);
 239     dumpmp("  result", d, dl);
 240     ok = 0;
 241   }
 242
 243   free(a); free(b); free(c); free(d); free(s);
 244   return (ok);
 245 }
 246
 247 static test_chunk defs[] = {
 248   { "mul", mul, { &type_hex, &type_hex, &type_hex, 0 } },
 249   { 0, 0, { 0 } }
 250 };
 251
 252 int main(int argc, char *argv[])
 253 {
 254   test_run(argc, argv, defs, SRCDIR"/tests/gfx");
 255   return (0);
 256 }
 257
 258 #endif
 259
 260 /*----- That's all, folks -------------------------------------------------*/