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

   1 /* -*-c-*-
   2  *
   3  * Efficient reduction modulo sparse binary polynomials
   4  *
   5  * (c) 2004 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 <mLib/alloc.h>
  31 #include <mLib/darray.h>
  32 #include <mLib/macros.h>
  33
  34 #include "gf.h"
  35 #include "gfreduce.h"
  36 #include "gfreduce-exp.h"
  37 #include "fibrand.h"
  38 #include "mprand.h"
  39
  40 /*----- Data structures ---------------------------------------------------*/
  41
  42 DA_DECL(instr_v, gfreduce_instr);
  43
  44 /*----- Main code ---------------------------------------------------------*/
  45
  46 /* --- What's going on here? --- *
  47  *
  48  * Let's face it, @gfx_div@ sucks.  It works (I hope), but it's not in any
  49  * sense fast.  Here, we do efficient reduction modulo sparse polynomials.
  50  * (It works for arbitrary polynomials, but isn't efficient for dense ones.)
  51  *
  52  * Suppose that %$p(x) = x^n + p'(x) = \sum_{0\le i<n} p_i x^i$%, hopefully
  53  * with only a few other %$p_i \ne 0$%.  We're going to compile %$p$% into a
  54  * sequence of instructions which can be used to perform reduction modulo
  55  * %$p$%.  The important observation is that %$x^n \equiv p' \pmod p$%.
  56  *
  57  * Suppose we're working with %$w$%-bit words; let %$n = N w + n'$% with
  58  * %$0 \le n' < w$%.  Let %$u(x)$% be some arbitrary polynomial.  Write
  59  * %$u = z x^k + u'$% with %$\deg u' < k \ge n$%; then a reduction step uses
  60  * that %$u \equiv u' + z p' x^{k-n} \pmod p$%: the right hand side has
  61  * degree %$\max \{ \deg u', k + \deg p' - n + \deg z \} < \deg u$%, so this
  62  * makes progress towards a complete reduction.
  63  *
  64  * The compiled instruction sequence computes
  65  * %$u' + z p' x^{k-n} = u' + \sum_{0\le i<n} z x^{k-n+i}$%.
  66  */
  67
  68 /* --- @gfreduce_create@ --- *
  69  *
  70  * Arguments:   @gfreduce *r@ = structure to fill in
  71  *              @mp *x@ = a (hopefully sparse) polynomial
  72  *
  73  * Returns:     ---
  74  *
  75  * Use:         Initializes a context structure for reduction.
  76  */
  77
  78 struct gen {
  79   unsigned f;                           /* Flags */
  80 #define f_lsr 1u                        /*   Overflow from previous word */
  81 #define f_load 2u                       /*   Outstanding @LOAD@ */
  82   instr_v iv;                           /* Instruction vector */
  83   size_t w;                             /* Currently loaded target word */
  84   size_t wi;                            /* Left-shifts for current word */
  85 };
  86
  87 #define INSTR(g_, op_, arg_) do {                                       \
  88   struct gen *_g = (g_);                                                \
  89   instr_v *_iv = &_g->iv;                                               \
  90   size_t _i = DA_LEN(_iv);                                              \
  91                                                                         \
  92   DA_ENSURE(_iv, 1);                                                    \
  93   DA(_iv)[_i].op = (op_);                                               \
  94   DA(_iv)[_i].arg = (arg_);                                             \
  95   DA_EXTEND(_iv, 1);                                                    \
  96 } while (0)
  97
  98 static void emit_load(struct gen *g, size_t w)
  99 {
 100   INSTR(g, GFRI_LOAD, w);
 101   g->f |= f_load;
 102   g->w = w;
 103 }
 104
 105 static void emit_right_shifts(struct gen *g)
 106 {
 107   gfreduce_instr *ip;
 108   size_t i, wl;
 109
 110   /* --- Close off the current word --- *
 111    *
 112    * If we shifted into this current word with a nonzero bit offset, then
 113    * we'll also need to arrange to perform a sequence of right shifts into
 114    * the following word, which we might as well do by scanning the
 115    * instruction sequence (which starts at @wi@).
 116    *
 117    * Either way, we leave a @LOAD@ unmatched if there was one before, in the
 118    * hope that callers have an easier time; @g->w@ is updated to reflect the
 119    * currently open word.
 120    */
 121
 122   if (!(g->f & f_lsr))
 123     return;
 124
 125   wl = DA_LEN(&g->iv);
 126   INSTR(g, GFRI_STORE, g->w);
 127   emit_load(g, g->w - 1);
 128   for (i = g->wi; i < wl; i++) {
 129     ip = &DA(&g->iv)[i];
 130     assert(ip->op == GFRI_LSL);
 131     if (ip->arg)
 132       INSTR(g, GFRI_LSR, MPW_BITS - ip->arg);
 133   }
 134   g->f &= ~f_lsr;
 135 }
 136
 137 static void ensure_loaded(struct gen *g, size_t w)
 138 {
 139   if (!(g->f & f_load)) {
 140     emit_load(g, w);
 141     g->wi = DA_LEN(&g->iv);
 142   } else if (w != g->w) {
 143     emit_right_shifts(g);
 144     if (w != g->w) {
 145       INSTR(g, GFRI_STORE, g->w);
 146       emit_load(g, w);
 147     }
 148     g->wi = DA_LEN(&g->iv);
 149   }
 150 }
 151
 152 void gfreduce_create(gfreduce *r, mp *p)
 153 {
 154   struct gen g = { 0, DA_INIT };
 155   unsigned long d;
 156   unsigned dw;
 157   mpscan sc;
 158   unsigned long i;
 159   size_t w, bb;
 160
 161   /* --- Sort out the easy stuff --- */
 162
 163   d = mp_bits(p); assert(d); d--;
 164   r->lim = d/MPW_BITS;
 165   dw = d%MPW_BITS;
 166   if (!dw)
 167     r->mask = 0;
 168   else {
 169     r->mask = MPW(((mpw)-1) << dw);
 170     r->lim++;
 171   }
 172   r->p = mp_copy(p);
 173
 174   /* --- How this works --- *
 175    *
 176    * The instruction sequence is run with two ambient parameters: a pointer
 177    * (usually) just past the most significant word of the polynomial to be
 178    * reduced; and a word %$z$% which is the multiple of %$p'$% we are meant
 179    * to add.
 180    *
 181    * The sequence visits each word of the polynomial at most once.  Suppose
 182    * %$u = z x^{w N} + u'$%; our pointer points just past the end of %$u'$%.
 183    * Word %$I$% of %$u'$% will be affected by modulus bits %$p_i$% where
 184    * %$(N - I - 1) w + 1 \le i \le (N - I + 1) w - 1$%, so %$p_i$% affects
 185    * word %$I = \lceil (n - i + 1)/w \rceil$% and (if %$i$% is not a multiple
 186    * of %$w$%) also word %$I - 1$%.
 187    *
 188    * We have four instructions: @LOAD@ reads a specified word of %$u$% into an
 189    * accumulator, and @STORE@ stores it back (we'll always store back to the
 190    * same word we most recently read, but this isn't a requirement); and
 191    * @LSL@ and @LSR@, which XOR in appropriately shifted copies of %$z$% into
 192    * the accumulator.  So a typical program will contain sequences of @LSR@
 193    * and @LSL@ instructions sandwiched between @LOAD@/@STORE@ pairs.
 194    *
 195    * We do a single right-to-left pass across %$p$%.
 196    */
 197
 198   bb = MPW_BITS - dw;
 199
 200   for (i = 0, mp_scan(&sc, p); mp_step(&sc) && i < d; i++) {
 201     if (!mp_bit(&sc))
 202       continue;
 203
 204     /* --- We've found a set bit, so work out which word it affects --- *
 205      *
 206      * In general, a bit affects two words: it needs to be shifted left into
 207      * one, and shifted right into the next.  We find the former here.
 208      */
 209
 210     w = (d - i + MPW_BITS - 1)/MPW_BITS;
 211
 212     /* --- Concentrate on the appropriate word --- */
 213
 214     ensure_loaded(&g, w);
 215
 216     /* --- Accumulate a new @LSL@ instruction --- *
 217      *
 218      * If this was a nonzero shift, then we'll need to arrange to do right
 219      * shifts into the following word.
 220      */
 221
 222     INSTR(&g, GFRI_LSL, (bb + i)%MPW_BITS);
 223     if ((bb + i)%MPW_BITS)
 224       g.f |= f_lsr;
 225   }
 226
 227   /* --- Wrapping up --- *
 228    *
 229    * We probably need a final @STORE@, and maybe a sequence of right shifts.
 230    */
 231
 232   if (g.f & f_load) {
 233     emit_right_shifts(&g);
 234     INSTR(&g, GFRI_STORE, g.w);
 235   }
 236
 237   r->in = DA_LEN(&g.iv);
 238   r->iv = xmalloc(r->in * sizeof(gfreduce_instr));
 239   memcpy(r->iv, DA(&g.iv), r->in * sizeof(gfreduce_instr));
 240   DA_DESTROY(&g.iv);
 241 }
 242
 243 #undef INSTR
 244
 245 #undef f_lsr
 246 #undef f_load
 247
 248 /* --- @gfreduce_destroy@ --- *
 249  *
 250  * Arguments:   @gfreduce *r@ = structure to free
 251  *
 252  * Returns:     ---
 253  *
 254  * Use:         Reclaims the resources from a reduction context.
 255  */
 256
 257 void gfreduce_destroy(gfreduce *r)
 258 {
 259   mp_drop(r->p);
 260   xfree(r->iv);
 261 }
 262
 263 /* --- @gfreduce_dump@ --- *
 264  *
 265  * Arguments:   @gfreduce *r@ = structure to dump
 266  *              @FILE *fp@ = file to dump on
 267  *
 268  * Returns:     ---
 269  *
 270  * Use:         Dumps a reduction context.
 271  */
 272
 273 void gfreduce_dump(gfreduce *r, FILE *fp)
 274 {
 275   size_t i;
 276
 277   fprintf(fp, "poly = "); mp_writefile(r->p, fp, 16);
 278   fprintf(fp, "\n  lim = %lu; mask = %lx\n",
 279           (unsigned long)r->lim, (unsigned long)r->mask);
 280   for (i = 0; i < r->in; i++) {
 281     static const char *opname[] = { "load", "lsl", "lsr", "store" };
 282     assert(r->iv[i].op < N(opname));
 283     fprintf(fp, "  %s %lu\n",
 284             opname[r->iv[i].op],
 285             (unsigned long)r->iv[i].arg);
 286   }
 287 }
 288
 289 /* --- @gfreduce_do@ --- *
 290  *
 291  * Arguments:   @gfreduce *r@ = reduction context
 292  *              @mp *d@ = destination
 293  *              @mp *x@ = source
 294  *
 295  * Returns:     Destination, @x@ reduced modulo the reduction poly.
 296  */
 297
 298 static void run(const gfreduce_instr *i, const gfreduce_instr *il,
 299                 mpw *v, mpw z)
 300 {
 301   mpw w = 0;
 302
 303   for (; i < il; i++) {
 304     switch (i->op) {
 305       case GFRI_LOAD: w = *(v - i->arg); break;
 306       case GFRI_LSL: w ^= z << i->arg; break;
 307       case GFRI_LSR: w ^= z >> i->arg; break;
 308       case GFRI_STORE: *(v - i->arg) = MPW(w); break;
 309       default: abort();
 310     }
 311   }
 312 }
 313
 314 mp *gfreduce_do(gfreduce *r, mp *d, mp *x)
 315 {
 316   mpw *v, *vl;
 317   const gfreduce_instr *il;
 318   mpw z;
 319
 320   /* --- Try to reuse the source's space --- */
 321
 322   MP_COPY(x);
 323   if (d) MP_DROP(d);
 324   MP_DEST(x, MP_LEN(x), x->f);
 325
 326   /* --- Do the reduction --- */
 327
 328   il = r->iv + r->in;
 329   if (MP_LEN(x) >= r->lim) {
 330     v = x->v + r->lim;
 331     vl = x->vl;
 332     while (vl-- > v) {
 333       while (*vl) {
 334         z = *vl;
 335         *vl = 0;
 336         run(r->iv, il, vl, z);
 337       }
 338     }
 339     if (r->mask) {
 340       while (*vl & r->mask) {
 341         z = *vl & r->mask;
 342         *vl &= ~r->mask;
 343         run(r->iv, il, vl, z);
 344       }
 345     }
 346   }
 347
 348   /* --- Done --- */
 349
 350   MP_SHRINK(x);
 351   return (x);
 352 }
 353
 354 /* --- @gfreduce_sqrt@ --- *
 355  *
 356  * Arguments:   @gfreduce *r@ = pointer to reduction context
 357  *              @mp *d@ = destination
 358  *              @mp *x@ = some polynomial
 359  *
 360  * Returns:     The square root of @x@ modulo @r->p@, or null.
 361  */
 362
 363 mp *gfreduce_sqrt(gfreduce *r, mp *d, mp *x)
 364 {
 365   mp *y = MP_COPY(x);
 366   mp *z, *spare = MP_NEW;
 367   unsigned long m = mp_bits(r->p) - 1;
 368   unsigned long i;
 369
 370   for (i = 0; i < m - 1; i++) {
 371     mp *t = gf_sqr(spare, y);
 372     spare = y;
 373     y = gfreduce_do(r, t, t);
 374   }
 375   z = gf_sqr(spare, y);
 376   z = gfreduce_do(r, z, z);
 377   if (!MP_EQ(x, z)) {
 378     mp_drop(y);
 379     y = 0;
 380   }
 381   mp_drop(z);
 382   mp_drop(d);
 383   return (y);
 384 }
 385
 386 /* --- @gfreduce_trace@ --- *
 387  *
 388  * Arguments:   @gfreduce *r@ = pointer to reduction context
 389  *              @mp *x@ = some polynomial
 390  *
 391  * Returns:     The trace of @x@. (%$\Tr(x)=x + x^2 + \cdots + x^{2^{m-1}}$%
 392  *              if %$x \in \gf{2^m}$%).
 393  */
 394
 395 int gfreduce_trace(gfreduce *r, mp *x)
 396 {
 397   mp *y = MP_COPY(x);
 398   mp *spare = MP_NEW;
 399   unsigned long m = mp_bits(r->p) - 1;
 400   unsigned long i;
 401   int rc;
 402
 403   for (i = 0; i < m - 1; i++) {
 404     mp *t = gf_sqr(spare, y);
 405     spare = y;
 406     y = gfreduce_do(r, t, t);
 407     y = gf_add(y, y, x);
 408   }
 409   rc = !MP_ZEROP(y);
 410   mp_drop(spare);
 411   mp_drop(y);
 412   return (rc);
 413 }
 414
 415 /* --- @gfreduce_halftrace@ --- *
 416  *
 417  * Arguments:   @gfreduce *r@ = pointer to reduction context
 418  *              @mp *d@ = destination
 419  *              @mp *x@ = some polynomial
 420  *
 421  * Returns:     The half-trace of @x@.
 422  *              (%$\HfTr(x)= x + x^{2^2} + \cdots + x^{2^{m-1}}$%
 423  *              if %$x \in \gf{2^m}$% with %$m$% odd).
 424  */
 425
 426 mp *gfreduce_halftrace(gfreduce *r, mp *d, mp *x)
 427 {
 428   mp *y = MP_COPY(x);
 429   mp *spare = MP_NEW;
 430   unsigned long m = mp_bits(r->p) - 1;
 431   unsigned long i;
 432
 433   mp_drop(d);
 434   for (i = 0; i < m - 1; i += 2) {
 435     mp *t = gf_sqr(spare, y);
 436     spare = y;
 437     y = gfreduce_do(r, t, t);
 438     t = gf_sqr(spare, y);
 439     spare = y;
 440     y = gfreduce_do(r, t, t);
 441     y = gf_add(y, y, x);
 442   }
 443   mp_drop(spare);
 444   return (y);
 445 }
 446
 447 /* --- @gfreduce_quadsolve@ --- *
 448  *
 449  * Arguments:   @gfreduce *r@ = pointer to reduction context
 450  *              @mp *d@ = destination
 451  *              @mp *x@ = some polynomial
 452  *
 453  * Returns:     A polynomial @y@ such that %$y^2 + y = x$%, or null.
 454  */
 455
 456 mp *gfreduce_quadsolve(gfreduce *r, mp *d, mp *x)
 457 {
 458   unsigned long m = mp_bits(r->p) - 1;
 459   mp *t;
 460
 461   MP_COPY(x);
 462   if (m & 1)
 463     d = gfreduce_halftrace(r, d, x);
 464   else {
 465     mp *z, *w, *rho = MP_NEW;
 466     mp *spare = MP_NEW;
 467     grand *fr = fibrand_create(0);
 468     unsigned long i;
 469
 470     for (;;) {
 471       rho = mprand(rho, m, fr, 0);
 472       z = MP_ZERO;
 473       w = MP_COPY(rho);
 474       for (i = 0; i < m - 1; i++) {
 475         t = gf_sqr(spare, z); spare = z; z = gfreduce_do(r, t, t);
 476         t = gf_sqr(spare, w); spare = w; w = gfreduce_do(r, t, t);
 477         t = gf_mul(spare, w, x); t = gfreduce_do(r, t, t); spare = t;
 478         z = gf_add(z, z, t);
 479         w = gf_add(w, w, rho);
 480       }
 481       if (!MP_ZEROP(w))
 482         break;
 483       MP_DROP(z);
 484       MP_DROP(w);
 485     }
 486     if (d) MP_DROP(d);
 487     MP_DROP(w);
 488     MP_DROP(spare);
 489     MP_DROP(rho);
 490     fr->ops->destroy(fr);
 491     d = z;
 492   }
 493
 494   t = gf_sqr(MP_NEW, d); t = gfreduce_do(r, t, t); t = gf_add(t, t, d);
 495   if (!MP_EQ(t, x)) {
 496     MP_DROP(d);
 497     d = 0;
 498   }
 499   MP_DROP(t);
 500   MP_DROP(x);
 501   if (d) d->v[0] &= ~(mpw)1;
 502   return (d);
 503 }
 504
 505 /* --- @gfreduce_exp@ --- *
 506  *
 507  * Arguments:   @gfreduce *gr@ = pointer to reduction context
 508  *              @mp *d@ = fake destination
 509  *              @mp *a@ = base
 510  *              @mp *e@ = exponent
 511  *
 512  * Returns:     Result, %$a^e \bmod m$%.
 513  */
 514
 515 mp *gfreduce_exp(gfreduce *gr, mp *d, mp *a, mp *e)
 516 {
 517   mp *x = MP_ONE;
 518   mp *spare = (e->f & MP_BURN) ? MP_NEWSEC : MP_NEW;
 519
 520   MP_SHRINK(e);
 521   MP_COPY(a);
 522   if (MP_ZEROP(e))
 523     ;
 524   else {
 525     if (MP_NEGP(e))
 526       a = gf_modinv(a, a, gr->p);
 527     if (MP_LEN(e) < EXP_THRESH)
 528       EXP_SIMPLE(x, a, e);
 529     else
 530       EXP_WINDOW(x, a, e);
 531   }
 532   mp_drop(d);
 533   mp_drop(a);
 534   mp_drop(spare);
 535   return (x);
 536 }
 537
 538 /*----- Test rig ----------------------------------------------------------*/
 539
 540 #ifdef TEST_RIG
 541
 542 #define MP(x) mp_readstring(MP_NEW, #x, 0, 0)
 543
 544 static int vreduce(dstr *v)
 545 {
 546   mp *d = *(mp **)v[0].buf;
 547   mp *n = *(mp **)v[1].buf;
 548   mp *r = *(mp **)v[2].buf;
 549   mp *c;
 550   int ok = 1;
 551   gfreduce rr;
 552
 553   gfreduce_create(&rr, d);
 554   c = gfreduce_do(&rr, MP_NEW, n);
 555   if (!MP_EQ(c, r)) {
 556     fprintf(stderr, "\n*** reduction failed\n*** ");
 557     gfreduce_dump(&rr, stderr);
 558     fprintf(stderr, "\n*** n = "); mp_writefile(n, stderr, 16);
 559     fprintf(stderr, "\n*** r = "); mp_writefile(r, stderr, 16);
 560     fprintf(stderr, "\n*** c = "); mp_writefile(c, stderr, 16);
 561     fprintf(stderr, "\n");
 562     ok = 0;
 563   }
 564   gfreduce_destroy(&rr);
 565   mp_drop(n); mp_drop(d); mp_drop(r); mp_drop(c);
 566   assert(mparena_count(MPARENA_GLOBAL) == 0);
 567   return (ok);
 568 }
 569
 570 static int vmodexp(dstr *v)
 571 {
 572   mp *p = *(mp **)v[0].buf;
 573   mp *g = *(mp **)v[1].buf;
 574   mp *x = *(mp **)v[2].buf;
 575   mp *r = *(mp **)v[3].buf;
 576   mp *c;
 577   int ok = 1;
 578   gfreduce rr;
 579
 580   gfreduce_create(&rr, p);
 581   c = gfreduce_exp(&rr, MP_NEW, g, x);
 582   if (!MP_EQ(c, r)) {
 583     fprintf(stderr, "\n*** modexp failed\n*** ");
 584     fprintf(stderr, "\n*** p = "); mp_writefile(p, stderr, 16);
 585     fprintf(stderr, "\n*** g = "); mp_writefile(g, stderr, 16);
 586     fprintf(stderr, "\n*** x = "); mp_writefile(x, stderr, 16);
 587     fprintf(stderr, "\n*** c = "); mp_writefile(c, stderr, 16);
 588     fprintf(stderr, "\n*** r = "); mp_writefile(r, stderr, 16);
 589     fprintf(stderr, "\n");
 590     ok = 0;
 591   }
 592   gfreduce_destroy(&rr);
 593   mp_drop(p); mp_drop(g); mp_drop(r); mp_drop(x); mp_drop(c);
 594   assert(mparena_count(MPARENA_GLOBAL) == 0);
 595   return (ok);
 596 }
 597
 598 static int vsqrt(dstr *v)
 599 {
 600   mp *p = *(mp **)v[0].buf;
 601   mp *x = *(mp **)v[1].buf;
 602   mp *r = *(mp **)v[2].buf;
 603   mp *c;
 604   int ok = 1;
 605   gfreduce rr;
 606
 607   gfreduce_create(&rr, p);
 608   c = gfreduce_sqrt(&rr, MP_NEW, x);
 609   if (!MP_EQ(c, r)) {
 610     fprintf(stderr, "\n*** sqrt failed\n*** ");
 611     fprintf(stderr, "\n*** p = "); mp_writefile(p, stderr, 16);
 612     fprintf(stderr, "\n*** x = "); mp_writefile(x, stderr, 16);
 613     fprintf(stderr, "\n*** c = "); mp_writefile(c, stderr, 16);
 614     fprintf(stderr, "\n*** r = "); mp_writefile(r, stderr, 16);
 615     fprintf(stderr, "\n");
 616     ok = 0;
 617   }
 618   gfreduce_destroy(&rr);
 619   mp_drop(p); mp_drop(r); mp_drop(x); mp_drop(c);
 620   assert(mparena_count(MPARENA_GLOBAL) == 0);
 621   return (ok);
 622 }
 623
 624 static int vtr(dstr *v)
 625 {
 626   mp *p = *(mp **)v[0].buf;
 627   mp *x = *(mp **)v[1].buf;
 628   int r = *(int *)v[2].buf, c;
 629   int ok = 1;
 630   gfreduce rr;
 631
 632   gfreduce_create(&rr, p);
 633   c = gfreduce_trace(&rr, x);
 634   if (c != r) {
 635     fprintf(stderr, "\n*** trace failed\n*** ");
 636     fprintf(stderr, "\n*** p = "); mp_writefile(p, stderr, 16);
 637     fprintf(stderr, "\n*** x = "); mp_writefile(x, stderr, 16);
 638     fprintf(stderr, "\n*** c = %d", c);
 639     fprintf(stderr, "\n*** r = %d", r);
 640     fprintf(stderr, "\n");
 641     ok = 0;
 642   }
 643   gfreduce_destroy(&rr);
 644   mp_drop(p); mp_drop(x);
 645   assert(mparena_count(MPARENA_GLOBAL) == 0);
 646   return (ok);
 647 }
 648
 649 static int vhftr(dstr *v)
 650 {
 651   mp *p = *(mp **)v[0].buf;
 652   mp *x = *(mp **)v[1].buf;
 653   mp *r = *(mp **)v[2].buf;
 654   mp *c;
 655   int ok = 1;
 656   gfreduce rr;
 657
 658   gfreduce_create(&rr, p);
 659   c = gfreduce_halftrace(&rr, MP_NEW, x);
 660   if (!MP_EQ(c, r)) {
 661     fprintf(stderr, "\n*** halftrace failed\n*** ");
 662     fprintf(stderr, "\n*** p = "); mp_writefile(p, stderr, 16);
 663     fprintf(stderr, "\n*** x = "); mp_writefile(x, stderr, 16);
 664     fprintf(stderr, "\n*** c = "); mp_writefile(c, stderr, 16);
 665     fprintf(stderr, "\n*** r = "); mp_writefile(r, stderr, 16);
 666     fprintf(stderr, "\n");
 667     ok = 0;
 668   }
 669   gfreduce_destroy(&rr);
 670   mp_drop(p); mp_drop(r); mp_drop(x); mp_drop(c);
 671   assert(mparena_count(MPARENA_GLOBAL) == 0);
 672   return (ok);
 673 }
 674
 675 static int vquad(dstr *v)
 676 {
 677   mp *p = *(mp **)v[0].buf;
 678   mp *x = *(mp **)v[1].buf;
 679   mp *r = *(mp **)v[2].buf;
 680   mp *c;
 681   int ok = 1;
 682   gfreduce rr;
 683
 684   gfreduce_create(&rr, p);
 685   c = gfreduce_quadsolve(&rr, MP_NEW, x);
 686   if (!MP_EQ(c, r)) {
 687     fprintf(stderr, "\n*** quadsolve failed\n*** ");
 688     fprintf(stderr, "\n*** p = "); mp_writefile(p, stderr, 16);
 689     fprintf(stderr, "\n*** x = "); mp_writefile(x, stderr, 16);
 690     fprintf(stderr, "\n*** c = "); mp_writefile(c, stderr, 16);
 691     fprintf(stderr, "\n*** r = "); mp_writefile(r, stderr, 16);
 692     fprintf(stderr, "\n");
 693     ok = 0;
 694   }
 695   gfreduce_destroy(&rr);
 696   mp_drop(p); mp_drop(r); mp_drop(x); mp_drop(c);
 697   assert(mparena_count(MPARENA_GLOBAL) == 0);
 698   return (ok);
 699 }
 700
 701 static test_chunk defs[] = {
 702   { "reduce", vreduce, { &type_mp, &type_mp, &type_mp, 0 } },
 703   { "modexp", vmodexp, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
 704   { "sqrt", vsqrt, { &type_mp, &type_mp, &type_mp, 0 } },
 705   { "trace", vtr, { &type_mp, &type_mp, &type_int, 0 } },
 706   { "halftrace", vhftr, { &type_mp, &type_mp, &type_mp, 0 } },
 707   { "quadsolve", vquad, { &type_mp, &type_mp, &type_mp, 0 } },
 708   { 0, 0, { 0 } }
 709 };
 710
 711 int main(int argc, char *argv[])
 712 {
 713   test_run(argc, argv, defs, SRCDIR"/t/gfreduce");
 714   return (0);
 715 }
 716
 717 #endif
 718
 719 /*----- That's all, folks -------------------------------------------------*/