Problem 229

229 Four Representations using Squares

この問題は難しかった．

なんてったって，問題サイズが 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 = ;
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 = ;
while ( n % *p ==  ) {
n /= *p;
c ++;
}
if ( c >  ) 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[] = 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[] = 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 = ; i <= n / 168 *3 + 3; i++)
if ( isBasePrime[i] ) switch ( i % 3 ) {
case : 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 = ;
FOR (p, basePrimes) c+=sqrt( n / *p ); // K^2 p
for (int p = ; 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 = , 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;
}

コンテンツ