Annotation of OpenXM_contrib/gmp/mpn/generic/mul_fft.c, Revision 1.1
1.1 ! maekawa 1: /* An implementation in GMP of Scho"nhage's fast multiplication algorithm
! 2: modulo 2^N+1, by Paul Zimmermann, INRIA Lorraine, February 1998.
! 3:
! 4: THE CONTENTS OF THIS FILE ARE FOR INTERNAL USE AND THE FUNCTIONS HAVE
! 5: MUTABLE INTERFACES. IT IS ONLY SAFE TO REACH THEM THROUGH DOCUMENTED
! 6: INTERFACES. IT IS ALMOST GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN
! 7: A FUTURE GNU MP RELEASE.
! 8:
! 9: Copyright (C) 1998, 1999, 2000 Free Software Foundation, Inc.
! 10:
! 11: This file is part of the GNU MP Library.
! 12:
! 13: The GNU MP Library is free software; you can redistribute it and/or modify
! 14: it under the terms of the GNU Lesser General Public License as published by
! 15: the Free Software Foundation; either version 2.1 of the License, or (at your
! 16: option) any later version.
! 17:
! 18: The GNU MP Library is distributed in the hope that it will be useful, but
! 19: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
! 20: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
! 21: License for more details.
! 22:
! 23: You should have received a copy of the GNU Lesser General Public License
! 24: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
! 25: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
! 26: MA 02111-1307, USA. */
! 27:
! 28:
! 29: /* References:
! 30:
! 31: Schnelle Multiplikation grosser Zahlen, by Arnold Scho"nhage and Volker
! 32: Strassen, Computing 7, p. 281-292, 1971.
! 33:
! 34: Asymptotically fast algorithms for the numerical multiplication
! 35: and division of polynomials with complex coefficients, by Arnold Scho"nhage,
! 36: Computer Algebra, EUROCAM'82, LNCS 144, p. 3-15, 1982.
! 37:
! 38: Tapes versus Pointers, a study in implementing fast algorithms,
! 39: by Arnold Scho"nhage, Bulletin of the EATCS, 30, p. 23-32, 1986.
! 40:
! 41: See also http://www.loria.fr/~zimmerma/bignum
! 42:
! 43:
! 44: Future:
! 45:
! 46: K==2 isn't needed in the current uses of this code and the bits specific
! 47: for that could be dropped.
! 48:
! 49: It might be possible to avoid a small number of MPN_COPYs by using a
! 50: rotating temporary or two.
! 51:
! 52: Multiplications of unequal sized operands can be done with this code, but
! 53: it needs a tighter test for identifying squaring (same sizes as well as
! 54: same pointers). */
! 55:
! 56:
! 57: #include <stdio.h>
! 58: #include "gmp.h"
! 59: #include "gmp-impl.h"
! 60:
! 61:
! 62: /* Change this to "#define TRACE(x) x" for some traces. */
! 63: #define TRACE(x)
! 64:
! 65:
! 66:
! 67: FFT_TABLE_ATTRS mp_size_t mpn_fft_table[2][MPN_FFT_TABLE_SIZE] = {
! 68: FFT_MUL_TABLE,
! 69: FFT_SQR_TABLE
! 70: };
! 71:
! 72:
! 73: static void mpn_mul_fft_internal
! 74: _PROTO ((mp_limb_t *op, mp_srcptr n, mp_srcptr m, mp_size_t pl,
! 75: int k, int K,
! 76: mp_limb_t **Ap, mp_limb_t **Bp,
! 77: mp_limb_t *A, mp_limb_t *B,
! 78: mp_size_t nprime, mp_size_t l, mp_size_t Mp, int **_fft_l,
! 79: mp_limb_t *T, int rec));
! 80:
! 81:
! 82: /* Find the best k to use for a mod 2^(n*BITS_PER_MP_LIMB)+1 FFT.
! 83: sqr==0 if for a multiply, sqr==1 for a square */
! 84: int
! 85: #if __STDC__
! 86: mpn_fft_best_k (mp_size_t n, int sqr)
! 87: #else
! 88: mpn_fft_best_k (n, sqr)
! 89: mp_size_t n;
! 90: int sqr;
! 91: #endif
! 92: {
! 93: mp_size_t t;
! 94: int i;
! 95:
! 96: for (i = 0; mpn_fft_table[sqr][i] != 0; i++)
! 97: if (n < mpn_fft_table[sqr][i])
! 98: return i + FFT_FIRST_K;
! 99:
! 100: /* treat 4*last as one further entry */
! 101: if (i == 0 || n < 4*mpn_fft_table[sqr][i-1])
! 102: return i + FFT_FIRST_K;
! 103: else
! 104: return i + FFT_FIRST_K + 1;
! 105: }
! 106:
! 107:
! 108: /* Returns smallest possible number of limbs >= pl for a fft of size 2^k.
! 109: FIXME: Is this simply pl rounded up to the next multiple of 2^k ? */
! 110:
! 111: mp_size_t
! 112: #if __STDC__
! 113: mpn_fft_next_size (mp_size_t pl, int k)
! 114: #else
! 115: mpn_fft_next_size (pl, k)
! 116: mp_size_t pl;
! 117: int k;
! 118: #endif
! 119: {
! 120: mp_size_t N, M;
! 121: int K;
! 122:
! 123: /* if (k==0) k = mpn_fft_best_k (pl, sqr); */
! 124: N = pl*BITS_PER_MP_LIMB;
! 125: K = 1<<k;
! 126: if (N%K) N=(N/K+1)*K;
! 127: M = N/K;
! 128: if (M%BITS_PER_MP_LIMB) N=((M/BITS_PER_MP_LIMB)+1)*BITS_PER_MP_LIMB*K;
! 129: return (N/BITS_PER_MP_LIMB);
! 130: }
! 131:
! 132:
! 133: static void
! 134: #if __STDC__
! 135: mpn_fft_initl(int **l, int k)
! 136: #else
! 137: mpn_fft_initl(l, k)
! 138: int **l;
! 139: int k;
! 140: #endif
! 141: {
! 142: int i,j,K;
! 143:
! 144: l[0][0] = 0;
! 145: for (i=1,K=2;i<=k;i++,K*=2) {
! 146: for (j=0;j<K/2;j++) {
! 147: l[i][j] = 2*l[i-1][j];
! 148: l[i][K/2+j] = 1+l[i][j];
! 149: }
! 150: }
! 151: }
! 152:
! 153:
! 154: /* a <- -a mod 2^(n*BITS_PER_MP_LIMB)+1 */
! 155: static void
! 156: #if __STDC__
! 157: mpn_fft_neg_modF(mp_limb_t *ap, mp_size_t n)
! 158: #else
! 159: mpn_fft_neg_modF(ap, n)
! 160: mp_limb_t *ap;
! 161: mp_size_t n;
! 162: #endif
! 163: {
! 164: mp_limb_t c;
! 165:
! 166: c = ap[n]+2;
! 167: mpn_com_n (ap, ap, n);
! 168: ap[n]=0; mpn_incr_u(ap, c);
! 169: }
! 170:
! 171:
! 172: /* a <- a*2^e mod 2^(n*BITS_PER_MP_LIMB)+1 */
! 173: static void
! 174: #if __STDC__
! 175: mpn_fft_mul_2exp_modF(mp_limb_t *ap, int e, mp_size_t n, mp_limb_t *tp)
! 176: #else
! 177: mpn_fft_mul_2exp_modF(ap, e, n, tp)
! 178: mp_limb_t *ap;
! 179: int e;
! 180: mp_size_t n;
! 181: mp_limb_t *tp;
! 182: #endif
! 183: {
! 184: int d, sh, i; mp_limb_t cc;
! 185:
! 186: d = e%(n*BITS_PER_MP_LIMB); /* 2^e = (+/-) 2^d */
! 187: sh = d % BITS_PER_MP_LIMB;
! 188: if (sh) mpn_lshift(tp, ap, n+1, sh); /* no carry here */
! 189: else MPN_COPY(tp, ap, n+1);
! 190: d /= BITS_PER_MP_LIMB; /* now shift of d limbs to the left */
! 191: if (d) {
! 192: /* ap[d..n-1] = tp[0..n-d-1], ap[0..d-1] = -tp[n-d..n-1] */
! 193: /* mpn_xor would be more efficient here */
! 194: for (i=d-1;i>=0;i--) ap[i] = ~tp[n-d+i];
! 195: cc = 1-mpn_add_1(ap, ap, d, 1);
! 196: if (cc) cc=mpn_sub_1(ap+d, tp, n-d, 1);
! 197: else MPN_COPY(ap+d, tp, n-d);
! 198: if (cc+=mpn_sub_1(ap+d, ap+d, n-d, tp[n]))
! 199: ap[n]=mpn_add_1(ap, ap, n, cc);
! 200: else ap[n]=0;
! 201: }
! 202: else if ((ap[n]=mpn_sub_1(ap, tp, n, tp[n]))) {
! 203: ap[n]=mpn_add_1(ap, ap, n, 1);
! 204: }
! 205: if ((e/(n*BITS_PER_MP_LIMB))%2) mpn_fft_neg_modF(ap, n);
! 206: }
! 207:
! 208:
! 209: /* a <- a+b mod 2^(n*BITS_PER_MP_LIMB)+1 */
! 210: static void
! 211: #if __STDC__
! 212: mpn_fft_add_modF (mp_limb_t *ap, mp_limb_t *bp, int n)
! 213: #else
! 214: mpn_fft_add_modF (ap, bp, n)
! 215: mp_limb_t *ap,*bp;
! 216: int n;
! 217: #endif
! 218: {
! 219: mp_limb_t c;
! 220:
! 221: c = ap[n] + bp[n] + mpn_add_n(ap, ap, bp, n);
! 222: if (c>1) c -= 1+mpn_sub_1(ap,ap,n,1);
! 223: ap[n]=c;
! 224: }
! 225:
! 226:
! 227: /* input: A[0] ... A[inc*(K-1)] are residues mod 2^N+1 where
! 228: N=n*BITS_PER_MP_LIMB
! 229: 2^omega is a primitive root mod 2^N+1
! 230: output: A[inc*l[k][i]] <- \sum (2^omega)^(ij) A[inc*j] mod 2^N+1 */
! 231:
! 232: static void
! 233: #if __STDC__
! 234: mpn_fft_fft_sqr (mp_limb_t **Ap, mp_size_t K, int **ll,
! 235: mp_size_t omega, mp_size_t n, mp_size_t inc, mp_limb_t *tp)
! 236: #else
! 237: mpn_fft_fft_sqr(Ap,K,ll,omega,n,inc,tp)
! 238: mp_limb_t **Ap,*tp;
! 239: mp_size_t K,omega,n,inc;
! 240: int **ll;
! 241: #endif
! 242: {
! 243: if (K==2) {
! 244: #ifdef ADDSUB
! 245: if (mpn_addsub_n(Ap[0], Ap[inc], Ap[0], Ap[inc], n+1) & 1)
! 246: #else
! 247: MPN_COPY(tp, Ap[0], n+1);
! 248: mpn_add_n(Ap[0], Ap[0], Ap[inc],n+1);
! 249: if (mpn_sub_n(Ap[inc], tp, Ap[inc],n+1))
! 250: #endif
! 251: Ap[inc][n] = mpn_add_1(Ap[inc], Ap[inc], n, 1);
! 252: }
! 253: else {
! 254: int j, inc2=2*inc;
! 255: int *lk = *ll;
! 256: mp_limb_t *tmp;
! 257: TMP_DECL(marker);
! 258:
! 259: TMP_MARK(marker);
! 260: tmp = TMP_ALLOC_LIMBS (n+1);
! 261: mpn_fft_fft_sqr(Ap, K/2,ll-1,2*omega,n,inc2, tp);
! 262: mpn_fft_fft_sqr(Ap+inc, K/2,ll-1,2*omega,n,inc2, tp);
! 263: /* A[2*j*inc] <- A[2*j*inc] + omega^l[k][2*j*inc] A[(2j+1)inc]
! 264: A[(2j+1)inc] <- A[2*j*inc] + omega^l[k][(2j+1)inc] A[(2j+1)inc] */
! 265: for (j=0;j<K/2;j++,lk+=2,Ap+=2*inc) {
! 266: MPN_COPY(tp, Ap[inc], n+1);
! 267: mpn_fft_mul_2exp_modF(Ap[inc], lk[1]*omega, n, tmp);
! 268: mpn_fft_add_modF(Ap[inc], Ap[0], n);
! 269: mpn_fft_mul_2exp_modF(tp,lk[0]*omega, n, tmp);
! 270: mpn_fft_add_modF(Ap[0], tp, n);
! 271: }
! 272: TMP_FREE(marker);
! 273: }
! 274: }
! 275:
! 276:
! 277: /* input: A[0] ... A[inc*(K-1)] are residues mod 2^N+1 where
! 278: N=n*BITS_PER_MP_LIMB
! 279: 2^omega is a primitive root mod 2^N+1
! 280: output: A[inc*l[k][i]] <- \sum (2^omega)^(ij) A[inc*j] mod 2^N+1 */
! 281:
! 282: static void
! 283: #if __STDC__
! 284: mpn_fft_fft (mp_limb_t **Ap, mp_limb_t **Bp, mp_size_t K, int **ll,
! 285: mp_size_t omega, mp_size_t n, mp_size_t inc, mp_limb_t *tp)
! 286: #else
! 287: mpn_fft_fft(Ap,Bp,K,ll,omega,n,inc,tp)
! 288: mp_limb_t **Ap,**Bp,*tp;
! 289: mp_size_t K,omega,n,inc;
! 290: int **ll;
! 291: #endif
! 292: {
! 293: if (K==2) {
! 294: #ifdef ADDSUB
! 295: if (mpn_addsub_n(Ap[0], Ap[inc], Ap[0], Ap[inc], n+1) & 1)
! 296: #else
! 297: MPN_COPY(tp, Ap[0], n+1);
! 298: mpn_add_n(Ap[0], Ap[0], Ap[inc],n+1);
! 299: if (mpn_sub_n(Ap[inc], tp, Ap[inc],n+1))
! 300: #endif
! 301: Ap[inc][n] = mpn_add_1(Ap[inc], Ap[inc], n, 1);
! 302: #ifdef ADDSUB
! 303: if (mpn_addsub_n(Bp[0], Bp[inc], Bp[0], Bp[inc], n+1) & 1)
! 304: #else
! 305: MPN_COPY(tp, Bp[0], n+1);
! 306: mpn_add_n(Bp[0], Bp[0], Bp[inc],n+1);
! 307: if (mpn_sub_n(Bp[inc], tp, Bp[inc],n+1))
! 308: #endif
! 309: Bp[inc][n] = mpn_add_1(Bp[inc], Bp[inc], n, 1);
! 310: }
! 311: else {
! 312: int j, inc2=2*inc;
! 313: int *lk=*ll;
! 314: mp_limb_t *tmp;
! 315: TMP_DECL(marker);
! 316:
! 317: TMP_MARK(marker);
! 318: tmp = TMP_ALLOC_LIMBS (n+1);
! 319: mpn_fft_fft(Ap, Bp, K/2,ll-1,2*omega,n,inc2, tp);
! 320: mpn_fft_fft(Ap+inc, Bp+inc, K/2,ll-1,2*omega,n,inc2, tp);
! 321: /* A[2*j*inc] <- A[2*j*inc] + omega^l[k][2*j*inc] A[(2j+1)inc]
! 322: A[(2j+1)inc] <- A[2*j*inc] + omega^l[k][(2j+1)inc] A[(2j+1)inc] */
! 323: for (j=0;j<K/2;j++,lk+=2,Ap+=2*inc,Bp+=2*inc) {
! 324: MPN_COPY(tp, Ap[inc], n+1);
! 325: mpn_fft_mul_2exp_modF(Ap[inc], lk[1]*omega, n, tmp);
! 326: mpn_fft_add_modF(Ap[inc], Ap[0], n);
! 327: mpn_fft_mul_2exp_modF(tp,lk[0]*omega, n, tmp);
! 328: mpn_fft_add_modF(Ap[0], tp, n);
! 329: MPN_COPY(tp, Bp[inc], n+1);
! 330: mpn_fft_mul_2exp_modF(Bp[inc], lk[1]*omega, n, tmp);
! 331: mpn_fft_add_modF(Bp[inc], Bp[0], n);
! 332: mpn_fft_mul_2exp_modF(tp,lk[0]*omega, n, tmp);
! 333: mpn_fft_add_modF(Bp[0], tp, n);
! 334: }
! 335: TMP_FREE(marker);
! 336: }
! 337: }
! 338:
! 339:
! 340: /* a[i] <- a[i]*b[i] mod 2^(n*BITS_PER_MP_LIMB)+1 for 0 <= i < K */
! 341: static void
! 342: #if __STDC__
! 343: mpn_fft_mul_modF_K (mp_limb_t **ap, mp_limb_t **bp, mp_size_t n, int K)
! 344: #else
! 345: mpn_fft_mul_modF_K(ap, bp, n, K)
! 346: mp_limb_t **ap, **bp;
! 347: mp_size_t n;
! 348: int K;
! 349: #endif
! 350: {
! 351: int i;
! 352: int sqr = (ap == bp);
! 353: TMP_DECL(marker);
! 354:
! 355: TMP_MARK(marker);
! 356:
! 357: if (n >= (sqr ? FFT_MODF_SQR_THRESHOLD : FFT_MODF_MUL_THRESHOLD)) {
! 358: int k, K2,nprime2,Nprime2,M2,maxLK,l,Mp2;
! 359: int **_fft_l;
! 360: mp_limb_t **Ap,**Bp,*A,*B,*T;
! 361:
! 362: k = mpn_fft_best_k (n, sqr);
! 363: K2 = 1<<k;
! 364: maxLK = (K2>BITS_PER_MP_LIMB) ? K2 : BITS_PER_MP_LIMB;
! 365: M2 = n*BITS_PER_MP_LIMB/K2;
! 366: l = n/K2;
! 367: Nprime2 = ((2*M2+k+2+maxLK)/maxLK)*maxLK; /* ceil((2*M2+k+3)/maxLK)*maxLK*/
! 368: nprime2 = Nprime2/BITS_PER_MP_LIMB;
! 369: Mp2 = Nprime2/K2;
! 370:
! 371: Ap = TMP_ALLOC_MP_PTRS (K2);
! 372: Bp = TMP_ALLOC_MP_PTRS (K2);
! 373: A = TMP_ALLOC_LIMBS (2*K2*(nprime2+1));
! 374: T = TMP_ALLOC_LIMBS (nprime2+1);
! 375: B = A + K2*(nprime2+1);
! 376: _fft_l = TMP_ALLOC_TYPE (k+1, int*);
! 377: for (i=0;i<=k;i++)
! 378: _fft_l[i] = TMP_ALLOC_TYPE (1<<i, int);
! 379: mpn_fft_initl(_fft_l, k);
! 380:
! 381: TRACE (printf("recurse: %dx%d limbs -> %d times %dx%d (%1.2f)\n", n,
! 382: n, K2, nprime2, nprime2, 2.0*(double)n/nprime2/K2));
! 383:
! 384: for (i=0;i<K;i++,ap++,bp++)
! 385: mpn_mul_fft_internal(*ap, *ap, *bp, n, k, K2, Ap, Bp, A, B, nprime2,
! 386: l, Mp2, _fft_l, T, 1);
! 387: }
! 388: else {
! 389: mp_limb_t *a, *b, cc, *tp, *tpn; int n2=2*n;
! 390: tp = TMP_ALLOC_LIMBS (n2);
! 391: tpn = tp+n;
! 392: TRACE (printf (" mpn_mul_n %d of %d limbs\n", K, n));
! 393: for (i=0;i<K;i++) {
! 394: a = *ap++; b=*bp++;
! 395: if (sqr)
! 396: mpn_sqr_n(tp, a, n);
! 397: else
! 398: mpn_mul_n(tp, b, a, n);
! 399: if (a[n]) cc=mpn_add_n(tpn, tpn, b, n); else cc=0;
! 400: if (b[n]) cc += mpn_add_n(tpn, tpn, a, n) + a[n];
! 401: if (cc) {
! 402: cc = mpn_add_1(tp, tp, n2, cc);
! 403: ASSERT_NOCARRY (mpn_add_1(tp, tp, n2, cc));
! 404: }
! 405: a[n] = mpn_sub_n(a, tp, tpn, n) && mpn_add_1(a, a, n, 1);
! 406: }
! 407: }
! 408: TMP_FREE(marker);
! 409: }
! 410:
! 411:
! 412: /* input: A^[l[k][0]] A^[l[k][1]] ... A^[l[k][K-1]]
! 413: output: K*A[0] K*A[K-1] ... K*A[1] */
! 414:
! 415: static void
! 416: #if __STDC__
! 417: mpn_fft_fftinv (mp_limb_t **Ap, int K, mp_size_t omega, mp_size_t n,
! 418: mp_limb_t *tp)
! 419: #else
! 420: mpn_fft_fftinv(Ap,K,omega,n,tp)
! 421: mp_limb_t **Ap, *tp;
! 422: int K;
! 423: mp_size_t omega, n;
! 424: #endif
! 425: {
! 426: if (K==2) {
! 427: #ifdef ADDSUB
! 428: if (mpn_addsub_n(Ap[0], Ap[1], Ap[0], Ap[1], n+1) & 1)
! 429: #else
! 430: MPN_COPY(tp, Ap[0], n+1);
! 431: mpn_add_n(Ap[0], Ap[0], Ap[1], n+1);
! 432: if (mpn_sub_n(Ap[1], tp, Ap[1], n+1))
! 433: #endif
! 434: Ap[1][n] = mpn_add_1(Ap[1], Ap[1], n, 1);
! 435: }
! 436: else {
! 437: int j, K2=K/2; mp_limb_t **Bp=Ap+K2, *tmp;
! 438: TMP_DECL(marker);
! 439:
! 440: TMP_MARK(marker);
! 441: tmp = TMP_ALLOC_LIMBS (n+1);
! 442: mpn_fft_fftinv(Ap, K2, 2*omega, n, tp);
! 443: mpn_fft_fftinv(Bp, K2, 2*omega, n, tp);
! 444: /* A[j] <- A[j] + omega^j A[j+K/2]
! 445: A[j+K/2] <- A[j] + omega^(j+K/2) A[j+K/2] */
! 446: for (j=0;j<K2;j++,Ap++,Bp++) {
! 447: MPN_COPY(tp, Bp[0], n+1);
! 448: mpn_fft_mul_2exp_modF(Bp[0], (j+K2)*omega, n, tmp);
! 449: mpn_fft_add_modF(Bp[0], Ap[0], n);
! 450: mpn_fft_mul_2exp_modF(tp, j*omega, n, tmp);
! 451: mpn_fft_add_modF(Ap[0], tp, n);
! 452: }
! 453: TMP_FREE(marker);
! 454: }
! 455: }
! 456:
! 457:
! 458: /* A <- A/2^k mod 2^(n*BITS_PER_MP_LIMB)+1 */
! 459: static void
! 460: #if __STDC__
! 461: mpn_fft_div_2exp_modF (mp_limb_t *ap, int k, mp_size_t n, mp_limb_t *tp)
! 462: #else
! 463: mpn_fft_div_2exp_modF(ap,k,n,tp)
! 464: mp_limb_t *ap,*tp;
! 465: int k;
! 466: mp_size_t n;
! 467: #endif
! 468: {
! 469: int i;
! 470:
! 471: i = 2*n*BITS_PER_MP_LIMB;
! 472: i = (i-k) % i;
! 473: mpn_fft_mul_2exp_modF(ap,i,n,tp);
! 474: /* 1/2^k = 2^(2nL-k) mod 2^(n*BITS_PER_MP_LIMB)+1 */
! 475: /* normalize so that A < 2^(n*BITS_PER_MP_LIMB)+1 */
! 476: if (ap[n]==1) {
! 477: for (i=0;i<n && ap[i]==0;i++);
! 478: if (i<n) {
! 479: ap[n]=0;
! 480: mpn_sub_1(ap, ap, n, 1);
! 481: }
! 482: }
! 483: }
! 484:
! 485:
! 486: /* R <- A mod 2^(n*BITS_PER_MP_LIMB)+1, n<=an<=3*n */
! 487: static void
! 488: #if __STDC__
! 489: mpn_fft_norm_modF(mp_limb_t *rp, mp_limb_t *ap, mp_size_t n, mp_size_t an)
! 490: #else
! 491: mpn_fft_norm_modF(rp, ap, n, an)
! 492: mp_limb_t *rp;
! 493: mp_limb_t *ap;
! 494: mp_size_t n;
! 495: mp_size_t an;
! 496: #endif
! 497: {
! 498: mp_size_t l;
! 499:
! 500: if (an>2*n) {
! 501: l = n;
! 502: rp[n] = mpn_add_1(rp+an-2*n, ap+an-2*n, 3*n-an,
! 503: mpn_add_n(rp,ap,ap+2*n,an-2*n));
! 504: }
! 505: else {
! 506: l = an-n;
! 507: MPN_COPY(rp, ap, n);
! 508: rp[n]=0;
! 509: }
! 510: if (mpn_sub_n(rp,rp,ap+n,l)) {
! 511: if (mpn_sub_1(rp+l,rp+l,n+1-l,1))
! 512: rp[n]=mpn_add_1(rp,rp,n,1);
! 513: }
! 514: }
! 515:
! 516:
! 517: static void
! 518: #if __STDC__
! 519: mpn_mul_fft_internal(mp_limb_t *op, mp_srcptr n, mp_srcptr m, mp_size_t pl,
! 520: int k, int K,
! 521: mp_limb_t **Ap, mp_limb_t **Bp,
! 522: mp_limb_t *A, mp_limb_t *B,
! 523: mp_size_t nprime, mp_size_t l, mp_size_t Mp,
! 524: int **_fft_l,
! 525: mp_limb_t *T, int rec)
! 526: #else
! 527: mpn_mul_fft_internal(op,n,m,pl,k,K,Ap,Bp,A,B,nprime,l,Mp,_fft_l,T,rec)
! 528: mp_limb_t *op;
! 529: mp_srcptr n, m;
! 530: mp_limb_t **Ap,**Bp,*A,*B,*T;
! 531: mp_size_t pl,nprime;
! 532: int **_fft_l;
! 533: int k,K,l,Mp,rec;
! 534: #endif
! 535: {
! 536: int i, sqr, pla, lo, sh, j;
! 537: mp_limb_t *p;
! 538:
! 539: sqr = (n==m);
! 540:
! 541: TRACE (printf ("pl=%d k=%d K=%d np=%d l=%d Mp=%d rec=%d sqr=%d\n",
! 542: pl,k,K,nprime,l,Mp,rec,sqr));
! 543:
! 544: /* decomposition of inputs into arrays Ap[i] and Bp[i] */
! 545: if (rec) for (i=0;i<K;i++) {
! 546: Ap[i] = A+i*(nprime+1); Bp[i] = B+i*(nprime+1);
! 547: /* store the next M bits of n into A[i] */
! 548: /* supposes that M is a multiple of BITS_PER_MP_LIMB */
! 549: MPN_COPY(Ap[i], n, l); n+=l; MPN_ZERO(Ap[i]+l, nprime+1-l);
! 550: /* set most significant bits of n and m (important in recursive calls) */
! 551: if (i==K-1) Ap[i][l]=n[0];
! 552: mpn_fft_mul_2exp_modF(Ap[i], i*Mp, nprime, T);
! 553: if (!sqr) {
! 554: MPN_COPY(Bp[i], m, l); m+=l; MPN_ZERO(Bp[i]+l, nprime+1-l);
! 555: if (i==K-1) Bp[i][l]=m[0];
! 556: mpn_fft_mul_2exp_modF(Bp[i], i*Mp, nprime, T);
! 557: }
! 558: }
! 559:
! 560: /* direct fft's */
! 561: if (sqr) mpn_fft_fft_sqr(Ap,K,_fft_l+k,2*Mp,nprime,1, T);
! 562: else mpn_fft_fft(Ap,Bp,K,_fft_l+k,2*Mp,nprime,1, T);
! 563:
! 564: /* term to term multiplications */
! 565: mpn_fft_mul_modF_K(Ap, (sqr) ? Ap : Bp, nprime, K);
! 566:
! 567: /* inverse fft's */
! 568: mpn_fft_fftinv(Ap, K, 2*Mp, nprime, T);
! 569:
! 570: /* division of terms after inverse fft */
! 571: for (i=0;i<K;i++) mpn_fft_div_2exp_modF(Ap[i],k+((K-i)%K)*Mp,nprime, T);
! 572:
! 573: /* addition of terms in result p */
! 574: MPN_ZERO(T,nprime+1);
! 575: pla = l*(K-1)+nprime+1; /* number of required limbs for p */
! 576: p = B; /* B has K*(n'+1) limbs, which is >= pla, i.e. enough */
! 577: MPN_ZERO(p, pla);
! 578: sqr=0; /* will accumulate the (signed) carry at p[pla] */
! 579: for (i=K-1,lo=l*i+nprime,sh=l*i;i>=0;i--,lo-=l,sh-=l) {
! 580: mp_ptr n = p+sh;
! 581: j = (K-i)%K;
! 582: if (mpn_add_n(n,n,Ap[j],nprime+1))
! 583: sqr += mpn_add_1(n+nprime+1,n+nprime+1,pla-sh-nprime-1,1);
! 584: T[2*l]=i+1; /* T = (i+1)*2^(2*M) */
! 585: if (mpn_cmp(Ap[j],T,nprime+1)>0) { /* subtract 2^N'+1 */
! 586: sqr -= mpn_sub_1(n,n,pla-sh,1);
! 587: sqr -= mpn_sub_1(p+lo,p+lo,pla-lo,1);
! 588: }
! 589: }
! 590: if (sqr==-1) {
! 591: if ((sqr=mpn_add_1(p+pla-pl,p+pla-pl,pl,1))) {
! 592: /* p[pla-pl]...p[pla-1] are all zero */
! 593: mpn_sub_1(p+pla-pl-1,p+pla-pl-1,pl+1,1);
! 594: mpn_sub_1(p+pla-1,p+pla-1,1,1);
! 595: }
! 596: }
! 597: else if (sqr==1) {
! 598: if (pla>=2*pl)
! 599: while ((sqr=mpn_add_1(p+pla-2*pl,p+pla-2*pl,2*pl,sqr)));
! 600: else {
! 601: sqr = mpn_sub_1(p+pla-pl,p+pla-pl,pl,sqr);
! 602: ASSERT (sqr == 0);
! 603: }
! 604: }
! 605: else
! 606: ASSERT (sqr == 0);
! 607:
! 608: /* here p < 2^(2M) [K 2^(M(K-1)) + (K-1) 2^(M(K-2)) + ... ]
! 609: < K 2^(2M) [2^(M(K-1)) + 2^(M(K-2)) + ... ]
! 610: < K 2^(2M) 2^(M(K-1))*2 = 2^(M*K+M+k+1) */
! 611: mpn_fft_norm_modF(op,p,pl,pla);
! 612: }
! 613:
! 614:
! 615: /* op <- n*m mod 2^N+1 with fft of size 2^k where N=pl*BITS_PER_MP_LIMB
! 616: n and m have respectively nl and ml limbs
! 617: op must have space for pl+1 limbs
! 618: One must have pl = mpn_fft_next_size(pl, k).
! 619: */
! 620:
! 621: void
! 622: #if __STDC__
! 623: mpn_mul_fft (mp_ptr op, mp_size_t pl,
! 624: mp_srcptr n, mp_size_t nl,
! 625: mp_srcptr m, mp_size_t ml,
! 626: int k)
! 627: #else
! 628: mpn_mul_fft (op, pl, n, nl, m, ml, k)
! 629: mp_ptr op;
! 630: mp_size_t pl;
! 631: mp_srcptr n;
! 632: mp_size_t nl;
! 633: mp_srcptr m;
! 634: mp_size_t ml;
! 635: int k;
! 636: #endif
! 637: {
! 638: int K,maxLK,i,j;
! 639: mp_size_t N,Nprime,nprime,M,Mp,l;
! 640: mp_limb_t **Ap,**Bp,*A,*T,*B;
! 641: int **_fft_l;
! 642: int sqr = (n==m && nl==ml);
! 643: TMP_DECL(marker);
! 644:
! 645: TRACE (printf ("\nmpn_mul_fft pl=%ld nl=%ld ml=%ld k=%d\n",
! 646: pl, nl, ml, k));
! 647: ASSERT_ALWAYS (mpn_fft_next_size(pl, k) == pl);
! 648:
! 649: TMP_MARK(marker);
! 650: N = pl*BITS_PER_MP_LIMB;
! 651: _fft_l = TMP_ALLOC_TYPE (k+1, int*);
! 652: for (i=0;i<=k;i++)
! 653: _fft_l[i] = TMP_ALLOC_TYPE (1<<i, int);
! 654: mpn_fft_initl(_fft_l, k);
! 655: K = 1<<k;
! 656: M = N/K; /* N = 2^k M */
! 657: l = M/BITS_PER_MP_LIMB;
! 658: maxLK = (K>BITS_PER_MP_LIMB) ? K : BITS_PER_MP_LIMB;
! 659:
! 660: Nprime = ((2*M+k+2+maxLK)/maxLK)*maxLK; /* ceil((2*M+k+3)/maxLK)*maxLK; */
! 661: nprime = Nprime/BITS_PER_MP_LIMB;
! 662: TRACE (printf ("N=%d K=%d, M=%d, l=%d, maxLK=%d, Np=%d, np=%d\n",
! 663: N, K, M, l, maxLK, Nprime, nprime));
! 664: if (nprime >= (sqr ? FFT_MODF_SQR_THRESHOLD : FFT_MODF_MUL_THRESHOLD)) {
! 665: maxLK = (1<<mpn_fft_best_k(nprime,n==m))*BITS_PER_MP_LIMB;
! 666: if (Nprime % maxLK) {
! 667: Nprime=((Nprime/maxLK)+1)*maxLK;
! 668: nprime = Nprime/BITS_PER_MP_LIMB;
! 669: }
! 670: TRACE (printf ("new maxLK=%d, Np=%d, np=%d\n", maxLK, Nprime, nprime));
! 671: }
! 672:
! 673: T = TMP_ALLOC_LIMBS (nprime+1);
! 674: Mp = Nprime/K;
! 675:
! 676: TRACE (printf("%dx%d limbs -> %d times %dx%d limbs (%1.2f)\n",
! 677: pl,pl,K,nprime,nprime,2.0*(double)N/Nprime/K);
! 678: printf(" temp space %ld\n", 2*K*(nprime+1)));
! 679:
! 680: A = _MP_ALLOCATE_FUNC_LIMBS (2*K*(nprime+1));
! 681: B = A+K*(nprime+1);
! 682: Ap = TMP_ALLOC_MP_PTRS (K);
! 683: Bp = TMP_ALLOC_MP_PTRS (K);
! 684: /* special decomposition for main call */
! 685: for (i=0;i<K;i++) {
! 686: Ap[i] = A+i*(nprime+1); Bp[i] = B+i*(nprime+1);
! 687: /* store the next M bits of n into A[i] */
! 688: /* supposes that M is a multiple of BITS_PER_MP_LIMB */
! 689: if (nl>0) {
! 690: j = (nl>=l) ? l : nl; /* limbs to store in Ap[i] */
! 691: MPN_COPY(Ap[i], n, j); n+=l; MPN_ZERO(Ap[i]+j, nprime+1-j);
! 692: mpn_fft_mul_2exp_modF(Ap[i], i*Mp, nprime, T);
! 693: }
! 694: else MPN_ZERO(Ap[i], nprime+1);
! 695: nl -= l;
! 696: if (n!=m) {
! 697: if (ml>0) {
! 698: j = (ml>=l) ? l : ml; /* limbs to store in Bp[i] */
! 699: MPN_COPY(Bp[i], m, j); m+=l; MPN_ZERO(Bp[i]+j, nprime+1-j);
! 700: mpn_fft_mul_2exp_modF(Bp[i], i*Mp, nprime, T);
! 701: }
! 702: else MPN_ZERO(Bp[i], nprime+1);
! 703: }
! 704: ml -= l;
! 705: }
! 706: mpn_mul_fft_internal(op,n,m,pl,k,K,Ap,Bp,A,B,nprime,l,Mp,_fft_l,T,0);
! 707: TMP_FREE(marker);
! 708: _MP_FREE_FUNC_LIMBS (A, 2*K*(nprime+1));
! 709: }
! 710:
! 711:
! 712: #if WANT_ASSERT
! 713: static int
! 714: #if __STDC__
! 715: mpn_zero_p (mp_ptr p, mp_size_t n)
! 716: #else
! 717: mpn_zero_p (p, n)
! 718: mp_ptr p;
! 719: mp_size_t n;
! 720: #endif
! 721: {
! 722: mp_size_t i;
! 723:
! 724: for (i = 0; i < n; i++)
! 725: {
! 726: if (p[i] != 0)
! 727: return 0;
! 728: }
! 729:
! 730: return 1;
! 731: }
! 732: #endif
! 733:
! 734:
! 735: /* Multiply {n,nl}*{m,ml} and write the result to {op,nl+ml}.
! 736:
! 737: FIXME: Duplicating the result like this is wasteful, do something better
! 738: perhaps at the norm_modF stage above. */
! 739:
! 740: void
! 741: #if __STDC__
! 742: mpn_mul_fft_full (mp_ptr op,
! 743: mp_srcptr n, mp_size_t nl,
! 744: mp_srcptr m, mp_size_t ml)
! 745: #else
! 746: mpn_mul_fft_full (op, n, nl, m, ml)
! 747: mp_ptr op;
! 748: mp_srcptr n;
! 749: mp_size_t nl;
! 750: mp_srcptr m;
! 751: mp_size_t ml;
! 752: #endif
! 753: {
! 754: mp_ptr pad_op;
! 755: mp_size_t pl;
! 756: int k;
! 757: int sqr = (n==m && nl==ml);
! 758:
! 759: k = mpn_fft_best_k (nl+ml, sqr);
! 760: pl = mpn_fft_next_size (nl+ml, k);
! 761:
! 762: TRACE (printf ("mpn_mul_fft_full nl=%ld ml=%ld -> pl=%ld k=%d\n",
! 763: nl, ml, pl, k));
! 764:
! 765: pad_op = _MP_ALLOCATE_FUNC_LIMBS (pl+1);
! 766: mpn_mul_fft (pad_op, pl, n, nl, m, ml, k);
! 767:
! 768: ASSERT (mpn_zero_p (pad_op+nl+ml, pl+1-(nl+ml)));
! 769: MPN_COPY (op, pad_op, nl+ml);
! 770:
! 771: _MP_FREE_FUNC_LIMBS (pad_op, pl+1);
! 772: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to mul_fft.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpn / generic