Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux1.fr machine (x86_64 architecture), and agregate them in the final report:

The target is 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - ifactor1.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 16781-61614ca) Lines: 1381 1685 82.0 %
Date: 2014-09-14 Functions: 75 84 89.3 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 1043 1519 68.7 %

           Branch data     Line data    Source code
       1                 :            : /* Copyright (C) 2000  The PARI group.
       2                 :            : 
       3                 :            : This file is part of the PARI/GP package.
       4                 :            : 
       5                 :            : PARI/GP is free software; you can redistribute it and/or modify it under the
       6                 :            : terms of the GNU General Public License as published by the Free Software
       7                 :            : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8                 :            : ANY WARRANTY WHATSOEVER.
       9                 :            : 
      10                 :            : Check the License for details. You should have received a copy of it, along
      11                 :            : with the package; see the file 'COPYING'. If not, write to the Free Software
      12                 :            : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13                 :            : #include "pari.h"
      14                 :            : #include "paripriv.h"
      15                 :            : 
      16                 :            : /***********************************************************************/
      17                 :            : /**                                                                   **/
      18                 :            : /**                       PRIMES IN SUCCESSION                        **/
      19                 :            : /** (abstracted by GN 1998Aug21 mainly for use in ellfacteur() below) **/
      20                 :            : /**                                                                   **/
      21                 :            : /***********************************************************************/
      22                 :            : 
      23                 :            : /* map from prime residue classes mod 210 to their numbers in {0...47}.
      24                 :            :  * Subscripts into this array take the form ((k-1)%210)/2, ranging from
      25                 :            :  * 0 to 104.  Unused entries are */
      26                 :            : #define NPRC 128                /* non-prime residue class */
      27                 :            : 
      28                 :            : static unsigned char prc210_no[] = {
      29                 :            :   0, NPRC, NPRC, NPRC, NPRC, 1, 2, NPRC, 3, 4, NPRC, /* 21 */
      30                 :            :   5, NPRC, NPRC, 6, 7, NPRC, NPRC, 8, NPRC, 9, /* 41 */
      31                 :            :   10, NPRC, 11, NPRC, NPRC, 12, NPRC, NPRC, 13, 14, NPRC, /* 63 */
      32                 :            :   NPRC, 15, NPRC, 16, 17, NPRC, NPRC, 18, NPRC, 19, /* 83 */
      33                 :            :   NPRC, NPRC, 20, NPRC, NPRC, NPRC, 21, NPRC, 22, 23, NPRC, /* 105 */
      34                 :            :   24, 25, NPRC, 26, NPRC, NPRC, NPRC, 27, NPRC, NPRC, /* 125 */
      35                 :            :   28, NPRC, 29, NPRC, NPRC, 30, 31, NPRC, 32, NPRC, NPRC, /* 147 */
      36                 :            :   33, 34, NPRC, NPRC, 35, NPRC, NPRC, 36, NPRC, 37, /* 167 */
      37                 :            :   38, NPRC, 39, NPRC, NPRC, 40, 41, NPRC, NPRC, 42, NPRC, /* 189 */
      38                 :            :   43, 44, NPRC, 45, 46, NPRC, NPRC, NPRC, NPRC, 47, /* 209 */
      39                 :            : };
      40                 :            : 
      41                 :            : #if 0
      42                 :            : /* map from prime residue classes mod 210 (by number) to their smallest
      43                 :            :  * positive representatives */
      44                 :            : static unsigned char prc210_rp[] = { /* 19 + 15 + 14 = [0..47] */
      45                 :            :   1, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79,
      46                 :            :   83, 89, 97, 101, 103, 107, 109, 113, 121, 127, 131, 137, 139, 143, 149,
      47                 :            :   151, 157, 163, 167, 169, 173, 179, 181, 187, 191, 193, 197, 199, 209,
      48                 :            : };
      49                 :            : #endif
      50                 :            : 
      51                 :            : /* first differences of the preceding */
      52                 :            : static unsigned char prc210_d1[] = {
      53                 :            :   10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6,
      54                 :            :   4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6,
      55                 :            :   2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 2,
      56                 :            : };
      57                 :            : 
      58                 :            : /* return 0 for overflow */
      59                 :            : ulong
      60                 :    4514589 : unextprime(ulong n)
      61                 :            : {
      62                 :            :   long rc, rc0, rcd, rcn;
      63                 :            : 
      64   [ +  +  +  +  :    4514589 :   switch(n) {
                      + ]
      65                 :        140 :     case 0: case 1: case 2: return 2;
      66                 :        147 :     case 3: return 3;
      67                 :        420 :     case 4: case 5: return 5;
      68                 :        406 :     case 6: case 7: return 7;
      69                 :            :   }
      70                 :            : #ifdef LONG_IS_64BIT
      71         [ +  + ]:    3868734 :   if (n > (ulong)-59) return 0;
      72                 :            : #else
      73         [ +  + ]:     644742 :   if (n > (ulong)-5) return 0;
      74                 :            : #endif
      75                 :            :   /* here n > 7 */
      76                 :    4513461 :   n |= 1; /* make it odd */
      77                 :    4513461 :   rc = rc0 = n % 210;
      78                 :            :   /* find next prime residue class mod 210 */
      79                 :            :   for(;;)
      80                 :            :   {
      81                 :    9172942 :     rcn = (long)(prc210_no[rc>>1]);
      82         [ +  + ]:    9172942 :     if (rcn != NPRC) break;
      83                 :    4659481 :     rc += 2; /* cannot wrap since 209 is coprime and rc odd */
      84                 :    4659481 :   }
      85         [ +  + ]:    4513461 :   if (rc > rc0) n += rc - rc0;
      86                 :            :   /* now find an actual (pseudo)prime */
      87                 :            :   for(;;)
      88                 :            :   {
      89         [ +  + ]:   12668981 :     if (uisprime(n)) break;
      90                 :    8155520 :     rcd = prc210_d1[rcn];
      91         [ +  + ]:    8155520 :     if (++rcn > 47) rcn = 0;
      92                 :    8155520 :     n += rcd;
      93                 :    8155520 :   }
      94                 :    4514589 :   return n;
      95                 :            : }
      96                 :            : 
      97                 :            : GEN
      98                 :    2026612 : nextprime(GEN n)
      99                 :            : {
     100                 :            :   long rc, rc0, rcd, rcn;
     101                 :    2026612 :   pari_sp av = avma;
     102                 :            : 
     103         [ -  + ]:    2026612 :   if (typ(n) != t_INT)
     104                 :            :   {
     105                 :          0 :     n = gceil(n);
     106         [ #  # ]:          0 :     if (typ(n) != t_INT) pari_err_TYPE("nextprime",n);
     107                 :            :   }
     108         [ +  + ]:    2026612 :   if (signe(n) <= 0) { avma = av; return gen_2; }
     109         [ +  + ]:    2026521 :   if (lgefint(n) == 3)
     110                 :            :   {
     111                 :    2026289 :     ulong k = unextprime(uel(n,2));
     112                 :    2026289 :     avma = av;
     113         [ +  - ]:    2026289 :     if (k) return utoipos(k);
     114                 :            : #ifdef LONG_IS_64BIT
     115                 :          0 :     return uutoi(1,13);
     116                 :            : #else
     117                 :          0 :     return uutoi(1,15);
     118                 :            : #endif
     119                 :            :   }
     120                 :            :   /* here n > 7 */
     121         [ +  + ]:        232 :   if (!mod2(n)) n = addsi(1,n);
     122                 :        232 :   rc = rc0 = smodis(n, 210);
     123                 :            :   /* find next prime residue class mod 210 */
     124                 :            :   for(;;)
     125                 :            :   {
     126                 :        459 :     rcn = (long)(prc210_no[rc>>1]);
     127         [ +  + ]:        459 :     if (rcn != NPRC) break;
     128                 :        227 :     rc += 2; /* cannot wrap since 209 is coprime and rc odd */
     129                 :        227 :   }
     130         [ +  + ]:        232 :   if (rc > rc0) n = addsi(rc - rc0, n);
     131                 :            :   /* now find an actual (pseudo)prime */
     132                 :            :   for(;;)
     133                 :            :   {
     134         [ +  + ]:       5625 :     if (BPSW_psp(n)) break;
     135                 :       5393 :     rcd = prc210_d1[rcn];
     136         [ +  + ]:       5393 :     if (++rcn > 47) rcn = 0;
     137                 :       5393 :     n = addsi(rcd, n);
     138                 :       5393 :   }
     139         [ -  + ]:        232 :   if (avma == av) return icopy(n);
     140                 :    2026612 :   return gerepileuptoint(av, n);
     141                 :            : }
     142                 :            : 
     143                 :            : ulong
     144                 :         25 : uprecprime(ulong n)
     145                 :            : {
     146                 :            :   long rc, rc0, rcd, rcn;
     147                 :            :   { /* check if n <= 10 */
     148         [ +  + ]:         25 :     if (n <= 1)  return 0;
     149         [ -  + ]:         18 :     if (n == 2)  return 2;
     150         [ -  + ]:         18 :     if (n <= 4)  return 3;
     151         [ -  + ]:         18 :     if (n <= 6)  return 5;
     152         [ -  + ]:         18 :     if (n <= 10) return 7;
     153                 :            :   }
     154                 :            :   /* here n >= 11 */
     155         [ +  - ]:         18 :   if (!(n % 2)) n--;
     156                 :         18 :   rc = rc0 = n % 210;
     157                 :            :   /* find previous prime residue class mod 210 */
     158                 :            :   for(;;)
     159                 :            :   {
     160                 :         36 :     rcn = (long)(prc210_no[rc>>1]);
     161         [ +  + ]:         36 :     if (rcn != NPRC) break;
     162                 :         18 :     rc -= 2; /* cannot wrap since 1 is coprime and rc odd */
     163                 :         18 :   }
     164         [ +  - ]:         18 :   if (rc < rc0) n += rc - rc0;
     165                 :            :   /* now find an actual (pseudo)prime */
     166                 :            :   for(;;)
     167                 :            :   {
     168         [ +  + ]:         36 :     if (uisprime(n)) break;
     169         [ -  + ]:         18 :     if (--rcn < 0) rcn = 47;
     170                 :         18 :     rcd = prc210_d1[rcn];
     171                 :         18 :     n -= rcd;
     172                 :         18 :   }
     173                 :         25 :   return n;
     174                 :            : }
     175                 :            : 
     176                 :            : GEN
     177                 :         35 : precprime(GEN n)
     178                 :            : {
     179                 :            :   long rc, rc0, rcd, rcn;
     180                 :         35 :   pari_sp av = avma;
     181                 :            : 
     182         [ -  + ]:         35 :   if (typ(n) != t_INT)
     183                 :            :   {
     184                 :          0 :     n = gfloor(n);
     185         [ #  # ]:          0 :     if (typ(n) != t_INT) pari_err_TYPE("nextprime",n);
     186                 :            :   }
     187         [ -  + ]:         35 :   if (signe(n) <= 0) { avma = av; return gen_0; }
     188         [ +  + ]:         35 :   if (lgefint(n) <= 3)
     189                 :            :   {
     190                 :         25 :     ulong k = uel(n,2);
     191                 :         25 :     avma = av;
     192                 :         25 :     return utoi(uprecprime(k));
     193                 :            :   }
     194         [ +  - ]:         10 :   if (!mod2(n)) n = addsi(-1,n);
     195                 :         10 :   rc = rc0 = smodis(n, 210);
     196                 :            :   /* find previous prime residue class mod 210 */
     197                 :            :   for(;;)
     198                 :            :   {
     199                 :         20 :     rcn = (long)(prc210_no[rc>>1]);
     200         [ +  + ]:         20 :     if (rcn != NPRC) break;
     201                 :         10 :     rc -= 2; /* cannot wrap since 1 is coprime and rc odd */
     202                 :         10 :   }
     203         [ +  - ]:         10 :   if (rc < rc0) n = addsi(rc - rc0, n);
     204                 :            :   /* now find an actual (pseudo)prime */
     205                 :            :   for(;;)
     206                 :            :   {
     207         [ +  + ]:         48 :     if (BPSW_psp(n)) break;
     208         [ -  + ]:         38 :     if (--rcn < 0) rcn = 47;
     209                 :         38 :     rcd = prc210_d1[rcn];
     210                 :         38 :     n = addsi(-rcd, n);
     211                 :         38 :   }
     212         [ -  + ]:         10 :   if (avma == av) return icopy(n);
     213                 :         35 :   return gerepileuptoint(av, n);
     214                 :            : }
     215                 :            : 
     216                 :            : /* Find next single-word prime strictly larger than p.
     217                 :            :  * If **d is non-NULL (somewhere in a diffptr), this is p + *(*d)++;
     218                 :            :  * otherwise imitate nextprime().
     219                 :            :  * *rcn = NPRC or the correct residue class for the current p; we'll use this
     220                 :            :  * to track the current prime residue class mod 210 once we're out of range of
     221                 :            :  * the diffptr table, and we'll update it before that if it isn't NPRC.
     222                 :            :  *
     223                 :            :  * *q is incremented whenever q!=NULL and we wrap from 209 mod 210 to
     224                 :            :  * 1 mod 210
     225                 :            :  * k =  second argument for MR_Jaeschke(). --GN1998Aug22 */
     226                 :            : ulong
     227                 :     235956 : snextpr(ulong p, byteptr *d, long *rcn, long *q, long k)
     228                 :            : {
     229                 :            :   ulong n;
     230         [ +  - ]:     235956 :   if (**d)
     231                 :            :   {
     232                 :     235956 :     byteptr dd = *d;
     233                 :     235956 :     long d1 = 0;
     234                 :            : 
     235                 :     235956 :     NEXT_PRIME_VIADIFF(d1,dd);
     236                 :            :     /* d1 = nextprime(p+1) - p */
     237         [ +  + ]:     235956 :     if (*rcn != NPRC)
     238                 :            :     {
     239                 :     233184 :       long rcn0 = *rcn;
     240         [ +  + ]:     745850 :       while (d1 > 0)
     241                 :            :       {
     242                 :     512666 :         d1 -= prc210_d1[*rcn];
     243 [ +  + ][ +  - ]:     512666 :         if (++*rcn > 47) { *rcn = 0; if (q) (*q)++; }
     244                 :            :       }
     245         [ -  + ]:     233184 :       if (d1 < 0)
     246                 :            :       {
     247                 :          0 :         char *s=stack_sprintf("snextpr: %lu!=prc210_rp[%ld] mod 210\n",p,rcn0);
     248                 :          0 :         pari_err_BUG(s);
     249                 :            :       }
     250                 :            :     }
     251                 :     235956 :     NEXT_PRIME_VIADIFF(p,*d);
     252                 :     235956 :     return p;
     253                 :            :   }
     254                 :            :   /* we are beyond the diffptr table */
     255         [ #  # ]:          0 :   if (*rcn == NPRC)
     256                 :            :   { /* initialize */
     257                 :          0 :     *rcn = prc210_no[(p % 210) >> 1];
     258         [ #  # ]:          0 :     if (*rcn == NPRC)
     259                 :            :     {
     260                 :          0 :       char *s = stack_sprintf("snextpr: %lu should have been prime\n", p);
     261                 :          0 :       pari_err_BUG(s);
     262                 :            :     }
     263                 :            :   }
     264                 :            :   /* look for the next one */
     265                 :          0 :   n = p + prc210_d1[*rcn];
     266         [ #  # ]:          0 :   if (++*rcn > 47) *rcn = 0;
     267         [ #  # ]:          0 :   while (!Fl_MR_Jaeschke(n, k))
     268                 :            :   {
     269                 :          0 :     n += prc210_d1[*rcn];
     270         [ #  # ]:          0 :     if (n <= 11) pari_err_OVERFLOW("snextpr");
     271 [ #  # ][ #  # ]:          0 :     if (++*rcn > 47) { *rcn = 0; if (q) (*q)++; }
     272                 :            :   }
     273                 :     235956 :   return n;
     274                 :            : }
     275                 :            : 
     276                 :            : /********************************************************************/
     277                 :            : /**                                                                **/
     278                 :            : /**                     INTEGER FACTORIZATION                      **/
     279                 :            : /**                                                                **/
     280                 :            : /********************************************************************/
     281                 :            : int factor_add_primes = 0, factor_proven = 0;
     282                 :            : 
     283                 :            : /***********************************************************************/
     284                 :            : /**                                                                   **/
     285                 :            : /**                 FACTORIZATION (ECM) -- GN Jul-Aug 1998            **/
     286                 :            : /**   Integer factorization using the elliptic curves method (ECM).   **/
     287                 :            : /**   ellfacteur() returns a non trivial factor of N, assuming N>0,   **/
     288                 :            : /**   is composite, and has no prime divisor below 2^14 or so.        **/
     289                 :            : /**   Thanks to Paul Zimmermann for much helpful advice and to        **/
     290                 :            : /**   Guillaume Hanrot and Igor Schein for intensive testing          **/
     291                 :            : /**                                                                   **/
     292                 :            : /***********************************************************************/
     293                 :            : #define nbcmax 64                /* max number of simultaneous curves */
     294                 :            : #define bstpmax 1024                /* max number of baby step table entries */
     295                 :            : 
     296                 :            : /* addition/doubling/multiplication of a point on an 'elliptic curve mod N'
     297                 :            :  * may result in one of three things:
     298                 :            :  * - a new bona fide point
     299                 :            :  * - a point at infinity (denominator divisible by N)
     300                 :            :  * - a point at infinity mod some p | N but finite mod q | N betraying itself
     301                 :            :  *   by a denominator which has nontrivial gcd with N.
     302                 :            :  *
     303                 :            :  * In the second case, addition/doubling aborts, copying one of the summands
     304                 :            :  * to the destination array of points unless they coincide.
     305                 :            :  * Multiplication will stop at some unpredictable intermediate stage:  The
     306                 :            :  * destination will contain _some_ multiple of the input point, but not
     307                 :            :  * necessarily the desired one, which doesn't matter.  As long as we're
     308                 :            :  * multiplying (B1 phase) we simply carry on with the next multiplier.
     309                 :            :  * During the B2 phase, the only additions are the giant steps, and the
     310                 :            :  * worst that can happen here is that we lose one residue class mod 210
     311                 :            :  * of prime multipliers on 4 of the curves, so again, we ignore the problem
     312                 :            :  * and just carry on.)
     313                 :            :  *
     314                 :            :  * Idea: select nbc curves mod N and one point P on each of them. For each
     315                 :            :  * such P, compute [M]P = Q where M is the product of all powers <= B2 of
     316                 :            :  * primes <= nextprime(B1). Then check whether [p]Q for p < nextprime(B2)
     317                 :            :  * betrays a factor. This second stage looks separately at the primes in
     318                 :            :  * each residue class mod 210, four curves at a time, and steps additively
     319                 :            :  * to ever larger multipliers, by comparing X coordinates of points which we
     320                 :            :  * would need to add in order to reach another prime multiplier in the same
     321                 :            :  * residue class. 'Comparing' means that we accumulate a product of
     322                 :            :  * differences of X coordinates, and from time to time take a gcd of this
     323                 :            :  * product with N. Montgomery's multi-inverse trick is used heavily. */
     324                 :            : 
     325                 :            : /* *** auxiliary functions for ellfacteur: *** */
     326                 :            : /* (Rx,Ry) <- (Px,Py)+(Qx,Qy) over Z/NZ, z=1/(Px-Qx). If Ry = NULL, don't set */
     327                 :            : static void
     328                 :     331520 : FpE_add_i(GEN N, GEN z, GEN Px, GEN Py, GEN Qx, GEN Qy, GEN *Rx, GEN *Ry)
     329                 :            : {
     330                 :     331520 :   GEN slope = modii(mulii(subii(Py, Qy), z), N);
     331                 :     331520 :   GEN t = subii(sqri(slope), addii(Qx, Px));
     332                 :     331520 :   affii(modii(t, N), *Rx);
     333         [ +  + ]:     331520 :   if (Ry) {
     334                 :     326340 :     t = subii(mulii(slope, subii(Px, *Rx)), Py);
     335                 :     326340 :     affii(modii(t, N), *Ry);
     336                 :            :   }
     337                 :     331520 : }
     338                 :            : /* X -> Z; cannot add on one of the curves: make sure Z contains
     339                 :            :  * something useful before letting caller proceed */
     340                 :            : static void
     341                 :       6104 : ZV_aff(long n, GEN *X, GEN *Z)
     342                 :            : {
     343         [ +  + ]:       6104 :   if (X != Z) {
     344                 :            :     long k;
     345         [ +  + ]:      92666 :     for (k = n; k--; ) affii(X[k],Z[k]);
     346                 :            :   }
     347                 :       6104 : }
     348                 :            : 
     349                 :            : /* Parallel addition on nbc curves, assigning the result to locations at and
     350                 :            :  * following *X3, *Y3. (If Y-coords of result not desired, set Y=NULL.)
     351                 :            :  * Safe even if (X3,Y3) = (X2,Y2), _not_ if (X1,Y1). It is also safe to
     352                 :            :  * overwrite Y2 with X3. If nbc1 < nbc, the first summand is
     353                 :            :  * assumed to hold only nbc1 distinct points, repeated as often as we need
     354                 :            :  * them  (to add one point on each of a few curves to several other points on
     355                 :            :  * the same curves): only used with nbc1 = nbc or nbc1 = 4 | nbc.
     356                 :            :  *
     357                 :            :  * Return 0 [SUCCESS], 1 [N | den], 2 [gcd(den, N) is a factor of N, preserved
     358                 :            :  * in gl.
     359                 :            :  * Stack space is bounded by a constant multiple of lgefint(N)*nbc:
     360                 :            :  * - Phase 2 creates 12 items on the stack per iteration, of which 4 are twice
     361                 :            :  *   as long and 1 is thrice as long as N, i.e. 18 units per iteration.
     362                 :            :  * - Phase  1 creates 4 units.
     363                 :            :  * Total can be as large as 4*nbcmax + 18*8 units; ecm_elladd2() is
     364                 :            :  * just as bad, and elldouble() comes to 3*nbcmax + 29*8 units. */
     365                 :            : static int
     366                 :      37457 : ecm_elladd0(GEN N, GEN *gl, long nbc, long nbc1,
     367                 :            :             GEN *X1, GEN *Y1, GEN *X2, GEN *Y2, GEN *X3, GEN *Y3)
     368                 :            : {
     369         [ +  + ]:      37457 :   const ulong mask = (nbc1 == 4)? 3: ~0UL; /*nbc1 = 4 or nbc*/
     370                 :      37457 :   GEN W[2*nbcmax], *A = W+nbc; /* W[0],A[0] unused */
     371                 :            :   long i;
     372                 :      37457 :   pari_sp av = avma;
     373                 :            : 
     374                 :      37457 :   W[1] = subii(X1[0], X2[0]);
     375         [ +  + ]:     323512 :   for (i=1; i<nbc; i++)
     376                 :            :   { /*prepare for multi-inverse*/
     377                 :     286055 :     A[i] = subii(X1[i&mask], X2[i]); /* don't waste time reducing mod N */
     378                 :     286055 :     W[i+1] = modii(mulii(A[i], W[i]), N);
     379                 :            :   }
     380         [ +  + ]:      37457 :   if (!invmod(W[nbc], N, gl))
     381                 :            :   {
     382         [ +  + ]:        476 :     if (!equalii(N,*gl)) return 2;
     383                 :        462 :     ZV_aff(nbc, X2,X3);
     384         [ +  + ]:        462 :     if (Y3) ZV_aff(nbc, Y2,Y3);
     385                 :        462 :     avma = av; return 1;
     386                 :            :   }
     387                 :            : 
     388         [ +  - ]:     316064 :   while (i--) /* nbc times */
     389                 :            :   {
     390                 :     316064 :     pari_sp av2 = avma;
     391                 :     316064 :     GEN Px = X1[i&mask], Py = Y1[i&mask], Qx = X2[i], Qy = Y2[i];
     392         [ +  + ]:     316064 :     GEN z = i? mulii(*gl,W[i]): *gl; /*1/(Px-Qx)*/
     393         [ +  + ]:     316064 :     FpE_add_i(N,z,  Px,Py,Qx,Qy, X3+i, Y3? Y3+i: NULL);
     394         [ +  + ]:     316064 :     if (!i) break;
     395                 :     279083 :     avma = av2; *gl = modii(mulii(*gl, A[i]), N);
     396                 :            :   }
     397                 :      37457 :   avma = av; return 0;
     398                 :            : }
     399                 :            : 
     400                 :            : /* Shortcut, for use in cases where Y coordinates follow their corresponding
     401                 :            :  * X coordinates, and first summand doesn't need to be repeated */
     402                 :            : static int
     403                 :      36925 : ecm_elladd(GEN N, GEN *gl, long nbc, GEN *X1, GEN *X2, GEN *X3) {
     404                 :      36925 :   return ecm_elladd0(N, gl, nbc, nbc, X1, X1+nbc, X2, X2+nbc, X3, X3+nbc);
     405                 :            : }
     406                 :            : 
     407                 :            : /* As ecm_elladd except it does twice as many additions (and hides even more
     408                 :            :  * of the cost of the modular inverse); the net effect is the same as
     409                 :            :  * ecm_elladd(nbc,X1,X2,X3) && ecm_elladd(nbc,X4,X5,X6). Safe to
     410                 :            :  * have X2=X3, X5=X6, or X1,X2 coincide with X4,X5 in any order. */
     411                 :            : static int
     412                 :        980 : ecm_elladd2(GEN N, GEN *gl, long nbc,
     413                 :            :             GEN *X1, GEN *X2, GEN *X3, GEN *X4, GEN *X5, GEN *X6)
     414                 :            : {
     415                 :        980 :   GEN *Y1 = X1+nbc, *Y2 = X2+nbc, *Y3 = X3+nbc;
     416                 :        980 :   GEN *Y4 = X4+nbc, *Y5 = X5+nbc, *Y6 = X6+nbc;
     417                 :        980 :   GEN W[4*nbcmax], *A = W+2*nbc; /* W[0],A[0] unused */
     418                 :            :   long i, j;
     419                 :        980 :   pari_sp av = avma;
     420                 :            : 
     421                 :        980 :   W[1] = subii(X1[0], X2[0]);
     422         [ +  + ]:       7840 :   for (i=1; i<nbc; i++)
     423                 :            :   {
     424                 :       6860 :     A[i] = subii(X1[i], X2[i]); /* don't waste time reducing mod N here */
     425                 :       6860 :     W[i+1] = modii(mulii(A[i], W[i]), N);
     426                 :            :   }
     427         [ +  + ]:       8820 :   for (j=0; j<nbc; i++,j++)
     428                 :            :   {
     429                 :       7840 :     A[i] = subii(X4[j], X5[j]);
     430                 :       7840 :     W[i+1] = modii(mulii(A[i], W[i]), N);
     431                 :            :   }
     432         [ +  + ]:        980 :   if (!invmod(W[2*nbc], N, gl))
     433                 :            :   {
     434         [ -  + ]:         14 :     if (!equalii(N,*gl)) return 2;
     435                 :         14 :     ZV_aff(2*nbc, X2,X3); /* hack: 2*nbc => copy Y2->Y3 */
     436                 :         14 :     ZV_aff(2*nbc, X5,X6); /* also copy Y5->Y6 */
     437                 :         14 :     avma = av; return 1;
     438                 :            :   }
     439                 :            : 
     440         [ +  + ]:       8694 :   while (j--) /* nbc times */
     441                 :            :   {
     442                 :       7728 :     pari_sp av2 = avma;
     443                 :       7728 :     GEN Px = X4[j], Py = Y4[j], Qx = X5[j], Qy = Y5[j];
     444                 :       7728 :     GEN z = mulii(*gl,W[--i]); /*1/(Px-Qx)*/
     445                 :       7728 :     FpE_add_i(N,z, Px,Py, Qx,Qy, X6+j,Y6+j);
     446                 :       7728 :     avma = av2; *gl = modii(mulii(*gl, A[i]), N);
     447                 :            :   }
     448         [ +  - ]:       7728 :   while (i--) /* nbc times */
     449                 :            :   {
     450                 :       7728 :     pari_sp av2 = avma;
     451                 :       7728 :     GEN Px = X1[i], Py = Y1[i], Qx = X2[i], Qy = Y2[i];
     452         [ +  + ]:       7728 :     GEN z = i? mulii(*gl, W[i]): *gl; /*1/(Px-Qx)*/
     453                 :       7728 :     FpE_add_i(N,z, Px,Py, Qx,Qy, X3+i,Y3+i);
     454         [ +  + ]:       7728 :     if (!i) break;
     455                 :       6762 :     avma = av2; *gl = modii(mulii(*gl, A[i]), N);
     456                 :            :   }
     457                 :        980 :   avma = av; return 0;
     458                 :            : }
     459                 :            : 
     460                 :            : /* Parallel doubling on nbc curves, assigning the result to locations at
     461                 :            :  * and following *X2.  Safe to be called with X2 equal to X1.  Return
     462                 :            :  * value as for ecm_elladd.  If we find a point at infinity mod N,
     463                 :            :  * and if X1 != X2, we copy the points at X1 to X2. */
     464                 :            : static int
     465                 :       7714 : elldouble(GEN N, GEN *gl, long nbc, GEN *X1, GEN *X2)
     466                 :            : {
     467                 :       7714 :   GEN *Y1 = X1+nbc, *Y2 = X2+nbc;
     468                 :            :   GEN W[nbcmax+1]; /* W[0] unused */
     469                 :            :   long i;
     470                 :       7714 :   pari_sp av = avma;
     471                 :       7714 :   /*W[0] = gen_1;*/ W[1] = Y1[0];
     472         [ +  + ]:      61712 :   for (i=1; i<nbc; i++) W[i+1] = modii(mulii(Y1[i], W[i]), N);
     473         [ -  + ]:       7714 :   if (!invmod(W[nbc], N, gl))
     474                 :            :   {
     475         [ #  # ]:          0 :     if (!equalii(N,*gl)) return 2;
     476                 :          0 :     ZV_aff(2*nbc,X1,X2); /* also copies Y1->Y2 */
     477                 :          0 :     avma = av; return 1;
     478                 :            :   }
     479         [ +  + ]:      69426 :   while (i--) /* nbc times */
     480                 :            :   {
     481                 :            :     pari_sp av2;
     482         [ +  + ]:      61712 :     GEN v, w, L, z = i? mulii(*gl,W[i]): *gl;
     483         [ +  + ]:      61712 :     if (i) *gl = modii(mulii(*gl, Y1[i]), N);
     484                 :      61712 :     av2 = avma;
     485                 :      61712 :     L = modii(mulii(addsi(1, mului(3, Fp_sqr(X1[i],N))), z), N);
     486         [ +  - ]:      61712 :     if (signe(L)) /* half of zero is still zero */
     487         [ +  + ]:      61712 :       L = shifti(mod2(L)? addii(L, N): L, -1);
     488                 :      61712 :     v = modii(subii(sqri(L), shifti(X1[i],1)), N);
     489                 :      61712 :     w = modii(subii(mulii(L, subii(X1[i], v)), Y1[i]), N);
     490                 :      61712 :     affii(v, X2[i]);
     491                 :      61712 :     affii(w, Y2[i]);
     492                 :      61712 :     avma = av2;
     493                 :            :   }
     494                 :       7714 :   avma = av; return 0;
     495                 :            : }
     496                 :            : 
     497                 :            : /* Parallel multiplication by an odd prime k on nbc curves, storing the
     498                 :            :  * result to locations at and following *X2. Safe to be called with X2 = X1.
     499                 :            :  * Return values as ecm_elladd. Uses (a simplified variant of) Montgomery's
     500                 :            :  * PRAC algorithm; see ftp://ftp.cwi.nl/pub/pmontgom/Lucas.ps.gz .
     501                 :            :  * With thanks to Paul Zimmermann for the reference.  --GN1998Aug13 */
     502                 :            : static int
     503                 :      30513 : get_rule(ulong d, ulong e)
     504                 :            : {
     505         [ +  + ]:      30513 :   if (d <= e + (e>>2)) /* floor(1.25*e) */
     506                 :            :   {
     507         [ +  + ]:       2359 :     if ((d+e)%3 == 0) return 0; /* rule 1 */
     508         [ -  + ]:       1421 :     if ((d-e)%6 == 0) return 1;  /* rule 2 */
     509                 :            :   }
     510                 :            :   /* d <= 4*e but no ofl */
     511         [ +  + ]:      29575 :   if ((d+3)>>2 <= e) return 2; /* rule 3, common case */
     512         [ +  + ]:       1582 :   if ((d&1)==(e&1))  return 1; /* rule 4 = rule 2 */
     513         [ +  + ]:        700 :   if (!(d&1))        return 3; /* rule 5 */
     514         [ +  + ]:        217 :   if (d%3 == 0)      return 4; /* rule 6 */
     515         [ +  - ]:         21 :   if ((d+e)%3 == 0)  return 5; /* rule 7 */
     516         [ #  # ]:          0 :   if ((d-e)%3 == 0)  return 6; /* rule 8 */
     517                 :            :   /* when we get here, e is even, otherwise one of rules 4,5 would apply */
     518                 :      30513 :   return 7; /* rule 9 */
     519                 :            : }
     520                 :            : 
     521                 :            : /* k>2 assumed prime, XAUX = scratchpad */
     522                 :            : static int
     523                 :       5166 : ellmult(GEN N, GEN *gl, long nbc, ulong k, GEN *X1, GEN *X2, GEN *XAUX)
     524                 :            : {
     525                 :            :   ulong r, d, e, e1;
     526                 :            :   int res;
     527                 :       5166 :   GEN *A = X2, *B = XAUX, *T = XAUX + 2*nbc;
     528                 :            : 
     529                 :       5166 :   ZV_aff(2*nbc,X1,XAUX);
     530                 :            :   /* first doubling picks up X1;  after this we'll be working in XAUX and
     531                 :            :    * X2 only, mostly via A and B and T */
     532         [ -  + ]:       5166 :   if ((res = elldouble(N, gl, nbc, X1, X2)) != 0) return res;
     533                 :            : 
     534                 :            :   /* split the work at the golden ratio */
     535                 :       5166 :   r = (ulong)(k*0.61803398875 + .5);
     536                 :       5166 :   d = k - r;
     537                 :       5166 :   e = r - d; /* d+e == r, so no danger of ofl below */
     538         [ +  + ]:      35483 :   while (d != e)
     539                 :            :   { /* apply one of the nine transformations from PM's Table 4. */
     540   [ +  +  +  +  :      30513 :     switch(get_rule(d,e))
             +  +  -  -  
                      - ]
     541                 :            :     {
     542                 :            :     case 0: /* rule 1 */
     543         [ +  + ]:        938 :       if ( (res = ecm_elladd(N, gl, nbc, A, B, T)) ) return res;
     544         [ +  + ]:        924 :       if ( (res = ecm_elladd2(N, gl, nbc, T, A, A, T, B, B)) != 0) return res;
     545                 :        910 :       e1 = d - e; d = (d + e1)/3; e = (e - e1)/3; break;
     546                 :            :     case 1: /* rules 2 and 4 */
     547         [ +  + ]:        882 :       if ( (res = ecm_elladd(N, gl, nbc, A, B, B)) ) return res;
     548         [ -  + ]:        861 :       if ( (res = elldouble(N, gl, nbc, A, A)) ) return res;
     549                 :        861 :       d = (d-e)>>1; break;
     550                 :            :     case 3: /* rule 5 */
     551         [ -  + ]:        483 :       if ( (res = elldouble(N, gl, nbc, A, A)) ) return res;
     552                 :        483 :       d >>= 1; break;
     553                 :            :     case 4: /* rule 6 */
     554         [ -  + ]:        196 :       if ( (res = elldouble(N, gl, nbc, A, T)) ) return res;
     555         [ -  + ]:        196 :       if ( (res = ecm_elladd(N, gl, nbc, T, A, A)) ) return res;
     556         [ -  + ]:        196 :       if ( (res = ecm_elladd(N, gl, nbc, A, B, B)) ) return res;
     557                 :        196 :       d = d/3 - e; break;
     558                 :            :     case 2: /* rule 3 */
     559         [ +  + ]:      27993 :       if ( (res = ecm_elladd(N, gl, nbc, A, B, B)) ) return res;
     560                 :      27846 :       d -= e; break;
     561                 :            :     case 5: /* rule 7 */
     562         [ -  + ]:         21 :       if ( (res = elldouble(N, gl, nbc, A, T)) ) return res;
     563         [ -  + ]:         21 :       if ( (res = ecm_elladd2(N, gl, nbc, T, A, A, T, B, B)) != 0) return res;
     564                 :         21 :       d = (d - 2*e)/3; break;
     565                 :            :     case 6: /* rule 8 */
     566         [ #  # ]:          0 :       if ( (res = ecm_elladd(N, gl, nbc, A, B, B)) ) return res;
     567         [ #  # ]:          0 :       if ( (res = elldouble(N, gl, nbc, A, T)) ) return res;
     568         [ #  # ]:          0 :       if ( (res = ecm_elladd(N, gl, nbc, T, A, A)) ) return res;
     569                 :          0 :       d = (d - e)/3; break;
     570                 :            :     case 7: /* rule 9 */
     571         [ #  # ]:          0 :       if ( (res = elldouble(N, gl, nbc, B, B)) ) return res;
     572                 :          0 :       e >>= 1; break;
     573                 :            :     }
     574                 :            :     /* swap d <-> e and A <-> B if necessary */
     575         [ +  + ]:      30317 :     if (d < e) { lswap(d,e); pswap(A,B); }
     576                 :            :   }
     577                 :       5166 :   return ecm_elladd(N, gl, nbc, XAUX, X2, X2);
     578                 :            : }
     579                 :            : 
     580                 :            : /* Auxiliary routines need < (3*nbc+240)*tf words on the PARI stack, in
     581                 :            :  * addition to the spc*(tf+1) words occupied by our main table.
     582                 :            :  * If stack space is already tight, use the heap & newblock(). */
     583                 :            : static GEN*
     584                 :         35 : alloc_scratch(long nbc, long spc, long tf)
     585                 :            : {
     586                 :         35 :   pari_sp bot = pari_mainstack->bot;
     587                 :         35 :   long i, tw = evallg(tf) | evaltyp(t_INT), len = spc + 385 + spc*tf;
     588                 :            :   GEN *X, w;
     589         [ -  + ]:         35 :   if ((long)((GEN)avma - (GEN)bot) < len + (3*nbc + 240)*tf)
     590                 :            :   {
     591         [ #  # ]:          0 :     if (DEBUGLEVEL>4) err_printf("ECM: stack tight, using heap space\n");
     592                 :          0 :     X = (GEN*)newblock(len);
     593                 :            :   } else
     594                 :         35 :     X = (GEN*)new_chunk(len);
     595                 :            :   /* hack for X[i] = cgeti(tf). X = current point in B1 phase */
     596                 :         35 :   w = (GEN)(X + spc + 385);
     597         [ +  + ]:     177555 :   for (i = spc-1; i >= 0; i--) { X[i] = w; *w = tw; w += tf; }
     598                 :         35 :   return X;
     599                 :            : }
     600                 :            : 
     601                 :            : /* PRAC implementation notes - main changes against the paper version:
     602                 :            :  * (1) The general function [m+n]P = f([m]P,[n]P,[m-n]P) collapses (for m!=n)
     603                 :            :  * to an ecm_elladd() which does not depend on the third argument; thus
     604                 :            :  * references to the third variable (C in the paper) can be eliminated.
     605                 :            :  * (2) Since our multipliers are prime, the outer loop of the paper
     606                 :            :  * version executes only once, and thus is invisible above.
     607                 :            :  * (3) The first step in the inner loop of the paper version will always be
     608                 :            :  * rule 3, but the addition requested by this rule amounts to a doubling, and
     609                 :            :  * will always be followed by a swap, so we have unrolled this first iteration.
     610                 :            :  * (4) Simplifications in rules 6 and 7 are possible given the above, and we
     611                 :            :  * save one addition in each of the two cases.  NB none of the other
     612                 :            :  * ecm_elladd()s in the loop can ever degenerate into an elldouble.
     613                 :            :  * (5) I tried to optimize for rule 3, which is used more frequently than all
     614                 :            :  * others together, but it didn't improve things, so I removed the nested
     615                 :            :  * tight loop again.  --GN */
     616                 :            : 
     617                 :            : /* The main loop body of ellfacteur() runs _slower_ under PRAC than under a
     618                 :            :  * straightforward left-shift binary multiplication when N has <30 digits and
     619                 :            :  * B1 is small;  PRAC wins when N and B1 get larger.  Weird. --GN */
     620                 :            : 
     621                 :            : /* memory layout in ellfacteur():  a large array of GEN pointers, and one
     622                 :            :  * huge chunk of memory containing all the actual GEN (t_INT) objects.
     623                 :            :  * nbc is constant throughout the invocation:
     624                 :            :  * - The B1 stage of each iteration through the main loop needs little
     625                 :            :  * space:  enough for the X and Y coordinates of the current points,
     626                 :            :  * and twice as much again as scratchpad for ellmult().
     627                 :            :  * - The B2 stage, starting from some current set of points Q, needs, in
     628                 :            :  * succession:
     629                 :            :  *   + space for [2]Q, [4]Q, ..., [10]Q, and [p]Q for building the helix;
     630                 :            :  *   + space for 48*nbc X and Y coordinates to hold the helix.  This could
     631                 :            :  *   re-use [2]Q,...,[8]Q, but only with difficulty, since we don't
     632                 :            :  *   know in advance which residue class mod 210 our p is going to be in.
     633                 :            :  *   It can and should re-use [p]Q, though;
     634                 :            :  *   + space for (temporarily [30]Q and then) [210]Q, [420]Q, and several
     635                 :            :  *   further doublings until the giant step multiplier is reached.  This
     636                 :            :  *   can re-use the remaining cells from above.  The computation of [210]Q
     637                 :            :  *   will have been the last call to ellmult() within this iteration of the
     638                 :            :  *   main loop, so the scratchpad is now also free to be re-used. We also
     639                 :            :  *   compute [630]Q by a parallel addition;  we'll need it later to get the
     640                 :            :  *   baby-step table bootstrapped a little faster.
     641                 :            :  *   + Finally, for no more than 4 curves at a time, room for up to 1024 X
     642                 :            :  *   coordinates only: the Y coordinates needed whilst setting up this baby
     643                 :            :  *   step table are temporarily stored in the upper half, and overwritten
     644                 :            :  *   during the last series of additions.
     645                 :            :  *
     646                 :            :  * Graphically:  after end of B1 stage (X,Y are the coords of Q):
     647                 :            :  * +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--
     648                 :            :  * | X Y |  scratch  | [2]Q| [4]Q| [6]Q| [8]Q|[10]Q|    ...    | ...
     649                 :            :  * +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--
     650                 :            :  * *X    *XAUX *XT   *XD                                       *XB
     651                 :            :  *
     652                 :            :  * [30]Q is computed from [10]Q.  [210]Q can go into XY, etc:
     653                 :            :  * +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--
     654                 :            :  * |[210]|[420]|[630]|[840]|[1680,3360,6720,...,2048*210]      |bstp table...
     655                 :            :  * +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--
     656                 :            :  * *X    *XAUX *XT   *XD      [*XG, somewhere here]            *XB .... *XH
     657                 :            :  *
     658                 :            :  * So we need (13 + 48) * 2 * nbc slots here + 4096 slots for the baby step
     659                 :            :  * table (not all of which will be used when we start with a small B1, but
     660                 :            :  * better to allocate and initialize ahead of time all the slots that might
     661                 :            :  * be needed later).
     662                 :            :  *
     663                 :            :  * Note on memory locality:  During the B2 phase, accesses to the helix
     664                 :            :  * (once it is set up) will be clustered by curves (4 out of nbc at a time).
     665                 :            :  * Accesses to the baby steps table will wander from one end of the array to
     666                 :            :  * the other and back, one such cycle per giant step, and during a full cycle
     667                 :            :  * we would expect on the order of 2E4 accesses when using the largest giant
     668                 :            :  * step size.  Thus we shouldn't be doing too bad with respect to thrashing
     669                 :            :  * a 512KBy L2 cache.  However, we don't want the baby step table to grow
     670                 :            :  * larger than this, even if it would reduce the number of EC operations by a
     671                 :            :  * few more per cent for very large B2, lest cache thrashing slow down
     672                 :            :  * everything disproportionally. --GN */
     673                 :            : 
     674                 :            : /* parameters for MR_Jaeschke() via snextpr(), for use by ellfacteur() */
     675                 :            : static const long MR_Jaeschke_k1 = 16;/* B1 phase, foolproof below 10^12 */
     676                 :            : static const long MR_Jaeschke_k2 = 1; /* B2 phase, not foolproof, 2xfaster */
     677                 :            : /* MR_Jaeschke_k2 will let thousands of composites slip through, which doesn't
     678                 :            :  * harm ECM, but ellmult() during the B1 phase should only be fed primes
     679                 :            :  * which really are prime */
     680                 :            : /* ellfacteur() has been re-tuned to be useful as a first stage before
     681                 :            :  * MPQS, especially for large arguments, when 'insist' is false, and now
     682                 :            :  * also for the case when 'insist' is true, vaguely following suggestions
     683                 :            :  * by Paul Zimmermann (http://www.loria.fr/~zimmerma/records/ecmnet.html).
     684                 :            :  * --GN 1998Jul,Aug */
     685                 :            : GEN
     686                 :        574 : ellfacteur(GEN N, int insist)
     687                 :            : {
     688                 :        574 :   const ulong TB1[] = {
     689                 :            :     142,172,208,252,305,370,450,545,661,801,972,1180,1430,
     690                 :            :     1735,2100,2550,3090,3745,4540,5505,6675,8090,9810,11900,
     691                 :            :     14420,17490,21200,25700,31160,37780UL,45810UL,55550UL,67350UL,
     692                 :            :     81660UL,99010UL,120050UL,145550UL,176475UL,213970UL,259430UL,
     693                 :            :     314550UL,381380UL,462415UL,560660UL,679780UL,824220UL,999340UL,
     694                 :            :     1211670UL,1469110UL,1781250UL,2159700UL,2618600UL,3175000UL,
     695                 :            :     3849600UL,4667500UL,5659200UL,6861600UL,8319500UL,10087100UL,
     696                 :            :     12230300UL,14828900UL,17979600UL,21799700UL,26431500UL,
     697                 :            :     32047300UL,38856400UL, /* 110 times that still fits into 32bits */
     698                 :            : #ifdef LONG_IS_64BIT
     699                 :            :     47112200UL,57122100UL,69258800UL,83974200UL,101816200UL,
     700                 :            :     123449000UL,149678200UL,181480300UL,220039400UL,266791100UL,
     701                 :            :     323476100UL,392204900UL,475536500UL,576573500UL,699077800UL,
     702                 :            :     847610500UL,1027701900UL,1246057200UL,1510806400UL,1831806700UL,
     703                 :            :     2221009800UL,2692906700UL,3265067200UL,3958794400UL,4799917500UL,
     704                 :            :     /* Someone can extend this table when the hardware gets faster */
     705                 :            : #endif
     706                 :            :     };
     707                 :        574 :   const ulong TB1_for_stage[] = {
     708                 :            :    /* Start below the optimal B1 for finding factors which would just have been
     709                 :            :     * missed by pollardbrent(), and escalate, changing curves to give good
     710                 :            :     * coverage of the small factor ranges. Entries grow faster than what would
     711                 :            :     * be optimal but a table instead of a 2D array keeps the code simple */
     712                 :            :     500,520,560,620,700,800,900,1000,1150,1300,1450,1600,1800,2000,
     713                 :            :     2200,2450,2700,2950,3250,3600,4000,4400,4850,5300,5800,6400,
     714                 :            :     7100,7850,8700,9600,10600,11700,12900,14200,15700,17300,
     715                 :            :     19000,21000,23200,25500,28000,31000,34500UL,38500UL,43000UL,
     716                 :            :     48000UL,53800UL,60400UL,67750UL,76000UL,85300UL,95700UL,
     717                 :            :     107400UL,120500UL,135400UL,152000UL,170800UL,191800UL,215400UL,
     718                 :            :     241800UL,271400UL,304500UL,341500UL,383100UL,429700UL,481900UL,
     719                 :            :     540400UL,606000UL,679500UL,761800UL,854100UL,957500UL,1073500UL,
     720                 :            :   };
     721                 :            :   long nbc,nbc2,dsn,dsnmax,rep,spc,gse,gss,rcn,rcn0,bstp,bstp0;
     722                 :        574 :   long a, i, j, k, size = expi(N) + 1, tf = lgefint(N);
     723                 :            :   ulong B1,B2,B2_p,B2_rt,m,p,p0,dp;
     724                 :            :   GEN *X,*XAUX,*XT,*XD,*XG,*YG,*XH,*XB,*XB2,*Xh,*Yh,*Xb;
     725                 :        574 :   GEN res = cgeti(tf), gl;
     726                 :        574 :   pari_sp av1, avtmp, av = avma;
     727                 :            :   pari_timer T;
     728                 :            :   int rflag;
     729                 :            : 
     730                 :            :   /* Determine where to start, how long to persist, and how many curves to
     731                 :            :    * use in parallel */
     732         [ +  + ]:        574 :   if (insist)
     733                 :            :   {
     734                 :            : #ifdef LONG_IS_64BIT
     735                 :         12 :     const long DSNMAX = 90;
     736                 :            : #else
     737                 :          2 :     const long DSNMAX = 65;
     738                 :            : #endif
     739                 :         14 :     dsnmax = (size >> 2) - 10;
     740         [ +  - ]:         14 :     if (dsnmax < 0) dsnmax = 0;
     741         [ #  # ]:          0 :     else if (dsnmax > DSNMAX) dsnmax = DSNMAX;
     742                 :         14 :     dsn = (size >> 3) - 5;
     743         [ +  - ]:         14 :     if (dsn < 0) dsn = 0;
     744         [ #  # ]:          0 :     else if (dsn > 47) dsn = 47;
     745                 :            :     /* pick up the torch where non-insistent stage would have given up */
     746                 :         14 :     nbc = dsn + (dsn >> 2) + 9; /* 8 or more curves in parallel */
     747                 :         14 :     nbc &= ~3; /* 4 | nbc */
     748         [ -  + ]:         14 :     if (nbc > nbcmax) nbc = nbcmax;
     749                 :         14 :     a = 1 + (nbcmax<<7)*(size&0xffff); /* seed for choice of curves */
     750                 :         14 :     rep = 0; /* gcc -Wall */
     751                 :            :   }
     752                 :            :   else
     753                 :            :   {
     754                 :        560 :     dsn = (size - 140) >> 3;
     755         [ -  + ]:        560 :     if (dsn > 12) dsn = 12;
     756                 :        560 :     dsnmax = 72;
     757         [ +  + ]:        560 :     if (dsn < 0) /* < 140 bits: decline the task */
     758                 :            :     {
     759                 :            : #ifdef __EMX__
     760                 :            :       /* unless DOS/EMX: MPQS's disk access is abysmally slow */
     761                 :            :       dsn = 0; rep = 20; nbc = 8;
     762                 :            : #else
     763         [ -  + ]:        539 :       if (DEBUGLEVEL >= 4)
     764                 :          0 :         err_printf("ECM: number too small to justify this stage\n");
     765                 :        539 :       avma = av; return NULL;
     766                 :            : #endif
     767                 :            :     }
     768                 :            :     else
     769                 :            :     {
     770                 :         21 :       rep = (size <= 248 ?
     771 [ +  - ][ +  - ]:         21 :              (size <= 176 ? (size - 124) >> 4 : (size - 148) >> 3) :
     772                 :          0 :              (size - 224) >> 1);
     773                 :         21 :       nbc = ((size >> 3) << 2) - 80;
     774         [ +  - ]:         21 :       if (nbc < 8) nbc = 8;
     775         [ #  # ]:          0 :       else if (nbc > nbcmax) nbc = nbcmax;
     776                 :            : #ifdef __EMX__
     777                 :            :       rep += 20;
     778                 :            : #endif
     779                 :            :     }
     780                 :            : 
     781                 :            :     /* Use disjoint sets of curves for non-insist and insist phases; moreover,
     782                 :            :      * repeated calls acting on factors of the same original number should try
     783                 :            :      * to use fresh curves. The following achieves this */
     784                 :         21 :     a = 1 + (nbcmax<<3)*(size & 0xf);
     785                 :            :   }
     786         [ -  + ]:         35 :   if (dsn > dsnmax) dsn = dsnmax;
     787                 :            : 
     788         [ -  + ]:         35 :   if (DEBUGLEVEL >= 4)
     789                 :            :   {
     790                 :          0 :     timer_start(&T);
     791                 :          0 :     err_printf("ECM: working on %ld curves at a time; initializing", nbc);
     792         [ #  # ]:          0 :     if (!insist)
     793                 :            :     {
     794         [ #  # ]:          0 :       if (rep == 1) err_printf(" for one round");
     795                 :          0 :       else          err_printf(" for up to %ld rounds", rep);
     796                 :            :     }
     797                 :          0 :     err_printf("...\n");
     798                 :            :   }
     799                 :            : 
     800                 :         35 :   nbc2 = nbc << 1;
     801                 :         35 :   spc = (13 + 48) * nbc2 + bstpmax * 4;
     802                 :         35 :   X = alloc_scratch(nbc, spc, tf);
     803                 :         35 :   XAUX = X    + nbc2; /* scratchpad for ellmult() */
     804                 :         35 :   XT   = XAUX + nbc2; /* ditto, will later hold [3*210]Q */
     805                 :         35 :   XD   = XT   + nbc2; /* room for various multiples */
     806                 :         35 :   XB   = XD   + 10*nbc2; /* start of baby steps table */
     807                 :         35 :   XB2  = XB   + 2 * bstpmax; /* middle of baby steps table */
     808                 :         35 :   XH   = XB2  + 2 * bstpmax; /* end of bstps table, start of helix */
     809                 :         35 :   Xh   = XH   + 48*nbc2; /* little helix, X coords */
     810                 :         35 :   Yh   = XH   + 192;     /* ditto, Y coords */
     811                 :            :   /* XG will be set inside the main loop, since it depends on B2 */
     812                 :            : 
     813                 :            :   /* Xh range of 384 pointers not set; these will later duplicate the pointers
     814                 :            :    * in the XH range, 4 curves at a time. Some of the cells reserved here for
     815                 :            :    * the XB range will never be used, instead, we'll warp the pointers to
     816                 :            :    * connect to (read-only) GENs in the X/XD range */
     817                 :            :   for(;;)
     818                 :            :   {
     819                 :         49 :     byteptr d0, d = diffptr;
     820                 :            : 
     821                 :         49 :     rcn = NPRC; /* multipliers begin at the beginning */
     822                 :            :     /* pick curves & bounds */
     823         [ +  + ]:        833 :     for (i = nbc2; i--; ) affui(a++, X[i]);
     824         [ +  + ]:         49 :     B1 = insist ? TB1[dsn] : TB1_for_stage[dsn];
     825                 :         49 :     B2 = 110*B1;
     826                 :         49 :     B2_rt = (ulong)(sqrt((double)B2));
     827                 :            :     /* pick giant step exponent and size.
     828                 :            :      * With 32 baby steps, a giant step corresponds to 32*420 = 13440,
     829                 :            :      * appropriate for the smallest B2s. With 1024, a giant step will be 430080;
     830                 :            :      * appropriate for B1 >~ 42000, where 512 baby steps would imply roughly
     831                 :            :      * the same number of E.C. additions. */
     832 [ +  - ][ +  + ]:         49 :     gse = B1 < 656
         [ #  # ][ #  # ]
                 [ #  # ]
     833                 :            :             ? (B1 < 200? 5: 6)
     834                 :            :             : (B1 < 10500
     835                 :            :               ? (B1 < 2625? 7: 8)
     836                 :            :               : (B1 < 42000? 9: 10));
     837                 :         49 :     gss = 1UL << gse;
     838                 :         49 :     XG = XT + gse*nbc2; /* will later hold [2^(gse+1)*210]Q */
     839                 :         49 :     YG = XG + nbc;
     840                 :            : 
     841         [ -  + ]:         49 :     if (DEBUGLEVEL >= 4) {
     842                 :          0 :       err_printf("ECM: time = %6ld ms\nECM: dsn = %2ld,\tB1 = %4lu,",
     843                 :            :                  timer_delay(&T), dsn, B1);
     844                 :          0 :       err_printf("\tB2 = %6lu,\tgss = %4ld*420\n", B2, gss);
     845                 :            :     }
     846                 :         49 :     p = 0;
     847                 :         49 :     NEXT_PRIME_VIADIFF(p,d);
     848                 :            : 
     849                 :            :     /* ---B1 PHASE--- */
     850                 :            :     /* treat p=2 separately */
     851                 :         49 :     B2_p = B2 >> 1;
     852         [ +  + ]:        735 :     for (m=1; m<=B2_p; m<<=1)
     853                 :            :     {
     854         [ +  - ]:        686 :       if ((rflag = elldouble(N, &gl, nbc, X, X)) > 1) goto fin;
     855         [ -  + ]:        686 :       else if (rflag) break;
     856                 :            :     }
     857                 :            :     /* p=3,...,nextprime(B1) */
     858 [ +  - ][ +  + ]:       1687 :     while (p < B1 && p <= B2_rt)
     859                 :            :     {
     860                 :       1652 :       pari_sp av = avma;
     861                 :       1652 :       p = snextpr(p, &d, &rcn, NULL, MR_Jaeschke_k1);
     862                 :       1652 :       B2_p = B2/p; /* beware integer overflow on 32-bit CPUs */
     863         [ +  + ]:       5446 :       for (m=1; m<=B2_p; m*=p)
     864                 :            :       {
     865         [ +  + ]:       3976 :         if ((rflag = ellmult(N, &gl, nbc, p, X, X, XAUX)) > 1) goto fin;
     866         [ +  + ]:       3962 :         else if (rflag) break;
     867                 :       3794 :         avma = av;
     868                 :            :       }
     869                 :       1638 :       avma = av;
     870                 :            :     }
     871                 :            :     /* primes p larger than sqrt(B2) appear only to the 1st power */
     872         [ +  + ]:       1120 :     while (p < B1)
     873                 :            :     {
     874                 :       1085 :       pari_sp av = avma;
     875                 :       1085 :       p = snextpr(p, &d, &rcn, NULL, MR_Jaeschke_k1);
     876         [ +  - ]:       1085 :       if (ellmult(N, &gl, nbc, p, X, X, XAUX) > 1) goto fin;
     877                 :       1085 :       avma = av;
     878                 :            :     }
     879         [ -  + ]:         35 :     if (DEBUGLEVEL >= 4) {
     880                 :          0 :       err_printf("ECM: time = %6ld ms, B1 phase done, ", timer_delay(&T));
     881                 :          0 :       err_printf("p = %lu, setting up for B2\n", p);
     882                 :            :     }
     883                 :            : 
     884                 :            :     /* ---B2 PHASE--- */
     885                 :            :     /* compute [2]Q,...,[10]Q, needed to build the helix */
     886         [ +  - ]:         35 :     if (elldouble(N, &gl, nbc, X, XD) > 1) goto fin; /* [2]Q */
     887         [ +  - ]:         35 :     if (elldouble(N, &gl, nbc, XD, XD + nbc2) > 1) goto fin; /* [4]Q */
     888         [ +  - ]:         35 :     if (ecm_elladd(N, &gl, nbc, XD, XD + nbc2, XD + (nbc<<2)) > 1)
     889                 :            :       goto fin; /* [6]Q */
     890         [ +  - ]:         35 :     if (ecm_elladd2(N, &gl, nbc,
     891                 :        140 :                 XD, XD + (nbc<<2), XT + (nbc<<3),
     892                 :        175 :                 XD + nbc2, XD + (nbc<<2), XD + (nbc<<3)) > 1)
     893                 :            :       goto fin; /* [8]Q and [10]Q */
     894         [ -  + ]:         35 :     if (DEBUGLEVEL >= 7) err_printf("\t(got [2]Q...[10]Q)\n");
     895                 :            : 
     896                 :            :     /* get next prime (still using the foolproof test) */
     897                 :         35 :     p = snextpr(p, &d, &rcn, NULL, MR_Jaeschke_k1);
     898                 :            :     /* make sure we have the residue class number (mod 210) */
     899         [ +  - ]:         35 :     if (rcn == NPRC)
     900                 :            :     {
     901                 :         35 :       rcn = prc210_no[(p % 210) >> 1];
     902         [ -  + ]:         35 :       if (rcn == NPRC)
     903                 :            :       {
     904                 :          0 :         err_printf("ECM: %lu should have been prime but isn\'t\n", p);
     905                 :          0 :         pari_err_BUG("ellfacteur");
     906                 :            :       }
     907                 :            :     }
     908                 :            : 
     909                 :            :     /* compute [p]Q and put it into its place in the helix */
     910         [ +  - ]:         35 :     if (ellmult(N, &gl, nbc, p, X, XH + rcn*nbc2, XAUX) > 1) goto fin;
     911         [ -  + ]:         35 :     if (DEBUGLEVEL >= 7)
     912                 :          0 :       err_printf("\t(got [p]Q, p = %lu = prc210_rp[%ld] mod 210)\n", p, rcn);
     913                 :            : 
     914                 :            :     /* save current p, d, and rcn;  we'll need them more than once below */
     915                 :         35 :     p0 = p;
     916                 :         35 :     d0 = d;
     917                 :         35 :     rcn0 = rcn; /* remember where the helix wraps */
     918                 :         35 :     bstp0 = 0; /* p is at baby-step offset 0 from itself */
     919                 :            : 
     920                 :            :     /* fill up the helix, stepping forward through the prime residue classes
     921                 :            :      * mod 210 until we're back at the r'class of p0.  Keep updating p so
     922                 :            :      * that we can print meaningful diagnostics if a factor shows up; don't
     923                 :            :      * bother checking which of these p's are in fact prime */
     924         [ +  + ]:       1680 :     for (i = 47; i; i--) /* 47 iterations */
     925                 :            :     {
     926                 :       1645 :       p += (dp = (ulong)prc210_d1[rcn]);
     927         [ +  + ]:       1645 :       if (rcn == 47)
     928                 :            :       { /* wrap mod 210 */
     929         [ +  - ]:         35 :         if (ecm_elladd(N, &gl, nbc, XT+dp*nbc, XH+rcn*nbc2, XH) > 1) goto fin;
     930                 :         35 :         rcn = 0; continue;
     931                 :            :       }
     932         [ +  - ]:       1610 :       if (ecm_elladd(N, &gl, nbc, XT+dp*nbc, XH+rcn*nbc2, XH+rcn*nbc2+nbc2) > 1)
     933                 :            :         goto fin;
     934                 :       1610 :       rcn++;
     935                 :            :     }
     936         [ -  + ]:         35 :     if (DEBUGLEVEL >= 7) err_printf("\t(got initial helix)\n");
     937                 :            :     /* compute [210]Q etc, needed for the baby step table */
     938         [ +  - ]:         35 :     if (ellmult(N, &gl, nbc, 3, XD + (nbc<<3), X, XAUX) > 1) goto fin;
     939         [ +  - ]:         35 :     if (ellmult(N, &gl, nbc, 7, X, X, XAUX) > 1) goto fin; /* [210]Q */
     940                 :            :     /* this was the last call to ellmult() in the main loop body; may now
     941                 :            :      * overwrite XAUX and slots XD and following */
     942         [ +  - ]:         35 :     if (elldouble(N, &gl, nbc, X, XAUX) > 1) goto fin; /* [420]Q */
     943         [ +  - ]:         35 :     if (ecm_elladd(N, &gl, nbc, X, XAUX, XT) > 1) goto fin;/* [630]Q */
     944         [ +  - ]:         35 :     if (ecm_elladd(N, &gl, nbc, X, XT, XD) > 1) goto fin;  /* [840]Q */
     945         [ +  + ]:        231 :     for (i=1; i <= gse; i++)
     946         [ +  - ]:        196 :       if (elldouble(N, &gl, nbc, XT + i*nbc2, XD + i*nbc2) > 1) goto fin;
     947                 :            :     /* (the last iteration has initialized XG to [210*2^(gse+1)]Q) */
     948                 :            : 
     949         [ -  + ]:         35 :     if (DEBUGLEVEL >= 4)
     950                 :          0 :       err_printf("ECM: time = %6ld ms, entering B2 phase, p = %lu\n",
     951                 :            :                  timer_delay(&T), p);
     952                 :            : 
     953                 :            :     /* inner loop over small sets of 4 curves at a time */
     954         [ +  + ]:         91 :     for (i = nbc - 4; i >= 0; i -= 4)
     955                 :            :     {
     956         [ -  + ]:         63 :       if (DEBUGLEVEL >= 6)
     957                 :          0 :         err_printf("ECM: finishing curves %ld...%ld\n", i, i+3);
     958                 :            :       /* Copy relevant pointers from XH to Xh. Memory layout in XH:
     959                 :            :        * nbc X coordinates, nbc Y coordinates for residue class
     960                 :            :        * 1 mod 210, then the same for r.c. 11 mod 210, etc. Memory layout for
     961                 :            :        * Xh is: four X coords for 1 mod 210, four for 11 mod 210, ..., four
     962                 :            :        * for 209 mod 210, then the corresponding Y coordinates in the same
     963                 :            :        * order. This allows a giant step on Xh using just three calls to
     964                 :            :        * ecm_elladd0() each acting on 64 points in parallel */
     965         [ +  + ]:       3087 :       for (j = 48; j--; )
     966                 :            :       {
     967                 :       3024 :         k = nbc2*j + i;
     968                 :       3024 :         m = j << 2; /* X coordinates */
     969                 :       3024 :         Xh[m]   = XH[k];   Xh[m+1] = XH[k+1];
     970                 :       3024 :         Xh[m+2] = XH[k+2]; Xh[m+3] = XH[k+3];
     971                 :       3024 :         k += nbc; /* Y coordinates */
     972                 :       3024 :         Yh[m]   = XH[k];   Yh[m+1] = XH[k+1];
     973                 :       3024 :         Yh[m+2] = XH[k+2]; Yh[m+3] = XH[k+3];
     974                 :            :       }
     975                 :            :       /* Build baby step table of X coords of multiples of [210]Q.  XB[4*j]
     976                 :            :        * will point at X coords on four curves from [(j+1)*210]Q.  Until
     977                 :            :        * we're done, we need some Y coords as well, which we keep in the
     978                 :            :        * second half of the table, overwriting them at the end when gse=10.
     979                 :            :        * Multiples which we already have  (by 1,2,3,4,8,16,...,2^gse) are
     980                 :            :        * entered simply by copying the pointers, ignoring the few slots in w
     981                 :            :        * that were initially reserved for them. Here are the initial entries */
     982         [ +  + ]:        189 :       for (Xb=XB,k=2,j=i; k--; Xb=XB2,j+=nbc) /* do first X, then Y coords */
     983                 :            :       {
     984                 :        126 :         Xb[0]  = X[j];      Xb[1]  = X[j+1]; /* [210]Q */
     985                 :        126 :         Xb[2]  = X[j+2];    Xb[3]  = X[j+3];
     986                 :        126 :         Xb[4]  = XAUX[j];   Xb[5]  = XAUX[j+1]; /* [420]Q */
     987                 :        126 :         Xb[6]  = XAUX[j+2]; Xb[7]  = XAUX[j+3];
     988                 :        126 :         Xb[8]  = XT[j];     Xb[9]  = XT[j+1]; /* [630]Q */
     989                 :        126 :         Xb[10] = XT[j+2];   Xb[11] = XT[j+3];
     990                 :        126 :         Xb += 4; /* points at [420]Q */
     991                 :            :         /* ... entries at powers of 2 times 210 .... */
     992         [ +  + ]:        637 :         for (m = 2; m < (ulong)gse+k; m++) /* omit Y coords of [2^gse*210]Q */
     993                 :            :         {
     994                 :        511 :           long m2 = m*nbc2 + j;
     995                 :        511 :           Xb += (2UL<<m); /* points at [2^m*210]Q */
     996                 :        511 :           Xb[0] = XAUX[m2];   Xb[1] = XAUX[m2+1];
     997                 :        511 :           Xb[2] = XAUX[m2+2]; Xb[3] = XAUX[m2+3];
     998                 :            :         }
     999                 :            :       }
    1000         [ -  + ]:         63 :       if (DEBUGLEVEL >= 7)
    1001                 :          0 :         err_printf("\t(extracted precomputed helix / baby step entries)\n");
    1002                 :            :       /* ... glue in between, up to 16*210 ... */
    1003         [ +  - ]:         63 :       if (ecm_elladd0(N, &gl, 12, 4, /* 12 pts + (4 pts replicated thrice) */
    1004                 :            :                   XB + 12, XB2 + 12,
    1005                 :            :                   XB,      XB2,
    1006                 :            :                   XB + 16, XB2 + 16) > 1) goto fin;  /*4+{1,2,3} = {5,6,7}*/
    1007         [ +  - ]:         63 :       if (ecm_elladd0(N, &gl, 28, 4, /* 28 pts + (4 pts replicated 7fold) */
    1008                 :            :                   XB + 28, XB2 + 28,
    1009                 :            :                   XB,      XB2,
    1010                 :            :                   XB + 32, XB2 + 32) > 1) goto fin; /*8+{1,...,7} = {9,...,15}*/
    1011                 :            :       /* ... and the remainder of the lot */
    1012         [ +  + ]:        161 :       for (m = 5; m <= (ulong)gse; m++)
    1013                 :            :       { /* fill in from 2^(m-1)+1 to 2^m-1 in chunks of 64 and 60 points */
    1014                 :         98 :         ulong m2 = 2UL << m; /* will point at 2^(m-1)+1 */
    1015         [ +  + ]:        133 :         for (j = 0; (ulong)j < m2-64; j+=64) /* executed 0 times when m = 5 */
    1016                 :            :         {
    1017 [ -  + ][ +  - ]:         35 :           if (ecm_elladd0(N, &gl, 64, 4,
    1018                 :            :                       XB + m2-4, XB2 + m2-4,
    1019                 :         70 :                       XB + j,    XB2 + j,
    1020                 :         70 :                       XB + m2+j, (m<(ulong)gse? XB2+m2+j: NULL)) > 1) goto fin;
    1021                 :            :         } /* j = m2-64 here, 60 points left */
    1022 [ +  + ][ +  - ]:         98 :         if (ecm_elladd0(N, &gl, 60, 4,
    1023                 :            :                     XB + m2-4, XB2 + m2-4,
    1024                 :        196 :                     XB + j,    XB2 + j,
    1025                 :        231 :                     XB + m2+j, (m<(ulong)gse? XB2+m2+j: NULL)) > 1) goto fin;
    1026                 :            :         /* when m=gse, drop Y coords of result, and when both equal 1024,
    1027                 :            :          * overwrite Y coords of second argument with X coords of result */
    1028                 :            :       }
    1029         [ -  + ]:         63 :       if (DEBUGLEVEL >= 7) err_printf("\t(baby step table complete)\n");
    1030                 :            :       /* initialize a few other things */
    1031                 :         63 :       bstp = bstp0;
    1032                 :         63 :       p = p0; d = d0; rcn = rcn0;
    1033                 :         63 :       gl = gen_1; av1 = avma;
    1034                 :            :       /* scratchspace for prod (x_i-x_j) */
    1035                 :         63 :       avtmp = (pari_sp)new_chunk(8 * lgefint(N));
    1036                 :            :       /* The correct entry in XB to use depends on bstp and on where we are
    1037                 :            :        * on the helix. As we skip from prime to prime, bstp is incremented
    1038                 :            :        * by snextpr each time we wrap around through residue class number 0
    1039                 :            :        * (1 mod 210), but the baby step should not be taken until rcn>=rcn0,
    1040                 :            :        * i.e. until we pass again the residue class of p0.
    1041                 :            :        *
    1042                 :            :        * The correct signed multiplier is thus k = bstp - (rcn < rcn0),
    1043                 :            :        * and the offset from XB is four times (|k| - 1).  When k=0, we ignore
    1044                 :            :        * the current prime: if it had led to a factorization, this
    1045                 :            :        * would have been noted during the last giant step, or -- when we
    1046                 :            :        * first get here -- whilst initializing the helix.  When k > gss,
    1047                 :            :        * we must do a giant step and bump bstp back by -2*gss.
    1048                 :            :        *
    1049                 :            :        * The gcd of the product of X coord differences against N is taken just
    1050                 :            :        * before we do a giant step. */
    1051         [ +  + ]:     233240 :       while (p < B2)
    1052                 :            :       {/* loop over probable primes p0 < p <= nextprime(B2), inserting giant
    1053                 :            :         * steps as necessary */
    1054                 :            :         /* get next probable prime */
    1055                 :     233184 :         p = snextpr(p, &d, &rcn, &bstp, MR_Jaeschke_k2);
    1056                 :            :         /* work out the corresponding baby-step multiplier */
    1057         [ +  + ]:     233184 :         k = bstp - (rcn < rcn0 ? 1 : 0);
    1058                 :            :         /* check whether it's giant-step time */
    1059         [ +  + ]:     233184 :         if (k > gss)
    1060                 :            :         { /* take gcd */
    1061                 :         98 :           gl = gcdii(gl, N);
    1062 [ +  + ][ +  + ]:         98 :           if (!is_pm1(gl) && !equalii(gl, N)) goto fin;
    1063                 :         91 :           gl = gen_1; avma = av1;
    1064         [ +  + ]:        182 :           while (k > gss)
    1065                 :            :           { /* giant step */
    1066         [ -  + ]:         91 :             if (DEBUGLEVEL >= 7) err_printf("\t(giant step at p = %lu)\n", p);
    1067         [ +  - ]:         91 :             if (ecm_elladd0(N, &gl, 64, 4, XG + i, YG + i,
    1068                 :            :                         Xh, Yh, Xh, Yh) > 1) goto fin;
    1069         [ +  - ]:         91 :             if (ecm_elladd0(N, &gl, 64, 4, XG + i, YG + i,
    1070                 :            :                         Xh + 64, Yh + 64, Xh + 64, Yh + 64) > 1) goto fin;
    1071         [ +  - ]:         91 :             if (ecm_elladd0(N, &gl, 64, 4, XG + i, YG + i,
    1072                 :            :                         Xh + 128, Yh + 128, Xh + 128, Yh + 128) > 1) goto fin;
    1073                 :         91 :             bstp -= (gss << 1);
    1074         [ -  + ]:         91 :             k = bstp - (rcn < rcn0? 1: 0); /* recompute multiplier */
    1075                 :            :           }
    1076                 :            :         }
    1077         [ +  + ]:     233177 :         if (!k) continue; /* point of interest is already in Xh */
    1078         [ +  + ]:     229096 :         if (k < 0) k = -k;
    1079                 :     229096 :         m = ((ulong)k - 1) << 2;
    1080                 :            :         /* accumulate product of differences of X coordinates */
    1081                 :     229096 :         j = rcn<<2;
    1082                 :     229096 :         avma = avtmp; /* go to garbage zone */
    1083                 :     229096 :         gl = modii(mulii(gl, subii(XB[m],   Xh[j])), N);
    1084                 :     229096 :         gl = modii(mulii(gl, subii(XB[m+1], Xh[j+1])), N);
    1085                 :     229096 :         gl = modii(mulii(gl, subii(XB[m+2], Xh[j+2])), N);
    1086                 :     229096 :         gl = mulii(gl, subii(XB[m+3], Xh[j+3]));
    1087                 :     229096 :         avma = av1;
    1088                 :     229096 :         gl = modii(gl, N);
    1089                 :            :       } /* loop over p */
    1090                 :         56 :       avma = av1;
    1091                 :            :     } /* for i (loop over sets of 4 curves) */
    1092                 :            : 
    1093                 :            :     /* continuation part of main loop */
    1094         [ +  + ]:         28 :     if (dsn < dsnmax)
    1095                 :            :     {
    1096         [ -  + ]:         14 :       dsn += insist ? 1 : 2;
    1097         [ -  + ]:         14 :       if (dsn > dsnmax) dsn = dsnmax;
    1098                 :            :     }
    1099                 :            : 
    1100 [ +  + ][ +  - ]:         28 :     if (!insist && !--rep)
    1101                 :            :     {
    1102         [ -  + ]:         14 :       if (DEBUGLEVEL >= 4) {
    1103                 :          0 :         err_printf("ECM: time = %6ld ms,\tellfacteur giving up.\n",
    1104                 :            :                    timer_delay(&T));
    1105                 :          0 :         err_flush();
    1106                 :            :       }
    1107                 :         14 :       res = NULL; goto ret;
    1108                 :            :     }
    1109                 :         35 :   }
    1110                 :            : fin:
    1111                 :         21 :   affii(gl, res);
    1112         [ -  + ]:         21 :   if (DEBUGLEVEL >= 4) {
    1113                 :          0 :     err_printf("ECM: time = %6ld ms,\tp <= %6lu,\n\tfound factor = %Ps\n",
    1114                 :            :                timer_delay(&T), p, res);
    1115                 :          0 :     err_flush();
    1116                 :            :   }
    1117                 :            : ret:
    1118         [ -  + ]:         35 :   if (!isonstack((GEN)X)) killblock((GEN)X);
    1119                 :        574 :   avma = av; return res;
    1120                 :            : }
    1121                 :            : 
    1122                 :            : /***********************************************************************/
    1123                 :            : /**                                                                   **/
    1124                 :            : /**                FACTORIZATION (Pollard-Brent rho) --GN1998Jun18-26 **/
    1125                 :            : /**  pollardbrent() returns a nontrivial factor of n, assuming n is   **/
    1126                 :            : /**  composite and has no small prime divisor, or NULL if going on    **/
    1127                 :            : /**  would take more time than we want to spend.  Sometimes it finds  **/
    1128                 :            : /**  more than one factor, and returns a structure suitable for       **/
    1129                 :            : /**  interpretation by ifac_crack. (Cf Algo 8.5.2 in ACiCNT)          **/
    1130                 :            : /**                                                                   **/
    1131                 :            : /***********************************************************************/
    1132                 :            : #define VALUE(x) gel(x,0)
    1133                 :            : #define EXPON(x) gel(x,1)
    1134                 :            : #define CLASS(x) gel(x,2)
    1135                 :            : 
    1136                 :            : INLINE void
    1137                 :      35443 : INIT(GEN x, GEN v, GEN e, GEN c) {
    1138                 :      35443 :   VALUE(x) = v;
    1139                 :      35443 :   EXPON(x) = e;
    1140                 :      35443 :   CLASS(x) = c;
    1141                 :      35443 : }
    1142                 :            : static void
    1143                 :      31221 : ifac_delete(GEN x) { INIT(x,NULL,NULL,NULL); }
    1144                 :            : 
    1145                 :            : static void
    1146                 :          0 : rho_dbg(pari_timer *T, long c, long msg_mask)
    1147                 :            : {
    1148         [ #  # ]:          0 :   if (c & msg_mask) return;
    1149         [ #  # ]:          0 :   err_printf("Rho: time = %6ld ms,\t%3ld round%s\n",
    1150                 :            :              timer_delay(T), c, (c==1?"":"s"));
    1151                 :          0 :   err_flush();
    1152                 :            : }
    1153                 :            : 
    1154                 :            : /* Tuning parameter:  for input up to 64 bits long, we must not spend more
    1155                 :            :  * than a very short time, for fear of slowing things down on average.
    1156                 :            :  * With the current tuning formula, increase our efforts somewhat at 49 bit
    1157                 :            :  * input (an extra round for each bit at first),  and go up more and more
    1158                 :            :  * rapidly after we pass 80 bits.-- Changed this to adjust for the presence of
    1159                 :            :  * squfof, which will finish input up to 59 bits quickly. */
    1160                 :            : 
    1161                 :            : /* Return NULL when we run out of time, or a single t_INT containing a
    1162                 :            :  * nontrivial factor of n, or a vector of t_INTs, each triple of successive
    1163                 :            :  * entries containing a factor, an exponent (equal to one),  and a factor
    1164                 :            :  * class (NULL for unknown or zero for known composite),  matching the
    1165                 :            :  * internal representation used by the ifac_*() routines below.  Repeated
    1166                 :            :  * factors may arise;  the caller will sort the factors anyway. */
    1167                 :            : GEN
    1168                 :       3753 : pollardbrent(GEN n)
    1169                 :            : {
    1170                 :       3753 :   const long tune_pb_min = 14; /* even 15 seems too much. */
    1171                 :       3753 :   long tf = lgefint(n), size = 0, delta, retries = 0, msg_mask;
    1172                 :            :   long c0, c, k, k1, l;
    1173                 :       3753 :   pari_sp GGG, avP, avx, av = avma;
    1174                 :            :   GEN x, x1, y, P, g, g1, res;
    1175                 :            :   pari_timer T;
    1176                 :            : 
    1177         [ -  + ]:       3753 :   if (DEBUGLEVEL >= 4) timer_start(&T);
    1178                 :            : 
    1179         [ +  + ]:       3753 :   if (tf >= 4)
    1180                 :       1212 :     size = expi(n) + 1;
    1181         [ +  - ]:       2541 :   else if (tf == 3)                /* try to keep purify happy...  */
    1182                 :       2541 :     size = 1 + expu(uel(n,2));
    1183                 :            : 
    1184         [ -  + ]:       3753 :   if (size <= 28)
    1185                 :          0 :     c0 = 32;/* amounts very nearly to 'insist'. Now that we have squfof(), we
    1186                 :            :              * don't insist any more when input is 2^29 ... 2^32 */
    1187         [ +  + ]:       3753 :   else if (size <= 42)
    1188                 :       1093 :     c0 = tune_pb_min;
    1189         [ +  + ]:       2660 :   else if (size <= 59) /* match squfof() cutoff point */
    1190                 :       1533 :     c0 = tune_pb_min + ((size - 42)<<1);
    1191         [ +  + ]:       1127 :   else if (size <= 72)
    1192                 :        455 :     c0 = tune_pb_min + size - 24;
    1193         [ +  + ]:        672 :   else if (size <= 301)
    1194                 :            :     /* nonlinear increase in effort, kicking in around 80 bits */
    1195                 :            :     /* 301 gives 48121 + tune_pb_min */
    1196                 :        665 :     c0 = tune_pb_min + size - 60 +
    1197                 :        665 :       ((size-73)>>1)*((size-70)>>3)*((size-56)>>4);
    1198                 :            :   else
    1199                 :          7 :     c0 = 49152;        /* ECM is faster when it'd take longer */
    1200                 :            : 
    1201                 :       3753 :   c = c0 << 5; /* 2^5 iterations per round */
    1202         [ +  - ]:       7506 :   msg_mask = (size >= 448? 0x1fff:
    1203         [ +  + ]:       3753 :                            (size >= 192? (256L<<((size-128)>>6))-1: 0xff));
    1204                 :            : PB_RETRY:
    1205                 :            :  /* trick to make a 'random' choice determined by n.  Don't use x^2+0 or
    1206                 :            :   * x^2-2, ever.  Don't use x^2-3 or x^2-7 with a starting value of 2.
    1207                 :            :   * x^2+4, x^2+9 are affine conjugate to x^2+1, so don't use them either.
    1208                 :            :   *
    1209                 :            :   * (the point being that when we get called again on a composite cofactor
    1210                 :            :   * of something we've already seen, we had better avoid the same delta) */
    1211   [ +  +  +  +  :       3753 :   switch ((size + retries) & 7)
             +  +  +  + ]
    1212                 :            :   {
    1213                 :        722 :     case 0:  delta=  1; break;
    1214                 :        420 :     case 1:  delta= -1; break;
    1215                 :        651 :     case 2:  delta=  3; break;
    1216                 :        287 :     case 3:  delta=  5; break;
    1217                 :        455 :     case 4:  delta= -5; break;
    1218                 :        413 :     case 5:  delta=  7; break;
    1219                 :        364 :     case 6:  delta= 11; break;
    1220                 :            :     /* case 7: */
    1221                 :        441 :     default: delta=-11; break;
    1222                 :            :   }
    1223         [ -  + ]:       3753 :   if (DEBUGLEVEL >= 4)
    1224                 :            :   {
    1225         [ #  # ]:          0 :     if (!retries)
    1226                 :          0 :       err_printf("Rho: searching small factor of %ld-bit integer\n", size);
    1227                 :            :     else
    1228                 :          0 :       err_printf("Rho: restarting for remaining rounds...\n");
    1229                 :          0 :     err_printf("Rho: using X^2%+1ld for up to %ld rounds of 32 iterations\n",
    1230                 :            :                delta, c >> 5);
    1231                 :          0 :     err_flush();
    1232                 :            :   }
    1233                 :       3753 :   x = gen_2; P = gen_1; g1 = NULL; k = 1; l = 1;
    1234                 :       3753 :   (void)new_chunk(10 + 6 * tf); /* enough for cgetg(10) + 3 modii */
    1235                 :       3753 :   y = cgeti(tf); affsi(2, y);
    1236                 :       3753 :   x1= cgeti(tf); affsi(2, x1);
    1237                 :       3753 :   avx = avma;
    1238                 :       3753 :   avP = (pari_sp)new_chunk(2 * tf); /* enough for x = addsi(tf+1) */
    1239                 :       3753 :   GGG = (pari_sp)new_chunk(4 * tf); /* enough for P = modii(2tf+1, tf) */
    1240                 :            : 
    1241                 :            :   for (;;)                        /* terminated under the control of c */
    1242                 :            :   {
    1243                 :            :     /* use the polynomial  x^2 + delta */
    1244                 :            : #define one_iter() STMT_START {\
    1245                 :            :     avma = GGG; x = remii(sqri(x), n); /* to garbage zone */\
    1246                 :            :     avma = avx; x = addsi(delta,x);    /* erase garbage */\
    1247                 :            :     avma = GGG; P = mulii(P, subii(x1, x));\
    1248                 :            :     avma = avP; P = modii(P,n); } STMT_END
    1249                 :            : 
    1250                 :    3258048 :     one_iter();
    1251                 :            : 
    1252         [ +  + ]:    3258048 :     if ((--c & 0x1f)==0)
    1253                 :            :     { /* one round complete */
    1254         [ +  + ]:      98061 :       g = gcdii(n, P); if (!is_pm1(g)) goto fin;
    1255         [ +  + ]:      96584 :       if (c <= 0)
    1256                 :            :       {        /* getting bored */
    1257         [ -  + ]:       1576 :         if (DEBUGLEVEL >= 4)
    1258                 :            :         {
    1259                 :          0 :           err_printf("Rho: time = %6ld ms,\tPollard-Brent giving up.\n",
    1260                 :            :                      timer_delay(&T));
    1261                 :          0 :           err_flush();
    1262                 :            :         }
    1263                 :       1576 :         avma = av; return NULL;
    1264                 :            :       }
    1265                 :      95008 :       P = gen_1;                        /* not necessary, but saves 1 mulii/round */
    1266         [ -  + ]:      95008 :       if (DEBUGLEVEL >= 4) rho_dbg(&T, c0-(c>>5), msg_mask);
    1267                 :      95008 :       affii(x,y);
    1268                 :            :     }
    1269                 :            : 
    1270         [ +  + ]:    3254995 :     if (--k) continue;                /* normal end of loop body */
    1271                 :            : 
    1272         [ +  + ]:      33846 :     if (c & 0x1f) /* otherwise, we already checked */
    1273                 :            :     {
    1274         [ +  + ]:      22518 :       g = gcdii(n, P); if (!is_pm1(g)) goto fin;
    1275                 :      22497 :       P = gen_1;
    1276                 :            :     }
    1277                 :            : 
    1278                 :            :    /* Fast forward phase, doing l inner iterations without computing gcds.
    1279                 :            :     * Check first whether it would take us beyond the alloted time.
    1280                 :            :     * Fast forward rounds count only half (although they're taking
    1281                 :            :     * more like 2/3 the time of normal rounds).  This to counteract the
    1282                 :            :     * nuisance that all c0 between 4096 and 6144 would act exactly as
    1283                 :            :     * 4096;  with the halving trick only the range 4096..5120 collapses
    1284                 :            :     * (similarly for all other powers of two)
    1285                 :            :     */
    1286         [ +  + ]:      33825 :     if ((c -= (l>>1)) <= 0)
    1287                 :            :     {                                /* got bored */
    1288         [ -  + ]:        679 :       if (DEBUGLEVEL >= 4)
    1289                 :            :       {
    1290                 :          0 :         err_printf("Rho: time = %6ld ms,\tPollard-Brent giving up.\n",
    1291                 :            :                    timer_delay(&T));
    1292                 :          0 :         err_flush();
    1293                 :            :       }
    1294                 :        679 :       avma = av; return NULL;
    1295                 :            :     }
    1296                 :      33146 :     c &= ~0x1f;                        /* keep it on multiples of 32 */
    1297                 :            : 
    1298                 :            :     /* Fast forward loop */
    1299                 :      33146 :     affii(x, x1); k = l; l <<= 1;
    1300                 :            :     /* don't show this for the first several (short) fast forward phases. */
    1301 [ -  + ][ #  # ]:      33146 :     if (DEBUGLEVEL >= 4 && (l>>7) > msg_mask)
    1302                 :            :     {
    1303                 :          0 :       err_printf("Rho: fast forward phase (%ld rounds of 64)...\n", l>>7);
    1304                 :          0 :       err_flush();
    1305                 :            :     }
    1306         [ +  + ]:    4226737 :     for (k1=k; k1; k1--) one_iter();
    1307 [ -  + ][ #  # ]:      33146 :     if (DEBUGLEVEL >= 4 && (l>>7) > msg_mask)
    1308                 :            :     {
    1309                 :          0 :       err_printf("Rho: time = %6ld ms,\t%3ld rounds, back to normal mode\n",
    1310                 :          0 :                  timer_delay(&T), c0-(c>>5));
    1311                 :          0 :       err_flush();
    1312                 :            :     }
    1313                 :            : 
    1314                 :      33146 :     affii(x,y);
    1315                 :    3254295 :   } /* forever */
    1316                 :            : 
    1317                 :            : fin:
    1318                 :            :   /* An accumulated gcd was > 1 */
    1319                 :            :   /* if it isn't n, and looks prime, return it */
    1320         [ +  + ]:       1498 :   if  (!equalii(g,n))
    1321                 :            :   {
    1322         [ +  + ]:       1372 :     if (MR_Jaeschke(g,17))
    1323                 :            :     {
    1324         [ -  + ]:       1358 :       if (DEBUGLEVEL >= 4)
    1325                 :            :       {
    1326                 :          0 :         rho_dbg(&T, c0-(c>>5), 0);
    1327                 :          0 :         err_printf("\tfound factor = %Ps\n",g);
    1328                 :          0 :         err_flush();
    1329                 :            :       }
    1330                 :       1358 :       avma = av; return icopy(g);
    1331                 :            :     }
    1332                 :         14 :     avma = avx; g1 = icopy(g);  /* known composite, keep it safe */
    1333                 :         14 :     avx = avma;
    1334                 :            :   }
    1335                 :        126 :   else g1 = n;                        /* and work modulo g1 for backtracking */
    1336                 :            : 
    1337                 :            :   /* Here g1 is known composite */
    1338 [ -  + ][ #  # ]:        140 :   if (DEBUGLEVEL >= 4 && size > 192)
    1339                 :            :   {
    1340                 :          0 :     err_printf("Rho: hang on a second, we got something here...\n");
    1341                 :          0 :     err_flush();
    1342                 :            :   }
    1343                 :            :   for(;;) /* backtrack until period recovered. Must terminate */
    1344                 :            :   {
    1345                 :      10108 :     avma = GGG; y = remii(sqri(y), g1);
    1346                 :      10108 :     avma = avx; y = addsi(delta,y);
    1347         [ +  + ]:      10108 :     g = gcdii(subii(x1, y), g1); if (!is_pm1(g)) break;
    1348                 :            : 
    1349 [ -  + ][ #  # ]:       9968 :     if (DEBUGLEVEL >= 4 && (--c & 0x1f) == 0) rho_dbg(&T, c0-(c>>5), msg_mask);
    1350                 :       9968 :   }
    1351                 :            : 
    1352                 :        140 :   avma = av; /* safe */
    1353 [ +  + ][ -  + ]:        140 :   if (g1 == n || equalii(g,g1))
    1354                 :            :   {
    1355 [ +  - ][ -  + ]:        126 :     if (g1 == n && equalii(g,g1))
    1356                 :            :     { /* out of luck */
    1357         [ #  # ]:          0 :       if (DEBUGLEVEL >= 4)
    1358                 :            :       {
    1359                 :          0 :         rho_dbg(&T, c0-(c>>5), 0);
    1360                 :          0 :         err_printf("\tPollard-Brent failed.\n"); err_flush();
    1361                 :            :       }
    1362         [ #  # ]:          0 :       if (++retries >= 4) return NULL;
    1363                 :          0 :       goto PB_RETRY;
    1364                 :            :     }
    1365                 :            :     /* half lucky: we've split n, but g1 equals either g or n */
    1366         [ -  + ]:        126 :     if (DEBUGLEVEL >= 4)
    1367                 :            :     {
    1368                 :          0 :       rho_dbg(&T, c0-(c>>5), 0);
    1369         [ #  # ]:          0 :       err_printf("\tfound %sfactor = %Ps\n", (g1!=n ? "composite " : ""), g);
    1370                 :          0 :       err_flush();
    1371                 :            :     }
    1372                 :        126 :     res = cgetg(7, t_VEC);
    1373                 :            :     /* g^1: known composite when g1!=n */
    1374         [ -  + ]:        126 :     INIT(res+1, icopy(g), gen_1, (g1!=n? gen_0: NULL));
    1375                 :            :     /* cofactor^1: status unknown */
    1376                 :        126 :     INIT(res+4, diviiexact(n,g), gen_1, NULL);
    1377                 :        126 :     return res;
    1378                 :            :   }
    1379                 :            :   /* g < g1 < n : our lucky day -- we've split g1, too */
    1380                 :         14 :   res = cgetg(10, t_VEC);
    1381                 :            :   /* unknown status for all three factors */
    1382                 :         14 :   INIT(res+1, icopy(g),         gen_1, NULL);
    1383                 :         14 :   INIT(res+4, diviiexact(g1,g), gen_1, NULL);
    1384                 :         14 :   INIT(res+7, diviiexact(n,g1), gen_1, NULL);
    1385         [ -  + ]:         14 :   if (DEBUGLEVEL >= 4)
    1386                 :            :   {
    1387                 :          0 :     rho_dbg(&T, c0-(c>>5), 0);
    1388                 :          0 :     err_printf("\tfound factors = %Ps, %Ps,\n\tand %Ps\n", res[1], res[4], res[7]);
    1389                 :          0 :     err_flush();
    1390                 :            :   }
    1391                 :       3753 :   return res;
    1392                 :            : }
    1393                 :            : 
    1394                 :            : /***********************************************************************/
    1395                 :            : /**                                                                   **/
    1396                 :            : /**              FACTORIZATION (Shanks' SQUFOF) --GN2000Sep30-Oct01   **/
    1397                 :            : /**  squfof() returns a nontrivial factor of n, assuming n is odd,    **/
    1398                 :            : /**  composite, not a pure square, and has no small prime divisor,    **/
    1399                 :            : /**  or NULL if it fails to find one.  It works on two discriminants  **/
    1400                 :            : /**  simultaneously  (n and 5n for n=1(4), 3n and 4n for n=3(4)).     **/
    1401                 :            : /**  Present implementation is limited to input <2^59, and works most **/
    1402                 :            : /**  of the time in signed arithmetic on integers <2^31 in absolute   **/
    1403                 :            : /**  size. (Cf. Algo 8.7.2 in ACiCNT)                                 **/
    1404                 :            : /**                                                                   **/
    1405                 :            : /***********************************************************************/
    1406                 :            : 
    1407                 :            : /* The following is invoked to walk back along the ambiguous cycle* until we
    1408                 :            :  * hit an ambiguous form and thus the desired factor, which it returns.  If it
    1409                 :            :  * fails for any reason, it returns 0.  It doesn't interfere with timing and
    1410                 :            :  * diagnostics, which it leaves to squfof().
    1411                 :            :  *
    1412                 :            :  * Before we invoke this, we've found a form (A, B, -C) with A = a^2, where a
    1413                 :            :  * isn't blacklisted and where gcd(a, B) = 1.  According to ACiCANT, we should
    1414                 :            :  * now proceed reducing the form (a, -B, -aC), but it is easy to show that the
    1415                 :            :  * first reduction step always sends this to (-aC, B, a), and the next one,
    1416                 :            :  * with q computed as usual from B and a (occupying the c position), gives a
    1417                 :            :  * reduced form, whose third member is easiest to recover by going back to D.
    1418                 :            :  * From this point onwards, we're once again working with single-word numbers.
    1419                 :            :  * No need to track signs, just work with the abs values of the coefficients. */
    1420                 :            : static long
    1421                 :       2255 : squfof_ambig(long a, long B, long dd, GEN D)
    1422                 :            : {
    1423                 :            :   long b, c, q, qc, qcb, a0, b0, b1, c0;
    1424                 :       2255 :   long cnt = 0; /* count reduction steps on the cycle */
    1425                 :            : 
    1426                 :       2255 :   q = (dd + (B>>1)) / a;
    1427                 :       2255 :   b = ((q*a) << 1) - B;
    1428                 :            :   {
    1429                 :       2255 :     pari_sp av = avma;
    1430                 :       2255 :     c = itos(divis(shifti(subii(D, sqrs(b)), -2), a));
    1431                 :       2255 :     avma = av;
    1432                 :            :   }
    1433                 :            : #ifdef DEBUG_SQUFOF
    1434                 :            :   err_printf("SQUFOF: ambigous cycle of discriminant %Ps\n", D);
    1435                 :            :   err_printf("SQUFOF: Form on ambigous cycle (%ld, %ld, %ld)\n", a, b, c);
    1436                 :            : #endif
    1437                 :            : 
    1438                 :       2255 :   a0 = a; b0 = b1 = b;        /* end of loop detection and safeguard */
    1439                 :            : 
    1440                 :            :   for (;;) /* reduced cycles are finite */
    1441                 :            :   { /* reduction step */
    1442                 :    3858576 :     c0 = c;
    1443         [ +  + ]:    3858576 :     if (c0 > dd)
    1444                 :    1076531 :       q = 1;
    1445                 :            :     else
    1446                 :    2782045 :       q = (dd + (b>>1)) / c0;
    1447         [ +  + ]:    3858576 :     if (q == 1)
    1448                 :            :     {
    1449                 :    1604072 :       qcb = c0 - b; b = c0 + qcb; c = a - qcb;
    1450                 :            :     }
    1451                 :            :     else
    1452                 :            :     {
    1453                 :    2254504 :       qc = q*c0; qcb = qc - b; b = qc + qcb; c = a - q*qcb;
    1454                 :            :     }
    1455                 :    3858576 :     a = c0;
    1456                 :            : 
    1457         [ +  + ]:    3858576 :     cnt++; if (b == b1) break;
    1458                 :            : 
    1459                 :            :     /* safeguard against infinite loop: recognize when we've walked the entire
    1460                 :            :      * cycle in vain. (I don't think this can actually happen -- exercise.) */
    1461 [ +  + ][ -  + ]:    3856321 :     if (b == b0 && a == a0) return 0;
    1462                 :            : 
    1463                 :    3856321 :     b1 = b;
    1464                 :    3856321 :   }
    1465         [ +  + ]:       2255 :   q = a&1 ? a : a>>1;
    1466         [ -  + ]:       2255 :   if (DEBUGLEVEL >= 4)
    1467                 :            :   {
    1468         [ #  # ]:          0 :     if (q > 1)
    1469                 :          0 :       err_printf("SQUFOF: found factor %ld from ambiguous form\n"
    1470                 :            :                  "\tafter %ld steps on the ambiguous cycle\n",
    1471                 :          0 :                  q / ugcd(q,15), cnt);
    1472                 :            :     else
    1473                 :          0 :       err_printf("SQUFOF: ...found nothing on the ambiguous cycle\n"
    1474                 :            :                  "\tafter %ld steps there\n", cnt);
    1475         [ #  # ]:          0 :     if (DEBUGLEVEL >= 6) err_printf("SQUFOF: squfof_ambig returned %ld\n", q);
    1476                 :            :   }
    1477                 :       2255 :   return q;
    1478                 :            : }
    1479                 :            : 
    1480                 :            : #define SQUFOF_BLACKLIST_SZ 64
    1481                 :            : 
    1482                 :            : /* assume 2,3,5 do not divide n */
    1483                 :            : GEN
    1484                 :       2255 : squfof(GEN n)
    1485                 :            : {
    1486                 :            :   ulong d1, d2;
    1487                 :       2255 :   long tf = lgefint(n), nm4, cnt = 0;
    1488                 :            :   long a1, b1, c1, dd1, L1, a2, b2, c2, dd2, L2, a, q, c, qc, qcb;
    1489                 :            :   GEN D1, D2;
    1490                 :       2255 :   pari_sp av = avma;
    1491                 :            :   long blacklist1[SQUFOF_BLACKLIST_SZ], blacklist2[SQUFOF_BLACKLIST_SZ];
    1492                 :       2255 :   long blp1 = 0, blp2 = 0;
    1493                 :       2255 :   int act1 = 1, act2 = 1;
    1494                 :            : 
    1495                 :            : #ifdef LONG_IS_64BIT
    1496 [ +  + ][ +  - ]:       1956 :   if (tf > 3 || (tf == 3 && uel(n,2)             >= (1UL << (BITS_IN_LONG-5))))
                 [ +  + ]
    1497                 :            : #else  /* 32 bits */
    1498 [ +  + ][ +  + ]:        299 :   if (tf > 4 || (tf == 4 && (ulong)(*int_MSW(n)) >= (1UL << (BITS_IN_LONG-5))))
                 [ +  + ]
    1499                 :            : #endif
    1500                 :        546 :     return NULL; /* n too large */
    1501                 :            : 
    1502                 :            :   /* now we have 5 < n < 2^59 */
    1503                 :       1709 :   nm4 = mod4(n);
    1504         [ +  + ]:       1709 :   if (nm4 == 1)
    1505                 :            :   { /* n = 1 (mod4):  run one iteration on D1 = n, another on D2 = 5n */
    1506                 :        749 :     D1 = n;
    1507                 :        749 :     D2 = mului(5,n); d2 = itou(sqrti(D2)); dd2 = (long)((d2>>1) + (d2&1));
    1508                 :        749 :     b2 = (long)((d2-1) | 1);        /* b1, b2 will always stay odd */
    1509                 :            :   }
    1510                 :            :   else
    1511                 :            :   { /* n = 3 (mod4):  run one iteration on D1 = 3n, another on D2 = 4n */
    1512                 :        960 :     D1 = mului(3,n);
    1513                 :        960 :     D2 = shifti(n,2); dd2 = itou(sqrti(n)); d2 =  dd2 << 1;
    1514                 :        960 :     b2 = (long)(d2 & (~1UL)); /* largest even below d2, will stay even */
    1515                 :            :   }
    1516                 :       1709 :   d1 = itou(sqrti(D1));
    1517                 :       1709 :   b1 = (long)((d1-1) | 1); /* largest odd number not exceeding d1 */
    1518                 :       1709 :   c1 = itos(shifti(subii(D1, sqru((ulong)b1)), -2));
    1519         [ -  + ]:       1709 :   if (!c1) pari_err_BUG("squfof [caller of] (n or 3n is a square)");
    1520                 :       1709 :   c2 = itos(shifti(subii(D2, sqru((ulong)b2)), -2));
    1521         [ -  + ]:       1709 :   if (!c2) pari_err_BUG("squfof [caller of] (5n is a square)");
    1522                 :       1709 :   L1 = (long)usqrt(d1);
    1523                 :       1709 :   L2 = (long)usqrt(d2);
    1524                 :            :   /* dd1 used to compute floor((d1+b1)/2) as dd1+floor(b1/2), without
    1525                 :            :    * overflowing the 31bit signed integer size limit. Same for dd2. */
    1526                 :       1709 :   dd1 = (long) ((d1>>1) + (d1&1));
    1527                 :       1709 :   a1 = a2 = 1;
    1528                 :            : 
    1529                 :            :   /* The two (identity) forms (a1,b1,-c1) and (a2,b2,-c2) are now set up.
    1530                 :            :    *
    1531                 :            :    * a1 and c1 represent the absolute values of the a,c coefficients; we keep
    1532                 :            :    * track of the sign separately, via the iteration counter cnt: when cnt is
    1533                 :            :    * even, c is understood to be negative, else c is positive and a < 0.
    1534                 :            :    *
    1535                 :            :    * L1, L2 are the limits for blacklisting small leading coefficients
    1536                 :            :    * on the principal cycle, to guarantee that when we find a square form,
    1537                 :            :    * its square root will belong to an ambiguous cycle  (i.e. won't be an
    1538                 :            :    * earlier form on the principal cycle).
    1539                 :            :    *
    1540                 :            :    * When n = 3(mod 4), D2 = 12(mod 16), and b^2 is always 0 or 4 mod 16.
    1541                 :            :    * It follows that 4*a*c must be 4 or 8 mod 16, respectively, so at most
    1542                 :            :    * one of a,c can be divisible by 2 at most to the first power.  This fact
    1543                 :            :    * is used a couple of times below.
    1544                 :            :    *
    1545                 :            :    * The flags act1, act2 remain true while the respective cycle is still
    1546                 :            :    * active;  we drop them to false when we return to the identity form with-
    1547                 :            :    * out having found a square form  (or when the blacklist overflows, which
    1548                 :            :    * shouldn't happen). */
    1549         [ -  + ]:       1709 :   if (DEBUGLEVEL >= 4)
    1550                 :          0 :     err_printf("SQUFOF: entering main loop with forms\n"
    1551                 :            :                "\t(1, %ld, %ld) and (1, %ld, %ld)\n\tof discriminants\n"
    1552                 :            :                "\t%Ps and %Ps, respectively\n", b1, -c1, b2, -c2, D1, D2);
    1553                 :            : 
    1554                 :            :   /* MAIN LOOP: walk around the principal cycle looking for a square form.
    1555                 :            :    * Blacklist small leading coefficients.
    1556                 :            :    *
    1557                 :            :    * The reduction operator can be computed entirely in 32-bit arithmetic:
    1558                 :            :    * Let q = floor(floor((d1+b1)/2)/c1)  (when c1>dd1, q=1, which happens
    1559                 :            :    * often enough to special-case it).  Then the new b1 = (q*c1-b1) + q*c1,
    1560                 :            :    * which does not overflow, and the new c1 = a1 - q*(q*c1-b1), which is
    1561                 :            :    * bounded by d1 in abs size since both the old and the new a1 are positive
    1562                 :            :    * and bounded by d1. */
    1563 [ +  + ][ +  + ]:    5384794 :   while (act1 || act2)
    1564                 :            :   {
    1565         [ +  + ]:    5384780 :     if (act1)
    1566                 :            :     { /* send first form through reduction operator if active */
    1567                 :    5384696 :       c = c1;
    1568         [ +  + ]:    5384696 :       q = (c > dd1)? 1: (dd1 + (b1>>1)) / c;
    1569         [ +  + ]:    5384696 :       if (q == 1)
    1570                 :    2233714 :       { qcb = c - b1; b1 = c + qcb; c1 = a1 - qcb; }
    1571                 :            :       else
    1572                 :    3150982 :       { qc = q*c; qcb = qc - b1; b1 = qc + qcb; c1 = a1 - q*qcb; }
    1573                 :    5384696 :       a1 = c;
    1574                 :            : 
    1575         [ +  + ]:    5384696 :       if (a1 <= L1)
    1576                 :            :       { /* blacklist this */
    1577         [ -  + ]:       1561 :         if (blp1 >= SQUFOF_BLACKLIST_SZ) /* overflows: shouldn't happen */
    1578                 :          0 :           act1 = 0;                /* silently */
    1579                 :            :         else
    1580                 :            :         {
    1581         [ -  + ]:       1561 :           if (DEBUGLEVEL >= 6)
    1582                 :          0 :             err_printf("SQUFOF: blacklisting a = %ld on first cycle\n", a1);
    1583                 :       1561 :           blacklist1[blp1++] = a1;
    1584                 :            :         }
    1585                 :            :       }
    1586                 :            :     }
    1587         [ +  + ]:    5384780 :     if (act2)
    1588                 :            :     { /* send second form through reduction operator if active */
    1589                 :    5383590 :       c = c2;
    1590         [ +  + ]:    5383590 :       q = (c > dd2)? 1: (dd2 + (b2>>1)) / c;
    1591         [ +  + ]:    5383590 :       if (q == 1)
    1592                 :    2236284 :       { qcb = c - b2; b2 = c + qcb; c2 = a2 - qcb; }
    1593                 :            :       else
    1594                 :    3147306 :       { qc = q*c; qcb = qc - b2; b2 = qc + qcb; c2 = a2 - q*qcb; }
    1595                 :    5383590 :       a2 = c;
    1596                 :            : 
    1597         [ +  + ]:    5383590 :       if (a2 <= L2)
    1598                 :            :       { /* blacklist this */
    1599         [ -  + ]:       1366 :         if (blp2 >= SQUFOF_BLACKLIST_SZ) /* overflows: shouldn't happen */
    1600                 :          0 :           act2 = 0;                /* silently */
    1601                 :            :         else
    1602                 :            :         {
    1603         [ -  + ]:       1366 :           if (DEBUGLEVEL >= 6)
    1604                 :          0 :             err_printf("SQUFOF: blacklisting a = %ld on second cycle\n", a2);
    1605                 :       1366 :           blacklist2[blp2++] = a2;
    1606                 :            :         }
    1607                 :            :       }
    1608                 :            :     }
    1609                 :            : 
    1610                 :            :     /* bump counter, loop if this is an odd iteration (i.e. if the real
    1611                 :            :      * leading coefficients are negative) */
    1612         [ +  + ]:    5384780 :     if (++cnt & 1) continue;
    1613                 :            : 
    1614                 :            :     /* second half of main loop entered only when the leading coefficients
    1615                 :            :      * are positive (i.e., during even-numbered iterations) */
    1616                 :            : 
    1617                 :            :     /* examine first form if active */
    1618 [ +  + ][ +  + ]:    2692390 :     if (act1 && a1 == 1) /* back to identity */
    1619                 :            :     { /* drop this discriminant */
    1620                 :         14 :       act1 = 0;
    1621         [ -  + ]:         14 :       if (DEBUGLEVEL >= 4)
    1622                 :          0 :         err_printf("SQUFOF: first cycle exhausted after %ld iterations,\n"
    1623                 :            :                    "\tdropping it\n", cnt);
    1624                 :            :     }
    1625         [ +  + ]:    2692390 :     if (act1)
    1626                 :            :     {
    1627         [ +  + ]:    2692334 :       if (uissquareall((ulong)a1, (ulong*)&a))
    1628                 :            :       { /* square form */
    1629         [ -  + ]:       1652 :         if (DEBUGLEVEL >= 4)
    1630                 :          0 :           err_printf("SQUFOF: square form (%ld^2, %ld, %ld) on first cycle\n"
    1631                 :            :                      "\tafter %ld iterations\n", a, b1, -c1, cnt);
    1632         [ +  - ]:       1652 :         if (a <= L1)
    1633                 :            :         { /* blacklisted? */
    1634                 :            :           long j;
    1635         [ +  + ]:       3262 :           for (j = 0; j < blp1; j++)
    1636         [ +  + ]:       2303 :             if (a == blacklist1[j]) { a = 0; break; }
    1637                 :            :         }
    1638         [ +  + ]:       1652 :         if (a > 0)
    1639                 :            :         { /* not blacklisted */
    1640                 :        959 :           q = ugcd(a, b1); /* imprimitive form? */
    1641         [ -  + ]:        959 :           if (q > 1)
    1642                 :            :           { /* q^2 divides D1 hence n [ assuming n % 3 != 0 ] */
    1643                 :          0 :             avma = av;
    1644         [ #  # ]:          0 :             if (DEBUGLEVEL >= 4) err_printf("SQUFOF: found factor %ld^2\n", q);
    1645                 :          0 :             return mkvec3(utoipos(q), gen_2, NULL);/* exponent 2, unknown status */
    1646                 :            :           }
    1647                 :            :           /* chase the inverse root form back along the ambiguous cycle */
    1648                 :        959 :           q = squfof_ambig(a, b1, dd1, D1);
    1649 [ +  + ][ +  + ]:        959 :           if (nm4 == 3 && q % 3 == 0) q /= 3;
    1650         [ +  + ]:        959 :           if (q > 1) { avma = av; return utoipos(q); } /* SUCCESS! */
    1651                 :            :         }
    1652         [ -  + ]:        693 :         else if (DEBUGLEVEL >= 4) /* blacklisted */
    1653                 :          0 :           err_printf("SQUFOF: ...but the root form seems to be on the "
    1654                 :            :                      "principal cycle\n");
    1655                 :            :       }
    1656                 :            :     }
    1657                 :            : 
    1658                 :            :     /* examine second form if active */
    1659 [ +  + ][ +  + ]:    2691690 :     if (act2 && a2 == 1) /* back to identity form */
    1660                 :            :     { /* drop this discriminant */
    1661                 :         21 :       act2 = 0;
    1662         [ -  + ]:         21 :       if (DEBUGLEVEL >= 4)
    1663                 :          0 :         err_printf("SQUFOF: second cycle exhausted after %ld iterations,\n"
    1664                 :            :                    "\tdropping it\n", cnt);
    1665                 :            :     }
    1666         [ +  + ]:    2691690 :     if (act2)
    1667                 :            :     {
    1668         [ +  + ]:    2691081 :       if (uissquareall((ulong)a2, (ulong*)&a))
    1669                 :            :       { /* square form */
    1670         [ -  + ]:       1569 :         if (DEBUGLEVEL >= 4)
    1671                 :          0 :           err_printf("SQUFOF: square form (%ld^2, %ld, %ld) on second cycle\n"
    1672                 :            :                      "\tafter %ld iterations\n", a, b2, -c2, cnt);
    1673         [ +  - ]:       1569 :         if (a <= L2)
    1674                 :            :         { /* blacklisted? */
    1675                 :            :           long j;
    1676         [ +  + ]:       2830 :           for (j = 0; j < blp2; j++)
    1677         [ +  + ]:       1534 :             if (a == blacklist2[j]) { a = 0; break; }
    1678                 :            :         }
    1679         [ +  + ]:       1569 :         if (a > 0)
    1680                 :            :         { /* not blacklisted */
    1681                 :       1296 :           q = ugcd(a, b2); /* imprimitive form? */
    1682                 :            :           /* NB if b2 is even, a is odd, so the gcd is always odd */
    1683         [ -  + ]:       1296 :           if (q > 1)
    1684                 :            :           { /* q^2 divides D2 hence n [ assuming n % 5 != 0 ] */
    1685                 :          0 :             avma = av;
    1686         [ #  # ]:          0 :             if (DEBUGLEVEL >= 4) err_printf("SQUFOF: found factor %ld^2\n", q);
    1687                 :          0 :             return mkvec3(utoipos(q), gen_2, NULL);/* exponent 2, unknown status */
    1688                 :            :           }
    1689                 :            :           /* chase the inverse root form along the ambiguous cycle */
    1690                 :       1296 :           q = squfof_ambig(a, b2, dd2, D2);
    1691 [ +  + ][ +  + ]:       1296 :           if (nm4 == 1 && q % 5 == 0) q /= 5;
    1692         [ +  + ]:       1296 :           if (q > 1) { avma = av; return utoipos(q); } /* SUCCESS! */
    1693                 :            :         }
    1694         [ -  + ]:        273 :         else if (DEBUGLEVEL >= 4)        /* blacklisted */
    1695                 :          0 :           err_printf("SQUFOF: ...but the root form seems to be on the "
    1696                 :            :                      "principal cycle\n");
    1697                 :            :       }
    1698                 :            :     }
    1699                 :            :   } /* end main loop */
    1700                 :            : 
    1701                 :            :   /* both discriminants turned out to be useless. */
    1702         [ -  + ]:         14 :   if (DEBUGLEVEL>=4) err_printf("SQUFOF: giving up\n");
    1703                 :       2255 :   avma = av; return NULL;
    1704                 :            : }
    1705                 :            : 
    1706                 :            : /***********************************************************************/
    1707                 :            : /*                                                                     */
    1708                 :            : /*                    DETECTING ODD POWERS  --GN1998Jun28              */
    1709                 :            : /*   Factoring engines like MPQS which ultimately rely on computing    */
    1710                 :            : /*   gcd(N, x^2-y^2) to find a nontrivial factor of N can't split      */
    1711                 :            : /*   N = p^k for an odd prime p, since (Z/p^k)^* is then cyclic. Here  */
    1712                 :            : /*   is an analogue of Z_issquareall() for 3rd, 5th and 7th powers.    */
    1713                 :            : /*   The general case is handled by is_kth_power                       */
    1714                 :            : /*                                                                     */
    1715                 :            : /***********************************************************************/
    1716                 :            : 
    1717                 :            : /* Multistage sieve. First stages work mod 211, 209, 61, 203 in this order
    1718                 :            :  * (first reduce mod the product of these and then take the remainder apart).
    1719                 :            :  * Second stages use 117, 31, 43, 71. Moduli which are no longer interesting
    1720                 :            :  * are skipped. Everything is encoded in a table of 106 24-bit masks. We only
    1721                 :            :  * need the first half of the residues.  Three bits per modulus indicate which
    1722                 :            :  * residues are 7th (bit 2), 5th (bit 1) or 3rd (bit 0) powers; the eight
    1723                 :            :  * moduli above are assigned right-to-left. The table was generated using: */
    1724                 :            : 
    1725                 :            : #if 0
    1726                 :            : L = [71, 43, 31, [O(3^2),O(13)], [O(7),O(29)], 61, [O(11),O(19)], 211];
    1727                 :            : ispow(x, N, k)=
    1728                 :            : {
    1729                 :            :   if (type(N) == "t_INT", return (ispower(Mod(x,N), k)));
    1730                 :            :   for (i = 1, #N, if (!ispower(x + N[i], k), return (0))); 1
    1731                 :            : }
    1732                 :            : check(r) =
    1733                 :            : {
    1734                 :            :   print1("  0");
    1735                 :            :   for (i=1,#L,
    1736                 :            :     N = 0;
    1737                 :            :     if (ispow(r, L[i], 3), N += 1);
    1738                 :            :     if (ispow(r, L[i], 5), N += 2);
    1739                 :            :     if (ispow(r, L[i], 7), N += 4);
    1740                 :            :     print1(N);
    1741                 :            :   ); print("ul,  /* ", r, " */")
    1742                 :            : }
    1743                 :            : for (r = 0, 105, check(r))
    1744                 :            : #endif
    1745                 :            : static ulong powersmod[106] = {
    1746                 :            :   077777777ul,  /* 0 */
    1747                 :            :   077777777ul,  /* 1 */
    1748                 :            :   013562440ul,  /* 2 */
    1749                 :            :   012402540ul,  /* 3 */
    1750                 :            :   013562440ul,  /* 4 */
    1751                 :            :   052662441ul,  /* 5 */
    1752                 :            :   016603440ul,  /* 6 */
    1753                 :            :   016463450ul,  /* 7 */
    1754                 :            :   013573551ul,  /* 8 */
    1755                 :            :   012462540ul,  /* 9 */
    1756                 :            :   012462464ul,  /* 10 */
    1757                 :            :   013462771ul,  /* 11 */
    1758                 :            :   012406473ul,  /* 12 */
    1759                 :            :   012463641ul,  /* 13 */
    1760                 :            :   052463646ul,  /* 14 */
    1761                 :            :   012503446ul,  /* 15 */
    1762                 :            :   013562440ul,  /* 16 */
    1763                 :            :   052466440ul,  /* 17 */
    1764                 :            :   012472451ul,  /* 18 */
    1765                 :            :   012462454ul,  /* 19 */
    1766                 :            :   032463550ul,  /* 20 */
    1767                 :            :   013403664ul,  /* 21 */
    1768                 :            :   013463460ul,  /* 22 */
    1769                 :            :   032562565ul,  /* 23 */
    1770                 :            :   012402540ul,  /* 24 */
    1771                 :            :   052662441ul,  /* 25 */
    1772                 :            :   032672452ul,  /* 26 */
    1773                 :            :   013573551ul,  /* 27 */
    1774                 :            :   012467541ul,  /* 28 */
    1775                 :            :   012567640ul,  /* 29 */
    1776                 :            :   032706450ul,  /* 30 */
    1777                 :            :   012762452ul,  /* 31 */
    1778                 :            :   033762662ul,  /* 32 */
    1779                 :            :   012502562ul,  /* 33 */
    1780                 :            :   032463562ul,  /* 34 */
    1781                 :            :   013563440ul,  /* 35 */
    1782                 :            :   016663440ul,  /* 36 */
    1783                 :            :   036662550ul,  /* 37 */
    1784                 :            :   012462552ul,  /* 38 */
    1785                 :            :   033502450ul,  /* 39 */
    1786                 :            :   012462643ul,  /* 40 */
    1787                 :            :   033467540ul,  /* 41 */
    1788                 :            :   017403441ul,  /* 42 */
    1789                 :            :   017463462ul,  /* 43 */
    1790                 :            :   017472460ul,  /* 44 */
    1791                 :            :   033462470ul,  /* 45 */
    1792                 :            :   052566450ul,  /* 46 */
    1793                 :            :   013562640ul,  /* 47 */
    1794                 :            :   032403640ul,  /* 48 */
    1795                 :            :   016463450ul,  /* 49 */
    1796                 :            :   016463752ul,  /* 50 */
    1797                 :            :   033402440ul,  /* 51 */
    1798                 :            :   012462540ul,  /* 52 */
    1799                 :            :   012472540ul,  /* 53 */
    1800                 :            :   053562462ul,  /* 54 */
    1801                 :            :   012463465ul,  /* 55 */
    1802                 :            :   012663470ul,  /* 56 */
    1803                 :            :   052607450ul,  /* 57 */
    1804                 :            :   012566553ul,  /* 58 */
    1805                 :            :   013466440ul,  /* 59 */
    1806                 :            :   012502741ul,  /* 60 */
    1807                 :            :   012762744ul,  /* 61 */
    1808                 :            :   012763740ul,  /* 62 */
    1809                 :            :   012763443ul,  /* 63 */
    1810                 :            :   013573551ul,  /* 64 */
    1811                 :            :   013462471ul,  /* 65 */
    1812                 :            :   052502460ul,  /* 66 */
    1813                 :            :   012662463ul,  /* 67 */
    1814                 :            :   012662451ul,  /* 68 */
    1815                 :            :   012403550ul,  /* 69 */
    1816                 :            :   073567540ul,  /* 70 */
    1817                 :            :   072463445ul,  /* 71 */
    1818                 :            :   072462740ul,  /* 72 */
    1819                 :            :   012472442ul,  /* 73 */
    1820                 :            :   012462644ul,  /* 74 */
    1821                 :            :   013406650ul,  /* 75 */
    1822                 :            :   052463471ul,  /* 76 */
    1823                 :            :   012563474ul,  /* 77 */
    1824                 :            :   013503460ul,  /* 78 */
    1825                 :            :   016462441ul,  /* 79 */
    1826                 :            :   016462440ul,  /* 80 */
    1827                 :            :   012462540ul,  /* 81 */
    1828                 :            :   013462641ul,  /* 82 */
    1829                 :            :   012463454ul,  /* 83 */
    1830                 :            :   013403550ul,  /* 84 */
    1831                 :            :   057563540ul,  /* 85 */
    1832                 :            :   017466441ul,  /* 86 */
    1833                 :            :   017606471ul,  /* 87 */
    1834                 :            :   053666573ul,  /* 88 */
    1835                 :            :   012562561ul,  /* 89 */
    1836                 :            :   013473641ul,  /* 90 */
    1837                 :            :   032573440ul,  /* 91 */
    1838                 :            :   016763440ul,  /* 92 */
    1839                 :            :   016702640ul,  /* 93 */
    1840                 :            :   033762552ul,  /* 94 */
    1841                 :            :   012562550ul,  /* 95 */
    1842                 :            :   052402451ul,  /* 96 */
    1843                 :            :   033563441ul,  /* 97 */
    1844                 :            :   012663561ul,  /* 98 */
    1845                 :            :   012677560ul,  /* 99 */
    1846                 :            :   012462464ul,  /* 100 */
    1847                 :            :   032562642ul,  /* 101 */
    1848                 :            :   013402551ul,  /* 102 */
    1849                 :            :   032462450ul,  /* 103 */
    1850                 :            :   012467445ul,  /* 104 */
    1851                 :            :   032403440ul,  /* 105 */
    1852                 :            : };
    1853                 :            : 
    1854                 :            : static int
    1855                 :    1837404 : check_res(ulong x, ulong N, int shift, ulong *mask)
    1856                 :            : {
    1857         [ +  + ]:    1837404 :   long r = x%N; if ((ulong)r> (N>>1)) r = N - r;
    1858                 :    1837404 :   *mask &= (powersmod[r] >> shift);
    1859                 :    1837404 :   return *mask;
    1860                 :            : }
    1861                 :            : 
    1862                 :            : /* is x mod 211*209*61*203*117*31*43*71 a 3rd, 5th or 7th power ? */
    1863                 :            : int
    1864                 :    1016235 : uis_357_powermod(ulong x, ulong *mask)
    1865                 :            : {
    1866         [ +  + ]:    1016235 :   if (             !check_res(x, 211UL, 0, mask)) return 0;
    1867 [ +  + ][ +  + ]:     521013 :   if (*mask & 3 && !check_res(x, 209UL, 3, mask)) return 0;
    1868 [ +  + ][ +  + ]:     252139 :   if (*mask & 3 && !check_res(x,  61UL, 6, mask)) return 0;
    1869 [ +  + ][ +  + ]:     166146 :   if (*mask & 5 && !check_res(x, 203UL, 9, mask)) return 0;
    1870 [ +  + ][ +  + ]:      37622 :   if (*mask & 1 && !check_res(x, 117UL,12, mask)) return 0;
    1871 [ +  + ][ +  + ]:      27045 :   if (*mask & 3 && !check_res(x,  31UL,15, mask)) return 0;
    1872 [ +  + ][ +  + ]:      20621 :   if (*mask & 5 && !check_res(x,  43UL,18, mask)) return 0;
    1873 [ +  + ][ +  + ]:       4486 :   if (*mask & 6 && !check_res(x,  71UL,21, mask)) return 0;
    1874                 :    1016235 :   return 1;
    1875                 :            : }
    1876                 :            : /* asume x > 0 and pt != NULL */
    1877                 :            : int
    1878                 :     982608 : uis_357_power(ulong x, ulong *pt, ulong *mask)
    1879                 :            : {
    1880                 :            :   double logx;
    1881         [ +  + ]:     982608 :   if (!odd(x))
    1882                 :            :   {
    1883                 :        210 :     long v = vals(x);
    1884         [ +  - ]:        210 :     if (v % 7) *mask &= ~4;
    1885         [ +  + ]:        210 :     if (v % 5) *mask &= ~2;
    1886         [ +  + ]:        210 :     if (v % 3) *mask &= ~1;
    1887         [ +  + ]:        210 :     if (!*mask) return 0;
    1888                 :            :   }
    1889         [ +  + ]:     982503 :   if (!uis_357_powermod(x, mask)) return 0;
    1890                 :       1037 :   logx = log((double)x);
    1891         [ +  + ]:       1821 :   while (*mask)
    1892                 :            :   {
    1893                 :            :     long e, b;
    1894                 :            :     ulong y, ye;
    1895         [ +  + ]:       1037 :     if (*mask & 1)      { b = 1; e = 3; }
    1896         [ +  + ]:        644 :     else if (*mask & 2) { b = 2; e = 5; }
    1897                 :        343 :     else                { b = 4; e = 7; }
    1898                 :       1037 :     y = (ulong)(exp(logx / e) + 0.5);
    1899                 :       1037 :     ye = upowuu(y,e);
    1900         [ +  + ]:       1037 :     if (ye == x) { *pt = y; return e; }
    1901                 :            : #ifdef LONG_IS_64BIT
    1902         [ +  + ]:        672 :     if (ye > x) y--; else y++;
    1903                 :        672 :     ye = upowuu(y,e);
    1904         [ -  + ]:        672 :     if (ye == x) { *pt = y; return e; }
    1905                 :            : #endif
    1906                 :        784 :     *mask &= ~b; /* turn the bit off */
    1907                 :            :   }
    1908                 :     982608 :   return 0;
    1909                 :            : }
    1910                 :            : 
    1911                 :            : #ifndef LONG_IS_64BIT
    1912                 :            : /* as above, split in two functions */
    1913                 :            : /* is x mod 211*209*61*203 a 3rd, 5th or 7th power ? */
    1914                 :            : static int
    1915                 :       8805 : uis_357_powermod_32bit_1(ulong x, ulong *mask)
    1916                 :            : {
    1917         [ +  + ]:       8805 :   if (             !check_res(x, 211UL, 0, mask)) return 0;
    1918 [ +  + ][ +  + ]:       4845 :   if (*mask & 3 && !check_res(x, 209UL, 3, mask)) return 0;
    1919 [ +  + ][ +  + ]:       2463 :   if (*mask & 3 && !check_res(x,  61UL, 6, mask)) return 0;
    1920 [ +  + ][ +  + ]:       1699 :   if (*mask & 5 && !check_res(x, 203UL, 9, mask)) return 0;
    1921                 :       8805 :   return 1;
    1922                 :            : }
    1923                 :            : /* is x mod 117*31*43*71 a 3rd, 5th or 7th power ? */
    1924                 :            : static int
    1925                 :        425 : uis_357_powermod_32bit_2(ulong x, ulong *mask)
    1926                 :            : {
    1927 [ +  + ][ +  + ]:        425 :   if (*mask & 1 && !check_res(x, 117UL,12, mask)) return 0;
    1928 [ +  + ][ +  + ]:        327 :   if (*mask & 3 && !check_res(x,  31UL,15, mask)) return 0;
    1929 [ +  + ][ +  + ]:        245 :   if (*mask & 5 && !check_res(x,  43UL,18, mask)) return 0;
    1930 [ +  + ][ +  + ]:         83 :   if (*mask & 6 && !check_res(x,  71UL,21, mask)) return 0;
    1931                 :        425 :   return 1;
    1932                 :            : }
    1933                 :            : #endif
    1934                 :            : 
    1935                 :            : /* Returns 3, 5, or 7 if x is a cube (but not a 5th or 7th power),  a 5th
    1936                 :            :  * power (but not a 7th),  or a 7th power, and in this case creates the
    1937                 :            :  * base on the stack and assigns its address to *pt.  Otherwise returns 0.
    1938                 :            :  * x must be of type t_INT and positive;  this is not checked.  The *mask
    1939                 :            :  * argument tells us which things to check -- bit 0: 3rd, bit 1: 5th,
    1940                 :            :  * bit 2: 7th pwr;  set a bit to have the corresponding power examined --
    1941                 :            :  * and is updated appropriately for a possible follow-up call */
    1942                 :            : int
    1943                 :    1382352 : is_357_power(GEN x, GEN *pt, ulong *mask)
    1944                 :            : {
    1945                 :    1382352 :   long lx = lgefint(x);
    1946                 :            :   ulong r;
    1947                 :            :   pari_sp av;
    1948                 :            :   GEN y;
    1949                 :            : 
    1950         [ +  + ]:    1382352 :   if (!*mask) return 0; /* useful when running in a loop */
    1951         [ -  + ]:    1012437 :   if (DEBUGLEVEL>4) err_printf("OddPwrs: examining %ld-bit integer\n", expi(x));
    1952         [ +  + ]:    1012437 :   if (lgefint(x) == 3) {
    1953                 :            :     ulong t;
    1954                 :     969900 :     long e = uis_357_power(x[2], &t, mask);
    1955         [ +  + ]:     969900 :     if (e)
    1956                 :            :     {
    1957         [ +  + ]:        253 :       if (pt) *pt = utoi(t);
    1958                 :        253 :       return e;
    1959                 :            :     }
    1960                 :     969900 :     return 0;
    1961                 :            :   }
    1962                 :            : #ifdef LONG_IS_64BIT
    1963         [ -  + ]:      33732 :   r = (lx == 3)? uel(x,2): umodiu(x, 6046846918939827UL);
    1964         [ +  + ]:      33732 :   if (!uis_357_powermod(r, mask)) return 0;
    1965                 :            : #else
    1966         [ -  + ]:       8805 :   r = (lx == 3)? uel(x,2): umodiu(x, 211*209*61*203);
    1967         [ +  + ]:       8805 :   if (!uis_357_powermod_32bit_1(r, mask)) return 0;
    1968         [ -  + ]:        425 :   r = (lx == 3)? uel(x,2): umodiu(x, 117*31*43*71);
    1969         [ +  + ]:        425 :   if (!uis_357_powermod_32bit_2(r, mask)) return 0;
    1970                 :            : #endif
    1971                 :        302 :   av = avma;
    1972         [ +  + ]:        360 :   while (*mask)
    1973                 :            :   {
    1974                 :            :     long e, b;
    1975                 :            :     /* priority to higher powers: if we have a 21st, it is easier to rediscover
    1976                 :            :      * that its 7th root is a cube than that its cube root is a 7th power */
    1977         [ +  + ]:        302 :          if (*mask & 4) { b = 4; e = 7; }
    1978         [ +  + ]:        244 :     else if (*mask & 2) { b = 2; e = 5; }
    1979                 :        118 :     else                { b = 1; e = 3; }
    1980                 :        302 :     y = mpround( sqrtnr(itor(x, nbits2prec(64 + bit_accuracy(lx) / e)), e) );
    1981         [ +  + ]:        302 :     if (equalii(powiu(y,e), x))
    1982                 :            :     {
    1983         [ -  + ]:        244 :       if (!pt) { avma = av; return e; }
    1984                 :        244 :       avma = (pari_sp)y; *pt = gerepileuptoint(av, y);
    1985                 :        244 :       return e;
    1986                 :            :     }
    1987         [ -  + ]:         58 :     if (DEBUGLEVEL>4)
    1988                 :          0 :       err_printf("\tBut it nevertheless wasn't a %ld%s power.\n", e,eng_ord(e));
    1989                 :         58 :     *mask &= ~b; /* turn the bit off */
    1990                 :         58 :     avma = av;
    1991                 :            :   }
    1992                 :    1382352 :   return 0;
    1993                 :            : }
    1994                 :            : 
    1995                 :            : /* Is x a n-th power ?
    1996                 :            :  * if d = NULL, n not necessarily prime, otherwise, n prime and d the
    1997                 :            :  * corresponding diffptr to go on looping over primes.
    1998                 :            :  * If pt != NULL, it receives the n-th root */
    1999                 :            : ulong
    2000                 :      30359 : is_kth_power(GEN x, ulong n, GEN *pt)
    2001                 :            : {
    2002                 :            :   forprime_t T;
    2003                 :            :   long j;
    2004                 :            :   ulong q, residue;
    2005                 :            :   GEN y;
    2006                 :      30359 :   pari_sp av = avma;
    2007                 :            : 
    2008         [ +  + ]:      30359 :   (void)u_forprime_arith_init(&T, odd(n)? 2*n+1: n+1, ULONG_MAX, 1,n);
    2009                 :            :   /* we'll start at q, smallest prime >= n */
    2010                 :            : 
    2011                 :            :   /* Modular checks, use small primes q congruent 1 mod n */
    2012                 :            :   /* A non n-th power nevertheless passes the test with proba n^(-#checks),
    2013                 :            :    * We'd like this < 1e-6 but let j = floor(log(1e-6) / log(n)) which
    2014                 :            :    * ensures much less. */
    2015         [ +  + ]:      30359 :   if (n < 16)
    2016                 :       3892 :     j = 5;
    2017         [ +  + ]:      26467 :   else if (n < 32)
    2018                 :       4928 :     j = 4;
    2019         [ +  + ]:      21539 :   else if (n < 101)
    2020                 :       3570 :     j = 3;
    2021         [ +  + ]:      17969 :   else if (n < 1001)
    2022                 :       1666 :     j = 2;
    2023         [ +  - ]:      16303 :   else if (n < 17886697) /* smallest such that smallest suitable q is > 2^32 */
    2024                 :      16303 :     j = 1;
    2025                 :            :   else
    2026                 :          0 :     j = 0;
    2027         [ +  + ]:      31423 :   for (; j > 0; j--)
    2028                 :            :   {
    2029         [ -  + ]:      31332 :     if (!(q = u_forprime_next(&T))) break;
    2030                 :            :     /* q a prime = 1 mod n */
    2031                 :      31332 :     residue = umodiu(x, q);
    2032         [ +  + ]:      31332 :     if (residue == 0)
    2033                 :            :     {
    2034         [ +  - ]:         35 :       if (Z_lval(x,q) % n) { avma = av; return 0; }
    2035                 :          0 :       continue;
    2036                 :            :     }
    2037                 :            :     /* n-th power mod q ? */
    2038         [ +  + ]:      31297 :     if (Fl_powu(residue, (q-1)/n, q) != 1) { avma = av; return 0; }
    2039                 :            :   }
    2040                 :         91 :   avma = av;
    2041                 :            : 
    2042         [ -  + ]:         91 :   if (DEBUGLEVEL>4) err_printf("\nOddPwrs: [%lu] passed modular checks\n",n);
    2043                 :            :   /* go to the horse's mouth... */
    2044                 :         91 :   y = roundr( sqrtnr(itor(x, nbits2prec((expi(x)+16*n)/n)), n) );
    2045         [ -  + ]:         91 :   if (!equalii(powiu(y, n), x)) {
    2046         [ #  # ]:          0 :     if (DEBUGLEVEL>4) err_printf("\tBut it wasn't a pure power.\n");
    2047                 :          0 :     avma = av; return 0;
    2048                 :            :   }
    2049         [ -  + ]:         91 :   if (!pt) avma = av; else { avma = (pari_sp)y; *pt = gerepileuptoint(av, y); }
    2050                 :      30359 :   return 1;
    2051                 :            : }
    2052                 :            : 
    2053                 :            : /* is x a p^i-th power, p >= 11 prime ? Similar to is_357_power(), but instead
    2054                 :            :  * of the mask, we keep the current test exponent around. Cut off when
    2055                 :            :  * log_2 x^(1/k) < cutoffbits since we would have found it by trial division.
    2056                 :            :  * Everything needed here (primitive roots etc.) is computed from scratch on
    2057                 :            :  * the fly; compared to the size of numbers under consideration, these
    2058                 :            :  * word-sized computations take negligible time.
    2059                 :            :  * Any cutoffbits > 0 is safe, but direct root extraction attempts are faster
    2060                 :            :  * when trial division has been used to discover very small bases. We become
    2061                 :            :  * competitive at cutoffbits ~ 10 */
    2062                 :            : int
    2063                 :      28452 : is_pth_power(GEN x, GEN *pt, forprime_t *T, ulong cutoffbits)
    2064                 :            : {
    2065                 :      28452 :   long cnt=0, size = expi(x) /* not +1 */;
    2066                 :            :   ulong p;
    2067                 :      28452 :   pari_sp av = avma;
    2068 [ +  + ][ +  + ]:      58594 :   while ((p = u_forprime_next(T)) && size/p >= cutoffbits) {
    2069                 :      30163 :     long v = 1;
    2070 [ -  + ][ #  # ]:      30163 :     if (DEBUGLEVEL>5 && cnt++==2000)
    2071                 :          0 :       { cnt=0; err_printf("%lu%% ", 100*p*cutoffbits/size); }
    2072         [ +  + ]:      30198 :     while (is_kth_power(x, p, pt)) {
    2073                 :         35 :       v *= p; x = *pt;
    2074                 :         35 :       size = expi(x);
    2075                 :            :     }
    2076         [ +  + ]:      30163 :     if (v > 1)
    2077                 :            :     {
    2078         [ -  + ]:         21 :       if (DEBUGLEVEL>5) err_printf("\nOddPwrs: is a %ld power\n",v);
    2079                 :         21 :       return v;
    2080                 :            :     }
    2081                 :            :   }
    2082         [ -  + ]:      28431 :   if (DEBUGLEVEL>5) err_printf("\nOddPwrs: not a power\n",p);
    2083                 :      28452 :   avma = av; return 0; /* give up */
    2084                 :            : }
    2085                 :            : 
    2086                 :            : /***********************************************************************/
    2087                 :            : /**                                                                   **/
    2088                 :            : /**                FACTORIZATION  (master iteration)                  **/
    2089                 :            : /**      Driver for the various methods of finding large factors      **/
    2090                 :            : /**      (after trial division has cast out the very small ones).     **/
    2091                 :            : /**                        GN1998Jun24--30                            **/
    2092                 :            : /**                                                                   **/
    2093                 :            : /***********************************************************************/
    2094                 :            : 
    2095                 :            : /* Direct use:
    2096                 :            :  *  ifac_start_hint(n,moebius,hint) registers with the iterative factorizer
    2097                 :            :  *  - an integer n (without prime factors  < tridiv_bound(n))
    2098                 :            :  *  - registers whether or not we should terminate early if we find a square
    2099                 :            :  *    factor,
    2100                 :            :  *  - a hint about which method(s) to use.
    2101                 :            :  *  This must always be called first. If input is not composite, oo loop.
    2102                 :            :  *  The routine decomposes n nontrivially into a product of two factors except
    2103                 :            :  *  in squarefreeness ('Moebius') mode.
    2104                 :            :  *
    2105                 :            :  *  ifac_start(n,moebius) same using default hint.
    2106                 :            :  *
    2107                 :            :  *  ifac_primary_factor()  returns a prime divisor (not necessarily the
    2108                 :            :  *    smallest) and the corresponding exponent.
    2109                 :            :  *
    2110                 :            :  * Encapsulated user interface: Many arithmetic functions have a 'contributor'
    2111                 :            :  * ifac_xxx, to be called on any large composite cofactor left over after trial
    2112                 :            :  * division by small primes: xxx is one of moebius, issquarefree, totient, etc.
    2113                 :            :  *
    2114                 :            :  * We never test whether the input number is prime or composite, since
    2115                 :            :  * presumably it will have come out of the small factors finder stage
    2116                 :            :  * (which doesn't really exist yet but which will test the left-over
    2117                 :            :  * cofactor for primality once it does). */
    2118                 :            : 
    2119                 :            : /* The data structure in which we preserve whatever we know about our number N
    2120                 :            :  * is kept on the PARI stack, and updated as needed.
    2121                 :            :  * This makes the machinery re-entrant, and avoids memory leaks when a lengthy
    2122                 :            :  * factorization is interrupted. We try to keep the whole affair connected,
    2123                 :            :  * and the parent object is always older than its children.  This may in
    2124                 :            :  * rare cases lead to some extra copying around, and knowing what is garbage
    2125                 :            :  * at any given time is not trivial. See below for examples how to do it right.
    2126                 :            :  * (Connectedness is destroyed if callers of ifac_main() create stuff on the
    2127                 :            :  * stack in between calls. This is harmless as long as ifac_realloc() is used
    2128                 :            :  * to re-create a connected object at the head of the stack just before
    2129                 :            :  * collecting garbage.)
    2130                 :            :  * A t_INT may well have > 10^6 distinct prime factors larger than 2^16. Since
    2131                 :            :  * we need not find factors in order of increasing size, we must be prepared to
    2132                 :            :  * drag a very large amount of data around.  We start with a small structure
    2133                 :            :  * and extend it when necessary. */
    2134                 :            : 
    2135                 :            : /* The idea of the algorithm is:
    2136                 :            :  * Let N0 be whatever is currently left of N after dividing off all the
    2137                 :            :  * prime powers we have already returned to the caller.  Then we maintain
    2138                 :            :  * N0 as a product
    2139                 :            :  * (1) N0 = \prod_i P_i^{e_i} * \prod_j Q_j^{f_j} * \prod_k C_k^{g_k}
    2140                 :            :  * where the P_i and Q_j are distinct primes, each C_k is known composite,
    2141                 :            :  * none of the P_i divides any C_k, and we also know the total ordering
    2142                 :            :  * of all the P_i, Q_j and C_k; in particular, we will never try to divide
    2143                 :            :  * a C_k by a larger Q_j.  Some of the C_k may have common factors.
    2144                 :            :  *
    2145                 :            :  * Caveat implementor:  Taking gcds among C_k's is very likely to cost at
    2146                 :            :  * least as much time as dividing off any primes as we find them, and book-
    2147                 :            :  * keeping would be tough (since D=gcd(C_1,C_2) can still have common factors
    2148                 :            :  * with both C_1/D and C_2/D, and so on...).
    2149                 :            :  *
    2150                 :            :  * At startup, we just initialize the structure to
    2151                 :            :  * (2) N = C_1^1   (composite).
    2152                 :            :  *
    2153                 :            :  * Whenever ifac_primary_factor() or one of the arithmetic user interface
    2154                 :            :  * routines needs a primary factor, and the smallest thing in our list is P_1,
    2155                 :            :  * we return that and its exponent, and remove it from our list. (When nothing
    2156                 :            :  * is left, we return a sentinel value -- gen_1.  And in Moebius mode, when we
    2157                 :            :  * see something with exponent > 1, whether prime or composite, we return gen_0
    2158                 :            :  * or 0, depending on the function). In all other cases, ifac_main() iterates
    2159                 :            :  * the following steps until we have a P_1 in the smallest position.
    2160                 :            :  *
    2161                 :            :  * When the smallest item is C_1, as it is initially:
    2162                 :            :  * (3.1) Crack C_1 into a nontrivial product  U_1 * U_2  by whatever method
    2163                 :            :  * comes to mind for this size. (U for 'unknown'.)  Cracking will detect
    2164                 :            :  * perfect powers, so we may instead see a power of some U_1 here, or even
    2165                 :            :  * something of the form U_1^k*U_2^k; of course the exponent already attached
    2166                 :            :  * to C_1 is taken into account in the following.
    2167                 :            :  * (3.2) If we have U_1*U_2, sort the two factors (distinct: squares are caught
    2168                 :            :  * in stage 3.1). N.B. U_1 and U_2 are smaller than anything else in our list.
    2169                 :            :  * (3.3) Check U_1 and U_2 for primality, and flag them accordingly.
    2170                 :            :  * (3.4) Iterate.
    2171                 :            :  *
    2172                 :            :  * When the smallest item is Q_1:
    2173                 :            :  * This is the unpleasant case.  We go through the entire list and try to
    2174                 :            :  * divide Q_1 off each of the current C_k's, which usually fails, but may
    2175                 :            :  * succeed several times. When a division was successful, the corresponding
    2176                 :            :  * C_k is removed from our list, and the cofactor becomes a U_l for the moment
    2177                 :            :  * unless it is 1 (which happens when C_k was a power of Q_1).  When we're
    2178                 :            :  * through we upgrade Q_1 to P_1 status, then do a primality check on each U_l
    2179                 :            :  * and sort it back into the list either as a Q_j or as a C_k.  If during the
    2180                 :            :  * insertion sort we discover that some U_l equals some P_i or Q_j or C_k we
    2181                 :            :  * already have, we just add U_l's exponent to that of its twin. (The sorting
    2182                 :            :  * therefore happens before the primality test). Since this may produce one or
    2183                 :            :  * more elements smaller than the P_1 we just confirmed, we may have to repeat
    2184                 :            :  * the iteration.
    2185                 :            :  * A trick avoids some Q_1 instances: just after the sweep classifying
    2186                 :            :  * all current unknowns as either composites or primes, we do another downward
    2187                 :            :  * sweep beginning with the largest current factor and stopping just above the
    2188                 :            :  * largest current composite.  Every Q_j we pass is turned into a P_i.
    2189                 :            :  * (Different primes are automatically coprime among each other, and primes do
    2190                 :            :  * not divide smaller composites.)
    2191                 :            :  * NB: We have no use for comparing the square of a prime to N0.  Normally
    2192                 :            :  * we will get called after casting out only the smallest primes, and
    2193                 :            :  * since we cannot guarantee that we see the large prime factors in as-
    2194                 :            :  * cending order, we cannot stop when we find one larger than sqrt(N0). */
    2195                 :            : 
    2196                 :            : /* Data structure: We keep everything in a single t_VEC of t_INTs.  The
    2197                 :            :  * first 2 components are read-only:
    2198                 :            :  * 1) the first records whether we're doing full (NULL) or Moebius (gen_1)
    2199                 :            :  * factorization; in the latter case subroutines return a sentinel value as
    2200                 :            :  * soon as they spot an exponent > 1.
    2201                 :            :  * 2) the second records the hint from factorint()'s optional flag, for use by
    2202                 :            :  * ifac_crack().
    2203                 :            :  *
    2204                 :            :  * The remaining components (initially 15) are used in groups of three:
    2205                 :            :  * [ factor (t_INT), exponent (t_INT), factor class ], where factor class is
    2206                 :            :  *  NULL : unknown
    2207                 :            :  *  gen_0: known composite C_k
    2208                 :            :  *  gen_1: known prime Q_j awaiting trial division
    2209                 :            :  *  gen_2: finished prime P_i.
    2210                 :            :  * When during the division stage we re-sort a C_k-turned-U_l to a lower
    2211                 :            :  * position, we rotate any intervening material upward towards its old
    2212                 :            :  * slot.  When a C_k was divided down to 1, its slot is left empty at
    2213                 :            :  * first; similarly when the re-sorting detects a repeated factor.
    2214                 :            :  * After the sorting phase, we de-fragment the list and squeeze all the
    2215                 :            :  * occupied slots together to the high end, so that ifac_crack() has room
    2216                 :            :  * for new factors.  When this doesn't suffice, we abandon the current vector
    2217                 :            :  * and allocate a somewhat larger one, defragmenting again while copying.
    2218                 :            :  *
    2219                 :            :  * For internal use: note that all exponents will fit into C longs, given
    2220                 :            :  * PARI's lgefint field size.  When we work with them, we sometimes read
    2221                 :            :  * out the GEN pointer, and sometimes do an itos, whatever is more con-
    2222                 :            :  * venient for the task at hand. */
    2223                 :            : 
    2224                 :            : /*** Overview ***/
    2225                 :            : 
    2226                 :            : /* The '*where' argument in the following points into *partial at the first of
    2227                 :            :  * the three fields of the first occupied slot.  It's there because the caller
    2228                 :            :  * would already know where 'here' is, so we don't want to search for it again.
    2229                 :            :  * We do not preserve this from one user-interface call to the next. */
    2230                 :            : 
    2231                 :            : /* In the most common cases, control flows from the user interface to
    2232                 :            :  * ifac_main() and then to a succession of ifac_crack()s and ifac_divide()s,
    2233                 :            :  * with (typically) none of the latter finding anything. */
    2234                 :            : 
    2235                 :            : static long ifac_insert_multiplet(GEN *, GEN *, GEN, long);
    2236                 :            : 
    2237                 :            : #define LAST(x) x+lg(x)-3
    2238                 :            : #define FIRST(x) x+3
    2239                 :            : 
    2240                 :            : #define MOEBIUS(x) gel(x,1)
    2241                 :            : #define HINT(x) gel(x,2)
    2242                 :            : 
    2243                 :            : /* y <- x */
    2244                 :            : INLINE void
    2245                 :          0 : SHALLOWCOPY(GEN x, GEN y) {
    2246                 :          0 :   VALUE(y) = VALUE(x);
    2247                 :          0 :   EXPON(y) = EXPON(x);
    2248                 :          0 :   CLASS(y) = CLASS(x);
    2249                 :          0 : }
    2250                 :            : /* y <- x */
    2251                 :            : INLINE void
    2252                 :          0 : COPY(GEN x, GEN y) {
    2253 [ #  # ][ #  # ]:          0 :   icopyifstack(VALUE(x), VALUE(y));
    2254 [ #  # ][ #  # ]:          0 :   icopyifstack(EXPON(x), EXPON(y));
    2255                 :          0 :   CLASS(y) = CLASS(x);
    2256                 :          0 : }
    2257                 :            : 
    2258                 :            : /* Diagnostics */
    2259                 :            : static void
    2260                 :          0 : ifac_factor_dbg(GEN x)
    2261                 :            : {
    2262                 :          0 :   GEN c = CLASS(x), v = VALUE(x);
    2263         [ #  # ]:          0 :   if (c == gen_2) err_printf("IFAC: factor %Ps\n\tis prime (finished)\n", v);
    2264         [ #  # ]:          0 :   else if (c == gen_1) err_printf("IFAC: factor %Ps\n\tis prime\n", v);
    2265         [ #  # ]:          0 :   else if (c == gen_0) err_printf("IFAC: factor %Ps\n\tis composite\n", v);
    2266                 :          0 : }
    2267                 :            : static void
    2268                 :          0 : ifac_check(GEN partial, GEN where)
    2269                 :            : {
    2270 [ #  # ][ #  # ]:          0 :   if (!where || where < FIRST(partial) || where > LAST(partial))
                 [ #  # ]
    2271                 :          0 :     pari_err_BUG("ifac_check ['where' out of bounds]");
    2272                 :          0 : }
    2273                 :            : static void
    2274                 :          0 : ifac_print(GEN part, GEN where)
    2275                 :            : {
    2276                 :          0 :   long l = lg(part);
    2277                 :            :   GEN p;
    2278                 :            : 
    2279                 :          0 :   err_printf("ifac partial factorization structure: %ld slots, ", (l-3)/3);
    2280         [ #  # ]:          0 :   if (MOEBIUS(part)) err_printf("Moebius mode, ");
    2281                 :          0 :   err_printf("hint = %ld\n", itos(HINT(part)));
    2282                 :          0 :   ifac_check(part, where);
    2283         [ #  # ]:          0 :   for (p = part+3; p < part + l; p += 3)
    2284                 :            :   {
    2285                 :          0 :     GEN v = VALUE(p), e = EXPON(p), c = CLASS(p);
    2286                 :          0 :     const char *s = "";
    2287         [ #  # ]:          0 :     if (!v) { err_printf("[empty slot]\n"); continue; }
    2288         [ #  # ]:          0 :     if (c == NULL) s = "unknown";
    2289         [ #  # ]:          0 :     else if (c == gen_0) s = "composite";
    2290         [ #  # ]:          0 :     else if (c == gen_1) s = "unfinished prime";
    2291         [ #  # ]:          0 :     else if (c == gen_2) s = "prime";
    2292                 :          0 :     else pari_err_BUG("unknown factor class");
    2293                 :          0 :     err_printf("[%Ps, %Ps, %s]\n", v, e, s);
    2294                 :            :   }
    2295                 :          0 :   err_printf("Done.\n");
    2296                 :          0 : }
    2297                 :            : 
    2298                 :            : static const long decomp_default_hint = 0;
    2299                 :            : /* assume n a non-zero t_INT */
    2300                 :            : /* return initial data structure, see ifac_crack() for the hint argument */
    2301                 :            : static GEN
    2302                 :       3823 : ifac_start_hint(GEN n, int moebius, long hint)
    2303                 :            : {
    2304                 :       3823 :   const long ifac_initial_length = 3 + 7*3;
    2305                 :            :   /* codeword, moebius, hint, 7 slots -- a 512-bit product of distinct 8-bit
    2306                 :            :    * primes needs at most 7 slots at a time) */
    2307                 :       3823 :   GEN here, part = cgetg(ifac_initial_length, t_VEC);
    2308                 :            : 
    2309         [ +  + ]:       3823 :   MOEBIUS(part) = moebius? gen_1 : NULL;
    2310                 :       3823 :   HINT(part) = stoi(hint);
    2311         [ +  + ]:       3823 :   if (isonstack(n)) n = absi(n);
    2312                 :            :   /* make copy, because we'll later want to replace it in place.
    2313                 :            :    * If it's not on stack, then we assume it is a clone made for us by
    2314                 :            :    * ifactor, and we assume the sign has already been set positive */
    2315                 :            :   /* fill first slot at the top end */
    2316                 :       3823 :   here = part + ifac_initial_length - 3; /* LAST(part) */
    2317                 :       3823 :   INIT(here, n,gen_1,gen_0); /* n^1: composite */
    2318         [ +  + ]:      26761 :   while ((here -= 3) > part) ifac_delete(here);
    2319                 :       3823 :   return part;
    2320                 :            : }
    2321                 :            : GEN
    2322                 :       1826 : ifac_start(GEN n, int moebius)
    2323                 :       1826 : { return ifac_start_hint(n,moebius,decomp_default_hint); }
    2324                 :            : 
    2325                 :            : /* Return next nonempty slot after 'here', NULL if none exist */
    2326                 :            : static GEN
    2327                 :      10653 : ifac_find(GEN partial)
    2328                 :            : {
    2329                 :      10653 :   GEN scan, end = partial + lg(partial);
    2330                 :            : 
    2331                 :            : #ifdef IFAC_DEBUG
    2332                 :            :   ifac_check(partial, partial);
    2333                 :            : #endif
    2334         [ +  + ]:      77862 :   for (scan = partial+3; scan < end; scan += 3)
    2335         [ +  + ]:      74508 :     if (VALUE(scan)) return scan;
    2336                 :      10653 :   return NULL;
    2337                 :            : }
    2338                 :            : 
    2339                 :            : /* Defragment: squeeze out unoccupied slots above *where. Unoccupied slots
    2340                 :            :  * arise when a composite factor dissolves completely whilst dividing off a
    2341                 :            :  * prime, or when ifac_resort() spots a coincidence and merges two factors.
    2342                 :            :  * Update *where */
    2343                 :            : static void
    2344                 :        210 : ifac_defrag(GEN *partial, GEN *where)
    2345                 :            : {
    2346                 :        210 :   GEN scan_new = LAST(*partial), scan_old;
    2347                 :            : 
    2348         [ +  + ]:        644 :   for (scan_old = scan_new; scan_old >= *where; scan_old -= 3)
    2349                 :            :   {
    2350         [ -  + ]:        434 :     if (!VALUE(scan_old)) continue; /* empty slot */
    2351         [ -  + ]:        434 :     if (scan_old < scan_new) SHALLOWCOPY(scan_old, scan_new);
    2352                 :        434 :     scan_new -= 3; /* point at next slot to be written */
    2353                 :            :   }
    2354                 :        210 :   scan_new += 3; /* back up to last slot written */
    2355                 :        210 :   *where = scan_new;
    2356         [ +  + ]:       1246 :   while ((scan_new -= 3) > *partial) ifac_delete(scan_new); /* erase junk */
    2357                 :        210 : }
    2358                 :            : 
    2359                 :            : /* Move to a larger main vector, updating *where if it points into it, and
    2360                 :            :  * *partial in any case. Can be used as a specialized gcopy before
    2361                 :            :  * a gerepileupto() (pass 0 as the new length). Normally, one would pass
    2362                 :            :  * new_lg=1 to let this function guess the new size.  To be used sparingly.
    2363                 :            :  * Complex version of ifac_defrag(), combined with reallocation.  If new_lg
    2364                 :            :  * is 0, use the old length, so this acts just like gcopy except that the
    2365                 :            :  * 'where' pointer is carried along; if it is 1, we make an educated guess.
    2366                 :            :  * Exception:  If new_lg is 0, the vector is full to the brim, and the first
    2367                 :            :  * entry is composite, we make it longer to avoid being called again a
    2368                 :            :  * microsecond later. It is safe to call this with *where = NULL:
    2369                 :            :  * if it doesn't point anywhere within the old structure, it is left alone */
    2370                 :            : static void
    2371                 :          0 : ifac_realloc(GEN *partial, GEN *where, long new_lg)
    2372                 :            : {
    2373                 :          0 :   long old_lg = lg(*partial);
    2374                 :            :   GEN newpart, scan_new, scan_old;
    2375                 :            : 
    2376         [ #  # ]:          0 :   if (new_lg == 1)
    2377                 :          0 :     new_lg = 2*old_lg - 6;        /* from 7 slots to 13 to 25... */
    2378         [ #  # ]:          0 :   else if (new_lg <= old_lg)        /* includes case new_lg == 0 */
    2379                 :            :   {
    2380                 :          0 :     GEN first = *partial + 3;
    2381                 :          0 :     new_lg = old_lg;
    2382                 :            :     /* structure full and first entry composite or unknown */
    2383 [ #  # ][ #  # ]:          0 :     if (VALUE(first) && (CLASS(first) == gen_0 || CLASS(first)==NULL))
                 [ #  # ]
    2384                 :          0 :       new_lg += 6; /* give it a little more breathing space */
    2385                 :            :   }
    2386                 :          0 :   newpart = cgetg(new_lg, t_VEC);
    2387         [ #  # ]:          0 :   if (DEBUGMEM >= 3)
    2388                 :          0 :     err_printf("IFAC: new partial factorization structure (%ld slots)\n",
    2389                 :          0 :                (new_lg - 3)/3);
    2390                 :          0 :   MOEBIUS(newpart) = MOEBIUS(*partial);
    2391 [ #  # ][ #  # ]:          0 :   icopyifstack(HINT(*partial), HINT(newpart));
    2392                 :            :   /* Downward sweep through the old *partial. Pick up 'where' and carry it
    2393                 :            :    * over if we pass it. (Only useful if it pointed at a non-empty slot.)
    2394                 :            :    * Factors are COPY'd so that we again have a nice object (parent older
    2395                 :            :    * than children, connected), except the one factor that may still be living
    2396                 :            :    * in a clone where n originally was; exponents are similarly copied if they
    2397                 :            :    * aren't global constants; class-of-factor fields are global constants so we
    2398                 :            :    * need only copy them as pointers. Caller may then do a gerepileupto() */
    2399                 :          0 :   scan_new = newpart + new_lg - 3; /* LAST(newpart) */
    2400                 :          0 :   scan_old = *partial + old_lg - 3; /* LAST(*partial) */
    2401         [ #  # ]:          0 :   for (; scan_old > *partial + 2; scan_old -= 3)
    2402                 :            :   {
    2403         [ #  # ]:          0 :     if (*where == scan_old) *where = scan_new;
    2404         [ #  # ]:          0 :     if (!VALUE(scan_old)) continue; /* skip empty slots */
    2405                 :          0 :     COPY(scan_old, scan_new); scan_new -= 3;
    2406                 :            :   }
    2407                 :          0 :   scan_new += 3; /* back up to last slot written */
    2408         [ #  # ]:          0 :   while ((scan_new -= 3) > newpart) ifac_delete(scan_new);
    2409                 :          0 :   *partial = newpart;
    2410                 :          0 : }
    2411                 :            : 
    2412                 :            : /* Re-sort one (typically unknown) entry from washere to a new position,
    2413                 :            :  * rotating intervening entries upward to fill the vacant space. If the new
    2414                 :            :  * position is the same as the old one, or the new value of the entry coincides
    2415                 :            :  * with a value already occupying a lower slot, then we just add exponents (and
    2416                 :            :  * use the 'more known' class, and return 1 immediately when in Moebius mode).
    2417                 :            :  * Slots between *where and washere must be in sorted order, so a sweep using
    2418                 :            :  * this to re-sort several unknowns must proceed upward, see ifac_resort().
    2419                 :            :  * Bubble-sort-of-thing sort. Won't be exercised frequently, so this is ok */
    2420                 :            : static void
    2421                 :        105 : ifac_sort_one(GEN *where, GEN washere)
    2422                 :            : {
    2423                 :        105 :   GEN old, scan = washere - 3;
    2424                 :            :   GEN value, exponent, class0, class1;
    2425                 :            :   long cmp_res;
    2426                 :            : 
    2427         [ -  + ]:        105 :   if (scan < *where) return; /* nothing to do, washere==*where */
    2428                 :        105 :   value    = VALUE(washere);
    2429                 :        105 :   exponent = EXPON(washere);
    2430                 :        105 :   class0 = CLASS(washere);
    2431                 :        105 :   cmp_res = -1; /* sentinel */
    2432         [ +  - ]:        105 :   while (scan >= *where) /* at least once */
    2433                 :            :   {
    2434         [ +  - ]:        105 :     if (VALUE(scan))
    2435                 :            :     { /* current slot nonempty, check against where */
    2436                 :        105 :       cmp_res = cmpii(value, VALUE(scan));
    2437         [ +  - ]:        105 :       if (cmp_res >= 0) break; /* have found where to stop */
    2438                 :            :     }
    2439                 :            :     /* copy current slot upward by one position and move pointers down */
    2440                 :          0 :     SHALLOWCOPY(scan, scan+3);
    2441                 :          0 :     scan -= 3;
    2442                 :            :   }
    2443                 :        105 :   scan += 3;
    2444                 :            :   /* At this point there are the following possibilities:
    2445                 :            :    * 1) cmp_res == -1. Either value is less than that at *where, or *where was
    2446                 :            :    * pointing at vacant slots and any factors we saw en route were larger than
    2447                 :            :    * value. At any rate, scan == *where now, and scan is pointing at an empty
    2448                 :            :    * slot, into which we'll stash our entry.
    2449                 :            :    * 2) cmp_res == 0. The entry at scan-3 is the one, we compare class0
    2450                 :            :    * fields and add exponents, and put it all into the vacated scan slot,
    2451                 :            :    * NULLing the one at scan-3 (and possibly updating *where).
    2452                 :            :    * 3) cmp_res == 1. The slot at scan is the one to store our entry into. */
    2453         [ +  - ]:        105 :   if (cmp_res)
    2454                 :            :   {
    2455 [ -  + ][ #  # ]:        105 :     if (cmp_res < 0 && scan != *where)
    2456                 :          0 :       pari_err_BUG("ifact_sort_one [misaligned partial]");
    2457                 :        105 :     INIT(scan, value, exponent, class0); return;
    2458                 :            :   }
    2459                 :            :   /* case cmp_res == 0: repeated factor detected */
    2460         [ #  # ]:          0 :   if (DEBUGLEVEL >= 4)
    2461                 :          0 :     err_printf("IFAC: repeated factor %Ps\n\tin ifac_sort_one\n", value);
    2462                 :          0 :   old = scan - 3;
    2463                 :            :   /* if old class0 was composite and new is prime, or vice versa, complain
    2464                 :            :    * (and if one class0 was unknown and the other wasn't, use the known one) */
    2465                 :          0 :   class1 = CLASS(old);
    2466         [ #  # ]:          0 :   if (class0) /* should never be used */
    2467                 :            :   {
    2468         [ #  # ]:          0 :     if (class1)
    2469                 :            :     {
    2470 [ #  # ][ #  # ]:          0 :       if (class0 == gen_0 && class1 != gen_0)
    2471                 :          0 :         pari_err_BUG("ifac_sort_one (composite = prime)");
    2472 [ #  # ][ #  # ]:          0 :       else if (class0 != gen_0 && class1 == gen_0)
    2473                 :          0 :         pari_err_BUG("ifac_sort_one (prime = composite)");
    2474         [ #  # ]:          0 :       else if (class0 == gen_2)
    2475                 :          0 :         CLASS(scan) = class0;
    2476                 :            :     }
    2477                 :            :     else
    2478                 :          0 :       CLASS(scan) = class0;
    2479                 :            :   }
    2480                 :            :   /* else stay with the existing known class0 */
    2481                 :          0 :   CLASS(scan) = class1;
    2482                 :            :   /* in any case, add exponents */
    2483 [ #  # ][ #  # ]:          0 :   if (EXPON(old) == gen_1 && exponent == gen_1)
    2484                 :          0 :     EXPON(scan) = gen_2;
    2485                 :            :   else
    2486                 :          0 :     EXPON(scan) = addii(EXPON(old), exponent);
    2487                 :            :   /* move the value over and null out the vacated slot below */
    2488                 :          0 :   old = scan - 3;
    2489                 :          0 :   *scan = *old;
    2490                 :          0 :   ifac_delete(old);
    2491                 :            :   /* finally, see whether *where should be pulled in */
    2492         [ #  # ]:        105 :   if (old == *where) *where += 3;
    2493                 :            : }
    2494                 :            : 
    2495                 :            : /* Sort all current unknowns downward to where they belong. Sweeps in the
    2496                 :            :  * upward direction. Not needed after ifac_crack(), only when ifac_divide()
    2497                 :            :  * returned true. Update *where. */
    2498                 :            : static void
    2499                 :        105 : ifac_resort(GEN *partial, GEN *where)
    2500                 :            : {
    2501                 :            :   GEN scan, end;
    2502                 :        105 :   ifac_defrag(partial, where); end = LAST(*partial);
    2503         [ +  + ]:        322 :   for (scan = *where; scan <= end; scan += 3)
    2504 [ +  - ][ +  + ]:        217 :     if (VALUE(scan) && !CLASS(scan)) ifac_sort_one(where, scan); /*unknown*/
    2505                 :        105 :   ifac_defrag(partial, where); /* remove newly created gaps */
    2506                 :        105 : }
    2507                 :            : 
    2508                 :            : /* Let x be a t_INT known not to have small divisors (< 2^14). Return 0 if x
    2509                 :            :  * is a proven composite. Return 1 if we believe it to be prime (fully proven
    2510                 :            :  * prime if factor_proven is set).  */
    2511                 :            : int
    2512                 :      11607 : ifac_isprime(GEN x)
    2513                 :            : {
    2514         [ +  + ]:      11607 :   if (!BPSW_psp_nosmalldiv(x)) return 0; /* composite */
    2515 [ -  + ][ #  # ]:       9270 :   if (factor_proven && ! BPSW_isprime(x))
    2516                 :            :   {
    2517                 :          0 :     pari_warn(warner,
    2518                 :            :               "IFAC: pseudo-prime %Ps\n\tis not prime. PLEASE REPORT!\n", x);
    2519                 :          0 :     return 0;
    2520                 :            :   }
    2521                 :      11607 :   return 1;
    2522                 :            : }
    2523                 :            : 
    2524                 :            : static int
    2525                 :       8023 : ifac_checkprime(GEN x)
    2526                 :            : {
    2527                 :       8023 :   int res = ifac_isprime(VALUE(x));
    2528         [ +  + ]:       8023 :   CLASS(x) = res? gen_1: gen_0;
    2529         [ -  + ]:       8023 :   if (DEBUGLEVEL>2) ifac_factor_dbg(x);
    2530                 :       8023 :   return res;
    2531                 :            : }
    2532                 :            : 
    2533                 :            : /* Determine primality or compositeness of all current unknowns, and set
    2534                 :            :  * class Q primes to finished (class P) if everything larger is already
    2535                 :            :  * known to be prime.  When after_crack >= 0, only look at the
    2536                 :            :  * first after_crack things in the list (do nothing when it's 0) */
    2537                 :            : static void
    2538                 :       4130 : ifac_whoiswho(GEN *partial, GEN *where, long after_crack)
    2539                 :            : {
    2540                 :       4130 :   GEN scan, scan_end = LAST(*partial);
    2541                 :            : 
    2542                 :            : #ifdef IFAC_DEBUG
    2543                 :            :   ifac_check(*partial, *where);
    2544                 :            : #endif
    2545         [ +  + ]:       7988 :   if (after_crack == 0) return;
    2546         [ +  + ]:       3858 :   if (after_crack > 0) /* check at most after_crack entries */
    2547                 :       3753 :     scan = *where + 3*(after_crack - 1); /* assert(scan <= scan_end) */
    2548                 :            :   else
    2549         [ +  + ]:        301 :     for (scan = scan_end; scan >= *where; scan -= 3)
    2550                 :            :     {
    2551         [ +  + ]:        203 :       if (CLASS(scan))
    2552                 :            :       { /* known class of factor */
    2553         [ +  + ]:        105 :         if (CLASS(scan) == gen_0) break;
    2554         [ -  + ]:         98 :         if (CLASS(scan) == gen_1)
    2555                 :            :         {
    2556         [ #  # ]:          0 :           if (DEBUGLEVEL>=3)
    2557                 :            :           {
    2558                 :          0 :             err_printf("IFAC: factor %Ps\n\tis prime (no larger composite)\n",
    2559                 :          0 :                        VALUE(*where));
    2560                 :          0 :             err_printf("IFAC: prime %Ps\n\tappears with exponent = %ld\n",
    2561                 :          0 :                        VALUE(*where), itos(EXPON(*where)));
    2562                 :            :           }
    2563                 :          0 :           CLASS(scan) = gen_2;
    2564                 :            :         }
    2565                 :         98 :         continue;
    2566                 :            :       }
    2567         [ -  + ]:         98 :       if (!ifac_checkprime(scan)) break; /* must disable Q-to-P */
    2568                 :         98 :       CLASS(scan) = gen_2; /* P_i, finished prime */
    2569         [ -  + ]:         98 :       if (DEBUGLEVEL>2) ifac_factor_dbg(scan);
    2570                 :            :     }
    2571                 :            :   /* go on, Q-to-P trick now disabled */
    2572         [ +  + ]:      11469 :   for (; scan >= *where; scan -= 3)
    2573                 :            :   {
    2574         [ +  + ]:       7611 :     if (CLASS(scan)) continue;
    2575                 :       7597 :     (void)ifac_checkprime(scan); /* Qj | Ck */
    2576                 :            :   }
    2577                 :            : }
    2578                 :            : 
    2579                 :            : /* Divide all current composites by first (prime, class Q) entry, updating its
    2580                 :            :  * exponent, and turning it into a finished prime (class P).  Return 1 if any
    2581                 :            :  * such divisions succeeded  (in Moebius mode, the update may then not have
    2582                 :            :  * been completed), or 0 if none of them succeeded.  Doesn't modify *where.
    2583                 :            :  * Here we normally do not check that the first entry is a not-finished
    2584                 :            :  * prime.  Stack management: we may allocate a new exponent */
    2585                 :            : static long
    2586                 :       7148 : ifac_divide(GEN *partial, GEN *where, long moebius_mode)
    2587                 :            : {
    2588                 :       7148 :   GEN scan, scan_end = LAST(*partial);
    2589                 :       7148 :   long res = 0, exponent, newexp, otherexp;
    2590                 :            : 
    2591                 :            : #ifdef IFAC_DEBUG
    2592                 :            :   ifac_check(*partial, *where);
    2593                 :            :   if (CLASS(*where) != gen_1)
    2594                 :            :     pari_err_BUG("ifac_divide [division by composite or finished prime]");
    2595                 :            :   if (!VALUE(*where)) pari_err_BUG("ifac_divide [division by nothing]");
    2596                 :            : #endif
    2597                 :       7148 :   newexp = exponent = itos(EXPON(*where));
    2598 [ +  + ][ -  + ]:       7148 :   if (exponent > 1 && moebius_mode) return 1;
    2599                 :            :   /* should've been caught by caller */
    2600                 :            : 
    2601         [ +  + ]:      11041 :   for (scan = *where+3; scan <= scan_end; scan += 3)
    2602                 :            :   {
    2603         [ +  + ]:       3900 :     if (CLASS(scan) != gen_0) continue; /* the other thing ain't composite */
    2604                 :        364 :     otherexp = 0;
    2605                 :            :     /* divide in place to keep stack clutter minimal */
    2606         [ +  + ]:        469 :     while (dvdiiz(VALUE(scan), VALUE(*where), VALUE(scan)))
    2607                 :            :     {
    2608         [ +  + ]:        112 :       if (moebius_mode) return 1; /* immediately */
    2609         [ +  - ]:        105 :       if (!otherexp) otherexp = itos(EXPON(scan));
    2610                 :        105 :       newexp += otherexp;
    2611                 :            :     }
    2612         [ +  + ]:        357 :     if (newexp > exponent)        /* did anything happen? */
    2613                 :            :     {
    2614         [ -  + ]:        105 :       EXPON(*where) = (newexp == 2 ? gen_2 : utoipos(newexp));
    2615                 :        105 :       exponent = newexp;
    2616         [ -  + ]:        105 :       if (is_pm1((GEN)*scan)) /* factor dissolved completely */
    2617                 :            :       {
    2618                 :          0 :         ifac_delete(scan);
    2619         [ #  # ]:          0 :         if (DEBUGLEVEL >= 4)
    2620                 :          0 :           err_printf("IFAC: a factor was a power of another prime factor\n");
    2621                 :            :       } else {
    2622                 :        105 :         CLASS(scan) = NULL;        /* at any rate it's Unknown now */
    2623         [ -  + ]:        105 :         if (DEBUGLEVEL >= 4)
    2624                 :          0 :           err_printf("IFAC: a factor was divisible by another prime factor,\n"
    2625                 :            :                      "\tleaving a cofactor = %Ps\n", VALUE(scan));
    2626                 :            :       }
    2627                 :        105 :       res = 1;
    2628         [ -  + ]:        105 :       if (DEBUGLEVEL >= 5)
    2629                 :          0 :         err_printf("IFAC: prime %Ps\n\tappears at least to the power %ld\n",
    2630                 :          0 :                    VALUE(*where), newexp);
    2631                 :            :     }
    2632                 :            :   } /* for */
    2633                 :       7141 :   CLASS(*where) = gen_2; /* make it a finished prime */
    2634         [ -  + ]:       7141 :   if (DEBUGLEVEL >= 3)
    2635                 :          0 :     err_printf("IFAC: prime %Ps\n\tappears with exponent = %ld\n",
    2636                 :          0 :                VALUE(*where), newexp);
    2637                 :       7148 :   return res;
    2638                 :            : }
    2639                 :            : 
    2640                 :            : /* found out our integer was factor^exp. Update */
    2641                 :            : static void
    2642                 :        377 : update_pow(GEN where, GEN factor, long exp, pari_sp *av)
    2643                 :            : {
    2644                 :        377 :   GEN ex = EXPON(where);
    2645         [ -  + ]:        377 :   if (DEBUGLEVEL>3)
    2646                 :          0 :     err_printf("IFAC: found %Ps =\n\t%Ps ^%ld\n", *where, factor, exp);
    2647                 :        377 :   affii(factor, VALUE(where)); avma = *av;
    2648         [ +  + ]:        377 :   if (ex == gen_1)
    2649         [ +  + ]:        335 :   { EXPON(where) = exp == 2? gen_2: utoipos(exp); *av = avma; }
    2650         [ +  + ]:         42 :   else if (ex == gen_2)
    2651                 :         28 :   { EXPON(where) = utoipos(exp<<1); *av = avma; }
    2652                 :            :   else
    2653                 :         14 :     affsi(exp * itos(ex), EXPON(where));
    2654                 :        377 : }
    2655                 :            : /* hint == 0 : Use a default strategy
    2656                 :            :  * hint & 1  : Avoid mpqs(), use ellfacteur() after pollardbrent()
    2657                 :            :  * hint & 2  : Avoid first-stage ellfacteur() in favour of mpqs()
    2658                 :            :  * (may still fall back to ellfacteur() if mpqs() is not installed or gives up)
    2659                 :            :  * hint & 4  : Avoid even the pollardbrent() and squfof() stages. Put under
    2660                 :            :  *  the same governing  bit, for no good reason other than avoiding a
    2661                 :            :  *  proliferation of bits.
    2662                 :            :  * hint & 8  : Avoid final ellfacteur(); this may declare a composite to be
    2663                 :            :  *  prime.  */
    2664                 :            : #define get_hint(partial) (itos(HINT(*partial)) & 15)
    2665                 :            : 
    2666                 :            : /* Split the first (composite) entry.  There _must_ already be room for another
    2667                 :            :  * factor below *where, and *where is updated. Two cases:
    2668                 :            :  * - entry = factor^k is a pure power: factor^k is inserted, leaving *where
    2669                 :            :  *   unchanged;
    2670                 :            :  * - entry = factor * cofactor (not necessarily coprime): both factors are
    2671                 :            :  *   inserted in the correct order, updating *where
    2672                 :            :  * The inserted factors class is set to unknown, they inherit the exponent
    2673                 :            :  * (or a multiple thereof) of their ancestor.
    2674                 :            :  *
    2675                 :            :  * Returns number of factors written into the structure, normally 2 (1 if pure
    2676                 :            :  * power, maybe > 2 if a factoring engine returned a vector of factors instead
    2677                 :            :  * of a single factor). Can reallocate the data structure in the
    2678                 :            :  * vector-of-factors case, not in the most common single-factor case.
    2679                 :            :  * Stack housekeeping:  this routine may create one or more objects  (a new
    2680                 :            :  * factor, or possibly several, and perhaps one or more new exponents > 2) */
    2681                 :            : static long
    2682                 :       4032 : ifac_crack(GEN *partial, GEN *where, long moebius_mode)
    2683                 :            : {
    2684                 :       4032 :   long cmp_res, hint = get_hint(partial);
    2685                 :            :   GEN factor, exponent;
    2686                 :            : 
    2687                 :            : #ifdef IFAC_DEBUG
    2688                 :            :   ifac_check(*partial, *where);
    2689                 :            :   if (*where < *partial + 6)
    2690                 :            :     pari_err_BUG("ifac_crack ['*where' out of bounds]");
    2691                 :            :   if (!(VALUE(*where)) || typ(VALUE(*where)) != t_INT)
    2692                 :            :     pari_err_BUG("ifac_crack [incorrect VALUE(*where)]");
    2693                 :            :   if (CLASS(*where) != gen_0)
    2694                 :            :     pari_err_BUG("ifac_crack [operand not known composite]");
    2695                 :            : #endif
    2696                 :            : 
    2697         [ -  + ]:       4032 :   if (DEBUGLEVEL>2) {
    2698                 :          0 :     err_printf("IFAC: cracking composite\n\t%Ps\n", **where);
    2699         [ #  # ]:          0 :     if (DEBUGLEVEL>3) err_printf("IFAC: checking for pure square\n");
    2700                 :            :   }
    2701                 :            :   /* MPQS cannot factor prime powers. Look for pure powers even if MPQS is
    2702                 :            :    * blocked by hint: fast and useful in bounded factorization */
    2703                 :            :   {
    2704                 :            :     forprime_t T;
    2705                 :       4032 :     ulong exp = 1, mask = 7;
    2706                 :       4032 :     long good = 0;
    2707                 :       4032 :     pari_sp av = avma;
    2708                 :       4032 :     (void)u_forprime_init(&T, 11, ULONG_MAX);
    2709                 :            :     /* crack squares */
    2710         [ +  + ]:       4318 :     while (Z_issquareall(VALUE(*where), &factor))
    2711                 :            :     {
    2712                 :        293 :       good = 1; /* remember we succeeded once */
    2713                 :        293 :       update_pow(*where, factor, 2, &av);
    2714         [ +  + ]:        293 :       if (moebius_mode) return 0; /* no need to carry on */
    2715                 :            :     }
    2716         [ +  + ]:       4109 :     while ( (exp = is_357_power(VALUE(*where), &factor, &mask)) )
    2717                 :            :     {
    2718                 :         84 :       good = 1; /* remember we succeeded once */
    2719                 :         84 :       update_pow(*where, factor, exp, &av);
    2720         [ -  + ]:         84 :       if (moebius_mode) return 0; /* no need to carry on */
    2721                 :            :     }
    2722                 :            :     /* cutoff at 14 bits as trial division must have found everything below */
    2723         [ -  + ]:       4025 :     while ( (exp = is_pth_power(VALUE(*where), &factor, &T, 15)) )
    2724                 :            :     {
    2725                 :          0 :       good = 1; /* remember we succeeded once */
    2726                 :          0 :       update_pow(*where, factor, exp, &av);
    2727         [ #  # ]:          0 :       if (moebius_mode) return 0; /* no need to carry on */
    2728                 :            :     }
    2729                 :            : 
    2730 [ +  + ][ +  - ]:       4025 :     if (good && hint != 15 && ifac_checkprime(*where))
                 [ +  + ]
    2731                 :            :     { /* our composite was a prime power */
    2732         [ -  + ]:        272 :       if (DEBUGLEVEL>3)
    2733                 :          0 :         err_printf("IFAC: factor %Ps\n\tis prime\n", VALUE(*where));
    2734                 :        272 :       return 0; /* bypass subsequent ifac_whoiswho() call */
    2735                 :            :     }
    2736                 :            :   } /* pure power stage */
    2737                 :            : 
    2738                 :       3753 :   factor = NULL;
    2739         [ +  - ]:       4032 :   if (!(hint & 4))
    2740                 :            :   { /* pollardbrent() Rho usually gets a first chance */
    2741         [ -  + ]:       3753 :     if (DEBUGLEVEL >= 4) err_printf("IFAC: trying Pollard-Brent rho method\n");
    2742                 :       3753 :     factor = pollardbrent(VALUE(*where));
    2743         [ +  + ]:       3753 :     if (!factor)
    2744                 :            :     { /* Shanks' squfof() */
    2745         [ -  + ]:       2255 :       if (DEBUGLEVEL >= 4)
    2746                 :          0 :         err_printf("IFAC: trying Shanks' SQUFOF, will fail silently if input\n"
    2747                 :            :                    "      is too large for it.\n");
    2748                 :       2255 :       factor = squfof(VALUE(*where));
    2749                 :            :     }
    2750                 :            :   }
    2751 [ +  + ][ +  - ]:       3753 :   if (!factor && !(hint & 2))
    2752                 :            :   { /* First ECM stage */
    2753         [ -  + ]:        560 :     if (DEBUGLEVEL >= 4) err_printf("IFAC: trying Lenstra-Montgomery ECM\n");
    2754                 :        560 :     factor = ellfacteur(VALUE(*where), 0); /* do not insist */
    2755                 :            :   }
    2756 [ +  + ][ +  + ]:       3753 :   if (!factor && !(hint & 1))
    2757                 :            :   { /* MPQS stage */
    2758         [ -  + ]:        539 :     if (DEBUGLEVEL >= 4) err_printf("IFAC: trying MPQS\n");
    2759                 :        539 :     factor = mpqs(VALUE(*where));
    2760                 :            :   }
    2761         [ +  + ]:       3753 :   if (!factor)
    2762                 :            :   {
    2763         [ +  - ]:         14 :     if (!(hint & 8))
    2764                 :            :     { /* still no luck? Final ECM stage, guaranteed to succeed */
    2765         [ -  + ]:         14 :       if (DEBUGLEVEL >= 4)
    2766                 :          0 :         err_printf("IFAC: forcing ECM, may take some time\n");
    2767                 :         14 :       factor = ellfacteur(VALUE(*where), 1);
    2768                 :            :     }
    2769                 :            :     else
    2770                 :            :     { /* limited factorization */
    2771         [ #  # ]:          0 :       if (DEBUGLEVEL >= 2)
    2772                 :            :       {
    2773         [ #  # ]:          0 :         if (hint != 15)
    2774                 :          0 :           pari_warn(warner, "IFAC: unfactored composite declared prime");
    2775                 :            :         else
    2776                 :          0 :           pari_warn(warner, "IFAC: untested integer declared prime");
    2777                 :            : 
    2778                 :            :         /* don't print it out at level 3 or above, where it would appear
    2779                 :            :          * several times before and after this message already */
    2780         [ #  # ]:          0 :         if (DEBUGLEVEL == 2) err_printf("\t%Ps\n", VALUE(*where));
    2781                 :            :       }
    2782                 :          0 :       CLASS(*where) = gen_1; /* might as well trial-divide by it... */
    2783                 :          0 :       return 1;
    2784                 :            :     }
    2785                 :            :   }
    2786         [ +  + ]:       3753 :   if (typ(factor) == t_VEC) /* delegate this case */
    2787                 :        679 :     return ifac_insert_multiplet(partial, where, factor, moebius_mode);
    2788                 :            :   /* typ(factor) == t_INT */
    2789                 :            :   /* got single integer back:  work out the cofactor (in place) */
    2790         [ -  + ]:       3074 :   if (!dvdiiz(VALUE(*where), factor, VALUE(*where)))
    2791                 :            :   {
    2792                 :          0 :     err_printf("IFAC: factoring %Ps\n", VALUE(*where));
    2793                 :          0 :     err_printf("\tyielded 'factor' %Ps\n\twhich isn't!\n", factor);
    2794                 :          0 :     pari_err_BUG("factoring");
    2795                 :            :   }
    2796                 :            :   /* factoring engines report the factor found; tell about the cofactor */
    2797         [ -  + ]:       3074 :   if (DEBUGLEVEL >= 4) err_printf("IFAC: cofactor = %Ps\n", VALUE(*where));
    2798                 :            : 
    2799                 :            :   /* The two factors are 'factor' and VALUE(*where), find out which is larger */
    2800                 :       3074 :   cmp_res = cmpii(factor, VALUE(*where));
    2801                 :       3074 :   CLASS(*where) = NULL; /* mark factor /cofactor 'unknown' */
    2802                 :       3074 :   exponent = EXPON(*where);
    2803                 :       3074 :   *where -= 3;
    2804                 :       3074 :   CLASS(*where) = NULL; /* mark factor /cofactor 'unknown' */
    2805         [ +  + ]:       3074 :   EXPON(*where) = isonstack(exponent)? icopy(exponent): exponent;
    2806         [ +  + ]:       3074 :   if (cmp_res < 0)
    2807                 :       2892 :     VALUE(*where) = factor; /* common case */
    2808         [ +  - ]:        182 :   else if (cmp_res > 0)
    2809                 :            :   { /* factor > cofactor, rearrange */
    2810                 :        182 :     GEN old = *where + 3;
    2811                 :        182 :     VALUE(*where) = VALUE(old); /* move cofactor pointer to lowest slot */
    2812                 :        182 :     VALUE(old) = factor; /* save factor */
    2813                 :            :   }
    2814                 :          0 :   else pari_err_BUG("ifac_crack [Z_issquareall miss]");
    2815                 :       4032 :   return 2;
    2816                 :            : }
    2817                 :            : 
    2818                 :            : /* Gets called to complete ifac_crack's job when a factoring engine splits
    2819                 :            :  * the current factor into a product of three or more new factors. Makes room
    2820                 :            :  * for them if necessary, sorts them, gives them the right exponents and class.
    2821                 :            :  * Also returns the number of factors actually written, which may be less than
    2822                 :            :  * the number of components in facvec if there are duplicates.--- Vectors of
    2823                 :            :  * factors  (cf pollardbrent()) actually contain 'slots' of three GENs per
    2824                 :            :  * factor with the three fields interpreted as in our partial factorization
    2825                 :            :  * data structure.  Thus 'engines' can tell us what they already happen to
    2826                 :            :  * know about factors being prime or composite and/or appearing to a power
    2827                 :            :  * larger than the first.
    2828                 :            :  * Don't collect garbage.  No diagnostics: the factoring engine should have
    2829                 :            :  * printed what it found. facvec contains slots of three components per factor;
    2830                 :            :  * repeated factors are allowed  (and their classes shouldn't contradict each
    2831                 :            :  * other whereas their exponents will be added up) */
    2832                 :            : static long
    2833                 :        679 : ifac_insert_multiplet(GEN *partial, GEN *where, GEN facvec, long moebius_mode)
    2834                 :            : {
    2835                 :        679 :   long j,k=1, lfv=lg(facvec)-1, nf=lfv/3, room=(long)(*where-*partial);
    2836                 :            :   /* one of the factors will go into the *where slot, so room is now 3 times
    2837                 :            :    * the number of slots we can use */
    2838                 :        679 :   long needroom = lfv - room;
    2839                 :        679 :   GEN e, newexp, cur, sorted, auxvec = cgetg(nf+1, t_VEC), factor;
    2840                 :        679 :   long exponent = itos(EXPON(*where)); /* the old exponent */
    2841                 :            : 
    2842         [ -  + ]:        679 :   if (DEBUGLEVEL >= 5) /* squfof may return a single squared factor as a set */
    2843                 :          0 :     err_printf("IFAC: incorporating set of %ld factor(s)\n", nf);
    2844         [ -  + ]:        679 :   if (needroom > 0) /* one extra slot for paranoia, errm, future use */
    2845                 :          0 :     ifac_realloc(partial, where, lg(*partial) + needroom + 3);
    2846                 :            : 
    2847                 :            :   /* create sort permutation from the values of the factors */
    2848         [ +  + ]:       2128 :   for (j=nf; j; j--) auxvec[j] = facvec[3*j-2]; /* just the pointers */
    2849                 :        679 :   sorted = indexsort(auxvec);
    2850                 :            :   /* and readjust the result for the triple spacing */
    2851         [ +  + ]:       2128 :   for (j=nf; j; j--) sorted[j] = 3*sorted[j]-2;
    2852                 :            : 
    2853                 :            :   /* store factors, beginning at *where, and catching any duplicates */
    2854                 :        679 :   cur = facvec + sorted[nf];
    2855                 :        679 :   VALUE(*where) = VALUE(cur);
    2856                 :        679 :   newexp = EXPON(cur);
    2857         [ -  + ]:        679 :   if (newexp != gen_1) /* new exponent > 1 */
    2858                 :            :   {
    2859         [ #  # ]:          0 :     if (exponent == 1)
    2860         [ #  # ]:          0 :       e = isonstack(newexp)? icopy(newexp): newexp;
    2861                 :            :     else
    2862                 :          0 :       e = mului(exponent, newexp);
    2863                 :          0 :     EXPON(*where) = e;
    2864                 :            :   } /* if new exponent is 1, the old exponent already in place will do */
    2865                 :        679 :   CLASS(*where) = CLASS(cur);
    2866         [ -  + ]:        679 :   if (DEBUGLEVEL >= 6) err_printf("\tstored (largest) factor no. %ld...\n", nf);
    2867                 :            : 
    2868         [ +  + ]:       1449 :   for (j=nf-1; j; j--)
    2869                 :            :   {
    2870                 :        770 :     cur = facvec + sorted[j];
    2871                 :        770 :     factor = VALUE(cur);
    2872         [ +  + ]:        770 :     if (equalii(factor, VALUE(*where)))
    2873                 :            :     {
    2874         [ -  + ]:          7 :       if (DEBUGLEVEL >= 6)
    2875         [ #  # ]:          0 :         err_printf("\tfactor no. %ld is a duplicate%s\n", j, (j>1? "...": ""));
    2876                 :            :       /* update exponent, ignore class which would already have been set,
    2877                 :            :        * then forget current factor */
    2878                 :          7 :       newexp = EXPON(cur);
    2879         [ -  + ]:          7 :       if (newexp != gen_1) /* new exp > 1 */
    2880                 :          0 :         e = addis(EXPON(*where), exponent * itos(newexp));
    2881 [ -  + ][ #  # ]:          7 :       else if (EXPON(*where) == gen_1 && exponent == 1)
    2882                 :          0 :         e = gen_2;
    2883                 :            :       else
    2884                 :          7 :         e = addis(EXPON(*where), exponent);
    2885                 :          7 :       EXPON(*where) = e;
    2886                 :            : 
    2887         [ -  + ]:          7 :       if (moebius_mode) return 0; /* stop now, but with exponent updated */
    2888                 :          7 :       continue;
    2889                 :            :     }
    2890                 :            : 
    2891                 :        763 :     *where -= 3;
    2892                 :        763 :     CLASS(*where) = CLASS(cur);        /* class as given */
    2893                 :        763 :     newexp = EXPON(cur);
    2894         [ -  + ]:        763 :     if (newexp != gen_1) /* new exp > 1 */
    2895                 :            :     {
    2896 [ #  # ][ #  # ]:          0 :       if (exponent == 1 && newexp == gen_2)
    2897                 :          0 :         e = gen_2;
    2898                 :            :       else /* exponent*newexp > 2 */
    2899                 :          0 :         e = mului(exponent, newexp);
    2900                 :            :     }
    2901                 :            :     else
    2902         [ +  + ]:        819 :       e = (exponent == 1 ? gen_1 :
    2903         [ +  + ]:         56 :             (exponent == 2 ? gen_2 :
    2904                 :         14 :                utoipos(exponent))); /* inherit parent's exponent */
    2905                 :        763 :     EXPON(*where) = e;
    2906                 :            :     /* keep components younger than *partial */
    2907         [ +  - ]:        763 :     VALUE(*where) = isonstack(factor) ? icopy(factor) : factor;
    2908                 :        763 :     k++;
    2909         [ -  + ]:        763 :     if (DEBUGLEVEL >= 6)
    2910         [ #  # ]:          0 :       err_printf("\tfactor no. %ld was unique%s\n", j, j>1? " (so far)...": "");
    2911                 :            :   }
    2912                 :            :   /* make the 'sorted' object safe for garbage collection (it should be in the
    2913                 :            :    * garbage zone from everybody's perspective, but it's easy to do it) */
    2914                 :        679 :   *sorted = evaltyp(t_INT) | evallg(nf+1);
    2915                 :        679 :   return k;
    2916                 :            : }
    2917                 :            : 
    2918                 :            : /* main loop:  iterate until smallest entry is a finished prime;  returns
    2919                 :            :  * a 'where' pointer, or NULL if nothing left, or gen_0 in Moebius mode if
    2920                 :            :  * we aren't squarefree */
    2921                 :            : static GEN
    2922                 :      10563 : ifac_main(GEN *partial)
    2923                 :            : {
    2924                 :      10563 :   const long moebius_mode = !!MOEBIUS(*partial);
    2925                 :      10563 :   GEN here = ifac_find(*partial);
    2926                 :            :   long nf;
    2927                 :            : 
    2928         [ +  + ]:      10563 :   if (!here) return NULL; /* nothing left */
    2929                 :            :   /* loop until first entry is a finished prime.  May involve reallocations,
    2930                 :            :    * thus updates of *partial */
    2931         [ +  + ]:      18405 :   while (CLASS(here) != gen_2)
    2932                 :            :   {
    2933         [ +  + ]:      11180 :     if (CLASS(here) == gen_0) /* composite: crack it */
    2934                 :            :     { /* make sure there's room for another factor */
    2935         [ -  + ]:       4032 :       if (here < *partial + 6)
    2936                 :            :       {
    2937                 :          0 :         ifac_defrag(partial, &here);
    2938         [ #  # ]:          0 :         if (here < *partial + 6) ifac_realloc(partial, &here, 1); /* no luck */
    2939                 :            :       }
    2940                 :       4032 :       nf = ifac_crack(partial, &here, moebius_mode);
    2941 [ +  + ][ +  + ]:       4032 :       if (moebius_mode && EXPON(here) != gen_1) /* that was a power */
    2942                 :            :       {
    2943         [ -  + ]:          7 :         if (DEBUGLEVEL >= 3)
    2944                 :          0 :           err_printf("IFAC: main loop: repeated new factor\n\t%Ps\n", *here);
    2945                 :          7 :         return gen_0;
    2946                 :            :       }
    2947                 :            :       /* deal with the new unknowns.  No sort: ifac_crack did it */
    2948                 :       4025 :       ifac_whoiswho(partial, &here, nf);
    2949                 :       4025 :       continue;
    2950                 :            :     }
    2951         [ +  - ]:       7148 :     if (CLASS(here) == gen_1) /* prime but not yet finished: finish it */
    2952                 :            :     {
    2953         [ +  + ]:       7148 :       if (ifac_divide(partial, &here, moebius_mode))
    2954                 :            :       {
    2955         [ +  + ]:        112 :         if (moebius_mode)
    2956                 :            :         {
    2957         [ -  + ]:          7 :           if (DEBUGLEVEL >= 3)
    2958                 :          0 :             err_printf("IFAC: main loop: another factor was divisible by\n"
    2959                 :            :                        "\t%Ps\n", *here);
    2960                 :          7 :           return gen_0;
    2961                 :            :         }
    2962                 :        105 :         ifac_resort(partial, &here); /* sort new cofactors down */
    2963                 :        105 :         ifac_whoiswho(partial, &here, -1);
    2964                 :            :       }
    2965                 :       7141 :       continue;
    2966                 :            :     }
    2967                 :          0 :     pari_err_BUG("ifac_main [non-existent factor class]");
    2968                 :            :   } /* while */
    2969 [ +  + ][ -  + ]:       7225 :   if (moebius_mode && EXPON(here) != gen_1)
    2970                 :            :   {
    2971         [ #  # ]:          0 :     if (DEBUGLEVEL >= 3)
    2972                 :          0 :       err_printf("IFAC: after main loop: repeated old factor\n\t%Ps\n", *here);
    2973                 :          0 :     return gen_0;
    2974                 :            :   }
    2975         [ -  + ]:       7225 :   if (DEBUGLEVEL >= 4)
    2976                 :            :   {
    2977                 :          0 :     nf = (*partial + lg(*partial) - here - 3)/3;
    2978         [ #  # ]:          0 :     if (nf)
    2979         [ #  # ]:          0 :       err_printf("IFAC: main loop: %ld factor%s left\n", nf, (nf>1)? "s": "");
    2980                 :            :     else
    2981                 :          0 :       err_printf("IFAC: main loop: this was the last factor\n");
    2982                 :            :   }
    2983 [ -  + ][ #  # ]:       7225 :   if (factor_add_primes && !(get_hint(partial) & 8))
    2984                 :            :   {
    2985                 :          0 :     GEN p = VALUE(here);
    2986 [ #  # ][ #  # ]:          0 :     if (lgefint(p)>3 || uel(p,2) > 0x1000000UL) (void)addprimes(p);
    2987                 :            :   }
    2988                 :      10563 :   return here;
    2989                 :            : }
    2990                 :            : 
    2991                 :            : /* Encapsulated routines */
    2992                 :            : 
    2993                 :            : /* prime/exponent pairs need to appear contiguously on the stack, but we also
    2994                 :            :  * need our data structure somewhere, and we don't know in advance how many
    2995                 :            :  * primes will turn up.  The following discipline achieves this:  When
    2996                 :            :  * ifac_decomp() is called, n should point at an object older than the oldest
    2997                 :            :  * small prime/exponent pair  (ifactor() guarantees this).
    2998                 :            :  * We allocate sufficient space to accommodate several pairs -- eleven pairs
    2999                 :            :  * ought to fit in a space not much larger than n itself -- before calling
    3000                 :            :  * ifac_start().  If we manage to complete the factorization before we run out
    3001                 :            :  * of space, we free the data structure and cull the excess reserved space
    3002                 :            :  * before returning.  When we do run out, we have to leapfrog to generate more
    3003                 :            :  * (guesstimating the requirements from what is left in the partial
    3004                 :            :  * factorization structure);  room for fresh pairs is allocated at the head of
    3005                 :            :  * the stack, followed by an ifac_realloc() to reconnect the data structure and
    3006                 :            :  * move it out of the way, followed by a few pointer tweaks to connect the new
    3007                 :            :  * pairs space to the old one. This whole affair translates into a surprisingly
    3008                 :            :  * compact routine. */
    3009                 :            : 
    3010                 :            : /* find primary factors of n */
    3011                 :            : static long
    3012                 :       1131 : ifac_decomp(GEN n, long hint)
    3013                 :            : {
    3014                 :       1131 :   pari_sp av = avma, lim = stack_lim(av, 1);
    3015                 :       1131 :   long nb = 0;
    3016                 :       1131 :   GEN part, here, workspc, pairs = (GEN)av;
    3017                 :            : 
    3018                 :            :   /* workspc will be doled out in pairs of smaller t_INTs. For n = prod p^{e_p}
    3019                 :            :    * (p not necessarily prime), need room to store all p and e_p [ cgeti(3) ],
    3020                 :            :    * bounded by
    3021                 :            :    *    sum_{p | n} ( log_{2^BIL} (p) + 6 ) <= log_{2^BIL} n + 6 log_2 n */
    3022                 :       1131 :   workspc = new_chunk((expi(n) + 1) * 7);
    3023                 :       1131 :   part = ifac_start_hint(n, 0, hint);
    3024                 :            :   for (;;)
    3025                 :            :   {
    3026                 :       3455 :     here = ifac_main(&part);
    3027         [ +  + ]:       3455 :     if (!here) break;
    3028         [ -  + ]:       2324 :     if (low_stack(lim, stack_lim(av,1)))
    3029                 :            :     {
    3030                 :            :       long offset;
    3031         [ #  # ]:          0 :       if(DEBUGMEM>1)
    3032                 :            :       {
    3033                 :          0 :         pari_warn(warnmem,"[2] ifac_decomp");
    3034                 :          0 :         ifac_print(part, here);
    3035                 :            :       }
    3036                 :          0 :       ifac_realloc(&part, &here, 0);
    3037                 :          0 :       offset = here - part;
    3038                 :          0 :       part = gerepileupto((pari_sp)workspc, part);
    3039                 :          0 :       here = part + offset;
    3040                 :            :     }
    3041                 :       2324 :     nb++;
    3042                 :       2324 :     pairs = icopy_avma(VALUE(here), (pari_sp)pairs);
    3043                 :       2324 :     pairs = icopy_avma(EXPON(here), (pari_sp)pairs);
    3044                 :       2324 :     ifac_delete(here);
    3045                 :       2324 :   }
    3046                 :       1131 :   avma = (pari_sp)pairs;
    3047         [ -  + ]:       1131 :   if (DEBUGLEVEL >= 3)
    3048         [ #  # ]:          0 :     err_printf("IFAC: found %ld large prime (power) factor%s.\n",
    3049                 :            :                nb, (nb>1? "s": ""));
    3050                 :       1131 :   return nb;
    3051                 :            : }
    3052                 :            : 
    3053                 :            : /***********************************************************************/
    3054                 :            : /**            ARITHMETIC FUNCTIONS WITH EARLY-ABORT                  **/
    3055                 :            : /**  needing direct access to the factoring machinery to avoid work:  **/
    3056                 :            : /**  e.g. if we find a square factor, moebius returns 0, core doesn't **/
    3057                 :            : /**  need to factor it, etc.                                          **/
    3058                 :            : /***********************************************************************/
    3059                 :            : /* memory management */
    3060                 :            : static void
    3061                 :          0 : ifac_GC(pari_sp av, GEN *part)
    3062                 :            : {
    3063                 :          0 :   GEN here = NULL;
    3064         [ #  # ]:          0 :   if(DEBUGMEM>1) pari_warn(warnmem,"ifac_xxx");
    3065                 :          0 :   ifac_realloc(part, &here, 0);
    3066                 :          0 :   *part = gerepileupto(av, *part);
    3067                 :          0 : }
    3068                 :            : 
    3069                 :            : static long
    3070                 :        189 : ifac_moebius(GEN n)
    3071                 :            : {
    3072                 :        189 :   long mu = 1;
    3073                 :        189 :   pari_sp av = avma, lim = stack_lim(av,1);
    3074                 :        189 :   GEN part = ifac_start(n, 1);
    3075                 :            :   for(;;)
    3076                 :            :   {
    3077                 :            :     long v;
    3078                 :            :     GEN p;
    3079 [ +  + ][ +  + ]:        539 :     if (!ifac_next(&part,&p,&v)) return v? 0: mu;
    3080                 :        350 :     mu = -mu;
    3081         [ -  + ]:        350 :     if (low_stack(lim, stack_lim(av,1))) ifac_GC(av,&part);
    3082                 :        350 :   }
    3083                 :            : }
    3084                 :            : 
    3085                 :            : int
    3086                 :         68 : ifac_read(GEN part, GEN *p, long *e)
    3087                 :            : {
    3088                 :         68 :   GEN here = ifac_find(part);
    3089         [ +  + ]:         68 :   if (!here) return 0;
    3090                 :         38 :   *p = VALUE(here);
    3091                 :         38 :   *e = EXPON(here)[2];
    3092                 :         68 :   return 1;
    3093                 :            : }
    3094                 :            : void
    3095                 :         22 : ifac_skip(GEN part)
    3096                 :            : {
    3097                 :         22 :   GEN here = ifac_find(part);
    3098         [ +  - ]:         22 :   if (here) ifac_delete(here);
    3099                 :         22 : }
    3100                 :            : 
    3101                 :            : static int
    3102                 :          7 : ifac_ispowerful(GEN n)
    3103                 :            : {
    3104                 :          7 :   pari_sp av = avma, lim = stack_lim(av,1);
    3105                 :          7 :   GEN part = ifac_start(n, 0);
    3106                 :            :   for(;;)
    3107                 :            :   {
    3108                 :            :     long e;
    3109                 :            :     GEN p;
    3110         [ +  + ]:         14 :     if (!ifac_read(part,&p,&e)) return 1;
    3111                 :            :     /* power: skip */
    3112 [ +  - ][ +  - ]:          7 :     if (e != 1 || Z_isanypower(p,NULL)) { ifac_skip(part); continue; }
    3113         [ #  # ]:          0 :     if (!ifac_next(&part,&p,&e)) return 1;
    3114         [ #  # ]:          0 :     if (e == 1) return 0;
    3115         [ #  # ]:          7 :     if (low_stack(lim, stack_lim(av,1))) ifac_GC(av,&part);
    3116                 :         14 :   }
    3117                 :            : }
    3118                 :            : static GEN
    3119                 :         23 : ifac_core(GEN n)
    3120                 :            : {
    3121                 :         23 :   GEN m = gen_1, c = cgeti(lgefint(n));
    3122                 :         23 :   pari_sp av = avma, lim = stack_lim(av,1);
    3123                 :         23 :   GEN part = ifac_start(n, 0);
    3124                 :            :   for(;;)
    3125                 :            :   {
    3126                 :            :     long e;
    3127                 :            :     GEN p;
    3128         [ +  + ]:         54 :     if (!ifac_read(part,&p,&e)) return m;
    3129                 :            :     /* square: skip */
    3130 [ +  - ][ +  + ]:         31 :     if (!odd(e) || Z_issquare(p)) { ifac_skip(part); continue; }
    3131         [ -  + ]:         16 :     if (!ifac_next(&part,&p,&e)) return m;
    3132         [ +  + ]:         16 :     if (odd(e)) m = mulii(m, p);
    3133         [ -  + ]:         39 :     if (low_stack(lim, stack_lim(av,1))) { affii(m,c); m=c; ifac_GC(av,&part); }
    3134                 :         54 :   }
    3135                 :            : }
    3136                 :            : 
    3137                 :            : /* Where to stop trial dividing in factorization. Guaranteed >= 2^14 */
    3138                 :            : ulong
    3139                 :      22830 : tridiv_bound(GEN n)
    3140                 :            : {
    3141                 :      22830 :   ulong l = (ulong)expi(n) + 1;
    3142         [ +  + ]:      22830 :   if (l <= 32)  return 1UL<<14;
    3143         [ +  + ]:      22732 :   if (l <= 512) return (l-16) << 10;
    3144                 :      22830 :   return 1UL<<19; /* Rho is generally faster above this */
    3145                 :            : }
    3146                 :            : 
    3147                 :            : /* return a value <= (48 << 10) = 49152 < primelinit */
    3148                 :            : static ulong
    3149                 :    1853743 : utridiv_bound(ulong n)
    3150                 :            : {
    3151                 :            : #ifdef LONG_IS_64BIT
    3152         [ +  + ]:    1600776 :   if (n & HIGHMASK)
    3153                 :      83976 :     return ((ulong)expu(n) + 1 - 16) << 10;
    3154                 :            : #else
    3155                 :            :   (void)n;
    3156                 :            : #endif
    3157                 :    1853743 :   return 1UL<<14;
    3158                 :            : }
    3159                 :            : 
    3160                 :            : static void
    3161                 :        866 : ifac_factoru(GEN n, long hint, GEN P, GEN E, long *pi)
    3162                 :            : {
    3163                 :        866 :   GEN part = ifac_start_hint(n, 0, hint);
    3164                 :            :   for(;;)
    3165                 :            :   {
    3166                 :            :     long v;
    3167                 :            :     GEN p;
    3168         [ +  + ]:       3410 :     if (!ifac_next(&part,&p,&v)) return;
    3169                 :       1678 :     P[*pi] = itou(p);
    3170                 :       1678 :     E[*pi] = v;
    3171                 :       1678 :     (*pi)++;
    3172                 :       2544 :   }
    3173                 :            : }
    3174                 :            : static long
    3175                 :       1068 : ifac_moebiusu(GEN n)
    3176                 :            : {
    3177                 :       1068 :   GEN part = ifac_start(n, 1);
    3178                 :       1068 :   long s = 1;
    3179                 :            :   for(;;)
    3180                 :            :   {
    3181                 :            :     long v;
    3182                 :            :     GEN p;
    3183 [ +  + ][ +  - ]:       3204 :     if (!ifac_next(&part,&p,&v)) return v? 0: s;
    3184                 :       2136 :     s = -s;
    3185                 :       2136 :   }
    3186                 :            : }
    3187                 :            : 
    3188                 :            : INLINE ulong
    3189                 :  340816680 : u_forprime_next_fast(forprime_t *T)
    3190                 :            : {
    3191         [ +  + ]:  340816680 :   if (*(T->d))
    3192                 :            :   {
    3193                 :  340799445 :     NEXT_PRIME_VIADIFF(T->p, T->d);
    3194         [ +  + ]:  340799445 :     return T->p > T->b ? 0: T->p;
    3195                 :            :   }
    3196                 :  340816680 :   return u_forprime_next(T);
    3197                 :            : }
    3198                 :            : 
    3199                 :            : /* Factor n and output [p,e] where
    3200                 :            :  * p, e are vecsmall with n = prod{p[i]^e[i]} */
    3201                 :            : static GEN
    3202                 :    2051261 : factoru_sign(ulong n, ulong all, long hint)
    3203                 :            : {
    3204                 :            :   GEN f, E, E2, P, P2;
    3205                 :            :   pari_sp av;
    3206                 :            :   ulong p, lim;
    3207                 :            :   long i;
    3208                 :            :   forprime_t S;
    3209                 :            : 
    3210         [ -  + ]:    2051261 :   if (n == 0) retmkvec2(mkvecsmall(0), mkvecsmall(1));
    3211         [ +  + ]:    2051261 :   if (n == 1) retmkvec2(cgetg(1,t_VECSMALL), cgetg(1,t_VECSMALL));
    3212                 :            : 
    3213                 :    1849507 :   f = cgetg(3,t_VEC); av = avma;
    3214         [ +  + ]:    1849507 :   lim = all; if (!lim) lim = utridiv_bound(n);
    3215                 :            :   /* enough room to store <= 15 primes and exponents (OK if n < 2^64) */
    3216                 :    1849507 :   (void)new_chunk(16*2);
    3217                 :    1849507 :   P = cgetg(16, t_VECSMALL); i = 1;
    3218                 :    1849507 :   E = cgetg(16, t_VECSMALL);
    3219         [ +  - ]:    1849507 :   if (lim > 2)
    3220                 :            :   {
    3221                 :    1849507 :     long v = vals(n), oldi;
    3222         [ +  + ]:    1849507 :     if (v)
    3223                 :            :     {
    3224                 :    1080480 :       P[1] = 2; E[1] = v; i = 2;
    3225         [ +  + ]:    1080480 :       n >>= v; if (n == 1) goto END;
    3226                 :            :     }
    3227                 :    1428688 :     u_forprime_init(&S, 3, lim);
    3228                 :    1428688 :     oldi = i;
    3229         [ +  + ]:   14701243 :     while ( (p = u_forprime_next_fast(&S)) )
    3230                 :            :     {
    3231                 :            :       int stop;
    3232                 :            :       /* tiny integers without small factors are often primes */
    3233         [ +  + ]:   14700167 :       if (p == 673)
    3234                 :            :       {
    3235                 :      10620 :         oldi = i;
    3236         [ +  + ]:      10620 :         if (uisprime_661(n)) { P[i] = n; E[i] = 1; i++; goto END; }
    3237                 :            :       }
    3238                 :   14695893 :       v = u_lvalrem_stop(&n, p, &stop);
    3239         [ +  + ]:   14695893 :       if (v) {
    3240                 :    1612093 :         P[i] = p;
    3241                 :    1612093 :         E[i] = v; i++;
    3242                 :            :       }
    3243         [ +  + ]:   14695893 :       if (stop) {
    3244         [ +  + ]:   14700167 :         if (n != 1) { P[i] = n; E[i] = 1; i++; }
    3245                 :            :         goto END;
    3246                 :            :       }
    3247                 :            :     }
    3248 [ +  + ][ +  + ]:       1076 :     if (oldi != i && uisprime_661(n)) { P[i] = n; E[i] = 1; i++; goto END; }
    3249                 :            :   }
    3250         [ +  + ]:        878 :   if (all)
    3251                 :            :   { /* smallfact: look for easy pure powers then stop */
    3252                 :            : #ifdef LONG_IS_64BIT
    3253 [ +  - ][ +  - ]:         12 :     ulong mask = all > 563 ? (all > 7129 ? 1: 3): 7;
    3254                 :            : #else
    3255 [ #  # ][ #  # ]:          0 :     ulong mask = all > 22 ? (all > 83 ? 1: 3): 7;
    3256                 :            : #endif
    3257                 :         12 :     long k = 1, ex;
    3258         [ +  + ]:         18 :     while (uissquareall(n, &n)) k <<= 1;
    3259         [ -  + ]:         12 :     while ( (ex = uis_357_power(n, &n, &mask)) ) k *= ex;
    3260                 :         12 :     P[i] = n; E[i] = k; i++; goto END;
    3261                 :            :   }
    3262                 :            :   {
    3263                 :            :     GEN perm;
    3264                 :        866 :     ifac_factoru(utoipos(n), hint, P, E, &i);
    3265                 :        866 :     setlg(P, i);
    3266                 :        866 :     perm = vecsmall_indexsort(P);
    3267                 :        866 :     P = vecsmallpermute(P, perm);
    3268                 :        866 :     E = vecsmallpermute(E, perm);
    3269                 :            :   }
    3270                 :            : END:
    3271                 :    1849507 :   avma = av;
    3272                 :    1849507 :   P2 = cgetg(i, t_VECSMALL); gel(f,1) = P2;
    3273                 :    1849507 :   E2 = cgetg(i, t_VECSMALL); gel(f,2) = E2;
    3274         [ +  + ]:    5401606 :   while (--i >= 1) { P2[i] = P[i]; E2[i] = E[i]; }
    3275                 :    2051261 :   return f;
    3276                 :            : }
    3277                 :            : GEN
    3278                 :     421994 : factoru(ulong n)
    3279                 :     421994 : { return factoru_sign(n, 0, decomp_default_hint); }
    3280                 :            : 
    3281                 :            : long
    3282                 :      11655 : moebiusu(ulong n)
    3283                 :            : {
    3284                 :            :   pari_sp av;
    3285                 :            :   ulong p;
    3286                 :            :   long s, v, test_prime;
    3287                 :            :   forprime_t S;
    3288                 :            : 
    3289   [ -  +  +  + ]:      11655 :   switch(n)
    3290                 :            :   {
    3291                 :          0 :     case 0: (void)check_arith_non0(gen_0,"moebius");/*error*/
    3292                 :       1288 :     case 1: return  1;
    3293                 :        917 :     case 2: return -1;
    3294                 :            :   }
    3295                 :       9450 :   v = vals(n);
    3296         [ +  + ]:       9450 :   if (v == 0)
    3297                 :       7816 :     s = 1;
    3298                 :            :   else
    3299                 :            :   {
    3300         [ +  + ]:       1634 :     if (v > 1) return 0;
    3301                 :       1256 :     n >>= 1;
    3302                 :       1256 :     s = -1;
    3303                 :            :   }
    3304                 :       9072 :   av = avma;
    3305                 :       9072 :   u_forprime_init(&S, 3, utridiv_bound(n));
    3306                 :       9072 :   test_prime = 0;
    3307         [ +  + ]:    6922627 :   while ((p = u_forprime_next_fast(&S)))
    3308                 :            :   {
    3309                 :            :     int stop;
    3310                 :            :     /* tiny integers without small factors are often primes */
    3311         [ +  + ]:    6921115 :     if (p == 673)
    3312                 :            :     {
    3313                 :       3728 :       test_prime = 0;
    3314         [ +  + ]:       3728 :       if (uisprime_661(n)) { avma = av; return -s; }
    3315                 :            :     }
    3316                 :    6919679 :     v = u_lvalrem_stop(&n, p, &stop);
    3317         [ +  + ]:    6919679 :     if (v) {
    3318         [ +  + ]:      10908 :       if (v > 1) { avma = av; return 0; }
    3319                 :       9792 :       test_prime = 1;
    3320                 :       9792 :       s = -s;
    3321                 :            :     }
    3322 [ +  + ][ +  + ]:    6921115 :     if (stop) { avma = av; return n == 1? s: -s; }
    3323                 :            :   }
    3324                 :       1512 :   avma = av;
    3325 [ +  + ][ +  + ]:       1512 :   if (test_prime && uisprime_661(n)) return -s;
    3326                 :            :   else
    3327                 :            :   {
    3328                 :       1068 :     long t = ifac_moebiusu(utoipos(n));
    3329                 :       1068 :     avma = av;
    3330         [ -  + ]:       1068 :     if (t == 0) return 0;
    3331         [ +  + ]:      11655 :     return (s == t)? 1: -1;
    3332                 :            :   }
    3333                 :            : }
    3334                 :            : 
    3335                 :            : long
    3336                 :       1018 : moebius(GEN n)
    3337                 :            : {
    3338                 :       1018 :   pari_sp av = avma;
    3339                 :            :   GEN F;
    3340                 :            :   ulong p;
    3341                 :            :   long i, l, s, v;
    3342                 :            :   forprime_t S;
    3343                 :            : 
    3344         [ +  + ]:       1018 :   if ((F = check_arith_non0(n,"moebius")))
    3345                 :            :   {
    3346                 :            :     GEN E;
    3347                 :         21 :     F = clean_Z_factor(F);
    3348                 :         21 :     E = gel(F,2);
    3349                 :         21 :     l = lg(E);
    3350         [ +  + ]:         28 :     for(i = 1; i < l; i++)
    3351         [ +  + ]:         14 :       if (!equali1(gel(E,1))) { avma = av; return 0; }
    3352         [ +  + ]:         14 :     avma = av; return odd(l)? 1: -1;
    3353                 :            :   }
    3354         [ +  + ]:        983 :   if (lgefint(n) == 3) return moebiusu(uel(n,2));
    3355         [ -  + ]:        784 :   p = mod4(n); if (!p) return 0;
    3356         [ +  + ]:        784 :   if (p == 2) { s = -1; n = shifti(n, -1); } else { s = 1; n = icopy(n); }
    3357                 :        784 :   setabssign(n);
    3358                 :            : 
    3359                 :        784 :   u_forprime_init(&S, 3, tridiv_bound(n));
    3360         [ +  + ]:    2194014 :   while ((p = u_forprime_next_fast(&S)))
    3361                 :            :   {
    3362                 :            :     int stop;
    3363                 :    2193459 :     v = Z_lvalrem_stop(&n, p, &stop);
    3364         [ +  + ]:    2193459 :     if (v)
    3365                 :            :     {
    3366         [ +  + ]:       1426 :       if (v > 1) { avma = av; return 0; }
    3367                 :       1282 :       s = -s;
    3368 [ +  + ][ +  - ]:    2193459 :       if (stop) { avma = av; return is_pm1(n)? s: -s; }
    3369                 :            :     }
    3370                 :            :   }
    3371                 :        555 :   l = lg(primetab);
    3372         [ +  + ]:        565 :   for (i = 1; i < l; i++)
    3373                 :            :   {
    3374                 :         25 :     v = Z_pvalrem(n, gel(primetab,i), &n);
    3375         [ +  - ]:         25 :     if (v)
    3376                 :            :     {
    3377         [ +  + ]:         25 :       if (v > 1) { avma = av; return 0; }
    3378                 :         11 :       s = -s;
    3379         [ +  + ]:         11 :       if (is_pm1(n)) { avma = av; return s; }
    3380                 :            :     }
    3381                 :            :   }
    3382         [ +  + ]:        540 :   if (ifac_isprime(n)) { avma = av; return -s; }
    3383                 :            :   /* large composite without small factors */
    3384                 :        189 :   v = ifac_moebius(n);
    3385         [ +  + ]:       1004 :   avma = av; return (s<0 ? -v : v); /* correct also if v==0 */
    3386                 :            : }
    3387                 :            : 
    3388                 :            : long
    3389                 :        294 : ispowerful(GEN n)
    3390                 :            : {
    3391                 :        294 :   pari_sp av = avma;
    3392                 :            :   GEN F;
    3393                 :            :   ulong p, bound;
    3394                 :            :   long i, l, v;
    3395                 :            :   forprime_t S;
    3396                 :            : 
    3397         [ +  + ]:        294 :   if ((F = check_arith_all(n, "ispowerful")))
    3398                 :            :   {
    3399                 :         35 :     GEN p, P = gel(F,1), E = gel(F,2);
    3400         [ +  + ]:         35 :     if (lg(P) == 1) return 1; /* 1 */
    3401                 :         28 :     p = gel(P,1);
    3402         [ +  + ]:         28 :     if (!signe(p)) return 1; /* 0 */
    3403                 :         14 :     l = lg(E);
    3404         [ +  - ]:         35 :     for (i = 1; i < l; i++)
    3405         [ +  + ]:         35 :       if (equali1(gel(E,i))) return 0;
    3406                 :          0 :     return 1;
    3407                 :            :   }
    3408         [ +  + ]:        259 :   if (!signe(n)) return 1;
    3409                 :            : 
    3410         [ +  + ]:        252 :   if (mod4(n) == 2) return 0;
    3411                 :        105 :   n = shifti(n, -vali(n));
    3412         [ +  + ]:        105 :   if (is_pm1(n)) return 1;
    3413                 :         98 :   setabssign(n);
    3414                 :         98 :   bound = tridiv_bound(n);
    3415                 :         98 :   u_forprime_init(&S, 3, bound);
    3416         [ +  + ]:     305921 :   while ((p = u_forprime_next_fast(&S)))
    3417                 :            :   {
    3418                 :            :     int stop;
    3419                 :     305914 :     v = Z_lvalrem_stop(&n, p, &stop);
    3420         [ +  + ]:     305914 :     if (v)
    3421                 :            :     {
    3422         [ +  + ]:        126 :       if (v == 1) { avma = av; return 0; }
    3423         [ +  + ]:     305914 :       if (stop) { avma = av; return is_pm1(n); }
    3424                 :            :     }
    3425                 :            :   }
    3426                 :          7 :   l = lg(primetab);
    3427         [ -  + ]:          7 :   for (i = 1; i < l; i++)
    3428                 :            :   {
    3429                 :          0 :     v = Z_pvalrem(n, gel(primetab,i), &n);
    3430         [ #  # ]:          0 :     if (v)
    3431                 :            :     {
    3432         [ #  # ]:          0 :       if (v == 1) { avma = av; return 0; }
    3433         [ #  # ]:          0 :       if (is_pm1(n)) { avma = av; return 1; }
    3434                 :            :     }
    3435                 :            :   }
    3436                 :            :   /* no need to factor: must be p^2 or not powerful */
    3437         [ -  + ]:          7 :   if(cmpii(powuu(bound+1, 3), n) > 0) {
    3438                 :          0 :     long res = Z_issquare(n);
    3439                 :          0 :     avma = av; return res;
    3440                 :            :   }
    3441                 :            : 
    3442         [ -  + ]:          7 :   if (ifac_isprime(n)) { avma=av; return 0; }
    3443                 :            :   /* large composite without small factors */
    3444                 :          7 :   v = ifac_ispowerful(n);
    3445                 :        294 :   avma = av; return v;
    3446                 :            : }
    3447                 :            : 
    3448                 :            : ulong
    3449                 :        176 : coreu(ulong n)
    3450                 :            : {
    3451         [ -  + ]:        176 :   if (n == 0) return 0;
    3452                 :            :   else
    3453                 :            :   {
    3454                 :        176 :     pari_sp av = avma;
    3455                 :        176 :     GEN f = factoru(n), P = gel(f,1), E = gel(f,2);
    3456                 :        176 :     long i, l = lg(P), m = 1;
    3457                 :            : 
    3458                 :        176 :     avma = av;
    3459         [ +  + ]:        846 :     for (i = 1; i < l; i++)
    3460                 :            :     {
    3461                 :        670 :       ulong p = P[i], e = E[i];
    3462         [ +  + ]:        670 :       if (e & 1) m *= p;
    3463                 :            :     }
    3464                 :        176 :     return m;
    3465                 :            :   }
    3466                 :            : }
    3467                 :            : GEN
    3468                 :        308 : core(GEN n)
    3469                 :            : {
    3470                 :        308 :   pari_sp av = avma;
    3471                 :            :   GEN m, F;
    3472                 :            :   ulong p;
    3473                 :            :   long i, l, v;
    3474                 :            :   forprime_t S;
    3475                 :            : 
    3476         [ +  + ]:        308 :   if ((F = check_arith_all(n, "core")))
    3477                 :            :   {
    3478                 :         77 :     GEN p, x, P = gel(F,1), E = gel(F,2);
    3479                 :         77 :     long j = 1;
    3480         [ +  + ]:         77 :     if (lg(P) == 1) return gen_1;
    3481                 :         63 :     p = gel(P,1);
    3482         [ +  + ]:         63 :     if (!signe(p)) return gen_0;
    3483                 :         35 :     l = lg(P); x = cgetg(l, t_VEC);
    3484         [ +  + ]:        140 :     for (i = 1; i < l; i++)
    3485         [ +  + ]:        105 :       if (mpodd(gel(E,i))) gel(x,j++) = gel(P,i);
    3486                 :         35 :     setlg(x, j); return ZV_prod(x);
    3487                 :            :   }
    3488      [ +  +  + ]:        231 :   switch(lgefint(n))
    3489                 :            :   {
    3490                 :         14 :     case 2: return gen_0;
    3491                 :            :     case 3:
    3492                 :        176 :       p = coreu(uel(n,2));
    3493         [ +  + ]:        176 :       return signe(n) > 0? utoipos(p): utoineg(p);
    3494                 :            :   }
    3495                 :            : 
    3496         [ -  + ]:         41 :   m = signe(n) < 0? gen_m1: gen_1;
    3497                 :         41 :   n = absi(n);
    3498                 :         41 :   u_forprime_init(&S, 2, tridiv_bound(n));
    3499         [ +  + ]:     195447 :   while ((p = u_forprime_next_fast(&S)))
    3500                 :            :   {
    3501                 :            :     int stop;
    3502                 :     195407 :     v = Z_lvalrem_stop(&n, p, &stop);
    3503         [ +  + ]:     195407 :     if (v)
    3504                 :            :     {
    3505         [ +  + ]:        139 :       if (v & 1) m = muliu(m, p);
    3506         [ +  + ]:        139 :       if (stop)
    3507                 :            :       {
    3508         [ +  - ]:          1 :         if (!is_pm1(n)) m = mulii(m, n);
    3509                 :     195407 :         return gerepileuptoint(av, m);
    3510                 :            :       }
    3511                 :            :     }
    3512                 :            :   }
    3513                 :         40 :   l = lg(primetab);
    3514         [ +  + ]:         80 :   for (i = 1; i < l; i++)
    3515                 :            :   {
    3516                 :         48 :     GEN q = gel(primetab,i);
    3517                 :         48 :     v = Z_pvalrem(n, q, &n);
    3518         [ +  + ]:         48 :     if (v)
    3519                 :            :     {
    3520         [ +  + ]:         32 :       if (v & 1) m = mulii(m, q);
    3521         [ +  + ]:         32 :       if (is_pm1(n)) return gerepileuptoint(av, m);
    3522                 :            :     }
    3523                 :            :   }
    3524         [ +  + ]:         32 :   if (ifac_isprime(n)) { m = mulii(m, n); return gerepileuptoint(av, m); }
    3525                 :            :   /* large composite without small factors */
    3526                 :        308 :   return gerepileuptoint(av, mulii(m, ifac_core(n)));
    3527                 :            : }
    3528                 :            : 
    3529                 :            : long
    3530                 :          0 : Z_issmooth(GEN m, ulong lim)
    3531                 :            : {
    3532                 :          0 :   pari_sp av=avma;
    3533                 :          0 :   ulong p = 2;
    3534                 :            :   forprime_t S;
    3535                 :          0 :   u_forprime_init(&S, 2, lim);
    3536         [ #  # ]:          0 :   while ((p = u_forprime_next_fast(&S)))
    3537                 :            :   {
    3538                 :            :     int stop;
    3539                 :          0 :     (void)Z_lvalrem_stop(&m, p, &stop);
    3540         [ #  # ]:          0 :     if (stop) { avma = av; return cmpiu(m,lim)<=0; }
    3541                 :            :   }
    3542                 :          0 :   avma = av; return 0;
    3543                 :            : }
    3544                 :            : 
    3545                 :            : GEN
    3546                 :     270599 : Z_issmooth_fact(GEN m, ulong lim)
    3547                 :            : {
    3548                 :     270599 :   pari_sp av=avma;
    3549                 :            :   GEN F, P, E;
    3550                 :            :   ulong p;
    3551                 :     270599 :   long i = 1, l = expi(m)+1;
    3552                 :            :   forprime_t S;
    3553                 :     270599 :   P = cgetg(l, t_VECSMALL);
    3554                 :     270599 :   E = cgetg(l, t_VECSMALL);
    3555                 :     270599 :   F = mkmat2(P,E);
    3556                 :     270599 :   u_forprime_init(&S, 2, lim);
    3557         [ +  + ]:   86956254 :   while ((p = u_forprime_next_fast(&S)))
    3558                 :            :   {
    3559                 :            :     long v;
    3560                 :            :     int stop;
    3561         [ +  + ]:   86903880 :     if ((v = Z_lvalrem_stop(&m, p, &stop)))
    3562                 :            :     {
    3563                 :    1154993 :       P[i] = p;
    3564                 :    1154993 :       E[i] = v; i++;
    3565         [ +  + ]:    1154993 :       if (stop)
    3566                 :            :       {
    3567         [ +  + ]:     218225 :         if (cmpiu(m,lim) > 0) break;
    3568                 :      76734 :         P[i] = m[2];
    3569                 :      76734 :         E[i] = 1; i++;
    3570                 :      76734 :         setlg(P, i);
    3571                 :   86903880 :         setlg(E, i); avma = (pari_sp)F; return F;
    3572                 :            :       }
    3573                 :            :     }
    3574                 :            :   }
    3575                 :     270599 :   avma = av; return NULL;
    3576                 :            : }
    3577                 :            : 
    3578                 :            : /***********************************************************************/
    3579                 :            : /**                                                                   **/
    3580                 :            : /**       COMPUTING THE MATRIX OF PRIME DIVISORS AND EXPONENTS        **/
    3581                 :            : /**                                                                   **/
    3582                 :            : /***********************************************************************/
    3583                 :            : static GEN
    3584                 :      30567 : aux_end(GEN M, GEN n, long nb)
    3585                 :            : {
    3586                 :      30567 :   GEN P,E, z = (GEN)avma;
    3587                 :            :   long i;
    3588                 :            : 
    3589         [ +  - ]:      30567 :   if (n) gunclone(n);
    3590                 :      30567 :   P = cgetg(nb+1,t_COL);
    3591                 :      30567 :   E = cgetg(nb+1,t_COL);
    3592         [ +  + ]:     182262 :   for (i=nb; i; i--)
    3593                 :            :   { /* allow a stackdummy in the middle */
    3594         [ -  + ]:     151695 :     while (typ(z) != t_INT) z += lg(z);
    3595                 :     151695 :     gel(E,i) = z; z += lg(z);
    3596                 :     151695 :     gel(P,i) = z; z += lg(z);
    3597                 :            :   }
    3598                 :      30567 :   gel(M,1) = P;
    3599                 :      30567 :   gel(M,2) = E;
    3600                 :      30567 :   return sort_factor(M, (void*)&absi_cmp, cmp_nodata);
    3601                 :            : }
    3602                 :            : 
    3603                 :            : static void
    3604                 :     149371 : STORE(long *nb, GEN x, long e) { (*nb)++; (void)x; (void)utoipos(e); }
    3605                 :            : static void
    3606                 :     132869 : STOREu(long *nb, ulong x, long e) { STORE(nb, utoipos(x), e); }
    3607                 :            : static void
    3608                 :      16481 : STOREi(long *nb, GEN x, long e) { STORE(nb, icopy(x), e); }
    3609                 :            : /* no prime less than p divides n */
    3610                 :            : static int
    3611                 :       9855 : special_primes(GEN n, ulong p, long *nb, GEN T)
    3612                 :            : {
    3613                 :       9855 :   long i, l = lg(T);
    3614         [ +  + ]:       9855 :   if (l > 1)
    3615                 :            :   { /* pp = square of biggest p tried so far */
    3616                 :        207 :     long pp[] = { evaltyp(t_INT)|_evallg(4), 0,0,0 };
    3617                 :        207 :     pari_sp av = avma; affii(sqru(p), pp); avma = av;
    3618                 :            : 
    3619         [ +  + ]:        334 :     for (i = 1; i < l; i++)
    3620         [ +  + ]:        246 :       if (dvdiiz(n,gel(T,i), n))
    3621                 :            :       {
    3622         [ +  + ]:        266 :         long k = 1; while (dvdiiz(n,gel(T,i), n)) k++;
    3623                 :        168 :         STOREi(nb, gel(T,i), k);
    3624         [ +  + ]:        168 :         if (absi_cmp(pp, n) > 0) return 1;
    3625                 :            :       }
    3626                 :            :   }
    3627                 :       9855 :   return 0;
    3628                 :            : }
    3629                 :            : 
    3630                 :            : /* factor(sn*|n|), where sn = -1,1 or 0.
    3631                 :            :  * all != 0 : only look for prime divisors < all */
    3632                 :            : static GEN
    3633                 :    1659841 : ifactor_sign(GEN n, ulong all, long hint, long sn)
    3634                 :            : {
    3635                 :            :   GEN M, N;
    3636                 :            :   pari_sp av;
    3637                 :    1659841 :   long nb = 0, i;
    3638                 :            :   ulong lim;
    3639                 :            :   forprime_t T;
    3640                 :            : 
    3641         [ +  + ]:    1659841 :   if (!sn) retmkmat2(mkcol(gen_0), mkcol(gen_1));
    3642         [ +  + ]:    1659834 :   if (lgefint(n) == 3)
    3643                 :            :   { /* small integer */
    3644                 :    1629267 :     GEN f, Pf, Ef, P, E, F = cgetg(3, t_MAT);
    3645                 :            :     long l;
    3646                 :    1629267 :     av = avma;
    3647                 :            :     /* enough room to store <= 15 primes and exponents (OK if n < 2^64) */
    3648                 :    1629267 :     (void)new_chunk((15*3 + 15 + 1) * 2);
    3649                 :    1629267 :     f = factoru_sign(uel(n,2), all, hint);
    3650                 :    1629267 :     avma = av;
    3651                 :    1629267 :     Pf = gel(f,1);
    3652                 :    1629267 :     Ef = gel(f,2);
    3653                 :    1629267 :     l = lg(Pf);
    3654         [ +  + ]:    1629267 :     if (sn < 0)
    3655                 :            :     { /* add sign */
    3656                 :         35 :       long L = l+1;
    3657                 :         35 :       gel(F,1) = P = cgetg(L, t_COL);
    3658                 :         35 :       gel(F,2) = E = cgetg(L, t_COL);
    3659                 :         35 :       gel(P,1) = gen_m1; P++;
    3660                 :         35 :       gel(E,1) = gen_1;  E++;
    3661                 :            :     }
    3662                 :            :     else
    3663                 :            :     {
    3664                 :    1629232 :       gel(F,1) = P = cgetg(l, t_COL);
    3665                 :    1629232 :       gel(F,2) = E = cgetg(l, t_COL);
    3666                 :            :     }
    3667         [ +  + ]:    4378326 :     for (i = 1; i < l; i++)
    3668                 :            :     {
    3669                 :    2749059 :       gel(P,i) = utoipos(Pf[i]);
    3670                 :    2749059 :       gel(E,i) = utoipos(Ef[i]);
    3671                 :            :     }
    3672                 :    1629267 :     return F;
    3673                 :            :   }
    3674                 :      30567 :   M = cgetg(3,t_MAT);
    3675         [ +  + ]:      30567 :   if (sn < 0) STORE(&nb, utoineg(1), 1);
    3676         [ -  + ]:      30567 :   if (is_pm1(n)) return aux_end(M,NULL,nb);
    3677                 :            : 
    3678                 :      30567 :   n = N = gclone(n); setabssign(n);
    3679                 :            :   /* trial division bound */
    3680         [ +  + ]:      30567 :   lim = all; if (!lim) lim = tridiv_bound(n);
    3681         [ +  + ]:      30567 :   if (lim > 2)
    3682                 :            :   {
    3683                 :            :     ulong maxp, p;
    3684                 :            :     pari_sp av2;
    3685                 :      30532 :     i = vali(n);
    3686         [ +  + ]:      30532 :     if (i)
    3687                 :            :     {
    3688                 :      22638 :       STOREu(&nb, 2, i);
    3689                 :      22638 :       av = avma; affii(shifti(n,-i), n); avma = av;
    3690                 :            :     }
    3691         [ +  + ]:      30532 :     if (is_pm1(n)) return aux_end(M,n,nb);
    3692                 :            :     /* trial division */
    3693                 :      30452 :     maxp = maxprime();
    3694                 :      30452 :     av = avma; u_forprime_init(&T, 3, minss(lim, maxp)); av2 = avma;
    3695                 :            :     /* first pass: known to fit in private prime table */
    3696         [ +  + ]:  229541174 :     while ((p = u_forprime_next_fast(&T)))
    3697                 :            :     {
    3698                 :  229531347 :       pari_sp av3 = avma;
    3699                 :            :       int stop;
    3700                 :  229531347 :       long k = Z_lvalrem_stop(&n, p, &stop);
    3701         [ +  + ]:  229531347 :       if (k)
    3702                 :            :       {
    3703                 :     110182 :         affii(n, N); n = N; avma = av3;
    3704                 :     110182 :         STOREu(&nb, p, k);
    3705                 :            :       }
    3706         [ +  + ]:  229531347 :       if (stop)
    3707                 :            :       {
    3708         [ +  + ]:      20625 :         if (!is_pm1(n)) STOREi(&nb, n, 1);
    3709                 :      20625 :         stackdummy(av, av2);
    3710                 :  229531347 :         return aux_end(M,n,nb);
    3711                 :            :       }
    3712                 :            :     }
    3713                 :       9827 :     stackdummy(av, av2);
    3714         [ +  + ]:       9827 :     if (lim > maxp)
    3715                 :            :     { /* second pass, usually empty: outside private prime table */
    3716                 :        737 :       av = avma; u_forprime_init(&T, maxp+1, lim); av2 = avma;
    3717         [ +  + ]:     130748 :       while ((p = u_forprime_next(&T)))
    3718                 :            :       {
    3719                 :     130018 :         pari_sp av3 = avma;
    3720                 :            :         int stop;
    3721                 :     130018 :         long k = Z_lvalrem_stop(&n, p, &stop);
    3722         [ +  + ]:     130018 :         if (k)
    3723                 :            :         {
    3724                 :         49 :           affii(n, N); n = N; avma = av3;
    3725                 :         49 :           STOREu(&nb, p, k);
    3726                 :            :         }
    3727         [ +  + ]:     130018 :         if (stop)
    3728                 :            :         {
    3729         [ -  + ]:          7 :           if (!is_pm1(n)) STOREi(&nb, n, 1);
    3730                 :          7 :           stackdummy(av, av2);
    3731                 :     130018 :           return aux_end(M,n,nb);
    3732                 :            :         }
    3733                 :            :       }
    3734                 :        730 :       stackdummy(av, av2);
    3735                 :            :     }
    3736                 :            :   }
    3737                 :            :   /* trial divide by the special primes */
    3738         [ +  + ]:       9855 :   if (special_primes(n, lim, &nb, primetab))
    3739                 :            :   {
    3740         [ +  - ]:        119 :     if (!is_pm1(n)) STOREi(&nb, n, 1);
    3741                 :        119 :     return aux_end(M,n,nb);
    3742                 :            :   }
    3743                 :            : 
    3744         [ +  + ]:       9736 :   if (all)
    3745                 :            :   { /* smallfact: look for easy pure powers then stop. Cf Z_isanypower */
    3746                 :            :     GEN x;
    3747                 :            :     long k;
    3748                 :       7592 :     av = avma;
    3749                 :       7592 :     k = isanypower_nosmalldiv(n, &x);
    3750         [ +  + ]:       7592 :     if (k > 1) affii(x, n);
    3751                 :       7592 :     avma = av; STOREi(&nb, n, k);
    3752         [ -  + ]:       7592 :     if (DEBUGLEVEL >= 2) {
    3753                 :          0 :       pari_warn(warner,
    3754                 :            :         "IFAC: untested %ld-bit integer declared prime", expi(n));
    3755         [ #  # ]:          0 :       if (expi(n) <= 256)
    3756                 :          0 :         err_printf("\t%Ps\n", n);
    3757                 :            :     }
    3758                 :       7592 :     return aux_end(M,n,nb);
    3759                 :            :   }
    3760         [ +  + ]:       2144 :   if (ifac_isprime(n)) { STOREi(&nb, n, 1); return aux_end(M,n,nb); }
    3761                 :       1131 :   nb += ifac_decomp(n, hint);
    3762                 :    1659841 :   return aux_end(M,n, nb);
    3763                 :            : }
    3764                 :            : 
    3765                 :            : static GEN
    3766                 :    1179863 : ifactor(GEN n, ulong all, long hint)
    3767                 :    1179863 : { return ifactor_sign(n, all, hint, signe(n)); }
    3768                 :            : 
    3769                 :            : int
    3770                 :       7108 : ifac_next(GEN *part, GEN *p, long *e)
    3771                 :            : {
    3772                 :       7108 :   GEN here = ifac_main(part);
    3773         [ +  + ]:       7108 :   if (here == gen_0) { *p = NULL; *e = 1; return 0; }
    3774         [ +  + ]:       7094 :   if (!here) { *p = NULL; *e = 0; return 0; }
    3775                 :       4901 :   *p = VALUE(here);
    3776                 :       4901 :   *e = EXPON(here)[2];
    3777                 :       7108 :   ifac_delete(here); return 1;
    3778                 :            : }
    3779                 :            : 
    3780                 :            : /* see before ifac_crack for current semantics of 'hint' (factorint's 'flag') */
    3781                 :            : GEN
    3782                 :         14 : factorint(GEN n, long flag)
    3783                 :            : {
    3784                 :            :   GEN F;
    3785         [ -  + ]:         14 :   if ((F = check_arith_all(n,"factorint"))) return gcopy(F);
    3786                 :         14 :   return ifactor(n,0,flag);
    3787                 :            : }
    3788                 :            : 
    3789                 :            : GEN
    3790                 :       4893 : Z_factor_limit(GEN n, ulong all)
    3791                 :            : {
    3792         [ +  + ]:       4893 :   if (!all) all = GP_DATA->primelimit + 1;
    3793                 :       4893 :   return ifactor(n,all,decomp_default_hint);
    3794                 :            : }
    3795                 :            : GEN
    3796                 :       8974 : absi_factor_limit(GEN n, ulong all)
    3797                 :            : {
    3798         [ +  + ]:       8974 :   if (!all) all = GP_DATA->primelimit + 1;
    3799                 :       8974 :   return ifactor_sign(n,all,decomp_default_hint, signe(n)?1 : 0);
    3800                 :            : }
    3801                 :            : GEN
    3802                 :    1174277 : Z_factor(GEN n)
    3803                 :    1174277 : { return ifactor(n,0,decomp_default_hint); }
    3804                 :            : GEN
    3805                 :     471004 : absi_factor(GEN n)
    3806                 :     471004 : { return ifactor_sign(n, 0, decomp_default_hint, signe(n)? 1: 0); }
    3807                 :            : 
    3808                 :            : /* Factor until the unfactored part is smaller than limit. Return the
    3809                 :            :  * factored part. Hence factorback(output) may be smaller than n */
    3810                 :            : GEN
    3811                 :        679 : Z_factor_until(GEN n, GEN limit)
    3812                 :            : {
    3813                 :        679 :   pari_sp av2, av = avma;
    3814                 :        679 :   ulong B = tridiv_bound(n);
    3815                 :        679 :   GEN q, part, F = ifactor(n, B, decomp_default_hint);
    3816                 :        679 :   GEN P = gel(F,1), E = gel(F,2), P2, E2, F2;
    3817                 :        679 :   long l = lg(P), l2;
    3818                 :            : 
    3819                 :        679 :   av2 = avma;
    3820                 :        679 :   q = gel(P,l-1);
    3821 [ +  + ][ +  - ]:        679 :   if (cmpiu(q, B) <= 0 || cmpii(q, sqru(B)) < 0 || ifac_isprime(q))
                 [ +  + ]
    3822                 :            :   {
    3823                 :        210 :     avma = av2; return F;
    3824                 :            :   }
    3825                 :            :   /* q = composite unfactored part, remove from P/E */
    3826                 :        469 :   setlg(E,l-1);
    3827                 :        469 :   setlg(P,l-1);
    3828         [ +  + ]:        469 :   if (cmpii(q, limit) > 0)
    3829                 :            :   { /* factor further */
    3830                 :        441 :     l2 = expi(q)+1;
    3831                 :        441 :     P2 = vectrunc_init(l2);
    3832                 :        441 :     E2 = vectrunc_init(l2);
    3833                 :        441 :     F2 = mkmat2(P2,E2);
    3834                 :        441 :     part = ifac_start(q, 0);
    3835                 :            :     for(;;)
    3836                 :            :     {
    3837                 :            :       long e;
    3838                 :            :       GEN p;
    3839         [ +  - ]:        581 :       if (!ifac_next(&part,&p,&e)) break;
    3840                 :        581 :       vectrunc_append(P2, p);
    3841                 :        581 :       vectrunc_append(E2, utoipos(e));
    3842                 :        581 :       q = diviiexact(q, powiu(p, e));
    3843         [ +  + ]:        581 :       if (cmpii(q, limit) <= 0) break;
    3844                 :        140 :     }
    3845                 :        441 :     F2 = sort_factor(F2, (void*)&absi_cmp, cmp_nodata);
    3846                 :        441 :     F = merge_factor(F, F2, (void*)&absi_cmp, cmp_nodata);
    3847                 :            :   }
    3848                 :        679 :   return gerepilecopy(av, F);
    3849                 :            : }

Generated by: LCOV version 1.9