はじめに
今回は Project Euler の問題 142 を Scheme で解いてみます。
問題 142 「完全平方数コレクション」
x + y, x - y, x + z, x - z, y + z, y - z が全て平方数となる整数の組 x > y > z > 0 で、最小の x + y + z を求めよ。
考え方
闇雲に x, y, z を走査して総当りで解を求めようとしても膨大な時間がかかってしまいます。この問題を解くにはいかに短時間で処理するかが肝となります。
問題として与えられた計算式が、すべて平方数 a, b, c で表せることを示し、平方数 a, b, c から x, y, z を求めます。これにより注目すべき変数の数を大幅に減らします。また、得られた計算式から各変数の走査範囲を適切に絞り込み、さらに余分な計算を省くようにします。
x + y = a^2 ... <式1>
x - y = b^2 ... <式2>
とおきます。
<式1> と <式2> を足すと、y の項が相殺されて
2x = a^2 + b^2
x = \frac{a^2 + b^2}{2} ... <式3>
を得ます。
同様に、<式1> から <式2> を引くと、x の項が相殺されて
2y = a^2 - b^2
y = \frac{a^2 - b^2}{2} ... <式4>
を得ます。
次に、
x + z = c^2 ...<式5>
x - z = d^2 ...<式6>
とおきます。
同様に <式5> と <式6> を加減算することにより
x = \frac{c^2 + d^2}{2} ...<式7>
z = \frac{c^2 - d^2}{2} ...<式8>
を得ます。
ここで、<式3> と <式7> は左辺が同じ x ですから
a^2 + b^2 = c^2 + d^2
が成り立ち、
d^2 = a^2 + b^2 - c^2
を得ます。さらに <式6> から、
x - z = d^2
としたので、
x - z = a^2 + b^2 - c^2 ... <式9>
つまり、変数 d は使わずとも x - z は計算できます。
<式5> から <式9> を引くと、x の項が相殺されて
2z = c^2 - (a^2 + b^2 - c^2)
2z = 2c^2 - (a^2 + b^2)
z = c^2 - \frac{a^2 + b^2}{2} ... <式10>
を得ます。
<式4> と <式10> から y + z は、
y + z = \frac{a^2 - b^2}{2} + c^2 - \frac{a^2 + b^2}{2}
と表すことができ、この式を整理すると、
y + z = c^2 - b^2 ... <式11>
を得ます。
同様に y - z は、
y - z = \frac{a^2 - b^2}{2} - (c^2 - \frac{a^2 + b^2}{2})
と表すことができ、この式を整理すると、
y - z = a^2 - c^2 ... <式12>
を得ます。
以上をまとめると x, y, z はそれぞれ
x = \frac{a^2 + b^2}{2}
y = \frac{a^2 - b^2}{2}
z = c^2 - \frac{a^2 - b^2}{2}
と表せます。また問題文の
x + y, x - y, x + z, x - z, y + z, y - z
はすべて a, b, c の 3 つの変数を使って表すことが出来るということもわかりました。
各変数の大小関係
変数 a、b、c の大小関係は、
a^2 = x + y
b^2 = x - y
c^2 = x + z
ですから、問題文より
x > y > z > 0
ですので、
a > c > b
となることがわかります。
走査範囲
変数 a, c, b の走査範囲についてです。
問題文より各平方数は整数でなければなりません。
変数 a の範囲は
10000 >= a >= 3
(無限ループはまずいので上限を 10000 と仮に設定します。)
変数 c の範囲は、
a > c >= 2
変数 b の範囲は、
c > b >= 1
とすれば良いことがわかります。
コード実装に向けて
さらにコードの実装するにあたって検討しておくべき事を整理しておきます。
プログラムコードは a, c, b の 3 つの変数を使用した 3 重ループとなります。
(a, b, c のアルファベット順ではない事に注意してください。)
次に以下の式を使って、プログラムを高速化します。
y + z = c^2 - b^2 ... <式A>
y - z = a^2 - c^2 ... <式B>
a 変数のループではスキップ条件はありません。c 変数のループ内では、<式B> を、b 変数のループ内では、<式A> を使用して平方数になるかをチェックします。平方数にならない場合は処理を直ちにスキップして無駄な計算を省きます。
また、x と y は常にプラスの値を取りますが、式10 から z はマイナスの値を取る可能性があります。問題文より与えられる式は
x > y > z > 0
ですので、マイナス値は除外します。最も内側となる b 変数のループにおいて z がプラスであることをチェックすれば良いです。
プログラムコードは a, b, c から x, y, z を求め、問題文にある以下の式がすべて平方数になる事をチェックします。
x + y, x - y, x + z, x - z, y + z, y - z
ただし、y + z, y - z の 2 つはループのスキップ条件で使用するため、チェック時は必ず平方数になりますので、除外することが出来ます。つまり、以下の 4 つ
x + y, x - y, x + z, x - z
が平方数であるかをチェックするだけで十分です。
加えて問題文より x, y, z は整数でなければならないので x が整数でないものは解にはなれません。これも除外します。
すべての条件を満たした x, y, z が見つかったら x + y + z を計算して処理を中断します。これが本問題の解となります。
プログラムコード
Scheme では以下のように書くことができます。
(define (square n)
(* n n))
(define (square? n)
(integer? (sqrt n)))
(define (pe142 limit)
(call/cc
(lambda (cc)
(do ((a 3 (+ a 1)))
((> a limit))
(do ((c 2 (+ c 1)))
((>= c a))
(and (square? (- (square a) (square c)))
(do ((b 1 (+ b 1)))
((>= b c))
(and (square? (- (square c) (square b)))
(let ((x (/ (+ (square a) (square b)) 2))
(y (/ (- (square a) (square b)) 2))
(z (- (square c) (/ (+ (square a) (square b)) 2))))
(and (integer? x)
(> z 0)
(square? (+ x y))
(square? (- x y))
(square? (+ x z))
(square? (- x z))
(cc (+ x y z))))))))))))
実行
gosh> (time (pe142 10000))
;(time (pe142 10000))
; real 0.152
; user 0.141
; sys 0.016
1006193
解は、1,006,193 です。
おわりに
処理時間は、私の PC では 152 ミリ秒でした。
そんなに新しい PC ではないので、みなさんの PC で実行するほうが短時間になると思います。