Project Euler 131

プロジェクトオイラー
http://projecteuler.net/index.php

Q131.
n3 + pn2が立方数となりうる素数pの100万以下の個数

n3 + pn2 = (n + k)3
(p-3k)n2 - 3k2n - k3 = 0
D = 4pk - 3k2

判別式が平方数だから、

4pk - 3k2 = m2
3m2 + (3k-2p)2 = 4p2

ここで挫折した。Pythonで書いて話にならなかったので、C++で書いた。32ビットの範囲で計算できないので、doubleを使って適当にごまかした。



#include
#include

using namespace std;

class gen_primes {
int p;

public:
gen_primes() { p = 2; }

int operator()() {
if(p == 2) {
p = 3;
return 2;
}
else {
while(!is_prime(p))
p += 2;
p += 2;
return p - 2;
}
}

private:
bool is_prime(int n) const { // n must be greater than 2.
for(int p = 3; p * p <= n; p += 2) {
if(n % p == 0)
return false;
}
return true;
}
};

int calc_m(int l, int q) {
double ratio = (double)q / l;
double m_f = l * sqrt((1 - ratio * ratio) / 3);

// is m_f almost integer?
int m = (int)(m_f + 0.5);
if(abs(m_f - m) < 1e-8)
return m;
else
return -1;
}

bool find_cube(int p) {
if(p == 3)
return false;

int limit = 2 * p;
for(int k = 1; k <= 2 * p / 3; k++) {
int q = 3 * k - limit;
int m = calc_m(limit, q);
if(m == -1)
continue;
double n_f = (3 * k + m) * (double)k / (2 * (p - 3 * k));
int n = (int)(n_f + 0.5);
if(abs(n_f - n) < 1e-8) {
cout << p << " " << n << endl;
return true;
}
}
return false;
}

const int N = (int)1e6;

int main() {
#if 0
find_cube(19);
#else
gen_primes g;
int counter = 0;
while(true) {
int p = g();
if(p >= N)
break;
if(find_cube(p))
counter += 1;
}
cout << counter << endl;
#endif
}