Problem 143 - もうカツ丼はいいよな

Problem 143 - Project Euler

ABCをどの内角も120°未満である三角形とする。Xを三角形の内点とし、XA = p, XB = q, XC = rとする。
フェルマー(ピエール・ド・フェルマー - Wikipedia)はトリチェリ(エヴァンジェリスタ・トリチェリ - Wikipedia)に対し、p + q + rを最小化するようなXを見つけられるかという問を出した。
トリチェリは正三角形AOB, BNC, AMCをABCの各辺の外側に作図すると、それらの外接円はABCの内側の一点で交わることを示し、さらにその交点Tがp + q + rを最小化することを示した。点Tはトリチェリ・フェルマー点と呼ばれる。加えて興味深いことに、p + q + rが最小化されているとき、AN = BM = CO = p + q + rであり、AN, BM, COは点Tで交わる。

p + q + rが最小化されており、a, b, c, p, q, rの全てが正整数であるとき、三角形ABCをトリチェリ三角形と呼ぼう。例えばa = 399, b=455, c=511はトリチェリ三角形であり、p + q + r = 784である。
トリチェリ三角形について、全ての相異なる値p + q + r ≦ 120000の総和を求めよ。

解説

Tを中心とした3つの角は全て120°になる(らしい)。
よって、余弦定理から
$a^2 = q^2 + r^2 + qr$
$b^2 = p^2 + r^2 + pr$
$c^2 = p^2 + q^2 + pq$
ここから、
$x^2 + xy + y^2$
が平方数となるようなとき、(x, y)は(p, q, r)のうち2つであり得る、ということがわかる。
方針としては、

(x, y)のペアをx + y < 120000の範囲で全て集める
(x, y)をエッジとみなしてグラフを作成
3つでループになっているノードの値の和がp + q + r

という感じになる。
ところで $x^2 + xy + y^2$ が平方数というのはピタゴラス数に似ているが、ピタゴラス数を生成するのと同様の考え方で(x, y)を生成することができる(らしい)。
詳細はFermat points and parameterizing the 120 degree integer triangle (Project Euler 143) – Lucky's Notesに詳しく解説されたものがある。
ここから先はこの解説を読んで解いたので解説はそちらに譲るとして、要するに数のペアは互いに素である正整数m > nについて

$2mn + n^2$
$m^2 - n^2$

として作成できる。また生成したペアの整数倍も条件を満たす。
ただし、m - nが3で割り切れるときに生成される数のペアは3で割り切ることができ、ペアを3で割ってできるペアは別の(m, n)を用いて生成可能。よってm - nが3で割り切れる場合はスキップしても問題ない。

実装

問題は実装。forumとか除く限り、ここまで問題を整理しておけば大抵の言語では素直に実装してしまえば数秒以内に答えが出る。
ただRだとなかなかそうはいかない。最初Fermat points and parameterizing the 120 degree integer triangle (Project Euler 143) – Lucky's NotesにあるC++のコード(2秒で答えが出るらしい)をそのままRに書き換えたら4時間くらいかかった。話にならない。
問題は(多分)Rの添字アクセスが遅いという点にある。R in a Nutshellには次のようなことが書いてある。

looking up a value in a vector with $n$ elements or a list can take $O(n)$ time.

なのでベクトルの要素を添字でなぞるだけで $O(n^2)$ かかってしまう(多分)。
ところでCRANにhashというパッケージがある(CRAN - Package hash)。
これを使えばなんとか耐えられる時間内に終了するようになるものの、それでもまだ7分前後かかる。
グラフということでigraph使ったりもしてみたけどやっぱり話にならなかった。
ペアを生成するところまでは高速に出来てるので、あとはグラフを効率よく扱えれば1分以内に解けるとは思うんだけど…

library(hash)
library(gmp)

limit <- 120000
pair <- matrix(numeric(1e6), ncol=2)

## 10sec
cnt <- 1
for(n in 1:(sqrt(limit)+1)){
  for(m in 1:max(1, n-1)){
    q <- 2*m*n + m^2
    r <- n^2 - m^2
    if(q + r > limit) break
    if(n-m %% 3 == 0) next
    if(max(gcd(m, n)) != 1) next
    l <- 1:((limit/(q+r)) + 1)
    pair[cnt:(cnt+max(l)-1),] <- matrix(c(q*l, r*l), ncol=2)
    cnt <- cnt + max(l)
  }
}
pair <- pair[pair[,1]!=0,]
pair <- rbind(pair, pair[,2:1])

## 80sec
h <- hash()
for(i in 1:length(pair[,1])){
  h[[as.character(pair[i,1])]] <-
    c(h[[as.character(pair[i,1])]], pair[i,2])
}

## 5min
sums <- logical(limit)

for(i in 1:length(pair[,1])){
  a <- pair[i,1]
  b <- pair[i,2]
  for(c in intersect(as.vector(values(h,a)),
                     as.vector(values(h,b)) )){
    if(a + b + c < limit) sums[a + b + c] <- TRUE
  }
}

sum((1:length(sums))[sums])