LCOV - code coverage report
Current view: top level - ecm - stage2.c (source / functions) Hit Total Coverage
Test: unnamed Lines: 288 310 92.9 %
Date: 2022-03-21 11:19:20 Functions: 6 6 100.0 %

          Line data    Source code
       1             : /* Common stage 2 for ECM, P-1 and P+1 (improved standard continuation
       2             :    with subquadratic polynomial arithmetic).
       3             : 
       4             : Copyright 2001-2015 Paul Zimmermann, Alexander Kruppa, Dave Newman.
       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 "config.h"
      24             : #include <stdio.h>
      25             : #include <stdlib.h>
      26             : #include <math.h> /* for floor */
      27             : #include <string.h> /* for strlen */
      28             : 
      29             : #ifdef HAVE_UNISTD_H
      30             : #include <unistd.h> /* for unlink */
      31             : #endif
      32             : 
      33             : #include "ecm-impl.h"
      34             : #include "sp.h"
      35             : 
      36             : extern unsigned int Fermat;
      37             : 
      38             : /* r <- Dickson(n,a)(x) */
      39             : static void 
      40       15176 : dickson (mpz_t r, mpz_t x, unsigned int n, int a)
      41             : {
      42       15176 :   unsigned int i, b = 0;
      43             :   mpz_t t, u;
      44             : 
      45       15176 :   ASSERT_ALWAYS(n > 0);
      46             : 
      47       22968 :   while (n > 2 && (n & 1) == 0)
      48             :     {
      49        7792 :       b++;
      50        7792 :       n >>= 1;
      51             :     }
      52             :   
      53       15176 :   mpz_set (r, x);
      54             :   
      55       15176 :   mpz_init (t);
      56       15176 :   mpz_init (u);
      57             : 
      58       15176 :   if (n > 1)
      59             :     {
      60       15176 :       mpz_set (r, x);
      61       15176 :       mpz_mul (r, r, r);
      62       15176 :       mpz_sub_si (r, r, a);
      63       15176 :       mpz_sub_si (r, r, a); /* r = dickson(x, 2, a) */
      64             :       
      65       15176 :       mpz_set (t, x);    /* t = dickson(x, 1, a) */
      66             :       
      67       45232 :       for (i = 2; i < n; i++)
      68             :         {
      69       30056 :           mpz_mul_si (u, t, a);
      70       30056 :           mpz_set (t, r);     /* t = dickson(x, i, a) */
      71       30056 :           mpz_mul (r, r, x);
      72       30056 :           mpz_sub (r, r, u);  /* r = dickson(x, i+1, a) */
      73             :         }
      74             :     }
      75             :   
      76       22968 :   for ( ; b > 0; b--)
      77             :     {
      78        7792 :       mpz_mul (t, r, r); /* t = dickson(x, n, a) ^ 2 */
      79        7792 :       mpz_ui_pow_ui (u, abs (a), n);
      80        7792 :       if (n & 1 && a < 0)
      81        7792 :         mpz_neg (u, u);
      82        7792 :       mpz_mul_2exp (u, u, 1); /* u = 2 * a^n */
      83        7792 :       mpz_sub (r, t, u); /* r = dickson(x, 2*n, a) */
      84        7792 :       n <<= 1;
      85             :     }
      86             :   
      87       15176 :   mpz_clear (t);
      88       15176 :   mpz_clear (u);
      89       15176 : }
      90             : 
      91             : 
      92             : /* Init table to allow computation of
      93             : 
      94             :    Dickson_{E, a} (s + n*D), 
      95             : 
      96             :    for successive n, where Dickson_{E, a} is the Dickson polynomial 
      97             :    of degree E with parameter a. For a == 0, Dickson_{E, a} (x) = x^E .
      98             : 
      99             :    See Knuth, TAOCP vol.2, 4.6.4 and exercise 7 in 4.6.4, and
     100             :    "An FFT Extension of the Elliptic Curve Method of Factorization",
     101             :    Peter Montgomery, Dissertation, 1992, Chapter 5.
     102             : 
     103             :    Ternary return value.
     104             : */
     105             : 
     106             : static void
     107       20564 : fin_diff_coeff (listz_t coeffs, mpz_t s, mpz_t D, unsigned int E, 
     108             :                 int dickson_a)
     109             : {
     110             :   unsigned int i, k;
     111             :   mpz_t t;
     112             :   
     113       20564 :   mpz_init (t);
     114       20564 :   mpz_set (t, s);
     115             :   
     116       75802 :   for (i = 0; i <= E; i++)
     117             :     {
     118       55238 :       if (dickson_a != 0)         /* fd[i] = dickson_{E,a} (s+i*D) */
     119       15176 :         dickson (coeffs[i], t, E, dickson_a); 
     120             :       else                        /* fd[i] = (s+i*D)^E */
     121       40062 :         mpz_pow_ui (coeffs[i], t, E);
     122       55238 :       mpz_add (t, t, D);          /* t = s + i * D */
     123             :     }
     124             :   
     125       55238 :   for (k = 1; k <= E; k++)
     126      110981 :     for (i = E; i >= k; i--)
     127       76307 :       mpz_sub (coeffs[i], coeffs[i], coeffs[i-1]);
     128             :   
     129       20564 :   mpz_clear (t);
     130       20564 : }
     131             : 
     132             : 
     133             : /* Init several disjoint progressions for the computation of 
     134             : 
     135             :    Dickson_{E,a} (e * (i0 + i + n * d * k)), for 0 <= i < d * k   (1)
     136             :                   with gcd(e * (i0 + i), d) == 1, i == 1 (mod m),
     137             :                   where m divides d
     138             :    
     139             :    for successive n (the variable n does not appear here, it is the 
     140             :    application that called this function that wants to evaluate (1)
     141             :    for n = 0, 1, 2, ...
     142             :    
     143             :    This means there will be k sets of progressions, where each set contains
     144             :    eulerphi(d) progressions that generate the values of Dickson_{E,a} (x)
     145             :    with x coprime to d and 
     146             :    with i == 1 (mod m), where x == e * (i0 + i) (mod m).
     147             : 
     148             :    i0 may be a NULL pointer, in this case i0 = 0 is assumed.
     149             : 
     150             :    Return NULL if an error occurred.
     151             : */
     152             : 
     153             : listz_t
     154        2086 : init_progression_coeffs (mpz_t i0, const unsigned long d, 
     155             :                          const unsigned long e, const unsigned int k, 
     156             :                          const unsigned int m, const unsigned int E, 
     157             :                          const int dickson_a)
     158             : {
     159             :   unsigned int i, j, size_fd;
     160             :   mpz_t t, dke, em;
     161             :   listz_t fd;
     162             : 
     163             :   ASSERT (d % m == 0);
     164             : 
     165        2086 :   size_fd = k * (eulerphi(d) / eulerphi(m)) * (E + 1);
     166        2086 :   fd = (listz_t) malloc (size_fd * sizeof (mpz_t));
     167        2086 :   ASSERT_ALWAYS(fd != NULL);
     168       57324 :   for (i = 0; i < size_fd; i++)
     169       55238 :     mpz_init (fd[i]);
     170             : 
     171        2086 :   mpz_init (t);
     172        2086 :   if (i0 != NULL)
     173        2086 :     mpz_set (t, i0);
     174             :   
     175        2086 :   outputf (OUTPUT_TRACE, "init_progression_coeffs: i0 = %Zd, d = %u, e = %u, "
     176             :            "k = %u, m = %u, E = %u, a = %d, size_fd = %u\n", 
     177             :            t, d, e, k, m, E, dickson_a, size_fd);
     178             : 
     179             :   /* Due to the condition i == 1 (mod m) we start at i = 1 or i = 0,
     180             :      depending on whether m > 1 or m == 1 */
     181        2086 :   i = (m > 1) ? 1 : 0;
     182        2086 :   mpz_add_ui (t, t, (unsigned long) i);
     183        2086 :   mpz_mul_ui (t, t, e);
     184             :   /* Now t = e * (i0 + i + n * d * k), for n = 0 */
     185             :   
     186             :   /* dke = d * k * e, the common difference of the arithmetic progressions
     187             :      (it is the same for all arithmetic progressions we initialise) */
     188        2086 :   mpz_init (dke);
     189        2086 :   mpz_set_ui (dke, d);
     190        2086 :   mpz_mul_ui (dke, dke, k);
     191        2086 :   mpz_mul_ui (dke, dke, e);
     192             :   /* em = e * m, the value by which t advances if we increase i by m */
     193        2086 :   mpz_init (em);
     194        2086 :   mpz_set_ui (em, e);
     195        2086 :   mpz_mul_ui (em, em, (unsigned long) m);
     196             :   
     197       27647 :   for (j = 0; i < k * d; i += m)
     198             :     {
     199       25561 :       if (mpz_gcd_ui (NULL, t, d) == 1)
     200             :         {
     201       20564 :           outputf (OUTPUT_TRACE, "init_progression_coeffs: initing a "
     202             :                    "progression for Dickson_{%d,%d}(%Zd + n * %Zd)\n", 
     203             :                    E, dickson_a, t, dke);
     204             :           /* Initialise for the evaluation of Dickson_{E,a} (t + n*dke)
     205             :              for n = 0, 1, 2, ... */
     206       20564 :           fin_diff_coeff (fd + j, t, dke, E, dickson_a);
     207       20564 :           j += E + 1;
     208             :         } else
     209        4997 :           if (test_verbose (OUTPUT_TRACE))
     210         110 :             outputf (OUTPUT_TRACE, "init_progression_coeffs: NOT initing a "
     211             :                      "progression for Dickson_{%d,%d}(%Zd + n * %Zd), "
     212             :                      "gcd (%Zd, %u) == %u)\n", E, dickson_a, t, dke, t, d,
     213             :                      mpz_gcd_ui (NULL, t, d));
     214             :       /* We increase i by m, so we increase t by e*m */
     215       25561 :       mpz_add (t, t, em);
     216             :     }
     217             : 
     218        2086 :   mpz_clear (em);
     219        2086 :   mpz_clear (dke);
     220        2086 :   mpz_clear (t);
     221        2086 :   return fd;
     222             : }
     223             : 
     224             : void 
     225        1048 : init_roots_params (progression_params_t *params, const int S, 
     226             :                    const unsigned long d1, const unsigned long d2, 
     227             :                    const double cost)
     228             : {
     229             :   ASSERT (gcd (d1, d2) == 1);
     230             :   /* If S < 0, use degree |S| Dickson poly, otherwise use x^S */
     231        1048 :   params->S = abs (S);
     232        1048 :   params->dickson_a = (S < 0) ? -1 : 0;
     233             : 
     234             :   /* We only calculate Dickson_{S, a}(j * d2) * s where
     235             :      gcd (j, dsieve) == 1 and j == 1 (mod 6)
     236             :      by doing nr = eulerphi(dsieve)/2 separate progressions. */
     237             :   /* Now choose a value for dsieve. */
     238        1048 :   params->dsieve = 6;
     239        1048 :   params->nr = 1;
     240             : 
     241             :   /* Prospective saving by sieving out multiples of 5:
     242             :      d1 / params->dsieve * params->nr / 5 roots, each one costs S point adds
     243             :      Prospective cost increase:
     244             :      4 times as many progressions to init (that is, 3 * params->nr more),
     245             :      each costs ~ S * S * log_2(5 * dsieve * d2) / 2 point adds
     246             :      The params->nr and one S cancel.
     247             :   */
     248        1048 :   if (d1 % 5 == 0 &&
     249         998 :       d1 / params->dsieve / 5. * cost > 
     250         998 :       3. * params->S * log (5. * params->dsieve * d2) / 2.)
     251             :     {
     252         766 :       params->dsieve *= 5;
     253         766 :       params->nr *= 4;
     254             :     }
     255             : 
     256        1048 :   if (d1 % 7 == 0 &&
     257         692 :       d1 / params->dsieve / 7. * cost > 
     258         692 :       5. * params->S * log (7. * params->dsieve * d2) / 2.)
     259             :     {
     260         315 :       params->dsieve *= 7;
     261         315 :       params->nr *= 6;
     262             :     }
     263             : 
     264             : #if 0 /* commented out since not covered by unit tests */
     265             :   if (d1 % 11 == 0 &&
     266             :       d1 / params->dsieve / 11. * cost > 
     267             :       9. * params->S * log (11. * params->dsieve * d2) / 2.)
     268             :     {
     269             :       params->dsieve *= 11;
     270             :       params->nr *= 10;
     271             :     }
     272             : #endif
     273             : 
     274        1048 :   params->size_fd = params->nr * (params->S + 1);
     275        1048 :   params->next = 0;
     276        1048 :   params->rsieve = 1;
     277        1048 : }
     278             : 
     279             : double 
     280        1331 : memory_use (unsigned long dF, unsigned int sp_num, unsigned int Ftreelvl,
     281             :             mpmod_t modulus)
     282             : {
     283             :   double mem;
     284             :   
     285             :   /* printf ("memory_use (%lu, %d, %d, )\n", dF, sp_num, Ftreelvl); */
     286             : 
     287        1331 :   mem = 9.0; /* F:1, T:3*2, invF:1, G:1 */
     288        1331 :   mem += (double) Ftreelvl;
     289        1331 :   mem *= (double) dF;
     290        1331 :   mem += 2. * list_mul_mem (dF); /* Also in T */
     291             :   /* estimated memory for list_mult_n /
     292             :      wrap-case in PrerevertDivision respectively */
     293        1331 :   mem += (24.0 + 1.0) * (double) (sp_num ? MIN(MUL_NTT_THRESHOLD, dF) : dF);
     294        1331 :   mem *= (double) (mpz_size (modulus->orig_modulus)) * sizeof (mp_limb_t)
     295        1331 :          + sizeof (mpz_t);
     296             :   
     297        1331 :   if (sp_num)
     298        1129 :     mem += /* peak malloc in ecm_ntt.c */
     299        1129 :          (4.0 * dF * sp_num * sizeof (sp_t))
     300             :          
     301             :          /* mpzspv_normalise */
     302        1129 :          + (MPZSPV_NORMALISE_STRIDE * ((double) sp_num * 
     303        1129 :                 sizeof (sp_t) + 6.0 * sizeof (sp_t) + sizeof (float)))
     304             : 
     305             :          /* sp_F, sp_invF */
     306        1129 :          + ((1.0 + 2.0) * dF * sp_num * sizeof (sp_t));
     307             : 
     308        1331 :   return mem;
     309             : }
     310             : 
     311             : /* Input:  X is the point at end of stage 1
     312             :            modulus contains the number to factor
     313             :            B2min-B2 is the stage 2 range (we consider B2min is done)
     314             :            k0 is the number of blocks (if 0, use default)
     315             :            S is the exponent for Brent-Suyama's extension
     316             :            invtrick is non-zero iff one uses x+1/x instead of x.
     317             :            Cf "Speeding the Pollard and Elliptic Curve Methods
     318             :                of Factorization", Peter Montgomery, Math. of Comp., 1987,
     319             :                page 257: using x^(i^e)+1/x^(i^e) instead of x^(i^(2e))
     320             :                reduces the cost of Brent-Suyama's extension from 2*e
     321             :                to e+3 multiplications per value of i.
     322             :    Output: f is the factor found
     323             :    Return value: 2 (step number) iff a factor was found,
     324             :                  or ECM_ERROR if an error occurred.
     325             : */
     326             : int
     327        1048 : stage2 (mpz_t f, void *X, mpmod_t modulus, unsigned long dF, unsigned long k, 
     328             :         root_params_t *root_params, int use_ntt, char *TreeFilename,
     329             :         unsigned int curve_number, int (*stop_asap)(void))
     330             : {
     331             :   unsigned long i, sizeT;
     332             :   mpz_t n;
     333             :   listz_t F, G, H, T;
     334        1048 :   int youpi = ECM_NO_FACTOR_FOUND;
     335             :   long st, st0;
     336        1048 :   void *rootsG_state = NULL;
     337        1048 :   listz_t *Tree = NULL; /* stores the product tree for F */
     338             :   unsigned int lgk; /* ceil(log(k)/log(2)) */
     339        1048 :   listz_t invF = NULL;
     340             :   double mem;
     341        1048 :   mpzspm_t mpzspm = NULL;
     342        1048 :   mpzspv_t sp_F = NULL, sp_invF = NULL;
     343             :   
     344             :   /* check alloc. size of f */
     345        1048 :   mpres_realloc (f, modulus);
     346             : 
     347        1048 :   st0 = cputime ();
     348             : 
     349        1048 :   Fermat = 0;
     350        1048 :   if (modulus->repr == ECM_MOD_BASE2 && modulus->Fermat > 0)
     351             :     {
     352           4 :       Fermat = modulus->Fermat;
     353           4 :       use_ntt = 0; /* don't use NTT for Fermat numbers */
     354             :     }
     355             : 
     356        1048 :   if (use_ntt)
     357             :     {
     358         929 :       mpzspm = mpzspm_init (2 * dF, modulus->orig_modulus);
     359         929 :       ASSERT_ALWAYS(mpzspm != NULL);
     360             : 
     361         929 :       outputf (OUTPUT_VERBOSE,
     362             :           "Using %u small primes for NTT\n", mpzspm->sp_num);
     363             :     }
     364             : 
     365        1048 :   lgk = ceil_log2 (dF);
     366             : 
     367        1048 :   mem = memory_use (dF, use_ntt ? mpzspm->sp_num : 0,
     368             :       (TreeFilename == NULL) ? lgk : 0, modulus);
     369             : 
     370             :   /* we want at least two significant digits */
     371        1048 :   if (mem < 1048576.0)
     372         601 :     outputf (OUTPUT_VERBOSE, "Estimated memory usage: %1.0fKB\n", mem / 1024.);
     373         447 :   else if (mem < 1073741824.0)
     374         445 :     outputf (OUTPUT_VERBOSE, "Estimated memory usage: %1.2fMB\n", 
     375             :              mem / 1048576.);
     376             :   else
     377           2 :     outputf (OUTPUT_VERBOSE, "Estimated memory usage: %1.2fGB\n", 
     378             :              mem / 1073741824.);
     379             : 
     380        1048 :   F = init_list2 (dF + 1, mpz_sizeinbase (modulus->orig_modulus, 2) + 
     381             :                           3 * GMP_NUMB_BITS);
     382        1048 :   ASSERT_ALWAYS(F != NULL);
     383             : 
     384        1048 :   sizeT = 3 * dF + list_mul_mem (dF);
     385        1048 :   if (dF > 3)
     386         998 :     sizeT += dF;
     387        1048 :   T = init_list2 (sizeT, 2 * mpz_sizeinbase (modulus->orig_modulus, 2) + 
     388             :                          3 * GMP_NUMB_BITS);
     389        1048 :   ASSERT_ALWAYS(T != NULL);
     390        1048 :   H = T;
     391             : 
     392             :   /* needs dF+1 cells in T */
     393        1048 :   youpi = ecm_rootsF (f, F, root_params, dF, (curve*) X, modulus);
     394             : 
     395        1048 :   if (youpi != ECM_NO_FACTOR_FOUND)
     396             :     {
     397          10 :       if (youpi != ECM_ERROR)
     398          10 :         youpi = ECM_FACTOR_FOUND_STEP2;
     399          10 :       goto clear_T;
     400             :     }
     401        1038 :   if (stop_asap != NULL && (*stop_asap)())
     402           0 :     goto clear_T;
     403             : 
     404        1038 :   if (test_verbose (OUTPUT_TRACE))
     405             :     {
     406             :       unsigned long j;
     407        4978 :       for (j = 0; j < dF; j++)
     408        4968 :         outputf (OUTPUT_TRACE, "f_%lu = %Zd\n", j, F[j]);
     409             :     }
     410             : 
     411             :   /* ----------------------------------------------
     412             :      |   F    |  invF  |   G   |         T        |
     413             :      ----------------------------------------------
     414             :      | rootsF |  ???   |  ???  |      ???         |
     415             :      ---------------------------------------------- */
     416             : 
     417        1038 :   if (TreeFilename == NULL)
     418             :     {
     419         922 :       Tree = (listz_t*) malloc (lgk * sizeof (listz_t));
     420         922 :       ASSERT_ALWAYS(Tree != NULL);
     421        7674 :       for (i = 0; i < lgk; i++)
     422             :         {
     423        6752 :           Tree[i] = init_list2 (dF, mpz_sizeinbase (modulus->orig_modulus, 2) 
     424             :                                     + GMP_NUMB_BITS);
     425        6752 :           ASSERT_ALWAYS(Tree[i] != NULL);
     426             :         }
     427             :     }
     428             :   else
     429         116 :     Tree = NULL;
     430             :   
     431             : #ifdef TELLEGEN_DEBUG
     432             :   outputf (OUTPUT_ALWAYS, "Roots = ");
     433             :   print_list (os, F, dF);
     434             : #endif
     435        1038 :   mpz_init_set (n, modulus->orig_modulus);
     436        1038 :   st = cputime ();
     437        1038 :   if (TreeFilename != NULL)
     438             :     {
     439             :       FILE *TreeFile;
     440         116 :       char *fullname = (char *) malloc (strlen (TreeFilename) + 1 + 2 + 1);
     441             :       int ret;
     442         116 :       ASSERT_ALWAYS(fullname != NULL);
     443             :       
     444         967 :       for (i = lgk; i > 0; i--)
     445             :         {
     446         851 :           if (stop_asap != NULL && (*stop_asap)())
     447           0 :             goto free_Tree_i;
     448         851 :           sprintf (fullname, "%s.%lu", TreeFilename, i - 1);
     449             :           
     450         851 :           TreeFile = fopen (fullname, "wb");
     451         851 :           if (TreeFile == NULL)
     452             :             {
     453           0 :               outputf (OUTPUT_ERROR, 
     454             :                        "Error opening file for product tree of F\n");
     455           0 :               youpi = ECM_ERROR;
     456           0 :               goto free_Tree_i;
     457             :             }
     458             :           
     459         835 :           ret = (use_ntt) ? ntt_PolyFromRoots_Tree (F, F, dF, T, i - 1,
     460             :                                                     mpzspm, NULL, TreeFile)
     461         851 :             : PolyFromRoots_Tree (F, F, dF, T, i - 1, n, NULL, TreeFile, 0);
     462         851 :           if (ret == ECM_ERROR)
     463             :             {
     464           0 :               fclose (TreeFile);
     465           0 :               youpi = ECM_ERROR;
     466           0 :               goto free_Tree_i;
     467             :             }
     468         851 :           fclose (TreeFile);
     469             :         }
     470         116 :       free (fullname);
     471             :     }
     472             :   else
     473             :     {
     474             :       /* TODO: how to check for stop_asap() here? */
     475         922 :       if (use_ntt)
     476         806 :         ntt_PolyFromRoots_Tree (F, F, dF, T, -1, mpzspm, Tree, NULL);
     477             :       else
     478         116 :         PolyFromRoots_Tree (F, F, dF, T, -1, n, Tree, NULL, 0);
     479             :     }
     480             :   
     481             :   
     482        1038 :   if (test_verbose (OUTPUT_TRACE))
     483             :     {
     484             :       unsigned long j;
     485        4978 :       for (j = 0; j < dF; j++)
     486        4968 :         outputf (OUTPUT_TRACE, "F[%lu] = %Zd\n", j, F[j]);
     487             :     }
     488        1038 :   outputf (OUTPUT_VERBOSE, "Building F from its roots took %ldms\n", 
     489             :            elltime (st, cputime ()));
     490             : 
     491        1038 :   if (stop_asap != NULL && (*stop_asap)())
     492           0 :     goto free_Tree_i;
     493             : 
     494             : 
     495             :   /* needs dF+list_mul_mem(dF/2) cells in T */
     496             : 
     497        1038 :   mpz_set_ui (F[dF], 1); /* the leading monic coefficient needs to be stored
     498             :                              explicitly for PrerevertDivision */
     499             : 
     500             :   /* ----------------------------------------------
     501             :      |   F    |  invF  |   G   |         T        |
     502             :      ----------------------------------------------
     503             :      |  F(x)  |  ???   |  ???  |      ???         |
     504             :      ---------------------------------------------- */
     505             : 
     506             :   /* G*H has degree 2*dF-2, hence we must cancel dF-1 coefficients
     507             :      to get degree dF-1 */
     508        1038 :   if (dF > 1)
     509             :     {
     510             :       /* only dF-1 coefficients of 1/F are needed to reduce G*H,
     511             :          but we need one more for TUpTree */
     512        1038 :       invF = init_list2 (dF + 1, mpz_sizeinbase (modulus->orig_modulus, 2) + 
     513             :                                  2 * GMP_NUMB_BITS);
     514        1038 :       ASSERT_ALWAYS(invF != NULL);
     515        1038 :       st = cputime ();
     516             :       
     517        1038 :       if (use_ntt)
     518             :         {
     519         920 :           sp_F = mpzspv_init (dF, mpzspm);
     520         920 :           mpzspv_from_mpzv (sp_F, 0, F, dF, mpzspm);
     521         920 :           mpzspv_to_ntt (sp_F, 0, dF, dF, 1, mpzspm);
     522             :           
     523         920 :           ntt_PolyInvert (invF, F + 1, dF, T, mpzspm);
     524         920 :           sp_invF = mpzspv_init (2 * dF, mpzspm);
     525         920 :           mpzspv_from_mpzv (sp_invF, 0, invF, dF, mpzspm);
     526         920 :           mpzspv_to_ntt (sp_invF, 0, dF, 2 * dF, 0, mpzspm);
     527             :         }
     528             :       else
     529         118 :         PolyInvert (invF, F + 1, dF, T, n);
     530             :       
     531             :       /* now invF[0..dF-1] = Quo(x^(2dF-1), F) */
     532        1038 :       outputf (OUTPUT_VERBOSE, "Computing 1/F took %ldms\n",
     533             :                elltime (st, cputime ()));
     534             :       
     535             :       /* ----------------------------------------------
     536             :          |   F    |  invF  |   G   |         T        |
     537             :          ----------------------------------------------
     538             :          |  F(x)  | 1/F(x) |  ???  |      ???         |
     539             :          ---------------------------------------------- */
     540             :     }
     541             : 
     542        1038 :   if (stop_asap != NULL && (*stop_asap)())
     543           0 :     goto clear_invF;
     544             : 
     545             :   /* start computing G with dF roots.
     546             : 
     547             :      In the non CM case, roots are at i0*d, (i0+1)*d, (i0+2)*d, ... 
     548             :      where i0*d <= B2min < (i0+1)*d .
     549             :   */
     550        1038 :   G = init_list2 (dF, mpz_sizeinbase (modulus->orig_modulus, 2) + 
     551             :                       3 * GMP_NUMB_BITS);
     552        1038 :   ASSERT_ALWAYS(G != NULL);
     553             : 
     554        1038 :   st = cputime ();
     555        1038 :   rootsG_state = ecm_rootsG_init (f, (curve *) X, root_params, dF, k, 
     556             :                                   modulus);
     557             : 
     558             :   /* rootsG_state=NULL if an error occurred or (ecm only) a factor was found */
     559        1038 :   if (rootsG_state == NULL)
     560             :     {
     561             :       /* ecm: f = -1 if an error occurred */
     562           0 :       youpi = (mpz_cmp_si (f, -1)) ? ECM_FACTOR_FOUND_STEP2 : ECM_ERROR;
     563           0 :       goto clear_G;
     564             :     }
     565             : 
     566        1038 :   if (stop_asap != NULL && (*stop_asap)())
     567           0 :     goto clear_fd;
     568             : 
     569        4389 :   for (i = 0; i < k; i++)
     570             :     {
     571             :       /* needs dF+1 cells in T+dF */
     572        3362 :         youpi = ecm_rootsG (f, G, dF, (ecm_roots_state_t *) rootsG_state, 
     573             :                               modulus);
     574             : 
     575        3362 :       if (test_verbose (OUTPUT_TRACE))
     576             :         {
     577             :           unsigned long j;
     578        4978 :           for (j = 0; j < dF; j++)
     579        4968 :             outputf (OUTPUT_TRACE, "g_%lu = %Zd\n", j, G[j]);
     580             :         }
     581             : 
     582             :       ASSERT(youpi != ECM_ERROR); /* xxx_rootsG cannot fail */
     583        3362 :       if (youpi) /* factor found */
     584             :         {
     585          11 :           youpi = ECM_FACTOR_FOUND_STEP2;
     586          11 :           goto clear_fd;
     587             :         }
     588             : 
     589        3351 :     if (stop_asap != NULL && (*stop_asap)())
     590           0 :       goto clear_fd;
     591             : 
     592             :   /* -----------------------------------------------
     593             :      |   F    |  invF  |   G    |         T        |
     594             :      -----------------------------------------------
     595             :      |  F(x)  | 1/F(x) | rootsG |      ???         |
     596             :      ----------------------------------------------- */
     597             : 
     598        3351 :       st = cputime ();
     599             : 
     600        3351 :       if (use_ntt)
     601        3108 :         ntt_PolyFromRoots (G, G, dF, T + dF, mpzspm);
     602             :       else
     603         243 :         PolyFromRoots (G, G, dF, T + dF, n);
     604             : 
     605        3351 :       if (test_verbose (OUTPUT_TRACE))
     606             :         {
     607             :           unsigned long j;
     608          10 :           outputf (OUTPUT_TRACE, "G(x) = x^%lu ", dF);
     609        4978 :           for (j = 0; j < dF; j++)
     610        4968 :             outputf (OUTPUT_TRACE, "+ (%Zd * x^%lu)", G[j], j);
     611          10 :           outputf (OUTPUT_TRACE, "\n");
     612             :         }
     613             : 
     614             :       /* needs 2*dF+list_mul_mem(dF/2) cells in T */
     615        3351 :       outputf (OUTPUT_VERBOSE, "Building G from its roots took %ldms\n", 
     616             :                elltime (st, cputime ()));
     617             : 
     618        3351 :     if (stop_asap != NULL && (*stop_asap)())
     619           0 :       goto clear_fd;
     620             : 
     621             :   /* -----------------------------------------------
     622             :      |   F    |  invF  |   G    |         T        |
     623             :      -----------------------------------------------
     624             :      |  F(x)  | 1/F(x) |  G(x)  |      ???         |
     625             :      ----------------------------------------------- */
     626             : 
     627        3351 :       if (i == 0)
     628             :         {
     629        1036 :           list_sub (H, G, F, dF); /* coefficients 1 of degree cancel,
     630             :                                      thus T is of degree < dF */
     631        1036 :           list_mod (H, H, dF, n);
     632             :           /* ------------------------------------------------
     633             :              |   F    |  invF  |    G    |         T        |
     634             :              ------------------------------------------------
     635             :              |  F(x)  | 1/F(x) |  ???    |G(x)-F(x)|  ???   |
     636             :              ------------------------------------------------ */
     637             :         }
     638             :       else
     639             :         {
     640             :           /* since F and G are monic of same degree, G mod F = G - F */
     641        2315 :           list_sub (G, G, F, dF);
     642        2315 :           list_mod (G, G, dF, n);
     643             : 
     644             :           /* ------------------------------------------------
     645             :              |   F    |  invF  |    G    |         T        |
     646             :              ------------------------------------------------
     647             :              |  F(x)  | 1/F(x) |G(x)-F(x)|  H(x)  |         |
     648             :              ------------------------------------------------ */
     649             : 
     650        2315 :           st = cputime ();
     651             :           /* previous G mod F is in H, with degree < dF, i.e. dF coefficients:
     652             :              requires 3dF-1+list_mul_mem(dF) cells in T */
     653        2315 :           if (use_ntt)
     654             :             {
     655        2188 :               ntt_mul (T + dF, G, H, dF, T + 3 * dF, 0, mpzspm);
     656        2188 :               list_mod (H, T + dF, 2 * dF, n);
     657             :             }
     658             :           else
     659         127 :             list_mulmod (H, T + dF, G, H, dF, T + 3 * dF, n);
     660             : 
     661        2315 :           outputf (OUTPUT_VERBOSE, "Computing G * H took %ldms\n", 
     662             :                    elltime (st, cputime ()));
     663             : 
     664        2315 :           if (stop_asap != NULL && (*stop_asap)())
     665           0 :             goto clear_fd;
     666             : 
     667             :           /* ------------------------------------------------
     668             :              |   F    |  invF  |    G    |         T        |
     669             :              ------------------------------------------------
     670             :              |  F(x)  | 1/F(x) |G(x)-F(x)| G * H  |         |
     671             :              ------------------------------------------------ */
     672             : 
     673        2315 :           st = cputime ();
     674             : 
     675        2315 :           if (use_ntt)
     676             :             {
     677        2188 :               ntt_PrerevertDivision (H, F, invF + 1, sp_F, sp_invF, dF,
     678             :                   T + 2 * dF, mpzspm);
     679             :             }
     680             :           else
     681             :             {
     682         127 :               if (PrerevertDivision (H, F, invF + 1, dF, T + 2 * dF, n))
     683             :                 {
     684           0 :                   youpi = ECM_ERROR;
     685           0 :                   goto clear_fd;
     686             :                 }
     687             :             }
     688             :           
     689        2315 :           outputf (OUTPUT_VERBOSE, "Reducing  G * H mod F took %ldms\n", 
     690             :                    elltime (st, cputime ()));
     691             : 
     692        2315 :           if (stop_asap != NULL && (*stop_asap)())
     693           0 :             goto clear_fd;
     694             :         }
     695             :     }
     696             :   
     697        1027 :   clear_list (F, dF + 1);
     698        1027 :   F = NULL;
     699        1027 :   clear_list (G, dF);
     700        1027 :   G = NULL;
     701        1027 :   st = cputime ();
     702        1027 :   if (use_ntt)
     703         911 :     youpi = ntt_polyevalT (T, dF, Tree, T + dF + 1, sp_invF,
     704             :         mpzspm, TreeFilename);
     705             :   else
     706         116 :     youpi = polyeval_tellegen (T, dF, Tree, T + dF + 1, sizeT - dF - 1, invF,
     707             :         n, TreeFilename);
     708             : 
     709        1027 :   if (youpi)
     710             :     {
     711           0 :       outputf (OUTPUT_ERROR, "Error, not enough memory\n");
     712           0 :       goto clear_fd;
     713             :     }
     714             : 
     715        1027 :   if (test_verbose (OUTPUT_TRACE))
     716             :     {
     717             :       unsigned long j;
     718        4978 :       for (j = 0; j < dF; j++)
     719        4968 :         outputf (OUTPUT_TRACE, "G(x_%lu) = %Zd\n", j, T[j]);
     720             :     }
     721             : 
     722        1027 :   outputf (OUTPUT_VERBOSE, "Computing polyeval(F,G) took %ldms\n", 
     723             :            elltime (st, cputime ()));
     724             : 
     725        1027 :   st = cputime ();
     726        1027 :   list_mulup (T, dF, n, T[dF]);
     727        1027 :   outputf (OUTPUT_VERBOSE, "Computing product of all F(g_i) took %ldms\n", 
     728             :            elltime (st, cputime ()));
     729             : 
     730        1027 :   mpz_gcd (f, T[dF - 1], n);
     731        1027 :   if (mpz_cmp_ui (f, 1) > 0)
     732         784 :     youpi = ECM_FACTOR_FOUND_STEP2;
     733             :   else
     734             :     /* Here, mpz_cmp_ui (f, 1) == 0, i.e. no factor was found */
     735         243 :     outputf (OUTPUT_RESVERBOSE, "Product of G(f_i) = %Zd\n", T[0]);
     736             : 
     737        1038 : clear_fd:
     738        1038 :   ecm_rootsG_clear ((ecm_roots_state_t *) rootsG_state, modulus);
     739             : 
     740        1038 : clear_G:
     741        1038 :   clear_list (G, dF);
     742        1038 : clear_invF:
     743        1038 :   clear_list (invF, dF + 1);
     744             : 
     745        1038 :   if (use_ntt)
     746             :     {
     747         920 :       mpzspv_clear (sp_F, mpzspm);
     748         920 :       mpzspv_clear (sp_invF, mpzspm);
     749             :     }
     750         118 : free_Tree_i:
     751        1038 :   if (Tree != NULL)
     752             :     {
     753        7674 :       for (i = 0; i < lgk; i++)
     754        6752 :         clear_list (Tree[i], dF);
     755         922 :       free (Tree);
     756             :     }
     757             :   /* the trees are already cleared by ntt_polyevalT or polyeval_tellegen */
     758        1038 :   mpz_clear (n);
     759             : 
     760        1048 : clear_T:
     761        1048 :   clear_list (T, sizeT);
     762        1048 :   clear_list (F, dF + 1);
     763             : 
     764        1048 :   if (use_ntt)
     765         929 :     mpzspm_clear (mpzspm);
     766             :   
     767        1048 :   if (Fermat)
     768           4 :     F_clear ();
     769             :   
     770             : 
     771        1048 :   if (stop_asap == NULL || !(*stop_asap)())
     772             :     {
     773        1048 :       st0 = elltime (st0, cputime ());
     774        1048 :       if (curve_number > 0) {
     775           0 :         outputf (OUTPUT_NORMAL, "Curve %d ", curve_number);
     776             :       }
     777        1048 :       outputf (OUTPUT_NORMAL, "Step 2 took %ldms\n", st0);
     778             :     }
     779             : 
     780        1048 :   return youpi;
     781             : }

Generated by: LCOV version 1.14