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-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 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.18.1 lcov report (development 29884-9d10e5cf54) Lines: 1806 2210 81.7 %
Date: 2025-01-20 09:10:26 Functions: 104 119 87.4 %
Legend: Lines: hit not hit

          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; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : #define DEBUGLEVEL DEBUGLEVEL_factorint
      18             : 
      19             : /***********************************************************************/
      20             : /**                       PRIMES IN SUCCESSION                        **/
      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 /* nonprime 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             : /* first differences of the preceding */
      42             : static unsigned char prc210_d1[] = {
      43             :   10, 2, 4, 2, 4, 6, 2, 6, 4, 2, 4, 6, 6, 2, 6, 4, 2, 6,
      44             :   4, 6, 8, 4, 2, 4, 2, 4, 8, 6, 4, 6, 2, 4, 6,
      45             :   2, 6, 6, 4, 2, 4, 6, 2, 6, 4, 2, 4, 2, 10, 2,
      46             : };
      47             : 
      48             : static int
      49      984000 : unextprime_overflow(ulong n)
      50             : {
      51             : #ifdef LONG_IS_64BIT
      52      983414 :   return (n > (ulong)-59);
      53             : #else
      54         586 :   return (n > (ulong)-5);
      55             : #endif
      56             : }
      57             : 
      58             : /* return 0 for overflow */
      59             : ulong
      60     1121103 : unextprime(ulong n)
      61             : {
      62             :   long rc, rc0, rcn;
      63             : 
      64     1121103 :   switch(n) {
      65        6858 :     case 0: case 1: case 2: return 2;
      66        2434 :     case 3: return 3;
      67        1668 :     case 4: case 5: return 5;
      68        1162 :     case 6: case 7: return 7;
      69             :   }
      70     1108981 :   if (n <= maxprime())
      71             :   {
      72      124964 :     long i = PRIMES_search(n);
      73      124964 :     return i > 0? n: pari_PRIMES[-i];
      74             :   }
      75      984005 :   if (unextprime_overflow(n)) return 0;
      76             :   /* here n > 7 */
      77      983971 :   n |= 1; /* make it odd */
      78      983971 :   rc = rc0 = n % 210;
      79             :   /* find next prime residue class mod 210 */
      80             :   for(;;)
      81             :   {
      82     2166691 :     rcn = (long)(prc210_no[rc>>1]);
      83     2166691 :     if (rcn != NPRC) break;
      84     1182720 :     rc += 2; /* cannot wrap since 209 is coprime and rc odd */
      85             :   }
      86      983971 :   if (rc > rc0) n += rc - rc0;
      87             :   /* now find an actual (pseudo)prime */
      88             :   for(;;)
      89             :   {
      90    11530693 :     if (uisprime(n)) break;
      91    10546722 :     n += prc210_d1[rcn];
      92    10546722 :     if (++rcn > 47) rcn = 0;
      93             :   }
      94      983998 :   return n;
      95             : }
      96             : 
      97             : GEN
      98      126847 : nextprime(GEN n)
      99             : {
     100             :   long rc, rc0, rcn;
     101      126847 :   pari_sp av = avma;
     102             : 
     103      126847 :   if (typ(n) != t_INT)
     104             :   {
     105          14 :     n = gceil(n);
     106          14 :     if (typ(n) != t_INT) pari_err_TYPE("nextprime",n);
     107             :   }
     108      126840 :   if (signe(n) <= 0) { set_avma(av); return gen_2; }
     109      126840 :   if (lgefint(n) == 3)
     110             :   {
     111      118710 :     ulong k = unextprime(uel(n,2));
     112      118710 :     set_avma(av);
     113      118710 :     if (k) return utoipos(k);
     114             : #ifdef LONG_IS_64BIT
     115           6 :     return uutoi(1,13);
     116             : #else
     117           1 :     return uutoi(1,15);
     118             : #endif
     119             :   }
     120             :   /* here n > 7 */
     121        8130 :   if (!mod2(n)) n = addui(1,n);
     122        8130 :   rc = rc0 = umodiu(n, 210);
     123             :   /* find next prime residue class mod 210 */
     124             :   for(;;)
     125             :   {
     126       17701 :     rcn = (long)(prc210_no[rc>>1]);
     127       17701 :     if (rcn != NPRC) break;
     128        9571 :     rc += 2; /* cannot wrap since 209 is coprime and rc odd */
     129             :   }
     130        8130 :   if (rc > rc0) n = addui(rc - rc0, n);
     131             :   /* now find an actual (pseudo)prime */
     132             :   for(;;)
     133             :   {
     134       84466 :     if (BPSW_psp(n)) break;
     135       76336 :     n = addui(prc210_d1[rcn], n);
     136       76336 :     if (++rcn > 47) rcn = 0;
     137             :   }
     138        8130 :   if (avma == av) return icopy(n);
     139        8130 :   return gerepileuptoint(av, n);
     140             : }
     141             : 
     142             : ulong
     143          32 : uprecprime(ulong n)
     144             : {
     145             :   long rc, rc0, rcn;
     146             :   { /* check if n <= 10 */
     147          32 :     if (n <= 1)  return 0;
     148          25 :     if (n == 2)  return 2;
     149          18 :     if (n <= 4)  return 3;
     150          18 :     if (n <= 6)  return 5;
     151          18 :     if (n <= 10) return 7;
     152             :   }
     153          18 :   if (n <= maxprimelim())
     154             :   {
     155           0 :     long i = PRIMES_search(n);
     156           0 :     return i > 0? n: pari_PRIMES[-i-1];
     157             :   }
     158             :   /* here n >= 11 */
     159          18 :   if (!(n % 2)) n--;
     160          18 :   rc = rc0 = n % 210;
     161             :   /* find previous prime residue class mod 210 */
     162             :   for(;;)
     163             :   {
     164          36 :     rcn = (long)(prc210_no[rc>>1]);
     165          36 :     if (rcn != NPRC) break;
     166          18 :     rc -= 2; /* cannot wrap since 1 is coprime and rc odd */
     167             :   }
     168          18 :   if (rc < rc0) n += rc - rc0;
     169             :   /* now find an actual (pseudo)prime */
     170             :   for(;;)
     171             :   {
     172          36 :     if (uisprime(n)) break;
     173          18 :     if (--rcn < 0) rcn = 47;
     174          18 :     n -= prc210_d1[rcn];
     175             :   }
     176          18 :   return n;
     177             : }
     178             : 
     179             : GEN
     180          49 : precprime(GEN n)
     181             : {
     182             :   long rc, rc0, rcn;
     183          49 :   pari_sp av = avma;
     184             : 
     185          49 :   if (typ(n) != t_INT)
     186             :   {
     187          14 :     n = gfloor(n);
     188          14 :     if (typ(n) != t_INT) pari_err_TYPE("nextprime",n);
     189             :   }
     190          42 :   if (signe(n) <= 0) { set_avma(av); return gen_0; }
     191          42 :   if (lgefint(n) <= 3)
     192             :   {
     193          32 :     ulong k = uel(n,2);
     194          32 :     return gc_utoi(av, uprecprime(k));
     195             :   }
     196          10 :   if (!mod2(n)) n = subiu(n,1);
     197          10 :   rc = rc0 = umodiu(n, 210);
     198             :   /* find previous prime residue class mod 210 */
     199             :   for(;;)
     200             :   {
     201          20 :     rcn = (long)(prc210_no[rc>>1]);
     202          20 :     if (rcn != NPRC) break;
     203          10 :     rc -= 2; /* cannot wrap since 1 is coprime and rc odd */
     204             :   }
     205          10 :   if (rc0 > rc) n = subiu(n, rc0 - rc);
     206             :   /* now find an actual (pseudo)prime */
     207             :   for(;;)
     208             :   {
     209          48 :     if (BPSW_psp(n)) break;
     210          38 :     if (--rcn < 0) rcn = 47;
     211          38 :     n = subiu(n, prc210_d1[rcn]);
     212             :   }
     213          10 :   if (avma == av) return icopy(n);
     214          10 :   return gerepileuptoint(av, n);
     215             : }
     216             : 
     217             : /* Find next single-word prime strictly larger than p.
     218             :  * If *n < pari_PRIMES[0], p is *n-th prime, 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 prime 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             : static ulong
     226     4531317 : snextpr(ulong p, long *n, long *rcn, long *q, int (*ispsp)(ulong))
     227             : {
     228     4531317 :   if (*n < pari_PRIMES[0])
     229             :   {
     230     4531317 :     ulong t, p1 = t = pari_PRIMES[++*n]; /* nextprime(p + 1) */
     231     4531317 :     if (*rcn != NPRC)
     232             :     {
     233    15888894 :       while (t > p)
     234             :       {
     235    11373946 :         t -= prc210_d1[*rcn];
     236    11373946 :         if (++*rcn > 47) { *rcn = 0; if (q) (*q)++; }
     237             :       }
     238             :       /* assert(d1 == p) */
     239             :     }
     240     4531317 :     return p1;
     241             :   }
     242           0 :   if (unextprime_overflow(p)) pari_err_OVERFLOW("snextpr");
     243             :   /* we are beyond the prime table, initialize if needed */
     244           0 :   if (*rcn == NPRC) *rcn = prc210_no[(p % 210) >> 1]; /* != NPRC */
     245             :   /* look for the next one */
     246             :   do {
     247           0 :     p += prc210_d1[*rcn];
     248           0 :     if (++*rcn > 47) { *rcn = 0; if (q) (*q)++; }
     249           0 :   } while (!ispsp(p));
     250           0 :   return p;
     251             : }
     252             : 
     253             : /********************************************************************/
     254             : /**                                                                **/
     255             : /**                     INTEGER FACTORIZATION                      **/
     256             : /**                                                                **/
     257             : /********************************************************************/
     258             : int factor_add_primes = 0, factor_proven = 0;
     259             : 
     260             : /***********************************************************************/
     261             : /**                                                                   **/
     262             : /**                 FACTORIZATION (ECM) -- GN Jul-Aug 1998            **/
     263             : /**   Integer factorization using the elliptic curves method (ECM).   **/
     264             : /**   ellfacteur() returns a non trivial factor of N, assuming N>0,   **/
     265             : /**   is composite, and has no prime divisor below tridiv_bound(N)    **/
     266             : /**   Thanks to Paul Zimmermann for much helpful advice and to        **/
     267             : /**   Guillaume Hanrot and Igor Schein for intensive testing          **/
     268             : /**                                                                   **/
     269             : /***********************************************************************/
     270             : #define nbcmax 64 /* max number of simultaneous curves */
     271             : 
     272             : static const ulong TB1[] = {
     273             :   142,172,208,252,305,370,450,545,661,801,972,1180,1430,
     274             :   1735,2100,2550,3090,3745,4540,5505,6675,8090,9810,11900,
     275             :   14420,17490,21200,25700,31160,37780UL,45810UL,55550UL,67350UL,
     276             :   81660UL,99010UL,120050UL,145550UL,176475UL,213970UL,259430UL,
     277             :   314550UL,381380UL,462415UL,560660UL,679780UL,824220UL,999340UL,
     278             :   1211670UL,1469110UL,1781250UL,2159700UL,2618600UL,3175000UL,
     279             :   3849600UL,4667500UL,5659200UL,6861600UL,8319500UL,10087100UL,
     280             :   12230300UL,14828900UL,17979600UL,21799700UL,26431500UL,
     281             :   32047300UL,38856400UL, /* 110 times that still fits into 32bits */
     282             : #ifdef LONG_IS_64BIT
     283             :   47112200UL,57122100UL,69258800UL,83974200UL,101816200UL,
     284             :   123449000UL,149678200UL,181480300UL,220039400UL,266791100UL,
     285             :   323476100UL,392204900UL,475536500UL,576573500UL,699077800UL,
     286             :   847610500UL,1027701900UL,1246057200UL,1510806400UL,1831806700UL,
     287             :   2221009800UL,2692906700UL,3265067200UL,3958794400UL,4799917500UL
     288             : #endif
     289             : };
     290             : static const ulong TB1_for_stage[] = {
     291             :  /* Start below the optimal B1 for finding factors which would just have been
     292             :   * missed by pollardbrent(), and escalate, changing curves to give good
     293             :   * coverage of the small factor ranges. Entries grow faster than what would
     294             :   * be optimal but a table instead of a 2D array keeps the code simple */
     295             :   500,520,560,620,700,800,900,1000,1150,1300,1450,1600,1800,2000,
     296             :   2200,2450,2700,2950,3250,3600,4000,4400,4850,5300,5800,6400,
     297             :   7100,7850,8700,9600,10600,11700,12900,14200,15700,17300,
     298             :   19000,21000,23200,25500,28000,31000,34500UL,38500UL,43000UL,
     299             :   48000UL,53800UL,60400UL,67750UL,76000UL,85300UL,95700UL,
     300             :   107400UL,120500UL,135400UL,152000UL,170800UL,191800UL,215400UL,
     301             :   241800UL,271400UL,304500UL,341500UL,383100UL,429700UL,481900UL,
     302             :   540400UL,606000UL,679500UL,761800UL,854100UL,957500UL,1073500UL
     303             : };
     304             : 
     305             : /* addition/doubling/multiplication of a point on an 'elliptic curve mod N'
     306             :  * may result in one of three things:
     307             :  * - a new bona fide point
     308             :  * - a point at infinity (denominator divisible by N)
     309             :  * - a point at infinity mod some p | N but finite mod q | N betraying itself
     310             :  *   by a denominator which has nontrivial gcd with N.
     311             :  *
     312             :  * In the second case, addition/doubling aborts, copying one of the summands
     313             :  * to the destination array of points unless they coincide.
     314             :  * Multiplication will stop at some unpredictable intermediate stage:  The
     315             :  * destination will contain _some_ multiple of the input point, but not
     316             :  * necessarily the desired one, which doesn't matter.  As long as we're
     317             :  * multiplying (B1 phase) we simply carry on with the next multiplier.
     318             :  * During the B2 phase, the only additions are the giant steps, and the
     319             :  * worst that can happen here is that we lose one residue class mod 210
     320             :  * of prime multipliers on 4 of the curves, so again, we ignore the problem
     321             :  * and just carry on.)
     322             :  *
     323             :  * Idea: select nbc curves mod N and one point P on each of them. For each
     324             :  * such P, compute [M]P = Q where M is the product of all powers <= B2 of
     325             :  * primes <= nextprime(B1). Then check whether [p]Q for p < nextprime(B2)
     326             :  * betrays a factor. This second stage looks separately at the primes in
     327             :  * each residue class mod 210, four curves at a time, and steps additively
     328             :  * to ever larger multipliers, by comparing X coordinates of points which we
     329             :  * would need to add in order to reach another prime multiplier in the same
     330             :  * residue class. 'Comparing' means that we accumulate a product of
     331             :  * differences of X coordinates, and from time to time take a gcd of this
     332             :  * product with N. Montgomery's multi-inverse trick is used heavily. */
     333             : 
     334             : /* *** auxiliary functions for ellfacteur: *** */
     335             : /* (Rx,Ry) <- (Px,Py)+(Qx,Qy) over Z/NZ, z=1/(Px-Qx). If Ry = NULL, don't set */
     336             : static void
     337     8291496 : FpE_add_i(GEN N, GEN z, GEN Px, GEN Py, GEN Qx, GEN Qy, GEN *Rx, GEN *Ry)
     338             : {
     339     8291496 :   GEN slope = modii(mulii(subii(Py, Qy), z), N);
     340     8291496 :   GEN t = subii(sqri(slope), addii(Qx, Px));
     341     8291496 :   affii(modii(t, N), *Rx);
     342     8291496 :   if (Ry) {
     343     8219188 :     t = subii(mulii(slope, subii(Px, *Rx)), Py);
     344     8219188 :     affii(modii(t, N), *Ry);
     345             :   }
     346     8291496 : }
     347             : /* X -> Z; cannot add on one of the curves: make sure Z contains
     348             :  * something useful before letting caller proceed */
     349             : static void
     350       25650 : ZV_aff(long n, GEN *X, GEN *Z)
     351             : {
     352       25650 :   if (X != Z) {
     353             :     long k;
     354     1569330 :     for (k = n; k--; ) affii(X[k],Z[k]);
     355             :   }
     356       25650 : }
     357             : 
     358             : /* Parallel addition on nbc curves, assigning the result to locations at and
     359             :  * following *X3, *Y3. (If Y-coords of result not desired, set Y=NULL.)
     360             :  * Safe even if (X3,Y3) = (X2,Y2), _not_ if (X1,Y1). It is also safe to
     361             :  * overwrite Y2 with X3. If nbc1 < nbc, the first summand is
     362             :  * assumed to hold only nbc1 distinct points, repeated as often as we need
     363             :  * them  (to add one point on each of a few curves to several other points on
     364             :  * the same curves): only used with nbc1 = nbc or nbc1 = 4 | nbc.
     365             :  *
     366             :  * Return 0 [SUCCESS], 1 [N | den], 2 [gcd(den, N) is a factor of N, preserved
     367             :  * in gl.
     368             :  * Stack space is bounded by a constant multiple of lgefint(N)*nbc:
     369             :  * - Phase 2 creates 12 items on the stack per iteration, of which 4 are twice
     370             :  *   as long and 1 is thrice as long as N, i.e. 18 units per iteration.
     371             :  * - Phase  1 creates 4 units.
     372             :  * Total can be as large as 4*nbcmax + 18*8 units; ecm_elladd2() is
     373             :  * just as bad, and elldouble() comes to 3*nbcmax + 29*8 units. */
     374             : static int
     375      240431 : ecm_elladd0(GEN N, GEN *gl, long nbc, long nbc1,
     376             :             GEN *X1, GEN *Y1, GEN *X2, GEN *Y2, GEN *X3, GEN *Y3)
     377             : {
     378      240431 :   const ulong mask = (nbc1 == 4)? 3: ~0UL; /*nbc1 = 4 or nbc*/
     379      240431 :   GEN W[2*nbcmax], *A = W+nbc; /* W[0],A[0] unused */
     380             :   long i;
     381      240431 :   pari_sp av = avma;
     382             : 
     383      240431 :   W[1] = subii(X1[0], X2[0]);
     384     7825896 :   for (i=1; i<nbc; i++)
     385             :   { /*prepare for multi-inverse*/
     386     7585465 :     A[i] = subii(X1[i&mask], X2[i]); /* don't waste time reducing mod N */
     387     7585465 :     W[i+1] = modii(mulii(A[i], W[i]), N);
     388             :   }
     389      240431 :   if (!invmod(W[nbc], N, gl))
     390             :   {
     391          18 :     if (!equalii(N,*gl)) return 2;
     392           0 :     ZV_aff(nbc, X2,X3);
     393           0 :     if (Y3) ZV_aff(nbc, Y2,Y3);
     394           0 :     return gc_int(av,1);
     395             :   }
     396             : 
     397     7825032 :   while (i--) /* nbc times */
     398             :   {
     399     7825032 :     pari_sp av2 = avma;
     400     7825032 :     GEN Px = X1[i&mask], Py = Y1[i&mask], Qx = X2[i], Qy = Y2[i];
     401     7825032 :     GEN z = i? mulii(*gl,W[i]): *gl; /*1/(Px-Qx)*/
     402     7825032 :     FpE_add_i(N,z,  Px,Py,Qx,Qy, X3+i, Y3? Y3+i: NULL);
     403     7825032 :     if (!i) break;
     404     7584619 :     set_avma(av2); *gl = modii(mulii(*gl, A[i]), N);
     405             :   }
     406      240413 :   return gc_int(av,0);
     407             : }
     408             : 
     409             : /* Shortcut, for use in cases where Y coordinates follow their corresponding
     410             :  * X coordinates, and first summand doesn't need to be repeated */
     411             : static int
     412      234487 : ecm_elladd(GEN N, GEN *gl, long nbc, GEN *X1, GEN *X2, GEN *X3) {
     413      234487 :   return ecm_elladd0(N, gl, nbc, nbc, X1, X1+nbc, X2, X2+nbc, X3, X3+nbc);
     414             : }
     415             : 
     416             : /* As ecm_elladd except it does twice as many additions (and hides even more
     417             :  * of the cost of the modular inverse); the net effect is the same as
     418             :  * ecm_elladd(nbc,X1,X2,X3) && ecm_elladd(nbc,X4,X5,X6). Safe to
     419             :  * have X2=X3, X5=X6, or X1,X2 coincide with X4,X5 in any order. */
     420             : static int
     421        7194 : ecm_elladd2(GEN N, GEN *gl, long nbc,
     422             :             GEN *X1, GEN *X2, GEN *X3, GEN *X4, GEN *X5, GEN *X6)
     423             : {
     424        7194 :   GEN *Y1 = X1+nbc, *Y2 = X2+nbc, *Y3 = X3+nbc;
     425        7194 :   GEN *Y4 = X4+nbc, *Y5 = X5+nbc, *Y6 = X6+nbc;
     426        7194 :   GEN W[4*nbcmax], *A = W+2*nbc; /* W[0],A[0] unused */
     427             :   long i, j;
     428        7194 :   pari_sp av = avma;
     429             : 
     430        7194 :   W[1] = subii(X1[0], X2[0]);
     431      233232 :   for (i=1; i<nbc; i++)
     432             :   {
     433      226038 :     A[i] = subii(X1[i], X2[i]); /* don't waste time reducing mod N here */
     434      226038 :     W[i+1] = modii(mulii(A[i], W[i]), N);
     435             :   }
     436      240426 :   for (j=0; j<nbc; i++,j++)
     437             :   {
     438      233232 :     A[i] = subii(X4[j], X5[j]);
     439      233232 :     W[i+1] = modii(mulii(A[i], W[i]), N);
     440             :   }
     441        7194 :   if (!invmod(W[2*nbc], N, gl))
     442             :   {
     443           0 :     if (!equalii(N,*gl)) return 2;
     444           0 :     ZV_aff(2*nbc, X2,X3); /* hack: 2*nbc => copy Y2->Y3 */
     445           0 :     ZV_aff(2*nbc, X5,X6); /* also copy Y5->Y6 */
     446           0 :     return gc_int(av,1);
     447             :   }
     448             : 
     449      240426 :   while (j--) /* nbc times */
     450             :   {
     451      233232 :     pari_sp av2 = avma;
     452      233232 :     GEN Px = X4[j], Py = Y4[j], Qx = X5[j], Qy = Y5[j];
     453      233232 :     GEN z = mulii(*gl,W[--i]); /*1/(Px-Qx)*/
     454      233232 :     FpE_add_i(N,z, Px,Py, Qx,Qy, X6+j,Y6+j);
     455      233232 :     set_avma(av2); *gl = modii(mulii(*gl, A[i]), N);
     456             :   }
     457      233232 :   while (i--) /* nbc times */
     458             :   {
     459      233232 :     pari_sp av2 = avma;
     460      233232 :     GEN Px = X1[i], Py = Y1[i], Qx = X2[i], Qy = Y2[i];
     461      233232 :     GEN z = i? mulii(*gl, W[i]): *gl; /*1/(Px-Qx)*/
     462      233232 :     FpE_add_i(N,z, Px,Py, Qx,Qy, X3+i,Y3+i);
     463      233232 :     if (!i) break;
     464      226038 :     set_avma(av2); *gl = modii(mulii(*gl, A[i]), N);
     465             :   }
     466        7194 :   return gc_int(av,0);
     467             : }
     468             : 
     469             : /* Parallel doubling on nbc curves, assigning the result to locations at
     470             :  * and following *X2.  Safe to be called with X2 equal to X1.  Return
     471             :  * value as for ecm_elladd.  If we find a point at infinity mod N,
     472             :  * and if X1 != X2, we copy the points at X1 to X2. */
     473             : static int
     474       40073 : elldouble(GEN N, GEN *gl, long nbc, GEN *X1, GEN *X2)
     475             : {
     476       40073 :   GEN *Y1 = X1+nbc, *Y2 = X2+nbc;
     477             :   GEN W[nbcmax+1]; /* W[0] unused */
     478             :   long i;
     479       40073 :   pari_sp av = avma;
     480       40073 :   /*W[0] = gen_1;*/ W[1] = Y1[0];
     481     1233544 :   for (i=1; i<nbc; i++) W[i+1] = modii(mulii(Y1[i], W[i]), N);
     482       40073 :   if (!invmod(W[nbc], N, gl))
     483             :   {
     484           0 :     if (!equalii(N,*gl)) return 2;
     485           0 :     ZV_aff(2*nbc,X1,X2); /* also copies Y1->Y2 */
     486           0 :     return gc_int(av,1);
     487             :   }
     488     1273617 :   while (i--) /* nbc times */
     489             :   {
     490             :     pari_sp av2;
     491     1233544 :     GEN v, w, L, z = i? mulii(*gl,W[i]): *gl;
     492     1233544 :     if (i) *gl = modii(mulii(*gl, Y1[i]), N);
     493     1233544 :     av2 = avma;
     494     1233544 :     L = modii(mulii(addui(1, mului(3, Fp_sqr(X1[i],N))), z), N);
     495     1233544 :     if (signe(L)) /* half of zero is still zero */
     496     1233544 :       L = shifti(mod2(L)? addii(L, N): L, -1);
     497     1233544 :     v = modii(subii(sqri(L), shifti(X1[i],1)), N);
     498     1233544 :     w = modii(subii(mulii(L, subii(X1[i], v)), Y1[i]), N);
     499     1233544 :     affii(v, X2[i]);
     500     1233544 :     affii(w, Y2[i]);
     501     1233544 :     set_avma(av2);
     502             :   }
     503       40073 :   return gc_int(av,0);
     504             : }
     505             : 
     506             : /* Parallel multiplication by an odd prime k on nbc curves, storing the
     507             :  * result to locations at and following *X2. Safe to be called with X2 = X1.
     508             :  * Return values as ecm_elladd. Uses (a simplified variant of) Montgomery's
     509             :  * PRAC algorithm; see ftp://ftp.cwi.nl/pub/pmontgom/Lucas.ps.gz .
     510             :  * With thanks to Paul Zimmermann for the reference.  --GN1998Aug13 */
     511             : static int
     512      208727 : get_rule(ulong d, ulong e)
     513             : {
     514      208727 :   if (d <= e + (e>>2)) /* floor(1.25*e) */
     515             :   {
     516       16630 :     if ((d+e)%3 == 0) return 0; /* rule 1 */
     517        9928 :     if ((d-e)%6 == 0) return 1;  /* rule 2 */
     518             :   }
     519             :   /* d <= 4*e but no ofl */
     520      201971 :   if ((d+3)>>2 <= e) return 2; /* rule 3, common case */
     521       12148 :   if ((d&1)==(e&1))  return 1; /* rule 4 = rule 2 */
     522        6344 :   if (!(d&1))        return 3; /* rule 5 */
     523        1769 :   if (d%3 == 0)      return 4; /* rule 6 */
     524         417 :   if ((d+e)%3 == 0)  return 5; /* rule 7 */
     525           0 :   if ((d-e)%3 == 0)  return 6; /* rule 8 */
     526             :   /* when we get here, e is even, otherwise one of rules 4,5 would apply */
     527           0 :   return 7; /* rule 9 */
     528             : }
     529             : 
     530             : /* PRAC implementation notes - main changes against the paper version:
     531             :  * (1) The general function [m+n]P = f([m]P,[n]P,[m-n]P) collapses (for m!=n)
     532             :  * to an ecm_elladd() which does not depend on the third argument; thus
     533             :  * references to the third variable (C in the paper) can be eliminated.
     534             :  * (2) Since our multipliers are prime, the outer loop of the paper
     535             :  * version executes only once, and thus is invisible above.
     536             :  * (3) The first step in the inner loop of the paper version will always be
     537             :  * rule 3, but the addition requested by this rule amounts to a doubling, and
     538             :  * will always be followed by a swap, so we have unrolled this first iteration.
     539             :  * (4) Simplifications in rules 6 and 7 are possible given the above, and we
     540             :  * save one addition in each of the two cases.  NB none of the other
     541             :  * ecm_elladd()s in the loop can ever degenerate into an elldouble.
     542             :  * (5) I tried to optimize for rule 3, which is used more frequently than all
     543             :  * others together, but it didn't improve things, so I removed the nested
     544             :  * tight loop again.  --GN */
     545             : /* The main loop body of ellfacteur() runs _slower_ under PRAC than under a
     546             :  * straightforward left-shift binary multiplication when N has <30 digits and
     547             :  * B1 is small;  PRAC wins when N and B1 get larger.  Weird. --GN */
     548             : /* k>2 assumed prime, XAUX = scratchpad */
     549             : static int
     550       25650 : ellmult(GEN N, GEN *gl, long nbc, ulong k, GEN *X1, GEN *X2, GEN *XAUX)
     551             : {
     552             :   ulong r, d, e, e1;
     553             :   int res;
     554       25650 :   GEN *A = X2, *B = XAUX, *T = XAUX + 2*nbc;
     555             : 
     556       25650 :   ZV_aff(2*nbc,X1,XAUX);
     557             :   /* first doubling picks up X1;  after this we'll be working in XAUX and
     558             :    * X2 only, mostly via A and B and T */
     559       25650 :   if ((res = elldouble(N, gl, nbc, X1, X2)) != 0) return res;
     560             : 
     561             :   /* split the work at the golden ratio */
     562       25650 :   r = (ulong)(k*0.61803398875 + .5);
     563       25650 :   d = k - r;
     564       25650 :   e = r - d; /* d+e == r, so no danger of ofl below */
     565      234377 :   while (d != e)
     566             :   { /* apply one of the nine transformations from PM's Table 4. */
     567      208727 :     switch(get_rule(d,e))
     568             :     {
     569        6702 :     case 0: /* rule 1 */
     570        6702 :       if ( (res = ecm_elladd(N, gl, nbc, A, B, T)) ) return res;
     571        6702 :       if ( (res = ecm_elladd2(N, gl, nbc, T, A, A, T, B, B)) != 0) return res;
     572        6702 :       e1 = d - e; d = (d + e1)/3; e = (e - e1)/3; break;
     573        5858 :     case 1: /* rules 2 and 4 */
     574        5858 :       if ( (res = ecm_elladd(N, gl, nbc, A, B, B)) ) return res;
     575        5858 :       if ( (res = elldouble(N, gl, nbc, A, A)) ) return res;
     576        5858 :       d = (d-e)>>1; break;
     577        4575 :     case 3: /* rule 5 */
     578        4575 :       if ( (res = elldouble(N, gl, nbc, A, A)) ) return res;
     579        4575 :       d >>= 1; break;
     580        1352 :     case 4: /* rule 6 */
     581        1352 :       if ( (res = elldouble(N, gl, nbc, A, T)) ) return res;
     582        1352 :       if ( (res = ecm_elladd(N, gl, nbc, T, A, A)) ) return res;
     583        1352 :       if ( (res = ecm_elladd(N, gl, nbc, A, B, B)) ) return res;
     584        1352 :       d = d/3 - e; break;
     585      189823 :     case 2: /* rule 3 */
     586      189823 :       if ( (res = ecm_elladd(N, gl, nbc, A, B, B)) ) return res;
     587      189823 :       d -= e; break;
     588         417 :     case 5: /* rule 7 */
     589         417 :       if ( (res = elldouble(N, gl, nbc, A, T)) ) return res;
     590         417 :       if ( (res = ecm_elladd2(N, gl, nbc, T, A, A, T, B, B)) != 0) return res;
     591         417 :       d = (d - 2*e)/3; break;
     592           0 :     case 6: /* rule 8 */
     593           0 :       if ( (res = ecm_elladd(N, gl, nbc, A, B, B)) ) return res;
     594           0 :       if ( (res = elldouble(N, gl, nbc, A, T)) ) return res;
     595           0 :       if ( (res = ecm_elladd(N, gl, nbc, T, A, A)) ) return res;
     596           0 :       d = (d - e)/3; break;
     597           0 :     case 7: /* rule 9 */
     598           0 :       if ( (res = elldouble(N, gl, nbc, B, B)) ) return res;
     599           0 :       e >>= 1; break;
     600             :     }
     601             :     /* swap d <-> e and A <-> B if necessary */
     602      208727 :     if (d < e) { lswap(d,e); pswap(A,B); }
     603             :   }
     604       25650 :   return ecm_elladd(N, gl, nbc, XAUX, X2, X2);
     605             : }
     606             : 
     607             : struct ECM {
     608             :   pari_timer T;
     609             :   long nbc, nbc2, seed;
     610             :   GEN *X, *XAUX, *XT, *XD, *XB, *XB2, *XH, *Xh, *Yh;
     611             : };
     612             : 
     613             : /* memory layout in ellfacteur():  a large array of GEN pointers, and one
     614             :  * huge chunk of memory containing all the actual GEN (t_INT) objects.
     615             :  * nbc is constant throughout the invocation:
     616             :  * - The B1 stage of each iteration through the main loop needs little
     617             :  * space:  enough for the X and Y coordinates of the current points,
     618             :  * and twice as much again as scratchpad for ellmult().
     619             :  * - The B2 stage, starting from some current set of points Q, needs, in
     620             :  * succession:
     621             :  *   + space for [2]Q, [4]Q, ..., [10]Q, and [p]Q for building the helix;
     622             :  *   + space for 48*nbc X and Y coordinates to hold the helix.  This could
     623             :  *   re-use [2]Q,...,[8]Q, but only with difficulty, since we don't
     624             :  *   know in advance which residue class mod 210 our p is going to be in.
     625             :  *   It can and should re-use [p]Q, though;
     626             :  *   + space for (temporarily [30]Q and then) [210]Q, [420]Q, and several
     627             :  *   further doublings until the giant step multiplier is reached.  This
     628             :  *   can re-use the remaining cells from above.  The computation of [210]Q
     629             :  *   will have been the last call to ellmult() within this iteration of the
     630             :  *   main loop, so the scratchpad is now also free to be re-used. We also
     631             :  *   compute [630]Q by a parallel addition;  we'll need it later to get the
     632             :  *   baby-step table bootstrapped a little faster.
     633             :  *   + Finally, for no more than 4 curves at a time, room for up to 1024 X
     634             :  *   coordinates only: the Y coordinates needed whilst setting up this baby
     635             :  *   step table are temporarily stored in the upper half, and overwritten
     636             :  *   during the last series of additions.
     637             :  *
     638             :  * Graphically:  after end of B1 stage (X,Y are the coords of Q):
     639             :  * +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--
     640             :  * | X Y |  scratch  | [2]Q| [4]Q| [6]Q| [8]Q|[10]Q|    ...    | ...
     641             :  * +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--
     642             :  * *X    *XAUX *XT   *XD                                       *XB
     643             :  *
     644             :  * [30]Q is computed from [10]Q.  [210]Q can go into XY, etc:
     645             :  * +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--
     646             :  * |[210]|[420]|[630]|[840]|[1680,3360,6720,...,2048*210]      |bstp table...
     647             :  * +--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--
     648             :  * *X    *XAUX *XT   *XD      [*XG, somewhere here]            *XB .... *XH
     649             :  *
     650             :  * So we need (13 + 48) * 2 * nbc slots here + 4096 slots for the baby step
     651             :  * table (not all of which will be used when we start with a small B1, but
     652             :  * better to allocate and initialize ahead of time all the slots that might
     653             :  * be needed later).
     654             :  *
     655             :  * Note on memory locality:  During the B2 phase, accesses to the helix
     656             :  * (once it is set up) will be clustered by curves (4 out of nbc at a time).
     657             :  * Accesses to the baby steps table will wander from one end of the array to
     658             :  * the other and back, one such cycle per giant step, and during a full cycle
     659             :  * we would expect on the order of 2E4 accesses when using the largest giant
     660             :  * step size.  Thus we shouldn't be doing too bad with respect to thrashing
     661             :  * a 512KBy L2 cache.  However, we don't want the baby step table to grow
     662             :  * larger than this, even if it would reduce the number of EC operations by a
     663             :  * few more per cent for very large B2, lest cache thrashing slow down
     664             :  * everything disproportionally. --GN */
     665             : /* Auxiliary routines need < (3*nbc+240)*lN words on the PARI stack, in
     666             :  * addition to the spc*(lN+1) words occupied by our main table. */
     667             : static void
     668          57 : ECM_alloc(struct ECM *E, long lN)
     669             : {
     670          57 :   const long bstpmax = 1024; /* max number of baby step table entries */
     671          57 :   long spc = (13 + 48) * E->nbc2 + bstpmax * 4;
     672          57 :   long len = spc + 385 + spc*lN;
     673          57 :   long i, tw = _evallg(lN) | evaltyp(t_INT);
     674          57 :   GEN w, *X = (GEN*)new_chunk(len);
     675             :   /* hack for X[i] = cgeti(lN). X = current point in B1 phase */
     676          57 :   w = (GEN)(X + spc + 385);
     677      377001 :   for (i = spc-1; i >= 0; i--) { X[i] = w; *w = tw; w += lN; }
     678          57 :   E->X = X;
     679          57 :   E->XAUX = E->X    + E->nbc2; /* scratchpad for ellmult() */
     680          57 :   E->XT   = E->XAUX + E->nbc2; /* ditto, will later hold [3*210]Q */
     681          57 :   E->XD   = E->XT   + E->nbc2; /* room for various multiples */
     682          57 :   E->XB   = E->XD   + 10*E->nbc2; /* start of baby steps table */
     683          57 :   E->XB2  = E->XB   + 2 * bstpmax; /* middle of baby steps table */
     684          57 :   E->XH   = E->XB2  + 2 * bstpmax; /* end of bstps table, start of helix */
     685          57 :   E->Xh   = E->XH   + 48*E->nbc2; /* little helix, X coords */
     686          57 :   E->Yh   = E->XH   + 192;     /* ditto, Y coords */
     687             :   /* XG,YG set inside the main loop, since they depend on B2 */
     688             :   /* E.Xh range of 384 pointers not set; these will later duplicate the pointers
     689             :    * in the E.XH range, 4 curves at a time. Some of the cells reserved here for
     690             :    * the E.XB range will never be used, instead, we'll warp the pointers to
     691             :    * connect to (read-only) GENs in the X/E.XD range */
     692          57 : }
     693             : /* N.B. E->seed is not initialized here */
     694             : static void
     695          57 : ECM_init(struct ECM *E, GEN N, long nbc)
     696             : {
     697          57 :   if (nbc < 0)
     698             :   { /* choose a sensible default */
     699          57 :     const long size = expi(N) + 1;
     700          57 :     nbc = ((size >> 3) << 2) - 80;
     701          57 :     if (nbc < 8) nbc = 8;
     702             :   }
     703          57 :   if (nbc > nbcmax) nbc = nbcmax;
     704          57 :   E->nbc = nbc;
     705          57 :   E->nbc2 = nbc << 1;
     706          57 :   ECM_alloc(E, lgefint(N));
     707          57 : }
     708             : 
     709             : static GEN
     710          93 : ECM_loop(struct ECM *E, GEN N, ulong B1)
     711             : {
     712          93 :   const ulong B2 = 110 * B1, B2_rt = usqrt(B2);
     713          93 :   const ulong nbc = E->nbc, nbc2 = E->nbc2;
     714             :   pari_sp av1, avtmp;
     715             :   long i, np, np0, gse, gss, bstp, bstp0, rcn0, rcn;
     716             :   ulong B2_p, m, p, p0;
     717             :   GEN g, *XG, *YG;
     718          93 :   GEN *X = E->X, *XAUX = E->XAUX, *XT = E->XT, *XD = E->XD;
     719          93 :   GEN *XB = E->XB, *XB2 = E->XB2, *XH = E->XH, *Xh = E->Xh, *Yh = E->Yh;
     720             :   /* pick curves */
     721        4461 :   for (i = nbc2; i--; ) affui(E->seed++, X[i]);
     722             :   /* pick giant step exponent and size */
     723          93 :   gse = B1 < 656
     724             :           ? (B1 < 200? 5: 6)
     725          93 :           : (B1 < 10500
     726             :             ? (B1 < 2625? 7: 8)
     727             :             : (B1 < 42000? 9: 10));
     728          93 :   gss = 1UL << gse;
     729             :   /* With 32 baby steps, a giant step corresponds to 32*420 = 13440,
     730             :    * appropriate for the smallest B2s. With 1024, a giant step will be 430080;
     731             :    * appropriate for B1 >~ 42000, where 512 baby steps would imply roughly
     732             :    * the same number of curve additions. */
     733          93 :   XG = XT + gse*nbc2; /* will later hold [2^(gse+1)*210]Q */
     734          93 :   YG = XG + nbc;
     735             : 
     736          93 :   if (DEBUGLEVEL >= 4) {
     737           0 :     err_printf("ECM: time = %6ld ms\nECM: B1 = %4lu,", timer_delay(&E->T), B1);
     738           0 :     err_printf("\tB2 = %6lu,\tgss = %4ld*420\n", B2, gss);
     739             :   }
     740          93 :   p = 2; np = 1; /* p is np-th prime */
     741             : 
     742             :   /* ---B1 PHASE--- */
     743             :   /* treat p=2 separately */
     744          93 :   B2_p = B2 >> 1;
     745        1603 :   for (m=1; m<=B2_p; m<<=1)
     746             :   {
     747        1510 :     int fl = elldouble(N, &g, nbc, X, X);
     748        1510 :     if (fl > 1) return g; else if (fl) break;
     749             :   }
     750          93 :   rcn = NPRC; /* multipliers begin at the beginning */
     751             :   /* p=3,...,nextprime(B1) */
     752        6538 :   while (p < B1 && p <= B2_rt)
     753             :   {
     754        6445 :     pari_sp av2 = avma;
     755        6445 :     p = snextpr(p, &np, &rcn, NULL, uisprime);
     756        6445 :     B2_p = B2/p; /* beware integer overflow on 32-bit CPUs */
     757       22021 :     for (m=1; m<=B2_p; m*=p)
     758             :     {
     759       15576 :       int fl = ellmult(N, &g, nbc, p, X, X, XAUX);
     760       15576 :       if (fl > 1) return g; else if (fl) break;
     761       15576 :       set_avma(av2);
     762             :     }
     763        6445 :     set_avma(av2);
     764             :   }
     765             :   /* primes p larger than sqrt(B2) appear only to the 1st power */
     766        9924 :   while (p < B1)
     767             :   {
     768        9849 :     pari_sp av2 = avma;
     769        9849 :     p = snextpr(p, &np, &rcn, NULL, uisprime);
     770        9849 :     if (ellmult(N, &g, nbc, p, X, X, XAUX) > 1) return g;
     771        9831 :     set_avma(av2);
     772             :   }
     773          75 :   if (DEBUGLEVEL >= 4) {
     774           0 :     err_printf("ECM: time = %6ld ms, B1 phase done, ", timer_delay(&E->T));
     775           0 :     err_printf("p = %lu, setting up for B2\n", p);
     776             :   }
     777             : 
     778             :   /* ---B2 PHASE--- */
     779             :   /* compute [2]Q,...,[10]Q, needed to build the helix */
     780          75 :   if (elldouble(N, &g, nbc, X, XD) > 1) return g; /*[2]Q*/
     781          75 :   if (elldouble(N, &g, nbc, XD, XD + nbc2) > 1) return g; /*[4]Q*/
     782          75 :   if (ecm_elladd(N, &g, nbc,
     783          75 :         XD, XD + nbc2, XD + (nbc<<2)) > 1) return g; /* [6]Q */
     784          75 :   if (ecm_elladd2(N, &g, nbc,
     785             :         XD, XD + (nbc<<2), XT + (nbc<<3),
     786          75 :         XD + nbc2, XD + (nbc<<2), XD + (nbc<<3)) > 1)
     787           0 :     return g; /* [8]Q and [10]Q */
     788          75 :   if (DEBUGLEVEL >= 7) err_printf("\t(got [2]Q...[10]Q)\n");
     789             : 
     790             :   /* get next prime (still using the foolproof test) */
     791          75 :   p = snextpr(p, &np, &rcn, NULL, uisprime);
     792             :   /* make sure we have the residue class number (mod 210) */
     793          75 :   if (rcn == NPRC)
     794             :   {
     795          75 :     rcn = prc210_no[(p % 210) >> 1];
     796          75 :     if (rcn == NPRC)
     797             :     {
     798           0 :       err_printf("ECM: %lu should have been prime but isn\'t\n", p);
     799           0 :       pari_err_BUG("ellfacteur");
     800             :     }
     801             :   }
     802             : 
     803             :   /* compute [p]Q and put it into its place in the helix */
     804          75 :   if (ellmult(N, &g, nbc, p, X, XH + rcn*nbc2, XAUX) > 1)
     805           0 :     return g;
     806          75 :   if (DEBUGLEVEL >= 7)
     807           0 :     err_printf("\t(got [p]Q, p = %lu = prc210_rp[%ld] mod 210)\n", p, rcn);
     808             : 
     809             :   /* save current p, np, and rcn;  we'll need them more than once below */
     810          75 :   p0 = p; np0 = np; rcn0 = rcn;
     811          75 :   bstp0 = 0; /* p is at baby-step offset 0 from itself */
     812             : 
     813             :   /* fill up the helix, stepping forward through the prime residue classes
     814             :    * mod 210 until we're back at the r'class of p0.  Keep updating p so
     815             :    * that we can print meaningful diagnostics if a factor shows up; don't
     816             :    * bother checking which of these p's are in fact prime */
     817        3600 :   for (i = 47; i; i--) /* 47 iterations */
     818             :   {
     819        3525 :     ulong dp = (ulong)prc210_d1[rcn];
     820        3525 :     p += dp;
     821        3525 :     if (rcn == 47)
     822             :     { /* wrap mod 210 */
     823          75 :       if (ecm_elladd(N, &g, nbc, XT+dp*nbc, XH+rcn*nbc2, XH) > 1) return g;
     824          75 :       rcn = 0; continue;
     825             :     }
     826        3450 :     if (ecm_elladd(N, &g, nbc, XT+dp*nbc, XH+rcn*nbc2, XH+rcn*nbc2+nbc2) > 1)
     827           0 :       return g;
     828        3450 :     rcn++;
     829             :   }
     830          75 :   if (DEBUGLEVEL >= 7) err_printf("\t(got initial helix)\n");
     831             :   /* compute [210]Q etc, needed for the baby step table */
     832          75 :   if (ellmult(N, &g, nbc, 3, XD + (nbc<<3), X, XAUX) > 1) return g;
     833          75 :   if (ellmult(N, &g, nbc, 7, X, X, XAUX) > 1) return g; /* [210]Q */
     834             :   /* this was the last call to ellmult() in the main loop body; may now
     835             :    * overwrite XAUX and slots XD and following */
     836          75 :   if (elldouble(N, &g, nbc, X, XAUX) > 1) return g; /* [420]Q */
     837          75 :   if (ecm_elladd(N, &g, nbc, X, XAUX, XT) > 1) return g;/*[630]Q*/
     838          75 :   if (ecm_elladd(N, &g, nbc, X, XT, XD) > 1) return g;  /*[840]Q*/
     839         561 :   for (i=1; i <= gse; i++)
     840         486 :     if (elldouble(N, &g, nbc, XT + i*nbc2, XD + i*nbc2) > 1) return g;
     841             :   /* (the last iteration has initialized XG to [210*2^(gse+1)]Q) */
     842             : 
     843          75 :   if (DEBUGLEVEL >= 4)
     844           0 :     err_printf("ECM: time = %6ld ms, entering B2 phase, p = %lu\n",
     845             :                timer_delay(&E->T), p);
     846             : 
     847         391 :   for (i = nbc - 4; i >= 0; i -= 4)
     848             :   { /* loop over small sets of 4 curves at a time */
     849             :     GEN *Xb;
     850             :     long j, k;
     851         323 :     if (DEBUGLEVEL >= 6)
     852           0 :       err_printf("ECM: finishing curves %ld...%ld\n", i, i+3);
     853             :     /* Copy relevant pointers from XH to Xh. Memory layout in XH:
     854             :      * nbc X coordinates, nbc Y coordinates for residue class
     855             :      * 1 mod 210, then the same for r.c. 11 mod 210, etc. Memory layout for
     856             :      * Xh is: four X coords for 1 mod 210, four for 11 mod 210, ..., four
     857             :      * for 209 mod 210, then the corresponding Y coordinates in the same
     858             :      * order. This allows a giant step on Xh using just three calls to
     859             :      * ecm_elladd0() each acting on 64 points in parallel */
     860       15827 :     for (j = 48; j--; )
     861             :     {
     862       15504 :       k = nbc2*j + i;
     863       15504 :       m = j << 2; /* X coordinates */
     864       15504 :       Xh[m]   = XH[k];   Xh[m+1] = XH[k+1];
     865       15504 :       Xh[m+2] = XH[k+2]; Xh[m+3] = XH[k+3];
     866       15504 :       k += nbc; /* Y coordinates */
     867       15504 :       Yh[m]   = XH[k];   Yh[m+1] = XH[k+1];
     868       15504 :       Yh[m+2] = XH[k+2]; Yh[m+3] = XH[k+3];
     869             :     }
     870             :     /* Build baby step table of X coords of multiples of [210]Q.  XB[4*j]
     871             :      * will point at X coords on four curves from [(j+1)*210]Q.  Until
     872             :      * we're done, we need some Y coords as well, which we keep in the
     873             :      * second half of the table, overwriting them at the end when gse=10.
     874             :      * Multiples which we already have  (by 1,2,3,4,8,16,...,2^gse) are
     875             :      * entered simply by copying the pointers, ignoring the few slots in w
     876             :      * that were initially reserved for them. Here are the initial entries */
     877         969 :     for (Xb=XB,k=2,j=i; k--; Xb=XB2,j+=nbc) /* first X, then Y coords */
     878             :     {
     879         646 :       Xb[0]  = X[j];      Xb[1]  = X[j+1]; /* [210]Q */
     880         646 :       Xb[2]  = X[j+2];    Xb[3]  = X[j+3];
     881         646 :       Xb[4]  = XAUX[j];   Xb[5]  = XAUX[j+1]; /* [420]Q */
     882         646 :       Xb[6]  = XAUX[j+2]; Xb[7]  = XAUX[j+3];
     883         646 :       Xb[8]  = XT[j];     Xb[9]  = XT[j+1]; /* [630]Q */
     884         646 :       Xb[10] = XT[j+2];   Xb[11] = XT[j+3];
     885         646 :       Xb += 4; /* points at [420]Q */
     886             :       /* ... entries at powers of 2 times 210 .... */
     887        4057 :       for (m = 2; m < (ulong)gse+k; m++) /* omit Y coords of [2^gse*210]Q */
     888             :       {
     889        3411 :         long m2 = m*nbc2 + j;
     890        3411 :         Xb += (2UL<<m); /* points at [2^m*210]Q */
     891        3411 :         Xb[0] = XAUX[m2];   Xb[1] = XAUX[m2+1];
     892        3411 :         Xb[2] = XAUX[m2+2]; Xb[3] = XAUX[m2+3];
     893             :       }
     894             :     }
     895         323 :     if (DEBUGLEVEL >= 7)
     896           0 :       err_printf("\t(extracted precomputed helix / baby step entries)\n");
     897             :     /* ... glue in between, up to 16*210 ... */
     898         323 :     if (ecm_elladd0(N, &g, 12, 4, /* 12 pts + (4 pts replicated thrice) */
     899             :           XB + 12, XB2 + 12,
     900             :           XB,      XB2,
     901           0 :           XB + 16, XB2 + 16) > 1) return g; /*4+{1,2,3} = {5,6,7}*/
     902         323 :     if (ecm_elladd0(N, &g, 28, 4, /* 28 pts + (4 pts replicated 7fold) */
     903             :           XB + 28, XB2 + 28,
     904             :           XB,      XB2,
     905           0 :           XB + 32, XB2 + 32) > 1) return g;/*8+{1...7} = {9...15}*/
     906             :     /* ... and the remainder of the lot */
     907        1221 :     for (m = 5; m <= (ulong)gse; m++)
     908             :     { /* fill in from 2^(m-1)+1 to 2^m-1 in chunks of 64 and 60 points */
     909         898 :       ulong m2 = 2UL << m; /* will point at 2^(m-1)+1 */
     910        1977 :       for (j = 0; (ulong)j < m2-64; j+=64) /* executed 0 times when m = 5 */
     911             :       {
     912        1906 :         if (ecm_elladd0(N, &g, 64, 4,
     913        1079 :               XB + m2-4, XB2 + m2-4,
     914        1079 :               XB + j,    XB2 + j,
     915        1906 :               XB + m2+j, (m<(ulong)gse? XB2+m2+j: NULL)) > 1)
     916           0 :           return g;
     917             :       } /* j = m2-64 here, 60 points left */
     918        1221 :       if (ecm_elladd0(N, &g, 60, 4,
     919         898 :             XB + m2-4, XB2 + m2-4,
     920         898 :             XB + j,    XB2 + j,
     921        1221 :             XB + m2+j, (m<(ulong)gse? XB2+m2+j: NULL)) > 1)
     922           0 :         return g;
     923             :       /* when m=gse, drop Y coords of result, and when both equal 1024,
     924             :        * overwrite Y coords of second argument with X coords of result */
     925             :     }
     926         323 :     if (DEBUGLEVEL >= 7) err_printf("\t(baby step table complete)\n");
     927             :     /* initialize a few other things */
     928         323 :     bstp = bstp0; p = p0; np = np0; rcn = rcn0;
     929         323 :     g = gen_1; av1 = avma;
     930             :     /* scratchspace for prod (x_i-x_j) */
     931         323 :     avtmp = (pari_sp)new_chunk(8 * lgefint(N));
     932             :     /* The correct entry in XB to use depends on bstp and on where we are
     933             :      * on the helix. As we skip from prime to prime, bstp is incremented
     934             :      * by snextpr each time we wrap around through residue class number 0
     935             :      * (1 mod 210), but the baby step should not be taken until rcn>=rcn0,
     936             :      * i.e. until we pass again the residue class of p0.
     937             :      *
     938             :      * The correct signed multiplier is thus k = bstp - (rcn < rcn0),
     939             :      * and the offset from XB is four times (|k| - 1).  When k=0, we ignore
     940             :      * the current prime: if it had led to a factorization, this
     941             :      * would have been noted during the last giant step, or -- when we
     942             :      * first get here -- whilst initializing the helix.  When k > gss,
     943             :      * we must do a giant step and bump bstp back by -2*gss.
     944             :      *
     945             :      * The gcd of the product of X coord differences against N is taken just
     946             :      * before we do a giant step. */
     947     4515264 :     while (p < B2)
     948             :     {/* loop over probable primes p0 < p <= nextprime(B2), inserting giant
     949             :       * steps as necessary */
     950     4514948 :       p = snextpr(p, &np, &rcn, &bstp, uis2psp); /* next probable prime */
     951             :       /* work out the corresponding baby-step multiplier */
     952     4514948 :       k = bstp - (rcn < rcn0 ? 1 : 0);
     953     4514948 :       if (k > gss)
     954             :       { /* giant-step time, take gcd */
     955        1114 :         g = gcdii(g, N);
     956        1114 :         if (!is_pm1(g) && !equalii(g, N)) return g;
     957        1107 :         g = gen_1; set_avma(av1);
     958        2214 :         while (k > gss)
     959             :         { /* giant step */
     960        1107 :           if (DEBUGLEVEL >= 7) err_printf("\t(giant step at p = %lu)\n", p);
     961        1107 :           if (ecm_elladd0(N, &g, 64, 4, XG + i, YG + i,
     962           0 :                 Xh, Yh, Xh, Yh) > 1) return g;
     963        1107 :           if (ecm_elladd0(N, &g, 64, 4, XG + i, YG + i,
     964             :                 Xh + 64, Yh + 64, Xh + 64, Yh + 64) > 1)
     965           0 :             return g;
     966        1107 :           if (ecm_elladd0(N, &g, 64, 4, XG + i, YG + i,
     967             :                 Xh + 128, Yh + 128, Xh + 128, Yh + 128) > 1)
     968           0 :             return g;
     969        1107 :           bstp -= (gss << 1);
     970        1107 :           k = bstp - (rcn < rcn0? 1: 0); /* recompute multiplier */
     971             :         }
     972             :       }
     973     4514941 :       if (!k) continue; /* point of interest is already in Xh */
     974     4490276 :       if (k < 0) k = -k;
     975     4490276 :       m = ((ulong)k - 1) << 2;
     976             :       /* accumulate product of differences of X coordinates */
     977     4490276 :       j = rcn<<2;
     978     4490276 :       avma = avtmp; /* go to garbage zone; don't use set_avma */
     979     4490276 :       g = modii(mulii(g, subii(XB[m],   Xh[j])), N);
     980     4490276 :       g = modii(mulii(g, subii(XB[m+1], Xh[j+1])), N);
     981     4490276 :       g = modii(mulii(g, subii(XB[m+2], Xh[j+2])), N);
     982     4490276 :       g = mulii(g, subii(XB[m+3], Xh[j+3]));
     983     4490276 :       set_avma(av1);
     984     4490276 :       g = modii(g, N);
     985             :     }
     986         316 :     set_avma(av1);
     987             :   }
     988          68 :   return NULL;
     989             : }
     990             : 
     991             : /* ellfacteur() tuned to be useful as a first stage before MPQS, especially for
     992             :  * large arguments, when 'insist' is false, and now also for the case when
     993             :  * 'insist' is true, vaguely following suggestions by Paul Zimmermann
     994             :  * (http://www.loria.fr/~zimmerma/records/ecmnet.html). --GN 1998Jul,Aug */
     995             : static GEN
     996        3267 : ellfacteur(GEN N, int insist)
     997             : {
     998        3267 :   const long size = expi(N) + 1;
     999        3267 :   pari_sp av = avma;
    1000             :   struct ECM E;
    1001        3267 :   long nbc, dsn, dsnmax, rep = 0;
    1002        3267 :   if (insist)
    1003             :   {
    1004           0 :     const long DSNMAX = numberof(TB1)-1;
    1005           0 :     dsnmax = (size >> 2) - 10;
    1006           0 :     if (dsnmax < 0) dsnmax = 0;
    1007           0 :     else if (dsnmax > DSNMAX) dsnmax = DSNMAX;
    1008           0 :     E.seed = 1 + (nbcmax<<7)*(size&0xffff); /* seed for choice of curves */
    1009             : 
    1010           0 :     dsn = (size >> 3) - 5;
    1011           0 :     if (dsn < 0) dsn = 0; else if (dsn > 47) dsn = 47;
    1012             :     /* pick up the torch where noninsistent stage would have given up */
    1013           0 :     nbc = dsn + (dsn >> 2) + 9; /* 8 or more curves in parallel */
    1014           0 :     nbc &= ~3; /* 4 | nbc */
    1015             :   }
    1016             :   else
    1017             :   {
    1018        3267 :     dsn = (size - 140) >> 3;
    1019        3267 :     if (dsn < 0)
    1020             :     {
    1021             : #ifndef __EMX__ /* unless DOS/EMX: MPQS's disk access is abysmally slow */
    1022        3210 :       if (DEBUGLEVEL >= 4)
    1023           0 :         err_printf("ECM: number too small to justify this stage\n");
    1024        3210 :       return NULL; /* too small, decline the task */
    1025             : #endif
    1026             :       dsn = 0;
    1027          57 :     } else if (dsn > 12) dsn = 12;
    1028          57 :     rep = (size <= 248 ?
    1029          57 :            (size <= 176 ? (size - 124) >> 4 : (size - 148) >> 3) :
    1030          18 :            (size - 224) >> 1);
    1031             : #ifdef __EMX__ /* DOS/EMX: extra rounds (shun MPQS) */
    1032             :     rep += 20;
    1033             : #endif
    1034          57 :     dsnmax = 72;
    1035             :     /* Use disjoint sets of curves for non-insist and insist phases; moreover,
    1036             :      * repeated calls acting on factors of the same original number should try
    1037             :      * to use fresh curves. The following achieves this */
    1038          57 :     E.seed = 1 + (nbcmax<<3)*(size & 0xf);
    1039          57 :     nbc = -1;
    1040             :   }
    1041          57 :   ECM_init(&E, N, nbc);
    1042          57 :   if (DEBUGLEVEL >= 4)
    1043             :   {
    1044           0 :     timer_start(&E.T);
    1045           0 :     err_printf("ECM: working on %ld curves at a time; initializing", E.nbc);
    1046           0 :     if (!insist)
    1047             :     {
    1048           0 :       if (rep == 1) err_printf(" for one round");
    1049           0 :       else          err_printf(" for up to %ld rounds", rep);
    1050             :     }
    1051           0 :     err_printf("...\n");
    1052             :   }
    1053          57 :   if (dsn > dsnmax) dsn = dsnmax;
    1054             :   for(;;)
    1055          36 :   {
    1056          93 :     ulong B1 = insist? TB1[dsn]: TB1_for_stage[dsn];
    1057          93 :     GEN g = ECM_loop(&E, N, B1);
    1058          93 :     if (g)
    1059             :     {
    1060          25 :       if (DEBUGLEVEL >= 4)
    1061           0 :         err_printf("ECM: time = %6ld ms\n\tfound factor = %Ps\n",
    1062             :                    timer_delay(&E.T), g);
    1063          25 :       return gerepilecopy(av, g);
    1064             :     }
    1065          68 :     if (dsn < dsnmax)
    1066             :     {
    1067          68 :       if (insist) dsn++;
    1068          68 :       else { dsn += 2; if (dsn > dsnmax) dsn = dsnmax; }
    1069             :     }
    1070          68 :     if (!insist && !--rep)
    1071             :     {
    1072          32 :       if (DEBUGLEVEL >= 4)
    1073           0 :         err_printf("ECM: time = %6ld ms,\tellfacteur giving up.\n",
    1074             :                    timer_delay(&E.T));
    1075          32 :       return gc_NULL(av);
    1076             :     }
    1077             :   }
    1078             : }
    1079             : /* assume rounds >= 1, seed >= 1, B1 <= ULONG_MAX / 110 */
    1080             : GEN
    1081           0 : Z_ECM(GEN N, long rounds, long seed, ulong B1)
    1082             : {
    1083           0 :   pari_sp av = avma;
    1084             :   struct ECM E;
    1085             :   long i;
    1086           0 :   E.seed = seed;
    1087           0 :   ECM_init(&E, N, -1);
    1088           0 :   if (DEBUGLEVEL >= 4) timer_start(&E.T);
    1089           0 :   for (i = rounds; i--; )
    1090             :   {
    1091           0 :     GEN g = ECM_loop(&E, N, B1);
    1092           0 :     if (g) return gerepilecopy(av, g);
    1093             :   }
    1094           0 :   return gc_NULL(av);
    1095             : }
    1096             : 
    1097             : /***********************************************************************/
    1098             : /**                                                                   **/
    1099             : /**                FACTORIZATION (Pollard-Brent rho) --GN1998Jun18-26 **/
    1100             : /**  pollardbrent() returns a nontrivial factor of n, assuming n is   **/
    1101             : /**  composite and has no small prime divisor, or NULL if going on    **/
    1102             : /**  would take more time than we want to spend.  Sometimes it finds  **/
    1103             : /**  more than one factor, and returns a structure suitable for       **/
    1104             : /**  interpretation by ifac_crack. (Cf Algo 8.5.2 in ACiCNT)          **/
    1105             : /**                                                                   **/
    1106             : /***********************************************************************/
    1107             : #define VALUE(x) gel(x,0)
    1108             : #define EXPON(x) gel(x,1)
    1109             : #define CLASS(x) gel(x,2)
    1110             : 
    1111             : INLINE void
    1112       51817 : INIT(GEN x, GEN v, GEN e, GEN c) {
    1113       51817 :   VALUE(x) = v;
    1114       51817 :   EXPON(x) = e;
    1115       51817 :   CLASS(x) = c;
    1116       51817 : }
    1117             : static void
    1118       45877 : ifac_delete(GEN x) { INIT(x,NULL,NULL,NULL); }
    1119             : 
    1120             : static void
    1121           0 : rho_dbg(pari_timer *T, long c, long msg_mask)
    1122             : {
    1123           0 :   if (c & msg_mask) return;
    1124           0 :   err_printf("Rho: time = %6ld ms,\t%3ld round%s\n",
    1125             :              timer_delay(T), c, (c==1?"":"s"));
    1126             : }
    1127             : 
    1128             : static void
    1129    28246075 : one_iter(GEN *x, GEN *P, GEN x1, GEN n, long delta)
    1130             : {
    1131    28246075 :   *x = addis(remii(sqri(*x), n), delta);
    1132    28197686 :   *P = modii(mulii(*P, subii(x1, *x)), n);
    1133    28250691 : }
    1134             : /* Return NULL when we run out of time, or a single t_INT containing a
    1135             :  * nontrivial factor of n, or a vector of t_INTs, each triple of successive
    1136             :  * entries containing a factor, an exponent (equal to one),  and a factor
    1137             :  * class (NULL for unknown or zero for known composite),  matching the
    1138             :  * internal representation used by the ifac_*() routines below. Repeated
    1139             :  * factors may arise; the caller will sort the factors anyway. Result
    1140             :  * is not gerepile-able (contains NULL) */
    1141             : static GEN
    1142         802 : pollardbrent_i(GEN n, long size, long c0, long retries)
    1143             : {
    1144         802 :   long tf = lgefint(n), delta, msg_mask, c, k, k1, l;
    1145             :   pari_sp av;
    1146             :   GEN x, x1, y, P, g, g1, res;
    1147             :   pari_timer T;
    1148             : 
    1149         802 :   if (DEBUGLEVEL >= 4) timer_start(&T);
    1150         802 :   c = c0 << 5; /* 2^5 iterations per round */
    1151        1604 :   msg_mask = (size >= 448? 0x1fff:
    1152         802 :                            (size >= 192? (256L<<((size-128)>>6))-1: 0xff));
    1153         802 :   y = cgeti(tf);
    1154         802 :   x1= cgeti(tf);
    1155         802 :   av = avma;
    1156             : 
    1157         802 : PB_RETRY:
    1158             :  /* trick to make a 'random' choice determined by n.  Don't use x^2+0 or
    1159             :   * x^2-2, ever.  Don't use x^2-3 or x^2-7 with a starting value of 2.
    1160             :   * x^2+4, x^2+9 are affine conjugate to x^2+1, so don't use them either.
    1161             :   *
    1162             :   * (the point being that when we get called again on a composite cofactor
    1163             :   * of something we've already seen, we had better avoid the same delta) */
    1164         802 :   switch ((size + retries) & 7)
    1165             :   {
    1166         107 :     case 0:  delta=  1; break;
    1167         177 :     case 1:  delta= -1; break;
    1168          96 :     case 2:  delta=  3; break;
    1169          73 :     case 3:  delta=  5; break;
    1170          72 :     case 4:  delta= -5; break;
    1171          56 :     case 5:  delta=  7; break;
    1172         137 :     case 6:  delta= 11; break;
    1173             :     /* case 7: */
    1174          84 :     default: delta=-11; break;
    1175             :   }
    1176         802 :   if (DEBUGLEVEL >= 4)
    1177             :   {
    1178           0 :     if (!retries)
    1179           0 :       err_printf("Rho: searching small factor of %ld-bit integer\n", size);
    1180             :     else
    1181           0 :       err_printf("Rho: restarting for remaining rounds...\n");
    1182           0 :     err_printf("Rho: using X^2%+1ld for up to %ld rounds of 32 iterations\n",
    1183             :                delta, c >> 5);
    1184             :   }
    1185         802 :   x = gen_2; P = gen_1; g1 = NULL; k = 1; l = 1;
    1186         802 :   affui(2, y);
    1187         802 :   affui(2, x1);
    1188             :   for (;;) /* terminated under the control of c */
    1189             :   { /* use the polynomial  x^2 + delta */
    1190    13205185 :     one_iter(&x, &P, x1, n, delta);
    1191             : 
    1192    13206750 :     if ((--c & 0x1f)==0)
    1193             :     { /* one round complete */
    1194      412593 :       g = gcdii(n, P); if (!is_pm1(g)) goto fin;
    1195      412326 :       if (c <= 0)
    1196             :       { /* getting bored */
    1197         396 :         if (DEBUGLEVEL >= 4)
    1198           0 :           err_printf("Rho: time = %6ld ms,\tPollard-Brent giving up.\n",
    1199             :                      timer_delay(&T));
    1200         396 :         return NULL;
    1201             :       }
    1202      411930 :       P = gen_1;
    1203      411930 :       if (DEBUGLEVEL >= 4) rho_dbg(&T, c0-(c>>5), msg_mask);
    1204      411930 :       affii(x,y); x = y; set_avma(av);
    1205             :     }
    1206             : 
    1207    13204676 :     if (--k) continue; /* normal end of loop body */
    1208             : 
    1209        8612 :     if (c & 0x1f) /* otherwise, we already checked */
    1210             :     {
    1211        4812 :       g = gcdii(n, P); if (!is_pm1(g)) goto fin;
    1212        4812 :       P = gen_1;
    1213             :     }
    1214             : 
    1215             :    /* Fast forward phase, doing l inner iterations without computing gcds.
    1216             :     * Check first whether it would take us beyond the alloted time.
    1217             :     * Fast forward rounds count only half (although they're taking
    1218             :     * more like 2/3 the time of normal rounds).  This to counteract the
    1219             :     * nuisance that all c0 between 4096 and 6144 would act exactly as
    1220             :     * 4096;  with the halving trick only the range 4096..5120 collapses
    1221             :     * (similarly for all other powers of two) */
    1222        8612 :     if ((c -= (l>>1)) <= 0)
    1223             :     { /* got bored */
    1224         179 :       if (DEBUGLEVEL >= 4)
    1225           0 :         err_printf("Rho: time = %6ld ms,\tPollard-Brent giving up.\n",
    1226             :                    timer_delay(&T));
    1227         179 :       return NULL;
    1228             :     }
    1229        8433 :     c &= ~0x1f; /* keep it on multiples of 32 */
    1230             : 
    1231             :     /* Fast forward loop */
    1232        8433 :     affii(x, x1); set_avma(av); x = x1;
    1233        8434 :     k = l; l <<= 1;
    1234             :     /* don't show this for the first several (short) fast forward phases. */
    1235        8434 :     if (DEBUGLEVEL >= 4 && (l>>7) > msg_mask)
    1236           0 :       err_printf("Rho: fast forward phase (%ld rounds of 64)...\n", l>>7);
    1237    15080029 :     for (k1=k; k1; k1--)
    1238             :     {
    1239    15072194 :       one_iter(&x, &P, x1, n, delta);
    1240    15069711 :       if ((k1 & 0x1f) == 0) gerepileall(av, 2, &x, &P);
    1241             :     }
    1242        7835 :     if (DEBUGLEVEL >= 4 && (l>>7) > msg_mask)
    1243           0 :       err_printf("Rho: time = %6ld ms,\t%3ld rounds, back to normal mode\n",
    1244           0 :                  timer_delay(&T), c0-(c>>5));
    1245        7835 :     affii(x,y); P = gerepileuptoint(av, P); x = y;
    1246             :   } /* forever */
    1247             : 
    1248         227 : fin:
    1249             :   /* An accumulated gcd was > 1 */
    1250         227 :   if  (!equalii(g,n))
    1251             :   { /* if it isn't n, and looks prime, return it */
    1252         227 :     if (MR_Jaeschke(g))
    1253             :     {
    1254         227 :       if (DEBUGLEVEL >= 4)
    1255             :       {
    1256           0 :         rho_dbg(&T, c0-(c>>5), 0);
    1257           0 :         err_printf("\tfound factor = %Ps\n",g);
    1258             :       }
    1259         227 :       return g;
    1260             :     }
    1261           0 :     set_avma(av); g1 = icopy(g);  /* known composite, keep it safe */
    1262           0 :     av = avma;
    1263             :   }
    1264           0 :   else g1 = n; /* and work modulo g1 for backtracking */
    1265             : 
    1266             :   /* Here g1 is known composite */
    1267           0 :   if (DEBUGLEVEL >= 4 && size > 192)
    1268           0 :     err_printf("Rho: hang on a second, we got something here...\n");
    1269           0 :   x = y;
    1270             :   for(;;)
    1271             :   { /* backtrack until period recovered. Must terminate */
    1272           0 :     x = addis(remii(sqri(x), g1), delta);
    1273           0 :     g = gcdii(subii(x1, x), g1); if (!is_pm1(g)) break;
    1274             : 
    1275           0 :     if (DEBUGLEVEL >= 4 && (--c & 0x1f) == 0) rho_dbg(&T, c0-(c>>5), msg_mask);
    1276             :   }
    1277             : 
    1278           0 :   if (g1 == n || equalii(g,g1))
    1279             :   {
    1280           0 :     if (g1 == n && equalii(g,g1))
    1281             :     { /* out of luck */
    1282           0 :       if (DEBUGLEVEL >= 4)
    1283             :       {
    1284           0 :         rho_dbg(&T, c0-(c>>5), 0);
    1285           0 :         err_printf("\tPollard-Brent failed.\n");
    1286             :       }
    1287           0 :       if (++retries >= 4) pari_err_BUG("");
    1288           0 :       goto PB_RETRY;
    1289             :     }
    1290             :     /* half lucky: we've split n, but g1 equals either g or n */
    1291           0 :     if (DEBUGLEVEL >= 4)
    1292             :     {
    1293           0 :       rho_dbg(&T, c0-(c>>5), 0);
    1294           0 :       err_printf("\tfound %sfactor = %Ps\n", (g1!=n ? "composite " : ""), g);
    1295             :     }
    1296           0 :     res = cgetg(7, t_VEC);
    1297             :     /* g^1: known composite when g1!=n */
    1298           0 :     INIT(res+1, g, gen_1, (g1!=n? gen_0: NULL));
    1299             :     /* cofactor^1: status unknown */
    1300           0 :     INIT(res+4, diviiexact(n,g), gen_1, NULL);
    1301           0 :     return res;
    1302             :   }
    1303             :   /* g < g1 < n : our lucky day -- we've split g1, too */
    1304           0 :   res = cgetg(10, t_VEC);
    1305             :   /* unknown status for all three factors */
    1306           0 :   INIT(res+1, g,                gen_1, NULL);
    1307           0 :   INIT(res+4, diviiexact(g1,g), gen_1, NULL);
    1308           0 :   INIT(res+7, diviiexact(n,g1), gen_1, NULL);
    1309           0 :   if (DEBUGLEVEL >= 4)
    1310             :   {
    1311           0 :     rho_dbg(&T, c0-(c>>5), 0);
    1312           0 :     err_printf("\tfound factors = %Ps, %Ps,\n\tand %Ps\n",
    1313           0 :                gel(res,1), gel(res,4), gel(res,7));
    1314             :   }
    1315           0 :   return res;
    1316             : }
    1317             : /* decline if n < 2^96 */
    1318             : static GEN
    1319        3494 : pollardbrent(GEN n)
    1320             : {
    1321        3494 :   const long tune = 14; /* FIXME: retune this */
    1322        3494 :   long c0, size, tf = lgefint(n);
    1323             : #ifdef LONG_IS_64BIT
    1324        3045 :   if (tf < 4 || (tf == 4 && uel(n,2) < (1UL << 32))) return NULL;
    1325             : #else  /* 32 bits */
    1326         449 :   if (tf < 5) return NULL;
    1327             : #endif
    1328         802 :   size = expi(n) + 1;
    1329             :   /* nonlinear increase in effort, kicking in around 80 bits */
    1330         802 :   if (size <= 301) /* 301 gives 48121 + tune */
    1331         795 :     c0 = tune + size-60 + ((size-73)>>1)*((size-70)>>3)*((size-56)>>4);
    1332             :   else
    1333           7 :     c0 = 49152; /* ECM is faster when it'd take longer */
    1334         802 :   return pollardbrent_i(n, size, c0, 0);
    1335             : }
    1336             : GEN
    1337           0 : Z_pollardbrent(GEN n, long rounds, long seed)
    1338             : {
    1339           0 :   pari_sp av = avma;
    1340           0 :   GEN v = pollardbrent_i(n, expi(n)+1, rounds, seed);
    1341           0 :   if (!v) return NULL;
    1342           0 :   if (typ(v) == t_INT) v = mkvec2(v, diviiexact(n,v));
    1343           0 :   else if (lg(v) == 7) v = mkvec2(gel(v,1), gel(v,4));
    1344           0 :   else v = mkvec3(gel(v,1), gel(v,4), gel(v,7));
    1345           0 :   return gerepilecopy(av, v);
    1346             : }
    1347             : 
    1348             : /***********************************************************************/
    1349             : /**              FACTORIZATION (Shanks' SQUFOF) --GN2000Sep30-Oct01   **/
    1350             : /**  squfof() returns a proper factor of n, or NULL (failure). Assume **/
    1351             : /**  n is composite, not a square, and has no small prime divisors.   **/
    1352             : /**  Works on two discriminants at once: n and 5n or 3n and 4n        **/
    1353             : /**  Present implementation is limited to input <2^59, and works most **/
    1354             : /**  of the time in signed arithmetic on integers <2^31 in absolute   **/
    1355             : /**  size. (Cf. Algo 8.7.2 in ACiCNT)                                 **/
    1356             : /***********************************************************************/
    1357             : 
    1358             : /* squfof_ambig walks back along the ambiguous cycle until we hit an ambiguous
    1359             :  * form and thus the desired factor, which it returns. Returs 0 on failure.
    1360             :  *
    1361             :  * Input: a form (A, B, -C) with A = a^2, where a isn't blacklisted and
    1362             :  * (a, B) = 1. We should now proceed reducing the form (a, -B, -aC), but the
    1363             :  * first reduction step always sends this to (-aC, B, a), and the next one,
    1364             :  * with q computed as usual from B and a (occupying the c position), gives a
    1365             :  * reduced form, whose third member is easiest to recover by going back to D.
    1366             :  * From this point onwards, we're once again working with single-word numbers.
    1367             :  * No need to track signs, just work with the abs values of the coefficients.
    1368             :  * HACK: if LONG_IS_64BIT, D is actually a typecast long */
    1369             : static long
    1370        2091 : squfof_ambig(long a, long B, long dd, GEN D)
    1371             : {
    1372             :   long b, c, q, qa, a0, b0, b1;
    1373        2091 :   long cnt = 0; /* count reduction steps on the cycle */
    1374             : 
    1375        2091 :   q = (dd + (B>>1)) / a;
    1376        2091 :   qa = q * a;
    1377        2091 :   b = (qa - B) + qa; /* avoid overflow */
    1378             : #ifdef LONG_IS_64BIT
    1379        1620 :   c = (((long)D - b*b) >> 2) / a;
    1380             : #else
    1381             :   {
    1382         471 :     pari_sp av = avma;
    1383         471 :     c = itos(divis(shifti(subii(D, sqrs(b)), -2), a));
    1384         471 :     set_avma(av);
    1385             :   }
    1386             : #endif
    1387        2091 :   a0 = a; b0 = b1 = b; /* end of loop detection and safeguard */
    1388             :   for (;;)
    1389      957157 :   { /* reduction step */
    1390      959248 :     long c0 = c, qc, qcb;
    1391      959248 :     if (c0 > dd)
    1392      267802 :       q = 1;
    1393             :     else
    1394      691446 :       q = (dd + (b>>1)) / c0;
    1395      959248 :     if (q == 1)
    1396             :     {
    1397      396866 :       qcb = c0 - b; b = c0 + qcb; c = a - qcb;
    1398             :     }
    1399             :     else
    1400             :     {
    1401      562382 :       qc = q*c0; qcb = qc - b; b = qc + qcb; c = a - q*qcb;
    1402             :     }
    1403      959248 :     a = c0;
    1404             : 
    1405      959248 :     cnt++; if (b == b1) break;
    1406             : 
    1407             :     /* safeguard against infinite loop: we walked the cycle in vain.
    1408             :      * (I don't think this can actually happen.) */
    1409      957157 :     if (b == b0 && a == a0) return 0;
    1410             : 
    1411      957157 :     b1 = b;
    1412             :   }
    1413        2091 :   q = a&1 ? a : a>>1;
    1414        2091 :   if (DEBUGLEVEL >= 4)
    1415             :   {
    1416           0 :     if (q > 1)
    1417           0 :       err_printf("SQUFOF: found factor %ld from ambiguous form\n"
    1418             :                  "\tafter %ld steps on the ambiguous cycle\n",
    1419           0 :                  q / ugcd(q,15), cnt);
    1420             :     else
    1421           0 :       err_printf("SQUFOF: ...found nothing on the ambiguous cycle\n"
    1422             :                  "\tafter %ld steps there\n", cnt);
    1423           0 :     if (DEBUGLEVEL >= 6) err_printf("SQUFOF: squfof_ambig returned %ld\n", q);
    1424             :   }
    1425        2091 :   return q;
    1426             : }
    1427             : 
    1428             : #define SQUFOF_BLACKLIST_SZ 64
    1429             : 
    1430             : /* assume (n,30) = 1 */
    1431             : static GEN
    1432        5125 : squfof(GEN n)
    1433             : {
    1434             :   ulong d1, d2;
    1435             : #ifdef LONG_IS_64BIT
    1436             :   ulong uD1, uD2;
    1437             : #endif
    1438        5125 :   long tf = lgefint(n), nm4, cnt = 0;
    1439             :   long a1, b1, c1, dd1, L1, a2, b2, c2, dd2, L2, a, q, c, qc, qcb;
    1440             :   GEN D1, D2;
    1441        5125 :   pari_sp av = avma;
    1442             :   long blacklist1[SQUFOF_BLACKLIST_SZ], blacklist2[SQUFOF_BLACKLIST_SZ];
    1443        5125 :   long blp1 = 0, blp2 = 0;
    1444        5125 :   int act1 = 1, act2 = 1;
    1445             : 
    1446             : #ifdef LONG_IS_64BIT
    1447        4323 :   if (tf > 3 || (tf == 3 && uel(n,2)             >= (1UL << 46)))
    1448             : #else  /* 32 bits */
    1449         802 :   if (tf > 4 || (tf == 4 && (ulong)(*int_MSW(n)) >= (1UL << 17))) /* 49 */
    1450             : #endif
    1451        3494 :     return NULL; /* n too large */
    1452             : 
    1453             :   /* now we have 5 < n < 2^59 */
    1454        1631 :   nm4 = mod4(n);
    1455             : #ifdef LONG_IS_64BIT
    1456        1278 :   if (nm4 == 1)
    1457             :   { /* n = 1 (mod4):  run one iteration on D1 = n, another on D2 = 5n */
    1458         600 :     uD1 = n[2];
    1459         600 :     uD2 = 5 * n[2]; d2 = usqrt(uD2); dd2 = (long)((d2>>1) + (d2&1));
    1460         600 :     b2 = (long)((d2-1) | 1); /* b1, b2 will always stay odd */
    1461             :   }
    1462             :   else
    1463             :   { /* n = 3 (mod4):  run one iteration on D1 = 3n, another on D2 = 4n */
    1464         678 :     uD1 = 3 * n[2];
    1465         678 :     uD2 = 4 * n[2]; dd2 = usqrt(n[2]); d2 =  dd2 << 1;
    1466         678 :     b2 = (long)(d2 & (~1UL)); /* largest even below d2, will stay even */
    1467             :   }
    1468        1278 :   D1 = (GEN)uD1;
    1469        1278 :   D2 = (GEN)uD2;
    1470        1278 :   d1 = usqrt(uD1);
    1471        1278 :   b1 = (long)((d1-1) | 1); /* largest odd number not exceeding d1 */
    1472             :   /* c1 != 0 else n or 3n would be a square */
    1473        1278 :   c1 = (uD1 - b1*b1) / 4;
    1474             :   /* c2 != 0 else 5n would be a square */
    1475        1278 :   c2 = (uD2 - b2*b2) / 4;
    1476             : #else
    1477         353 :   if (nm4 == 1)
    1478             :   { /* n = 1 (mod4):  run one iteration on D1 = n, another on D2 = 5n */
    1479         179 :     D1 = n;
    1480         179 :     D2 = mului(5,n); d2 = itou(sqrti(D2)); dd2 = (long)((d2>>1) + (d2&1));
    1481         179 :     b2 = (long)((d2-1) | 1); /* b1, b2 will always stay odd */
    1482             :   }
    1483             :   else
    1484             :   { /* n = 3 (mod4):  run one iteration on D1 = 3n, another on D2 = 4n */
    1485         174 :     D1 = mului(3,n);
    1486         174 :     D2 = shifti(n,2); dd2 = itou(sqrti(n)); d2 =  dd2 << 1;
    1487         174 :     b2 = (long)(d2 & (~1UL)); /* largest even below d2, will stay even */
    1488             :   }
    1489         353 :   d1 = itou(sqrti(D1));
    1490         353 :   b1 = (long)((d1-1) | 1); /* largest odd number not exceeding d1 */
    1491             :   /* c1 != 0 else n or 3n would be a square */
    1492         353 :   c1 = itos(shifti(subii(D1, sqru((ulong)b1)), -2));
    1493             :   /* c2 != 0 else 5n would be a square */
    1494         353 :   c2 = itos(shifti(subii(D2, sqru((ulong)b2)), -2));
    1495             : #endif
    1496        1631 :   L1 = (long)usqrt(d1);
    1497        1631 :   L2 = (long)usqrt(d2);
    1498             :   /* dd1 used to compute floor((d1+b1)/2) as dd1+floor(b1/2), without
    1499             :    * overflowing the 31bit signed integer size limit. Same for dd2. */
    1500        1631 :   dd1 = (long) ((d1>>1) + (d1&1));
    1501        1631 :   a1 = a2 = 1;
    1502             : 
    1503             :   /* The two (identity) forms (a1,b1,-c1) and (a2,b2,-c2) are now set up.
    1504             :    *
    1505             :    * a1 and c1 represent the absolute values of the a,c coefficients; we keep
    1506             :    * track of the sign separately, via the iteration counter cnt: when cnt is
    1507             :    * even, c < 0 and a > 0, else c > 0 and a < 0.
    1508             :    *
    1509             :    * L1, L2 are the limits for blacklisting small leading coefficients
    1510             :    * on the principal cycle, to guarantee that when we find a square form,
    1511             :    * its square root will belong to an ambiguous cycle, i.e. won't be an
    1512             :    * earlier form on the principal cycle.
    1513             :    *
    1514             :    * When n = 3(mod 4), D2 = 12(mod 16), and b^2 is 0 or 4 mod 16.
    1515             :    * It follows that 4*a*c must be 4 or 8 mod 16, respectively, so at most
    1516             :    * one of a,c can be divisible by 2 at most to the first power. This fact
    1517             :    * is used a couple of times below.
    1518             :    *
    1519             :    * The flags act1, act2 remain true while the respective cycle is still
    1520             :    * active; we drop them to false when we return to the identity form with-
    1521             :    * out having found a square form (or when the blacklist overflows, which
    1522             :    * shouldn't happen). */
    1523        1631 :   if (DEBUGLEVEL >= 4)
    1524           0 :     err_printf("SQUFOF: entering main loop with forms\n"
    1525             :                "\t(1, %ld, %ld) and (1, %ld, %ld)\n", b1, -c1, b2, -c2);
    1526             : 
    1527             :   /* MAIN LOOP: walk around the principal cycle looking for a square form.
    1528             :    * Blacklist small leading coefficients.
    1529             :    *
    1530             :    * The reduction operator can be computed entirely in 32-bit arithmetic:
    1531             :    * Let q = floor(floor((d1+b1)/2)/c1)  (when c1>dd1, q=1, which happens
    1532             :    * often enough to special-case it). Then the new b1 = (q*c1-b1) + q*c1,
    1533             :    * which does not overflow, and the new c1 = a1 - q*(q*c1-b1), which is
    1534             :    * bounded by d1 in abs size since both the old and the new a1 are positive
    1535             :    * and bounded by d1. */
    1536     1384926 :   while (act1 || act2)
    1537             :   {
    1538     1384926 :     if (act1)
    1539             :     { /* send first form through reduction operator if active */
    1540     1384926 :       c = c1;
    1541     1384926 :       q = (c > dd1)? 1: (dd1 + (b1>>1)) / c;
    1542     1384926 :       if (q == 1)
    1543      574343 :       { qcb = c - b1; b1 = c + qcb; c1 = a1 - qcb; }
    1544             :       else
    1545      810583 :       { qc = q*c; qcb = qc - b1; b1 = qc + qcb; c1 = a1 - q*qcb; }
    1546     1384926 :       a1 = c;
    1547             : 
    1548     1384926 :       if (a1 <= L1)
    1549             :       { /* blacklist this */
    1550        1080 :         if (blp1 >= SQUFOF_BLACKLIST_SZ) /* overflows: shouldn't happen */
    1551           0 :           act1 = 0;
    1552             :         else
    1553             :         {
    1554        1080 :           if (DEBUGLEVEL >= 6)
    1555           0 :             err_printf("SQUFOF: blacklisting a = %ld on first cycle\n", a1);
    1556        1080 :           blacklist1[blp1++] = a1;
    1557             :         }
    1558             :       }
    1559             :     }
    1560     1384926 :     if (act2)
    1561             :     { /* send second form through reduction operator if active */
    1562     1384734 :       c = c2;
    1563     1384734 :       q = (c > dd2)? 1: (dd2 + (b2>>1)) / c;
    1564     1384734 :       if (q == 1)
    1565      574193 :       { qcb = c - b2; b2 = c + qcb; c2 = a2 - qcb; }
    1566             :       else
    1567      810541 :       { qc = q*c; qcb = qc - b2; b2 = qc + qcb; c2 = a2 - q*qcb; }
    1568     1384734 :       a2 = c;
    1569             : 
    1570     1384734 :       if (a2 <= L2)
    1571             :       { /* blacklist this */
    1572        1108 :         if (blp2 >= SQUFOF_BLACKLIST_SZ) /* overflows: shouldn't happen */
    1573           0 :           act2 = 0;
    1574             :         else
    1575             :         {
    1576        1108 :           if (DEBUGLEVEL >= 6)
    1577           0 :             err_printf("SQUFOF: blacklisting a = %ld on second cycle\n", a2);
    1578        1108 :           blacklist2[blp2++] = a2;
    1579             :         }
    1580             :       }
    1581             :     }
    1582             : 
    1583     1384926 :     if (++cnt & 1) continue; /* odd iteration */
    1584             :     /* even iteration: the leading coefficients are positive */
    1585             : 
    1586             :     /* examine first form if active */
    1587      692463 :     if (act1 && a1 == 1) /* back to identity */
    1588             :     { /* drop this discriminant */
    1589           0 :       act1 = 0;
    1590           0 :       if (DEBUGLEVEL >= 4)
    1591           0 :         err_printf("SQUFOF: first cycle exhausted after %ld iterations,\n"
    1592             :                    "\tdropping it\n", cnt);
    1593             :     }
    1594      692463 :     if (act1)
    1595             :     {
    1596      692463 :       if (uissquareall((ulong)a1, (ulong*)&a))
    1597             :       { /* square form */
    1598        1309 :         if (DEBUGLEVEL >= 4)
    1599           0 :           err_printf("SQUFOF: square form (%ld^2, %ld, %ld) on first cycle\n"
    1600             :                      "\tafter %ld iterations\n", a, b1, -c1, cnt);
    1601        1309 :         if (a <= L1)
    1602             :         { /* blacklisted? */
    1603             :           long j;
    1604        2436 :           for (j = 0; j < blp1; j++)
    1605        1613 :             if (a == blacklist1[j]) { a = 0; break; }
    1606             :         }
    1607        1309 :         if (a > 0)
    1608             :         { /* not blacklisted */
    1609         823 :           q = ugcd(a, b1); /* imprimitive form? */
    1610         823 :           if (q > 1)
    1611             :           { /* q^2 divides D1 hence n [ assuming n % 3 != 0 ] */
    1612           0 :             set_avma(av);
    1613           0 :             if (DEBUGLEVEL >= 4) err_printf("SQUFOF: found factor %ld^2\n", q);
    1614           0 :             return mkvec3(utoipos(q), gen_2, NULL);/* q^2, unknown status */
    1615             :           }
    1616             :           /* chase the inverse root form back along the ambiguous cycle */
    1617         823 :           q = squfof_ambig(a, b1, dd1, D1);
    1618         823 :           if (q > 3)
    1619             :           {
    1620         677 :             if (nm4 == 3 && q % 3 == 0) q /= 3;
    1621         677 :             return gc_utoipos(av, q); /* SUCCESS! */
    1622             :           }
    1623             :         }
    1624         486 :         else if (DEBUGLEVEL >= 4) /* blacklisted */
    1625           0 :           err_printf("SQUFOF: ...but the root form seems to be on the "
    1626             :                      "principal cycle\n");
    1627             :       }
    1628             :     }
    1629             : 
    1630             :     /* examine second form if active */
    1631      691786 :     if (act2 && a2 == 1) /* back to identity form */
    1632             :     { /* drop this discriminant */
    1633           2 :       act2 = 0;
    1634           2 :       if (DEBUGLEVEL >= 4)
    1635           0 :         err_printf("SQUFOF: second cycle exhausted after %ld iterations,\n"
    1636             :                    "\tdropping it\n", cnt);
    1637             :     }
    1638      691786 :     if (act2)
    1639             :     {
    1640      691690 :       if (uissquareall((ulong)a2, (ulong*)&a))
    1641             :       { /* square form */
    1642        1549 :         if (DEBUGLEVEL >= 4)
    1643           0 :           err_printf("SQUFOF: square form (%ld^2, %ld, %ld) on second cycle\n"
    1644             :                      "\tafter %ld iterations\n", a, b2, -c2, cnt);
    1645        1549 :         if (a <= L2)
    1646             :         { /* blacklisted? */
    1647             :           long j;
    1648        2729 :           for (j = 0; j < blp2; j++)
    1649        1461 :             if (a == blacklist2[j]) { a = 0; break; }
    1650             :         }
    1651        1549 :         if (a > 0)
    1652             :         { /* not blacklisted */
    1653        1268 :           q = ugcd(a, b2); /* imprimitive form? */
    1654             :           /* NB if b2 is even, a is odd, so the gcd is always odd */
    1655        1268 :           if (q > 1)
    1656             :           { /* q^2 divides D2 hence n [ assuming n % 5 != 0 ] */
    1657           0 :             set_avma(av);
    1658           0 :             if (DEBUGLEVEL >= 4) err_printf("SQUFOF: found factor %ld^2\n", q);
    1659           0 :             return mkvec3(utoipos(q), gen_2, NULL);/* q^2, unknown status */
    1660             :           }
    1661             :           /* chase the inverse root form along the ambiguous cycle */
    1662        1268 :           q = squfof_ambig(a, b2, dd2, D2);
    1663        1268 :           if (q > 5)
    1664             :           {
    1665         954 :             if (nm4 == 1 && q % 5 == 0) q /= 5;
    1666         954 :             return gc_utoipos(av, q); /* SUCCESS! */
    1667             :           }
    1668             :         }
    1669         281 :         else if (DEBUGLEVEL >= 4)        /* blacklisted */
    1670           0 :           err_printf("SQUFOF: ...but the root form seems to be on the "
    1671             :                      "principal cycle\n");
    1672             :       }
    1673             :     }
    1674             :   } /* end main loop */
    1675             : 
    1676           0 :   if (DEBUGLEVEL>=4) err_printf("SQUFOF: giving up\n");
    1677           0 :   return gc_NULL(av);
    1678             : }
    1679             : 
    1680             : /***********************************************************************/
    1681             : /*                    DETECTING ODD POWERS  --GN1998Jun28              */
    1682             : /*   Factoring engines like MPQS which ultimately rely on computing    */
    1683             : /*   gcd(N, x^2-y^2) to find a nontrivial factor of N can't split      */
    1684             : /*   N = p^k for an odd prime p, since (Z/p^k)^* is then cyclic. Here  */
    1685             : /*   is an analogue of Z_issquareall() for 3rd, 5th and 7th powers.    */
    1686             : /*   The general case is handled by is_kth_power                       */
    1687             : /***********************************************************************/
    1688             : 
    1689             : /* Multistage sieve. First stages work mod 211, 209, 61, 203 in this order
    1690             :  * (first reduce mod the product of these and then take the remainder apart).
    1691             :  * Second stages use 117, 31, 43, 71. Moduli which are no longer interesting
    1692             :  * are skipped. Everything is encoded in a table of 106 24-bit masks. We only
    1693             :  * need the first half of the residues.  Three bits per modulus indicate which
    1694             :  * residues are 7th (bit 2), 5th (bit 1) or 3rd (bit 0) powers; the eight
    1695             :  * moduli above are assigned right-to-left. The table was generated using: */
    1696             : 
    1697             : #if 0
    1698             : L = [71, 43, 31, [O(3^2),O(13)], [O(7),O(29)], 61, [O(11),O(19)], 211];
    1699             : ispow(x, N, k)=
    1700             : {
    1701             :   if (type(N) == "t_INT", return (ispower(Mod(x,N), k)));
    1702             :   for (i = 1, #N, if (!ispower(x + N[i], k), return (0))); 1
    1703             : }
    1704             : check(r) =
    1705             : {
    1706             :   print1("  0");
    1707             :   for (i=1,#L,
    1708             :     N = 0;
    1709             :     if (ispow(r, L[i], 3), N += 1);
    1710             :     if (ispow(r, L[i], 5), N += 2);
    1711             :     if (ispow(r, L[i], 7), N += 4);
    1712             :     print1(N);
    1713             :   ); print("ul,  /* ", r, " */")
    1714             : }
    1715             : for (r = 0, 105, check(r))
    1716             : #endif
    1717             : static ulong powersmod[106] = {
    1718             :   077777777ul,  /* 0 */
    1719             :   077777777ul,  /* 1 */
    1720             :   013562440ul,  /* 2 */
    1721             :   012402540ul,  /* 3 */
    1722             :   013562440ul,  /* 4 */
    1723             :   052662441ul,  /* 5 */
    1724             :   016603440ul,  /* 6 */
    1725             :   016463450ul,  /* 7 */
    1726             :   013573551ul,  /* 8 */
    1727             :   012462540ul,  /* 9 */
    1728             :   012462464ul,  /* 10 */
    1729             :   013462771ul,  /* 11 */
    1730             :   012406473ul,  /* 12 */
    1731             :   012463641ul,  /* 13 */
    1732             :   052463646ul,  /* 14 */
    1733             :   012503446ul,  /* 15 */
    1734             :   013562440ul,  /* 16 */
    1735             :   052466440ul,  /* 17 */
    1736             :   012472451ul,  /* 18 */
    1737             :   012462454ul,  /* 19 */
    1738             :   032463550ul,  /* 20 */
    1739             :   013403664ul,  /* 21 */
    1740             :   013463460ul,  /* 22 */
    1741             :   032562565ul,  /* 23 */
    1742             :   012402540ul,  /* 24 */
    1743             :   052662441ul,  /* 25 */
    1744             :   032672452ul,  /* 26 */
    1745             :   013573551ul,  /* 27 */
    1746             :   012467541ul,  /* 28 */
    1747             :   012567640ul,  /* 29 */
    1748             :   032706450ul,  /* 30 */
    1749             :   012762452ul,  /* 31 */
    1750             :   033762662ul,  /* 32 */
    1751             :   012502562ul,  /* 33 */
    1752             :   032463562ul,  /* 34 */
    1753             :   013563440ul,  /* 35 */
    1754             :   016663440ul,  /* 36 */
    1755             :   036662550ul,  /* 37 */
    1756             :   012462552ul,  /* 38 */
    1757             :   033502450ul,  /* 39 */
    1758             :   012462643ul,  /* 40 */
    1759             :   033467540ul,  /* 41 */
    1760             :   017403441ul,  /* 42 */
    1761             :   017463462ul,  /* 43 */
    1762             :   017472460ul,  /* 44 */
    1763             :   033462470ul,  /* 45 */
    1764             :   052566450ul,  /* 46 */
    1765             :   013562640ul,  /* 47 */
    1766             :   032403640ul,  /* 48 */
    1767             :   016463450ul,  /* 49 */
    1768             :   016463752ul,  /* 50 */
    1769             :   033402440ul,  /* 51 */
    1770             :   012462540ul,  /* 52 */
    1771             :   012472540ul,  /* 53 */
    1772             :   053562462ul,  /* 54 */
    1773             :   012463465ul,  /* 55 */
    1774             :   012663470ul,  /* 56 */
    1775             :   052607450ul,  /* 57 */
    1776             :   012566553ul,  /* 58 */
    1777             :   013466440ul,  /* 59 */
    1778             :   012502741ul,  /* 60 */
    1779             :   012762744ul,  /* 61 */
    1780             :   012763740ul,  /* 62 */
    1781             :   012763443ul,  /* 63 */
    1782             :   013573551ul,  /* 64 */
    1783             :   013462471ul,  /* 65 */
    1784             :   052502460ul,  /* 66 */
    1785             :   012662463ul,  /* 67 */
    1786             :   012662451ul,  /* 68 */
    1787             :   012403550ul,  /* 69 */
    1788             :   073567540ul,  /* 70 */
    1789             :   072463445ul,  /* 71 */
    1790             :   072462740ul,  /* 72 */
    1791             :   012472442ul,  /* 73 */
    1792             :   012462644ul,  /* 74 */
    1793             :   013406650ul,  /* 75 */
    1794             :   052463471ul,  /* 76 */
    1795             :   012563474ul,  /* 77 */
    1796             :   013503460ul,  /* 78 */
    1797             :   016462441ul,  /* 79 */
    1798             :   016462440ul,  /* 80 */
    1799             :   012462540ul,  /* 81 */
    1800             :   013462641ul,  /* 82 */
    1801             :   012463454ul,  /* 83 */
    1802             :   013403550ul,  /* 84 */
    1803             :   057563540ul,  /* 85 */
    1804             :   017466441ul,  /* 86 */
    1805             :   017606471ul,  /* 87 */
    1806             :   053666573ul,  /* 88 */
    1807             :   012562561ul,  /* 89 */
    1808             :   013473641ul,  /* 90 */
    1809             :   032573440ul,  /* 91 */
    1810             :   016763440ul,  /* 92 */
    1811             :   016702640ul,  /* 93 */
    1812             :   033762552ul,  /* 94 */
    1813             :   012562550ul,  /* 95 */
    1814             :   052402451ul,  /* 96 */
    1815             :   033563441ul,  /* 97 */
    1816             :   012663561ul,  /* 98 */
    1817             :   012677560ul,  /* 99 */
    1818             :   012462464ul,  /* 100 */
    1819             :   032562642ul,  /* 101 */
    1820             :   013402551ul,  /* 102 */
    1821             :   032462450ul,  /* 103 */
    1822             :   012467445ul,  /* 104 */
    1823             :   032403440ul,  /* 105 */
    1824             : };
    1825             : 
    1826             : static int
    1827     3928662 : check_res(ulong x, ulong N, int shift, ulong *mask)
    1828             : {
    1829     3928662 :   long r = x%N; if ((ulong)r> (N>>1)) r = N - r;
    1830     3928662 :   *mask &= (powersmod[r] >> shift);
    1831     3928662 :   return *mask;
    1832             : }
    1833             : 
    1834             : /* is x mod 211*209*61*203*117*31*43*71 a 3rd, 5th or 7th power ? */
    1835             : int
    1836     2422077 : uis_357_powermod(ulong x, ulong *mask)
    1837             : {
    1838     2422077 :   if (             !check_res(x, 211UL, 0, mask)) return 0;
    1839      978599 :   if (*mask & 3 && !check_res(x, 209UL, 3, mask)) return 0;
    1840      399592 :   if (*mask & 3 && !check_res(x,  61UL, 6, mask)) return 0;
    1841      222296 :   if (*mask & 5 && !check_res(x, 203UL, 9, mask)) return 0;
    1842       56862 :   if (*mask & 1 && !check_res(x, 117UL,12, mask)) return 0;
    1843       32886 :   if (*mask & 3 && !check_res(x,  31UL,15, mask)) return 0;
    1844       24697 :   if (*mask & 5 && !check_res(x,  43UL,18, mask)) return 0;
    1845        7452 :   if (*mask & 6 && !check_res(x,  71UL,21, mask)) return 0;
    1846        3786 :   return 1;
    1847             : }
    1848             : /* asume x > 0 and pt != NULL */
    1849             : int
    1850     2366983 : uis_357_power(ulong x, ulong *pt, ulong *mask)
    1851             : {
    1852             :   double logx;
    1853     2366983 :   if (!odd(x))
    1854             :   {
    1855        9081 :     long v = vals(x);
    1856        9081 :     if (v % 7) *mask &= ~4;
    1857        9081 :     if (v % 5) *mask &= ~2;
    1858        9081 :     if (v % 3) *mask &= ~1;
    1859        9081 :     if (!*mask) return 0;
    1860             :   }
    1861     2359390 :   if (!uis_357_powermod(x, mask)) return 0;
    1862        2981 :   logx = log((double)x);
    1863        3853 :   while (*mask)
    1864             :   {
    1865             :     long e, b;
    1866             :     ulong y, ye;
    1867        2981 :     if (*mask & 1)      { b = 1; e = 3; }
    1868         873 :     else if (*mask & 2) { b = 2; e = 5; }
    1869         355 :     else                { b = 4; e = 7; }
    1870        2981 :     y = (ulong)(exp(logx / e) + 0.5);
    1871        2981 :     ye = upowuu(y,e);
    1872        2981 :     if (ye == x) { *pt = y; return e; }
    1873             : #ifdef LONG_IS_64BIT
    1874         750 :     if (ye > x) y--; else y++;
    1875         750 :     ye = upowuu(y,e);
    1876         750 :     if (ye == x) { *pt = y; return e; }
    1877             : #endif
    1878         872 :     *mask &= ~b; /* turn the bit off */
    1879             :   }
    1880         872 :   return 0;
    1881             : }
    1882             : 
    1883             : #ifndef LONG_IS_64BIT
    1884             : /* as above, split in two functions */
    1885             : /* is x mod 211*209*61*203 a 3rd, 5th or 7th power ? */
    1886             : static int
    1887       14433 : uis_357_powermod_32bit_1(ulong x, ulong *mask)
    1888             : {
    1889       14433 :   if (             !check_res(x, 211UL, 0, mask)) return 0;
    1890        8008 :   if (*mask & 3 && !check_res(x, 209UL, 3, mask)) return 0;
    1891        4115 :   if (*mask & 3 && !check_res(x,  61UL, 6, mask)) return 0;
    1892        2914 :   if (*mask & 5 && !check_res(x, 203UL, 9, mask)) return 0;
    1893         828 :   return 1;
    1894             : }
    1895             : /* is x mod 117*31*43*71 a 3rd, 5th or 7th power ? */
    1896             : static int
    1897         828 : uis_357_powermod_32bit_2(ulong x, ulong *mask)
    1898             : {
    1899         828 :   if (*mask & 1 && !check_res(x, 117UL,12, mask)) return 0;
    1900         656 :   if (*mask & 3 && !check_res(x,  31UL,15, mask)) return 0;
    1901         538 :   if (*mask & 5 && !check_res(x,  43UL,18, mask)) return 0;
    1902         266 :   if (*mask & 6 && !check_res(x,  71UL,21, mask)) return 0;
    1903         176 :   return 1;
    1904             : }
    1905             : #endif
    1906             : 
    1907             : /* Returns 3, 5, or 7 if x is a cube (but not a 5th or 7th power),  a 5th
    1908             :  * power (but not a 7th),  or a 7th power, and in this case creates the
    1909             :  * base on the stack and assigns its address to *pt.  Otherwise returns 0.
    1910             :  * x must be of type t_INT and positive;  this is not checked.  The *mask
    1911             :  * argument tells us which things to check -- bit 0: 3rd, bit 1: 5th,
    1912             :  * bit 2: 7th pwr;  set a bit to have the corresponding power examined --
    1913             :  * and is updated appropriately for a possible follow-up call */
    1914             : int
    1915     2801660 : is_357_power(GEN x, GEN *pt, ulong *mask)
    1916             : {
    1917     2801660 :   long lx = lgefint(x);
    1918             :   ulong r;
    1919             :   pari_sp av;
    1920             :   GEN y;
    1921             : 
    1922     2801660 :   if (!*mask) return 0; /* useful when running in a loop */
    1923     2430273 :   if (DEBUGLEVEL>4)
    1924           0 :     err_printf("OddPwrs: examining %ld-bit integer\n", expi(x)+1);
    1925     2430273 :   if (lgefint(x) == 3) {
    1926             :     ulong t;
    1927     2353154 :     long e = uis_357_power(x[2], &t, mask);
    1928     2353154 :     if (e)
    1929             :     {
    1930        2083 :       if (pt) *pt = utoi(t);
    1931        2083 :       return e;
    1932             :     }
    1933     2351071 :     return 0;
    1934             :   }
    1935             : #ifdef LONG_IS_64BIT
    1936       62686 :   r = (lx == 3)? uel(x,2): umodiu(x, 6046846918939827UL);
    1937       62686 :   if (!uis_357_powermod(r, mask)) return 0;
    1938             : #else
    1939       14433 :   r = (lx == 3)? uel(x,2): umodiu(x, 211*209*61*203);
    1940       14433 :   if (!uis_357_powermod_32bit_1(r, mask)) return 0;
    1941         828 :   r = (lx == 3)? uel(x,2): umodiu(x, 117*31*43*71);
    1942         828 :   if (!uis_357_powermod_32bit_2(r, mask)) return 0;
    1943             : #endif
    1944         981 :   av = avma;
    1945        1602 :   while (*mask)
    1946             :   {
    1947             :     long e, b;
    1948             :     /* priority to higher powers: if we have a 21st, it is easier to rediscover
    1949             :      * that its 7th root is a cube than that its cube root is a 7th power */
    1950        1303 :          if (*mask & 4) { b = 4; e = 7; }
    1951         963 :     else if (*mask & 2) { b = 2; e = 5; }
    1952         366 :     else                { b = 1; e = 3; }
    1953        1303 :     y = mpround( sqrtnr(itor(x, nbits2prec(64 + bit_accuracy(lx) / e)), e) );
    1954        1303 :     if (equalii(powiu(y,e), x))
    1955             :     {
    1956         682 :       if (!pt) return gc_int(av,e);
    1957         668 :       set_avma((pari_sp)y); *pt = gerepileuptoint(av, y);
    1958         668 :       return e;
    1959             :     }
    1960         621 :     *mask &= ~b; /* turn the bit off */
    1961         621 :     set_avma(av);
    1962             :   }
    1963         299 :   return 0;
    1964             : }
    1965             : 
    1966             : /* Is x a n-th power ? If pt != NULL, it receives the n-th root */
    1967             : ulong
    1968      471038 : is_kth_power(GEN x, ulong n, GEN *pt)
    1969             : {
    1970             :   forprime_t T;
    1971             :   long j;
    1972             :   ulong q, residue;
    1973             :   GEN y;
    1974      471038 :   pari_sp av = avma;
    1975             : 
    1976      471038 :   (void)u_forprime_arith_init(&T, odd(n)? 2*n+1: n+1, ULONG_MAX, 1,n);
    1977             :   /* we'll start at q, smallest prime >= n */
    1978             : 
    1979             :   /* Modular checks, use small primes q congruent 1 mod n */
    1980             :   /* A non n-th power nevertheless passes the test with proba n^(-#checks),
    1981             :    * We'd like this < 1e-6 but let j = floor(log(1e-6) / log(n)) which
    1982             :    * ensures much less. */
    1983      470964 :   if (n < 16)
    1984      118068 :     j = 5;
    1985      352896 :   else if (n < 32)
    1986      154522 :     j = 4;
    1987      198374 :   else if (n < 101)
    1988      176423 :     j = 3;
    1989       21951 :   else if (n < 1001)
    1990        5390 :     j = 2;
    1991       16561 :   else if (n < 17886697) /* smallest such that smallest suitable q is > 2^32 */
    1992       16303 :     j = 1;
    1993             :   else
    1994         275 :     j = 0;
    1995      514422 :   for (; j > 0; j--)
    1996             :   {
    1997      511806 :     if (!(q = u_forprime_next(&T))) break;
    1998             :     /* q a prime = 1 mod n */
    1999      511899 :     residue = umodiu(x, q);
    2000      511893 :     if (residue == 0)
    2001             :     {
    2002         483 :       if (Z_lval(x,q) % n) return gc_ulong(av,0);
    2003          49 :       continue;
    2004             :     }
    2005             :     /* n-th power mod q ? */
    2006      511410 :     if (Fl_powu(residue, (q-1)/n, q) != 1) return gc_ulong(av,0);
    2007             :   }
    2008        2616 :   set_avma(av);
    2009             : 
    2010        2611 :   if (DEBUGLEVEL>4) err_printf("\nOddPwrs: [%lu] passed modular checks\n",n);
    2011             :   /* go to the horse's mouth... */
    2012        2611 :   y = roundr( sqrtnr(itor(x, nbits2prec(expi(x)/n + 16)), n) );
    2013        2611 :   if (!equalii(powiu(y, n), x)) {
    2014        1662 :     if (DEBUGLEVEL>4) err_printf("\tBut it wasn't a pure power.\n");
    2015        1662 :     return gc_ulong(av,0);
    2016             :   }
    2017         949 :   if (!pt) set_avma(av); else { set_avma((pari_sp)y); *pt = gerepileuptoint(av,y); }
    2018         949 :   return 1;
    2019             : }
    2020             : 
    2021             : /* is x a p^i-th power, p >= 11 prime ? Similar to is_357_power(), but instead
    2022             :  * of the mask, we keep the current test exponent around. Cut off when
    2023             :  * log_2 x^(1/k) < cutoffbits since we would have found it by trial division.
    2024             :  * Everything needed here (primitive roots etc.) is computed from scratch on
    2025             :  * the fly; compared to the size of numbers under consideration, these
    2026             :  * word-sized computations take negligible time.
    2027             :  * Any cutoffbits > 0 is safe, but direct root extraction attempts are faster
    2028             :  * when trial division has been used to discover very small bases. We become
    2029             :  * competitive at cutoffbits ~ 10 */
    2030             : int
    2031       69277 : is_pth_power(GEN x, GEN *pt, forprime_t *T, ulong cutoffbits)
    2032             : {
    2033       69277 :   long cnt=0, size = expi(x) /* not +1 */;
    2034             :   ulong p;
    2035       69277 :   pari_sp av = avma;
    2036      486016 :   while ((p = u_forprime_next(T)) && size/p >= cutoffbits) {
    2037      416770 :     long v = 1;
    2038      416770 :     if (DEBUGLEVEL>5 && cnt++==2000)
    2039           0 :       { cnt=0; err_printf("%lu%% ", 100*p*cutoffbits/size); }
    2040      416828 :     while (is_kth_power(x, p, pt)) {
    2041          56 :       v *= p; x = *pt;
    2042          56 :       size = expi(x);
    2043             :     }
    2044      416781 :     if (v > 1)
    2045             :     {
    2046          42 :       if (DEBUGLEVEL>5) err_printf("\nOddPwrs: is a %ld power\n",v);
    2047          42 :       return v;
    2048             :     }
    2049             :   }
    2050       69224 :   if (DEBUGLEVEL>5) err_printf("\nOddPwrs: not a power\n",p);
    2051       69224 :   return gc_int(av,0); /* give up */
    2052             : }
    2053             : 
    2054             : /***********************************************************************/
    2055             : /**                FACTORIZATION  (master iteration)                  **/
    2056             : /**      Driver for the various methods of finding large factors      **/
    2057             : /**      (after trial division has cast out the very small ones).     **/
    2058             : /**                        GN1998Jun24--30                            **/
    2059             : /***********************************************************************/
    2060             : 
    2061             : /* Direct use:
    2062             :  *  ifac_start_hint(n,moebius,hint) registers with the iterative factorizer
    2063             :  *  - an integer n (without prime factors  < tridiv_bound(n))
    2064             :  *  - registers whether or not we should terminate early if we find a square
    2065             :  *    factor,
    2066             :  *  - a hint about which method(s) to use.
    2067             :  *  This must always be called first. If input is not composite, oo loop.
    2068             :  *  The routine decomposes n nontrivially into a product of two factors except
    2069             :  *  in squarefreeness ('Moebius') mode.
    2070             :  *
    2071             :  *  ifac_start(n,moebius) same using default hint.
    2072             :  *
    2073             :  *  ifac_primary_factor()  returns a prime divisor (not necessarily the
    2074             :  *    smallest) and the corresponding exponent.
    2075             :  *
    2076             :  * Encapsulated user interface: Many arithmetic functions have a 'contributor'
    2077             :  * ifac_xxx, to be called on any large composite cofactor left over after trial
    2078             :  * division by small primes: xxx is one of moebius, issquarefree, totient, etc.
    2079             :  *
    2080             :  * We never test whether the input number is prime or composite, since
    2081             :  * presumably it will have come out of the small factors finder stage
    2082             :  * (which doesn't really exist yet but which will test the left-over
    2083             :  * cofactor for primality once it does). */
    2084             : 
    2085             : /* The data structure in which we preserve whatever we know about our number N
    2086             :  * is kept on the PARI stack, and updated as needed.
    2087             :  * This makes the machinery re-entrant, and avoids memory leaks when a lengthy
    2088             :  * factorization is interrupted. We try to keep the whole affair connected,
    2089             :  * and the parent object is always older than its children.  This may in
    2090             :  * rare cases lead to some extra copying around, and knowing what is garbage
    2091             :  * at any given time is not trivial. See below for examples how to do it right.
    2092             :  * (Connectedness is destroyed if callers of ifac_main() create stuff on the
    2093             :  * stack in between calls. This is harmless as long as ifac_realloc() is used
    2094             :  * to re-create a connected object at the head of the stack just before
    2095             :  * collecting garbage.)
    2096             :  * A t_INT may well have > 10^6 distinct prime factors larger than 2^16. Since
    2097             :  * we need not find factors in order of increasing size, we must be prepared to
    2098             :  * drag a very large amount of data around.  We start with a small structure
    2099             :  * and extend it when necessary. */
    2100             : 
    2101             : /* The idea of the algorithm is:
    2102             :  * Let N0 be whatever is currently left of N after dividing off all the
    2103             :  * prime powers we have already returned to the caller.  Then we maintain
    2104             :  * N0 as a product
    2105             :  * (1) N0 = \prod_i P_i^{e_i} * \prod_j Q_j^{f_j} * \prod_k C_k^{g_k}
    2106             :  * where the P_i and Q_j are distinct primes, each C_k is known composite,
    2107             :  * none of the P_i divides any C_k, and we also know the total ordering
    2108             :  * of all the P_i, Q_j and C_k; in particular, we will never try to divide
    2109             :  * a C_k by a larger Q_j.  Some of the C_k may have common factors.
    2110             :  *
    2111             :  * Caveat implementor:  Taking gcds among C_k's is very likely to cost at
    2112             :  * least as much time as dividing off any primes as we find them, and book-
    2113             :  * keeping would be tough (since D=gcd(C_1,C_2) can still have common factors
    2114             :  * with both C_1/D and C_2/D, and so on...).
    2115             :  *
    2116             :  * At startup, we just initialize the structure to
    2117             :  * (2) N = C_1^1   (composite).
    2118             :  *
    2119             :  * Whenever ifac_primary_factor() or one of the arithmetic user interface
    2120             :  * routines needs a primary factor, and the smallest thing in our list is P_1,
    2121             :  * we return that and its exponent, and remove it from our list. (When nothing
    2122             :  * is left, we return a sentinel value -- gen_1.  And in Moebius mode, when we
    2123             :  * see something with exponent > 1, whether prime or composite, we return gen_0
    2124             :  * or 0, depending on the function). In all other cases, ifac_main() iterates
    2125             :  * the following steps until we have a P_1 in the smallest position.
    2126             :  *
    2127             :  * When the smallest item is C_1, as it is initially:
    2128             :  * (3.1) Crack C_1 into a nontrivial product  U_1 * U_2  by whatever method
    2129             :  * comes to mind for this size. (U for 'unknown'.)  Cracking will detect
    2130             :  * perfect powers, so we may instead see a power of some U_1 here, or even
    2131             :  * something of the form U_1^k*U_2^k; of course the exponent already attached
    2132             :  * to C_1 is taken into account in the following.
    2133             :  * (3.2) If we have U_1*U_2, sort the two factors (distinct: squares are caught
    2134             :  * in stage 3.1). N.B. U_1 and U_2 are smaller than anything else in our list.
    2135             :  * (3.3) Check U_1 and U_2 for primality, and flag them accordingly.
    2136             :  * (3.4) Iterate.
    2137             :  *
    2138             :  * When the smallest item is Q_1:
    2139             :  * This is the unpleasant case.  We go through the entire list and try to
    2140             :  * divide Q_1 off each of the current C_k's, which usually fails, but may
    2141             :  * succeed several times. When a division was successful, the corresponding
    2142             :  * C_k is removed from our list, and the cofactor becomes a U_l for the moment
    2143             :  * unless it is 1 (which happens when C_k was a power of Q_1).  When we're
    2144             :  * through we upgrade Q_1 to P_1 status, then do a primality check on each U_l
    2145             :  * and sort it back into the list either as a Q_j or as a C_k.  If during the
    2146             :  * insertion sort we discover that some U_l equals some P_i or Q_j or C_k we
    2147             :  * already have, we just add U_l's exponent to that of its twin. (The sorting
    2148             :  * therefore happens before the primality test). Since this may produce one or
    2149             :  * more elements smaller than the P_1 we just confirmed, we may have to repeat
    2150             :  * the iteration.
    2151             :  * A trick avoids some Q_1 instances: just after the sweep classifying
    2152             :  * all current unknowns as either composites or primes, we do another downward
    2153             :  * sweep beginning with the largest current factor and stopping just above the
    2154             :  * largest current composite.  Every Q_j we pass is turned into a P_i.
    2155             :  * (Different primes are automatically coprime among each other, and primes do
    2156             :  * not divide smaller composites.)
    2157             :  * NB: We have no use for comparing the square of a prime to N0.  Normally
    2158             :  * we will get called after casting out only the smallest primes, and
    2159             :  * since we cannot guarantee that we see the large prime factors in as-
    2160             :  * cending order, we cannot stop when we find one larger than sqrt(N0). */
    2161             : 
    2162             : /* Data structure: We keep everything in a single t_VEC of t_INTs.  The
    2163             :  * first 2 components are read-only:
    2164             :  * 1) the first records whether we're doing full (NULL) or Moebius (gen_1)
    2165             :  * factorization; in the latter case subroutines return a sentinel value as
    2166             :  * soon as they spot an exponent > 1.
    2167             :  * 2) the second records the hint from factorint()'s optional flag, for use by
    2168             :  * ifac_crack().
    2169             :  *
    2170             :  * The remaining components (initially 15) are used in groups of three:
    2171             :  * [ factor (t_INT), exponent (t_INT), factor class ], where factor class is
    2172             :  *  NULL : unknown
    2173             :  *  gen_0: known composite C_k
    2174             :  *  gen_1: known prime Q_j awaiting trial division
    2175             :  *  gen_2: finished prime P_i.
    2176             :  * When during the division stage we re-sort a C_k-turned-U_l to a lower
    2177             :  * position, we rotate any intervening material upward towards its old
    2178             :  * slot.  When a C_k was divided down to 1, its slot is left empty at
    2179             :  * first; similarly when the re-sorting detects a repeated factor.
    2180             :  * After the sorting phase, we de-fragment the list and squeeze all the
    2181             :  * occupied slots together to the high end, so that ifac_crack() has room
    2182             :  * for new factors.  When this doesn't suffice, we abandon the current vector
    2183             :  * and allocate a somewhat larger one, defragmenting again while copying.
    2184             :  *
    2185             :  * For internal use: note that all exponents will fit into C longs, given
    2186             :  * PARI's lgefint field size.  When we work with them, we sometimes read
    2187             :  * out the GEN pointer, and sometimes do an itos, whatever is more con-
    2188             :  * venient for the task at hand. */
    2189             : 
    2190             : /*** Overview ***/
    2191             : 
    2192             : /* The '*where' argument in the following points into *partial at the first of
    2193             :  * the three fields of the first occupied slot.  It's there because the caller
    2194             :  * would already know where 'here' is, so we don't want to search for it again.
    2195             :  * We do not preserve this from one user-interface call to the next. */
    2196             : 
    2197             : /* In the most cases, control flows from the user interface to ifac_main() and
    2198             :  * then to a succession of ifac_crack()s and ifac_divide()s, with (typically)
    2199             :  * none of the latter finding anything. */
    2200             : 
    2201             : #define LAST(x) x+lg(x)-3
    2202             : #define FIRST(x) x+3
    2203             : 
    2204             : #define MOEBIUS(x) gel(x,1)
    2205             : #define HINT(x) gel(x,2)
    2206             : 
    2207             : /* y <- x */
    2208             : INLINE void
    2209           0 : SHALLOWCOPY(GEN x, GEN y) {
    2210           0 :   VALUE(y) = VALUE(x);
    2211           0 :   EXPON(y) = EXPON(x);
    2212           0 :   CLASS(y) = CLASS(x);
    2213           0 : }
    2214             : /* y <- x */
    2215             : INLINE void
    2216           0 : COPY(GEN x, GEN y) {
    2217           0 :   icopyifstack(VALUE(x), VALUE(y));
    2218           0 :   icopyifstack(EXPON(x), EXPON(y));
    2219           0 :   CLASS(y) = CLASS(x);
    2220           0 : }
    2221             : 
    2222             : /* Diagnostics */
    2223             : static void
    2224           0 : ifac_factor_dbg(GEN x)
    2225             : {
    2226           0 :   GEN c = CLASS(x), v = VALUE(x);
    2227           0 :   if (c == gen_2) err_printf("IFAC: factor %Ps\n\tis prime (finished)\n", v);
    2228           0 :   else if (c == gen_1) err_printf("IFAC: factor %Ps\n\tis prime\n", v);
    2229           0 :   else if (c == gen_0) err_printf("IFAC: factor %Ps\n\tis composite\n", v);
    2230           0 : }
    2231             : static void
    2232           0 : ifac_check(GEN partial, GEN where)
    2233             : {
    2234           0 :   if (!where || where < FIRST(partial) || where > LAST(partial))
    2235           0 :     pari_err_BUG("ifac_check ['where' out of bounds]");
    2236           0 : }
    2237             : static void
    2238           0 : ifac_print(GEN part, GEN where)
    2239             : {
    2240           0 :   long l = lg(part);
    2241             :   GEN p;
    2242             : 
    2243           0 :   err_printf("ifac partial factorization structure: %ld slots, ", (l-3)/3);
    2244           0 :   if (MOEBIUS(part)) err_printf("Moebius mode, ");
    2245           0 :   err_printf("hint = %ld\n", itos(HINT(part)));
    2246           0 :   ifac_check(part, where);
    2247           0 :   for (p = part+3; p < part + l; p += 3)
    2248             :   {
    2249           0 :     GEN v = VALUE(p), e = EXPON(p), c = CLASS(p);
    2250           0 :     const char *s = "";
    2251           0 :     if (!v) { err_printf("[empty slot]\n"); continue; }
    2252           0 :     if (c == NULL) s = "unknown";
    2253           0 :     else if (c == gen_0) s = "composite";
    2254           0 :     else if (c == gen_1) s = "unfinished prime";
    2255           0 :     else if (c == gen_2) s = "prime";
    2256           0 :     else pari_err_BUG("unknown factor class");
    2257           0 :     err_printf("[%Ps, %Ps, %s]\n", v, e, s);
    2258             :   }
    2259           0 :   err_printf("Done.\n");
    2260           0 : }
    2261             : 
    2262             : static const long decomp_default_hint = 0;
    2263             : /* assume n > 0, which we can assign to */
    2264             : /* return initial data structure, see ifac_crack() for the hint argument */
    2265             : static GEN
    2266        5940 : ifac_start_hint(GEN n, int moebius, long hint)
    2267             : {
    2268        5940 :   const long ifac_initial_length = 3 + 7*3;
    2269             :   /* codeword, moebius, hint, 7 slots -- a 512-bit product of distinct 8-bit
    2270             :    * primes needs at most 7 slots at a time) */
    2271        5940 :   GEN here, part = cgetg(ifac_initial_length, t_VEC);
    2272             : 
    2273        5940 :   MOEBIUS(part) = moebius? gen_1 : NULL;
    2274        5940 :   HINT(part) = stoi(hint);
    2275             :   /* fill first slot at the top end */
    2276        5940 :   here = part + ifac_initial_length - 3; /* LAST(part) */
    2277        5940 :   INIT(here, n,gen_1,gen_0); /* n^1: composite */
    2278       41580 :   while ((here -= 3) > part) ifac_delete(here);
    2279        5940 :   return part;
    2280             : }
    2281             : GEN
    2282        2509 : ifac_start(GEN n, int moebius)
    2283        2509 : { return ifac_start_hint(n,moebius,decomp_default_hint); }
    2284             : 
    2285             : /* Return next nonempty slot after 'here', NULL if none exist */
    2286             : static GEN
    2287       15293 : ifac_find(GEN partial)
    2288             : {
    2289       15293 :   GEN scan, end = partial + lg(partial);
    2290             : 
    2291             : #ifdef IFAC_DEBUG
    2292             :   ifac_check(partial, partial);
    2293             : #endif
    2294      111505 :   for (scan = partial+3; scan < end; scan += 3)
    2295      106781 :     if (VALUE(scan)) return scan;
    2296        4724 :   return NULL;
    2297             : }
    2298             : 
    2299             : /* Defragment: squeeze out unoccupied slots above *where. Unoccupied slots
    2300             :  * arise when a composite factor dissolves completely whilst dividing off a
    2301             :  * prime, or when ifac_resort() spots a coincidence and merges two factors.
    2302             :  * Update *where */
    2303             : static void
    2304          14 : ifac_defrag(GEN *partial, GEN *where)
    2305             : {
    2306          14 :   GEN scan_new = LAST(*partial), scan_old;
    2307             : 
    2308          42 :   for (scan_old = scan_new; scan_old >= *where; scan_old -= 3)
    2309             :   {
    2310          28 :     if (!VALUE(scan_old)) continue; /* empty slot */
    2311          28 :     if (scan_old < scan_new) SHALLOWCOPY(scan_old, scan_new);
    2312          28 :     scan_new -= 3; /* point at next slot to be written */
    2313             :   }
    2314          14 :   scan_new += 3; /* back up to last slot written */
    2315          14 :   *where = scan_new;
    2316          84 :   while ((scan_new -= 3) > *partial) ifac_delete(scan_new); /* erase junk */
    2317          14 : }
    2318             : 
    2319             : /* Move to a larger main vector, updating *where if it points into it, and
    2320             :  * *partial in any case. Can be used as a specialized gcopy before
    2321             :  * a gerepileupto() (pass 0 as the new length). Normally, one would pass
    2322             :  * new_lg=1 to let this function guess the new size.  To be used sparingly.
    2323             :  * Complex version of ifac_defrag(), combined with reallocation.  If new_lg
    2324             :  * is 0, use the old length, so this acts just like gcopy except that the
    2325             :  * 'where' pointer is carried along; if it is 1, we make an educated guess.
    2326             :  * Exception:  If new_lg is 0, the vector is full to the brim, and the first
    2327             :  * entry is composite, we make it longer to avoid being called again a
    2328             :  * microsecond later. It is safe to call this with *where = NULL:
    2329             :  * if it doesn't point anywhere within the old structure, it is left alone */
    2330             : static void
    2331           0 : ifac_realloc(GEN *partial, GEN *where, long new_lg)
    2332             : {
    2333           0 :   long old_lg = lg(*partial);
    2334             :   GEN newpart, scan_new, scan_old;
    2335             : 
    2336           0 :   if (new_lg == 1)
    2337           0 :     new_lg = 2*old_lg - 6;        /* from 7 slots to 13 to 25... */
    2338           0 :   else if (new_lg <= old_lg)        /* includes case new_lg == 0 */
    2339             :   {
    2340           0 :     GEN first = *partial + 3;
    2341           0 :     new_lg = old_lg;
    2342             :     /* structure full and first entry composite or unknown */
    2343           0 :     if (VALUE(first) && (CLASS(first) == gen_0 || CLASS(first)==NULL))
    2344           0 :       new_lg += 6; /* give it a little more breathing space */
    2345             :   }
    2346           0 :   newpart = cgetg(new_lg, t_VEC);
    2347           0 :   if (DEBUGMEM >= 3)
    2348           0 :     err_printf("IFAC: new partial factorization structure (%ld slots)\n",
    2349           0 :                (new_lg - 3)/3);
    2350           0 :   MOEBIUS(newpart) = MOEBIUS(*partial);
    2351           0 :   icopyifstack(HINT(*partial), HINT(newpart));
    2352             :   /* Downward sweep through the old *partial. Pick up 'where' and carry it
    2353             :    * over if we pass it. (Only useful if it pointed at a nonempty slot.)
    2354             :    * Factors are COPY'd so that we again have a nice object (parent older
    2355             :    * than children, connected), except the one factor that may still be living
    2356             :    * in a clone where n originally was; exponents are similarly copied if they
    2357             :    * aren't global constants; class-of-factor fields are global constants so we
    2358             :    * need only copy them as pointers. Caller may then do a gerepileupto() */
    2359           0 :   scan_new = newpart + new_lg - 3; /* LAST(newpart) */
    2360           0 :   scan_old = *partial + old_lg - 3; /* LAST(*partial) */
    2361           0 :   for (; scan_old > *partial + 2; scan_old -= 3)
    2362             :   {
    2363           0 :     if (*where == scan_old) *where = scan_new;
    2364           0 :     if (!VALUE(scan_old)) continue; /* skip empty slots */
    2365           0 :     COPY(scan_old, scan_new); scan_new -= 3;
    2366             :   }
    2367           0 :   scan_new += 3; /* back up to last slot written */
    2368           0 :   while ((scan_new -= 3) > newpart) ifac_delete(scan_new);
    2369           0 :   *partial = newpart;
    2370           0 : }
    2371             : 
    2372             : /* Re-sort one (typically unknown) entry from washere to a new position,
    2373             :  * rotating intervening entries upward to fill the vacant space. If the new
    2374             :  * position is the same as the old one, or the new value of the entry coincides
    2375             :  * with a value already occupying a lower slot, then we just add exponents (and
    2376             :  * use the 'more known' class, and return 1 immediately when in Moebius mode).
    2377             :  * Slots between *where and washere must be in sorted order, so a sweep using
    2378             :  * this to re-sort several unknowns must proceed upward, see ifac_resort().
    2379             :  * Bubble-sort-of-thing sort. Won't be exercised frequently, so this is ok */
    2380             : static void
    2381           7 : ifac_sort_one(GEN *where, GEN washere)
    2382             : {
    2383           7 :   GEN old, scan = washere - 3;
    2384             :   GEN value, exponent, class0, class1;
    2385             :   long cmp_res;
    2386             : 
    2387           7 :   if (scan < *where) return; /* nothing to do, washere==*where */
    2388           7 :   value    = VALUE(washere);
    2389           7 :   exponent = EXPON(washere);
    2390           7 :   class0 = CLASS(washere);
    2391           7 :   cmp_res = -1; /* sentinel */
    2392           7 :   while (scan >= *where) /* at least once */
    2393             :   {
    2394           7 :     if (VALUE(scan))
    2395             :     { /* current slot nonempty, check against where */
    2396           7 :       cmp_res = cmpii(value, VALUE(scan));
    2397           7 :       if (cmp_res >= 0) break; /* have found where to stop */
    2398             :     }
    2399             :     /* copy current slot upward by one position and move pointers down */
    2400           0 :     SHALLOWCOPY(scan, scan+3);
    2401           0 :     scan -= 3;
    2402             :   }
    2403           7 :   scan += 3;
    2404             :   /* At this point there are the following possibilities:
    2405             :    * 1) cmp_res == -1. Either value is less than that at *where, or *where was
    2406             :    * pointing at vacant slots and any factors we saw en route were larger than
    2407             :    * value. At any rate, scan == *where now, and scan is pointing at an empty
    2408             :    * slot, into which we'll stash our entry.
    2409             :    * 2) cmp_res == 0. The entry at scan-3 is the one, we compare class0
    2410             :    * fields and add exponents, and put it all into the vacated scan slot,
    2411             :    * NULLing the one at scan-3 (and possibly updating *where).
    2412             :    * 3) cmp_res == 1. The slot at scan is the one to store our entry into. */
    2413           7 :   if (cmp_res)
    2414             :   {
    2415           7 :     if (cmp_res < 0 && scan != *where)
    2416           0 :       pari_err_BUG("ifact_sort_one [misaligned partial]");
    2417           7 :     INIT(scan, value, exponent, class0); return;
    2418             :   }
    2419             :   /* case cmp_res == 0: repeated factor detected */
    2420           0 :   if (DEBUGLEVEL >= 4)
    2421           0 :     err_printf("IFAC: repeated factor %Ps\n\tin ifac_sort_one\n", value);
    2422           0 :   old = scan - 3;
    2423             :   /* if old class0 was composite and new is prime, or vice versa, complain
    2424             :    * (and if one class0 was unknown and the other wasn't, use the known one) */
    2425           0 :   class1 = CLASS(old);
    2426           0 :   if (class0) /* should never be used */
    2427             :   {
    2428           0 :     if (class1)
    2429             :     {
    2430           0 :       if (class0 == gen_0 && class1 != gen_0)
    2431           0 :         pari_err_BUG("ifac_sort_one (composite = prime)");
    2432           0 :       else if (class0 != gen_0 && class1 == gen_0)
    2433           0 :         pari_err_BUG("ifac_sort_one (prime = composite)");
    2434           0 :       else if (class0 == gen_2)
    2435           0 :         CLASS(scan) = class0;
    2436             :     }
    2437             :     else
    2438           0 :       CLASS(scan) = class0;
    2439             :   }
    2440             :   /* else stay with the existing known class0 */
    2441           0 :   CLASS(scan) = class1;
    2442             :   /* in any case, add exponents */
    2443           0 :   if (EXPON(old) == gen_1 && exponent == gen_1)
    2444           0 :     EXPON(scan) = gen_2;
    2445             :   else
    2446           0 :     EXPON(scan) = addii(EXPON(old), exponent);
    2447             :   /* move the value over and null out the vacated slot below */
    2448           0 :   old = scan - 3;
    2449           0 :   *scan = *old;
    2450           0 :   ifac_delete(old);
    2451             :   /* finally, see whether *where should be pulled in */
    2452           0 :   if (old == *where) *where += 3;
    2453             : }
    2454             : 
    2455             : /* Sort all current unknowns downward to where they belong. Sweeps in the
    2456             :  * upward direction. Not needed after ifac_crack(), only when ifac_divide()
    2457             :  * returned true. Update *where. */
    2458             : static void
    2459           7 : ifac_resort(GEN *partial, GEN *where)
    2460             : {
    2461             :   GEN scan, end;
    2462           7 :   ifac_defrag(partial, where); end = LAST(*partial);
    2463          21 :   for (scan = *where; scan <= end; scan += 3)
    2464          14 :     if (VALUE(scan) && !CLASS(scan)) ifac_sort_one(where, scan); /*unknown*/
    2465           7 :   ifac_defrag(partial, where); /* remove newly created gaps */
    2466           7 : }
    2467             : 
    2468             : /* Let x be a t_INT known not to have small divisors (< 661, and < 2^14 for huge
    2469             :  * x > 2^512). Return 0 if x is a proven composite. Return 1 if we believe it
    2470             :  * to be prime (fully proven prime if factor_proven is set).  */
    2471             : int
    2472       28148 : ifac_isprime(GEN x)
    2473             : {
    2474       28148 :   if (!BPSW_psp_nosmalldiv(x)) return 0; /* composite */
    2475       19818 :   if (factor_proven && ! BPSW_isprime(x))
    2476             :   {
    2477           0 :     pari_warn(warner,
    2478             :               "IFAC: pseudo-prime %Ps\n\tis not prime. PLEASE REPORT!\n", x);
    2479           0 :     return 0;
    2480             :   }
    2481       19818 :   return 1;
    2482             : }
    2483             : 
    2484             : static int
    2485       11291 : ifac_checkprime(GEN x)
    2486             : {
    2487       11291 :   int res = ifac_isprime(VALUE(x));
    2488       11291 :   CLASS(x) = res? gen_1: gen_0;
    2489       11291 :   if (DEBUGLEVEL>2) ifac_factor_dbg(x);
    2490       11291 :   return res;
    2491             : }
    2492             : 
    2493             : /* Determine primality or compositeness of all current unknowns, and set
    2494             :  * class Q primes to finished (class P) if everything larger is already
    2495             :  * known to be prime.  When after_crack >= 0, only look at the
    2496             :  * first after_crack things in the list (do nothing when it's 0) */
    2497             : static void
    2498        5794 : ifac_whoiswho(GEN *partial, GEN *where, long after_crack)
    2499             : {
    2500        5794 :   GEN scan, scan_end = LAST(*partial);
    2501             : 
    2502             : #ifdef IFAC_DEBUG
    2503             :   ifac_check(*partial, *where);
    2504             : #endif
    2505        5794 :   if (after_crack == 0) return;
    2506        5180 :   if (after_crack > 0) /* check at most after_crack entries */
    2507        5173 :     scan = *where + 3*(after_crack - 1); /* assert(scan <= scan_end) */
    2508             :   else
    2509           7 :     for (scan = scan_end; scan >= *where; scan -= 3)
    2510             :     {
    2511           7 :       if (CLASS(scan))
    2512             :       { /* known class of factor */
    2513           0 :         if (CLASS(scan) == gen_0) break;
    2514           0 :         if (CLASS(scan) == gen_1)
    2515             :         {
    2516           0 :           if (DEBUGLEVEL>=3)
    2517             :           {
    2518           0 :             err_printf("IFAC: factor %Ps\n\tis prime (no larger composite)\n",
    2519           0 :                        VALUE(*where));
    2520           0 :             err_printf("IFAC: prime %Ps\n\tappears with exponent = %ld\n",
    2521           0 :                        VALUE(*where), itos(EXPON(*where)));
    2522             :           }
    2523           0 :           CLASS(scan) = gen_2;
    2524             :         }
    2525           0 :         continue;
    2526             :       }
    2527           7 :       if (!ifac_checkprime(scan)) break; /* must disable Q-to-P */
    2528           0 :       CLASS(scan) = gen_2; /* P_i, finished prime */
    2529           0 :       if (DEBUGLEVEL>2) ifac_factor_dbg(scan);
    2530             :     }
    2531             :   /* go on, Q-to-P trick now disabled */
    2532       15809 :   for (; scan >= *where; scan -= 3)
    2533             :   {
    2534       10629 :     if (CLASS(scan)) continue;
    2535       10594 :     (void)ifac_checkprime(scan); /* Qj | Ck */
    2536             :   }
    2537             : }
    2538             : 
    2539             : /* Divide all current composites by first (prime, class Q) entry, updating its
    2540             :  * exponent, and turning it into a finished prime (class P).  Return 1 if any
    2541             :  * such divisions succeeded  (in Moebius mode, the update may then not have
    2542             :  * been completed), or 0 if none of them succeeded.  Doesn't modify *where.
    2543             :  * Here we normally do not check that the first entry is a not-finished
    2544             :  * prime.  Stack management: we may allocate a new exponent */
    2545             : static long
    2546        9847 : ifac_divide(GEN *partial, GEN *where, long moebius_mode)
    2547             : {
    2548        9847 :   GEN scan, scan_end = LAST(*partial);
    2549        9847 :   long res = 0, exponent, newexp, otherexp;
    2550             : 
    2551             : #ifdef IFAC_DEBUG
    2552             :   ifac_check(*partial, *where);
    2553             :   if (CLASS(*where) != gen_1)
    2554             :     pari_err_BUG("ifac_divide [division by composite or finished prime]");
    2555             :   if (!VALUE(*where)) pari_err_BUG("ifac_divide [division by nothing]");
    2556             : #endif
    2557        9847 :   newexp = exponent = itos(EXPON(*where));
    2558        9847 :   if (exponent > 1 && moebius_mode) return 1;
    2559             :   /* should've been caught by caller */
    2560             : 
    2561       15559 :   for (scan = *where+3; scan <= scan_end; scan += 3)
    2562             :   {
    2563        5712 :     if (CLASS(scan) != gen_0) continue; /* the other thing ain't composite */
    2564         205 :     otherexp = 0;
    2565             :     /* divide in place to keep stack clutter minimal */
    2566         219 :     while (dvdiiz(VALUE(scan), VALUE(*where), VALUE(scan)))
    2567             :     {
    2568          14 :       if (moebius_mode) return 1; /* immediately */
    2569          14 :       if (!otherexp) otherexp = itos(EXPON(scan));
    2570          14 :       newexp += otherexp;
    2571             :     }
    2572         205 :     if (newexp > exponent)        /* did anything happen? */
    2573             :     {
    2574           7 :       EXPON(*where) = (newexp == 2 ? gen_2 : utoipos(newexp));
    2575           7 :       exponent = newexp;
    2576           7 :       if (is_pm1((GEN)*scan)) /* factor dissolved completely */
    2577             :       {
    2578           0 :         ifac_delete(scan);
    2579           0 :         if (DEBUGLEVEL >= 4)
    2580           0 :           err_printf("IFAC: a factor was a power of another prime factor\n");
    2581             :       } else {
    2582           7 :         CLASS(scan) = NULL;        /* at any rate it's Unknown now */
    2583           7 :         if (DEBUGLEVEL >= 4)
    2584           0 :           err_printf("IFAC: a factor was divisible by another prime factor,\n"
    2585             :                      "\tleaving a cofactor = %Ps\n", VALUE(scan));
    2586             :       }
    2587           7 :       res = 1;
    2588           7 :       if (DEBUGLEVEL >= 5)
    2589           0 :         err_printf("IFAC: prime %Ps\n\tappears at least to the power %ld\n",
    2590           0 :                    VALUE(*where), newexp);
    2591             :     }
    2592             :   } /* for */
    2593        9847 :   CLASS(*where) = gen_2; /* make it a finished prime */
    2594        9847 :   if (DEBUGLEVEL >= 3)
    2595           0 :     err_printf("IFAC: prime %Ps\n\tappears with exponent = %ld\n",
    2596           0 :                VALUE(*where), newexp);
    2597        9847 :   return res;
    2598             : }
    2599             : 
    2600             : /* found out our integer was factor^exp. Update */
    2601             : static void
    2602         880 : update_pow(GEN where, GEN factor, long exp, pari_sp *av)
    2603             : {
    2604         880 :   GEN ex = EXPON(where);
    2605         880 :   if (DEBUGLEVEL>3)
    2606           0 :     err_printf("IFAC: found %Ps =\n\t%Ps ^%ld\n", *where, factor, exp);
    2607         880 :   affii(factor, VALUE(where)); set_avma(*av);
    2608         880 :   if (ex == gen_1)
    2609         677 :   { EXPON(where) = exp == 2? gen_2: utoipos(exp); *av = avma; }
    2610         203 :   else if (ex == gen_2)
    2611         182 :   { EXPON(where) = utoipos(exp<<1); *av = avma; }
    2612             :   else
    2613          21 :     affsi(exp * itos(ex), EXPON(where));
    2614         880 : }
    2615             : /* hint = 0 : Use a default strategy
    2616             :  * hint & 1 : avoid MPQS
    2617             :  * hint & 2 : avoid first-stage ECM (may fall back to ECM if MPQS gives up)
    2618             :  * hint & 4 : avoid Pollard and SQUFOF stages.
    2619             :  * hint & 8 : avoid final ECM; may flag a composite as prime. */
    2620             : #define get_hint(partial) (itos(HINT(*partial)) & 15)
    2621             : 
    2622             : /* Complete ifac_crack's job when a factoring engine splits the current factor
    2623             :  * into a product of three or more new factors. Makes room for them if
    2624             :  * necessary, sorts them, gives them the right exponents and class. Returns the
    2625             :  * number of factors actually written, which may be less than #facvec if there
    2626             :  * are duplicates. Vectors of factors (cf pollardbrent() or mpqs()) contain
    2627             :  * 'slots' of 3 GENs per factor, interpreted as in our partial factorization
    2628             :  * data structure. Thus engines can tell us what they already know about
    2629             :  * factors being prime/composite or appearing to a power larger than thefirst.
    2630             :  * Don't collect garbage.  No diagnostics: engine has printed what it found.
    2631             :  * facvec contains slots of three components per factor; repeated factors are
    2632             :  * allowed (their classes shouldn't contradict each other whereas their
    2633             :  * exponents will be added up) */
    2634             : static long
    2635        3284 : ifac_insert_multiplet(GEN *partial, GEN *where, GEN facvec, long moebius_mode)
    2636             : {
    2637        3284 :   long j,k=1, lfv=lg(facvec)-1, nf=lfv/3, room=(long)(*where-*partial);
    2638             :   /* one of the factors will go into the *where slot, so room is now 3 times
    2639             :    * the number of slots we can use */
    2640        3284 :   long needroom = lfv - room;
    2641        3284 :   GEN newexp, cur, sorted, auxvec = cgetg(nf+1, t_VEC), factor;
    2642        3284 :   long E = itos(EXPON(*where)); /* the old exponent */
    2643             : 
    2644        3284 :   if (DEBUGLEVEL >= 5) /* squfof may return a single squared factor as a set */
    2645           0 :     err_printf("IFAC: incorporating set of %ld factor(s)\n", nf);
    2646        3284 :   if (needroom > 0) ifac_realloc(partial, where, lg(*partial) + needroom);
    2647             : 
    2648             :   /* create sort permutation from the values of the factors */
    2649       10128 :   for (j=nf; j; j--) gel(auxvec,j) = gel(facvec,3*j-2);
    2650        3284 :   sorted = ZV_indexsort(auxvec);
    2651             :   /* store factors, beginning at *where, and catching any duplicates */
    2652        3284 :   cur = facvec + 3*sorted[nf]-2; /* adjust for triple spacing */
    2653        3284 :   VALUE(*where) = VALUE(cur);
    2654        3284 :   newexp = EXPON(cur);
    2655             :   /* if new exponent is 1, the old exponent already in place will do */
    2656        3284 :   if (newexp != gen_1) EXPON(*where) = mului(E,newexp);
    2657        3284 :   CLASS(*where) = CLASS(cur);
    2658        3284 :   if (DEBUGLEVEL >= 6) err_printf("\tstored (largest) factor no. %ld...\n", nf);
    2659             : 
    2660        6844 :   for (j=nf-1; j; j--)
    2661             :   {
    2662             :     GEN e;
    2663        3560 :     cur = facvec + 3*sorted[j]-2;
    2664        3560 :     factor = VALUE(cur);
    2665        3560 :     if (equalii(factor, VALUE(*where)))
    2666             :     {
    2667           0 :       if (DEBUGLEVEL >= 6)
    2668           0 :         err_printf("\tfactor no. %ld is a duplicate%s\n", j, (j>1? "...": ""));
    2669             :       /* update exponent, ignore class which would already have been set,
    2670             :        * then forget current factor */
    2671           0 :       newexp = EXPON(cur);
    2672           0 :       if (newexp != gen_1) /* new exp > 1 */
    2673           0 :         e = addiu(EXPON(*where), E * itou(newexp));
    2674           0 :       else if (EXPON(*where) == gen_1 && E == 1)
    2675           0 :         e = gen_2;
    2676             :       else
    2677           0 :         e = addiu(EXPON(*where), E);
    2678           0 :       EXPON(*where) = e;
    2679             : 
    2680           0 :       if (moebius_mode) return 0; /* stop now, with exponent updated */
    2681           0 :       continue;
    2682             :     }
    2683             : 
    2684        3560 :     *where -= 3;
    2685        3560 :     CLASS(*where) = CLASS(cur); /* class as given */
    2686        3560 :     newexp = EXPON(cur);
    2687        3560 :     if (newexp != gen_1) /* new exp > 1 */
    2688          99 :       e = (E == 1 && newexp == gen_2)? gen_2: mului(E, newexp);
    2689             :     else /* inherit parent's exponent */
    2690        3461 :       e = (E == 1 ? gen_1 : (E == 2 ? gen_2 : utoipos(E)));
    2691        3560 :     EXPON(*where) = e;
    2692             :     /* keep components younger than *partial */
    2693        3560 :     icopyifstack(factor, VALUE(*where));
    2694        3560 :     k++;
    2695        3560 :     if (DEBUGLEVEL >= 6)
    2696           0 :       err_printf("\tfactor no. %ld was unique%s\n", j, j>1? " (so far)...": "");
    2697             :   }
    2698        3284 :   return k;
    2699             : }
    2700             : 
    2701             : /* x /= y; exact division */
    2702             : static void
    2703        1883 : diviiexact_inplace(GEN x, GEN y)
    2704        1883 : { pari_sp av = avma; affii(diviiexact(x, y), x); set_avma(av); }
    2705             : 
    2706             : /* Split the first (composite) entry.  There must already be room for another
    2707             :  * factor below *where, and *where is updated. Two cases:
    2708             :  * - entry is a pure power: factor^k is inserted, leaving *where unchanged;
    2709             :  * - entry = factor * cofactor (not necessarily coprime): both factors are
    2710             :  *   inserted in the correct order, updating *where
    2711             :  * The inserted factors class is set to unknown, they inherit the exponent
    2712             :  * (or a multiple thereof) of their ancestor.
    2713             :  *
    2714             :  * Returns number of factors written into the structure, usually 2; only 1
    2715             :  * if pure power, and > 2 if a factoring engine returned a vector of factors.
    2716             :  * Can reallocate the data structure in the rare > 2 case; may create one or
    2717             :  * more objects: new factors or exponents > 2 */
    2718             : static long
    2719        5789 : ifac_crack(GEN *partial, GEN *where, long moebius_mode)
    2720             : {
    2721        5789 :   long hint = get_hint(partial);
    2722             :   GEN cofactor, factor, exponent;
    2723             : 
    2724             : #ifdef IFAC_DEBUG
    2725             :   ifac_check(*partial, *where);
    2726             :   if (*where < *partial + 6)
    2727             :     pari_err_BUG("ifac_crack ['*where' out of bounds]");
    2728             :   if (!(VALUE(*where)) || typ(VALUE(*where)) != t_INT)
    2729             :     pari_err_BUG("ifac_crack [incorrect VALUE(*where)]");
    2730             :   if (CLASS(*where) != gen_0)
    2731             :     pari_err_BUG("ifac_crack [operand not known composite]");
    2732             : #endif
    2733             : 
    2734        5789 :   if (DEBUGLEVEL>2) {
    2735           0 :     err_printf("IFAC: cracking composite\n\t%Ps\n", **where);
    2736           0 :     if (DEBUGLEVEL>3) err_printf("IFAC: checking for pure square\n");
    2737             :   }
    2738             :   /* MPQS cannot factor prime powers. Look for pure powers even if MPQS is
    2739             :    * blocked by hint: fast and useful in bounded factorization */
    2740             :   {
    2741             :     forprime_t T;
    2742        5789 :     ulong exp = 1, mask = 7;
    2743        5789 :     long good = 0;
    2744        5789 :     pari_sp av = avma;
    2745        5789 :     (void)u_forprime_init(&T, 11, ULONG_MAX);
    2746             :     /* crack squares */
    2747        6626 :     while (Z_issquareall(VALUE(*where), &factor))
    2748             :     {
    2749         838 :       good = 1; /* remember we succeeded once */
    2750         838 :       update_pow(*where, factor, 2, &av);
    2751        1452 :       if (moebius_mode) return 0; /* no need to carry on */
    2752             :     }
    2753        5830 :     while ( (exp = is_357_power(VALUE(*where), &factor, &mask)) )
    2754             :     {
    2755          42 :       good = 1; /* remember we succeeded once */
    2756          42 :       update_pow(*where, factor, exp, &av);
    2757          42 :       if (moebius_mode) return 0; /* no need to carry on */
    2758             :     }
    2759             :     /* cutoff at 14 bits: OK if tridiv_bound >= 2^14 OR if >= 661 for
    2760             :      * an integer < 701^11 (103 bits). */
    2761        5788 :     while ( (exp = is_pth_power(VALUE(*where), &factor, &T, 15)) )
    2762             :     {
    2763           0 :       good = 1; /* remember we succeeded once */
    2764           0 :       update_pow(*where, factor, exp, &av);
    2765           0 :       if (moebius_mode) return 0; /* no need to carry on */
    2766             :     }
    2767             : 
    2768        5788 :     if (good && hint != 15 && ifac_checkprime(*where))
    2769             :     { /* our composite was a prime power */
    2770         614 :       if (DEBUGLEVEL>3)
    2771           0 :         err_printf("IFAC: factor %Ps\n\tis prime\n", VALUE(*where));
    2772         614 :       return 0; /* bypass subsequent ifac_whoiswho() call */
    2773             :     }
    2774             :   } /* pure power stage */
    2775             : 
    2776        5174 :   factor = NULL;
    2777        5174 :   if (!(hint & 4))
    2778             :   { /* SQUFOF then Rho */
    2779        5125 :     if (DEBUGLEVEL >= 4)
    2780           0 :       err_printf("IFAC: trying Shanks' SQUFOF, will fail silently if input\n"
    2781             :                  "      is too large for it.\n");
    2782        5125 :     factor = squfof(VALUE(*where));
    2783        5125 :     if (!factor)
    2784             :     {
    2785        3494 :       if (DEBUGLEVEL >= 4) err_printf("IFAC: trying Pollard-Brent rho\n");
    2786        3494 :       factor = pollardbrent(VALUE(*where));
    2787             :     }
    2788             :   }
    2789        5174 :   if (!factor && !(hint & 2))
    2790             :   { /* First ECM stage */
    2791        3267 :     if (DEBUGLEVEL >= 4) err_printf("IFAC: trying Lenstra-Montgomery ECM\n");
    2792        3267 :     factor = ellfacteur(VALUE(*where), 0); /* do not insist */
    2793             :   }
    2794        5174 :   if (!factor && !(hint & 1))
    2795             :   { /* MPQS stage */
    2796        3291 :     if (DEBUGLEVEL >= 4) err_printf("IFAC: trying MPQS\n");
    2797        3291 :     factor = mpqs(VALUE(*where));
    2798             :   }
    2799        5174 :   if (!factor && !(hint & 8))
    2800             :   { /* Final ECM stage, guaranteed to succeed */
    2801           0 :     if (DEBUGLEVEL >= 4)
    2802           0 :       err_printf("IFAC: forcing ECM, may take some time\n");
    2803           0 :     factor = ellfacteur(VALUE(*where), 1);
    2804             :   }
    2805        5174 :   if (!factor)
    2806             :   { /* limited factorization */
    2807           7 :     if (DEBUGLEVEL >= 2)
    2808             :     {
    2809           0 :       pari_warn(warner, hint==15? "IFAC: untested integer declared prime"
    2810             :                                 : "IFAC: unfactored composite declared prime");
    2811             :       /* don't print it out at level 3 or above, where it would appear
    2812             :        * several times before and after this message already */
    2813           0 :       if (DEBUGLEVEL == 2) err_printf("\t%Ps\n", VALUE(*where));
    2814             :     }
    2815           7 :     CLASS(*where) = gen_1; /* might as well trial-divide by it... */
    2816           7 :     return 1;
    2817             :   }
    2818             :   /* At least two factors */
    2819        5167 :   if (typ(factor) == t_VEC)
    2820        3284 :     return ifac_insert_multiplet(partial, where, factor, moebius_mode);
    2821             : 
    2822             :   /* Single factor (t_INT): work out cofactor in place */
    2823        1883 :   diviiexact_inplace(VALUE(*where), factor);
    2824        1883 :   cofactor = VALUE(*where);
    2825             :   /* factoring engines reported factor; tell about the cofactor */
    2826        1883 :   if (DEBUGLEVEL >= 4) err_printf("IFAC: cofactor = %Ps\n", cofactor);
    2827        1883 :   CLASS(*where) = NULL; /* mark factor /cofactor 'unknown' */
    2828        1883 :   exponent = EXPON(*where); /* common exponent */
    2829             : 
    2830        1883 :   *where -= 3;
    2831        1883 :   CLASS(*where) = NULL; /* mark factor /cofactor 'unknown' */
    2832        1883 :   icopyifstack(exponent, EXPON(*where)); /* copy common exponent */
    2833             : 
    2834        1883 :   if (cmpii(factor, cofactor) < 0)
    2835        1705 :     VALUE(*where) = factor; /* common case */
    2836             :   else
    2837             :   { /* factor > cofactor, swap. Exponent is the same, so no need to swap. */
    2838         178 :     GEN old = *where + 3;
    2839         178 :     VALUE(*where) = VALUE(old); /* move cofactor pointer to lowest slot */
    2840         178 :     VALUE(old) = factor; /* save factor */
    2841             :   }
    2842        1883 :   return 2;
    2843             : }
    2844             : 
    2845             : /* main loop:  iterate until smallest entry is a finished prime;  returns
    2846             :  * a 'where' pointer, or NULL if nothing left, or gen_0 in Moebius mode if
    2847             :  * we aren't squarefree */
    2848             : static GEN
    2849       14213 : ifac_main(GEN *partial)
    2850             : {
    2851       14213 :   const long moebius_mode = !!MOEBIUS(*partial);
    2852       14213 :   GEN here = ifac_find(*partial);
    2853             :   long nf;
    2854             : 
    2855       14213 :   if (!here) return NULL; /* nothing left */
    2856             :   /* loop until first entry is a finished prime.  May involve reallocations,
    2857             :    * thus updates of *partial */
    2858       25483 :   while (CLASS(here) != gen_2)
    2859             :   {
    2860       15636 :     if (CLASS(here) == gen_0) /* composite: crack it */
    2861             :     { /* make sure there's room for another factor */
    2862        5789 :       if (here < *partial + 6)
    2863             :       {
    2864           0 :         ifac_defrag(partial, &here);
    2865           0 :         if (here < *partial + 6) ifac_realloc(partial, &here, 1); /* no luck */
    2866             :       }
    2867        5789 :       nf = ifac_crack(partial, &here, moebius_mode);
    2868        5789 :       if (moebius_mode && EXPON(here) != gen_1) /* that was a power */
    2869             :       {
    2870           2 :         if (DEBUGLEVEL >= 3)
    2871           0 :           err_printf("IFAC: main loop: repeated new factor\n\t%Ps\n", *here);
    2872           2 :         return gen_0;
    2873             :       }
    2874             :       /* deal with the new unknowns.  No sort: ifac_crack did it */
    2875        5787 :       ifac_whoiswho(partial, &here, nf);
    2876        5787 :       continue;
    2877             :     }
    2878        9847 :     if (CLASS(here) == gen_1) /* prime but not yet finished: finish it */
    2879             :     {
    2880        9847 :       if (ifac_divide(partial, &here, moebius_mode))
    2881             :       {
    2882           7 :         if (moebius_mode)
    2883             :         {
    2884           0 :           if (DEBUGLEVEL >= 3)
    2885           0 :             err_printf("IFAC: main loop: another factor was divisible by\n"
    2886             :                        "\t%Ps\n", *here);
    2887           0 :           return gen_0;
    2888             :         }
    2889           7 :         ifac_resort(partial, &here); /* sort new cofactors down */
    2890           7 :         ifac_whoiswho(partial, &here, -1);
    2891             :       }
    2892        9847 :       continue;
    2893             :     }
    2894           0 :     pari_err_BUG("ifac_main [nonexistent factor class]");
    2895             :   } /* while */
    2896        9847 :   if (moebius_mode && EXPON(here) != gen_1)
    2897             :   {
    2898           0 :     if (DEBUGLEVEL >= 3)
    2899           0 :       err_printf("IFAC: after main loop: repeated old factor\n\t%Ps\n", *here);
    2900           0 :     return gen_0;
    2901             :   }
    2902        9847 :   if (DEBUGLEVEL >= 4)
    2903             :   {
    2904           0 :     nf = (*partial + lg(*partial) - here - 3)/3;
    2905           0 :     if (nf)
    2906           0 :       err_printf("IFAC: main loop: %ld factor%s left\n", nf, (nf>1)? "s": "");
    2907             :     else
    2908           0 :       err_printf("IFAC: main loop: this was the last factor\n");
    2909             :   }
    2910        9847 :   if (factor_add_primes && !(get_hint(partial) & 8))
    2911             :   {
    2912           0 :     GEN p = VALUE(here);
    2913           0 :     if (lgefint(p)>3 || uel(p,2) > 0x1000000UL) (void)addprimes(p);
    2914             :   }
    2915        9847 :   return here;
    2916             : }
    2917             : 
    2918             : /* Encapsulated routines */
    2919             : 
    2920             : /* prime/exponent pairs need to appear contiguously on the stack, but we also
    2921             :  * need our data structure somewhere, and we don't know in advance how many
    2922             :  * primes will turn up.  The following discipline achieves this:  When
    2923             :  * ifac_decomp() is called, n should point at an object older than the oldest
    2924             :  * small prime/exponent pair  (ifactor() guarantees this).
    2925             :  * We allocate sufficient space to accommodate several pairs -- eleven pairs
    2926             :  * ought to fit in a space not much larger than n itself -- before calling
    2927             :  * ifac_start().  If we manage to complete the factorization before we run out
    2928             :  * of space, we free the data structure and cull the excess reserved space
    2929             :  * before returning.  When we do run out, we have to leapfrog to generate more
    2930             :  * (guesstimating the requirements from what is left in the partial
    2931             :  * factorization structure);  room for fresh pairs is allocated at the head of
    2932             :  * the stack, followed by an ifac_realloc() to reconnect the data structure and
    2933             :  * move it out of the way, followed by a few pointer tweaks to connect the new
    2934             :  * pairs space to the old one. This whole affair translates into a surprisingly
    2935             :  * compact routine. */
    2936             : 
    2937             : /* find primary factors of n; destroy n */
    2938             : static long
    2939        2568 : ifac_decomp(GEN n, long hint)
    2940             : {
    2941        2568 :   pari_sp av = avma;
    2942        2568 :   long nb = 0;
    2943        2568 :   GEN part, here, workspc, pairs = (GEN)av;
    2944             : 
    2945             :   /* workspc will be doled out in pairs of smaller t_INTs. For n = prod p^{e_p}
    2946             :    * (p not necessarily prime), need room to store all p and e_p [ cgeti(3) ],
    2947             :    * bounded by
    2948             :    *    sum_{p | n} ( log_{2^BIL} (p) + 6 ) <= log_{2^BIL} n + 6 log_2 n */
    2949        2568 :   workspc = new_chunk((expi(n) + 1) * 7);
    2950        2568 :   part = ifac_start_hint(n, 0, hint);
    2951             :   for (;;)
    2952             :   {
    2953        7562 :     here = ifac_main(&part);
    2954        7562 :     if (!here) break;
    2955        4994 :     if (gc_needed(av,1))
    2956             :     {
    2957             :       long offset;
    2958           0 :       if(DEBUGMEM>1)
    2959             :       {
    2960           0 :         pari_warn(warnmem,"[2] ifac_decomp");
    2961           0 :         ifac_print(part, here);
    2962             :       }
    2963           0 :       ifac_realloc(&part, &here, 0);
    2964           0 :       offset = here - part;
    2965           0 :       part = gerepileupto((pari_sp)workspc, part);
    2966           0 :       here = part + offset;
    2967             :     }
    2968        4994 :     nb++;
    2969        4994 :     pairs = icopy_avma(VALUE(here), (pari_sp)pairs);
    2970        4994 :     pairs = icopy_avma(EXPON(here), (pari_sp)pairs);
    2971        4994 :     ifac_delete(here);
    2972             :   }
    2973        2568 :   set_avma((pari_sp)pairs);
    2974        2568 :   if (DEBUGLEVEL >= 3)
    2975           0 :     err_printf("IFAC: found %ld large prime (power) factor%s.\n",
    2976             :                nb, (nb>1? "s": ""));
    2977        2568 :   return nb;
    2978             : }
    2979             : 
    2980             : /***********************************************************************/
    2981             : /**            ARITHMETIC FUNCTIONS WITH EARLY-ABORT                  **/
    2982             : /**  needing direct access to the factoring machinery to avoid work:  **/
    2983             : /**  e.g. if we find a square factor, moebius returns 0, core doesn't **/
    2984             : /**  need to factor it, etc.                                          **/
    2985             : /***********************************************************************/
    2986             : /* memory management */
    2987             : static void
    2988           0 : ifac_GC(pari_sp av, GEN *part)
    2989             : {
    2990           0 :   GEN here = NULL;
    2991           0 :   if(DEBUGMEM>1) pari_warn(warnmem,"ifac_xxx");
    2992           0 :   ifac_realloc(part, &here, 0);
    2993           0 :   *part = gerepileupto(av, *part);
    2994           0 : }
    2995             : 
    2996             : /* destroys n */
    2997             : static long
    2998         236 : ifac_moebius(GEN n)
    2999             : {
    3000         236 :   long mu = 1;
    3001         236 :   pari_sp av = avma;
    3002         236 :   GEN part = ifac_start(n, 1);
    3003             :   for(;;)
    3004         468 :   {
    3005             :     long v;
    3006             :     GEN p;
    3007         704 :     if (!ifac_next(&part,&p,&v)) return v? 0: mu;
    3008         468 :     mu = -mu;
    3009         468 :     if (gc_needed(av,1)) ifac_GC(av,&part);
    3010             :   }
    3011             : }
    3012             : 
    3013             : int
    3014         760 : ifac_read(GEN part, GEN *p, long *e)
    3015             : {
    3016         760 :   GEN here = ifac_find(part);
    3017         760 :   if (!here) return 0;
    3018         400 :   *p = VALUE(here);
    3019         400 :   *e = EXPON(here)[2];
    3020         400 :   return 1;
    3021             : }
    3022             : void
    3023         320 : ifac_skip(GEN part)
    3024             : {
    3025         320 :   GEN here = ifac_find(part);
    3026         320 :   if (here) ifac_delete(here);
    3027         320 : }
    3028             : 
    3029             : /* destroys n */
    3030             : static int
    3031           7 : ifac_ispowerful(GEN n)
    3032             : {
    3033           7 :   pari_sp av = avma;
    3034           7 :   GEN part = ifac_start(n, 0);
    3035             :   for(;;)
    3036           7 :   {
    3037             :     long e;
    3038             :     GEN p;
    3039          14 :     if (!ifac_read(part,&p,&e)) return 1;
    3040             :     /* power: skip */
    3041           7 :     if (e != 1 || Z_isanypower(p,NULL)) { ifac_skip(part); continue; }
    3042           0 :     if (!ifac_next(&part,&p,&e)) return 1;
    3043           0 :     if (e == 1) return 0;
    3044           0 :     if (gc_needed(av,1)) ifac_GC(av,&part);
    3045             :   }
    3046             : }
    3047             : /* destroys n; assume n != 0 is composite */
    3048             : static GEN
    3049         353 : ifac_core(GEN n)
    3050             : {
    3051         353 :   GEN m = gen_1, c = cgeti(lgefint(n));
    3052         353 :   pari_sp av = avma;
    3053         353 :   GEN part = ifac_start(n, 0);
    3054             :   for(;;)
    3055         393 :   {
    3056             :     long e;
    3057             :     GEN p;
    3058         746 :     if (!ifac_read(part,&p,&e)) return m;
    3059             :     /* square: skip */
    3060         393 :     if (!odd(e) || Z_issquare(p)) { ifac_skip(part); continue; }
    3061          80 :     if (!ifac_next(&part,&p,&e)) return m;
    3062          80 :     if (odd(e)) m = mulii(m, p);
    3063          80 :     if (gc_needed(av,1)) { affii(m,c); m=c; ifac_GC(av,&part); }
    3064             :   }
    3065             : }
    3066             : 
    3067             : /* must be >= 661 (various functions assume it in order to call uisprime_661
    3068             :  * instead of uisprime, and Z_isanypower_nosmalldiv instead of Z_isanypower) */
    3069             : ulong
    3070     4369152 : tridiv_boundu(ulong n)
    3071             : {
    3072     4369152 :   long e = expu(n);
    3073     4369017 :   if(e<30) return 1UL<<12;
    3074             : #ifdef LONG_IS_64BIT
    3075      182925 :   if(e<34) return 1UL<<13;
    3076      100188 :   if(e<37) return 1UL<<14;
    3077       60054 :   if(e<42) return 1UL<<15;
    3078       26547 :   if(e<47) return 1UL<<16;
    3079       14882 :   if(e<56) return 1UL<<17;
    3080        4589 :   if(e<56) return 1UL<<18;
    3081        4589 :   if(e<62) return 1UL<<19;
    3082        1220 :   return 1UL<<18;
    3083             : #else
    3084        6945 :   return 1UL<<13;
    3085             : #endif
    3086             : }
    3087             : 
    3088             : /* Where to stop trial dividing in factorization. Must be >= 661.
    3089             :  * If further n > 2^512, must be >= 2^14 */
    3090             : ulong
    3091      869536 : tridiv_bound(GEN n)
    3092             : {
    3093      869536 :   if (lgefint(n)==3) return tridiv_boundu(n[2]);
    3094             :   else
    3095             :   {
    3096       76318 :     ulong l = (ulong)expi(n) + 1;
    3097       76318 :     if (l <= 512) return (l-16) << 10;
    3098        1196 :     return 1UL<<19; /* Rho is generally faster above this */
    3099             :   }
    3100             : }
    3101             : 
    3102             : /* destroys n */
    3103             : static void
    3104         863 : ifac_factoru(GEN n, long hint, GEN P, GEN E, long *pi)
    3105             : {
    3106         863 :   GEN part = ifac_start_hint(n, 0, hint);
    3107         863 :   long i = *pi;
    3108             :   for(;;)
    3109        1634 :   {
    3110             :     long v;
    3111             :     GEN p;
    3112        2497 :     if (!ifac_next(&part,&p,&v)) { *pi = i; return; }
    3113        1634 :     P[i] = itou(p); E[i] = v; i++;
    3114             :   }
    3115             : }
    3116             : /* destroys n */
    3117             : static long
    3118         663 : ifac_moebiusu(GEN n)
    3119             : {
    3120         663 :   GEN part = ifac_start(n, 1);
    3121         663 :   long s = 1;
    3122             :   for(;;)
    3123        1326 :   {
    3124             :     long v;
    3125             :     GEN p;
    3126        1989 :     if (!ifac_next(&part,&p,&v)) return v? 0: s;
    3127        1326 :     s = -s;
    3128             :   }
    3129             : }
    3130             : 
    3131             : INLINE ulong
    3132   143807561 : u_forprime_next_fast(forprime_t *T)
    3133             : {
    3134   143807561 :   if (++T->n <= pari_PRIMES[0])
    3135             :   {
    3136   143808292 :     T->p = pari_PRIMES[T->n];
    3137   143808292 :     return T->p > T->b ? 0: T->p;
    3138             :   }
    3139           0 :   return u_forprime_next(T);
    3140             : }
    3141             : 
    3142             : /* uisprime(n) knowing n has no prime divisor <= lim */
    3143             : static int
    3144        7866 : uisprime_nosmall(ulong n, ulong lim)
    3145        7866 : { return (lim >= 661)? uisprime_661(n): uisprime(n); }
    3146             : 
    3147             : static GEN factoru_sign(ulong n, ulong all, long hint, ulong *pU1, ulong *pU2);
    3148             : static GEN ifactor_sign(GEN n, ulong all, long hint, long sn, GEN *pU);
    3149             : 
    3150             : /* simplified version of factoru_sign, to be called on squarefree n whose
    3151             :  * prime divisors are in [minp, maxp]. In practice called with
    3152             :  * maxp <= maxprimelim() */
    3153             : static GEN
    3154      917212 : factoru_primes(ulong n, ulong minp, ulong maxp)
    3155             : {
    3156             :   forprime_t S;
    3157             :   ulong p;
    3158             :   long i;
    3159             :   GEN P;
    3160             : 
    3161      917212 :   if (n < minp) return NULL;
    3162      898337 :   if (n <= maxp && PRIMES_search(n) > 0) return mkvecsmall(n);
    3163      761722 :   P = cgetg(16, t_VECSMALL); i = 1;
    3164      761722 :   u_forprime_init(&S, minp, maxp);
    3165    39059687 :   while ( (p = u_forprime_next_fast(&S)) )
    3166             :   {
    3167    39059687 :     ulong q = n / p;
    3168    39059687 :     if (n % p == 0)
    3169             :     {
    3170     1246799 :       P[i++] = p; n = q;
    3171     1246799 :       if (q <= p || (n <= maxp && PRIMES_search(n) > 0)) { P[i++] = n; break; }
    3172             :     }
    3173    37812888 :     else if (q <= p) { P[i++] = n; break; } /* n <= p^2: n is now prime */
    3174             :   }
    3175      761722 :   if (i == 1) return NULL;
    3176      761722 :   setlg(P, i); return P;
    3177             : }
    3178             : static GEN
    3179      127536 : Z_factor_primes(GEN N, ulong minp, ulong maxp)
    3180             : {
    3181             :   forprime_t S;
    3182      127536 :   ulong p, n = 0;
    3183             :   long i;
    3184             :   GEN P;
    3185      127536 :   if (lgefint(N) == 3) return factoru_primes(uel(N,2), minp, maxp);
    3186       19086 :   u_forprime_init(&S, minp, maxp);
    3187       19086 :   P = cgetg(expi(N) + 1, t_VECSMALL); i = 1;
    3188    11437046 :   while ( (p = u_forprime_next_fast(&S)) )
    3189             :   {
    3190             :     int stop;
    3191    11437050 :     long v = Z_lvalrem_stop(&N, p, &stop);
    3192    11437046 :     if (v) P[i++] = p;
    3193    11437046 :     if (stop)
    3194             :     {
    3195         320 :       if (!equali1(N)) P[i++] = uel(N,2);
    3196         320 :       goto END;
    3197             :     }
    3198    11436726 :     if (v && lgefint(N) == 3) { n = uel(N,2); break; }
    3199             :   }
    3200    23972052 :   if (n) while ( (p = u_forprime_next_fast(&S)) )
    3201             :   {
    3202    23972052 :     ulong q = n / p;
    3203    23972052 :     if (n % p == 0)
    3204             :     {
    3205       64015 :       P[i++] = p; n = q;
    3206       64015 :       if (q <= p || (n <= maxp && PRIMES_search(n) > 0)) { P[i++] = n; break; }
    3207             :     }
    3208    23908037 :     else if (q <= p) { P[i++] = n; break; } /* n <= p^2: n is now prime */
    3209             :   }
    3210           0 : END:
    3211       19086 :   if (i == 1) return NULL;
    3212       19086 :   setlg(P, i); return P;
    3213             : }
    3214             : 
    3215             : /* N != 0. Product of odd prime divisors less than
    3216             :  *   min(*pLIM, factorlimit) [WARNING!];
    3217             :  * with lim <= *pLIM < 2*lim and *pLIM prime
    3218             :  * Assume lim >= 128. Better for efficiency if N >= lim^2. */
    3219             : static ulong
    3220      809154 : u_oddprimedivisors_gcd(ulong N, ulong lim, ulong *pLIM)
    3221             : {
    3222      809154 :   GEN PR = prodprimes(), LIM = prodprimeslim();
    3223      809154 :   long b = minss(lg(PR)-1, expu(lim)-6);
    3224             :   /* 2^{b+6} <= lim < 2^{b+7}, b >= 1 */
    3225      809154 :   *pLIM = LIM[b]; return ugcd(N, umodiu(gel(PR,b), N));
    3226             : }
    3227             : /* not GC-clean */
    3228             : static GEN
    3229      128213 : Z_oddprimedivisors_gcd(GEN N, ulong lim, ulong *pLIM)
    3230             : {
    3231      128213 :   GEN PR = prodprimes(), LIM = prodprimeslim();
    3232      128213 :   long b = minss(lg(PR)-1, expu(lim)-6);
    3233      128213 :   *pLIM = LIM[b]; return gcdii(N, gel(PR,b));
    3234             : }
    3235             : 
    3236             : /* Assume lim >= 128 and N odd. */
    3237             : static GEN
    3238      126933 : Z_oddprimedivisors_fast(GEN N, ulong lim)
    3239             : {
    3240      126933 :   pari_sp av = avma;
    3241      126933 :   GEN Nr = Z_oddprimedivisors_gcd(N, lim, &lim);
    3242      126933 :   GEN P = Z_factor_primes(Nr, 3, lim);
    3243      126933 :   return P? P: gc_NULL(av);
    3244             : }
    3245             : /* return mask with bit 0, 1, 2 set if respectively 3, 5, 7 divide n */
    3246             : static int
    3247    15807641 : u_357_divides(ulong n)
    3248             : { /* vector (105, i, n = i-1; !(n%3) + 2 * !(n%5) + 4 * !(n%7)) */
    3249    15807641 :   const unsigned int tab[] = {
    3250             :   7,0,0,1,0,2,1,4,0,1,2,0,1,0,4,3,0,0,1,0,2,5,0,0,1,2,0,1,4,0,3,0,0,1,0,6,1,0,
    3251             :   0,1,2,0,5,0,0,3,0,0,1,4,2,1,0,0,1,2,4,1,0,0,3,0,0,5,0,2,1,0,0,1,6,0,1,0,0,3,
    3252             :   0,4,1,0,2,1,0,0,5,2,0,1,0,0,3,4,0,1,0,2,1,0,4,1,2,0,1,0,0};
    3253    15807641 :   return tab[n % 105UL];
    3254             : }
    3255             : 
    3256             : static GEN
    3257    17181630 : factoru_result(GEN P, GEN E, long i)
    3258             : {
    3259    17181630 :   GEN P2, E2, f = cgetg(3,t_VEC);
    3260    17181261 :   gel(f,1) = P2 = cgetg(i, t_VECSMALL);
    3261    17180670 :   gel(f,2) = E2 = cgetg(i, t_VECSMALL);
    3262    52490595 :   while (--i >= 1) { P2[i] = P[i]; E2[i] = E[i]; }
    3263    17180555 :   return f;
    3264             : }
    3265             : 
    3266             : /* Factor n and output [p,e] where
    3267             :  * p, e are vecsmall with n = prod{p[i]^e[i]}. If all != 0:
    3268             :  * if pU1 is not NULL, set *pU1 and *pU2 so that unfactored composite is
    3269             :  * U1^U2 with U1 not a pure power; else include it in factorization */
    3270             : static GEN
    3271    17613518 : factoru_sign(ulong n, ulong all, long hint, ulong *pU1, ulong *pU2)
    3272             : {
    3273    17613518 :   pari_sp av = avma;
    3274    17613518 :   ulong ALL, p, lim = 0;
    3275    17613518 :   long i, oldi = -1;
    3276             :   forprime_t S;
    3277             :   GEN E, P;
    3278             : 
    3279    17613518 :   if (pU1) *pU1 = *pU2 = 1;
    3280    17613518 :   if (n == 0) retmkvec2(mkvecsmall(0), mkvecsmall(1));
    3281    17613518 :   if (n == 1) retmkvec2(cgetg(1,t_VECSMALL), cgetg(1,t_VECSMALL));
    3282             : 
    3283             :   /* enough room to store <= 15 primes and exponents (OK if n < 2^64) */
    3284    17181045 :   (void)new_chunk(3 + 16*2);
    3285    17181009 :   P = cgetg(16, t_VECSMALL); i = 1;
    3286    17180916 :   E = cgetg(16, t_VECSMALL);
    3287    17180551 :   ALL = all? all: ULONG_MAX; /* (!all || all > n) == ALL > n */
    3288    17180551 :   if (ALL > 2)
    3289             :   {
    3290             :     ulong maxp;
    3291    17180517 :     long v = vals(n);
    3292    17180477 :     if (v)
    3293             :     {
    3294    11638672 :       P[1] = 2; E[1] = v; i = 2;
    3295    11638672 :       n >>= v; if (n == 1) goto END;
    3296             :     }
    3297    14053473 :     if (ALL > 3)
    3298             :     {
    3299    14053932 :       int mask = u_357_divides(n);
    3300    14054408 :       if (mask)
    3301             :       {
    3302     9458704 :         if (mask & 1)
    3303     6001299 :         { P[i] = 3; E[i] = 1 + u_lvalrem(n / 3, 3, &n); i++;
    3304     6001355 :           if (n == 1) goto END; }
    3305     7144607 :         if ((mask & 2) && ALL > 5)
    3306     3468839 :         { P[i] = 5; E[i] = 1 + u_lvalrem(n / 5, 5, &n); i++;
    3307     3468855 :           if (n == 1) goto END; }
    3308     5578039 :         if ((mask & 4) && ALL > 7)
    3309     2079364 :         { P[i] = 7; E[i] = 1 + u_lvalrem(n / 7, 7, &n); i++;
    3310     2079371 :           if (n == 1) goto END; }
    3311             :       }
    3312             :     }
    3313     9013178 :     maxp = maxprime();
    3314     9014411 :     if (n <= maxp && PRIMES_search(n) > 0) { P[i] = n; E[i] = 1; i++; goto END; }
    3315     3254154 :     lim = all? all-1: tridiv_boundu(n);
    3316     3254152 :     if (lim >= 128 && n >= 691 * 691) /* expu(lim) >= 7 */
    3317        7777 :     { /* fast trial division */
    3318      804400 :       ulong gcdlim, gcd, sqrtn = usqrt(n);
    3319             :       GEN Q;
    3320      804400 :       lim = minuu(lim, sqrtn);
    3321      804400 :       gcd = u_oddprimedivisors_gcd(n, lim, &gcdlim);
    3322      804400 :       Q = factoru_primes(gcd, 11, gcdlim);
    3323      804400 :       maxp = GP_DATA->factorlimit;
    3324      804400 :       if (Q)
    3325             :       {
    3326      789362 :         long j, l = lg(Q);
    3327      789362 :         int stop = 0;
    3328     2578953 :         for (j = 1; j < l; j++)
    3329             :         {
    3330     1789597 :           ulong p = uel(Q,j);
    3331     1789597 :           if (all && p >= all) { stop = 1; break; }
    3332     1789591 :           E[i] = u_lvalrem_stop(&n, p, &stop); /* > 0 */
    3333     1789591 :           P[i] = p; i++;
    3334             :         }
    3335      825110 :         if (n == 1) goto END;
    3336       39113 :         if (stop || (n <= maxp && PRIMES_search(n) > 0))
    3337       35748 :         { P[i] = n; E[i] = 1; i++; goto END; }
    3338             :       }
    3339       15038 :       else if (lim == sqrtn && lim <= maxp)
    3340       10626 :         { P[i] = n; E[i] = 1; i++; goto END; }
    3341             :     }
    3342             :     else
    3343             :     { /* naive trial division */
    3344     2449752 :       maxp = lim;
    3345     2449752 :       u_forprime_init(&S, 11, lim);
    3346    17695257 :       while ( (p = u_forprime_next_fast(&S)) )
    3347             :       {
    3348             :         int stop;
    3349             :         /* tiny integers without small factors are often primes */
    3350    17694893 :         if (p == 673)
    3351             :         {
    3352     2449566 :           if (uisprime_661(n)) { P[i] = n; E[i] = 1; i++; goto END; }
    3353         190 :           oldi = i;
    3354             :         }
    3355    17694893 :         v = u_lvalrem_stop(&n, p, &stop);
    3356    17694886 :         if (v) { P[i] = p; E[i] = v; i++; }
    3357    17694886 :         if (stop)
    3358             :         {
    3359     2449376 :           if (n != 1) { P[i] = n; E[i] = 1; i++; }
    3360     2449376 :           goto END;
    3361             :         }
    3362             :       }
    3363             :     }
    3364        8168 :     if (lim > maxp)
    3365             :     { /* second pass usually empty, outside fast trial division range */
    3366             :       long v;
    3367           6 :       u_forprime_init(&S, maxp+1, lim);
    3368     5296866 :       while ((p = u_forprime_next(&S)))
    3369             :       {
    3370             :         int stop;
    3371     5296900 :         v = u_lvalrem_stop(&n, p, &stop);
    3372     5296866 :         if (v) { P[i] = p; E[i] = v; i++; }
    3373     5296866 :         if (stop)
    3374             :         {
    3375           6 :           if (n != 1) { P[i] = n; E[i] = 1; i++; }
    3376           6 :           goto END;
    3377             :         }
    3378             :       }
    3379             :     }
    3380             :   }
    3381             :   /* if i > oldi (includes oldi = -1) we don't know that n is composite */
    3382        8162 :   if (all)
    3383             :   { /* smallfact: look for easy pure powers then stop */
    3384             : #ifdef LONG_IS_64BIT
    3385        1080 :     ulong mask = all > 563 ? (all > 7129 ? 1: 3): 7;
    3386             : #else
    3387          15 :     ulong mask = all > 22 ? (all > 83 ? 1: 3): 7;
    3388             : #endif
    3389        1095 :     long k = 1, ex;
    3390        1588 :     while (uissquareall(n, &n)) k <<= 1;
    3391        1108 :     while ( (ex = uis_357_power(n, &n, &mask)) ) k *= ex;
    3392        1095 :     if (pU1 && (i == oldi || !uisprime_nosmall(n, lim)))
    3393         266 :     { *pU1 = n; *pU2 = (ulong)k; }
    3394             :     else
    3395         829 :     { P[i] = n; E[i] = k; i++; }
    3396        1095 :     goto END;
    3397             :   }
    3398             :   /* we don't known that n is composite ? */
    3399        7067 :   if (oldi != i && uisprime_nosmall(n, lim)) { P[i]=n; E[i]=1; i++; goto END; }
    3400             : 
    3401             :   {
    3402             :     GEN perm;
    3403         863 :     ifac_factoru(utoipos(n), hint, P, E, &i);
    3404         863 :     setlg(P, i);
    3405         863 :     perm = vecsmall_indexsort(P);
    3406         863 :     P = vecsmallpermute(P, perm);
    3407         863 :     E = vecsmallpermute(E, perm);
    3408             :   }
    3409    17182456 : END:
    3410    17182456 :   set_avma(av); return factoru_result(P, E, i);
    3411             : }
    3412             : GEN
    3413     3741503 : factoru(ulong n)
    3414     3741503 : { return factoru_sign(n, 0, decomp_default_hint, NULL, NULL); }
    3415             : 
    3416             : ulong
    3417           0 : radicalu(ulong n)
    3418             : {
    3419           0 :   pari_sp av = avma;
    3420           0 :   return gc_long(av, zv_prod(gel(factoru(n),1)));
    3421             : }
    3422             : 
    3423             : long
    3424       54194 : moebiusu_fact(GEN f)
    3425             : {
    3426       54194 :   GEN E = gel(f,2);
    3427       54194 :   long i, l = lg(E);
    3428       93569 :   for (i = 1; i < l; i++)
    3429       57834 :     if (E[i] > 1) return 0;
    3430       35735 :   return odd(l)? 1: -1;
    3431             : }
    3432             : 
    3433             : long
    3434     2487299 : moebiusu(ulong n)
    3435             : {
    3436             :   pari_sp av;
    3437             :   long s, v, test_prime;
    3438             :   ulong p, lim;
    3439             :   int mask;
    3440             : 
    3441     2487299 :   switch(n)
    3442             :   {
    3443           0 :     case 0: (void)check_arith_non0(gen_0,"moebius");/*error*/
    3444      568237 :     case 1: return  1;
    3445      106667 :     case 2: return -1;
    3446             :   }
    3447             :   /* n > 2 */
    3448     1832112 :   p = n & 3; if (!p) return 0;
    3449     1743874 :   if (p == 2) { n >>= 1; s = -1; } else s = 1;
    3450     1743874 :   mask = u_357_divides(n);
    3451     1762802 :   if (mask)
    3452             :   {
    3453      677844 :     if (mask & 1) { n /= 3; s = -s; if (n % 3 == 0) return 0; }
    3454      613916 :     if (mask & 2) { n /= 5; s = -s; if (n % 5 == 0) return 0; }
    3455      572553 :     if (mask & 4) { n /= 7; s = -s; if (n % 7 == 0) return 0; }
    3456             :   }
    3457     1622366 :   if (n <= maxprimelim() && PRIMES_search(n) > 0) return -s;
    3458      661399 :   av = avma; lim = tridiv_boundu(n);
    3459      672117 :   if (n >= 691 * 691)
    3460             :   {
    3461        4507 :     ulong gcdlim, gcd, sqrtn = usqrt(n);
    3462             :     GEN P;
    3463        4507 :     lim = minuu(sqrtn, lim);
    3464        4507 :     gcd = u_oddprimedivisors_gcd(n, lim, &gcdlim);
    3465        4507 :     if (gcd != 1)
    3466             :     {
    3467        3523 :       n /= gcd;
    3468        3846 :       if (ugcd(gcd, n) != 1) return 0;
    3469             :     }
    3470        4362 :     P = factoru_primes(gcd, 11, gcdlim);
    3471        4362 :     if (P && odd(lg(P) - 1)) s = -s;
    3472        4362 :     if (n == 1) return gc_long(av, s);
    3473        4060 :     if (lim == sqrtn && lim <= GP_DATA->factorlimit) return gc_long(av, -s);
    3474        4039 :     test_prime = 1;
    3475             :   }
    3476             :   else
    3477             :   {
    3478             :     forprime_t S;
    3479      667610 :     u_forprime_init(&S, 3, lim);
    3480      666022 :     test_prime = 0;
    3481     7417388 :     while ((p = u_forprime_next_fast(&S)))
    3482             :     {
    3483             :       int stop;
    3484             :       /* tiny integers without small factors are often primes */
    3485     7416328 :       if (p == 673)
    3486             :       {
    3487           0 :         test_prime = 0;
    3488      665023 :         if (uisprime_661(n)) return gc_long(av,-s);
    3489             :       }
    3490     7416328 :       v = u_lvalrem_stop(&n, p, &stop);
    3491     7415431 :       if (v) {
    3492      635317 :         if (v > 1) return gc_long(av,0);
    3493      593329 :         test_prime = 1;
    3494      593329 :         s = -s;
    3495             :       }
    3496     7373443 :       if (stop) return gc_long(av, n==1? s: -s);
    3497             :     }
    3498             :   }
    3499        4039 :   set_avma(av);
    3500        4039 :   if (test_prime && uisprime_661(n)) return -s;
    3501             :   else
    3502             :   {
    3503         663 :     long t = ifac_moebiusu(utoipos(n));
    3504         663 :     set_avma(av);
    3505         663 :     if (t == 0) return 0;
    3506         663 :     return (s == t)? 1: -1;
    3507             :   }
    3508             : }
    3509             : 
    3510             : long
    3511       58815 : moebius(GEN n)
    3512             : {
    3513       58815 :   pari_sp av = avma;
    3514             :   GEN F;
    3515             :   ulong p, lim, n357;
    3516             :   long i, l, s, v, copy;
    3517             :   int mask;
    3518             : 
    3519       58815 :   if ((F = check_arith_non0(n,"moebius")))
    3520             :   {
    3521             :     GEN E;
    3522         728 :     F = clean_Z_factor(F);
    3523         728 :     E = gel(F,2);
    3524         728 :     l = lg(E);
    3525        1428 :     for(i = 1; i < l; i++)
    3526         980 :       if (!equali1(gel(E,i))) return gc_long(av,0);
    3527         448 :     return gc_long(av, odd(l)? 1: -1);
    3528             :   }
    3529       58050 :   if (lgefint(n) == 3) return moebiusu(uel(n,2));
    3530        1608 :   p = mod4(n); if (!p) return 0;
    3531        1408 :   copy = s = 1;
    3532        1408 :   if (p == 2)
    3533             :   {
    3534         358 :     n = shifti(n, -1);
    3535         358 :     copy = 0; s = -1;
    3536             :   }
    3537        1408 :   n357 = umodiu(n, 9 * 25 * 49);
    3538        1408 :   mask = u_357_divides(n357);
    3539        1408 :   if (mask)
    3540             :   {
    3541         764 :     ulong m = 1;
    3542         764 :     if (mask & 1) { m = 3;  s = -s; }
    3543         764 :     if (mask & 2) { m *= 5; s = -s; }
    3544         764 :     if (mask & 4) { m *= 7; s = -s; }
    3545         764 :     if (u_357_divides(n357 / m)) return gc_long(av, 0);
    3546         549 :     copy = 0; n = diviuexact(n, m);
    3547             :   }
    3548        1193 :   if (copy) n = icopy(n);
    3549         701 :   else if (lgefint(n) == 3) return gc_long(av, s * moebiusu(uel(n,2)));
    3550         883 :   setabssign(n); lim = tridiv_bound(n);
    3551         883 :   if (lim >= 128)
    3552             :   {
    3553             :     ulong gcdlim;
    3554         883 :     GEN gcd = Z_oddprimedivisors_gcd(n, lim, &gcdlim);
    3555         883 :     if (!equali1(gcd))
    3556             :     {
    3557             :       GEN P;
    3558         622 :       n = diviiexact(n, gcd);
    3559         808 :       if (!equali1(gcdii(gcd, n))) return gc_long(av, 0);
    3560         603 :       P = Z_factor_primes(gcd, 11, gcdlim);
    3561         603 :       if (P)
    3562             :       {
    3563         603 :         if (odd(lg(P) - 1)) s = -s;
    3564         603 :         if (is_pm1(n)) return gc_long(av, s);
    3565        1158 :         if (lim <= GP_DATA->factorlimit &&
    3566         741 :             cmpii(sqru(lim), n) >= 0) return gc_long(av, -s); /* n prime */
    3567             :       }
    3568             :     }
    3569             :   }
    3570             :   else
    3571             :   {
    3572             :     forprime_t S;
    3573           0 :     u_forprime_init(&S, 3, lim);
    3574           0 :     while ((p = u_forprime_next_fast(&S)))
    3575             :     {
    3576             :       int stop;
    3577           0 :       v = Z_lvalrem_stop(&n, p, &stop);
    3578           0 :       if (v)
    3579             :       {
    3580           0 :         if (v > 1) return gc_long(av,0);
    3581           0 :         s = -s;
    3582             :       }
    3583           0 :       if (stop) return gc_long(av, is_pm1(n)? s: -s);
    3584             :     }
    3585             :   }
    3586         678 :   l = lg(primetab);
    3587         682 :   for (i = 1; i < l; i++)
    3588             :   {
    3589           7 :     v = Z_pvalrem(n, gel(primetab,i), &n);
    3590           7 :     if (v)
    3591             :     {
    3592           7 :       if (v > 1) return gc_long(av,0);
    3593           5 :       s = -s;
    3594           5 :       if (is_pm1(n)) return gc_long(av,s);
    3595             :     }
    3596             :   }
    3597         675 :   if (ifac_isprime(n)) return gc_long(av,-s);
    3598             :   /* large composite without small factors */
    3599         236 :   v = ifac_moebius(n);
    3600         236 :   return gc_long(av, s < 0? -v: v); /* correct also if v==0 */
    3601             : }
    3602             : 
    3603             : long
    3604        1708 : ispowerful(GEN n)
    3605             : {
    3606        1708 :   pari_sp av = avma;
    3607             :   GEN F;
    3608             :   ulong p, lim, n357;
    3609             :   long i, l, v;
    3610        1708 :   int mask, copy = 1;
    3611             : 
    3612        1708 :   if ((F = check_arith_all(n, "ispowerful")))
    3613             :   {
    3614         742 :     GEN p, P = gel(F,1), E = gel(F,2);
    3615         742 :     if (lg(P) == 1) return 1; /* 1 */
    3616         728 :     p = gel(P,1);
    3617         728 :     if (!signe(p)) return 1; /* 0 */
    3618         707 :     i = is_pm1(p)? 2: 1; /* skip -1 */
    3619         707 :     l = lg(E);
    3620         980 :     for (; i < l; i++)
    3621         847 :       if (equali1(gel(E,i))) return 0;
    3622         133 :     return 1;
    3623             :   }
    3624         966 :   if (!signe(n)) return 1;
    3625         952 :   if (mod4(n) == 2) return 0;
    3626         623 :   n357 = umodiu(n, 9 * 25 * 49);
    3627         623 :   mask = u_357_divides(n357);
    3628         623 :   if (mask)
    3629             :   {
    3630         315 :     if ((mask & 1) && n357 % 9)  return 0;
    3631         203 :     if ((mask & 2) && n357 % 25) return 0;
    3632         126 :     if ((mask & 4) && n357 % 49) return 0;
    3633          98 :     if (mask & 1) (void)Z_lvalrem(diviuexact(n,9),  3, &n);
    3634          98 :     if (mask & 2) (void)Z_lvalrem(diviuexact(n,25), 5, &n);
    3635          98 :     if (mask & 4) (void)Z_lvalrem(diviuexact(n,49), 7, &n);
    3636          98 :     copy = 0;
    3637             :   }
    3638         406 :   if (!mod2(n)) { n = shifti(n, -vali(n)); copy = 0; }
    3639         406 :   if (is_pm1(n)) return gc_long(av, 1);
    3640         238 :   if (copy) n = icopy(n);
    3641         238 :   setabssign(n); lim = tridiv_bound(n);
    3642         238 :   if (cmpiu(n, 691 * 691) >= 0)
    3643             :   {
    3644          70 :     ulong gcdlim, sqrtn = 0;
    3645             :     GEN gcd;
    3646          70 :     if (lgefint(n) == 3)
    3647             :     {
    3648           6 :       sqrtn = usqrt(n[2]);
    3649           6 :       lim = minuu(sqrtn, lim);
    3650             :     }
    3651          70 :     gcd = Z_oddprimedivisors_gcd(n, lim, &gcdlim);
    3652          70 :     if (!equali1(gcd))
    3653             :     {
    3654             :       GEN r;
    3655          70 :       n = diviiexact(n, gcd);
    3656          70 :       n = dvmdii(n, gcd, &r);
    3657          70 :       if (r != gen_0) return gc_long(av, 0);
    3658          70 :       n = Z_ppo(n, gcd);
    3659             :     }
    3660             :     /* prime divisors > gcdlim */
    3661          70 :     if (equali1(n)) return gc_long(av, 1);
    3662           7 :     if (sqrtn && gcdlim >= sqrtn) return gc_long(av, 0); /* prime */
    3663             :   }
    3664             :   else
    3665             :   {
    3666             :     forprime_t S;
    3667         168 :     u_forprime_init(&S, 3, lim);
    3668         378 :     while ((p = u_forprime_next_fast(&S)))
    3669             :     {
    3670             :       int stop;
    3671         378 :       v = Z_lvalrem_stop(&n, p, &stop);
    3672         546 :       if (v == 1) return gc_long(av,0);
    3673         378 :       if (stop) return gc_long(av, is_pm1(n)); /* n > 1 is now prime */
    3674             :     }
    3675             :   }
    3676           7 :   l = lg(primetab);
    3677           7 :   for (i = 1; i < l; i++)
    3678             :   {
    3679           0 :     v = Z_pvalrem(n, gel(primetab,i), &n);
    3680           0 :     if (v)
    3681             :     {
    3682           0 :       if (v == 1) return gc_long(av,0);
    3683           0 :       if (is_pm1(n)) return gc_long(av,1);
    3684             :     }
    3685             :   }
    3686             :   /* no need to factor: must be p^2 or not powerful */
    3687           7 :   if (cmpii(powuu(lim+1, 3), n) > 0) return gc_long(av,  Z_issquare(n));
    3688             : 
    3689           7 :   if (ifac_isprime(n)) return gc_long(av,0);
    3690             :   /* large composite without small factors */
    3691           7 :   return gc_long(av, ifac_ispowerful(n));
    3692             : }
    3693             : 
    3694             : ulong
    3695           0 : coreu_fact(GEN f)
    3696             : {
    3697           0 :   GEN P = gel(f,1), E = gel(f,2);
    3698           0 :   long i, l = lg(P), m = 1;
    3699           0 :   for (i = 1; i < l; i++)
    3700             :   {
    3701           0 :     ulong p = P[i], e = E[i];
    3702           0 :     if (e & 1) m *= p;
    3703             :   }
    3704           0 :   return m;
    3705             : }
    3706             : 
    3707             : /* d = a squarefree divisor of n. Return n / (n, d^oo)
    3708             :  * and set *pcore = \prod_{p | (n,d), v_p(n) odd} p
    3709             :  * Simpified form of Z_cba. */
    3710             : static GEN
    3711         327 : core_init_from_squarefree(GEN n, GEN d, GEN *pcore)
    3712             : {
    3713         327 :   GEN c = gen_1;
    3714             :   long v;
    3715             : 
    3716         327 :   if (equali1(d)) { *pcore = c; return n; }
    3717         260 :   v = Z_pvalrem(n, d, &n);
    3718          59 :   for (;; v++)
    3719          59 :   { /* d^v divides "original n" */
    3720         319 :     GEN newd = gcdii(n, d); /* newd^{v+1} divides original n */
    3721         319 :     if (!equalii(d, newd))
    3722             :     { /* new d loses primes dividing original n to exact power v */
    3723         310 :       if (odd(v)) c = mulii(c, diviiexact(d, newd)); /* lost primes */
    3724         310 :       d = newd; if (equali1(d)) break;
    3725             :     }
    3726          63 :     if (equalii(d, n))
    3727             :     {
    3728           4 :       if (odd(v + 1)) c = mulii(c, d);
    3729           4 :       *pcore = c; return gen_1;
    3730             :     }
    3731          59 :     n = diviiexact(n, d);
    3732             :   }
    3733         256 :   *pcore = c; return n;
    3734             : }
    3735             : static ulong
    3736         247 : coreu_init_from_squarefree(ulong n, ulong d, ulong *pcore)
    3737             : {
    3738         247 :   ulong c = 1;
    3739             :   long v;
    3740             : 
    3741         247 :   if (d == 1) { *pcore = c; return n; }
    3742         205 :   v = u_lvalrem(n, d, &n);
    3743          18 :   for (;; v++)
    3744          18 :   { /* d^v divides "original n" */
    3745         223 :     ulong newd = ugcd(n, d); /* newd^{v+1} divides original n */
    3746         223 :     if (d != newd)
    3747             :     { /* new d loses primes dividing original n to exact power v */
    3748         211 :       if (odd(v)) c *= d / newd; /* lost primes */
    3749         211 :       d = newd; if (d == 1) break;
    3750             :     }
    3751          90 :     if (d == n)
    3752             :     {
    3753          72 :       if (odd(v + 1)) c *= d;
    3754          72 :       *pcore = c; return 1;
    3755             :     }
    3756          18 :     n /= d;
    3757             :   }
    3758         133 :   *pcore = c; return n;
    3759             : }
    3760             : 
    3761             : ulong
    3762       62112 : coreu(ulong n)
    3763             : {
    3764             :   ulong p, lim, m;
    3765             :   long v;
    3766             : 
    3767       62112 :   if (!n) return 0;
    3768       62112 :   m = 1;
    3769       62112 :   v = vals(n);
    3770       62112 :   if (v)
    3771             :   {
    3772        4133 :     n >>= v;
    3773        4133 :     if (odd(v)) m = 2;
    3774             :   }
    3775       62112 :   v = u_lvalrem(n, 3, &n); if (odd(v)) m *= 3;
    3776       62112 :   v = u_lvalrem(n, 5, &n); if (odd(v)) m *= 5;
    3777       62112 :   v = u_lvalrem(n, 7, &n); if (odd(v)) m *= 7;
    3778       62112 :   if (n == 1) return m;
    3779        3320 :   if (n <= maxprimelim() && PRIMES_search(n) > 0) return m * n;
    3780         408 :   lim = tridiv_boundu(n);
    3781         408 :   if (n >= 691 * 691)
    3782             :   {
    3783         247 :     ulong mpart, gcd, gcdlim, sqrtn = usqrt(n);
    3784         247 :     lim = minuu(sqrtn, lim);
    3785         247 :     gcd = u_oddprimedivisors_gcd(n, lim, &gcdlim);
    3786         247 :     n = coreu_init_from_squarefree(n, gcd, &mpart);
    3787         247 :     m *= mpart;
    3788         266 :     if (n == 1) return m;
    3789             :     /* n has no prime divisor <= gcdlim */
    3790         103 :     if ((lim == sqrtn && lim <= GP_DATA->factorlimit)
    3791          96 :         || (gcdlim + 1) * (gcdlim + 1) > n)
    3792          19 :       return m * n; /* prime */
    3793             :   }
    3794             :   else
    3795             :   {
    3796             :     forprime_t S;
    3797         161 :     u_forprime_init(&S, 11, lim);
    3798        3633 :     while ((p = u_forprime_next_fast(&S)))
    3799             :     {
    3800             :       int stop;
    3801        3633 :       v = u_lvalrem_stop(&n, p, &stop);
    3802        3633 :       if (v & 1) m *= p;
    3803        3633 :       if (stop) return n == 1? m: m * n; /* n > 1 is now prime */
    3804             :     }
    3805             :   }
    3806          84 :   if (uisprime_661(n)) return m * n;
    3807             :   else
    3808             :   {
    3809          84 :     pari_sp av = avma;
    3810          84 :     m *= itou(ifac_core(utoipos(n)));
    3811          84 :     return gc_ulong(av, m);
    3812             :   }
    3813             : }
    3814             : 
    3815             : GEN
    3816      708659 : core(GEN n)
    3817             : {
    3818      708659 :   pari_sp av = avma;
    3819             :   GEN m, mpart, gcd, F;
    3820             :   ulong lim, gcdlim, mask, m0;
    3821             :   long i, l, v, s;
    3822      708659 :   int copy = 1;
    3823             : 
    3824      708659 :   if ((F = check_arith_all(n, "core")))
    3825             :   {
    3826      646236 :     GEN p, x, P = gel(F,1), E = gel(F,2);
    3827      646236 :     long j = 1;
    3828      646236 :     if (lg(P) == 1) return gen_1;
    3829      646208 :     p = gel(P,1);
    3830      646208 :     if (!signe(p)) return gen_0;
    3831      646166 :     l = lg(P); x = cgetg(l, t_VEC);
    3832     2283030 :     for (i = 1; i < l; i++)
    3833     1636865 :       if (mpodd(gel(E,i))) gel(x,j++) = gel(P,i);
    3834      646165 :     setlg(x, j); return ZV_prod(x);
    3835             :   }
    3836       62410 :   s = signe(n);
    3837       62410 :   if (!s) return gen_0;
    3838       62382 :   if (lgefint(n) == 3)
    3839             :   {
    3840       61975 :     ulong c = coreu(uel(n,2));
    3841       61975 :     return s < 0? utoineg(c): utoipos(c);
    3842             :   }
    3843         407 :   v = vali(n); m0 = 1;
    3844         408 :   if (v)
    3845             :   {
    3846         123 :     n = shifti(n, -v); if (odd(v)) m0 *= 2;
    3847         123 :     copy = 0;
    3848             :   }
    3849         408 :   if ((mask = u_357_divides(umodiu(n, 3 * 5 * 7))))
    3850             :   {
    3851         275 :     if (mask & 1) { v = Z_lvalrem(n, 3, &n); if (odd(v)) m0 *= 3; }
    3852         275 :     if (mask & 2) { v = Z_lvalrem(n, 5, &n); if (odd(v)) m0 *= 5; }
    3853         275 :     if (mask & 4) { v = Z_lvalrem(n, 7, &n); if (odd(v)) m0 *= 7; ; }
    3854         275 :     copy = 0;
    3855             :   }
    3856         408 :   if (copy) n = absi(n); /* ifac_core destroys n */
    3857         289 :   else if (lgefint(n) == 3)
    3858             :   {
    3859          81 :     ulong c = coreu(uel(n,2));
    3860          81 :     m = muluu(m0, c); if (s < 0) setsigne(m, -1);
    3861          81 :     return gerepileuptoint(av, m);
    3862             :   }
    3863         327 :   setabssign(n); lim = tridiv_bound(n);
    3864             :   /* n >= 691^2 */
    3865         327 :   gcd = Z_oddprimedivisors_gcd(n, lim, &gcdlim);
    3866         327 :   n = core_init_from_squarefree(n, gcd, &mpart);
    3867         327 :   m = mului(m0, mpart); if (s < 0) setsigne(m, -1);
    3868         327 :   if (equali1(n)) return gerepileuptoint(av, m);
    3869             :   /* n has no prime divisor <= gcdlim */
    3870         276 :   if (cmpii(sqru(gcdlim + 1), n) > 0)
    3871           2 :     return gerepileuptoint(av, mulii(m, n)); /* prime */
    3872         274 :   l = lg(primetab);
    3873         750 :   for (i = 1; i < l; i++)
    3874             :   {
    3875         478 :     GEN q = gel(primetab,i);
    3876         478 :     v = Z_pvalrem(n, q, &n);
    3877         478 :     if (v)
    3878             :     {
    3879           8 :       if (v & 1) m = mulii(m, q);
    3880           8 :       if (is_pm1(n)) return gerepileuptoint(av, m);
    3881             :     }
    3882             :   }
    3883         272 :   if (!ifac_isprime(n)) n = ifac_core(n); /* composite without small factors */
    3884         272 :   return gerepileuptoint(av, mulii(m, n));
    3885             : }
    3886             : 
    3887             : long
    3888           0 : Z_issmooth(GEN m, ulong lim)
    3889             : {
    3890           0 :   pari_sp av = avma;
    3891           0 :   ulong p = 2;
    3892             :   forprime_t S;
    3893           0 :   u_forprime_init(&S, 2, lim);
    3894           0 :   while ((p = u_forprime_next_fast(&S)))
    3895             :   {
    3896             :     int stop;
    3897           0 :     (void)Z_lvalrem_stop(&m, p, &stop);
    3898           0 :     if (stop) return gc_long(av, abscmpiu(m,lim) <= 0);
    3899             :   }
    3900           0 :   return gc_long(av,0);
    3901             : }
    3902             : 
    3903             : GEN
    3904      178649 : Z_issmooth_fact(GEN m, ulong lim)
    3905             : {
    3906      178649 :   pari_sp av = avma;
    3907             :   GEN F, P, E;
    3908             :   ulong p;
    3909      178649 :   long i = 1, l = expi(m)+1;
    3910             :   forprime_t S;
    3911      178641 :   P = cgetg(l, t_VECSMALL);
    3912      178576 :   E = cgetg(l, t_VECSMALL); F = mkmat2(P,E);
    3913      178566 :   if (l == 1) return F; /* m == 1 */
    3914      178524 :   u_forprime_init(&S, 2, lim);
    3915    44193391 :   while ((p = u_forprime_next_fast(&S)))
    3916             :   {
    3917             :     int stop;
    3918    44151856 :     long v = Z_lvalrem_stop(&m, p, &stop);
    3919    44152639 :     if (v) { P[i] = p; E[i] = v; i++; }
    3920    44152639 :     if (stop)
    3921             :     {
    3922      137802 :       if (abscmpiu(m,lim) > 0) break;
    3923      111610 :       if (m[2] > 1) { P[i] = m[2]; E[i] = 1; i++; }
    3924      111610 :       setlg(P, i);
    3925      111692 :       setlg(E, i); return gc_const((pari_sp)F, F);
    3926             :     }
    3927             :   }
    3928       66057 :   return gc_NULL(av);
    3929             : }
    3930             : 
    3931             : /* Is (a mod p^e) a K-th power ? p is prime and e > 0 */
    3932             : static int
    3933         798 : Zp_ispower(GEN a, GEN L, GEN K, GEN p, long e)
    3934             : {
    3935         798 :   GEN t = gen_0;
    3936         798 :   long v = Z_pvalrem(a, p, &a), d = e - v;
    3937         798 :   if (d > 0)
    3938             :   { /* is a mod p^d a K-th power ? a a p-unit */
    3939             :     ulong r;
    3940         763 :     v = uabsdivui_rem(v, K, &r); if (r) return 0;
    3941         651 :     if (d == 1)
    3942             :     { /* mod p: faster */
    3943         567 :       if (!Fp_ispower(a, K, p)) return 0;
    3944         560 :       if (L) t = Fp_sqrtn(a, K, p, NULL);
    3945             :     }
    3946             :     else
    3947             :     { /* mod p^{2 +} */
    3948          84 :       if (!ispower(cvtop(a, p, d), K, L? &t: NULL)) return 0;
    3949          63 :       if (L) t = gtrunc(t);
    3950             :     }
    3951         623 :     if (L && v) t = mulii(t, powiu(p, v));
    3952             :   }
    3953         686 :   if (L) vectrunc_append(L, mkintmod(t, powiu(p, e)));
    3954         686 :   return 1;
    3955             : }
    3956             : long
    3957         756 : Zn_ispower(GEN a, GEN q, GEN K, GEN *pt)
    3958             : {
    3959             :   GEN L, N;
    3960             :   pari_sp av;
    3961             :   long e, i, l;
    3962             :   ulong pp, lim;
    3963             : 
    3964         756 :   if (!signe(a))
    3965             :   {
    3966          91 :     if (pt) {
    3967          91 :       GEN t = cgetg(3, t_INTMOD);
    3968          91 :       gel(t,1) = icopy(q); gel(t,2) = gen_0; *pt = t;
    3969             :     }
    3970          91 :     return 1;
    3971             :   }
    3972             :   /* a != 0 */
    3973         665 :   av = avma;
    3974             : 
    3975         665 :   if (typ(q) != t_INT) /* integer factorization */
    3976             :   {
    3977           0 :     GEN P = gel(q,1), E = gel(q,2);
    3978           0 :     l = lg(P);
    3979           0 :     L = pt? vectrunc_init(l): NULL;
    3980           0 :     for (i = 1; i < l; i++)
    3981             :     {
    3982           0 :       GEN p = gel(P,i);
    3983           0 :       long e = itos(gel(E,i));
    3984           0 :       if (!Zp_ispower(a, L, K, p, e)) return gc_long(av,0);
    3985             :     }
    3986           0 :     goto END;
    3987             :   }
    3988         665 :   e = vali(q); if (e) q = shifti(q, -e);
    3989         665 :   if (!mod2(K) && kronecker(a, q) == -1) return gc_long(av,0);
    3990         658 :   L = pt? vectrunc_init(expi(q)+2): NULL;
    3991         658 :   if (e)
    3992             :   {
    3993         469 :     if (!Zp_ispower(a, L, K, gen_2, e)) return gc_long(av,0);
    3994         455 :     a = modii(a, q);
    3995             :   }
    3996         644 :   lim = tridiv_bound(q);
    3997         644 :   if (cmpiu(q, 691 * 691) >= 0)
    3998             :   {
    3999         161 :     ulong sqrtq = lgefint(q) == 3? usqrt(q[2]): 0;
    4000             :     GEN P;
    4001         161 :     if (sqrtq) lim = minuu(sqrtq, lim);
    4002         161 :     P = Z_oddprimedivisors_fast(q, lim);
    4003         161 :     if (P)
    4004             :     {
    4005         103 :       long nP = lg(P) - 1;
    4006         103 :       int stop = 0;
    4007         206 :       for (i = 1; i <= nP; i++)
    4008             :       {
    4009         151 :         ulong pp = uel(P,i);
    4010         151 :         e = Z_lvalrem_stop(&q, pp, &stop);
    4011         151 :         if (!Zp_ispower(a, L, K, utoipos(pp), e)) return gc_long(av,0);
    4012         103 :         a = modii(a, q);
    4013             :       }
    4014          55 :       if (stop)
    4015             :       {
    4016          48 :         if (!is_pm1(q) && !Zp_ispower(a, L, K, q, 1)) return gc_long(av,0);
    4017          48 :         goto END;
    4018             :       }
    4019             :     }
    4020          58 :     else if (lim == sqrtq && lim <= GP_DATA->factorlimit)
    4021             :     {
    4022           0 :       if (!Zp_ispower(a, L, K, q, 1)) return gc_long(av,0);
    4023           0 :       goto END;
    4024             :     }
    4025             :   }
    4026             :   else
    4027             :   {
    4028             :     forprime_t S;
    4029         483 :     u_forprime_init(&S, 3, lim);
    4030      237174 :     while ((pp = u_forprime_next(&S)))
    4031             :     {
    4032             :       int stop;
    4033      236754 :       e = Z_lvalrem_stop(&q, pp, &stop);
    4034      236754 :       if (!e) continue;
    4035          63 :       if (!Zp_ispower(a, L, K, utoipos(pp), e)) return gc_long(av,0);
    4036          56 :       a = modii(a, q);
    4037          56 :       if (stop)
    4038             :       {
    4039          56 :         if (!is_pm1(q) && !Zp_ispower(a, L, K, q, 1)) return gc_long(av,0);
    4040          56 :         goto END;
    4041             :       }
    4042             :     }
    4043             :   }
    4044         485 :   l = lg(primetab);
    4045         485 :   for (i = 1; i < l; i++)
    4046             :   {
    4047           0 :     GEN p = gel(primetab,i);
    4048           0 :     e = Z_pvalrem(q, p, &q);
    4049           0 :     if (!e) continue;
    4050           0 :     if (!Zp_ispower(a, L, K, p, e)) return gc_long(av,0);
    4051           0 :     if (is_pm1(q)) goto END;
    4052           0 :     a = modii(a, q);
    4053             :   }
    4054         485 :   N = gcdii(a,q);
    4055         485 :   if (!is_pm1(N))
    4056             :   {
    4057          52 :     if (ifac_isprime(N))
    4058             :     {
    4059          34 :       e = Z_pvalrem(q, N, &q);
    4060          34 :       if (!Zp_ispower(a, L, K, N, e)) return gc_long(av,0);
    4061           0 :       a = modii(a, q);
    4062             :     }
    4063             :     else
    4064             :     {
    4065          18 :       GEN part = ifac_start(N, 0);
    4066             :       for(;;)
    4067          18 :       {
    4068             :         long e;
    4069             :         GEN p;
    4070          36 :         if (!ifac_next(&part, &p, &e)) break;
    4071          18 :         e = Z_pvalrem(q, p, &q);
    4072          18 :         if (!Zp_ispower(a, L, K, p, e)) return gc_long(av,0);
    4073          18 :         a = modii(a, q);
    4074             :       }
    4075             :     }
    4076             :   }
    4077         451 :   if (!is_pm1(q))
    4078             :   {
    4079          31 :     if (ifac_isprime(q))
    4080             :     {
    4081           4 :       if (!Zp_ispower(a, L, K, q, 1)) return gc_long(av,0);
    4082             :     }
    4083             :     else
    4084             :     { /* icopy needed: ifac_next would destroy q */
    4085          27 :       GEN part = ifac_start(icopy(q), 0);
    4086             :       for(;;)
    4087          43 :       {
    4088             :         long e;
    4089             :         GEN p;
    4090          70 :         if (!ifac_next(&part, &p, &e)) break;
    4091          52 :         if (!Zp_ispower(a, L, K, p, e)) return gc_long(av,0);
    4092          43 :         a = modii(a, q);
    4093             :       }
    4094             :     }
    4095             :   }
    4096         420 : END:
    4097         546 :   if (!pt) return gc_long(av, 1);
    4098         476 :   *pt = gerepileupto(av, chinese1_coprime_Z(L)); return 1;
    4099             : }
    4100             : 
    4101             : 
    4102             : /***********************************************************************/
    4103             : /**                                                                   **/
    4104             : /**       COMPUTING THE MATRIX OF PRIME DIVISORS AND EXPONENTS        **/
    4105             : /**                                                                   **/
    4106             : /***********************************************************************/
    4107             : static GEN
    4108      128365 : aux_end(GEN M, GEN n, long nb)
    4109             : {
    4110      128365 :   GEN P,E, z = (GEN)avma;
    4111             :   long i;
    4112             : 
    4113      128365 :   guncloneNULL(n);
    4114      128365 :   P = cgetg(nb+1,t_COL);
    4115      128365 :   E = cgetg(nb+1,t_COL);
    4116      846197 :   for (i=nb; i; i--)
    4117             :   { /* allow a stackdummy in the middle */
    4118      792964 :     while (typ(z) != t_INT) z += lg(z);
    4119      717832 :     gel(E,i) = z; z += lg(z);
    4120      717832 :     gel(P,i) = z; z += lg(z);
    4121             :   }
    4122      128365 :   gel(M,1) = P;
    4123      128365 :   gel(M,2) = E;
    4124      128365 :   return sort_factor(M, (void*)&abscmpii, cmp_nodata);
    4125             : }
    4126             : 
    4127             : static void
    4128      712838 : STORE(long *nb, GEN x, long e) { (*nb)++; (void)x; (void)utoipos(e); }
    4129             : static void
    4130      681523 : STOREu(long *nb, ulong x, long e) { STORE(nb, utoipos(x), e); }
    4131             : static void
    4132       31163 : STOREi(long *nb, GEN x, long e) { STORE(nb, icopy(x), e); }
    4133             : /* no prime less than p divides n; return 1 if factored completely */
    4134             : static int
    4135       39425 : special_primes(GEN n, ulong p, long *nb, GEN T)
    4136             : {
    4137       39425 :   long i, l = lg(T);
    4138       39425 :   if (l > 1)
    4139             :   { /* pp = square of biggest p tried so far */
    4140         540 :     long pp[] = { evaltyp(t_INT)|_evallg(4), 0,0,0 };
    4141         540 :     pari_sp av = avma; affii(sqru(p), pp); set_avma(av);
    4142             : 
    4143        1184 :     for (i = 1; i < l; i++)
    4144         777 :       if (dvdiiz(n, gel(T,i), n))
    4145             :       {
    4146         329 :         long k = 1; while (dvdiiz(n, gel(T,i), n)) k++;
    4147         231 :         STOREi(nb, gel(T,i), k);
    4148         231 :         if (abscmpii(pp, n) > 0)
    4149             :         {
    4150         133 :           if (!is_pm1(n)) STOREi(nb, n, 1);
    4151         133 :           return 1;
    4152             :         }
    4153             :       }
    4154             :   }
    4155       39292 :   return 0;
    4156             : }
    4157             : 
    4158             : /* factor(sn*|n|), where sn = -1 or 1.
    4159             :  * all != 0 : only look for prime divisors < all. If pU is not NULL,
    4160             :  * set it to unfactored composite */
    4161             : static GEN
    4162    13999600 : ifactor_sign(GEN n, ulong all, long hint, long sn, GEN *pU)
    4163             : {
    4164             :   GEN M, N;
    4165             :   pari_sp av;
    4166    13999600 :   long nb = 0, nb0 = -1, i;
    4167             :   ulong lim;
    4168             :   forprime_t T;
    4169             : 
    4170    13999600 :   if (lgefint(n) == 3)
    4171             :   { /* small integer */
    4172    13871250 :     GEN f, Pf, Ef, P, E, F = cgetg(3, t_MAT);
    4173             :     ulong U1, U2;
    4174             :     long l;
    4175    13871752 :     av = avma;
    4176             :     /* enough room to store <= 15 primes and exponents (OK if n < 2^64) */
    4177    13871752 :     (void)new_chunk((15*3 + 15 + 1) * 2);
    4178    13871744 :     f = factoru_sign(uel(n,2), all, hint, pU? &U1: NULL, pU? &U2: NULL);
    4179    13871956 :     set_avma(av);
    4180    13871954 :     Pf = gel(f,1);
    4181    13871954 :     Ef = gel(f,2);
    4182    13871954 :     l = lg(Pf);
    4183    13871954 :     if (sn < 0)
    4184             :     { /* add sign */
    4185        6519 :       long L = l+1;
    4186        6519 :       gel(F,1) = P = cgetg(L, t_COL);
    4187        6519 :       gel(F,2) = E = cgetg(L, t_COL);
    4188        6519 :       gel(P,1) = gen_m1; P++;
    4189        6519 :       gel(E,1) = gen_1;  E++;
    4190             :     }
    4191             :     else
    4192             :     {
    4193    13865435 :       gel(F,1) = P = cgetg(l, t_COL);
    4194    13865455 :       gel(F,2) = E = cgetg(l, t_COL);
    4195             :     }
    4196    42749118 :     for (i = 1; i < l; i++)
    4197             :     {
    4198    28877210 :       gel(P,i) = utoipos(Pf[i]);
    4199    28877152 :       gel(E,i) = utoipos(Ef[i]);
    4200             :     }
    4201    13871908 :     if (pU) *pU = U1 == 1? NULL: mkvec2(utoipos(U1), utoipos(U2));
    4202    13871908 :     return F;
    4203             :   }
    4204      128350 :   if (pU) *pU = NULL;
    4205      128350 :   M = cgetg(3,t_MAT);
    4206      128365 :   if (sn < 0) STORE(&nb, utoineg(1), 1);
    4207      128365 :   if (is_pm1(n)) return aux_end(M,NULL,nb);
    4208             : 
    4209      128365 :   n = N = gclone(n); setabssign(n);
    4210             :   /* trial division bound; look for primes <= lim; nb is the number of
    4211             :    * distinct prime factors so far; if nb0 >= 0, it records the value of nb
    4212             :    * for which we made a successful compositeness test: if later nb = nb0,
    4213             :    * we know that n is composite */
    4214      128365 :   lim = 1;
    4215      128365 :   if (!all || all > 2)
    4216             :   { /* trial divide p < all if all != 0, else p <= tridiv_bound() */
    4217             :     ulong maxp, p;
    4218             :     pari_sp av2;
    4219      128351 :     i = vali(n);
    4220      128351 :     if (i)
    4221             :     {
    4222       76854 :       STOREu(&nb, 2, i);
    4223       76854 :       av = avma; affii(shifti(n,-i), n); set_avma(av);
    4224             :     }
    4225      128351 :     if (is_pm1(n)) return aux_end(M,n,nb);
    4226      128235 :     lim = all? all-1: tridiv_bound(n);
    4227      128235 :     av = avma;
    4228      128235 :     if (lim >= 128)
    4229             :     { /* fast trial division */
    4230      126772 :       GEN Q = Z_oddprimedivisors_fast(n, lim);
    4231      126772 :       av2 = avma;
    4232      126772 :       if (Q)
    4233             :       {
    4234      123977 :         long l = lg(Q);
    4235      724350 :         for (i = 1; i < l; i++)
    4236             :         {
    4237      601361 :           pari_sp av3 = avma;
    4238      601361 :           ulong p = uel(Q, i);
    4239             :           long k;
    4240      601361 :           if (all && p >= all) break;
    4241      600373 :           k = Z_lvalrem(n, p, &n); /* > 0 */
    4242      600373 :           affii(n, N); n = N; set_avma(av3);
    4243      600373 :           STOREu(&nb, p, k);
    4244             :         }
    4245      123977 :         if (is_pm1(n))
    4246             :         {
    4247       88502 :           stackdummy(av, av2);
    4248       88502 :           return aux_end(M,n,nb);
    4249             :         }
    4250             :       }
    4251       38270 :       maxp = GP_DATA->factorlimit;
    4252             :     }
    4253             :     else
    4254             :     { /* naive trial division */
    4255        1463 :       maxp = maxprime();
    4256        1463 :       u_forprime_init(&T, 3, minuu(lim, maxp)); av2 = avma;
    4257             :       /* first pass: known to fit in private prime table */
    4258       31133 :       while ((p = u_forprime_next_fast(&T)))
    4259             :       {
    4260       29991 :         pari_sp av3 = avma;
    4261             :         int stop;
    4262       29991 :         long k = Z_lvalrem_stop(&n, p, &stop);
    4263       29991 :         if (k)
    4264             :         {
    4265        4295 :           affii(n, N); n = N; set_avma(av3);
    4266        4295 :           STOREu(&nb, p, k);
    4267             :         }
    4268             :         /* prodeuler(p=2,16381,1-1/p) ~ 0.0578; if probability of being prime
    4269             :          * knowing P^-(n) > 16381 is at least 10%, try BPSW */
    4270       29991 :         if (!stop && p == 16381)
    4271             :         {
    4272           0 :           if (bit_accuracy_mul(lgefint(n), 0.0578 * M_LN2) < 10)
    4273           0 :           { nb0 = nb; stop = ifac_isprime(n); }
    4274             :         }
    4275       29991 :         if (stop)
    4276             :         {
    4277         321 :           if (!is_pm1(n)) STOREi(&nb, n, 1);
    4278         321 :           stackdummy(av, av2);
    4279         321 :           return aux_end(M,n,nb);
    4280             :         }
    4281             :       }
    4282             :     }
    4283       39412 :     stackdummy(av, av2);
    4284       39412 :     if (lim > maxp)
    4285             :     { /* second pass usually empty, outside fast trial division range */
    4286           1 :       av = avma; u_forprime_init(&T, maxp+1, lim); av2 = avma;
    4287      882811 :       while ((p = u_forprime_next(&T)))
    4288             :       {
    4289      882811 :         pari_sp av3 = avma;
    4290             :         int stop;
    4291      882811 :         long k = Z_lvalrem_stop(&n, p, &stop);
    4292      882811 :         if (k)
    4293             :         {
    4294           1 :           affii(n, N); n = N; set_avma(av3);
    4295           1 :           STOREu(&nb, p, k);
    4296             :         }
    4297      882811 :         if (stop)
    4298             :         {
    4299           1 :           if (!is_pm1(n)) STOREi(&nb, n, 1);
    4300           1 :           stackdummy(av, av2);
    4301           1 :           return aux_end(M,n,nb);
    4302             :         }
    4303             :       }
    4304           0 :       stackdummy(av, av2);
    4305             :     }
    4306             :   }
    4307       39425 :   if (special_primes(n, lim, &nb, primetab)) return aux_end(M,n, nb);
    4308             :   /* if nb > nb0 (includes nb0 = -1) we don't know that n is composite */
    4309       39292 :   if (all)
    4310             :   { /* smallfact: look for easy pure powers then stop. Cf Z_isanypower */
    4311             :     GEN x;
    4312       31290 :     long k, e = expu(lim);
    4313       31290 :     av = avma;
    4314       29976 :     k = e >= 10? Z_isanypower_nosmalldiv(n, e, &x)
    4315       31290 :                : Z_isanypower(n, &x);
    4316       31290 :     if (k > 1) { affii(x, n); nb0 = -1; } else if (k < 1) k = 1;
    4317       31290 :     if (pU)
    4318             :     {
    4319             :       GEN F;
    4320       13156 :       if (abscmpiu(n, lim) <= 0
    4321       13156 :           || cmpii(n, sqru(lim)) <= 0
    4322        8833 :           || ((e >= 14) &&
    4323        7880 :               (nb>nb0 && bit_accuracy(lgefint(n))<2048 && ifac_isprime(n))))
    4324       13156 :       { set_avma(av); STOREi(&nb, n, k); return aux_end(M,n, nb); }
    4325        5946 :       set_avma(av); F = aux_end(M, NULL, nb); /* don't destroy n */
    4326        5946 :       *pU = mkvec2(icopy(n), utoipos(k)); /* composite cofactor */
    4327        5946 :       gunclone(n); return F;
    4328             :     }
    4329       18134 :     set_avma(av); STOREi(&nb, n, k);
    4330       18134 :     if (DEBUGLEVEL >= 2) {
    4331           0 :       pari_warn(warner,
    4332           0 :         "IFAC: untested %ld-bit integer declared prime", expi(n)+1);
    4333           0 :       if (expi(n) <= 256) err_printf("\t%Ps\n", n);
    4334             :     }
    4335             :   }
    4336        8002 :   else if (nb > nb0 && ifac_isprime(n)) STOREi(&nb, n, 1);
    4337        2568 :   else nb += ifac_decomp(n, hint);
    4338       26136 :   return aux_end(M,n, nb);
    4339             : }
    4340             : 
    4341             : static GEN
    4342     9541902 : ifactor(GEN n, ulong all, long hint)
    4343             : {
    4344     9541902 :   long s = signe(n);
    4345     9541902 :   if (!s) retmkmat2(mkcol(gen_0), mkcol(gen_1));
    4346     9541853 :   return ifactor_sign(n, all, hint, s, NULL);
    4347             : }
    4348             : 
    4349             : int
    4350        6651 : ifac_next(GEN *part, GEN *p, long *e)
    4351             : {
    4352        6651 :   GEN here = ifac_main(part);
    4353        6651 :   if (here == gen_0) { *p = NULL; *e = 1; return 0; }
    4354        6649 :   if (!here) { *p = NULL; *e = 0; return 0; }
    4355        4853 :   *p = VALUE(here);
    4356        4853 :   *e = EXPON(here)[2];
    4357        4853 :   ifac_delete(here); return 1;
    4358             : }
    4359             : 
    4360             : /* see before ifac_crack for current semantics of 'hint' (factorint's 'flag') */
    4361             : GEN
    4362       10290 : factorint(GEN n, long flag)
    4363             : {
    4364             :   GEN F;
    4365       10290 :   if ((F = check_arith_all(n,"factorint"))) return gcopy(F);
    4366       10276 :   return ifactor(n,0,flag);
    4367             : }
    4368             : 
    4369             : GEN
    4370       49797 : Z_factor_limit(GEN n, ulong all)
    4371             : {
    4372       49797 :   if (!all) all = GP_DATA->factorlimit + 1;
    4373       49797 :   return ifactor(n, all, decomp_default_hint);
    4374             : }
    4375             : GEN
    4376      886941 : absZ_factor_limit_strict(GEN n, ulong all, GEN *pU)
    4377             : {
    4378             :   GEN F, U;
    4379      886941 :   if (!signe(n))
    4380             :   {
    4381           0 :     if (pU) *pU = NULL;
    4382           0 :     retmkmat2(mkcol(gen_0), mkcol(gen_1));
    4383             :   }
    4384      886941 :   if (!all) all = GP_DATA->factorlimit + 1;
    4385      886941 :   F = ifactor_sign(n, all, decomp_default_hint, 1, &U);
    4386      886972 :   if (pU) *pU = U;
    4387      886972 :   return F;
    4388             : }
    4389             : GEN
    4390      290811 : absZ_factor_limit(GEN n, ulong all)
    4391             : {
    4392      290811 :   if (!signe(n)) retmkmat2(mkcol(gen_0), mkcol(gen_1));
    4393      290811 :   if (!all) all = GP_DATA->factorlimit + 1;
    4394      290811 :   return ifactor_sign(n, all, decomp_default_hint, 1, NULL);
    4395             : }
    4396             : GEN
    4397     9481780 : Z_factor(GEN n)
    4398     9481780 : { return ifactor(n,0,decomp_default_hint); }
    4399             : GEN
    4400     3277039 : absZ_factor(GEN n)
    4401             : {
    4402     3277039 :   if (!signe(n)) retmkmat2(mkcol(gen_0), mkcol(gen_1));
    4403     3277025 :   return ifactor_sign(n, 0, decomp_default_hint, 1, NULL);
    4404             : }
    4405             : /* Factor until the unfactored part is smaller than limit. Return the
    4406             :  * factored part. Hence factorback(output) may be smaller than n */
    4407             : GEN
    4408        3045 : Z_factor_until(GEN n, GEN limit)
    4409             : {
    4410        3045 :   pari_sp av = avma;
    4411        3045 :   long s = signe(n), eq;
    4412             :   GEN q, F, U;
    4413             : 
    4414        3045 :   if (!s) retmkmat2(mkcol(gen_0), mkcol(gen_1));
    4415        3045 :   F = ifactor_sign(n, tridiv_bound(n), decomp_default_hint, s, &U);
    4416        3045 :   if (!U) return F;
    4417        1155 :   q = gel(U,1); /* composite, q^eq = unfactored part */
    4418        1155 :   eq = itou(gel(U,2));
    4419        1155 :   if (cmpii(eq == 1? q: powiu(q, eq), limit) > 0)
    4420             :   { /* factor further */
    4421        1022 :     long l2 = expi(q)+1;
    4422             :     GEN P2, E2, F2, part;
    4423        1022 :     if (eq > 1) limit = sqrtnint(limit, eq);
    4424        1022 :     P2 = coltrunc_init(l2);
    4425        1022 :     E2 = coltrunc_init(l2); F2 = mkmat2(P2,E2);
    4426        1022 :     part = ifac_start(icopy(q), 0); /* ifac_next would destroy q */
    4427             :     for(;;)
    4428          70 :     {
    4429             :       long e;
    4430             :       GEN p;
    4431        1092 :       if (!ifac_next(&part,&p,&e)) break;
    4432        1092 :       vectrunc_append(P2, p);
    4433        1092 :       vectrunc_append(E2, utoipos(e * eq));
    4434        1092 :       q = diviiexact(q, powiu(p, e));
    4435        1092 :       if (cmpii(q, limit) <= 0) break;
    4436             :     }
    4437        1022 :     F2 = sort_factor(F2, (void*)&abscmpii, cmp_nodata);
    4438        1022 :     F = merge_factor(F, F2, (void*)&abscmpii, cmp_nodata);
    4439             :   }
    4440        1155 :   return gerepilecopy(av, F);
    4441             : }
    4442             : 
    4443             : static void
    4444    96569823 : matsmalltrunc_append(GEN m, ulong p, ulong e)
    4445             : {
    4446    96569823 :   GEN P = gel(m,1), E = gel(m,2);
    4447    96569823 :   long l = lg(P);
    4448    96569823 :   P[l] = p; lg_increase(P);
    4449    96573968 :   E[l] = e; lg_increase(E);
    4450    96554374 : }
    4451             : static GEN
    4452    37767785 : matsmalltrunc_init(long l)
    4453             : {
    4454    37767785 :   GEN P = vecsmalltrunc_init(l);
    4455    37730662 :   GEN E = vecsmalltrunc_init(l); return mkvec2(P,E);
    4456             : }
    4457             : 
    4458             : /* return optimal N s.t. omega(b) <= N for all b <= x */
    4459             : long
    4460       71606 : maxomegau(ulong x)
    4461             : { /* P=primes(15); for(i=1,15, print([i, vecprod(P[1..i])])) */
    4462       71606 :   if (x < 30030UL)/* rare trivial cases */
    4463             :   {
    4464       37455 :     if (x < 2UL) return 0;
    4465       19199 :     if (x < 6UL) return 1;
    4466       13599 :     if (x < 30UL) return 2;
    4467       12892 :     if (x < 210UL) return 3;
    4468       12633 :     if (x < 2310UL) return 4;
    4469       11604 :     return 5;
    4470             :   }
    4471       34151 :   if (x < 510510UL) return 6; /* most frequent case */
    4472       18753 :   if (x < 9699690UL) return 7;
    4473           7 :   if (x < 223092870UL) return 8;
    4474             : #ifdef LONG_IS_64BIT
    4475           6 :   if (x < 6469693230UL) return 9;
    4476           0 :   if (x < 200560490130UL) return 10;
    4477           0 :   if (x < 7420738134810UL) return 11;
    4478           0 :   if (x < 304250263527210UL) return 12;
    4479           0 :   if (x < 13082761331670030UL) return 13;
    4480           0 :   if (x < 614889782588491410UL) return 14;
    4481           0 :   return 15;
    4482             : #else
    4483           1 :   return 9;
    4484             : #endif
    4485             : }
    4486             : /* return optimal N s.t. omega(b) <= N for all odd b <= x */
    4487             : long
    4488        2229 : maxomegaoddu(ulong x)
    4489             : { /* P=primes(15+1); for(i=1,15, print([i, vecprod(P[2..i+1])])) */
    4490        2229 :   if (x < 255255UL)/* rare trivial cases */
    4491             :   {
    4492        1355 :     if (x < 3UL) return 0;
    4493        1355 :     if (x < 15UL) return 1;
    4494        1355 :     if (x < 105UL) return 2;
    4495        1355 :     if (x < 1155UL) return 3;
    4496        1327 :     if (x < 15015UL) return 4;
    4497        1327 :     return 5;
    4498             :   }
    4499         874 :   if (x < 4849845UL) return 6; /* most frequent case */
    4500           0 :   if (x < 111546435UL) return 7;
    4501           0 :   if (x < 3234846615UL) return 8;
    4502             : #ifdef LONG_IS_64BIT
    4503           0 :   if (x < 100280245065UL) return 9;
    4504           0 :   if (x < 3710369067405UL) return 10;
    4505           0 :   if (x < 152125131763605UL) return 11;
    4506           0 :   if (x < 6541380665835015UL) return 12;
    4507           0 :   if (x < 307444891294245705UL) return 13;
    4508           0 :   if (x < 16294579238595022365UL) return 14;
    4509           0 :   return 15;
    4510             : #else
    4511           0 :   return 9;
    4512             : #endif
    4513             : }
    4514             : 
    4515             : /* If a <= c <= b , factoru(c) = L[c-a+1] */
    4516             : GEN
    4517       31925 : vecfactoru_i(ulong a, ulong b)
    4518             : {
    4519       31925 :   ulong k, p, n = b-a+1, N = maxomegau(b) + 1;
    4520       31925 :   GEN v = const_vecsmall(n, 1);
    4521       31925 :   GEN L = cgetg(n+1, t_VEC);
    4522             :   forprime_t T;
    4523    20783074 :   for (k = 1; k <= n; k++) gel(L,k) = matsmalltrunc_init(N);
    4524       31925 :   u_forprime_init(&T, 2, usqrt(b));
    4525      883441 :   while ((p = u_forprime_next(&T)))
    4526             :   { /* p <= sqrt(b) */
    4527      853170 :     ulong pk = p, K = ulogint(b, p);
    4528     2967666 :     for (k = 1; k <= K; k++)
    4529             :     {
    4530     2116150 :       ulong j, t = a / pk, ap = t * pk;
    4531     2116150 :       if (ap < a) { ap += pk; t++; }
    4532             :       /* t = (j+a-1) \ pk */
    4533     2116150 :       t %= p;
    4534    60953760 :       for (j = ap-a+1; j <= n; j += pk)
    4535             :       {
    4536    58839268 :         if (t) { v[j] *= pk; matsmalltrunc_append(gel(L,j), p,k); }
    4537    58837610 :         if (++t == p) t = 0;
    4538             :       }
    4539     2114492 :       pk *= p;
    4540             :     }
    4541             :   }
    4542             :   /* complete factorisation of non-sqrt(b)-smooth numbers */
    4543    21150542 :   for (k = 1, N = a; k <= n; k++, N++)
    4544    21118326 :     if (uel(v,k) != N) matsmalltrunc_append(gel(L,k), N/uel(v,k),1UL);
    4545       32216 :   return L;
    4546             : }
    4547             : GEN
    4548           0 : vecfactoru(ulong a, ulong b)
    4549             : {
    4550           0 :   pari_sp av = avma;
    4551           0 :   return gerepilecopy(av, vecfactoru_i(a,b));
    4552             : }
    4553             : 
    4554             : /* Assume a and b odd, return L s.t. L[k] = factoru(a + 2*(k-1))
    4555             :  * If a <= c <= b odd, factoru(c) = L[(c-a)>>1 + 1] */
    4556             : GEN
    4557        2229 : vecfactoroddu_i(ulong a, ulong b)
    4558             : {
    4559        2229 :   ulong k, p, n = ((b-a)>>1) + 1, N = maxomegaoddu(b) + 1;
    4560        2229 :   GEN v = const_vecsmall(n, 1);
    4561        2229 :   GEN L = cgetg(n+1, t_VEC);
    4562             :   forprime_t T;
    4563             : 
    4564    17321271 :   for (k = 1; k <= n; k++) gel(L,k) = matsmalltrunc_init(N);
    4565        2229 :   u_forprime_init(&T, 3, usqrt(b));
    4566      185277 :   while ((p = u_forprime_next(&T)))
    4567             :   { /* p <= sqrt(b) */
    4568      183064 :     ulong pk = p, K = ulogint(b, p);
    4569      623601 :     for (k = 1; k <= K; k++)
    4570             :     {
    4571      440553 :       ulong j, t = (a / pk) | 1UL, ap = t * pk;
    4572             :       /* t and ap are odd, ap multiple of pk = p^k */
    4573      440553 :       if (ap < a) { ap += pk<<1; t+=2; }
    4574             :       /* c=t*p^k by steps of 2*p^k; factorization of c*=p^k if (t,p)=1 */
    4575      440553 :       t %= p;
    4576    32504310 :       for (j = ((ap-a)>>1)+1; j <= n; j += pk)
    4577             :       {
    4578    32063881 :         if (t) { v[j] *= pk; matsmalltrunc_append(gel(L,j), p,k); }
    4579    32063757 :         t += 2; if (t >= p) t -= p;
    4580             :       }
    4581      440429 :       pk *= p;
    4582             :     }
    4583             :   }
    4584             :   /* complete factorisation of non-sqrt(b)-smooth numbers */
    4585    17688492 :   for (k = 1, N = a; k <= n; k++, N+=2)
    4586    17686437 :     if (uel(v,k) != N) matsmalltrunc_append(gel(L,k), N/uel(v,k),1UL);
    4587        2055 :   return L;
    4588             : }
    4589             : GEN
    4590           0 : vecfactoroddu(ulong a, ulong b)
    4591             : {
    4592           0 :   pari_sp av = avma;
    4593           0 :   return gerepilecopy(av, vecfactoroddu_i(a,b));
    4594             : }
    4595             : 
    4596             : /* If 0 <= a <= c <= b; L[c-a+1] = factoru(c)[,1] if c squarefree, else NULL */
    4597             : GEN
    4598        7014 : vecfactorsquarefreeu(ulong a, ulong b)
    4599             : {
    4600        7014 :   ulong k, p, n = b-a+1, N = maxomegau(b) + 1;
    4601        7014 :   GEN v = const_vecsmall(n, 1);
    4602        7014 :   GEN L = cgetg(n+1, t_VEC);
    4603             :   forprime_t T;
    4604    14007238 :   for (k = 1; k <= n; k++) gel(L,k) = vecsmalltrunc_init(N);
    4605        7014 :   u_forprime_init(&T, 2, usqrt(b));
    4606      838334 :   while ((p = u_forprime_next(&T)))
    4607             :   { /* p <= sqrt(b), kill nonsquarefree */
    4608      831320 :     ulong j, pk = p*p, t = a / pk, ap = t * pk;
    4609      831320 :     if (ap < a) ap += pk;
    4610     7160090 :     for (j = ap-a+1; j <= n; j += pk) gel(L,j) = NULL;
    4611             : 
    4612      831320 :     t = a / p; ap = t * p;
    4613      831320 :     if (ap < a) { ap += p; t++; }
    4614    30551556 :     for (j = ap-a+1; j <= n; j += p, t++)
    4615    29720236 :       if (gel(L,j)) { v[j] *= p; vecsmalltrunc_append(gel(L,j), p); }
    4616             :   }
    4617             :   /* complete factorisation of non-sqrt(b)-smooth numbers */
    4618    14007238 :   for (k = 1, N = a; k <= n; k++, N++)
    4619    14000224 :     if (gel(L,k) && uel(v,k) != N) vecsmalltrunc_append(gel(L,k), N/uel(v,k));
    4620        7014 :   return L;
    4621             : }
    4622             : /* If 0 <= a <= c <= b; L[c-a+1] = factoru(c)[,1] if c squarefree and coprime
    4623             :  * to all the primes in sorted zv P, else NULL */
    4624             : GEN
    4625       32590 : vecfactorsquarefreeu_coprime(ulong a, ulong b, GEN P)
    4626             : {
    4627       32590 :   ulong k, p, n = b-a+1, sqb = usqrt(b), N = maxomegau(b) + 1;
    4628       32590 :   GEN v = const_vecsmall(n, 1);
    4629       32591 :   GEN L = cgetg(n+1, t_VEC);
    4630             :   forprime_t T;
    4631    90637851 :   for (k = 1; k <= n; k++) gel(L,k) = vecsmalltrunc_init(N);
    4632       32592 :   u_forprime_init(&T, 2, sqb);
    4633     3680347 :   while ((p = u_forprime_next(&T)))
    4634             :   { /* p <= sqrt(b), kill nonsquarefree */
    4635     3648305 :     ulong j, t, ap, bad = zv_search(P, p), pk = bad ? p: p * p;
    4636     3648553 :     t = a / pk; ap = t * pk; if (ap < a) ap += pk;
    4637    80847095 :     for (j = ap-a+1; j <= n; j += pk) gel(L,j) = NULL;
    4638     3648553 :     if (bad) continue;
    4639             : 
    4640     3586400 :     t = a / p; ap = t * p;
    4641     3586400 :     if (ap < a) { ap += p; t++; }
    4642   116461758 :     for (j = ap-a+1; j <= n; j += p, t++)
    4643   112876156 :       if (gel(L,j)) { v[j] *= p; vecsmalltrunc_append(gel(L,j), p); }
    4644             :   }
    4645       32592 :   if (uel(P,lg(P)-1) <= sqb) P = NULL;
    4646             :   /* complete factorisation of non-sqrt(b)-smooth numbers */
    4647    90824890 :   for (k = 1, N = a; k <= n; k++, N++)
    4648    90792409 :     if (gel(L,k) && uel(v,k) != N)
    4649             :     {
    4650    25662448 :       ulong q = N / uel(v,k);
    4651    25662448 :       if (!P || !zv_search(P, q)) vecsmalltrunc_append(gel(L,k), q);
    4652             :     }
    4653       32481 :   return L;
    4654             : }
    4655             : 
    4656             : GEN
    4657          49 : vecsquarefreeu(ulong a, ulong b)
    4658             : {
    4659          49 :   ulong j, k, p, n = b-a+1;
    4660          49 :   GEN L = const_vecsmall(n, 1);
    4661             :   forprime_t T;
    4662          49 :   u_forprime_init(&T, 2, usqrt(b));
    4663         462 :   while ((p = u_forprime_next(&T)))
    4664             :   { /* p <= sqrt(b), kill nonsquarefree */
    4665         413 :     ulong pk = p*p, t = a / pk, ap = t * pk;
    4666         413 :     if (ap < a) { ap += pk; t++; }
    4667             :     /* t = (j+a-1) \ pk */
    4668       21777 :     for (j = ap-a+1; j <= n; j += pk, t++) L[j] = 0;
    4669             :   }
    4670       48258 :   for (k = j = 1; k <= n; k++)
    4671       48209 :     if (L[k]) L[j++] = a+k-1;
    4672          49 :   setlg(L,j); return L;
    4673             : }

Generated by: LCOV version 1.16