All pastes #454763 Raw Edit

Anonymous

public text v1 · immutable
·published 2007-04-23 20:17 UTC
rendered paste body
static GEN
auxdecomp1(GEN n, long (*ifac_break)(GEN n, GEN pairs, GEN here, GEN state),
                  GEN state, ulong all, long hint)
{
  pari_sp av;
  long pp[] = { evaltyp(t_INT)|_evallg(4), 0,0,0 };
  long nb = 0, i, lp;
  ulong p, k, lim;
  byteptr d = diffptr+1; /* start at p = 3 */
#define STORE(x,e) { nb++; (void)x; (void)utoipos(e); }
#define STOREu(x,e) STORE(utoipos(x), e)
#define STOREi(x,e) STORE(icopy(x),   e)

  if (typ(n) != t_INT) pari_err(arither1);
  i = signe(n); if (!i) pari_err(talker, "zero argument in factorint");
  (void)cgetg(3,t_MAT);
  if (i < 0) STORE(utoineg(1), 1);
  if (is_pm1(n)) return aux_end(NULL,nb);

  n = gclone(n); setsigne(n,1);
  i = vali(n);
  if (i)
  {
    STOREu(2, i);
    av = avma; affii(shifti(n,-i), n); avma = av;
  }
  if (is_pm1(n)) return aux_end(n,nb);
  lim = tridiv_bound(n, all);

  /* trial division */
  p = 2;
  for(;;)
  {
    int stop;
    NEXT_PRIME_VIADIFF(p,d);
    if (p >= lim) break;

    k = Z_lvalrem_stop(n, p, &stop);
    if (k) STOREu(p, k);
    if (stop)
    {
      if (!is_pm1(n)) STOREi(n, 1);
      return aux_end(n,nb);
    }
  }

  /* pp = square of biggest p tried so far */
  av = avma; affii(muluu(p,p), pp); avma = av;

  /* trial divide by the "special primes" (usually huge composites) */
  lp = lg(primetab);
  for (i=1; i<lp; i++)
    if (dvdiiz(n,gel(primetab,i), n))
    {
      long k = 1; while (dvdiiz(n,gel(primetab,i), n)) k++;
      STOREi(gel(primetab,i), k);
      if (absi_cmp(pp, n) > 0)
      {
        if (!is_pm1(n)) STOREi(n, 1);
        return aux_end(n,nb);
      }
    }

  if (all != 1) hint = 15; /* smallfact: turn off all except pure powers */
  else /* else test primality */
    if (BSW_psp(n)) { STOREi(n, 1); return aux_end(n,nb); }

  /* now we have a large composite */
  if (ifac_break && (*ifac_break)(n,NULL,NULL,state)) /*initialize ifac_break*/
  {
    if (DEBUGLEVEL>2)
      fprintferr("IFAC: (Partial fact.) Initial stop requested.\n");
  }
  else
    nb += ifac_decomp_break(n, ifac_break, state, hint);

  return aux_end(n, nb);
}