LCOV - code coverage report
Current view: top level - ecm - pm1.c (source / functions) Hit Total Coverage
Test: unnamed Lines: 243 261 93.1 %
Date: 2022-03-21 11:19:20 Functions: 8 8 100.0 %

          Line data    Source code
       1             : /* Pollard 'P-1' algorithm.
       2             : 
       3             : Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011,
       4             : 2012, Paul Zimmermann and Alexander Kruppa.
       5             : 
       6             : This file is part of the ECM Library.
       7             : 
       8             : The ECM Library is free software; you can redistribute it and/or modify
       9             : it under the terms of the GNU Lesser General Public License as published by
      10             : the Free Software Foundation; either version 3 of the License, or (at your
      11             : option) any later version.
      12             : 
      13             : The ECM Library is distributed in the hope that it will be useful, but
      14             : WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
      15             : or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
      16             : License for more details.
      17             : 
      18             : You should have received a copy of the GNU Lesser General Public License
      19             : along with the ECM Library; see the file COPYING.LIB.  If not, see
      20             : http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
      21             : 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
      22             : 
      23             : #include <math.h>
      24             : #include <stdlib.h>
      25             : #include "ecm-impl.h"
      26             : #include "getprime_r.h"
      27             : 
      28             : #define CASCADE_THRES 3
      29             : #define CASCADE_MAX 50000000.0
      30             : 
      31             : typedef struct {
      32             :   unsigned int size;
      33             :   mpz_t *val;
      34             : } mul_casc;
      35             : 
      36             : /******************************************************************************
      37             : *                                                                             *
      38             : *                                  Stage 1                                    *
      39             : *                                                                             *
      40             : ******************************************************************************/
      41             : 
      42             : /* prime powers are accumulated up to about n^L1 */
      43             : #define L1 16
      44             : 
      45             : /*** Cascaded multiply ***/
      46             : 
      47             : /* return NULL if an error occurred */
      48             : static mul_casc *
      49         267 : mulcascade_init (void)
      50             : {
      51             :   mul_casc *t;
      52             : 
      53         267 :   t = (mul_casc *) malloc (sizeof (mul_casc));
      54         267 :   ASSERT_ALWAYS(t != NULL);
      55         267 :   t->val = (mpz_t*) malloc (sizeof (mpz_t));
      56         267 :   ASSERT_ALWAYS(t->val != NULL);
      57         267 :   mpz_init (t->val[0]);
      58         267 :   t->size = 1;
      59         267 :   return t;
      60             : }
      61             : 
      62             : static void 
      63         267 : mulcascade_free (mul_casc *c)
      64             : {
      65             :   unsigned int i;
      66             : 
      67        2183 :   for (i = 0; i < c->size; i++)
      68        1916 :     mpz_clear (c->val[i]);
      69         267 :   free (c->val);
      70         267 :   free (c);
      71         267 : }
      72             : 
      73             : static void
      74     1996877 : mulcascade_mul_d (mul_casc *c, const double n, ATTRIBUTE_UNUSED mpz_t t)
      75             : {
      76             :   unsigned int i;
      77             : 
      78     1996877 :   if (mpz_sgn (c->val[0]) == 0)
      79             :     {
      80      164067 :       mpz_set_d (c->val[0], n);
      81      164067 :       return;
      82             :     }
      83             : 
      84     1832810 :   mpz_mul_d (c->val[0], c->val[0], n, t);
      85     1832810 :   if (mpz_size (c->val[0]) <= CASCADE_THRES)
      86     1668984 :     return;
      87             :   
      88      326742 :   for (i = 1; i < c->size; i++) 
      89             :     {
      90      325093 :       if (mpz_sgn (c->val[i]) == 0) 
      91             :         {
      92      162177 :           mpz_set (c->val[i], c->val[i-1]);
      93      162177 :           mpz_set_ui (c->val[i-1], 0);
      94      162177 :           return;
      95             :         }
      96             :       else
      97             :         {
      98      162916 :           mpz_mul (c->val[i], c->val[i], c->val[i-1]);
      99      162916 :           mpz_set_ui (c->val[i-1], 0);
     100             :         }
     101             :     }
     102             :   
     103             :   /* Allocate more space for cascade */
     104             :   
     105        1649 :   i = c->size++;
     106        1649 :   c->val = (mpz_t*) realloc (c->val, c->size * sizeof (mpz_t));
     107        1649 :   ASSERT_ALWAYS(c->val != NULL);
     108        1649 :   mpz_init (c->val[i]);
     109        1649 :   mpz_swap (c->val[i], c->val[i-1]);
     110             : }
     111             : 
     112             : /* initialize c to n (assumes c is empty) */
     113             : static void
     114          21 : mulcascade_set (mul_casc *c, mpz_t n)
     115             : {
     116             :   ASSERT(mpz_sgn (c->val[0]) == 0);
     117             : 
     118          21 :   mpz_set (c->val[0], n);
     119          21 : }
     120             : 
     121             : static void 
     122         267 : mulcascade_get_z (mpz_t r, mul_casc *c) 
     123             : {
     124             :   unsigned int i;
     125             : 
     126             :   ASSERT(c->size != 0);
     127             : 
     128         267 :   mpz_set_ui (r, 1);
     129             :   
     130        2183 :   for (i = 0; i < c->size; i++)
     131        1916 :     if (mpz_sgn (c->val[i]) != 0)
     132        1172 :       mpz_mul (r, r, c->val[i]);
     133         267 : }
     134             : 
     135             : 
     136             : /* Input:  a is the generator (sigma)
     137             :            n is the number to factor
     138             :            B1 is the stage 1 bound
     139             :            B1done: stage 1 was already done up to that limit
     140             :            go is the group order to preload
     141             :    Output: f is the factor found, a is the value at end of stage 1
     142             :            B1done is set to B1 if stage 1 completed normally,
     143             :            or to the largest prime processed if interrupted, but never
     144             :            to a smaller value than B1done was upon function entry.
     145             :    Return value: non-zero iff a factor was found (or an error occurred).
     146             : */
     147             : 
     148             : static int
     149         267 : pm1_stage1 (mpz_t f, mpres_t a, mpmod_t n, double B1, double *B1done, 
     150             :             mpz_t go, int (*stop_asap)(void), char *chkfilename)
     151             : {
     152             :   double p, q, r, cascade_limit, last_chkpnt_p;
     153             :   mpz_t g, d;
     154         267 :   int youpi = ECM_NO_FACTOR_FOUND;
     155             :   unsigned int size_n, max_size;
     156         267 :   unsigned int smallbase = 0;
     157             :   mul_casc *cascade;
     158             :   long last_chkpnt_time;
     159         267 :   const double B0 = sqrt (B1);
     160             :   prime_info_t prime_info;
     161             : 
     162         267 :   mpz_init (g);
     163         267 :   mpz_init (d);
     164             : 
     165         267 :   size_n = mpz_sizeinbase (n->orig_modulus, 2);
     166         267 :   max_size = L1 * size_n;
     167             : 
     168         267 :   mpres_get_z (g, a, n);
     169         267 :   if (mpz_fits_uint_p (g))
     170         247 :     smallbase = mpz_get_ui (g);
     171             : 
     172         267 :   mpz_set_ui (g, 1);
     173             : 
     174             :   /* Set a limit of roughly 10000 * log_10(N) for the primes that are 
     175             :      multiplied up in the exponent, i.e. 1M for a 100 digit number, 
     176             :      but limit to CASCADE_MAX to avoid problems with stack allocation */
     177             :   
     178         267 :   cascade_limit = 3000.0 * (double) size_n;
     179             : 
     180         267 :   if (cascade_limit > CASCADE_MAX)
     181           1 :     cascade_limit = CASCADE_MAX;
     182             :   
     183         267 :   if (cascade_limit > B1)
     184         216 :     cascade_limit = B1;
     185             : 
     186         267 :   cascade = mulcascade_init ();
     187         267 :   if (cascade == NULL)
     188             :     {
     189           0 :       youpi = ECM_ERROR;
     190           0 :       goto clear_pm1_stage1;
     191             :     }
     192             : 
     193             :   /* since B0 = sqrt(B1), we can have B0 > cascade_limit only when
     194             :      B1 > cascade_limit^2. This cannot happen when cascade_limit=B1,
     195             :      thus we need B1 > min(CASCADE_MAX, 3000*sizeinbase(n,2))^2.
     196             :      For sizeinbase(n,2) <= CASCADE_MAX/3000 (less than 5017 digits 
     197             :      for CASCADE_MAX=5e7) this means B1 > 9e6*sizeinbase(n,2)^2.
     198             :      For sizeinbase(n,2) > CASCADE_MAX/3000, this means B1 > CASCADE_MAX^2,
     199             :      i.e. B1 > 25e14 for CASCADE_MAX=5e7.
     200             :   */
     201             : 
     202             :   /* if the user knows that P-1 has a given divisor, he can supply it */
     203         267 :   if (mpz_cmp_ui (go, 1) > 0)
     204          21 :     mulcascade_set (cascade, go);
     205             :   
     206         267 :   last_chkpnt_time = cputime ();
     207         267 :   last_chkpnt_p = 2.;
     208             :   
     209             :   /* Fill the multiplication cascade with the product of small stage 1 
     210             :      primes */
     211             :   /* Add small primes <= MIN(sqrt(B1), cascade_limit) in the appropriate 
     212             :      power to the cascade */
     213         267 :   prime_info_init (prime_info);
     214       19332 :   for (p = 2.; p <= MIN(B0, cascade_limit); p = (double) getprime_mt (prime_info))
     215             :     {
     216       66097 :       for (q = 1., r = p; r <= B1; r *= p)
     217       47032 :         if (r > *B1done) q *= p;
     218       19065 :       mulcascade_mul_d (cascade, q, d);
     219             :     }
     220             : 
     221             :   /* If B0 < cascade_limit, we can add some primes > sqrt(B1) with 
     222             :      exponent 1 to the cascade */
     223     2120924 :   for ( ; p <= cascade_limit; p = (double) getprime_mt (prime_info))
     224     2120657 :     if (p > *B1done)
     225     1977812 :       mulcascade_mul_d (cascade, p, d);
     226             : 
     227             :   /* Now p > cascade_limit, flush cascade and exponentiate */
     228         267 :   mulcascade_get_z (g, cascade);
     229         267 :   mulcascade_free (cascade);
     230         267 :   outputf (OUTPUT_DEVVERBOSE, "Exponent has %u bits\n", 
     231             :            mpz_sizeinbase (g, 2));
     232             :   
     233         267 :   if (smallbase)
     234             :     {
     235         247 :       outputf (OUTPUT_DEVVERBOSE, "Using mpres_ui_pow, base %u\n", smallbase);
     236         247 :       mpres_ui_pow (a, smallbase, g, n);
     237             :     }
     238             :   else
     239             :     {
     240          20 :       mpres_pow (a, a, g, n);
     241             :     }
     242         267 :   mpz_set_ui (g, 1);
     243             : 
     244             :   /* If B0 > cascade_limit, we need to process the primes 
     245             :      cascade_limit < p < B0 in the appropriate exponent yet */
     246         273 :   for ( ; p <= B0; p = (double) getprime_mt (prime_info))
     247             :     {
     248          18 :       for (q = 1, r = p; r <= B1; r *= p)
     249          12 :         if (r > *B1done) q *= p;
     250           6 :       mpz_mul_d (g, g, q, d);
     251           6 :       if (mpz_sizeinbase (g, 2) >= max_size)
     252             :         {
     253           1 :           mpres_pow (a, a, g, n);
     254           1 :           mpz_set_ui (g, 1);
     255           1 :         if (stop_asap != NULL && (*stop_asap) ())
     256             :           {
     257           0 :             outputf (OUTPUT_NORMAL, "Interrupted at prime %.0f\n", p);
     258           0 :             if (p > *B1done)
     259           0 :               *B1done = p;
     260           0 :             goto clear_pm1_stage1;
     261             :           }
     262             :         }
     263             :     }
     264             : 
     265             :   /* All primes sqrt(B1) < p <= B1 appear with exponent 1. All primes <= B1done
     266             :      are already included with exponent at least 1, so it's safe to skip
     267             :      ahead to B1done+1. */
     268             : 
     269         267 :   while (p <= *B1done)
     270           0 :     p = (double) getprime_mt (prime_info);
     271             : 
     272             :   /* then remaining primes > max(sqrt(B1), cascade_limit) and taken 
     273             :      with exponent 1 */
     274    30123131 :   for (; p <= B1; p = (double) getprime_mt (prime_info))
     275             :   {
     276    30122864 :     mpz_mul_d (g, g, p, d);
     277    30122864 :     if (mpz_sizeinbase (g, 2) >= max_size)
     278             :       {
     279     5340447 :         mpres_pow (a, a, g, n);
     280     5340447 :         mpz_set_ui (g, 1);
     281     5340447 :         if (stop_asap != NULL && (*stop_asap) ())
     282             :           {
     283           0 :             outputf (OUTPUT_NORMAL, "Interrupted at prime %.0f\n", p);
     284           0 :              if (p > *B1done)
     285           0 :               *B1done = p;
     286           0 :             goto clear_pm1_stage1;
     287             :           }
     288     5340447 :         if (chkfilename != NULL && p > last_chkpnt_p + 10000. &&
     289           0 :             elltime (last_chkpnt_time, cputime ()) > CHKPNT_PERIOD)
     290             :           {
     291           0 :             writechkfile (chkfilename, ECM_PM1, p, n, NULL, a, NULL, NULL);
     292           0 :             last_chkpnt_p = p;
     293           0 :             last_chkpnt_time = cputime ();
     294             :           }
     295             :       }
     296             :   }
     297             : 
     298         267 :   mpres_pow (a, a, g, n);
     299             :   
     300             :   /* If stage 1 finished normally, p is the smallest prime >B1 here.
     301             :      In that case, set to B1 */
     302         267 :   if (p > B1)
     303         267 :       p = B1;
     304             :   
     305         267 :   if (p > *B1done)
     306         262 :     *B1done = p;
     307             :   
     308         267 :   mpres_sub_ui (a, a, 1, n);
     309         267 :   mpres_gcd (f, a, n);
     310         267 :   if (mpz_cmp_ui (f, 1) > 0)
     311          21 :     youpi = ECM_FACTOR_FOUND_STEP1;
     312         267 :   mpres_add_ui (a, a, 1, n);
     313             : 
     314         267 :  clear_pm1_stage1:
     315         267 :   if (chkfilename != NULL)
     316           5 :     writechkfile (chkfilename, ECM_PM1, *B1done, n, NULL, a, NULL, NULL);
     317         267 :   prime_info_clear (prime_info); /* free the prime table */
     318         267 :   mpz_clear (d);
     319         267 :   mpz_clear (g);
     320             : 
     321         267 :   return youpi;
     322             : }
     323             : 
     324             : 
     325             : void
     326          65 : print_prob (double B1, const mpz_t B2, unsigned long dF, unsigned long k, 
     327             :             int S, const mpz_t go)
     328             : {
     329             :   double prob;
     330             :   int i;
     331             :   char sep;
     332             : 
     333          65 :   outputf (OUTPUT_VERBOSE, "Probability of finding a factor of n digits (assuming one exists):\n");
     334          65 :   outputf (OUTPUT_VERBOSE, "20\t25\t30\t35\t40\t45\t50\t55\t60\t65\n");
     335         715 :   for (i = 20; i <= 65; i += 5)
     336             :     {
     337         650 :       sep = (i < 65) ? '\t' : '\n';
     338         650 :       prob = pm1prob (B1, mpz_get_d (B2),
     339         650 :                       pow (10., i - .5), (double) dF * dF * k, S, go);
     340         650 :       outputf (OUTPUT_VERBOSE, "%.2g%c", prob, sep);
     341             :     }
     342          65 : }
     343             : 
     344             : 
     345             : 
     346             : /******************************************************************************
     347             : *                                                                             *
     348             : *                                Pollard P-1                                  *
     349             : *                                                                             *
     350             : ******************************************************************************/
     351             : 
     352             : /* Input: p is the initial generator (sigma), if 0, generate it at random.
     353             :           N is the number to factor
     354             :           B1 is the stage 1 bound
     355             :           B2 is the stage 2 bound
     356             :           B1done is the stage 1 limit to which supplied residue has 
     357             :             already been computed
     358             :           k is the number of blocks for stage 2
     359             :           verbose is the verbosity level
     360             :    Output: f is the factor found, p is the residue at end of stage 1
     361             :    Return value: non-zero iff a factor is found (1 for stage 1, 2 for stage 2)
     362             : */
     363             : int
     364         272 : pm1 (mpz_t f, mpz_t p, mpz_t N, mpz_t go, double *B1done, double B1,
     365             :      mpz_t B2min_parm, mpz_t B2_parm, unsigned long k, 
     366             :      int verbose, int repr, int use_ntt, FILE *os, FILE *es, 
     367             :      char *chkfilename, char *TreeFilename, double maxmem, 
     368             :      gmp_randstate_t rng, int (*stop_asap)(void))
     369             : {
     370         272 :   int youpi = ECM_NO_FACTOR_FOUND;
     371             :   long st;
     372             :   mpmod_t modulus;
     373             :   mpres_t x;
     374             :   mpz_t B2min, B2; /* Local B2, B2min to avoid changing caller's values */
     375             :   faststage2_param_t params;
     376             : 
     377         272 :   set_verbose (verbose);
     378         272 :   ECM_STDOUT = (os == NULL) ? stdout : os;
     379         272 :   ECM_STDERR = (es == NULL) ? stdout : es;
     380             : 
     381             :   /* if n is even, return 2 */
     382         272 :   if (mpz_divisible_2exp_p (N, 1))
     383             :     {
     384           0 :       mpz_set_ui (f, 2);
     385           0 :       return ECM_FACTOR_FOUND_STEP1;
     386             :     }
     387             : 
     388         272 :   st = cputime ();
     389             : 
     390         272 :   if (mpz_cmp_ui (p, 0) == 0)
     391         212 :     pm1_random_seed (p, N, rng);
     392             : 
     393         272 :   mpz_init_set (B2min, B2min_parm);
     394         272 :   mpz_init_set (B2, B2_parm);
     395             : 
     396             :   /* Set default B2. See ecm.c for comments */
     397         272 :   if (ECM_IS_DEFAULT_B2(B2))
     398          61 :     mpz_set_d (B2, pow (B1 * PM1FS2_COST, PM1FS2_DEFAULT_B2_EXPONENT));
     399             : 
     400             :   /* set B2min */
     401         272 :   if (mpz_sgn (B2min) < 0)
     402         232 :     mpz_set_d (B2min, B1);
     403             : 
     404             :   /* choice of modular arithmetic: if default choice, choose mpzmod which
     405             :      is always faster, since mpz_powm uses base-k sliding window exponentiation
     406             :      and mpres_pow does not */
     407         272 :   if (repr == ECM_MOD_DEFAULT && isbase2 (N, BASE2_THRESHOLD) == 0)
     408          95 :     mpmod_init (modulus, N, ECM_MOD_MPZ);
     409             :   else
     410         177 :     mpmod_init (modulus, N, repr);
     411             : 
     412             :   /* Determine parameters (polynomial degree etc.) */
     413             : 
     414             :     {
     415             :       long P_ntt, P_nontt;
     416         272 :       const unsigned long lmax = 1UL << 30; /* An upper bound */
     417             :       unsigned long lmax_NTT, lmax_noNTT;
     418             :       faststage2_param_t params_ntt, params_nontt, *better_params;
     419             :       mpz_t effB2min_ntt, effB2_ntt, effB2min_nontt, effB2_nontt;
     420             : 
     421         272 :       mpz_init (params.m_1);
     422         272 :       params.l = 0;
     423         272 :       mpz_init (params_ntt.m_1);
     424         272 :       params_ntt.l = 0;
     425         272 :       mpz_init (params_nontt.m_1);
     426         272 :       params_nontt.l = 0;
     427         272 :       mpz_init (effB2min_ntt);
     428         272 :       mpz_init (effB2_ntt);
     429         272 :       mpz_init (effB2min_nontt);
     430         272 :       mpz_init (effB2_nontt);
     431             : 
     432             :       /* Find out what the longest transform length is we can do at all.
     433             :          If no maxmem is given, the non-NTT can theoretically do any length. */
     434             : 
     435         272 :       lmax_NTT = 0;
     436         272 :       if (use_ntt)
     437             :         {
     438             :           unsigned long t;
     439             :           /* See what transform length the NTT can handle (due to limited 
     440             :              primes and limited memory) */
     441         210 :           t = mpzspm_max_len (N);
     442         210 :           lmax_NTT = MIN (lmax, t);
     443         210 :           if (maxmem != 0.)
     444             :             {
     445           5 :               t = pm1fs2_maxlen (double_to_size (maxmem), N, use_ntt);
     446           5 :               lmax_NTT = MIN (lmax_NTT, t);
     447             :             }
     448         210 :           outputf (OUTPUT_DEVVERBOSE, "NTT can handle lmax <= %lu\n", lmax_NTT);
     449         210 :           P_ntt = choose_P (B2min, B2, lmax_NTT, k, &params_ntt, 
     450             :                             effB2min_ntt, effB2_ntt, 1, ECM_PM1);
     451         210 :           if (P_ntt != ECM_ERROR)
     452         205 :             outputf (OUTPUT_DEVVERBOSE,
     453             :                      "Parameters for NTT: P=%lu, l=%lu\n", 
     454             :                      params_ntt.P, params_ntt.l);
     455             :         }
     456             :       else
     457          62 :         P_ntt = 0; /* or GCC complains about uninitialized var */
     458             :       
     459             :       /* See what transform length the non-NTT code can handle */
     460         272 :       lmax_noNTT = lmax;
     461         272 :       if (maxmem != 0.)
     462             :         {
     463             :           unsigned long t;
     464          10 :           t = pm1fs2_maxlen (double_to_size (maxmem), N, 0);
     465          10 :           lmax_noNTT = MIN (lmax_noNTT, t);
     466          10 :           outputf (OUTPUT_DEVVERBOSE, "non-NTT can handle lmax <= %lu\n", 
     467             :                    lmax_noNTT);
     468             :         }
     469         272 :       if (use_ntt != 2)
     470         257 :         P_nontt = choose_P (B2min, B2, lmax_noNTT, k, &params_nontt, 
     471             :                             effB2min_nontt, effB2_nontt, 0, ECM_PM1);
     472             :       else
     473          15 :         P_nontt = ECM_ERROR;
     474         272 :       if (P_nontt != ECM_ERROR)
     475         257 :         outputf (OUTPUT_DEVVERBOSE,
     476             :                  "Parameters for non-NTT: P=%lu, l=%lu\n", 
     477             :                  params_nontt.P, params_nontt.l);
     478             :       
     479         272 :       if (((!use_ntt || P_ntt == ECM_ERROR) && P_nontt == ECM_ERROR) ||
     480          10 :           (use_ntt == 2 && P_ntt == ECM_ERROR))
     481             :         {
     482           5 :           outputf (OUTPUT_ERROR, 
     483             :                    "Error: cannot choose suitable P value for your stage 2 "
     484             :                    "parameters.\nTry a shorter B2min,B2 interval.\n");
     485           5 :           mpz_clear (params.m_1);
     486           5 :           mpz_clear (params_ntt.m_1);
     487           5 :           mpz_clear (params_nontt.m_1);
     488           5 :           return ECM_ERROR;
     489             :         }
     490             : 
     491             :       /* Now decide whether to take NTT or non-NTT. Since the non-NTT code
     492             :          uses more memory, we only use it when -no-ntt was given, or when
     493             :          we can't find good parameters for the NTT code.
     494             :          Warning: we only do that when B2 >= B2min. */
     495         267 :       if (mpz_cmp (B2, B2min) >= 0)
     496             :         {
     497         262 :           if (use_ntt == 0 || P_ntt == ECM_ERROR)
     498             :             {
     499          61 :               better_params = &params_nontt;
     500          61 :               mpz_set (B2min, effB2min_nontt);
     501          61 :               mpz_set (B2, effB2_nontt);
     502          61 :               use_ntt = 0;
     503             :             }
     504             :           else
     505             :             {
     506         201 :               better_params = &params_ntt;
     507         201 :               mpz_set (B2min, effB2min_ntt);
     508         201 :               mpz_set (B2, effB2_ntt);
     509         201 :               use_ntt = 1;
     510             :             }
     511             : 
     512         262 :           params.P = better_params->P;
     513         262 :           params.s_1 = better_params->s_1;
     514         262 :           params.s_2 = better_params->s_2;
     515         262 :           params.l = better_params->l;
     516         262 :           mpz_set (params.m_1, better_params->m_1);
     517         262 :           params.file_stem = TreeFilename;
     518         262 :           params.file_stem = TreeFilename;
     519             :         }
     520             : 
     521         267 :       mpz_clear (params_ntt.m_1);
     522         267 :       mpz_clear (params_nontt.m_1);
     523         267 :       mpz_clear (effB2min_ntt);
     524         267 :       mpz_clear (effB2_ntt);
     525         267 :       mpz_clear (effB2min_nontt);
     526         267 :       mpz_clear (effB2_nontt);
     527             :       
     528         267 :       outputf (OUTPUT_VERBOSE, "Using lmax = %lu with%s NTT which takes "
     529             :                "about %luMB of memory\n", params.l, 
     530             :                (use_ntt) ? "" : "out", 
     531         267 :                pm1fs2_memory_use (params.l, N, use_ntt) / 1048576);
     532             :     }
     533             :   
     534             :   /* Print B1, B2, polynomial and x0 */
     535         267 :   print_B1_B2_poly (OUTPUT_NORMAL, ECM_PM1, B1, *B1done, B2min_parm, B2min, 
     536             :                     B2, 1, p, 0, 0, NULL, 0, 0);
     537             : 
     538             :   /* If we do a stage 2, print its parameters */
     539         267 :   if (mpz_cmp (B2, B2min) >= 0)
     540             :     {
     541             :       /* can't mix 64-bit types and mpz_t on win32 for some reason */
     542         262 :       outputf (OUTPUT_VERBOSE, "P = %" PRId64 ", l = %lu"
     543             :                     ", s_1 = %" PRId64 ", k = s_2 = %" PRId64 ,
     544             :              params.P, params.l,
     545             :              params.s_1,params.s_2);
     546         262 :       outputf (OUTPUT_VERBOSE, ", m_1 = %Zd\n", params.m_1);
     547             :     }
     548             : 
     549         267 :   if (test_verbose (OUTPUT_VERBOSE))
     550             :     {
     551          45 :       if (mpz_sgn (B2min_parm) >= 0)
     552             :         {
     553          10 :           outputf (OUTPUT_VERBOSE, 
     554             :             "Can't compute success probabilities for B1 <> B2min\n");
     555             :         }
     556             :       else
     557             :         {
     558          35 :           rhoinit (256, 10);
     559          35 :           print_prob (B1, B2, 0, k, 1, go);
     560             :         }
     561             :     }
     562             : 
     563         267 :   mpres_init (x, modulus);
     564         267 :   mpres_set_z (x, p, modulus);
     565             : 
     566         267 :   st = cputime ();
     567             : 
     568         267 :   if (B1 > *B1done || mpz_cmp_ui (go, 1) > 0)
     569         267 :     youpi = pm1_stage1 (f, x, modulus, B1, B1done, go, stop_asap, chkfilename);
     570             : 
     571         267 :   st = elltime (st, cputime ());
     572             : 
     573         267 :   outputf (OUTPUT_NORMAL, "Step 1 took %ldms\n", st);
     574         267 :   if (test_verbose (OUTPUT_RESVERBOSE))
     575             :     {
     576             :       mpz_t tx;
     577          30 :       mpz_init (tx);
     578          30 :       mpres_get_z (tx, x, modulus);
     579          30 :       outputf (OUTPUT_RESVERBOSE, "x=%Zd\n", tx);
     580          30 :       mpz_clear (tx);
     581             :     }
     582             : 
     583         267 :   if (stop_asap != NULL && (*stop_asap) ())
     584           0 :     goto clear_and_exit;
     585             : 
     586         267 :   if (youpi == ECM_NO_FACTOR_FOUND && mpz_cmp (B2, B2min) >= 0)
     587             :     {
     588         246 :       if (use_ntt)
     589         187 :         youpi = pm1fs2_ntt (f, x, modulus, &params);
     590             :       else
     591          59 :         youpi = pm1fs2 (f, x, modulus, &params);
     592             :     }
     593             : 
     594         267 :   if (test_verbose (OUTPUT_VERBOSE))
     595             :     {
     596          45 :       if (mpz_sgn (B2min_parm) < 0)
     597          35 :         rhoinit (1, 0); /* Free memory of rhotable */
     598             :     }
     599             : 
     600         232 : clear_and_exit:
     601         267 :   mpres_get_z (p, x, modulus);
     602         267 :   mpres_clear (x, modulus);
     603         267 :   mpmod_clear (modulus);
     604         267 :   mpz_clear (params.m_1);
     605         267 :   mpz_clear (B2);
     606         267 :   mpz_clear (B2min);
     607             : 
     608         267 :   return youpi;
     609             : }

Generated by: LCOV version 1.14