cheng/mpir
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Added missing mpn_mul_fft function (now mpn_mulmod_Bexpp1_fft) and used it.

Browse Source
This commit is contained in:
wbhart 2012-07-17 21:33:24 +00:00
parent 0e06c37a20
commit b50f939a74
4 changed files with 47 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

31
fft/mulmod_2expp1.c
View File

@ -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;
}

5
gmp-impl.h
View File

@ -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,

13
mpn/generic/inv_div_qr_n.c
View File

@ -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);

7
mpn/generic/invert.c
View File

@ -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);