Track mLib change: symbols no longer need to include a terminating
[u/mdw/catacomb] / mpmont.c
1 /* -*-c-*-
2 *
3 * $Id: mpmont.c,v 1.12 2000/10/08 15:48:35 mdw Exp $
4 *
5 * Montgomery reduction
6 *
7 * (c) 1999 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: mpmont.c,v $
33 * Revision 1.12 2000/10/08 15:48:35 mdw
34 * Rename Karatsuba constants now that we have @gfx_kmul@ too.
35 *
36 * Revision 1.11 2000/10/08 12:04:27 mdw
37 * (mpmont_reduce, mpmont_mul): Cope with negative numbers.
38 *
39 * Revision 1.10 2000/07/29 17:05:43 mdw
40 * (mpmont_expr): Use sliding window exponentiation, with a drop-through
41 * for small exponents to use a simple left-to-right bitwise routine. This
42 * can reduce modexp times by up to a quarter.
43 *
44 * Revision 1.9 2000/06/17 11:45:09 mdw
45 * Major memory management overhaul. Added arena support. Use the secure
46 * arena for secret integers. Replace and improve the MP management macros
47 * (e.g., replace MP_MODIFY by MP_DEST).
48 *
49 * Revision 1.8 1999/12/22 15:55:00 mdw
50 * Adjust Karatsuba parameters.
51 *
52 * Revision 1.7 1999/12/11 01:51:14 mdw
53 * Use a Karatsuba-based reduction for large moduli.
54 *
55 * Revision 1.6 1999/12/10 23:18:39 mdw
56 * Change interface for suggested destinations.
57 *
58 * Revision 1.5 1999/11/22 13:58:40 mdw
59 * Add an option to disable Montgomery reduction, so that performance
60 * comparisons can be done.
61 *
62 * Revision 1.4 1999/11/21 12:27:06 mdw
63 * Remove a division from the Montgomery setup by calculating
64 * %$R^2 \bmod m$% first and then %$R \bmod m$% by Montgomery reduction of
65 * %$R^2$%.
66 *
67 * Revision 1.3 1999/11/21 11:35:10 mdw
68 * Performance improvement: use @mp_sqr@ and @mpmont_reduce@ instead of
69 * @mpmont_mul@ for squaring in exponentiation.
70 *
71 * Revision 1.2 1999/11/19 13:17:26 mdw
72 * Add extra interface to exponentiation which returns a Montgomerized
73 * result.
74 *
75 * Revision 1.1 1999/11/17 18:02:16 mdw
76 * New multiprecision integer arithmetic suite.
77 *
78 */
79
80 /*----- Header files ------------------------------------------------------*/
81
82 #include "mp.h"
83 #include "mpmont.h"
84
85 /*----- Tweakables --------------------------------------------------------*/
86
87 /* --- @MPMONT_DISABLE@ --- *
88 *
89 * Replace all the clever Montgomery reduction with good old-fashioned long
90 * division.
91 */
92
93 /* #define MPMONT_DISABLE */
94
95 /*----- Main code ---------------------------------------------------------*/
96
97 /* --- @mpmont_create@ --- *
98 *
99 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
100 * @mp *m@ = modulus to use
101 *
102 * Returns: ---
103 *
104 * Use: Initializes a Montgomery reduction context ready for use.
105 * The argument @m@ must be a positive odd integer.
106 */
107
108 #ifdef MPMONT_DISABLE
109
110 void mpmont_create(mpmont *mm, mp *m)
111 {
112 mp_shrink(m);
113 mm->m = MP_COPY(m);
114 mm->r = MP_ONE;
115 mm->r2 = MP_ONE;
116 mm->mi = MP_ONE;
117 }
118
119 #else
120
121 void mpmont_create(mpmont *mm, mp *m)
122 {
123 size_t n = MP_LEN(m);
124 mp *r2 = mp_new(2 * n + 1, 0);
125 mp r;
126
127 /* --- Validate the arguments --- */
128
129 assert(((void)"Montgomery modulus must be positive",
130 (m->f & MP_NEG) == 0));
131 assert(((void)"Montgomery modulus must be odd", m->v[0] & 1));
132
133 /* --- Take a copy of the modulus --- */
134
135 mp_shrink(m);
136 mm->m = MP_COPY(m);
137
138 /* --- Determine %$R^2$% --- */
139
140 mm->n = n;
141 MPX_ZERO(r2->v, r2->vl - 1);
142 r2->vl[-1] = 1;
143
144 /* --- Find the magic value @mi@ --- */
145
146 mp_build(&r, r2->v + n, r2->vl);
147 mm->mi = MP_NEW;
148 mp_gcd(0, 0, &mm->mi, &r, m);
149 mm->mi = mp_sub(mm->mi, &r, mm->mi);
150
151 /* --- Discover the values %$R \bmod m$% and %$R^2 \bmod m$% --- */
152
153 mm->r2 = MP_NEW;
154 mp_div(0, &mm->r2, r2, m);
155 mm->r = mpmont_reduce(mm, MP_NEW, mm->r2);
156 MP_DROP(r2);
157 }
158
159 #endif
160
161 /* --- @mpmont_destroy@ --- *
162 *
163 * Arguments: @mpmont *mm@ = pointer to a Montgomery reduction context
164 *
165 * Returns: ---
166 *
167 * Use: Disposes of a context when it's no longer of any use to
168 * anyone.
169 */
170
171 void mpmont_destroy(mpmont *mm)
172 {
173 MP_DROP(mm->m);
174 MP_DROP(mm->r);
175 MP_DROP(mm->r2);
176 MP_DROP(mm->mi);
177 }
178
179 /* --- @mpmont_reduce@ --- *
180 *
181 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
182 * @mp *d@ = destination
183 * @mp *a@ = source, assumed positive
184 *
185 * Returns: Result, %$a R^{-1} \bmod m$%.
186 */
187
188 #ifdef MPMONT_DISABLE
189
190 mp *mpmont_reduce(mpmont *mm, mp *d, mp *a)
191 {
192 mp_div(0, &d, a, mm->m);
193 return (d);
194 }
195
196 #else
197
198 mp *mpmont_reduce(mpmont *mm, mp *d, mp *a)
199 {
200 size_t n = mm->n;
201
202 /* --- Check for serious Karatsuba reduction --- */
203
204 if (n > MPK_THRESH * 3) {
205 mp al;
206 mpw *vl;
207 mp *u;
208
209 if (MP_LEN(a) >= n)
210 vl = a->v + n;
211 else
212 vl = a->vl;
213 mp_build(&al, a->v, vl);
214 u = mp_mul(MP_NEW, &al, mm->mi);
215 if (MP_LEN(u) > n)
216 u->vl = u->v + n;
217 u = mp_mul(u, u, mm->m);
218 d = mp_add(d, a, u);
219 mp_drop(u);
220 }
221
222 /* --- Otherwise do it the hard way --- */
223
224 else {
225 mpw *dv, *dvl;
226 mpw *mv, *mvl;
227 mpw mi;
228 size_t k = n;
229
230 /* --- Initial conditioning of the arguments --- */
231
232 a = MP_COPY(a);
233 if (d)
234 MP_DROP(d);
235 d = a;
236 MP_DEST(d, 2 * n + 1, a->f);
237
238 dv = d->v; dvl = d->vl;
239 mv = mm->m->v; mvl = mm->m->vl;
240
241 /* --- Let's go to work --- */
242
243 mi = mm->mi->v[0];
244 while (k--) {
245 mpw u = MPW(*dv * mi);
246 MPX_UMLAN(dv, dvl, mv, mvl, u);
247 dv++;
248 }
249 }
250
251 /* --- Wrap everything up --- */
252
253 memmove(d->v, d->v + n, MPWS(MP_LEN(d) - n));
254 d->vl -= n;
255 if (MPX_UCMP(d->v, d->vl, >=, mm->m->v, mm->m->vl))
256 mpx_usub(d->v, d->vl, d->v, d->vl, mm->m->v, mm->m->vl);
257 if (d->f & MP_NEG) {
258 mpx_usub(d->v, d->vl, mm->m->v, mm->m->vl, d->v, d->vl);
259 d->f &= ~MP_NEG;
260 }
261 MP_SHRINK(d);
262 return (d);
263 }
264
265 #endif
266
267 /* --- @mpmont_mul@ --- *
268 *
269 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
270 * @mp *d@ = destination
271 * @mp *a, *b@ = sources, assumed positive
272 *
273 * Returns: Result, %$a b R^{-1} \bmod m$%.
274 */
275
276 #ifdef MPMONT_DISABLE
277
278 mp *mpmont_mul(mpmont *mm, mp *d, mp *a, mp *b)
279 {
280 d = mp_mul(d, a, b);
281 mp_div(0, &d, d, mm->m);
282 return (d);
283 }
284
285 #else
286
287 mp *mpmont_mul(mpmont *mm, mp *d, mp *a, mp *b)
288 {
289 if (mm->n > MPK_THRESH * 3) {
290 d = mp_mul(d, a, b);
291 d = mpmont_reduce(mm, d, d);
292 } else {
293 mpw *dv, *dvl;
294 mpw *av, *avl;
295 mpw *bv, *bvl;
296 mpw *mv, *mvl;
297 mpw y;
298 size_t n, i;
299 mpw mi;
300
301 /* --- Initial conditioning of the arguments --- */
302
303 if (MP_LEN(a) > MP_LEN(b)) {
304 mp *t = a; a = b; b = t;
305 }
306 n = MP_LEN(mm->m);
307
308 a = MP_COPY(a);
309 b = MP_COPY(b);
310 MP_DEST(d, 2 * n + 1, a->f | b->f | MP_UNDEF);
311 dv = d->v; dvl = d->vl;
312 MPX_ZERO(dv, dvl);
313 av = a->v; avl = a->vl;
314 bv = b->v; bvl = b->vl;
315 mv = mm->m->v; mvl = mm->m->vl;
316 y = *bv;
317
318 /* --- Montgomery multiplication phase --- */
319
320 i = 0;
321 mi = mm->mi->v[0];
322 while (i < n && av < avl) {
323 mpw x = *av++;
324 mpw u = MPW((*dv + x * y) * mi);
325 MPX_UMLAN(dv, dvl, bv, bvl, x);
326 MPX_UMLAN(dv, dvl, mv, mvl, u);
327 dv++;
328 i++;
329 }
330
331 /* --- Simpler Montgomery reduction phase --- */
332
333 while (i < n) {
334 mpw u = MPW(*dv * mi);
335 MPX_UMLAN(dv, dvl, mv, mvl, u);
336 dv++;
337 i++;
338 }
339
340 /* --- Done --- */
341
342 memmove(d->v, dv, MPWS(dvl - dv));
343 d->vl -= dv - d->v;
344 if (MPX_UCMP(d->v, d->vl, >=, mm->m->v, mm->m->vl))
345 mpx_usub(d->v, d->vl, d->v, d->vl, mm->m->v, mm->m->vl);
346 if ((a->f ^ b->f) & MP_NEG)
347 mpx_usub(d->v, d->vl, mm->m->v, mm->m->vl, d->v, d->vl);
348 MP_SHRINK(d);
349 d->f = (a->f | b->f) & MP_BURN;
350 MP_DROP(a);
351 MP_DROP(b);
352 }
353
354 return (d);
355 }
356
357 #endif
358
359 /* --- @mpmont_expr@ --- *
360 *
361 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
362 * @mp *d@ = fake destination
363 * @mp *a@ = base
364 * @mp *e@ = exponent
365 *
366 * Returns: Result, %$a^e R \bmod m$%.
367 */
368
369 #define WINSZ 5
370 #define TABSZ (1 << (WINSZ - 1))
371
372 #define THRESH (((MPW_BITS / WINSZ) << 2) + 1)
373
374 static mp *exp_simple(mpmont *mm, mp *d, mp *a, mp *e)
375 {
376 mpscan sc;
377 mp *ar;
378 mp *x = MP_COPY(mm->r);
379 mp *spare = (e->f & MP_BURN) ? MP_NEWSEC : MP_NEW;
380 unsigned sq = 0;
381
382 mp_rscan(&sc, e);
383 if (!MP_RSTEP(&sc))
384 goto exit;
385 while (!MP_RBIT(&sc))
386 MP_RSTEP(&sc);
387
388 /* --- Do the main body of the work --- */
389
390 ar = mpmont_mul(mm, MP_NEW, a, mm->r2);
391 for (;;) {
392 sq++;
393 while (sq) {
394 mp *y;
395 y = mp_sqr(spare, x);
396 y = mpmont_reduce(mm, y, y);
397 spare = x; x = y;
398 sq--;
399 }
400 { mp *y = mpmont_mul(mm, spare, x, ar); spare = x; x = y; }
401 sq = 0;
402 for (;;) {
403 if (!MP_RSTEP(&sc))
404 goto done;
405 if (MP_RBIT(&sc))
406 break;
407 sq++;
408 }
409 }
410
411 /* --- Do a final round of squaring --- */
412
413 done:
414 while (sq) {
415 mp *y;
416 y = mp_sqr(spare, x);
417 y = mpmont_reduce(mm, y, y);
418 spare = x; x = y;
419 sq--;
420 }
421
422 /* --- Done --- */
423
424 MP_DROP(ar);
425 exit:
426 if (spare != MP_NEW)
427 MP_DROP(spare);
428 if (d != MP_NEW)
429 MP_DROP(d);
430 return (x);
431 }
432
433 mp *mpmont_expr(mpmont *mm, mp *d, mp *a, mp *e)
434 {
435 mp **tab;
436 mp *ar, *a2;
437 mp *spare = (e->f & MP_BURN) ? MP_NEWSEC : MP_NEW;
438 mp *x = MP_COPY(mm->r);
439 unsigned i, sq = 0;
440 mpscan sc;
441
442 /* --- Do we bother? --- */
443
444 MP_SHRINK(e);
445 if (MP_LEN(e) == 0)
446 goto exit;
447 if (MP_LEN(e) < THRESH) {
448 x->ref--;
449 return (exp_simple(mm, d, a, e));
450 }
451
452 /* --- Do the precomputation --- */
453
454 ar = mpmont_mul(mm, MP_NEW, a, mm->r2);
455 a2 = mp_sqr(MP_NEW, ar);
456 a2 = mpmont_reduce(mm, a2, a2);
457 tab = xmalloc(TABSZ * sizeof(mp *));
458 tab[0] = ar;
459 for (i = 1; i < TABSZ; i++)
460 tab[i] = mpmont_mul(mm, MP_NEW, tab[i - 1], a2);
461 mp_drop(a2);
462 mp_rscan(&sc, e);
463
464 /* --- Skip top-end zero bits --- *
465 *
466 * If the initial step worked, there must be a set bit somewhere, so keep
467 * stepping until I find it.
468 */
469
470 MP_RSTEP(&sc);
471 while (!MP_RBIT(&sc)) {
472 MP_RSTEP(&sc);
473 }
474
475 /* --- Now for the main work --- */
476
477 for (;;) {
478 unsigned l = 0;
479 unsigned z = 0;
480
481 /* --- The next bit is set, so read a window index --- *
482 *
483 * Reset @i@ to zero and increment @sq@. Then, until either I read
484 * @WINSZ@ bits or I run out of bits, scan in a bit: if it's clear, bump
485 * the @z@ counter; if it's set, push a set bit into @i@, shift it over
486 * by @z@ bits, bump @sq@ by @z + 1@ and clear @z@. By the end of this
487 * palaver, @i@ is an index to the precomputed value in @tab@.
488 */
489
490 i = 0;
491 sq++;
492 for (;;) {
493 l++;
494 if (l >= WINSZ || !MP_RSTEP(&sc))
495 break;
496 if (!MP_RBIT(&sc))
497 z++;
498 else {
499 i = ((i << 1) | 1) << z;
500 sq += z + 1;
501 z = 0;
502 }
503 }
504
505 /* --- Do the squaring --- *
506 *
507 * Remember that @sq@ carries over from the zero-skipping stuff below.
508 */
509
510 while (sq) {
511 mp *y;
512 y = mp_sqr(spare, x);
513 y = mpmont_reduce(mm, y, y);
514 spare = x; x = y;
515 sq--;
516 }
517
518 /* --- Do the multiply --- */
519
520 { mp *y = mpmont_mul(mm, spare, x, tab[i]); spare = x; x = y; }
521
522 /* --- Now grind along through the rest of the bits --- */
523
524 sq = z;
525 for (;;) {
526 if (!MP_RSTEP(&sc))
527 goto done;
528 if (MP_RBIT(&sc))
529 break;
530 sq++;
531 }
532 }
533
534 /* --- Do a final round of squaring --- */
535
536 done:
537 while (sq) {
538 mp *y;
539 y = mp_sqr(spare, x);
540 y = mpmont_reduce(mm, y, y);
541 spare = x; x = y;
542 sq--;
543 }
544
545 /* --- Done --- */
546
547 for (i = 0; i < TABSZ; i++)
548 mp_drop(tab[i]);
549 xfree(tab);
550 exit:
551 if (d != MP_NEW)
552 mp_drop(d);
553 if (spare)
554 mp_drop(spare);
555 return (x);
556 }
557
558 /* --- @mpmont_exp@ --- *
559 *
560 * Arguments: @mpmont *mm@ = pointer to Montgomery reduction context
561 * @mp *d@ = fake destination
562 * @mp *a@ = base
563 * @mp *e@ = exponent
564 *
565 * Returns: Result, %$a^e \bmod m$%.
566 */
567
568 mp *mpmont_exp(mpmont *mm, mp *d, mp *a, mp *e)
569 {
570 d = mpmont_expr(mm, d, a, e);
571 d = mpmont_reduce(mm, d, d);
572 return (d);
573 }
574
575 /*----- Test rig ----------------------------------------------------------*/
576
577 #ifdef TEST_RIG
578
579 static int tcreate(dstr *v)
580 {
581 mp *m = *(mp **)v[0].buf;
582 mp *mi = *(mp **)v[1].buf;
583 mp *r = *(mp **)v[2].buf;
584 mp *r2 = *(mp **)v[3].buf;
585
586 mpmont mm;
587 int ok = 1;
588
589 mpmont_create(&mm, m);
590
591 if (mm.mi->v[0] != mi->v[0]) {
592 fprintf(stderr, "\n*** bad mi: found %lu, expected %lu",
593 (unsigned long)mm.mi->v[0], (unsigned long)mi->v[0]);
594 fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
595 fputc('\n', stderr);
596 ok = 0;
597 }
598
599 if (!MP_EQ(mm.r, r)) {
600 fputs("\n*** bad r", stderr);
601 fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
602 fputs("\nexpected ", stderr); mp_writefile(r, stderr, 10);
603 fputs("\n found ", stderr); mp_writefile(mm.r, stderr, 10);
604 fputc('\n', stderr);
605 ok = 0;
606 }
607
608 if (!MP_EQ(mm.r2, r2)) {
609 fputs("\n*** bad r2", stderr);
610 fputs("\nm = ", stderr); mp_writefile(m, stderr, 10);
611 fputs("\nexpected ", stderr); mp_writefile(r2, stderr, 10);
612 fputs("\n found ", stderr); mp_writefile(mm.r2, stderr, 10);
613 fputc('\n', stderr);
614 ok = 0;
615 }
616
617 MP_DROP(m);
618 MP_DROP(mi);
619 MP_DROP(r);
620 MP_DROP(r2);
621 mpmont_destroy(&mm);
622 assert(mparena_count(MPARENA_GLOBAL) == 0);
623 return (ok);
624 }
625
626 static int tmul(dstr *v)
627 {
628 mp *m = *(mp **)v[0].buf;
629 mp *a = *(mp **)v[1].buf;
630 mp *b = *(mp **)v[2].buf;
631 mp *r = *(mp **)v[3].buf;
632 int ok = 1;
633
634 mpmont mm;
635 mpmont_create(&mm, m);
636
637 {
638 mp *qr = mp_mul(MP_NEW, a, b);
639 mp_div(0, &qr, qr, m);
640
641 if (!MP_EQ(qr, r)) {
642 fputs("\n*** classical modmul failed", stderr);
643 fputs("\n m = ", stderr); mp_writefile(m, stderr, 10);
644 fputs("\n a = ", stderr); mp_writefile(a, stderr, 10);
645 fputs("\n b = ", stderr); mp_writefile(b, stderr, 10);
646 fputs("\n r = ", stderr); mp_writefile(r, stderr, 10);
647 fputs("\nqr = ", stderr); mp_writefile(qr, stderr, 10);
648 fputc('\n', stderr);
649 ok = 0;
650 }
651
652 mp_drop(qr);
653 }
654
655 {
656 mp *ar = mpmont_mul(&mm, MP_NEW, a, mm.r2);
657 mp *br = mpmont_mul(&mm, MP_NEW, b, mm.r2);
658 mp *mr = mpmont_mul(&mm, MP_NEW, ar, br);
659 mr = mpmont_reduce(&mm, mr, mr);
660 if (!MP_EQ(mr, r)) {
661 fputs("\n*** montgomery modmul failed", stderr);
662 fputs("\n m = ", stderr); mp_writefile(m, stderr, 10);
663 fputs("\n a = ", stderr); mp_writefile(a, stderr, 10);
664 fputs("\n b = ", stderr); mp_writefile(b, stderr, 10);
665 fputs("\n r = ", stderr); mp_writefile(r, stderr, 10);
666 fputs("\nmr = ", stderr); mp_writefile(mr, stderr, 10);
667 fputc('\n', stderr);
668 ok = 0;
669 }
670 MP_DROP(ar); MP_DROP(br);
671 mp_drop(mr);
672 }
673
674
675 MP_DROP(m);
676 MP_DROP(a);
677 MP_DROP(b);
678 MP_DROP(r);
679 mpmont_destroy(&mm);
680 assert(mparena_count(MPARENA_GLOBAL) == 0);
681 return ok;
682 }
683
684 static int texp(dstr *v)
685 {
686 mp *m = *(mp **)v[0].buf;
687 mp *a = *(mp **)v[1].buf;
688 mp *b = *(mp **)v[2].buf;
689 mp *r = *(mp **)v[3].buf;
690 mp *mr;
691 int ok = 1;
692
693 mpmont mm;
694 mpmont_create(&mm, m);
695
696 mr = mpmont_exp(&mm, MP_NEW, a, b);
697
698 if (!MP_EQ(mr, r)) {
699 fputs("\n*** montgomery modexp failed", stderr);
700 fputs("\n m = ", stderr); mp_writefile(m, stderr, 10);
701 fputs("\n a = ", stderr); mp_writefile(a, stderr, 10);
702 fputs("\n e = ", stderr); mp_writefile(b, stderr, 10);
703 fputs("\n r = ", stderr); mp_writefile(r, stderr, 10);
704 fputs("\nmr = ", stderr); mp_writefile(mr, stderr, 10);
705 fputc('\n', stderr);
706 ok = 0;
707 }
708
709 MP_DROP(m);
710 MP_DROP(a);
711 MP_DROP(b);
712 MP_DROP(r);
713 MP_DROP(mr);
714 mpmont_destroy(&mm);
715 assert(mparena_count(MPARENA_GLOBAL) == 0);
716 return ok;
717 }
718
719
720 static test_chunk tests[] = {
721 { "create", tcreate, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
722 { "mul", tmul, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
723 { "exp", texp, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
724 { 0, 0, { 0 } },
725 };
726
727 int main(int argc, char *argv[])
728 {
729 sub_init();
730 test_run(argc, argv, tests, SRCDIR "/tests/mpmont");
731 return (0);
732 }
733
734 #endif
735
736 /*----- That's all, folks -------------------------------------------------*/