X-Git-Url: https://git.distorted.org.uk/u/mdw/catacomb/blobdiff_plain/4468424ee46294bb2073ca546d54cced2b9bdf02..2e8eb64a72472bea4faba1cee5edde4fe3616808:/mpx-kmul.c diff --git a/mpx-kmul.c b/mpx-kmul.c index d53622b..f6b0b4f 100644 --- a/mpx-kmul.c +++ b/mpx-kmul.c @@ -1,6 +1,6 @@ /* -*-c-*- * - * $Id: mpx-kmul.c,v 1.3 1999/12/13 15:35:01 mdw Exp $ + * $Id: mpx-kmul.c,v 1.7 2000/10/08 15:48:35 mdw Exp $ * * Karatsuba's multiplication algorithm * @@ -30,6 +30,19 @@ /*----- Revision history --------------------------------------------------* * * $Log: mpx-kmul.c,v $ + * Revision 1.7 2000/10/08 15:48:35 mdw + * Rename Karatsuba constants now that we have @gfx_kmul@ too. + * + * Revision 1.6 2000/10/08 12:11:01 mdw + * Use @mpx_ueq@ instead of @MPX_UCMP@. + * + * Revision 1.5 2000/07/29 17:04:02 mdw + * Remove useless header `mpscan.h'. + * + * Revision 1.4 2000/06/17 11:42:11 mdw + * Moved the Karatsuba macros into a separate file for better sharing. + * Fixed some comments. + * * Revision 1.3 1999/12/13 15:35:01 mdw * Simplify and improve. * @@ -47,84 +60,15 @@ #include #include "mpx.h" +#include "karatsuba.h" /*----- Tweakables --------------------------------------------------------*/ #ifdef TEST_RIG -# undef KARATSUBA_CUTOFF -# define KARATSUBA_CUTOFF 2 +# undef MPK_THRESH +# define MPK_THRESH 2 #endif -/*----- Addition macros ---------------------------------------------------*/ - -#define UADD(dv, av, avl) do { \ - mpw *_dv = (dv); \ - const mpw *_av = (av), *_avl = (avl); \ - mpw _c = 0; \ - \ - while (_av < _avl) { \ - mpw _a, _b; \ - mpd _x; \ - _a = *_av++; \ - _b = *_dv; \ - _x = (mpd)_a + (mpd)_b + _c; \ - *_dv++ = MPW(_x); \ - _c = _x >> MPW_BITS; \ - } \ - while (_c) { \ - mpd _x = (mpd)*_dv + (mpd)_c; \ - *_dv++ = MPW(_x); \ - _c = _x >> MPW_BITS; \ - } \ -} while (0) - -#define UADD2(dv, dvl, av, avl, bv, bvl) do { \ - mpw *_dv = (dv), *_dvl = (dvl); \ - const mpw *_av = (av), *_avl = (avl); \ - const mpw *_bv = (bv), *_bvl = (bvl); \ - mpw _c = 0; \ - \ - while (_av < _avl || _bv < _bvl) { \ - mpw _a, _b; \ - mpd _x; \ - _a = (_av < _avl) ? *_av++ : 0; \ - _b = (_bv < _bvl) ? *_bv++ : 0; \ - _x = (mpd)_a + (mpd)_b + _c; \ - *_dv++ = MPW(_x); \ - _c = _x >> MPW_BITS; \ - } \ - *_dv++ = _c; \ - while (_dv < _dvl) \ - *_dv++ = 0; \ -} while (0) - -#define USUB(dv, av, avl) do { \ - mpw *_dv = (dv); \ - const mpw *_av = (av), *_avl = (avl); \ - mpw _c = 0; \ - \ - while (_av < _avl) { \ - mpw _a, _b; \ - mpd _x; \ - _a = *_av++; \ - _b = *_dv; \ - _x = (mpd)_b - (mpd)_a - _c; \ - *_dv++ = MPW(_x); \ - if (_x >> MPW_BITS) \ - _c = 1; \ - else \ - _c = 0; \ - } \ - while (_c) { \ - mpd _x = (mpd)*_dv - (mpd)_c; \ - *_dv++ = MPW(_x); \ - if (_x >> MPW_BITS) \ - _c = 1; \ - else \ - _c = 0; \ - } \ -} while (0) - /*----- Main code ---------------------------------------------------------*/ /* --- @mpx_kmul@ --- * @@ -143,7 +87,7 @@ * * The destination must be twice as large as the larger * argument. The scratch space must be twice as large as the - * larger argument, plus the magic number @KARATSUBA_SLOP@. + * larger argument, plus the magic number @MPK_SLOP@. */ void mpx_kmul(mpw *dv, mpw *dvl, @@ -165,18 +109,18 @@ void mpx_kmul(mpw *dv, mpw *dvl, MPX_SHRINK(av, avl); MPX_SHRINK(bv, bvl); - if (avl - av <= KARATSUBA_CUTOFF || bvl - bv <= KARATSUBA_CUTOFF) { + if (avl - av <= MPK_THRESH || bvl - bv <= MPK_THRESH) { mpx_umul(dv, dvl, av, avl, bv, bvl); return; } /* --- How the algorithm works --- * * - * Let %$A = xb + y$% and %$B = ub + v$%. Then, simply by expanding, %$AB - * = x u b^2 + b(x v + y u) + y v$%. That's not helped any, because I've - * got four multiplications, each four times easier than the one I started - * with. However, note that I can rewrite the coefficient of %$b$% as - * %$xv + yu = (x + y)(u + v) - xu - yv$%. The terms %$xu$% and %$yv$% + * Let %$A = xb + y$% and %$B = ub + v$%. Then, simply by expanding, + * %$AB = x u b^2 + b(x v + y u) + y v$%. That's not helped any, because + * I've got four multiplications, each four times easier than the one I + * started with. However, note that I can rewrite the coefficient of %$b$% + * as %$xv + yu = (x + y)(u + v) - xu - yv$%. The terms %$xu$% and %$yv$% * I've already calculated, and that leaves only one more multiplication to * do. So now I have three multiplications, each four times easier, and * that's a win. @@ -185,7 +129,7 @@ void mpx_kmul(mpw *dv, mpw *dvl, /* --- First things --- * * * Sort out where to break the factors in half. I'll choose the midpoint - * of the largest one, since this minimizes the amount of work I have to do + * of the larger one, since this minimizes the amount of work I have to do * most effectively. */ @@ -218,7 +162,7 @@ void mpx_kmul(mpw *dv, mpw *dvl, UADD2(sv, bsv, av, avm, avm, avl); UADD2(bsv, ssv, bv, bvm, bvm, bvl); - if (m > KARATSUBA_CUTOFF) + if (m > MPK_THRESH) mpx_kmul(rdv, rdvl, sv, bsv, bsv, ssv, ssv, svl); else mpx_umul(rdv, rdvl, sv, bsv, bsv, ssv); @@ -234,7 +178,7 @@ void mpx_kmul(mpw *dv, mpw *dvl, if (avl == avm || bvl == bvm) MPX_ZERO(rdv + m + 1, dvl); else { - if (m > KARATSUBA_CUTOFF) + if (m > MPK_THRESH) mpx_kmul(sv, ssv, avm, avl, bvm, bvl, ssv, svl); else mpx_umul(sv, ssv, avm, avl, bvm, bvl); @@ -243,7 +187,7 @@ void mpx_kmul(mpw *dv, mpw *dvl, USUB(tdv, sv, svn); } - if (m > KARATSUBA_CUTOFF) + if (m > MPK_THRESH) mpx_kmul(sv, ssv, av, avm, bv, bvm, ssv, svl); else mpx_umul(sv, ssv, av, avm, bv, bvm); @@ -260,8 +204,6 @@ void mpx_kmul(mpw *dv, mpw *dvl, #include #include -#include "mpscan.h" - #define ALLOC(v, vl, sz) do { \ size_t _sz = (sz); \ mpw *_vv = xmalloc(MPWS(_sz)); \ @@ -308,7 +250,7 @@ static int umul(dstr *v) ALLOC(s, sl, 2 * m + 32); mpx_kmul(d, dl, a, al, b, bl, s, sl); - if (MPX_UCMP(d, dl, !=, c, cl)) { + if (!mpx_ueq(d, dl, c, cl)) { fprintf(stderr, "\n*** umul failed\n"); dumpmp(" a", a, al); dumpmp(" b", b, bl);