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

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