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 = 0;
  for (int a = 2; a <= N / 3; a++) {
    if (a % 10000 == 0 ) cout << "now count " << a << endl;
    int k = (a + 1) % 2 == 0 ? 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 == 0 ) 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 0.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

まぁ，速いけど，そんなに速くない．
解を全列挙しているし，解の個数も結構多いから，仕方がないのか？
もう少し速いことを期待していたから，残念．