C++ PrimeSeq::next方法代码示例

bool ProvePrime(const ZZ& _n) {
  ZZ n(_n);
  if (n<0)
    abs(n,n);
  if (n<=1)
    return 0;

  if (n<=1000000) {
    // n is small so use trial division to check primality
    long ln = to_long(n);
    long end = to_long(SqrRoot(n));
    PrimeSeq s;
    for (long p=s.next(); p<=end; p=s.next())
      if ((ln%p)==0)
	return 0;
    return 1;
  }

  // check small primes
  PrimeSeq s;
  for (long p=s.next(); p<1000; p=s.next())
    if (divide(n,p))
      return 0;

  // obviously, something is missing here!

  return ProbPrime(n);
}

long PrimeTest(const ZZ& n, long t)
{
   if (n <= 1) return 0;

   PrimeSeq s;  
   long p;

   p = s.next();
   while (p && p < 2000) {
      if ((n % p) == 0) return (n == p);
      p = s.next();  
   }

   ZZ x;
   long i;

   for (i = 0; i < t; i++)
   {
      x = RandomBnd(n); 

      if (witness(n, x)) 
         return 0;
   }

   return 1;
}

bool RSA::isPrime(const ZZ& prime, long t)
{

  if (prime <= 1) 
		return false;

	// First, just try out the division by the first 2000 primes
	// Source: http://www.shoup.net/ntl/doc/tour-ex1.html
	PrimeSeq primeSequence; 
	long tempPrime;

	tempPrime = primeSequence.next();
	while(tempPrime && tempPrime < 2000) 
	{
		if((prime % tempPrime) == 0) 
			return prime == tempPrime;
		tempPrime = primeSequence.next();
	}

	// Not it's time for real prime testing
	ZZ x;
	long i;

	for(i = 0; i < t; i++) 
	{
		x = bigRandom(this->numberOfBits / 4); // random number between 0 and n-1

		if (MRTest(prime, x)) 
			return false;
	}

	return true;
}

// trial division primitive
void TrialDivision(vec_pair_ZZ_long& factors, ZZ& q, const ZZ& n, long bnd) {
  factors.SetLength(0);

  if (&q!=&n) q=n;
  if (bnd==0) {
    bnd=10000;  // should probably be higher
  }

  PrimeSeq s;
  ZZ d;
  for (long p=s.next(); (p>0 && p<=bnd); p=s.next()) {
    if (DivRem(d,q,p)==0) {
      long e=1;
      q=d;
      while (DivRem(d,q,p)==0) {
	++e;
	q=d;
      }
      addFactor(factors,to_ZZ(p),e);
      if (IsOne(q))
	return;
    }
    if (d<=p) {
      // q must be prime
      addFactor(factors,q);
      set(q);
      return;
    }
  }
}

int main() {
    PrimeSeq gen;
    for (int i = 1; i < P; ++i)
        primes[i] = gen.next();
    print("{}\n", gen.next());

    dfs(1, 1, 1, ZZ(1));
    print("ans = {}\n", ans);

    return 0;
}

// q is largest prime factor of m
// return true if q>=threshold and prime
bool check_for_factor(ZZ& q, const ZZ& m, const ZZ& threshold) {
  static const long bound=100000;

  // m must be composite
  if (ProbPrime(m))
    return false;

  // remove small factors
  q=m;
  ZZ t;
  PrimeSeq s;
  for (long p=s.next(); p<bound; p=s.next())
    while (DivRem(t,q,p)==0)
      q=t;

  if (q<threshold)
    return false;

  // lower bound on size of smallest factor
  // we are assuming that ECM finds factors smallest first (more or less)
  ZZ small_factor;
  conv(small_factor,bound);

  do {
    if (ProbPrime(q))
      return true;

    // q is composite and has a factor greater than small_factor,
    // thus it must be at least threshold*small_factor to have a 
    // large factor >threshold
    if (q<threshold*small_factor)
      return false;

    // attempt to factor with ECM
    // we are looking for factor such that
    //    small_factor < factor < q/threshold
    ECM(t,q,to_ZZ(NTL_MAX_LONG),0.5,true);

    if (IsOne(t))
      return false;

    if (t<threshold) {
      // new lower bound on remaining factors
      if (small_factor<t)
	small_factor=t;
      q/=t;
    }
    else {
      // this will be rare, but we should update small_factor anyway
      q=t;
    }

  } while (true);
}

 void main() {
     PrimeSeq gen;
     for (int u; (u = gen.next()) < n; )
         pi[u] = 1;
     for (int i = 1; i <= n; ++i) {
         pi[i] += pi[i - 1];
         if (i % 1000000 == 0) {
             print("current = {}\n", i);
         }
         for (int w = i, cnt = 0; w; w = pi[w]) {
             cnt += pi[w] == pi[w - 1];
             if (w != i) {
                 if (cnt == 100) {
                     print("XXX\n");
                     return;
                 }
                 P[cnt] += 1;
             }
         }
     }
     long ans = 1;
     for (int i = 0; i <= 100; ++i) {
         if (P[i] == 0) {
             break;
         }
         ans = ans * (P[i] % MOD) % MOD;
     }
     print("ans = {}\n", ans);
 }

void factorize(Vec< Pair<long, long> > &factors, long N)
{
  factors.SetLength(0);

  if (N < 2) return;

  PrimeSeq s;
  long n = N;
  while (n > 1) {
    if (ProbPrime(n)) {
      append(factors, cons(n, 1L));
      return;
    }

    long p = s.next();
    if ((n % p) == 0) {
      long e = 1;
      n = n/p;
      while ((n % p) == 0) {
        n = n/p;
        e++;
      }
      append(factors, cons(p, e));
    }
  }
}

// returns a non-trivial factor q of ZZ_p::modulus(), or 1 otherwise
void ECM_stage_one(ZZ& q, EC_p& Q, PrimeSeq& seq, long bound) {
  long sbound = (long)sqrt((double)bound);
  seq.reset(0);
  long p = seq.next();
  for (; p<=sbound; p=seq.next()) {
    long pp,t=p;
    do { pp=t; t*=p; } while (t>pp && t<=bound); // we might overflow t here
    mul(Q,Q,pp);
  }
  for (; p<=bound; p=seq.next())
    mul(Q,Q,p);
  if (!IsZero(Q))
    GCD(q,ZZ_p::modulus(),rep(Q.Z));
  else
    set(q);
}

int step0(ZZ n){

    long p;

    PrimeSeq s;
    
    p = s.next();  
    while (p && p < 2000) {
        if ((n % p) == 0){
            cout << p << " divides " << n << '\n';
            return (1 + (n == p)); // 1: composite
        }                          // 2: prime
        p = s.next();              
    }
    return 0; // 0: continue
}

// Returns a list of prime factors and their multiplicity, 
// N = \prod_i factors[i].first^{factors[i].second}
void factorize(Vec< Pair<long, long> > &factors, long N)
{
  factors.SetLength(0);

  if (N < 2) return;

  PrimeSeq s;
  long n = N;
  while (n > 1) {
    if (ProbPrime(n)) { // n itself is a prime, add (n,1) to the list
      append(factors, cons(n, 1L));
      return;
    }

    long p = s.next();
    if ((n % p) == 0) { // p divides n, find its multiplicity
      long e = 1;
      n = n/p;
      while ((n % p) == 0) {
        n = n/p;
        e++;
      }
      append(factors, cons(p, e)); // add (p,e) to the list
    }
  }
}

template<class zp,class zz> void FindPrimRootT(zp &root, unsigned long e)
{
  zz qm1 = zp::modulus()-1;

  assert(qm1 % e == 0);
  
  vector<long> facts;
  factorize(facts,e); // factorization of e

  root = 1;

  for (unsigned long i = 0; i < facts.size(); i++) {
    long p = facts[i];
    long pp = p;
    long ee = e/p;
    while (ee % p == 0) {
      ee = ee/p;
      pp = pp*p;
    }
    // so now we have e = pp * ee, where pp is 
    // the power of p that divides e.
    // Our goal is to find an element of order pp

    PrimeSeq s;
    long q;
    zp qq, qq1;
    long iter = 0;
    do {
      iter++;
      if (iter > 1000000) 
        Error("FindPrimitiveRoot: possible infinite loop?");
      q = s.next();
      conv(qq, q);
      power(qq1, qq, qm1/p);
    } while (qq1 == 1);
    power(qq1, qq, qm1/pp); // qq1 has order pp

    mul(root, root, qq1);
  }

  // independent check that we have an e-th root of unity 
  {
    zp s;

    power(s, root, e);
    if (s != 1) Error("FindPrimitiveRoot: internal error (1)");

    // check that s^{e/p} != 1 for any prime divisor p of e
    for (unsigned long i=0; i<facts.size(); i++) {
      long e2 = e/facts[i];
      power(s, root, e2);   // s = root^{e/p}
      if (s == 1) 
        Error("FindPrimitiveRoot: internal error (2)");
    }
  }
}

/* Compute Phi(N) */
long phi_N(long N)
{
  long phiN=1,p,e;
  PrimeSeq s;
  while (N!=1)
    { p=s.next();
      e=0;
      while ((N%p)==0) { N=N/p; e++; }
      if (e!=0)
        { phiN=phiN*(p-1)*power_long(p,e-1); }
    }
  return phiN;
}

bool ProvePrime(const ZZ& N, vec_vec_ZZ& certs) {
  certs.SetLength(0);
  long i=0;
  ZZ num(N);
  HCP_generate HCP;
  while (NumBits(num)>30) {
    certs.SetLength(i+1);
    if (!ProvePrime_Atkin(num,certs[i],HCP))
      return false;
    num=certs[i][1];
    ++i;
  }
  // prove small prime by trial division
  long n,limit;
  conv(limit,SqrRoot(num));
  conv(n,num);
  PrimeSeq s;
  for (long p=s.next(); p<=limit; p=s.next())
    if (n%p==0)
      return false;
  return true;
}

int main(){
        long n = 183783600;
        long b = 18;
        long rootn=SqrRoot(n);
        long p;
        PrimeSeq s;
        vector<long> values;
        vector<long> logV;
        for(long m=n+1;m<1.66*n;m++){
                values.push_back(m);
                logV.push_back(log(m));
        }
                s.reset(b);
                p=s.next();
                int temp=log(p);
                while(p<rootn){
        for(int i=0;i<logV.size();i++){
                        long t=pow(10,(logV[i]/temp));
                        if(ceil(t)-t<0.05 || t-floor(t)<0.05){
//                        cout<<"here at "<<i<<endl;
                                values[i]=0;
                                
                                break;
                        }
                }
                        p=s.next();
        }
        int y=0;
        cout<<"arr=[";
        for(int i=0;i<values.size();++i){
                if(values[i]!=0)cout<<values[i]<<",";
//                if(temp[i]!=0)cout<<temp[i]<<endl;
        }
        cout<<"0]"<<endl;
        cout<<"ND = numdiv("<<n<<"); nd = 0; i = 1; while( nd<=ND, nd = numdiv(arr[i]); i = i+1); print(arr[i-1])"<<endl;

        return 0;
}