chiark - git - mdw - catacomb/blob - gfx.c

   1 /* -*-c-*-
   2  *
   3  * $Id$
   4  *
   5  * Low-level arithmetic 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 /*----- Header files ------------------------------------------------------*/
  31
  32 #include <assert.h>
  33
  34 #include "mpx.h"
  35 #include "mpscan.h"
  36
  37 /*----- Main code ---------------------------------------------------------*/
  38
  39 /* --- @gfx_add@ --- *
  40  *
  41  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
  42  *              @const mpw *av, *avl@ = first addend vector base and limit
  43  *              @const mpw *bv, *bvl@ = second addend vector base and limit
  44  *
  45  * Returns:     ---
  46  *
  47  * Use:         Adds two %$\gf{2}$% polynomials.  This is the same as
  48  *              subtraction.
  49  */
  50
  51 void gfx_add(mpw *dv, mpw *dvl,
  52              const mpw *av, const mpw *avl,
  53              const mpw *bv, const mpw *bvl)
  54 {
  55   MPX_SHRINK(av, avl);
  56   MPX_SHRINK(bv, bvl);
  57
  58   while (av < avl || bv < bvl) {
  59     mpw a, b;
  60     if (dv >= dvl)
  61       return;
  62     a = (av < avl) ? *av++ : 0;
  63     b = (bv < bvl) ? *bv++ : 0;
  64     *dv++ = a ^ b;
  65   }
  66   if (dv < dvl)
  67     MPX_ZERO(dv, dvl);
  68 }
  69
  70 /* --- @gfx_acc@ --- *
  71  *
  72  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
  73  *              @const mpw *av, *avl@ = addend vector base and limit
  74  *
  75  * Returns:     ---
  76  *
  77  * Use:         Adds the addend into the destination.  This is considerably
  78  *              faster than the three-address add call.
  79  */
  80
  81 void gfx_acc(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl)
  82 {
  83   size_t dlen, alen;
  84
  85   MPX_SHRINK(av, avl);
  86   dlen = dvl - dv;
  87   alen = avl - av;
  88   if (dlen < alen)
  89     avl = av + dlen;
  90   while (av < avl)
  91     *dv++ ^= *av++;
  92 }
  93
  94 /* --- @gfx_accshift@ --- *
  95  *
  96  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
  97  *              @const mpw *av, *avl@ = addend vector base and limit
  98  *              @size_t n@ = number of bits to shift
  99  *
 100  * Returns:     ---
 101  *
 102  * Use:         Shifts the argument left by %$n$% places and adds it to the
 103  *              destination.  This is a primitive used by multiplication and
 104  *              division.
 105  */
 106
 107 void gfx_accshift(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, size_t n)
 108 {
 109   size_t c = n / MPW_BITS;
 110   mpw x = 0, y;
 111   size_t dlen, alen;
 112
 113   /* --- Sort out the shift amounts --- */
 114
 115   if (dvl - dv < c)
 116     return;
 117   dv += c;
 118   n %= MPW_BITS;
 119   if (!n) {
 120     gfx_acc(dv, dvl, av, avl);
 121     return;
 122   }
 123   c = MPW_BITS - n;
 124
 125   /* --- Sort out vector lengths --- */
 126
 127   MPX_SHRINK(av, avl);
 128   dlen = dvl - dv;
 129   alen = avl - av;
 130   if (dlen < alen)
 131     avl = av + dlen;
 132
 133   /* --- Now do the hard work --- */
 134
 135   while (av < avl) {
 136     y = *av++;
 137     *dv++ ^= MPW((y << n) | (x >> c));
 138     x = y;
 139   }
 140   if (dv < dvl)
 141     *dv++ ^= x >> c;
 142 }
 143
 144 /* --- @gfx_mul@ --- *
 145  *
 146  * Arguments:   @mpw *dv, *dvl@ = destination vector base and limit
 147  *              @const mpw *av, *avl@ = first argument vector base and limit
 148  *              @const mpw *bv, *bvl@ = second argument vector base and limit
 149  *
 150  * Returns:     ---
 151  *
 152  * Use:         Does multiplication of polynomials over %$\gf{2}$%.
 153  */
 154
 155 void gfx_mul(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
 156              const mpw *bv, const mpw *bvl)
 157 {
 158   mpscan sc;
 159   const mpw *v;
 160   mpw *vv;
 161   mpw z;
 162   mpd x, y;
 163
 164   MPX_SHRINK(av, avl);
 165   MPX_SHRINK(bv, bvl);
 166   MPSCAN_INITX(&sc, av, avl);
 167   MPX_ZERO(dv, dvl);
 168
 169   while (bv < bvl && dv < dvl) {
 170     x = 0;
 171     for (v = av, vv = dv++; v < avl && vv < dvl; v++) {
 172       z = *bv; y = *v;
 173       while (z) {
 174         if (z & 1u) x ^= y;
 175         z >>= 1; y <<= 1;
 176       }
 177       *vv++ ^= MPW(x);
 178       x >>= MPW_BITS;
 179     }
 180     if (vv < dvl)
 181       *vv++ = MPW(x);
 182     bv++;
 183   }
 184 }
 185
 186 /* --- @gfx_div@ --- *
 187  *
 188  * Arguments:   @mpw *qv, *qvl@ = quotient vector base and limit
 189  *              @mpw *rv, *rvl@ = dividend/remainder vector base and limit
 190  *              @const mpw *dv, *dvl@ = divisor vector base and limit
 191  *
 192  * Returns:     ---
 193  *
 194  * Use:         Performs division on polynomials over %$\gf{2}$%.
 195  */
 196
 197 void gfx_div(mpw *qv, mpw *qvl, mpw *rv, mpw *rvl,
 198              const mpw *dv, const mpw *dvl)
 199 {
 200   size_t dlen, rlen, qlen;
 201   size_t dbits;
 202   mpw *rvv, *rvd;
 203   unsigned rvm, n, qi;
 204   mpw q;
 205
 206   MPX_SHRINK(rv, rvl);
 207   MPX_SHRINK(dv, dvl);
 208   assert(((void)"division by zero in gfx_div", dv < dvl));
 209   MPX_BITS(dbits, dv, dvl);
 210   dlen = dvl - dv;
 211   rlen = rvl - rv;
 212   qlen = qvl - qv;
 213
 214   MPX_ZERO(qv, qvl);
 215   if (dlen > rlen)
 216     return;
 217   rvd = rvl - dlen;
 218   rvv = rvl - 1;
 219   rvm = 1 << (MPW_BITS - 1);
 220   n = MPW_BITS - (dbits % MPW_BITS);
 221   if (n == MPW_BITS)
 222     n = 0;
 223   q = 0;
 224   qi = rvd - rv;
 225
 226   for (;;) {
 227     q <<= 1;
 228     if (*rvv & rvm) {
 229       q |= 1;
 230       gfx_accshift(rvd, rvl, dv, dvl, n);
 231     }
 232     rvm >>= 1;
 233     if (!rvm) {
 234       rvm = 1 << (MPW_BITS - 1);
 235       rvv--;
 236     }
 237     if (n)
 238       n--;
 239     else {
 240       if (qi < qlen)
 241         qv[qi] = q;
 242       q = 0;
 243       qi--;
 244       if (rvd == rv)
 245         break;
 246       n = MPW_BITS - 1;
 247       rvd--;
 248     }
 249   }
 250 }
 251
 252 /*----- Test rig ----------------------------------------------------------*/
 253
 254 #ifdef TEST_RIG
 255
 256 #include <mLib/alloc.h>
 257 #include <mLib/dstr.h>
 258 #include <mLib/quis.h>
 259 #include <mLib/testrig.h>
 260
 261 #define ALLOC(v, vl, sz) do {                                           \
 262   size_t _sz = (sz);                                                    \
 263   mpw *_vv = xmalloc(MPWS(_sz));                                        \
 264   mpw *_vvl = _vv + _sz;                                                \
 265   (v) = _vv;                                                            \
 266   (vl) = _vvl;                                                          \
 267 } while (0)
 268
 269 #define LOAD(v, vl, d) do {                                             \
 270   const dstr *_d = (d);                                                 \
 271   mpw *_v, *_vl;                                                        \
 272   ALLOC(_v, _vl, MPW_RQ(_d->len));                                      \
 273   mpx_loadb(_v, _vl, _d->buf, _d->len);                                 \
 274   (v) = _v;                                                             \
 275   (vl) = _vl;                                                           \
 276 } while (0)
 277
 278 #define MAX(x, y) ((x) > (y) ? (x) : (y))
 279
 280 static void dumpmp(const char *msg, const mpw *v, const mpw *vl)
 281 {
 282   fputs(msg, stderr);
 283   MPX_SHRINK(v, vl);
 284   while (v < vl)
 285     fprintf(stderr, " %08lx", (unsigned long)*--vl);
 286   fputc('\n', stderr);
 287 }
 288
 289 static int vadd(dstr *v)
 290 {
 291   mpw *a, *al;
 292   mpw *b, *bl;
 293   mpw *c, *cl;
 294   mpw *d, *dl;
 295   int ok = 1;
 296
 297   LOAD(a, al, &v[0]);
 298   LOAD(b, bl, &v[1]);
 299   LOAD(c, cl, &v[2]);
 300   ALLOC(d, dl, MAX(al - a, bl - b) + 1);
 301
 302   gfx_add(d, dl, a, al, b, bl);
 303   if (!mpx_ueq(d, dl, c, cl)) {
 304     fprintf(stderr, "\n*** vadd failed\n");
 305     dumpmp("       a", a, al);
 306     dumpmp("       b", b, bl);
 307     dumpmp("expected", c, cl);
 308     dumpmp("  result", d, dl);
 309     ok = 0;
 310   }
 311
 312   xfree(a); xfree(b); xfree(c); xfree(d);
 313   return (ok);
 314 }
 315
 316 static int vmul(dstr *v)
 317 {
 318   mpw *a, *al;
 319   mpw *b, *bl;
 320   mpw *c, *cl;
 321   mpw *d, *dl;
 322   int ok = 1;
 323
 324   LOAD(a, al, &v[0]);
 325   LOAD(b, bl, &v[1]);
 326   LOAD(c, cl, &v[2]);
 327   ALLOC(d, dl, (al - a) + (bl - b));
 328
 329   gfx_mul(d, dl, a, al, b, bl);
 330   if (!mpx_ueq(d, dl, c, cl)) {
 331     fprintf(stderr, "\n*** vmul failed\n");
 332     dumpmp("       a", a, al);
 333     dumpmp("       b", b, bl);
 334     dumpmp("expected", c, cl);
 335     dumpmp("  result", d, dl);
 336     ok = 0;
 337   }
 338
 339   xfree(a); xfree(b); xfree(c); xfree(d);
 340   return (ok);
 341 }
 342
 343 static int vdiv(dstr *v)
 344 {
 345   mpw *a, *al;
 346   mpw *b, *bl;
 347   mpw *q, *ql;
 348   mpw *r, *rl;
 349   mpw *qq, *qql;
 350   int ok = 1;
 351
 352   ALLOC(a, al, MPW_RQ(v[0].len) + 2); mpx_loadb(a, al, v[0].buf, v[0].len);
 353   LOAD(b, bl, &v[1]);
 354   LOAD(q, ql, &v[2]);
 355   LOAD(r, rl, &v[3]);
 356   ALLOC(qq, qql, al - a);
 357
 358   gfx_div(qq, qql, a, al, b, bl);
 359   if (!mpx_ueq(qq, qql, q, ql) ||
 360       !mpx_ueq(a, al, r, rl)) {
 361     fprintf(stderr, "\n*** vdiv failed\n");
 362     dumpmp(" divisor", b, bl);
 363     dumpmp("expect r", r, rl);
 364     dumpmp("result r", a, al);
 365     dumpmp("expect q", q, ql);
 366     dumpmp("result q", qq, qql);
 367     ok = 0;
 368   }
 369
 370   xfree(a); xfree(b); xfree(r); xfree(q); xfree(qq);
 371   return (ok);
 372 }
 373
 374 static test_chunk defs[] = {
 375   { "add", vadd, { &type_hex, &type_hex, &type_hex, 0 } },
 376   { "mul", vmul, { &type_hex, &type_hex, &type_hex, 0 } },
 377   { "div", vdiv, { &type_hex, &type_hex, &type_hex, &type_hex, 0 } },
 378   { 0, 0, { 0 } }
 379 };
 380
 381 int main(int argc, char *argv[])
 382 {
 383   test_run(argc, argv, defs, SRCDIR"/tests/gfx");
 384   return (0);
 385 }
 386
 387 #endif
 388
 389 /*----- That's all, folks -------------------------------------------------*/