chiark / gitweb /
cleanup: Big pile of whitespace fixes, all at once.
[catacomb] / gfreduce.c
1 /* -*-c-*-
2  *
3  * $Id$
4  *
5  * Efficient reduction modulo sparse binary polynomials
6  *
7  * (c) 2004 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 <mLib/alloc.h>
33 #include <mLib/darray.h>
34 #include <mLib/macros.h>
35
36 #include "gf.h"
37 #include "gfreduce.h"
38 #include "gfreduce-exp.h"
39 #include "fibrand.h"
40 #include "mprand.h"
41
42 /*----- Data structures ---------------------------------------------------*/
43
44 DA_DECL(instr_v, gfreduce_instr);
45
46 /*----- Main code ---------------------------------------------------------*/
47
48 /* --- What's going on here? --- *
49  *
50  * Let's face it, @gfx_div@ sucks.  It works (I hope), but it's not in any
51  * sense fast.  Here, we do efficient reduction modulo sparse polynomials.
52  *
53  * Suppose we have a polynomial @X@ we're trying to reduce mod @P@.  If we
54  * take the topmost nonzero word of @X@, call it @w@, then we can eliminate
55  * it by subtracting off @w P x^{k}@ for an appropriate value of @k@.  The
56  * trick is in observing that if @P@ is sparse we can do this multiplication
57  * and subtraction efficiently, just by XORing appropriate shifts of @w@ into
58  * @X@.
59  *
60  * The first tricky bit is in working out when to stop.  I'll use eight-bit
61  * words to demonstrate what I'm talking about.
62  *
63  *  xxxxxx xxxxxxxx xxxxxxxx xxxxxxxx xxxxxxxx xxxxxxxx
64  *                  001ppppp pppppppp pppppppp pppppppp
65  *                    |<rp>|
66  *                    |<------------ bp ------------->|
67  *                  |<------------ nw --------------->|
68  *
69  * The trick of taking whole words off of @X@ stops working when there are
70  * only @nw@ words left.  Then we have to mask off the bottom bits of @w@
71  * before continuing.
72  */
73
74 /* --- @gfreduce_create@ --- *
75  *
76  * Arguments:   @gfreduce *r@ = structure to fill in
77  *              @mp *x@ = a (hopefully sparse) polynomial
78  *
79  * Returns:     ---
80  *
81  * Use:         Initializes a context structure for reduction.
82  */
83
84 void gfreduce_create(gfreduce *r, mp *p)
85 {
86   instr_v iv = DA_INIT;
87   unsigned long d;
88   unsigned dw;
89   mpscan sc;
90   unsigned long i;
91   gfreduce_instr *ip;
92   unsigned f = 0;
93   size_t w, ww, wi, wl, ll, bb;
94
95   /* --- Sort out the easy stuff --- */
96
97   d = mp_bits(p); assert(d); d--;
98   r->lim = d/MPW_BITS;
99   dw = d%MPW_BITS;
100   if (!dw)
101     r->mask = 0;
102   else {
103     r->mask = MPW(((mpw)-1) << dw);
104     r->lim++;
105   }
106   r->p = mp_copy(p);
107
108   /* --- Stash a new instruction --- */
109
110 #define INSTR(op_, arg_) do {                                           \
111   DA_ENSURE(&iv, 1);                                                    \
112   DA(&iv)[DA_LEN(&iv)].op = (op_);                                      \
113   DA(&iv)[DA_LEN(&iv)].arg = (arg_);                                    \
114   DA_EXTEND(&iv, 1);                                                    \
115 } while (0)
116
117 #define f_lsr 1u
118
119   w = (d + MPW_BITS - 1)/MPW_BITS;
120   INSTR(GFRI_LOAD, w);
121   wi = DA_LEN(&iv);
122   f = 0;
123   ll = 0;
124   bb = MPW_BITS - dw;
125   for (i = 0, mp_scan(&sc, p); mp_step(&sc) && i < d; i++) {
126     if (!mp_bit(&sc))
127       continue;
128     ww = (d - i + MPW_BITS - 1)/MPW_BITS;
129     if (ww != w) {
130       wl = DA_LEN(&iv);
131       INSTR(GFRI_STORE, w);
132       if (!ll)
133         ll = DA_LEN(&iv);
134       if (!(f & f_lsr))
135         INSTR(GFRI_LOAD, ww);
136       else {
137         INSTR(GFRI_LOAD, w - 1);
138         for (; wi < wl; wi++) {
139           ip = &DA(&iv)[wi];
140           assert(ip->op == GFRI_LSL);
141           if (ip->arg)
142             INSTR(GFRI_LSR, MPW_BITS - ip->arg);
143         }
144         if (w - 1 != ww) {
145           INSTR(GFRI_STORE, w - 1);
146           INSTR(GFRI_LOAD, ww);
147         }
148         f &= ~f_lsr;
149       }
150       w = ww;
151       wi = DA_LEN(&iv);
152     }
153     INSTR(GFRI_LSL, (bb + i)%MPW_BITS);
154     if ((bb + i)%MPW_BITS)
155       f |= f_lsr;
156   }
157   wl = DA_LEN(&iv);
158   INSTR(GFRI_STORE, w);
159   if (!ll)
160     ll = DA_LEN(&iv);
161   if (f & f_lsr) {
162     INSTR(GFRI_LOAD, w - 1);
163     for (; wi < wl; wi++) {
164       ip = &DA(&iv)[wi];
165       assert(ip->op == GFRI_LSL);
166       if (ip->arg)
167         INSTR(GFRI_LSR, MPW_BITS - ip->arg);
168     }
169     INSTR(GFRI_STORE, w - 1);
170   }
171
172 #undef INSTR
173
174   r->in = DA_LEN(&iv);
175   r->iv = xmalloc(r->in * sizeof(gfreduce_instr));
176   r->liv = r->iv + ll;
177   memcpy(r->iv, DA(&iv), r->in * sizeof(gfreduce_instr));
178   DA_DESTROY(&iv);
179 }
180
181 /* --- @gfreduce_destroy@ --- *
182  *
183  * Arguments:   @gfreduce *r@ = structure to free
184  *
185  * Returns:     ---
186  *
187  * Use:         Reclaims the resources from a reduction context.
188  */
189
190 void gfreduce_destroy(gfreduce *r)
191 {
192   mp_drop(r->p);
193   xfree(r->iv);
194 }
195
196 /* --- @gfreduce_dump@ --- *
197  *
198  * Arguments:   @gfreduce *r@ = structure to dump
199  *              @FILE *fp@ = file to dump on
200  *
201  * Returns:     ---
202  *
203  * Use:         Dumps a reduction context.
204  */
205
206 void gfreduce_dump(gfreduce *r, FILE *fp)
207 {
208   size_t i;
209
210   fprintf(fp, "poly = "); mp_writefile(r->p, fp, 16);
211   fprintf(fp, "\n  lim = %lu; mask = %lx\n",
212           (unsigned long)r->lim, (unsigned long)r->mask);
213   for (i = 0; i < r->in; i++) {
214     static const char *opname[] = { "load", "lsl", "lsr", "store" };
215     assert(r->iv[i].op < N(opname));
216     fprintf(fp, "  %s %lu\n",
217             opname[r->iv[i].op],
218             (unsigned long)r->iv[i].arg);
219   }
220 }
221
222 /* --- @gfreduce_do@ --- *
223  *
224  * Arguments:   @gfreduce *r@ = reduction context
225  *              @mp *d@ = destination
226  *              @mp *x@ = source
227  *
228  * Returns:     Destination, @x@ reduced modulo the reduction poly.
229  */
230
231 static void run(const gfreduce_instr *i, const gfreduce_instr *il,
232                 mpw *v, mpw z)
233 {
234   mpw w = 0;
235
236   for (; i < il; i++) {
237     switch (i->op) {
238       case GFRI_LOAD: w = *(v - i->arg); break;
239       case GFRI_LSL: w ^= z << i->arg; break;
240       case GFRI_LSR: w ^= z >> i->arg; break;
241       case GFRI_STORE: *(v - i->arg) = MPW(w); break;
242       default: abort();
243     }
244   }
245 }
246
247 mp *gfreduce_do(gfreduce *r, mp *d, mp *x)
248 {
249   mpw *v, *vl;
250   const gfreduce_instr *il;
251   mpw z;
252
253   /* --- Try to reuse the source's space --- */
254
255   MP_COPY(x);
256   if (d) MP_DROP(d);
257   MP_DEST(x, MP_LEN(x), x->f);
258
259   /* --- Do the reduction --- */
260
261   il = r->iv + r->in;
262   if (MP_LEN(x) >= r->lim) {
263     v = x->v + r->lim;
264     vl = x->vl;
265     while (vl-- > v) {
266       while (*vl) {
267         z = *vl;
268         *vl = 0;
269         run(r->iv, il, vl, z);
270       }
271     }
272     if (r->mask) {
273       while (*vl & r->mask) {
274         z = *vl & r->mask;
275         *vl &= ~r->mask;
276         run(r->iv, il, vl, z);
277       }
278     }
279   }
280
281   /* --- Done --- */
282
283   MP_SHRINK(x);
284   return (x);
285 }
286
287 /* --- @gfreduce_sqrt@ --- *
288  *
289  * Arguments:   @gfreduce *r@ = pointer to reduction context
290  *              @mp *d@ = destination
291  *              @mp *x@ = some polynomial
292  *
293  * Returns:     The square root of @x@ modulo @r->p@, or null.
294  */
295
296 mp *gfreduce_sqrt(gfreduce *r, mp *d, mp *x)
297 {
298   mp *y = MP_COPY(x);
299   mp *z, *spare = MP_NEW;
300   unsigned long m = mp_bits(r->p) - 1;
301   unsigned long i;
302
303   for (i = 0; i < m - 1; i++) {
304     mp *t = gf_sqr(spare, y);
305     spare = y;
306     y = gfreduce_do(r, t, t);
307   }
308   z = gf_sqr(spare, y);
309   z = gfreduce_do(r, z, z);
310   if (!MP_EQ(x, z)) {
311     mp_drop(y);
312     y = 0;
313   }
314   mp_drop(z);
315   mp_drop(d);
316   return (y);
317 }
318
319 /* --- @gfreduce_trace@ --- *
320  *
321  * Arguments:   @gfreduce *r@ = pointer to reduction context
322  *              @mp *x@ = some polynomial
323  *
324  * Returns:     The trace of @x@. (%$\Tr(x)=x + x^2 + \cdots + x^{2^{m-1}}$%
325  *              if %$x \in \gf{2^m}$%).
326  */
327
328 int gfreduce_trace(gfreduce *r, mp *x)
329 {
330   mp *y = MP_COPY(x);
331   mp *spare = MP_NEW;
332   unsigned long m = mp_bits(r->p) - 1;
333   unsigned long i;
334   int rc;
335
336   for (i = 0; i < m - 1; i++) {
337     mp *t = gf_sqr(spare, y);
338     spare = y;
339     y = gfreduce_do(r, t, t);
340     y = gf_add(y, y, x);
341   }
342   rc = !MP_ZEROP(y);
343   mp_drop(spare);
344   mp_drop(y);
345   return (rc);
346 }
347
348 /* --- @gfreduce_halftrace@ --- *
349  *
350  * Arguments:   @gfreduce *r@ = pointer to reduction context
351  *              @mp *d@ = destination
352  *              @mp *x@ = some polynomial
353  *
354  * Returns:     The half-trace of @x@.
355  *              (%$\HfTr(x)= x + x^{2^2} + \cdots + x^{2^{m-1}}$%
356  *              if %$x \in \gf{2^m}$% with %$m$% odd).
357  */
358
359 mp *gfreduce_halftrace(gfreduce *r, mp *d, mp *x)
360 {
361   mp *y = MP_COPY(x);
362   mp *spare = MP_NEW;
363   unsigned long m = mp_bits(r->p) - 1;
364   unsigned long i;
365
366   mp_drop(d);
367   for (i = 0; i < m - 1; i += 2) {
368     mp *t = gf_sqr(spare, y);
369     spare = y;
370     y = gfreduce_do(r, t, t);
371     t = gf_sqr(spare, y);
372     spare = y;
373     y = gfreduce_do(r, t, t);
374     y = gf_add(y, y, x);
375   }
376   mp_drop(spare);
377   return (y);
378 }
379
380 /* --- @gfreduce_quadsolve@ --- *
381  *
382  * Arguments:   @gfreduce *r@ = pointer to reduction context
383  *              @mp *d@ = destination
384  *              @mp *x@ = some polynomial
385  *
386  * Returns:     A polynomial @y@ such that %$y^2 + y = x$%, or null.
387  */
388
389 mp *gfreduce_quadsolve(gfreduce *r, mp *d, mp *x)
390 {
391   unsigned long m = mp_bits(r->p) - 1;
392   mp *t;
393
394   MP_COPY(x);
395   if (m & 1)
396     d = gfreduce_halftrace(r, d, x);
397   else {
398     mp *z, *w, *rho = MP_NEW;
399     mp *spare = MP_NEW;
400     grand *fr = fibrand_create(0);
401     unsigned long i;
402
403     for (;;) {
404       rho = mprand(rho, m, fr, 0);
405       z = MP_ZERO;
406       w = MP_COPY(rho);
407       for (i = 0; i < m - 1; i++) {
408         t = gf_sqr(spare, z); spare = z; z = gfreduce_do(r, t, t);
409         t = gf_sqr(spare, w); spare = w; w = gfreduce_do(r, t, t);
410         t = gf_mul(spare, w, x); t = gfreduce_do(r, t, t); spare = t;
411         z = gf_add(z, z, t);
412         w = gf_add(w, w, rho);
413       }
414       if (!MP_ZEROP(w))
415         break;
416       MP_DROP(z);
417       MP_DROP(w);
418     }
419     if (d) MP_DROP(d);
420     MP_DROP(w);
421     MP_DROP(spare);
422     MP_DROP(rho);
423     fr->ops->destroy(fr);
424     d = z;
425   }
426
427   t = gf_sqr(MP_NEW, d); t = gfreduce_do(r, t, t); t = gf_add(t, t, d);
428   if (!MP_EQ(t, x)) {
429     MP_DROP(d);
430     d = 0;
431   }
432   MP_DROP(t);
433   MP_DROP(x);
434   if (d) d->v[0] &= ~(mpw)1;
435   return (d);
436 }
437
438 /* --- @gfreduce_exp@ --- *
439  *
440  * Arguments:   @gfreduce *gr@ = pointer to reduction context
441  *              @mp *d@ = fake destination
442  *              @mp *a@ = base
443  *              @mp *e@ = exponent
444  *
445  * Returns:     Result, %$a^e \bmod m$%.
446  */
447
448 mp *gfreduce_exp(gfreduce *gr, mp *d, mp *a, mp *e)
449 {
450   mp *x = MP_ONE;
451   mp *spare = (e->f & MP_BURN) ? MP_NEWSEC : MP_NEW;
452
453   MP_SHRINK(e);
454   MP_COPY(a);
455   if (MP_ZEROP(e))
456     ;
457   else {
458     if (MP_NEGP(e))
459       a = gf_modinv(a, a, gr->p);
460     if (MP_LEN(e) < EXP_THRESH)
461       EXP_SIMPLE(x, a, e);
462     else
463       EXP_WINDOW(x, a, e);
464   }
465   mp_drop(d);
466   mp_drop(a);
467   mp_drop(spare);
468   return (x);
469 }
470
471 /*----- Test rig ----------------------------------------------------------*/
472
473 #ifdef TEST_RIG
474
475 #define MP(x) mp_readstring(MP_NEW, #x, 0, 0)
476
477 static int vreduce(dstr *v)
478 {
479   mp *d = *(mp **)v[0].buf;
480   mp *n = *(mp **)v[1].buf;
481   mp *r = *(mp **)v[2].buf;
482   mp *c;
483   int ok = 1;
484   gfreduce rr;
485
486   gfreduce_create(&rr, d);
487   c = gfreduce_do(&rr, MP_NEW, n);
488   if (!MP_EQ(c, r)) {
489     fprintf(stderr, "\n*** reduction failed\n*** ");
490     gfreduce_dump(&rr, stderr);
491     fprintf(stderr, "\n*** n = "); mp_writefile(n, stderr, 16);
492     fprintf(stderr, "\n*** r = "); mp_writefile(r, stderr, 16);
493     fprintf(stderr, "\n*** c = "); mp_writefile(c, stderr, 16);
494     fprintf(stderr, "\n");
495     ok = 0;
496   }
497   gfreduce_destroy(&rr);
498   mp_drop(n); mp_drop(d); mp_drop(r); mp_drop(c);
499   assert(mparena_count(MPARENA_GLOBAL) == 0);
500   return (ok);
501 }
502
503 static int vmodexp(dstr *v)
504 {
505   mp *p = *(mp **)v[0].buf;
506   mp *g = *(mp **)v[1].buf;
507   mp *x = *(mp **)v[2].buf;
508   mp *r = *(mp **)v[3].buf;
509   mp *c;
510   int ok = 1;
511   gfreduce rr;
512
513   gfreduce_create(&rr, p);
514   c = gfreduce_exp(&rr, MP_NEW, g, x);
515   if (!MP_EQ(c, r)) {
516     fprintf(stderr, "\n*** modexp failed\n*** ");
517     fprintf(stderr, "\n*** p = "); mp_writefile(p, stderr, 16);
518     fprintf(stderr, "\n*** g = "); mp_writefile(g, stderr, 16);
519     fprintf(stderr, "\n*** x = "); mp_writefile(x, stderr, 16);
520     fprintf(stderr, "\n*** c = "); mp_writefile(c, stderr, 16);
521     fprintf(stderr, "\n*** r = "); mp_writefile(r, stderr, 16);
522     fprintf(stderr, "\n");
523     ok = 0;
524   }
525   gfreduce_destroy(&rr);
526   mp_drop(p); mp_drop(g); mp_drop(r); mp_drop(x); mp_drop(c);
527   assert(mparena_count(MPARENA_GLOBAL) == 0);
528   return (ok);
529 }
530
531 static int vsqrt(dstr *v)
532 {
533   mp *p = *(mp **)v[0].buf;
534   mp *x = *(mp **)v[1].buf;
535   mp *r = *(mp **)v[2].buf;
536   mp *c;
537   int ok = 1;
538   gfreduce rr;
539
540   gfreduce_create(&rr, p);
541   c = gfreduce_sqrt(&rr, MP_NEW, x);
542   if (!MP_EQ(c, r)) {
543     fprintf(stderr, "\n*** sqrt failed\n*** ");
544     fprintf(stderr, "\n*** p = "); mp_writefile(p, stderr, 16);
545     fprintf(stderr, "\n*** x = "); mp_writefile(x, stderr, 16);
546     fprintf(stderr, "\n*** c = "); mp_writefile(c, stderr, 16);
547     fprintf(stderr, "\n*** r = "); mp_writefile(r, stderr, 16);
548     fprintf(stderr, "\n");
549     ok = 0;
550   }
551   gfreduce_destroy(&rr);
552   mp_drop(p); mp_drop(r); mp_drop(x); mp_drop(c);
553   assert(mparena_count(MPARENA_GLOBAL) == 0);
554   return (ok);
555 }
556
557 static int vtr(dstr *v)
558 {
559   mp *p = *(mp **)v[0].buf;
560   mp *x = *(mp **)v[1].buf;
561   int r = *(int *)v[2].buf, c;
562   int ok = 1;
563   gfreduce rr;
564
565   gfreduce_create(&rr, p);
566   c = gfreduce_trace(&rr, x);
567   if (c != r) {
568     fprintf(stderr, "\n*** trace failed\n*** ");
569     fprintf(stderr, "\n*** p = "); mp_writefile(p, stderr, 16);
570     fprintf(stderr, "\n*** x = "); mp_writefile(x, stderr, 16);
571     fprintf(stderr, "\n*** c = %d", c);
572     fprintf(stderr, "\n*** r = %d", r);
573     fprintf(stderr, "\n");
574     ok = 0;
575   }
576   gfreduce_destroy(&rr);
577   mp_drop(p); mp_drop(x);
578   assert(mparena_count(MPARENA_GLOBAL) == 0);
579   return (ok);
580 }
581
582 static int vhftr(dstr *v)
583 {
584   mp *p = *(mp **)v[0].buf;
585   mp *x = *(mp **)v[1].buf;
586   mp *r = *(mp **)v[2].buf;
587   mp *c;
588   int ok = 1;
589   gfreduce rr;
590
591   gfreduce_create(&rr, p);
592   c = gfreduce_halftrace(&rr, MP_NEW, x);
593   if (!MP_EQ(c, r)) {
594     fprintf(stderr, "\n*** halftrace failed\n*** ");
595     fprintf(stderr, "\n*** p = "); mp_writefile(p, stderr, 16);
596     fprintf(stderr, "\n*** x = "); mp_writefile(x, stderr, 16);
597     fprintf(stderr, "\n*** c = "); mp_writefile(c, stderr, 16);
598     fprintf(stderr, "\n*** r = "); mp_writefile(r, stderr, 16);
599     fprintf(stderr, "\n");
600     ok = 0;
601   }
602   gfreduce_destroy(&rr);
603   mp_drop(p); mp_drop(r); mp_drop(x); mp_drop(c);
604   assert(mparena_count(MPARENA_GLOBAL) == 0);
605   return (ok);
606 }
607
608 static int vquad(dstr *v)
609 {
610   mp *p = *(mp **)v[0].buf;
611   mp *x = *(mp **)v[1].buf;
612   mp *r = *(mp **)v[2].buf;
613   mp *c;
614   int ok = 1;
615   gfreduce rr;
616
617   gfreduce_create(&rr, p);
618   c = gfreduce_quadsolve(&rr, MP_NEW, x);
619   if (!MP_EQ(c, r)) {
620     fprintf(stderr, "\n*** quadsolve failed\n*** ");
621     fprintf(stderr, "\n*** p = "); mp_writefile(p, stderr, 16);
622     fprintf(stderr, "\n*** x = "); mp_writefile(x, stderr, 16);
623     fprintf(stderr, "\n*** c = "); mp_writefile(c, stderr, 16);
624     fprintf(stderr, "\n*** r = "); mp_writefile(r, stderr, 16);
625     fprintf(stderr, "\n");
626     ok = 0;
627   }
628   gfreduce_destroy(&rr);
629   mp_drop(p); mp_drop(r); mp_drop(x); mp_drop(c);
630   assert(mparena_count(MPARENA_GLOBAL) == 0);
631   return (ok);
632 }
633
634 static test_chunk defs[] = {
635   { "reduce", vreduce, { &type_mp, &type_mp, &type_mp, 0 } },
636   { "modexp", vmodexp, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
637   { "sqrt", vsqrt, { &type_mp, &type_mp, &type_mp, 0 } },
638   { "trace", vtr, { &type_mp, &type_mp, &type_int, 0 } },
639   { "halftrace", vhftr, { &type_mp, &type_mp, &type_mp, 0 } },
640   { "quadsolve", vquad, { &type_mp, &type_mp, &type_mp, 0 } },
641   { 0, 0, { 0 } }
642 };
643
644 int main(int argc, char *argv[])
645 {
646   test_run(argc, argv, defs, SRCDIR"/tests/gfreduce");
647   return (0);
648 }
649
650 #endif
651
652 /*----- That's all, folks -------------------------------------------------*/