素数の解法:
17271 ワード
一、素朴な判断素数アルゴリズム
素数を判断するのは、実はとても簡単です.定義によれば、1つの整数nが素数であるか否かを判断するには、整数区間[2,n−1]内にある数mがあるか否かを判断し、n%m==0とする必要がある.コードは次のように書くことができます.
実際,このアルゴリズムはO(n)であり,感覚は速いが,依然として需要を満たすことができない.したがって,O(sqrt(n))のアルゴリズムがある.コードは次のように書くことができます.
原理は巧みで,コードのi二、Miller_Rabin素性試験
上のO(sqrt(n))のアルゴリズムは非常に優れているが,より高い数級の「大数」に直面すると力不足に見える.そしてこの時、Miller_Rabinの優位性が示された.Miller_Rabinの理論的基礎はFermaの小さな定理に由来する.Fermaの小さな定理は数論の四大定理の一つであることに値する.
Fermaの小さい定理:nは1つの奇素数で、aはいかなる整数(1≦a≦n-1)で、a^(n-1)≡1(mod n)
nが素数であるかどうかをテストするには、まずn−1を(2^s)*dに分解し、テスト開始のたびに[1,N−1]の整数aをランダムに選択し、その後、すべてのr∈[0,s−1]、a^d mod N≠1およびa^((2^r)*d)mod N≠−1であれば、nは合数であり、そうでなければnが3/4の確率で素数である.
これもなぜこのアルゴリズムが素性テストにすぎず、結果が正しいことを完全に保証することはできないが、テスト回数が約20の場合、基本的に正解率が97%以上に達することを確保することができるのか.コードは次のように書くことができます.
四、選別法
上で紹介したいくつかの素数判断のアルゴリズムは,いくつかの問題で基本的に適用できるようになった.しかし、別の問題については気まずい.例えば、整数区間[1,n]内の素数の個数を聞いてみましょう.このとき一つ一つ列挙し、一つ一つ判断すれば、素朴な方法にとってはO(nsqrt(n))が最適であり、Miller_RabinアルゴリズムもO(n)の時間内に結果を得ることができなかった.そこで、エラトスターニふるい法(エラトスターニふるい法)がありました.
フィルタ整数n以内の素数については,アルゴリズムは,素数2の倍数を全て削除し,残りの数は1番目が3であり,素数3の倍数を全て削除し,残りの1番目の数は5であり,素数5の倍数を全て削除する・・・nを内最後の素数の倍数で削除し,n以内のすべての素数を得る.
コードは次のように書くことができます.
実際、観察によると、上記の書き方は何度も計算を繰り返しているので、線形フィルタリングは明らかにできないが、別の書き方は線形フィルタリングを得ることができ、時間の複雑さはO(n)に達し、コードはこのように書くことができる.
kuangbinからのテンプレート.
四、反発原理
上記のコードから、このふるい法は1 e 7という数段の演算にしか対応できないことが明らかになった.線形のふるい法でも、ACMコンテストでは1 e 8のメモリがMemery Limit Exceedを得る可能性が高いため、満足できない.そこで,反発原理を考えることができる.AHUOJ 557を例にとると,1 e 8の場合はフィルタリング法では全く満足できないが,a*b=cの場合を考慮し,1 e 8は10000以内の素数p[10000]のみを考慮し,その後毎回n/p[i]を減算し,n/(p[i]*p[j])を加えてn/(p[i]*p[j]*p[k])を減算することで,このように推すだけで…と正確な結果が得られる.
コードは次のとおりです.
五、Meissel-Lehmerアルゴリズム
最後に紹介したこのアルゴリズムはブラックテクノロジーレベルと言えるが,3 s以内に1 e 11以内の素数個数を求めることができる.300 ms以内に1 e 11の個数を求めるものもあるそうです.wikiの原理を参考にすることができます.そしてCodeforces 665 Fテーマからのコードを与えます.
素数を判断するのは、実はとても簡単です.定義によれば、1つの整数nが素数であるか否かを判断するには、整数区間[2,n−1]内にある数mがあるか否かを判断し、n%m==0とする必要がある.コードは次のように書くことができます.
int isPrime(int n) {
int i;
for (i = 2; i < n; ++i) {
if (n % i == 0) return 0;
}
return 1;
}
実際,このアルゴリズムはO(n)であり,感覚は速いが,依然として需要を満たすことができない.したがって,O(sqrt(n))のアルゴリズムがある.コードは次のように書くことができます.
int isPrime(long long n) {
long long i;
for (i = 2; i * i <= n; ++i) {
if (n % i == 0) return 0;
}
return 1;
}
原理は巧みで,コードのi
上のO(sqrt(n))のアルゴリズムは非常に優れているが,より高い数級の「大数」に直面すると力不足に見える.そしてこの時、Miller_Rabinの優位性が示された.Miller_Rabinの理論的基礎はFermaの小さな定理に由来する.Fermaの小さな定理は数論の四大定理の一つであることに値する.
Fermaの小さい定理:nは1つの奇素数で、aはいかなる整数(1≦a≦n-1)で、a^(n-1)≡1(mod n)
nが素数であるかどうかをテストするには、まずn−1を(2^s)*dに分解し、テスト開始のたびに[1,N−1]の整数aをランダムに選択し、その後、すべてのr∈[0,s−1]、a^d mod N≠1およびa^((2^r)*d)mod N≠−1であれば、nは合数であり、そうでなければnが3/4の確率で素数である.
これもなぜこのアルゴリズムが素性テストにすぎず、結果が正しいことを完全に保証することはできないが、テスト回数が約20の場合、基本的に正解率が97%以上に達することを確保することができるのか.コードは次のように書くことができます.
const int S = 20;
LL mod_mul(LL a, LL b, LL n) {
LL res = 0;
while (b) {
if (b & 1) res = (res + a) % n;
a = (a + a) % n;
b >>= 1;
}
return res;
}
LL mod_exp(LL a, LL b, LL n) {
LL res = 1;
while (b) {
if (b & 1) res = mod_mul(res, a, n);
a = mod_mul(a, a, n);
b >>= 1;
}
return res;
}
bool miller_rabin(LL n) {
if (n == 2 || n == 3 || n == 5 || n == 7 || n == 11) return true;
if (n == 1 || !(n % 2) || !(n % 3) || !(n % 5) || !(n % 7) || !(n % 11)) return false;
LL x, pre, u = n - 1, k = 0;
while (!(u & 1)) {
++k;
u >>= 1;
}
srand((LL)time(NULL));
for (int i = 0; i < S; ++i) { // S ,S ,
x = rand() % (n - 2) + 2; // [2, n)
if (x % n == 0) continue;
x = mod_exp(x, u, n); // x^u % n
pre = x;
for (int j = 0; j < k; ++j) {
x = mod_mul(x, x, n);
if (x == 1 && pre != 1 && pre != n - 1)
return false;
pre = x;
}
if (x != 1) return false;
}
return true;
}
四、選別法
上で紹介したいくつかの素数判断のアルゴリズムは,いくつかの問題で基本的に適用できるようになった.しかし、別の問題については気まずい.例えば、整数区間[1,n]内の素数の個数を聞いてみましょう.このとき一つ一つ列挙し、一つ一つ判断すれば、素朴な方法にとってはO(nsqrt(n))が最適であり、Miller_RabinアルゴリズムもO(n)の時間内に結果を得ることができなかった.そこで、エラトスターニふるい法(エラトスターニふるい法)がありました.
フィルタ整数n以内の素数については,アルゴリズムは,素数2の倍数を全て削除し,残りの数は1番目が3であり,素数3の倍数を全て削除し,残りの1番目の数は5であり,素数5の倍数を全て削除する・・・nを内最後の素数の倍数で削除し,n以内のすべての素数を得る.
コードは次のように書くことができます.
const int maxn = 1e7 + 5;
int pri[maxn];
void getPrime(int n) {
for (int i = 0; i <= n; ++i) pri[i] = i;
pri[1] = 0;
for (int i = 2; i <= n; ++i) {
if (!pri[i]) continue;
pri[++pri[0]] = i;
for (int j = 2; i * j <= n && j < n; ++j) {
pri[i * j] = 0;
}
}
}
実際、観察によると、上記の書き方は何度も計算を繰り返しているので、線形フィルタリングは明らかにできないが、別の書き方は線形フィルタリングを得ることができ、時間の複雑さはO(n)に達し、コードはこのように書くことができる.
const int MAXN = 1e7 + 5;
void getPrime(){
memset(prime, 0, sizeof(prime));
for (int i = 2;i <= MAXN;i++) {
if (!prime[i])prime[++prime[0]] = i;
for (int j = 1;j <= prime[0] && prime[j] <= MAXN / i;j++) {
prime[prime[j] * i] = 1;
if (i%prime[j] == 0) break;
}
}
}
kuangbinからのテンプレート.
四、反発原理
上記のコードから、このふるい法は1 e 7という数段の演算にしか対応できないことが明らかになった.線形のふるい法でも、ACMコンテストでは1 e 8のメモリがMemery Limit Exceedを得る可能性が高いため、満足できない.そこで,反発原理を考えることができる.AHUOJ 557を例にとると,1 e 8の場合はフィルタリング法では全く満足できないが,a*b=cの場合を考慮し,1 e 8は10000以内の素数p[10000]のみを考慮し,その後毎回n/p[i]を減算し,n/(p[i]*p[j])を加えてn/(p[i]*p[j]*p[k])を減算することで,このように推すだけで…と正確な結果が得られる.
コードは次のとおりです.
#include
#include
using namespace std;
const int maxn = 10005;
int sqrn, n, ans = 0;
bool vis[maxn];
int pri[1500] = {0};
void init(){
vis[1] = true;
int k = 0;
for(int i = 2; i < maxn; i++){
if(!vis[i]) pri[k++] = i;
for(int j = 0; j < k && pri[j] * i < maxn; j++){
vis[pri[j] * i] = true;
if(i % pri[j] == 0) break;
}
}
}
void dfs(int num, int res, int index){
for(int i = index; pri[i] <= sqrn; i++){
if(1LL * res * pri[i] > n){
return;
}
dfs(num + 1, res * pri[i], i+1);
if(num % 2 == 1){
ans -= n / (res * pri[i]);
}else{
ans += n / (res * pri[i]);
}
if(num == 1) ans++;
}
}
int main(){
init();
while(~scanf("%d",&n) && n){
ans = n;
sqrn = sqrt((double)n);
dfs(1,1,0);
printf("%d
",ans-1);
}
return 0;
}
五、Meissel-Lehmerアルゴリズム
最後に紹介したこのアルゴリズムはブラックテクノロジーレベルと言えるが,3 s以内に1 e 11以内の素数個数を求めることができる.300 ms以内に1 e 11の個数を求めるものもあるそうです.wikiの原理を参考にすることができます.そしてCodeforces 665 Fテーマからのコードを与えます.
#define MAXN 100 // pre-calc max n for phi(m, n)
#define MAXM 10010 // pre-calc max m for phi(m, n)
#define MAXP 40000 // max primes counter
#define MAX 400010 // max prime
#define setbit(ar, i) (((ar[(i) >> 6]) |= (1 << (((i) >> 1) & 31))))
#define chkbit(ar, i) (((ar[(i) >> 6]) & (1 << (((i) >> 1) & 31))))
#define isprime(x) (( (x) && ((x)&1) && (!chkbit(ar, (x)))) || ((x) == 2))
namespace pcf {
long long dp[MAXN][MAXM];
unsigned int ar[(MAX >> 6) + 5] = { 0 };
int len = 0, primes[MAXP], counter[MAX];
void Sieve() {
setbit(ar, 0), setbit(ar, 1);
for (int i = 3; (i * i) < MAX; i++, i++) {
if (!chkbit(ar, i)) {
int k = i << 1;
for (int j = (i * i); j < MAX; j += k) setbit(ar, j);
}
}
for (int i = 1; i < MAX; i++) {
counter[i] = counter[i - 1];
if (isprime(i)) primes[len++] = i, counter[i]++;
}
}
void init() {
Sieve();
for (int n = 0; n < MAXN; n++) {
for (int m = 0; m < MAXM; m++) {
if (!n) dp[n][m] = m;
else dp[n][m] = dp[n - 1][m] - dp[n - 1][m / primes[n - 1]];
}
}
}
long long phi(long long m, int n) {
if (n == 0) return m;
if (primes[n - 1] >= m) return 1;
if (m < MAXM && n < MAXN) return dp[n][m];
return phi(m, n - 1) - phi(m / primes[n - 1], n - 1);
}
long long Lehmer(long long m) {
if (m < MAX) return counter[m];
long long w, res = 0;
int i, a, s, c, x, y;
s = sqrt(0.9 + m), y = c = cbrt(0.9 + m);
a = counter[y], res = phi(m, a) + a - 1;
for (i = a; primes[i] <= s; i++) res = res - Lehmer(m / primes[i]) + Lehmer(primes[i]) - 1;
return res;
}
}