diff --git a/lib/cm-impl.h b/lib/cm-impl.h index 78399645d96f64631a0a453406dc3ec2d33d8b61..c1da1594ed5ef6e4fe8c281d384c3ad7e5636e32 100644 --- a/lib/cm-impl.h +++ b/lib/cm-impl.h @@ -132,6 +132,8 @@ extern char* cm_pari_sprintf_hfactor (int_cl_t d); #ifdef HAVE_FLINT3 extern void cm_flint_mpz_powm (mpz_ptr z, mpz_srcptr a, mpz_srcptr e, mpz_srcptr p); +extern void cm_flint_tree_gcd (mpz_t *gcd, mpz_srcptr n, mpz_t *m, + int no_m); #endif extern void cm_flint_mpzx_xplusa_pow_modmod (mpzx_ptr g, unsigned long int a, mpz_srcptr e, mpzx_srcptr m, mpz_srcptr p); diff --git a/lib/flint.c b/lib/flint.c index 738065ccd2d4e80d9acf6f043226933c837cd283..548483d182f0555a3c43064a38fa5b0bb9bb071a 100644 --- a/lib/flint.c +++ b/lib/flint.c @@ -146,6 +146,82 @@ void cm_flint_mpz_powm (mpz_ptr z, mpz_srcptr a, mpz_srcptr e, mpz_srcptr p) } #endif +/*****************************************************************************/ +/* */ +/* Function for product tree relying on FLINT. */ +/* */ +/*****************************************************************************/ + +#ifdef HAVE_FLINT3 +void cm_flint_tree_gcd (mpz_t *gcd, mpz_srcptr n, mpz_t *m, int no_m) + /* This is a copy-paste of the static function tree_gcd from nt.c + with fmpz types. */ +{ + fmpz_t np; + mpz_t tmp; + fmpz_t **tree; + int *width; + int levels; + int i, j; + + fmpz_init (np); + fmpz_set_mpz (np, n); + mpz_init (tmp); + + /* Compute the height of the subproduct tree of the m in levels. */ + for (i = no_m, levels = 1; i > 1; i = (i+1) / 2, levels++); + + /* Compute bottom-up a subproduct tree with m on the leaves. */ + tree = (fmpz_t **) malloc (levels * sizeof (fmpz_t *)); + width = (int *) malloc (levels * sizeof (int)); + width [0] = no_m; + tree [0] = (fmpz_t *) malloc (no_m * sizeof (fmpz_t)); + for (j = 0; j < no_m; j++) { + fmpz_init (tree [0][j]); + fmpz_set_mpz (tree [0][j], m [j]); + } + for (i = 1; i < levels; i++) { + width [i] = (width [i-1] + 1) / 2; + tree [i] = (fmpz_t *) malloc (width [i] * sizeof (fmpz_t)); + for (j = 0; j < width [i-1] / 2; j++) { + fmpz_init (tree [i][j]); + fmpz_mul (tree [i][j], tree [i-1][2*j], tree [i-1][2*j+1]); + } + if (width [i-1] % 2 != 0) + fmpz_init_set (tree [i][j], tree [i-1][2*j]); + } + + /* Replace the tree tops by n modulo the entry. */ + for (j = 0; j < width [levels-1]; j++) + fmpz_mod (tree [levels-1][j], np, tree [levels-1][j]); + + /* Replace top-down the tree entries by n modulo the entry. */ + for (i = levels - 2; i >= 0; i--) { + for (j = 0; j < (width [i] / 2) * 2; j++) + fmpz_mod (tree [i][j], tree [i+1][j/2], tree [i][j]); + if (width [i] % 2 != 0) + fmpz_set (tree [i][j], tree [i+1][j/2]); + } + + /* Compute the gcd of n mod m [j] and m [j]. */ + for (j = 0; j < no_m; j++) { + fmpz_get_mpz (tmp, tree [0][j]); + mpz_gcd (gcd [j], tmp, m [j]); + } + + /* Clear the tree. */ + for (i = 0; i < levels; i++) { + for (j = 0; j < width [i]; j++) + fmpz_clear (tree [i][j]); + free (tree [i]); + } + free (tree); + free (width); + fmpz_clear (np); + mpz_clear (tmp); +} +#endif + /*****************************************************************************/ /* */ /* Functions for mpzx modulo p relying on FLINT. */ diff --git a/lib/nt.c b/lib/nt.c index a2a8c8abcd7dbb59404491c460a8325f090abe4a..e76a129de14ea3b5de564b50588a436f22dc0259 100644 --- a/lib/nt.c +++ b/lib/nt.c @@ -24,7 +24,9 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. #include "cm-impl.h" static void cm_mpz_powm (mpz_ptr z, mpz_srcptr a, mpz_srcptr e, mpz_srcptr p); +#ifndef HAVE_FLINT3 static void tree_gcd (mpz_t *gcd, mpz_srcptr n, mpz_t *m, int no_m); +#endif static int miller_rabin (mpz_srcptr n); static void cm_nt_mpz_tonelli_with_generator (mpz_ptr root, mpz_srcptr a, mpz_srcptr p, unsigned int e, mpz_srcptr q, mpz_srcptr z); @@ -455,6 +457,7 @@ void cm_nt_mpz_tonelli_si (mpz_ptr root, const long int a, mpz_srcptr p) /*****************************************************************************/ +#ifndef HAVE_FLINT3 static void tree_gcd (mpz_t *gcd, mpz_srcptr n, mpz_t *m, int no_m) /* This is the workhorse behind the cm_nt_mpz_tree_gcd function; it assumes that the total size of the no_m entries of m is less than @@ -513,6 +516,7 @@ static void tree_gcd (mpz_t *gcd, mpz_srcptr n, mpz_t *m, int no_m) free (tree); free (width); } +#endif /*****************************************************************************/ @@ -538,7 +542,11 @@ void cm_nt_mpz_tree_gcd (mpz_t *gcd, mpz_srcptr n, mpz_t *m, int no_m) && size_m + mpz_sizeinbase (m [offset + no_batch], 2) < size_n / 2; size_m += mpz_sizeinbase (m [offset + no_batch], 2), no_batch++); +#ifdef HAVE_FLINT3 + cm_flint_tree_gcd (gcd + offset, n, m + offset, no_batch); +#else tree_gcd (gcd + offset, n, m + offset, no_batch); +#endif } }