From b50f939a743376620ddd875a10586228f25fa07f Mon Sep 17 00:00:00 2001
From: wbhart <goodwillhart@googlemail.com>
Date: Tue, 17 Jul 2012 21:33:24 +0000
Subject: [PATCH] Added missing mpn_mul_fft function (now
 mpn_mulmod_Bexpp1_fft) and used it.

---
 fft/mulmod_2expp1.c        | 31 +++++++++++++++++++++++++++++++
 gmp-impl.h                 |  5 +++++
 mpn/generic/inv_div_qr_n.c | 13 +++++++------
 mpn/generic/invert.c       |  7 ++++---
 4 files changed, 47 insertions(+), 9 deletions(-)

diff --git a/fft/mulmod_2expp1.c b/fft/mulmod_2expp1.c
index 54789334..1e0aa2fc 100644
--- a/fft/mulmod_2expp1.c
+++ b/fft/mulmod_2expp1.c
@@ -210,3 +210,34 @@ mpir_si fft_adjust_limbs(mp_size_t limbs)
 
    return limbs;
 }
+
+int
+mpn_mulmod_Bexpp1_fft (mp_ptr op, mp_size_t pl,
+	     mp_srcptr n, mp_size_t nl,
+	     mp_srcptr m, mp_size_t ml)
+{
+   mp_ptr a, b, tt;
+   mp_limb_t cy;
+
+   TMP_DECL;
+   TMP_MARK;
+
+   /* temporary space */
+   tt = TMP_ALLOC_LIMBS(2*pl);
+   
+   /* make copies of inputs, padded out to pl limbs */
+   a = TMP_ALLOC_LIMBS(pl + 1);
+   mpn_copyi(a, n, nl);
+   MPN_ZERO(a + nl, pl + 1 - nl);
+   
+   b = TMP_ALLOC_LIMBS(pl + 1);
+   mpn_copyi(b, m, ml);
+   MPN_ZERO(b + ml, pl + 1 - ml);
+   
+   /* this function only cares about the product, limbs = pl*GMP_LIMB_BITS */
+   cy = mpn_mulmod_2expp1(op, a, b, pl, GMP_LIMB_BITS, tt);
+   
+   TMP_FREE;
+
+   return cy;
+}
\ No newline at end of file
diff --git a/gmp-impl.h b/gmp-impl.h
index 7115e534..5d1a3f76 100644
--- a/gmp-impl.h
+++ b/gmp-impl.h
@@ -1421,6 +1421,11 @@ int mpn_mul_fft _PROTO ((mp_ptr op, mp_size_t pl,
                           mp_srcptr m, mp_size_t ml,
                           int k));
 
+#define mpn_mulmod_Bexpp1_fft __MPN(mulmod_Bexpp1_fft)
+int mpn_mulmod_Bexpp1_fft _PROTO ((mp_ptr op, mp_size_t pl,
+	     mp_srcptr n, mp_size_t nl,
+	     mp_srcptr m, mp_size_t ml));
+
 #define mpn_mul_fft_full  __MPN(mul_fft_full)
 void mpn_mul_fft_full _PROTO ((mp_ptr op,
                                mp_srcptr n, mp_size_t nl,
diff --git a/mpn/generic/inv_div_qr_n.c b/mpn/generic/inv_div_qr_n.c
index 030f155f..052b7a2a 100644
--- a/mpn/generic/inv_div_qr_n.c
+++ b/mpn/generic/inv_div_qr_n.c
@@ -22,6 +22,7 @@ MA 02110-1301, USA. */
 #include <mpir.h>
 #include "gmp-impl.h"
 #include "longlong.h"
+#include "fft/fft_tuning.h"
 
 /* 
    Computes the quotient and remainder of { np, 2*dn } by { dp, dn }.
@@ -34,7 +35,7 @@ mpn_inv_div_qr_n(mp_ptr qp, mp_ptr np,
 {
    mp_limb_t cy, lo, ret = 0, ret2 = 0;
    mp_size_t m, i;
-   mp_ptr tp, tp2;
+   mp_ptr tp;
    TMP_DECL;
 
    TMP_MARK;
@@ -82,11 +83,11 @@ mpn_inv_div_qr_n(mp_ptr qp, mp_ptr np,
    } else
    {
       mp_limb_t cy, cy2;
-      mp_size_t k;
-      k = mpn_fft_best_k (m, 0);
-      m = mpn_fft_next_size (m, k);
-
-      cy = mpn_mul_fft (tp, m, qp, dn, dp, dn, k);
+      
+      if (m >= FFT_MULMOD_2EXPP1_CUTOFF)
+         m = fft_adjust_limbs (m);
+      cy = mpn_mulmod_Bexpp1_fft (tp, m, qp, dn, dp, dn);
+      
       /* cy, {tp, m} = qp * dp mod (B^m+1) */ 
       cy2 = mpn_add_n(tp, tp, np + m, 2*dn - m);
       mpn_add_1(tp + 2*dn - m, tp + 2*dn - m, 2*m - 2*dn, cy2);
diff --git a/mpn/generic/invert.c b/mpn/generic/invert.c
index 1dc9606f..46ab0957 100644
--- a/mpn/generic/invert.c
+++ b/mpn/generic/invert.c
@@ -28,6 +28,7 @@ MA 02110-1301, USA. */
 #include "mpir.h"
 #include "gmp-impl.h"
 #include "longlong.h"
+#include "fft/fft_tuning.h"
 
 #define ZERO (mp_limb_t) 0
 #define ONE  (mp_limb_t) 1
@@ -151,11 +152,11 @@ mpn_invert (mp_ptr xp, mp_srcptr ap, mp_size_t n)
           mpir_ui k;
           int cc;
 
-          k = mpn_fft_best_k (m, 0);
-          m = mpn_fft_next_size (m, k);
+          if (m >= FFT_MULMOD_2EXPP1_CUTOFF)
+             m = fft_adjust_limbs (m);
           /* we have m >= n + 1 by construction, thus m > h */
           ASSERT(m < n + h);
-          cy = mpn_mul_fft (tp, m, ap, n, xp + l, h, k);
+          cy = mpn_mulmod_Bexpp1_fft (tp, m, ap, n, xp + l, h);
           /* cy, {tp, m} = A * {xp + l, h} mod (B^m+1) */
           cy += mpn_add_n (tp + h, tp + h, ap, m - h);
           cc = mpn_sub_n (tp, tp, ap + m - h, n + h - m);