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