chiark / gitweb /
math/ptab.in: Include the correct Oakley 2048 group!
[catacomb] / math / gfx.c
CommitLineData
ae747c9b 1/* -*-c-*-
ae747c9b 2 *
3 * Low-level arithmetic on binary polynomials
4 *
5 * (c) 2000 Straylight/Edgeware
6 */
7
45c0fd36 8/*----- Licensing notice --------------------------------------------------*
ae747c9b 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.
45c0fd36 16 *
ae747c9b 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.
45c0fd36 21 *
ae747c9b 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
ae747c9b 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
49void 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
79void 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
105void 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
153void 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
45c0fd36
MW
187 * @mpw *rv, *rvl@ = dividend/remainder vector base and limit
188 * @const mpw *dv, *dvl@ = divisor vector base and limit
ae747c9b 189 *
45c0fd36 190 * Returns: ---
ae747c9b 191 *
45c0fd36 192 * Use: Performs division on polynomials over %$\gf{2}$%.
ae747c9b 193 */
194
195void 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))
45c0fd36 277
ae747c9b 278static 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
287static 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");
45c0fd36
MW
303 dumpmp(" a", a, al);
304 dumpmp(" b", b, bl);
ae747c9b 305 dumpmp("expected", c, cl);
306 dumpmp(" result", d, dl);
307 ok = 0;
308 }
309
12ed8a1f 310 xfree(a); xfree(b); xfree(c); xfree(d);
ae747c9b 311 return (ok);
312}
313
314static 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");
45c0fd36
MW
330 dumpmp(" a", a, al);
331 dumpmp(" b", b, bl);
ae747c9b 332 dumpmp("expected", c, cl);
333 dumpmp(" result", d, dl);
334 ok = 0;
335 }
336
12ed8a1f 337 xfree(a); xfree(b); xfree(c); xfree(d);
ae747c9b 338 return (ok);
339}
340
341static 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
12ed8a1f 368 xfree(a); xfree(b); xfree(r); xfree(q); xfree(qq);
ae747c9b 369 return (ok);
370}
371
372static 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
379int main(int argc, char *argv[])
380{
0f00dc4c 381 test_run(argc, argv, defs, SRCDIR"/t/gfx");
ae747c9b 382 return (0);
383}
384
385#endif
386
387/*----- That's all, folks -------------------------------------------------*/