Problem 229

229 Four Representations using Squares
Problem 229 - Project Euler
この問題は難しかった．
なんてったって，問題サイズが 2*10^9 と非常に大きい．
詳しくは，明日ぐらいに追記する．
遅いほう．Haskell, Java で実装したらメモリが足りなくて，うわーん，ってなった．
力技．マシンパワーに頼る作戦．

#include <iostream>
#include <vector>
#include <cmath>
#define FOR(x, xs) for (__typeof( (xs).begin() ) x = (xs).begin(); x != (xs).end(); x++)

using namespace std;

static const int n = 2000000000;
static const int u = sqrt(n);

void sieve(int k, vector<bool> &s, vector<bool> &t) {
  FOR(x, s) *x = false;
  for (int x = 1; x <= u; x++)
    for (int y = 1, u_y = sqrt( ( n - x * x ) / k); y <= u_y; y++)
      if( t[ x * x + k * y * y] )
	s[ x * x + k * y * y ] = true;
}

int main() {
  vector<bool> set1( n + 1 );
  vector<bool> set2( n + 1 );
  FOR(p, set2) *p = true;

  sieve(7, set1, set2);
  sieve(3, set2, set1);
  sieve(2, set1, set2);
  sieve(1, set2, set1);

  int count = 0;
  for (int i =1; i < n; i++)
    if( set2[i] ) count++;
  cout << count << endl;
}

速いほう．平方数とそれ以外で場合分けして考えれば，良いことに気がついた．

#include <iostream>
#include <vector>
#include <map>
#include <cmath>
#define FOR(x, xs) for (__typeof( (xs).begin() ) x = (xs).begin(); x != (xs).end(); x++)

using namespace std;

static const int n = 2000000000;
static const int u = sqrt(n);
static const int v = sqrt(u);

map<int, int> factors (int n, vector<int> &primes) {
  map<int, int> fs;
  FOR (p, primes) {
    if ( *p * *p > n ) break;
    int c = 0;
    while ( n % *p == 0 ) {
      n /= *p;
      c ++;
    }
    if ( c > 0 ) fs[*p] = c;
  }
  if ( n > 1 ) fs[n] = 1;
  return fs;
}

int main() {

  // find priems up to u, using sieve
  vector<bool> isPrime( u + 1);
  FOR (q, isPrime) *q = true;
  isPrime[0] = false;
  isPrime[1] = false;

  for (int q = 2; q <= v; q++)
    if ( isPrime[q] )
      for (int s = q * q; s <= u; s += q)
	isPrime[s] = false;

  vector<int> smallPrimes;
  smallPrimes.push_back(2);
  for (int q = 3; q <= u; q += 2)
    if ( isPrime[q] ) smallPrimes.push_back(q);

  // find base-priems up to n, using sieve
  // p = 1, 25, 121 ( mod 168 )
  vector<bool> isBasePrime  ( n / 168 * 3 + 3 );
  FOR (p , isBasePrime) *p = true;
  isBasePrime[0] = false;

  FOR (p, smallPrimes)
    for (int q = *p; q < *p * 168; q += *p)
      if ( q % 168 == 1 || q % 168 == 25 || q % 168 == 121 )
	for (int s = *p == q ? *p * 169: q; s < n; s += *p * 168) switch ( s % 168 ) {
	  case 1  : isBasePrime[ s / 168 * 3 ] = false; break;
	  case 25 : isBasePrime[ s / 168 * 3 + 1 ] = false; break;
	  case 121: isBasePrime[ s / 168 * 3 + 2 ] = false; break;
	  }

  vector<int> basePrimes;
  for (int i = 0; i <= n / 168 *3 + 3; i++)
    if ( isBasePrime[i] ) switch ( i % 3 ) {
      case 0: basePrimes.push_back( i / 3 * 168 + 1 ); break;
      case 1: basePrimes.push_back( i / 3 * 168 + 25 ); break;
      case 2: basePrimes.push_back( i / 3 * 168 + 121 ); break;
      }

  // counting non-sqaure-type numbers
  int c = 0;
  FOR (p, basePrimes) c+=sqrt( n / *p ); // K^2 p

  for (int p = 0; basePrimes[p] <= u; p++) // K^2 p q
    for (int q = p + 1, m = n / basePrimes[p] ; basePrimes[q] <= m; q++)
      c += sqrt ( n / ( basePrimes[p] * basePrimes[q] ) );

  for (int p = 0, mp = cbrt(n); basePrimes[p] <= mp; p++) // K^2 p q r
    for (int q = p + 1, mq = sqrt ( n / basePrimes[p] ); basePrimes[q] <= mq; q++)
      for (int r = q + 1, mr = n / ( basePrimes[p] * basePrimes[q] ); basePrimes[r] <= mr; r++)
  	c += sqrt ( n / ( basePrimes[p] * basePrimes[q] * basePrimes[r] ) );

  // counting sqaure-type numbers
  for (int x = 1; x <= u; x++) {
    map<int, int> fs = factors(x, smallPrimes);
    bool b1 = false, b2 = false, b3 = false, b7 = false;
    FOR (pn, fs) {
      int p = (*pn).first;
      b1 = b1 || p % 4 == 1;
      b2 = b2 || p % 8 == 1 || p % 8 == 3;
      b3 = b3 || fs[2] >= 1 || p % 6 == 1;
      b7 = b7 || fs[2] >= 2 || p % 14 == 1 || p % 14 == 9 || p % 14 == 11;
      if ( b1 && b2 && b3 && b7 ) {
	c++;
	break;
      }
    }
  }
  cout << c << endl;
}