プロジェクトオイラー
http://projecteuler.net/
Q223.
三角形の辺a≤b≤cが、a2 + b2 = c2 + 1を満たすとき、ぎりぎり鋭角三角形と呼ぶ。周が2500万以下のぎりぎり鋭角三角形はいくつあるか。
特に工夫もない。aについてエラトステネスのふるい的に素因数分解する。そして、(a - 1)(a + 1)の約数をすべて出す。そして条件を満たすものを数えていく。ただし、メモリが厳しいので、素因数分解は、素数と指数で一つの整数として格納した。また、(a - 1)(a + 1)が偶数・奇数で分けて、別々に計算した。それでも足りないので、結局C++で組みなおした。
#include
#include
#includeusing namespace std;
vector
primes; bool is_prime(int n) {
typedef vector::const_iterator PIT;
for(PIT 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 *a = new int[(n + 15) / 16];
for(int i = 0; i < (n + 15) / 16; i++)
a[i] = 0;
typedef vector::const_iterator PIT;
for(PIT 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);
}
delete [] a;
}vector
*factorize(int M) {
int *a = new int[M+1];
vector*b = new vector [M+1];
for(int n = 2; n <= M; n++) {
a[n] = n;
}
typedef vector::const_iterator PIT;
for(PIT q = primes.begin(); q != primes.end(); ++q) {
int p = *q;
long long p_pow = p;
int e = 1;
while(p_pow <= M) {
for(int n = p_pow; n <= M; n += p_pow) {
a[n] /= p;
if(e == 1)
b[n].push_back( (p << 5) + e);
else
b[n].back() += 1;
}
p_pow *= p;
e += 1;
}
}
for(int n = 2; n <= M; n++) {
if(a[n] != 1)
b[n].push_back( (a[n] << 5) + 1);
}
delete [] a;
return b;
}vector
*factorize_odd(int M) {
int *a = new int[M+1];
vector*b = new vector [M+1];
for(int n = 2; n <= M; n++) {
a[n] = n * 2 + 1;
}
typedef vector::const_iterator PIT;
for(PIT q = primes.begin(); q != primes.end(); ++q) {
int p = *q;
if(p == 2)
continue;
long long p_pow = p;
int e = 1;
while(p_pow <= M * 2 + 1) {
for(int n = (p_pow - 1) / 2; n <= M; n += p_pow) {
a[n] /= p;
if(e == 1)
b[n].push_back( (p << 5) + e);
else
b[n].back() += 1;
}
p_pow *= p;
e += 1;
}
}
for(int n = 1; n <= M; n++) {
if(a[n] != 1)
b[n].push_back( (a[n] << 5) + 1);
}
delete [] a;
return b;
}long long pow(int n, int e) {
if(e == 0)
return 1;
long long m = pow(n, e / 2);
m *= m;
if(e & 1)
return m * n;
else
return m;
}class cGenDivisors {
const vector*v;
int n;
int k;
public:
cGenDivisors(const vector& a) {
this->v = &a;
n = 1;
for(int i = 0; i < (int)v->size(); i++)
n *= (*(v->begin() + i) & 31) + 1;
k = 0;
}
long long operator ()() {
if(k == n)
return 0;
long long m = 1;
int l = k;
for(int i = 0; i < (int)v->size(); i++) {
int p = *(v->begin() + i) >> 5;
int q = (*(v->begin() + i) & 31) + 1;
int e = l % q;
m *= pow(p, e);
l /= q;
}
k++;
return m;
}
};void vec_add(vector
& v, const vector & a, const vector & b) {
v = a;
for(vector::const_iterator p = b.begin(); p != b.end(); ++p)
v.push_back(*p);
}const int N = (int)2.5e7;
int num_solutions(const vector
*b, int k) {
int counter = 0;
long long a = 2 * k + 1;
long long n = a * a - 1;
vectorv;
vec_add(v, b[k], b[k+1]);
cGenDivisors g(v);
while(true) {
long long d1 = g();
if(d1 == 0)
break;
long long d2 = (n / d1) / 4;
if(d2 * 2 <= N - a) {
long long b1 = d2 - d1;
if(a <= b1)
counter++;
}
}
return counter;
}int num_solutions_odd(const vector
*b, int k) {
int counter = 0;
long long a = 2 * k;
long long n = a * a - 1;
vectorv;
vec_add(v, b[k-1], b[k]);
cGenDivisors g(v);
while(true) {
long long d1 = g();
if(d1 == 0)
break;
long long d2 = n / d1;
if(a + d2 <= N) {
long long b1 = (d2 - d1) / 2;
if(a <= b1)
counter++;
}
}
return counter;
}int count_barely_acute_triangles() {
const int M = (N - 3) / 4;
const int L = (N - 1) / 4;
make_primes( (int)sqrt(M * 2 + 1.5));
int counter = (N - 1) / 2;
vector*b = factorize(M);
for(int k = 1; k < M; k++)
counter += num_solutions(b, k);
delete [] b;
b = factorize_odd(L);
for(int k = 1; k < L + 1; k++)
counter += num_solutions_odd(b, k);
delete [] b;
return counter;
}int main() {
primes.push_back(2);
cout << count_barely_acute_triangles() << endl;
}
