Problem 223

223 Almost right-angled triangles I

Problem 223 – Project Euler

なる整数の組を探す．

Problem 221と近いものを感じますが…

まず，思いつくのは，

と因数分解する方法．Haskellでは遅いので，C++で

#include <iostream>
#include <algorithm>
#include <vector>
#include <cmath>
#define N 25000000
#define A (N / 3 + 1)
#define FOR(x, xs) for (__typeof( (xs).begin() ) x = (xs).begin(); x != (xs).end(); x++)
#define ALL(xs) (xs).begin(), (xs).end()
using namespace std;
vector<long long> concat (vector< vector<long long> >& xxs) {
vector<long long> ys;
FOR (xs, xxs)
ys.insert(ys.end(), ALL(*xs));
return ys;
}
vector<int> joint (vector<int>& xs, vector<int>& ys) {
vector<int> zs;
zs.insert(zs.end(), ALL(ys));
zs.insert(zs.end(), ALL(xs));
sort(ALL(zs));
return zs;
}
vector<long long> divisors (vector<int> ps) {
vector< vector<long long> > xxs(1, vector<long long> (1,1));
int p = 1;
vector<long long> ys;
FOR (q, ps) {
if ( p != *q ) {
ys = concat(xxs);
xxs.clear();
xxs.push_back(ys);
}
ys = xxs.back();
FOR (y, ys) (*y) *= (*q);
xxs.push_back(ys);
p = *q;
}
ys = concat(xxs);
sort(ALL(ys));
return ys;
}
int main() {
cout << "N is " << N << endl;
cout << "sieve start." << endl;
vector< vector<int> > factors( A + 1 );
for (int i = 2; i <= A; i++) {
if ( !factors[i].empty() ) continue;
long long j = i; // j = i^n
while ( j <= A ) {
for (int k = j; k <= A; k += j)
factors[k].push_back(i);
j *= (long long) i;
}
}
cout << "sieve complete." << endl;
cout << "counting start." << endl;
long long c = ;
for (int a = 2; a <= N / 3; a++) {
if (a % 10000 ==  ) cout << "now count " << a << endl;
int k = (a + 1) % 2 ==  ? 2: 1;
vector<int> f1 = factors[ (a + 1) / k ], f2 = factors[ (a - 1) / k ];
vector<long long> ds = divisors( joint(f1, f2) );
if ( (a + 1) % 2 ==  ) FOR (d, ds) *d *= k;
int upp = N - a;
int low = ceil( sqrt(2LL*a*a - 1) ) + a;
FOR (d, ds) {
if ( *d < low ) continue;
if ( *d > upp ) break;
c++;
}
}
c += (N - 1) / 2;
cout << "result " << c << endl;
}

これでも，遅い．

でフォーラムでこんなページの存在を知った．

http://www.cut-the-knot.org/pythagoras/PT_matrix.shtml

でHaskellでコード書いてみた．

import Data.List (unfoldr, foldl')
trans :: (Num a) => [a] -> [a]
trans [a,b,c] = [ -a - 2*b + 2*c, -2*a - b + 2*c, -2*a - 2*b + 3*c]
next :: (Num a) => [a] -> [[a]]
next x@[a,b,c] | a == b    = map trans [ [-a, b, c], [-a, -b, c]]
| otherwise = map trans [ [-a, b, c], [a, -b, c], [-a, -b, c]]
p223 :: Integer -> Integer
p223 l = foldl' count .unfoldr step $ [ [1,1,1], [1,2,2] ]
where step [] = Nothing
step xs = Just (xs, filter feasible.concatMap next $ xs)
feasible xs = sum xs <= l
count c xs = c + toInteger (length xs)
main :: IO ()
main = print.p223 $ 25*10^6

まぁ，速いけど，そんなに速くない．

解を全列挙しているし，解の個数も結構多いから，仕方がないのか？

もう少し速いことを期待していたから，残念．