[BACK]Return to ifactor1.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / pari-2.2 / src / basemath

Diff for /OpenXM_contrib/pari-2.2/src/basemath/Attic/ifactor1.c between version 1.1 and 1.2

version 1.1, 2001/10/02 11:17:04 version 1.2, 2002/09/11 07:26:51
Line 19  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
Line 19  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 /**                                                                **/  /**                                                                **/
 /********************************************************************/  /********************************************************************/
 #include "pari.h"  #include "pari.h"
   extern GEN decomp_limit(GEN n, GEN limit);
   extern int BSW_isprime(GEN x);
   extern int BSW_isprime_small(GEN x);
   extern ulong ucarrecomplet(ulong A);
   extern GEN mpqs(GEN N);/* in src/modules/mpqs.c,
                           * returns a factor, a vector of factors, or NULL */
   
 /*********************************************************************/  /*********************************************************************/
 /**                                                                 **/  /**                                                                 **/
Line 47  init_miller(GEN n)
Line 53  init_miller(GEN n)
 static int  static int
 bad_for_base(GEN n, GEN a)  bad_for_base(GEN n, GEN a)
 {  {
   long r, av=avma, lim=stack_lim(av,1);    long r;
     gpmem_t av=avma, lim=stack_lim(av, 1);
   GEN c2, c = powmodulo(a,t1,n);    GEN c2, c = powmodulo(a,t1,n);
   
   if (!is_pm1(c) && !egalii(t,c)) /* go fishing for -1, not for 1 */    if (!is_pm1(c) && !egalii(t,c)) /* go fishing for -1, not for 1 */
Line 79  bad_for_base(GEN n, GEN a)
Line 86  bad_for_base(GEN n, GEN a)
 long  long
 millerrabin(GEN n, long k)  millerrabin(GEN n, long k)
 {  {
   long r,i,av2, av = avma;    long r, i;
     gpmem_t av2, av = avma;
   
   if (!signe(n)) return 0;    if (!signe(n)) return 0;
   /* If |n| <= 3, check if n = +- 1 */    /* If |n| <= 3, check if n = +- 1 */
Line 135  millerrabin(GEN n, long k)
Line 143  millerrabin(GEN n, long k)
 int                             /* no longer static -- needed in mpqs.c */  int                             /* no longer static -- needed in mpqs.c */
 miller(GEN n, long k)  miller(GEN n, long k)
 {  {
   long r,i,av2, av = avma;    long r, i;
     gpmem_t av2, av = avma;
   static long pr[] =    static long pr[] =
     { 0, 2,3,5,7,11,13,17,19,23,29, 31,73, 2,13,23,1662803UL, };      { 0, 2,3,5,7,11,13,17,19,23,29, 31,73, 2,13,23,1662803UL, };
   long *p;    long *p;
Line 163  miller(GEN n, long k)
Line 172  miller(GEN n, long k)
   }    }
   avma=av; return 1;    avma=av; return 1;
 }  }
   
   /* compute n-th term of Lucas sequence modulo N.
    * v_{k+2} = P v_{k+1} - v_k, v_0 = 2, v_1 = P.
    * Assume n > 0 */
   static GEN
   LucasMod(GEN n, long P, GEN N)
   {
     gpmem_t av = avma, lim = stack_lim(av, 1);
     GEN nd = n+2;
     long i, m = *nd, j = 1+bfffo((ulong)m);
     GEN v = stoi(P), v1 = stoi(P*P - 2);
   
     m <<= j; j = BITS_IN_LONG - j;
     for (i=lgefint(n)-2;;) /* cf. leftright_pow */
     {
       for (; j; m<<=1,j--)
       { /* v = v_k, v1 = v_{k+1} */
         if (m < 0)
         { /* set v = v_{2k+1}, v1 = v_{2k+2} */
           v = subis(mulii(v,v1), P);
           v1= subis(sqri(v1), 2);
         }
         else
         {/* set v = v_{2k}, v1 = v_{2k+1} */
           v1= subis(mulii(v,v1), P);
           v = subis(sqri(v), 2);
         }
         v = modii(v, N);
         v1= modii(v1,N);
         if (low_stack(lim,stack_lim(av,1)))
         {
           GEN *gptr[2]; gptr[0]=&v; gptr[1]=&v1;
           if(DEBUGMEM>1) err(warnmem,"LucasMod");
           gerepilemany(av,gptr,2);
         }
       }
       if (--i == 0) return v;
       j = BITS_IN_LONG;
       m = *++nd;
     }
   }
   
   /* check that N not a square first (taken care of here, but inefficient)
    * assume N > 3 */
   static int
   IsLucasPsP0(GEN N)
   {
     long b, i, v;
     GEN N_2, m, z;
   
     for (b=3, i=0;; b+=2, i++)
     {
       if (i == 64 && carreparfait(N)) return 0; /* avoid oo loop if N = m^2 */
       if (krosg(b*b - 4, N) < 0) break;
     }
   
     m = addis(N,1); v = vali(m); m = shifti(m,-v);
     z = LucasMod(m, b, N);
     if (egalii(z, gdeux)) return 1;
     N_2 = subis(N,2);
     if (egalii(z, N_2)) return 1;
     for (i=1; i<v; i++)
     {
       if (!signe(z)) return 1;
       z = modii(subis(sqri(z), 2), N);
       if (egalii(z, gdeux)) return 0;
     }
     return 0;
   }
   
   long
   BSW_psp(GEN N)
   {
     gpmem_t av = avma;
     int k;
     GEN T;
   
     if (typ(N) != t_INT) err(arither1);
     if (signe(N) <= 0) return 0;
     if (!is_bigint(N))
       switch(itos(N))
       {
         case 1: return 0;
         case 2:
         case 3: return 1;
       }
     if (!mod2(N)) return 0;
   
     T = init_miller(N);
     k = bad_for_base(T,gdeux);
     avma = av;
     k = (!k && IsLucasPsP0(N));
     avma = av; return k;
   }
   
 /***********************************************************************/  /***********************************************************************/
 /**                                                                   **/  /**                                                                   **/
 /**                       Pocklington-Lehmer                          **/  /**                       Pocklington-Lehmer                          **/
Line 173  miller(GEN n, long k)
Line 277  miller(GEN n, long k)
 /*assume n>=2*/  /*assume n>=2*/
 static long pl831(GEN N, GEN p)  static long pl831(GEN N, GEN p)
 {  {
   ulong ltop=avma,av;    gpmem_t ltop=avma, av;
   long a;    long a;
   GEN Nmun,Nmunp;    GEN Nmun,Nmunp;
   Nmun=addis(N,-1);    Nmun=addis(N,-1);
Line 187  static long pl831(GEN N, GEN p)
Line 291  static long pl831(GEN N, GEN p)
     {      {
       GEN g;        GEN g;
       g=mppgcd(addis(b,-1),N);        g=mppgcd(addis(b,-1),N);
       if (gcmp1(g))        if (gcmp1(g)) { avma=ltop; return a; }
       {        if (!gegal(g,N)) { avma=ltop; return 0; }
         avma=ltop;  
         return a;  
       }  
       if (!gegal(g,N))  
       {  
         avma=ltop;  
         return 0;  
       }  
     }      }
     else      else { avma=ltop; return 0; }
     {  
       avma=ltop;  
       return 0;  
     }  
     avma=av;      avma=av;
   }    }
 }  }
 /*  /* Assume N is a strong BSW pseudoprime
    *
  * flag 0: return gun (prime), gzero (composite)   * flag 0: return gun (prime), gzero (composite)
  * flag 1: return gzero (composite), gun (small prime), matrix (large prime)   * flag 1: return gzero (composite), gun (small prime), matrix (large prime)
  *   *
  * The matrix has 3 columns, [a,b,c] with   * The matrix has 3 columns, [a,b,c] with
  * a[i] prime factor of N-1,   * a[i] prime factor of N-1,
  * b[i] witness for a[i] as in pl831   * b[i] witness for a[i] as in pl831
  * c[i] plisprime(a[i])   * c[i] plisprime(a[i])
  */   */
 extern GEN decomp_limit(GEN n, GEN limit);  
 GEN  GEN
 plisprime(GEN N, long flag)  plisprime(GEN N, long flag)
 {  {
   ulong ltop=avma;    gpmem_t ltop = avma;
   long i;    long i, l, t = typ(N);
   int eps;    int eps;
   GEN C,F;    GEN C, F = NULL;
   if ( typ(N) != t_INT ) err(arither1);  
     if (t == t_VEC)
     { /* [ N, [p1,...,pk] ], pi list of pseudoprime divisors of N */
       F = (GEN)N[2];
       N = (GEN)N[1]; t = typ(N);
     }
     if (t != t_INT) err(arither1);
     if (DEBUGLEVEL>3) fprintferr("PL: proving primality of N = %Z\n", N);
   
   eps = absi_cmp(N,gdeux);    eps = absi_cmp(N,gdeux);
   if (eps<=0) return eps? gzero: gun;    if (eps<=0) return eps? gzero: gun;
   
   N = absi(N);    N = absi(N);
   /* Use Jaeschke results. See above */    if (!F)
   if (miller(N,7))    {
   { /* compare to 341550071728321 */      F = (GEN)decomp_limit(addis(N,-1), racine(N))[1];
     if (cmpii(N, u2toi(0x136a3, 0x52b2c8c1)) < 0) { avma=ltop; return gun; }      if (DEBUGLEVEL>3) fprintferr("PL: N-1 factored!\n");
   }    }
   else { avma=ltop; return gzero; }  
   F=(GEN)decomp_limit(addis(N,-1),racine(N))[1];    C = cgetg(4,t_MAT); l = lg(F);
   if (DEBUGLEVEL>=3) fprintferr("P.L.:factor O.K.\n");    C[1] = lgetg(l,t_COL);
   C=cgetg(4,t_MAT);    C[2] = lgetg(l,t_COL);
   C[1]=lgetg(lg(F),t_COL);    C[3] = lgetg(l,t_COL);
   C[2]=lgetg(lg(F),t_COL);    for(i=1; i<l; i++)
   C[3]=lgetg(lg(F),t_COL);  
   for(i=1;i<lg(F);i++)  
   {    {
     long witness;      GEN p = (GEN)F[i], r;
     GEN p;      long witness = pl831(N,p);
     p=(GEN)F[i];  
     witness=pl831(N,p);      if (!witness) { avma = ltop; return gzero; }
     if (!witness)      mael(C,1,i) = licopy(p);
       mael(C,2,i) = lstoi(witness);
       if (!flag) r = BSW_isprime(p)? gun: gzero;
       else
     {      {
       avma=ltop;        if (BSW_isprime_small(p)) r = gun;
       return gzero;        else if (expi(p) > 250) r = isprimeAPRCL(p)? gdeux: gzero;
         else r = plisprime(p,flag);
     }      }
     mael(C,1,i)=lcopy(p);      mael(C,3,i) = (long)r;
     mael(C,2,i)=lstoi(witness);      if (r == gzero) err(talker,"False prime number %Z in plisprime", p);
     mael(C,3,i)=(long)plisprime(p,flag);  
     if (gmael(C,3,i)==gzero)  
       err(talker,"Sorry false prime number %Z in plisprime",p);  
   }    }
   if (!flag)   { avma=ltop; return gun; }    if (!flag) { avma = ltop; return gun; }
   return gerepileupto(ltop,C);    return gerepileupto(ltop,C);
 }  }
   
Line 311  unsigned char prc210_d1[] =
Line 411  unsigned char prc210_d1[] =
 GEN  GEN
 nextprime(GEN n)  nextprime(GEN n)
 {  {
   long rc,rc0,rcd,rcn,av1,av2, av = avma;    long rc, rc0, rcd, rcn;
     gpmem_t av1, av2, av = avma;
   
   if (typ(n) != t_INT) n=gceil(n); /* accept arguments in R --GN */    if (typ(n) != t_INT) n=gceil(n); /* accept arguments in R --GN */
   if (typ(n) != t_INT) err(arither1);    if (typ(n) != t_INT) err(arither1);
Line 339  nextprime(GEN n)
Line 440  nextprime(GEN n)
   av2 = av1 = avma;    av2 = av1 = avma;
   for(;;)    for(;;)
   {    {
     if (miller(n,10)) break;      if (BSW_psp(n)) break;
     av1 = avma;      av1 = avma;
     rcd = prc210_d1[rcn];      rcd = prc210_d1[rcn];
     if (++rcn > 47) rcn = 0;      if (++rcn > 47) rcn = 0;
Line 352  nextprime(GEN n)
Line 453  nextprime(GEN n)
 GEN  GEN
 precprime(GEN n)  precprime(GEN n)
 {  {
   long rc,rc0,rcd,rcn,av1,av2, av = avma;    long rc, rc0, rcd, rcn;
     gpmem_t av1, av2, av = avma;
   
   if (typ(n) != t_INT) n=gfloor(n); /* accept arguments in R --GN */    if (typ(n) != t_INT) n=gfloor(n); /* accept arguments in R --GN */
   if (typ(n) != t_INT) err(arither1);    if (typ(n) != t_INT) err(arither1);
Line 381  precprime(GEN n)
Line 483  precprime(GEN n)
   av2 = av1 = avma;    av2 = av1 = avma;
   for(;;)    for(;;)
   {    {
     if (miller(n,10)) break;      if (BSW_psp(n)) break;
     av1 = avma;      av1 = avma;
     if (rcn == 0)      if (rcn == 0)
     { rcd = 2; rcn = 47; }      { rcd = 2; rcn = 47; }
Line 408  ulong
Line 510  ulong
 snextpr(ulong p, byteptr *d, long *rcn, long *q, long k)  snextpr(ulong p, byteptr *d, long *rcn, long *q, long k)
 {  {
   static ulong pp[] =    static ulong pp[] =
     { evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3), 0 };      { evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 0 };
   static ulong *pp2 = pp + 2;    static ulong *pp2 = pp + 2;
   static GEN gp = (GEN)pp;    static GEN gp = (GEN)pp;
   long d1 = **d, rcn0;    long rcn0;
   
   if (d1)    if (**d)
   {    {
       byteptr dd = *d;
       long d1 = 0;
   
       NEXT_PRIME_VIADIFF(d1,dd);
     if (*rcn != NPRC)      if (*rcn != NPRC)
     {      {
       rcn0 = *rcn;        rcn0 = *rcn;
Line 430  snextpr(ulong p, byteptr *d, long *rcn, long *q, long 
Line 536  snextpr(ulong p, byteptr *d, long *rcn, long *q, long 
         err(bugparier, "[caller of] snextpr");          err(bugparier, "[caller of] snextpr");
       }        }
     }      }
     return p + *(*d)++;      NEXT_PRIME_VIADIFF(p,*d);
       return p;
   }    }
   /* we are beyond the diffptr table */    /* we are beyond the diffptr table */
   if (*rcn == NPRC)             /* we need to initialize this now */    if (*rcn == NPRC)             /* we need to initialize this now */
Line 542  elladd0(long nbc, long nbc1,
Line 649  elladd0(long nbc, long nbc1,
 {  {
   GEN lambda;    GEN lambda;
   GEN W[2*nbcmax], *A=W+nbc;    /* W[0],A[0] never used */    GEN W[2*nbcmax], *A=W+nbc;    /* W[0],A[0] never used */
   long i, av=avma, tetpil;    long i;
     gpmem_t av=avma, tetpil;
   ulong mask = ~0UL;    ulong mask = ~0UL;
   
   /* actually, this is only ever called with nbc1==nbc or nbc1==4, so: */    /* actually, this is only ever called with nbc1==nbc or nbc1==4, so: */
Line 613  elladd2(long nbc, GEN *X1, GEN *X2, GEN *X3, GEN *X4, 
Line 721  elladd2(long nbc, GEN *X1, GEN *X2, GEN *X3, GEN *X4, 
   GEN lambda, *Y1 = X1+nbc, *Y2 = X2+nbc, *Y3 = X3+nbc;    GEN lambda, *Y1 = X1+nbc, *Y2 = X2+nbc, *Y3 = X3+nbc;
   GEN *Y4 = X4+nbc, *Y5 = X5+nbc, *Y6 = X6+nbc;    GEN *Y4 = X4+nbc, *Y5 = X5+nbc, *Y6 = X6+nbc;
   GEN W[4*nbcmax], *A=W+2*nbc;  /* W[0],A[0] never used */    GEN W[4*nbcmax], *A=W+2*nbc;  /* W[0],A[0] never used */
   long i,j, av=avma, tetpil;    long i, j;
     gpmem_t av=avma, tetpil;
   /* W[0] = gun; */    /* W[0] = gun; */
   W[1] = /* A[0] =*/ subii(X1[0], X2[0]);    W[1] = /* A[0] =*/ subii(X1[0], X2[0]);
   for (i=1; i<nbc; i++)    for (i=1; i<nbc; i++)
Line 685  elldouble(long nbc, GEN *X1, GEN *X2)
Line 794  elldouble(long nbc, GEN *X1, GEN *X2)
 {  {
   GEN lambda,v, *Y1 = X1+nbc, *Y2 = X2+nbc;    GEN lambda,v, *Y1 = X1+nbc, *Y2 = X2+nbc;
   GEN W[nbcmax+1];              /* W[0] never used */    GEN W[nbcmax+1];              /* W[0] never used */
   long i, av=avma, tetpil;    long i;
     gpmem_t av=avma, tetpil;
   /*W[0] = gun;*/ W[1] = Y1[0];    /*W[0] = gun;*/ W[1] = Y1[0];
   for (i=1; i<nbc; i++)    for (i=1; i<nbc; i++)
     W[i+1] = modii(mulii(Y1[i], W[i]), N);      W[i+1] = modii(mulii(Y1[i], W[i]), N);
Line 736  elldouble(long nbc, GEN *X1, GEN *X2)
Line 846  elldouble(long nbc, GEN *X1, GEN *X2)
 static int  static int
 ellmult(long nbc, ulong k, GEN *X1, GEN *X2) /* k>2 prime, not checked */  ellmult(long nbc, ulong k, GEN *X1, GEN *X2) /* k>2 prime, not checked */
 {  {
   long i,d,e,e1,r,av=avma,tetpil;    long i, d, e, e1, r;
     gpmem_t av=avma, tetpil;
   int res;    int res;
   GEN *A=X2, *B=XAUX, *S, *T=XAUX+2*nbc;    GEN *A=X2, *B=XAUX, *S, *T=XAUX+2*nbc;
   
Line 1039  ellfacteur(GEN n, int insist)
Line 1150  ellfacteur(GEN n, int insist)
       540400UL,606000UL,679500UL,761800UL,854100UL,957500UL,1073500UL,        540400UL,606000UL,679500UL,761800UL,854100UL,957500UL,1073500UL,
     };      };
   long nbc,nbc2,dsn,dsnmax,rep,spc,gse,gss,rcn,rcn0,bstp,bstp0;    long nbc,nbc2,dsn,dsnmax,rep,spc,gse,gss,rcn,rcn0,bstp,bstp0;
   long a,i,j,k, av,av1,avtmp, size = expi(n) + 1, tf = lgefint(n);    long a, i, j, k, size = expi(n) + 1, tf = lgefint(n);
     gpmem_t av, av1, avtmp;
   ulong B1,B2,B2_p,B2_rt,m,p,p0,p2,dp;    ulong B1,B2,B2_p,B2_rt,m,p,p0,p2,dp;
   GEN w,w0,x,*X,*XT,*XD,*XG,*YG,*XH,*XB,*XB2,*Xh,*Yh,*Xb, res = cgeti(tf);    GEN w,w0,x,*X,*XT,*XD,*XG,*YG,*XH,*XB,*XB2,*Xh,*Yh,*Xb, res = cgeti(tf);
   int rflag, use_clones = 0;    int rflag, use_clones = 0;
Line 1211  ellfacteur(GEN n, int insist)
Line 1323  ellfacteur(GEN n, int insist)
       fprintferr("\tB2 = %6lu,\tgss = %4ld*420\n", B2, gss);        fprintferr("\tB2 = %6lu,\tgss = %4ld*420\n", B2, gss);
       flusherr();        flusherr();
     }      }
     p = *d++;      p = 0;
       NEXT_PRIME_VIADIFF(p,d);
   
     /* ---B1 PHASE--- */      /* ---B1 PHASE--- */
     /* treat p=2 separately */      /* treat p=2 separately */
Line 1425  ellfacteur(GEN n, int insist)
Line 1538  ellfacteur(GEN n, int insist)
       gl = gun;        gl = gun;
       av1 = avma;        av1 = avma;
       /* scratchspace for prod (x_i-x_j) */        /* scratchspace for prod (x_i-x_j) */
       avtmp = (long)new_chunk(8 * lgefint(n));        avtmp = (gpmem_t)new_chunk(8 * lgefint(n));
       /* the correct entry in XB to use depends on bstp and on where we are        /* the correct entry in XB to use depends on bstp and on where we are
        * on the helix.  As we skip from prime to prime, bstp will be incre-         * on the helix.  As we skip from prime to prime, bstp will be incre-
        * mented by snextpr() each time we wrap around through residue class         * mented by snextpr() each time we wrap around through residue class
Line 1574  GEN
Line 1687  GEN
 pollardbrent(GEN n)  pollardbrent(GEN n)
 {  {
   long tf = lgefint(n), size = 0, delta, retries = 0, msg_mask;    long tf = lgefint(n), size = 0, delta, retries = 0, msg_mask;
   long c0, c, k, k1, l, avP, avx, GGG, av = avma;    long c0, c, k, k1, l;
     gpmem_t GGG, avP, avx, av = avma;
   GEN x, x1, y, P, g, g1, res;    GEN x, x1, y, P, g, g1, res;
   
   if (DEBUGLEVEL >= 4) (void)timer2(); /* clear timer */    if (DEBUGLEVEL >= 4) (void)timer2(); /* clear timer */
Line 1582  pollardbrent(GEN n)
Line 1696  pollardbrent(GEN n)
   if (tf >= 4)    if (tf >= 4)
     size = expi(n) + 1;      size = expi(n) + 1;
   else if (tf == 3)             /* try to keep purify happy...  */    else if (tf == 3)             /* try to keep purify happy...  */
     size = BITS_IN_LONG - bfffo(n[2]);      size = BITS_IN_LONG - bfffo((ulong)n[2]);
   
   if (size <= 28)    if (size <= 28)
     c0 = 32;                    /* amounts very nearly to `insist'.      c0 = 32;                    /* amounts very nearly to `insist'.
Line 1600  pollardbrent(GEN n)
Line 1714  pollardbrent(GEN n)
     /* 301 gives 48121 + tune_pb_min */      /* 301 gives 48121 + tune_pb_min */
     c0 = tune_pb_min + size - 60 +      c0 = tune_pb_min + size - 60 +
       ((size-73)>>1)*((size-70)>>3)*((size-56)>>4);        ((size-73)>>1)*((size-70)>>3)*((size-56)>>4);
   else    else
     c0 = 49152;                 /* ECM is faster when it'd take longer */      c0 = 49152;                 /* ECM is faster when it'd take longer */
   
   c = c0 << 5;                  /* 32 iterations per round */    c = c0 << 5;                  /* 32 iterations per round */
Line 1646  PB_RETRY:
Line 1760  PB_RETRY:
   y = cgeti(tf); affsi(2, y);    y = cgeti(tf); affsi(2, y);
   x1= cgeti(tf); affsi(2, x1);    x1= cgeti(tf); affsi(2, x1);
   avx = avma;    avx = avma;
   avP = (long)new_chunk(2 * tf); /* enough for x = addsi(tf+1) */    avP = (gpmem_t)new_chunk(2 * tf); /* enough for x = addsi(tf+1) */
   GGG = (long)new_chunk(4 * tf); /* enough for P = modii(2tf+1, tf) */    GGG = (gpmem_t)new_chunk(4 * tf); /* enough for P = modii(2tf+1, tf) */
   
   for (;;)                      /* terminated under the control of c */    for (;;)                      /* terminated under the control of c */
   {    {
Line 1826  fin:
Line 1940  fin:
 /**                 (Cf. Algorithm 8.7.2 in ACiCNT)                   **/  /**                 (Cf. Algorithm 8.7.2 in ACiCNT)                   **/
 /**                                                                   **/  /**                                                                   **/
 /***********************************************************************/  /***********************************************************************/
 extern ulong ucarrecomplet(ulong A);  
 static long squfof_ambig(long a, long B, long dd, GEN D, long *cntamb);  static long squfof_ambig(long a, long B, long dd, GEN D, long *cntamb);
   
 #define SQUFOF_BLACKLIST_SZ 64  #define SQUFOF_BLACKLIST_SZ 64
Line 1837  squfof(GEN n, long quiet)
Line 1950  squfof(GEN n, long quiet)
   long tf = lgefint(n), nm4, cnt = 0, cntamb;    long tf = lgefint(n), nm4, cnt = 0, cntamb;
   long a1, b1, c1, d1, dd1, L1, a2, b2, c2, d2, dd2, L2, a, q, c, qc, qcb;    long a1, b1, c1, d1, dd1, L1, a2, b2, c2, d2, dd2, L2, a, q, c, qc, qcb;
   GEN D1, D2, Q, res;    GEN D1, D2, Q, res;
   long av = avma;    gpmem_t av = avma;
   static long blacklist1[SQUFOF_BLACKLIST_SZ], blacklist2[SQUFOF_BLACKLIST_SZ];    static long blacklist1[SQUFOF_BLACKLIST_SZ], blacklist2[SQUFOF_BLACKLIST_SZ];
   long blp1 = 0, blp2 = 0;    long blp1 = 0, blp2 = 0;
   long mydebug = DEBUGLEVEL - quiet;    long mydebug = DEBUGLEVEL - quiet;
Line 1846  squfof(GEN n, long quiet)
Line 1959  squfof(GEN n, long quiet)
   if (cmpis(n,5) <= 0) return NULL; /* input n <= 5 */    if (cmpis(n,5) <= 0) return NULL; /* input n <= 5 */
   
 #ifdef LONG_IS_64BIT  #ifdef LONG_IS_64BIT
   if (tf > 3 || (tf == 3 && bfffo(n[2]) < 5)) /* n too large */    if (tf > 3 || (tf == 3 && bfffo((ulong)n[2]) < 5)) /* n too large */
     return NULL;      return NULL;
 #else  /* 32 bits */  #else  /* 32 bits */
   if (tf > 4 || (tf == 4 && bfffo(n[2]) < 5)) /* n too large */    if (tf > 4 || (tf == 4 && bfffo((ulong)n[2]) < 5)) /* n too large */
     return NULL;      return NULL;
 #endif  #endif
   /* now we have 5 < n < 2^59 */    /* now we have 5 < n < 2^59 */
Line 2288  static
Line 2401  static
 long  long
 squfof_ambig(long a, long B, long dd, GEN D, long *cntamb)  squfof_ambig(long a, long B, long dd, GEN D, long *cntamb)
 {  {
   long b, c, q, qc, qcb, av = avma;    long b, c, q, qc, qcb;
     gpmem_t av = avma;
   long a0, b0, b1, c0;    long a0, b0, b1, c0;
   
   q = (dd + (B>>1)) / a; qc = q*a; qcb = qc - B;    q = (dd + (B>>1)) / a; qc = q*a; qcb = qc - B;
Line 2490  static ulong powersmod[106] = {
Line 2604  static ulong powersmod[106] = {
 long                            /* no longer static -- used in mpqs.c */  long                            /* no longer static -- used in mpqs.c */
 is_odd_power(GEN x, GEN *pt, long *mask)  is_odd_power(GEN x, GEN *pt, long *mask)
 {  {
   long av=avma, tetpil, lgx=lgefint(x), exponent=0, residue, resbyte;    long lgx=lgefint(x), exponent=0, residue, resbyte;
     gpmem_t av=avma, tetpil;
   GEN y;    GEN y;
   
   *mask &= 7;                   /* paranoia */    *mask &= 7;                   /* paranoia */
Line 3325  ifac_whoiswho(GEN *partial, GEN *where, long after_cra
Line 3440  ifac_whoiswho(GEN *partial, GEN *where, long after_cra
       continue;        continue;
     }      }
     scan[2] =      scan[2] =
       (isprime((GEN)(*scan)) ?        (BSW_psp((GEN)(*scan)) ?
        (larger_compos ? un : deux) : /* un- or finished prime */         (larger_compos ? un : deux) : /* un- or finished prime */
        zero);                   /* composite */         zero);                   /* composite */
   
Line 3415  ifac_divide(GEN *partial, GEN *where)
Line 3530  ifac_divide(GEN *partial, GEN *where)
 }  }
   
   
 GEN mpqs(GEN N);                /* in src/modules/mpqs.c, maybe a dummy,  
                                  * returns a factor, or a vector of factors,  
                                  * or NULL  
                                  */  
   
 /* The following takes the place of 2.0.9.alpha's find_factor(). */  /* The following takes the place of 2.0.9.alpha's find_factor(). */
   
 /* The meaning of the hint changes against 2.0.9.alpha to:  /* The meaning of the hint changes against 2.0.9.alpha to:
Line 3440  GEN mpqs(GEN N);  /* in src/modules/mpqs.c, maybe a du
Line 3550  GEN mpqs(GEN N);  /* in src/modules/mpqs.c, maybe a du
 static long  static long
 ifac_crack(GEN *partial, GEN *where)  ifac_crack(GEN *partial, GEN *where)
 {  {
   long hint, cmp_res, exp1 = 1, exp2 = 1, av;    long hint, cmp_res, exp1 = 1, exp2 = 1;
     gpmem_t av;
   GEN factor = NULL, exponent;    GEN factor = NULL, exponent;
   
   if (DEBUGLEVEL >= 5)          /* none of these should ever happen */    if (DEBUGLEVEL >= 5)          /* none of these should ever happen */
Line 3486  ifac_crack(GEN *partial, GEN *where)
Line 3597  ifac_crack(GEN *partial, GEN *where)
   } /* while carrecomplet */    } /* while carrecomplet */
   
   /* check whether our composite hasn't become prime */    /* check whether our composite hasn't become prime */
   if (exp1 > 1 && hint != 15 && isprime((GEN)(**where)))    if (exp1 > 1 && hint != 15 && BSW_psp((GEN)(**where)))
   {    {
     (*where)[2] = un;      (*where)[2] = un;
     if (DEBUGLEVEL >= 4)      if (DEBUGLEVEL >= 4)
Line 3527  ifac_crack(GEN *partial, GEN *where)
Line 3638  ifac_crack(GEN *partial, GEN *where)
       if (moebius_mode) return 0; /* no need to carry on... */        if (moebius_mode) return 0; /* no need to carry on... */
     } /* while is_odd_power */      } /* while is_odd_power */
   
     if (exp2 > 1 && hint != 15 && isprime((GEN)(**where)))      if (exp2 > 1 && hint != 15 && BSW_psp((GEN)(**where)))
     { /* Something nice has happened and our composite has become prime */      { /* Something nice has happened and our composite has become prime */
       (*where)[2] = un;        (*where)[2] = un;
       if (DEBUGLEVEL >= 4)        if (DEBUGLEVEL >= 4)
Line 3949  ifac_primary_factor(GEN *partial, long *exponent)
Line 4060  ifac_primary_factor(GEN *partial, long *exponent)
  */   */
   
 #define ifac_overshoot 64       /* lgefint(n)+64 words reserved */  #define ifac_overshoot 64       /* lgefint(n)+64 words reserved */
 /* ifac_decomp_break:  /* ifac_decomp_break:
  *   *
  * Find primary factors of n until ifac_break return true, or n is   * Find primary factors of n until ifac_break return true, or n is
  * factored if ifac_break is NULL.   * factored if ifac_break is NULL.
Line 3968  long
Line 4079  long
 ifac_decomp_break(GEN n, long (*ifac_break)(GEN n,GEN pairs,GEN here,GEN state),  ifac_decomp_break(GEN n, long (*ifac_break)(GEN n,GEN pairs,GEN here,GEN state),
                   GEN state, long hint)                    GEN state, long hint)
 {  {
   long tf=lgefint(n), av=avma, lim=stack_lim(av,1);    long tf=lgefint(n);
     gpmem_t av=avma, lim=stack_lim(av, 1);
   long nb=0;    long nb=0;
   GEN part, here, workspc = new_chunk(tf + ifac_overshoot), pairs = (GEN)av;    GEN part, here, workspc = new_chunk(tf + ifac_overshoot), pairs = (GEN)av;
   /* workspc will be doled out by us in pairs of smaller t_INTs */    /* workspc will be doled out by us in pairs of smaller t_INTs */
   long tetpil = avma;           /* remember head of workspc zone */    gpmem_t tetpil = avma;                /* remember head of workspc zone */;
   
   if (!n || typ(n) != t_INT) err(typeer, "ifac_decomp");    if (!n || typ(n) != t_INT) err(typeer, "ifac_decomp");
   if (!signe(n) || tf < 3) err(talker, "factoring 0 in ifac_decomp");    if (!signe(n) || tf < 3) err(talker, "factoring 0 in ifac_decomp");
Line 3996  ifac_decomp_break(GEN n, long (*ifac_break)(GEN n,GEN 
Line 4108  ifac_decomp_break(GEN n, long (*ifac_break)(GEN n,GEN 
       workspc = new_chunk(lf + 3 + ifac_overshoot);        workspc = new_chunk(lf + 3 + ifac_overshoot);
       ifac_realloc(&part, &here, 0);        ifac_realloc(&part, &here, 0);
       here = ifac_find(&part, &part);        here = ifac_find(&part, &part);
       tetpil = (long)workspc;        tetpil = (gpmem_t)workspc;
     }      }
     /* room enough now */      /* room enough now */
     nb++;      nb++;
Line 4021  ifac_decomp_break(GEN n, long (*ifac_break)(GEN n,GEN 
Line 4133  ifac_decomp_break(GEN n, long (*ifac_break)(GEN n,GEN 
       part = gerepileupto(tetpil, part);        part = gerepileupto(tetpil, part);
     }      }
   }    }
   avma = (long)pairs;    avma = (gpmem_t)pairs;
   if (DEBUGLEVEL >= 3)    if (DEBUGLEVEL >= 3)
   {    {
     fprintferr("IFAC: found %ld large prime (power) factor%s.\n",      fprintferr("IFAC: found %ld large prime (power) factor%s.\n",
Line 4040  ifac_decomp(GEN n, long hint)
Line 4152  ifac_decomp(GEN n, long hint)
 long  long
 ifac_moebius(GEN n, long hint)  ifac_moebius(GEN n, long hint)
 {  {
   long mu=1, av=avma, lim=stack_lim(av,1);    long mu=1;
     gpmem_t av=avma, lim=stack_lim(av, 1);
   GEN part = ifac_start(n, 1, hint);    GEN part = ifac_start(n, 1, hint);
   GEN here = ifac_main(&part);    GEN here = ifac_main(&part);
   
Line 4065  ifac_moebius(GEN n, long hint)
Line 4178  ifac_moebius(GEN n, long hint)
 long  long
 ifac_issquarefree(GEN n, long hint)  ifac_issquarefree(GEN n, long hint)
 {  {
   long av=avma, lim=stack_lim(av,1);    gpmem_t av=avma, lim=stack_lim(av, 1);
   GEN part = ifac_start(n, 1, hint);    GEN part = ifac_start(n, 1, hint);
   GEN here = ifac_main(&part);    GEN here = ifac_main(&part);
   
Line 4089  ifac_issquarefree(GEN n, long hint)
Line 4202  ifac_issquarefree(GEN n, long hint)
 long  long
 ifac_omega(GEN n, long hint)  ifac_omega(GEN n, long hint)
 {  {
   long omega=0, av=avma, lim=stack_lim(av,1);    long omega=0;
     gpmem_t av=avma, lim=stack_lim(av, 1);
   GEN part = ifac_start(n, 0, hint);    GEN part = ifac_start(n, 0, hint);
   GEN here = ifac_main(&part);    GEN here = ifac_main(&part);
   
Line 4112  ifac_omega(GEN n, long hint)
Line 4226  ifac_omega(GEN n, long hint)
 long  long
 ifac_bigomega(GEN n, long hint)  ifac_bigomega(GEN n, long hint)
 {  {
   long Omega=0, av=avma, lim=stack_lim(av,1);    long Omega=0;
     gpmem_t av=avma, lim=stack_lim(av, 1);
   GEN part = ifac_start(n, 0, hint);    GEN part = ifac_start(n, 0, hint);
   GEN here = ifac_main(&part);    GEN here = ifac_main(&part);
   
Line 4136  GEN
Line 4251  GEN
 ifac_totient(GEN n, long hint)  ifac_totient(GEN n, long hint)
 {  {
   GEN res = cgeti(lgefint(n));    GEN res = cgeti(lgefint(n));
   long exponent, av=avma, tetpil, lim=stack_lim(av,1);    long exponent;
     gpmem_t av=avma, tetpil, lim=stack_lim(av, 1);
   GEN phi = gun;    GEN phi = gun;
   GEN part = ifac_start(n, 0, hint);    GEN part = ifac_start(n, 0, hint);
   GEN here = ifac_main(&part);    GEN here = ifac_main(&part);
Line 4186  ifac_numdiv(GEN n, long hint)
Line 4302  ifac_numdiv(GEN n, long hint)
    * size here     * size here
    */     */
   GEN res;    GEN res;
   long av=avma, tetpil, lim=stack_lim(av,1);    gpmem_t av=avma, tetpil, lim=stack_lim(av, 1);
   GEN exponent, tau = gun;    GEN exponent, tau = gun;
   GEN part = ifac_start(n, 0, hint);    GEN part = ifac_start(n, 0, hint);
   GEN here = ifac_main(&part);    GEN here = ifac_main(&part);
Line 4220  ifac_sumdiv(GEN n, long hint)
Line 4336  ifac_sumdiv(GEN n, long hint)
 {  {
   /* don't preallocate */    /* don't preallocate */
   GEN res;    GEN res;
   long exponent, av=avma, tetpil, lim=stack_lim(av,1);    long exponent;
     gpmem_t av=avma, tetpil, lim=stack_lim(av, 1);
   GEN contrib, sigma = gun;    GEN contrib, sigma = gun;
   GEN part = ifac_start(n, 0, hint);    GEN part = ifac_start(n, 0, hint);
   GEN here = ifac_main(&part);    GEN here = ifac_main(&part);
Line 4260  ifac_sumdivk(GEN n, long k, long hint)
Line 4377  ifac_sumdivk(GEN n, long k, long hint)
 {  {
   /* don't preallocate */    /* don't preallocate */
   GEN res;    GEN res;
   long exponent, av=avma, tetpil, lim=stack_lim(av,1);    long exponent;
     gpmem_t av=avma, tetpil, lim=stack_lim(av, 1);
   GEN contrib, q, sigma = gun;    GEN contrib, q, sigma = gun;
   GEN part = ifac_start(n, 0, hint);    GEN part = ifac_start(n, 0, hint);
   GEN here = ifac_main(&part);    GEN here = ifac_main(&part);
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>