chiark / gitweb /
rand/rand.c (rand_gate): Evolve r->ibits in a more sensible manner.
[catacomb] / 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 -------------------------------------------------*/