+++ /dev/null
-/* An implementation in GMP of Scho"nhage's fast multiplication algorithm
- modulo 2^N+1, by Paul Zimmermann, INRIA Lorraine, February 1998.
-
- THE CONTENTS OF THIS FILE ARE FOR INTERNAL USE AND THE FUNCTIONS HAVE
- MUTABLE INTERFACES. IT IS ONLY SAFE TO REACH THEM THROUGH DOCUMENTED
- INTERFACES. IT IS ALMOST GUARANTEED THAT THEY'LL CHANGE OR DISAPPEAR IN
- A FUTURE GNU MP RELEASE.
-
-Copyright (C) 1998, 1999, 2000 Free Software Foundation, Inc.
-
-This file is part of the GNU MP Library.
-
-The GNU MP Library is free software; you can redistribute it and/or modify
-it under the terms of the GNU Lesser General Public License as published by
-the Free Software Foundation; either version 2.1 of the License, or (at your
-option) any later version.
-
-The GNU MP Library is distributed in the hope that it will be useful, but
-WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
-or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
-License for more details.
-
-You should have received a copy of the GNU Lesser General Public License
-along with the GNU MP Library; see the file COPYING.LIB. If not, write to
-the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
-MA 02111-1307, USA. */
-
-
-/* References:
-
- Schnelle Multiplikation grosser Zahlen, by Arnold Scho"nhage and Volker
- Strassen, Computing 7, p. 281-292, 1971.
-
- Asymptotically fast algorithms for the numerical multiplication
- and division of polynomials with complex coefficients, by Arnold Scho"nhage,
- Computer Algebra, EUROCAM'82, LNCS 144, p. 3-15, 1982.
-
- Tapes versus Pointers, a study in implementing fast algorithms,
- by Arnold Scho"nhage, Bulletin of the EATCS, 30, p. 23-32, 1986.
-
- See also http://www.loria.fr/~zimmerma/bignum
-
-
- Future:
-
- K==2 isn't needed in the current uses of this code and the bits specific
- for that could be dropped.
-
- It might be possible to avoid a small number of MPN_COPYs by using a
- rotating temporary or two.
-
- Multiplications of unequal sized operands can be done with this code, but
- it needs a tighter test for identifying squaring (same sizes as well as
- same pointers). */
-
-
-#include <stdio.h>
-#include "gmp.h"
-#include "gmp-impl.h"
-
-
-/* Change this to "#define TRACE(x) x" for some traces. */
-#define TRACE(x)
-
-
-
-FFT_TABLE_ATTRS mp_size_t mpn_fft_table[2][MPN_FFT_TABLE_SIZE] = {
- FFT_MUL_TABLE,
- FFT_SQR_TABLE
-};
-
-
-static void mpn_mul_fft_internal
-_PROTO ((mp_limb_t *op, mp_srcptr n, mp_srcptr m, mp_size_t pl,
- int k, int K,
- mp_limb_t **Ap, mp_limb_t **Bp,
- mp_limb_t *A, mp_limb_t *B,
- mp_size_t nprime, mp_size_t l, mp_size_t Mp, int **_fft_l,
- mp_limb_t *T, int rec));
-
-
-/* Find the best k to use for a mod 2^(n*BITS_PER_MP_LIMB)+1 FFT.
- sqr==0 if for a multiply, sqr==1 for a square */
-int
-#if __STDC__
-mpn_fft_best_k (mp_size_t n, int sqr)
-#else
-mpn_fft_best_k (n, sqr)
- mp_size_t n;
- int sqr;
-#endif
-{
- mp_size_t t;
- int i;
-
- for (i = 0; mpn_fft_table[sqr][i] != 0; i++)
- if (n < mpn_fft_table[sqr][i])
- return i + FFT_FIRST_K;
-
- /* treat 4*last as one further entry */
- if (i == 0 || n < 4*mpn_fft_table[sqr][i-1])
- return i + FFT_FIRST_K;
- else
- return i + FFT_FIRST_K + 1;
-}
-
-
-/* Returns smallest possible number of limbs >= pl for a fft of size 2^k.
- FIXME: Is this simply pl rounded up to the next multiple of 2^k ? */
-
-mp_size_t
-#if __STDC__
-mpn_fft_next_size (mp_size_t pl, int k)
-#else
-mpn_fft_next_size (pl, k)
- mp_size_t pl;
- int k;
-#endif
-{
- mp_size_t N, M;
- int K;
-
- /* if (k==0) k = mpn_fft_best_k (pl, sqr); */
- N = pl*BITS_PER_MP_LIMB;
- K = 1<<k;
- if (N%K) N=(N/K+1)*K;
- M = N/K;
- if (M%BITS_PER_MP_LIMB) N=((M/BITS_PER_MP_LIMB)+1)*BITS_PER_MP_LIMB*K;
- return (N/BITS_PER_MP_LIMB);
-}
-
-
-static void
-#if __STDC__
-mpn_fft_initl(int **l, int k)
-#else
-mpn_fft_initl(l, k)
- int **l;
- int k;
-#endif
-{
- int i,j,K;
-
- l[0][0] = 0;
- for (i=1,K=2;i<=k;i++,K*=2) {
- for (j=0;j<K/2;j++) {
- l[i][j] = 2*l[i-1][j];
- l[i][K/2+j] = 1+l[i][j];
- }
- }
-}
-
-
-/* a <- -a mod 2^(n*BITS_PER_MP_LIMB)+1 */
-static void
-#if __STDC__
-mpn_fft_neg_modF(mp_limb_t *ap, mp_size_t n)
-#else
-mpn_fft_neg_modF(ap, n)
- mp_limb_t *ap;
- mp_size_t n;
-#endif
-{
- mp_limb_t c;
-
- c = ap[n]+2;
- mpn_com_n (ap, ap, n);
- ap[n]=0; mpn_incr_u(ap, c);
-}
-
-
-/* a <- a*2^e mod 2^(n*BITS_PER_MP_LIMB)+1 */
-static void
-#if __STDC__
-mpn_fft_mul_2exp_modF(mp_limb_t *ap, int e, mp_size_t n, mp_limb_t *tp)
-#else
-mpn_fft_mul_2exp_modF(ap, e, n, tp)
- mp_limb_t *ap;
- int e;
- mp_size_t n;
- mp_limb_t *tp;
-#endif
-{
- int d, sh, i; mp_limb_t cc;
-
- d = e%(n*BITS_PER_MP_LIMB); /* 2^e = (+/-) 2^d */
- sh = d % BITS_PER_MP_LIMB;
- if (sh) mpn_lshift(tp, ap, n+1, sh); /* no carry here */
- else MPN_COPY(tp, ap, n+1);
- d /= BITS_PER_MP_LIMB; /* now shift of d limbs to the left */
- if (d) {
- /* ap[d..n-1] = tp[0..n-d-1], ap[0..d-1] = -tp[n-d..n-1] */
- /* mpn_xor would be more efficient here */
- for (i=d-1;i>=0;i--) ap[i] = ~tp[n-d+i];
- cc = 1-mpn_add_1(ap, ap, d, 1);
- if (cc) cc=mpn_sub_1(ap+d, tp, n-d, 1);
- else MPN_COPY(ap+d, tp, n-d);
- if (cc+=mpn_sub_1(ap+d, ap+d, n-d, tp[n]))
- ap[n]=mpn_add_1(ap, ap, n, cc);
- else ap[n]=0;
- }
- else if ((ap[n]=mpn_sub_1(ap, tp, n, tp[n]))) {
- ap[n]=mpn_add_1(ap, ap, n, 1);
- }
- if ((e/(n*BITS_PER_MP_LIMB))%2) mpn_fft_neg_modF(ap, n);
-}
-
-
-/* a <- a+b mod 2^(n*BITS_PER_MP_LIMB)+1 */
-static void
-#if __STDC__
-mpn_fft_add_modF (mp_limb_t *ap, mp_limb_t *bp, int n)
-#else
-mpn_fft_add_modF (ap, bp, n)
- mp_limb_t *ap,*bp;
- int n;
-#endif
-{
- mp_limb_t c;
-
- c = ap[n] + bp[n] + mpn_add_n(ap, ap, bp, n);
- if (c>1) c -= 1+mpn_sub_1(ap,ap,n,1);
- ap[n]=c;
-}
-
-
-/* input: A[0] ... A[inc*(K-1)] are residues mod 2^N+1 where
- N=n*BITS_PER_MP_LIMB
- 2^omega is a primitive root mod 2^N+1
- output: A[inc*l[k][i]] <- \sum (2^omega)^(ij) A[inc*j] mod 2^N+1 */
-
-static void
-#if __STDC__
-mpn_fft_fft_sqr (mp_limb_t **Ap, mp_size_t K, int **ll,
- mp_size_t omega, mp_size_t n, mp_size_t inc, mp_limb_t *tp)
-#else
-mpn_fft_fft_sqr(Ap,K,ll,omega,n,inc,tp)
-mp_limb_t **Ap,*tp;
-mp_size_t K,omega,n,inc;
-int **ll;
-#endif
-{
- if (K==2) {
-#ifdef ADDSUB
- if (mpn_addsub_n(Ap[0], Ap[inc], Ap[0], Ap[inc], n+1) & 1)
-#else
- MPN_COPY(tp, Ap[0], n+1);
- mpn_add_n(Ap[0], Ap[0], Ap[inc],n+1);
- if (mpn_sub_n(Ap[inc], tp, Ap[inc],n+1))
-#endif
- Ap[inc][n] = mpn_add_1(Ap[inc], Ap[inc], n, 1);
- }
- else {
- int j, inc2=2*inc;
- int *lk = *ll;
- mp_limb_t *tmp;
- TMP_DECL(marker);
-
- TMP_MARK(marker);
- tmp = TMP_ALLOC_LIMBS (n+1);
- mpn_fft_fft_sqr(Ap, K/2,ll-1,2*omega,n,inc2, tp);
- mpn_fft_fft_sqr(Ap+inc, K/2,ll-1,2*omega,n,inc2, tp);
- /* A[2*j*inc] <- A[2*j*inc] + omega^l[k][2*j*inc] A[(2j+1)inc]
- A[(2j+1)inc] <- A[2*j*inc] + omega^l[k][(2j+1)inc] A[(2j+1)inc] */
- for (j=0;j<K/2;j++,lk+=2,Ap+=2*inc) {
- MPN_COPY(tp, Ap[inc], n+1);
- mpn_fft_mul_2exp_modF(Ap[inc], lk[1]*omega, n, tmp);
- mpn_fft_add_modF(Ap[inc], Ap[0], n);
- mpn_fft_mul_2exp_modF(tp,lk[0]*omega, n, tmp);
- mpn_fft_add_modF(Ap[0], tp, n);
- }
- TMP_FREE(marker);
- }
-}
-
-
-/* input: A[0] ... A[inc*(K-1)] are residues mod 2^N+1 where
- N=n*BITS_PER_MP_LIMB
- 2^omega is a primitive root mod 2^N+1
- output: A[inc*l[k][i]] <- \sum (2^omega)^(ij) A[inc*j] mod 2^N+1 */
-
-static void
-#if __STDC__
-mpn_fft_fft (mp_limb_t **Ap, mp_limb_t **Bp, mp_size_t K, int **ll,
- mp_size_t omega, mp_size_t n, mp_size_t inc, mp_limb_t *tp)
-#else
-mpn_fft_fft(Ap,Bp,K,ll,omega,n,inc,tp)
- mp_limb_t **Ap,**Bp,*tp;
- mp_size_t K,omega,n,inc;
- int **ll;
-#endif
-{
- if (K==2) {
-#ifdef ADDSUB
- if (mpn_addsub_n(Ap[0], Ap[inc], Ap[0], Ap[inc], n+1) & 1)
-#else
- MPN_COPY(tp, Ap[0], n+1);
- mpn_add_n(Ap[0], Ap[0], Ap[inc],n+1);
- if (mpn_sub_n(Ap[inc], tp, Ap[inc],n+1))
-#endif
- Ap[inc][n] = mpn_add_1(Ap[inc], Ap[inc], n, 1);
-#ifdef ADDSUB
- if (mpn_addsub_n(Bp[0], Bp[inc], Bp[0], Bp[inc], n+1) & 1)
-#else
- MPN_COPY(tp, Bp[0], n+1);
- mpn_add_n(Bp[0], Bp[0], Bp[inc],n+1);
- if (mpn_sub_n(Bp[inc], tp, Bp[inc],n+1))
-#endif
- Bp[inc][n] = mpn_add_1(Bp[inc], Bp[inc], n, 1);
- }
- else {
- int j, inc2=2*inc;
- int *lk=*ll;
- mp_limb_t *tmp;
- TMP_DECL(marker);
-
- TMP_MARK(marker);
- tmp = TMP_ALLOC_LIMBS (n+1);
- mpn_fft_fft(Ap, Bp, K/2,ll-1,2*omega,n,inc2, tp);
- mpn_fft_fft(Ap+inc, Bp+inc, K/2,ll-1,2*omega,n,inc2, tp);
- /* A[2*j*inc] <- A[2*j*inc] + omega^l[k][2*j*inc] A[(2j+1)inc]
- A[(2j+1)inc] <- A[2*j*inc] + omega^l[k][(2j+1)inc] A[(2j+1)inc] */
- for (j=0;j<K/2;j++,lk+=2,Ap+=2*inc,Bp+=2*inc) {
- MPN_COPY(tp, Ap[inc], n+1);
- mpn_fft_mul_2exp_modF(Ap[inc], lk[1]*omega, n, tmp);
- mpn_fft_add_modF(Ap[inc], Ap[0], n);
- mpn_fft_mul_2exp_modF(tp,lk[0]*omega, n, tmp);
- mpn_fft_add_modF(Ap[0], tp, n);
- MPN_COPY(tp, Bp[inc], n+1);
- mpn_fft_mul_2exp_modF(Bp[inc], lk[1]*omega, n, tmp);
- mpn_fft_add_modF(Bp[inc], Bp[0], n);
- mpn_fft_mul_2exp_modF(tp,lk[0]*omega, n, tmp);
- mpn_fft_add_modF(Bp[0], tp, n);
- }
- TMP_FREE(marker);
- }
-}
-
-
-/* a[i] <- a[i]*b[i] mod 2^(n*BITS_PER_MP_LIMB)+1 for 0 <= i < K */
-static void
-#if __STDC__
-mpn_fft_mul_modF_K (mp_limb_t **ap, mp_limb_t **bp, mp_size_t n, int K)
-#else
-mpn_fft_mul_modF_K(ap, bp, n, K)
- mp_limb_t **ap, **bp;
- mp_size_t n;
- int K;
-#endif
-{
- int i;
- int sqr = (ap == bp);
- TMP_DECL(marker);
-
- TMP_MARK(marker);
-
- if (n >= (sqr ? FFT_MODF_SQR_THRESHOLD : FFT_MODF_MUL_THRESHOLD)) {
- int k, K2,nprime2,Nprime2,M2,maxLK,l,Mp2;
- int **_fft_l;
- mp_limb_t **Ap,**Bp,*A,*B,*T;
-
- k = mpn_fft_best_k (n, sqr);
- K2 = 1<<k;
- maxLK = (K2>BITS_PER_MP_LIMB) ? K2 : BITS_PER_MP_LIMB;
- M2 = n*BITS_PER_MP_LIMB/K2;
- l = n/K2;
- Nprime2 = ((2*M2+k+2+maxLK)/maxLK)*maxLK; /* ceil((2*M2+k+3)/maxLK)*maxLK*/
- nprime2 = Nprime2/BITS_PER_MP_LIMB;
- Mp2 = Nprime2/K2;
-
- Ap = TMP_ALLOC_MP_PTRS (K2);
- Bp = TMP_ALLOC_MP_PTRS (K2);
- A = TMP_ALLOC_LIMBS (2*K2*(nprime2+1));
- T = TMP_ALLOC_LIMBS (nprime2+1);
- B = A + K2*(nprime2+1);
- _fft_l = TMP_ALLOC_TYPE (k+1, int*);
- for (i=0;i<=k;i++)
- _fft_l[i] = TMP_ALLOC_TYPE (1<<i, int);
- mpn_fft_initl(_fft_l, k);
-
- TRACE (printf("recurse: %dx%d limbs -> %d times %dx%d (%1.2f)\n", n,
- n, K2, nprime2, nprime2, 2.0*(double)n/nprime2/K2));
-
- for (i=0;i<K;i++,ap++,bp++)
- mpn_mul_fft_internal(*ap, *ap, *bp, n, k, K2, Ap, Bp, A, B, nprime2,
- l, Mp2, _fft_l, T, 1);
- }
- else {
- mp_limb_t *a, *b, cc, *tp, *tpn; int n2=2*n;
- tp = TMP_ALLOC_LIMBS (n2);
- tpn = tp+n;
- TRACE (printf (" mpn_mul_n %d of %d limbs\n", K, n));
- for (i=0;i<K;i++) {
- a = *ap++; b=*bp++;
- if (sqr)
- mpn_sqr_n(tp, a, n);
- else
- mpn_mul_n(tp, b, a, n);
- if (a[n]) cc=mpn_add_n(tpn, tpn, b, n); else cc=0;
- if (b[n]) cc += mpn_add_n(tpn, tpn, a, n) + a[n];
- if (cc) {
- cc = mpn_add_1(tp, tp, n2, cc);
- ASSERT_NOCARRY (mpn_add_1(tp, tp, n2, cc));
- }
- a[n] = mpn_sub_n(a, tp, tpn, n) && mpn_add_1(a, a, n, 1);
- }
- }
- TMP_FREE(marker);
-}
-
-
-/* input: A^[l[k][0]] A^[l[k][1]] ... A^[l[k][K-1]]
- output: K*A[0] K*A[K-1] ... K*A[1] */
-
-static void
-#if __STDC__
-mpn_fft_fftinv (mp_limb_t **Ap, int K, mp_size_t omega, mp_size_t n,
- mp_limb_t *tp)
-#else
-mpn_fft_fftinv(Ap,K,omega,n,tp)
- mp_limb_t **Ap, *tp;
- int K;
- mp_size_t omega, n;
-#endif
-{
- if (K==2) {
-#ifdef ADDSUB
- if (mpn_addsub_n(Ap[0], Ap[1], Ap[0], Ap[1], n+1) & 1)
-#else
- MPN_COPY(tp, Ap[0], n+1);
- mpn_add_n(Ap[0], Ap[0], Ap[1], n+1);
- if (mpn_sub_n(Ap[1], tp, Ap[1], n+1))
-#endif
- Ap[1][n] = mpn_add_1(Ap[1], Ap[1], n, 1);
- }
- else {
- int j, K2=K/2; mp_limb_t **Bp=Ap+K2, *tmp;
- TMP_DECL(marker);
-
- TMP_MARK(marker);
- tmp = TMP_ALLOC_LIMBS (n+1);
- mpn_fft_fftinv(Ap, K2, 2*omega, n, tp);
- mpn_fft_fftinv(Bp, K2, 2*omega, n, tp);
- /* A[j] <- A[j] + omega^j A[j+K/2]
- A[j+K/2] <- A[j] + omega^(j+K/2) A[j+K/2] */
- for (j=0;j<K2;j++,Ap++,Bp++) {
- MPN_COPY(tp, Bp[0], n+1);
- mpn_fft_mul_2exp_modF(Bp[0], (j+K2)*omega, n, tmp);
- mpn_fft_add_modF(Bp[0], Ap[0], n);
- mpn_fft_mul_2exp_modF(tp, j*omega, n, tmp);
- mpn_fft_add_modF(Ap[0], tp, n);
- }
- TMP_FREE(marker);
- }
-}
-
-
-/* A <- A/2^k mod 2^(n*BITS_PER_MP_LIMB)+1 */
-static void
-#if __STDC__
-mpn_fft_div_2exp_modF (mp_limb_t *ap, int k, mp_size_t n, mp_limb_t *tp)
-#else
-mpn_fft_div_2exp_modF(ap,k,n,tp)
- mp_limb_t *ap,*tp;
- int k;
- mp_size_t n;
-#endif
-{
- int i;
-
- i = 2*n*BITS_PER_MP_LIMB;
- i = (i-k) % i;
- mpn_fft_mul_2exp_modF(ap,i,n,tp);
- /* 1/2^k = 2^(2nL-k) mod 2^(n*BITS_PER_MP_LIMB)+1 */
- /* normalize so that A < 2^(n*BITS_PER_MP_LIMB)+1 */
- if (ap[n]==1) {
- for (i=0;i<n && ap[i]==0;i++);
- if (i<n) {
- ap[n]=0;
- mpn_sub_1(ap, ap, n, 1);
- }
- }
-}
-
-
-/* R <- A mod 2^(n*BITS_PER_MP_LIMB)+1, n<=an<=3*n */
-static void
-#if __STDC__
-mpn_fft_norm_modF(mp_limb_t *rp, mp_limb_t *ap, mp_size_t n, mp_size_t an)
-#else
-mpn_fft_norm_modF(rp, ap, n, an)
- mp_limb_t *rp;
- mp_limb_t *ap;
- mp_size_t n;
- mp_size_t an;
-#endif
-{
- mp_size_t l;
-
- if (an>2*n) {
- l = n;
- rp[n] = mpn_add_1(rp+an-2*n, ap+an-2*n, 3*n-an,
- mpn_add_n(rp,ap,ap+2*n,an-2*n));
- }
- else {
- l = an-n;
- MPN_COPY(rp, ap, n);
- rp[n]=0;
- }
- if (mpn_sub_n(rp,rp,ap+n,l)) {
- if (mpn_sub_1(rp+l,rp+l,n+1-l,1))
- rp[n]=mpn_add_1(rp,rp,n,1);
- }
-}
-
-
-static void
-#if __STDC__
-mpn_mul_fft_internal(mp_limb_t *op, mp_srcptr n, mp_srcptr m, mp_size_t pl,
- int k, int K,
- mp_limb_t **Ap, mp_limb_t **Bp,
- mp_limb_t *A, mp_limb_t *B,
- mp_size_t nprime, mp_size_t l, mp_size_t Mp,
- int **_fft_l,
- mp_limb_t *T, int rec)
-#else
-mpn_mul_fft_internal(op,n,m,pl,k,K,Ap,Bp,A,B,nprime,l,Mp,_fft_l,T,rec)
- mp_limb_t *op;
- mp_srcptr n, m;
- mp_limb_t **Ap,**Bp,*A,*B,*T;
- mp_size_t pl,nprime;
- int **_fft_l;
- int k,K,l,Mp,rec;
-#endif
-{
- int i, sqr, pla, lo, sh, j;
- mp_limb_t *p;
-
- sqr = (n==m);
-
- TRACE (printf ("pl=%d k=%d K=%d np=%d l=%d Mp=%d rec=%d sqr=%d\n",
- pl,k,K,nprime,l,Mp,rec,sqr));
-
- /* decomposition of inputs into arrays Ap[i] and Bp[i] */
- if (rec) for (i=0;i<K;i++) {
- Ap[i] = A+i*(nprime+1); Bp[i] = B+i*(nprime+1);
- /* store the next M bits of n into A[i] */
- /* supposes that M is a multiple of BITS_PER_MP_LIMB */
- MPN_COPY(Ap[i], n, l); n+=l; MPN_ZERO(Ap[i]+l, nprime+1-l);
- /* set most significant bits of n and m (important in recursive calls) */
- if (i==K-1) Ap[i][l]=n[0];
- mpn_fft_mul_2exp_modF(Ap[i], i*Mp, nprime, T);
- if (!sqr) {
- MPN_COPY(Bp[i], m, l); m+=l; MPN_ZERO(Bp[i]+l, nprime+1-l);
- if (i==K-1) Bp[i][l]=m[0];
- mpn_fft_mul_2exp_modF(Bp[i], i*Mp, nprime, T);
- }
- }
-
- /* direct fft's */
- if (sqr) mpn_fft_fft_sqr(Ap,K,_fft_l+k,2*Mp,nprime,1, T);
- else mpn_fft_fft(Ap,Bp,K,_fft_l+k,2*Mp,nprime,1, T);
-
- /* term to term multiplications */
- mpn_fft_mul_modF_K(Ap, (sqr) ? Ap : Bp, nprime, K);
-
- /* inverse fft's */
- mpn_fft_fftinv(Ap, K, 2*Mp, nprime, T);
-
- /* division of terms after inverse fft */
- for (i=0;i<K;i++) mpn_fft_div_2exp_modF(Ap[i],k+((K-i)%K)*Mp,nprime, T);
-
- /* addition of terms in result p */
- MPN_ZERO(T,nprime+1);
- pla = l*(K-1)+nprime+1; /* number of required limbs for p */
- p = B; /* B has K*(n'+1) limbs, which is >= pla, i.e. enough */
- MPN_ZERO(p, pla);
- sqr=0; /* will accumulate the (signed) carry at p[pla] */
- for (i=K-1,lo=l*i+nprime,sh=l*i;i>=0;i--,lo-=l,sh-=l) {
- mp_ptr n = p+sh;
- j = (K-i)%K;
- if (mpn_add_n(n,n,Ap[j],nprime+1))
- sqr += mpn_add_1(n+nprime+1,n+nprime+1,pla-sh-nprime-1,1);
- T[2*l]=i+1; /* T = (i+1)*2^(2*M) */
- if (mpn_cmp(Ap[j],T,nprime+1)>0) { /* subtract 2^N'+1 */
- sqr -= mpn_sub_1(n,n,pla-sh,1);
- sqr -= mpn_sub_1(p+lo,p+lo,pla-lo,1);
- }
- }
- if (sqr==-1) {
- if ((sqr=mpn_add_1(p+pla-pl,p+pla-pl,pl,1))) {
- /* p[pla-pl]...p[pla-1] are all zero */
- mpn_sub_1(p+pla-pl-1,p+pla-pl-1,pl+1,1);
- mpn_sub_1(p+pla-1,p+pla-1,1,1);
- }
- }
- else if (sqr==1) {
- if (pla>=2*pl)
- while ((sqr=mpn_add_1(p+pla-2*pl,p+pla-2*pl,2*pl,sqr)));
- else {
- sqr = mpn_sub_1(p+pla-pl,p+pla-pl,pl,sqr);
- ASSERT (sqr == 0);
- }
- }
- else
- ASSERT (sqr == 0);
-
- /* here p < 2^(2M) [K 2^(M(K-1)) + (K-1) 2^(M(K-2)) + ... ]
- < K 2^(2M) [2^(M(K-1)) + 2^(M(K-2)) + ... ]
- < K 2^(2M) 2^(M(K-1))*2 = 2^(M*K+M+k+1) */
- mpn_fft_norm_modF(op,p,pl,pla);
-}
-
-
-/* op <- n*m mod 2^N+1 with fft of size 2^k where N=pl*BITS_PER_MP_LIMB
- n and m have respectively nl and ml limbs
- op must have space for pl+1 limbs
- One must have pl = mpn_fft_next_size(pl, k).
-*/
-
-void
-#if __STDC__
-mpn_mul_fft (mp_ptr op, mp_size_t pl,
- mp_srcptr n, mp_size_t nl,
- mp_srcptr m, mp_size_t ml,
- int k)
-#else
-mpn_mul_fft (op, pl, n, nl, m, ml, k)
- mp_ptr op;
- mp_size_t pl;
- mp_srcptr n;
- mp_size_t nl;
- mp_srcptr m;
- mp_size_t ml;
- int k;
-#endif
-{
- int K,maxLK,i,j;
- mp_size_t N,Nprime,nprime,M,Mp,l;
- mp_limb_t **Ap,**Bp,*A,*T,*B;
- int **_fft_l;
- int sqr = (n==m && nl==ml);
- TMP_DECL(marker);
-
- TRACE (printf ("\nmpn_mul_fft pl=%ld nl=%ld ml=%ld k=%d\n",
- pl, nl, ml, k));
- ASSERT_ALWAYS (mpn_fft_next_size(pl, k) == pl);
-
- TMP_MARK(marker);
- N = pl*BITS_PER_MP_LIMB;
- _fft_l = TMP_ALLOC_TYPE (k+1, int*);
- for (i=0;i<=k;i++)
- _fft_l[i] = TMP_ALLOC_TYPE (1<<i, int);
- mpn_fft_initl(_fft_l, k);
- K = 1<<k;
- M = N/K; /* N = 2^k M */
- l = M/BITS_PER_MP_LIMB;
- maxLK = (K>BITS_PER_MP_LIMB) ? K : BITS_PER_MP_LIMB;
-
- Nprime = ((2*M+k+2+maxLK)/maxLK)*maxLK; /* ceil((2*M+k+3)/maxLK)*maxLK; */
- nprime = Nprime/BITS_PER_MP_LIMB;
- TRACE (printf ("N=%d K=%d, M=%d, l=%d, maxLK=%d, Np=%d, np=%d\n",
- N, K, M, l, maxLK, Nprime, nprime));
- if (nprime >= (sqr ? FFT_MODF_SQR_THRESHOLD : FFT_MODF_MUL_THRESHOLD)) {
- maxLK = (1<<mpn_fft_best_k(nprime,n==m))*BITS_PER_MP_LIMB;
- if (Nprime % maxLK) {
- Nprime=((Nprime/maxLK)+1)*maxLK;
- nprime = Nprime/BITS_PER_MP_LIMB;
- }
- TRACE (printf ("new maxLK=%d, Np=%d, np=%d\n", maxLK, Nprime, nprime));
- }
-
- T = TMP_ALLOC_LIMBS (nprime+1);
- Mp = Nprime/K;
-
- TRACE (printf("%dx%d limbs -> %d times %dx%d limbs (%1.2f)\n",
- pl,pl,K,nprime,nprime,2.0*(double)N/Nprime/K);
- printf(" temp space %ld\n", 2*K*(nprime+1)));
-
- A = _MP_ALLOCATE_FUNC_LIMBS (2*K*(nprime+1));
- B = A+K*(nprime+1);
- Ap = TMP_ALLOC_MP_PTRS (K);
- Bp = TMP_ALLOC_MP_PTRS (K);
- /* special decomposition for main call */
- for (i=0;i<K;i++) {
- Ap[i] = A+i*(nprime+1); Bp[i] = B+i*(nprime+1);
- /* store the next M bits of n into A[i] */
- /* supposes that M is a multiple of BITS_PER_MP_LIMB */
- if (nl>0) {
- j = (nl>=l) ? l : nl; /* limbs to store in Ap[i] */
- MPN_COPY(Ap[i], n, j); n+=l; MPN_ZERO(Ap[i]+j, nprime+1-j);
- mpn_fft_mul_2exp_modF(Ap[i], i*Mp, nprime, T);
- }
- else MPN_ZERO(Ap[i], nprime+1);
- nl -= l;
- if (n!=m) {
- if (ml>0) {
- j = (ml>=l) ? l : ml; /* limbs to store in Bp[i] */
- MPN_COPY(Bp[i], m, j); m+=l; MPN_ZERO(Bp[i]+j, nprime+1-j);
- mpn_fft_mul_2exp_modF(Bp[i], i*Mp, nprime, T);
- }
- else MPN_ZERO(Bp[i], nprime+1);
- }
- ml -= l;
- }
- mpn_mul_fft_internal(op,n,m,pl,k,K,Ap,Bp,A,B,nprime,l,Mp,_fft_l,T,0);
- TMP_FREE(marker);
- _MP_FREE_FUNC_LIMBS (A, 2*K*(nprime+1));
-}
-
-
-#if WANT_ASSERT
-static int
-#if __STDC__
-mpn_zero_p (mp_ptr p, mp_size_t n)
-#else
- mpn_zero_p (p, n)
- mp_ptr p;
- mp_size_t n;
-#endif
-{
- mp_size_t i;
-
- for (i = 0; i < n; i++)
- {
- if (p[i] != 0)
- return 0;
- }
-
- return 1;
-}
-#endif
-
-
-/* Multiply {n,nl}*{m,ml} and write the result to {op,nl+ml}.
-
- FIXME: Duplicating the result like this is wasteful, do something better
- perhaps at the norm_modF stage above. */
-
-void
-#if __STDC__
-mpn_mul_fft_full (mp_ptr op,
- mp_srcptr n, mp_size_t nl,
- mp_srcptr m, mp_size_t ml)
-#else
-mpn_mul_fft_full (op, n, nl, m, ml)
- mp_ptr op;
- mp_srcptr n;
- mp_size_t nl;
- mp_srcptr m;
- mp_size_t ml;
-#endif
-{
- mp_ptr pad_op;
- mp_size_t pl;
- int k;
- int sqr = (n==m && nl==ml);
-
- k = mpn_fft_best_k (nl+ml, sqr);
- pl = mpn_fft_next_size (nl+ml, k);
-
- TRACE (printf ("mpn_mul_fft_full nl=%ld ml=%ld -> pl=%ld k=%d\n",
- nl, ml, pl, k));
-
- pad_op = _MP_ALLOCATE_FUNC_LIMBS (pl+1);
- mpn_mul_fft (pad_op, pl, n, nl, m, ml, k);
-
- ASSERT (mpn_zero_p (pad_op+nl+ml, pl+1-(nl+ml)));
- MPN_COPY (op, pad_op, nl+ml);
-
- _MP_FREE_FUNC_LIMBS (pad_op, pl+1);
-}
projects / ghc-hetmet.git / blobdiff