Problem 233

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[0] = 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 = 0; c <= w; c++) multi[c] = 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] = 0;
    }
  for (int c = 1; c <= w; c++) multi[c] += multi[c-1];

  long long sum = 0;

  // 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 などのことを完全に忘れていた.
どうも,答えがあわない.で,いろいろと悩んだあげくに辿り着いた.
まぁ,詰めが甘いですな.