プロジェクトオイラー
http://projecteuler.net/index.php
またしてもPythonのメモリの使い方に苦しめられた。エラトステネスのふるい的な解法を試したが、多倍長整数がやたらとメモリを食う。そのため、1つの数で2つの整数を使うように配列の長さを2倍してみたのだが、それでもメモリを食う。結局どうしてメモリを食っているのか分からなかった。
C++で64ビット整数を使ったらあっさり解けた。
素数pについて、あるk0が2k02 ≡ 1(mod p)を満たせば、n = mp + k0もこれを満たす。また、k1 = p - k0も同様に満たす。よって、k0が分かれば、エラトステネスのふるいと同様のアルゴリズムが使える。
しかし、k0はすぐには分からないので、k0からpを求めるようにする。a[n] = 2n - 1となる配列を用意し、p = a[n]が素数なら、pごとにpで割っていく。というか、pは必ず素数となる。なぜなら、複数の素数からなる合成数なら、小さいほうの素数がその大きさからしてその前に現れていなければならないからである。また、2は素因数になりえないので、N = 50000000として、k0は、2√N/3まで調べればよい。
#include
#include
#includeusing namespace std;
vector
primes; bool is_prime(int n) {
if(n < 2)
return false;
typedef vector::const_iterator CIT;
for(CIT q = primes.begin(); q != primes.end(); ++q) {
int p = *q;
if(p * p > n)
return true;
else if(n % p == 0)
return false;
}
return true;
}void make_primes(int n) {
int m = (int)sqrt( (n + 0.5));
for(int p = 3; p <= m; p += 2) {
if(is_prime(p))
primes.push_back(p);
}
int size = (n + 15) / 16;
int *a = new int[size];
if(!a) {
cout << size << endl;
exit(1);
}
for(int i = 0; i < size; i++)
a[i] = 0;
typedef vector::const_iterator CIT;
for(CIT q = primes.begin(); q != primes.end(); ++q) {
int p = *q;
if(p * p > n)
break;
if(p == 2)
continue;
for(int k = 3 * p; k < n; k += 2 * p)
a[k>>4] |= 1 << (k & 15);
}
for(int k = (m + 1) / 2 * 2 + 1; k <= n; k += 2) {
if((a[k>>4] & (1 << (k & 15))) == 0)
primes.push_back(k);
}
}void init_a(long long *a, int size) {
for(int n = 0; n < size; n++)
a[n] = 2 * (long long)n * n - 1;
}void divide(long long *a, int n, long long p) {
do
a[n] /= p;
while(a[n] % p == 0);
}int main() {
const int N = (int)5e7;
const int M = (int)(sqrt(2.) * N / 3);
primes.push_back(2);
make_primes(N);
long long *a = new long long[N+1];
init_a(a, N + 1);
for(int k0 = 2; k0 <= M; k0++) {
if(a[k0] > 1) {
long long p = a[k0]; // 2k0^2 ≡ 1(mod p)となる最小のk
long long k1 = p - k0; // もう一つのk
for(long long n = k0; n <= N; n += p) {
if(n == k0 && p == 2 * (long long)k0 * k0 - 1)
continue;
divide(a, n, p);
}
for(long long n = k1; n <= N; n += p) {
divide(a, n, p);
}
}
}
int counter = 0;
for(int n = 2; n <= N; n++) {
if(a[n] == 2 * (long long)n * n - 1)
counter++;
}
cout << counter << endl;
delete [] a;
}
