Mentions légales du service

Skip to content
Snippets Groups Projects
Commit b7289136 authored by Andreas Enge's avatar Andreas Enge
Browse files

Use flint-3 for the subproduct trees in the trial division step.

* lib/cm-impl.h (cm_flint_tree_gcd): Declare new function.
* lib/flint.c (cm_flint_tree_gcd): New function.
* lib/nt.c (tree_gcd): Drop function when flint-3 is available.
  (cm_nt_mpz_tree_gcd): Use new function when flint-3 is available.
parent 1b4c9faa
No related branches found
No related tags found
No related merge requests found
......@@ -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);
......
......@@ -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. */
......
......@@ -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
}
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment