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

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