233 Lattice points on a circle

Problem 233 – Project Euler

円周上にある格子点の数についての問題．

難しい．数学的知識が必要．

ポイントは

(N-2x)^2 + （N-2y)^2 = 2N^2 という式と，

与えられた数を二つの平方数の和として表わす方法の数．

ふるいを2回使うので，haskellでは面倒．

C++で実装．

自分のデスクトップPCの環境では実行時間は 0.2 秒以下で十分速い． 

#include <iostream>
#include <vector>
#include <cmath>
#define FOR(x, xs) for (__typeof( (xs).begin() ) x = (xs).begin(); x != (xs).end(); x++)
#define WHILE(x, xs, p) for (__typeof( (xs).begin() ) x = (xs).begin(); (p) && x != (xs).end(); x++)
using namespace std;
static const long long u = pow(10, 11);
static const int v = u / ( pow(5, 3) * pow(13, 2) );
static const int w = u / ( pow(5, 3) * pow(13, 2) * 17);
int main() {
vector<bool> isPrime( v + 1); //primes up to v
FOR (q, isPrime) *q = true;
isPrime[] = false;
isPrime[1] = false;
for (int q = 2; q <= sqrt(v) ; q++)
if ( isPrime[q] )
for (int i = q * q; i <= v; i += q)
isPrime[i] = false;
vector<int> primes1; // primes form of 4k+1
vector<long long> multi( w + 1 );
for (int c = ; c <= w; c++) multi = c;
for (int i = 1; i <= v ; i++)
if ( isPrime[i] && i % 4 == 1 ) {
primes1.push_back(i);
for (int j = i; j <= w; j += i)
multi[j] = ;
}
for (int c = 1; c <= w; c++) multi += multi[c-1];
long long sum = ;
// n = p * q^2 * r^3
WHILE (r, primes1, *r <= pow( u / (13 * 5 * 5), 1.0 / 3 ) )
WHILE (q, primes1, *q <= sqrt ( u / ( 5 * pow(*r, 3) ) ) ) {
if ( *r == *q ) continue;
WHILE (p, primes1, *p <= u / ( pow(*q, 2) * pow(*r, 3) ) ) {
if ( *r == *p || *q == *p ) continue;
long long n =  *p * pow(*q, 2) * pow(*r, 3);
sum += n * multi[ u / n ];
}
}
// n = q^3 * r^7
WHILE (r, primes1, *r <= pow ( u / pow(5, 3), 1.0 / 7 ) )
WHILE (q, primes1, *q <= pow( u / pow(*r, 7), 1.0 / 3) ) {
if ( *r == *q ) continue;
long long n =   pow(*r, 7) * pow(*q, 3);
sum += n * multi[ u / n ];
}
// n = q^2 * r^10
WHILE (r, primes1, *r <= pow ( u / pow(5, 2), 0.1 ) )
WHILE (q, primes1, *q <= sqrt( u / pow(*r, 10) ) ){
if ( *r == *q ) continue;
long long n =   pow(*r, 10) * pow(*q, 2);
sum += n * multi[ u / n ];
}
cout << sum << endl;
}

はじめは，n = p*q^2*r^3 だけしか，考えておらず，

n ＝ q^3*r^7, q^2*r10 などのことを完全に忘れていた．

どうも，答えがあわない．で，いろいろと悩んだあげくに辿り着いた．

まぁ，詰めが甘いですな．