幾星霜の時を経て一問解く

問題 142 です。問題の内容は公式サイトをご覧ください。
公式サイトは英語ですが、有志の方々が日本語に翻訳しているページもあったりします。

最初に書いたプログラムはこちら。

(define nil '())

(define (pe142 low high)
  (define (iter-x x)
    (if (> x high)
        'none
        (let ((l (iter-y x (- x 1) nil)))
          (if (>= (length l) 2)
              (let* ((lst (sort (append l (cons x nil))))
                     (ret (pe142-b (sort (append l (cons x nil))))))
                (if (null? ret)
                    (iter-x (+ x 1))
                    lst))
              (iter-x (+ x 1))))))
  (define (iter-y x y lst)
    (if (= y 0)
        lst
        (iter-y x (- y 1)
              (if (and (integer? (sqrt (- x y))) (integer? (sqrt (+ x y))))
                  (cons y lst)
                  lst))))
  (iter-x low))

(define (pe142-b lst)
  (define (iter-y z l ans)
    (if (null? (cdr l))
        ans
        (let ((y (car l)))
          ;;(print "y=" y ", z=" z)
          (iter-y z (cdr l)
                  (if (and (integer? (sqrt (- y z))) (integer? (sqrt (+ y z))))
                      (append ans (list (cons y z)))
                      ans)))))
  (define (iter-z l ans)
    (if (null? (cddr l))
        ans
        (iter-z (cdr l) (append ans (iter-y (car l) (cdr l) nil)))))
  (iter-z lst nil))

なんとか解答は出力するのですが、とても遅いです。解を表示するまで 5-6 時間かかります。

次にいろいろ勉強して書いたプログラムがこちらです。

(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))))))))))))

後者のほうについてロジックを解説したいと思います。
ざっくりとした説明ですが、x と y と z を直接求めようとするのではなく、a^2, b^2, c^2 から x, y, z を求めます。

まず、

x + y = a^2 ... (a)
x - y = b^2 ... (b)

と定義します。

(a) と (b) を式のまま足し算します。
足し算するので、y は打ち消し合って無くなります。
両辺を 2 で割れば、式 (e) が得られます。

2x = a^2 + b^2
x = (a^2 + b^2) / 2 ... (e)

次に (a) と (b) を式のまま引き算します。
引き算するので x は打ち消し合って無くなります。
両辺を 2 で割れば、式 (f) が得られます。

2y = a^2 - b^2
y = (a^2 - b^2) / 2 ... (f)

ここまでで、x と y は a と b を使って求められることがわかりました。
次に、

x + z = c^2 ... (c)
x - z = d^2 ... (d)

と定義します。

先程と同様に、(c) と (d) を足し算します。
足し算するので、z は打ち消し合って無くなります。
両辺を 2 で割れば、式 (g) が得られます。

2x = c^2 + d^2
x = (c^2 + d^2) / 2 ... (g)

同様に、(c) と (d) を引き算します。
引き算するので、x は打ち消し合って無くなります。
両辺を 2 で割れば、式 (h) が得られます。

2z = c^2 - d^2
z = (c^2 - d^2) / 2 ... (h)

ここで、式 (e) と式 (g) を見てみると
左辺は x で同じなので、次の式を導くことができます。

a^2 + b^2 = c^2 + d^2 ... (i)

この式 (i) から、d^2 を a, b, c を使った式で表せることがわかります。

d^2 = a^2 + b^2 - c^2

d^2 を 式 (d) に代入して

x - z = a^2 + b^2 - c^2 ... (d')

ここで、式 (c) と式 (d') を引き算すると、x が打ち消されて

2z = c^2 - (a^2 + b^2 - c^2)
2z = 2c^2 - (a^2 + b^2)
z = c^2 - (a^2 + b^2) / 2 ... (j)

z を a, b, c を使って求める式を得ました。
また、式 (f) と式 (j) を使って y + z は、

y + z = (a^2 - b^2) / 2 + c^2 - (a^2 + b^2) / 2

このように書くことができ、a^2 は打ち消されて

y + z = c^2 - b^2

を得ます。
同様に y - z は、

y - z = (a^2 - b^2) / 2 - (c^2 - (a^2 + b^2) / 2)

このように書くことができ、b^2 は打ち消されて

y - z = a^2 - c^2

を得ます。

まとめると、x, y, z は a, b, c を使用して

x = (a * a + b * b) / 2
y = (a * a - b * b) / 2
z = c * c - (a * a + b * b) / 2

として求めることが出来ます。

また、a, b, c の大小関係を考えると、

a * a = x + y
b * b = x - y
c * c = x + z

と定義しており、問題文から

x > y > z > 0

なので

a > c > b

となります。

さて、次に a, b, c をどこからどこまで走査すれば良いかを考えます。

a > c > b

であることがわかったので、

a は、 3 以上 (ただし、無限に走査するプログラムは問題がある気がするので、一応 10000 以下、3 以上とします)
c は、a 未満 2 以上
b は、c 未満 1 以上

となります。

外側から a, c, b の三重ループとなります。a, b, c の順ではない事に注意してください。
さらに以下の式を使い、ループ処理をスキップすることでプログラムを高速化します。

y + z = c^2 - b^2
y - z = a^2 - c^2

a ループと c ループにおいては、a^2 - c^2 を使用し、
c ループと b ループにおいては、c^2 - b^2 を使用します。
この式が平方数であるかを判定し、平方数でないときはループ処理をスキップします。

蛇足ですが、Scheme ではループは再帰で書くのですが、
単純なループ処理を再帰で書くとプログラムが非常に汚くなるので、
今回は do を使用して書きます。

次に解の除外のしかたについてです。
まず、Scheme は静的型付けをしない言語なので、x が整数でない場合がありえますので除外します。
つぎに x と y は常にプラスの値を取りますが、z だけはマイナスの値を取る場合があります。
問題文から x > y > z > 0 なのでマイナス値は除外します。

三重ループの内側で、a, b, c から x, y, z を求めて、
問題文にある以下の式がすべて平方数になる事をチェックします。

x + y
x - y
x + z
x - z
y + z
y - z

ただし、ループの高速化判定に y + z と y - z が平方数であることを利用しているため、
最終的な判定するロジック部分では、y + z と y - z が平方数になることは明らかです。
ですので、以下の 4 つの式が平方数になることをチェックすれば十分です。

x + y
x - y
x + z
x - z

すべての条件を満す最初に見つかった x, y, z を足し算したものが解答となります。

gosh> (time (pe142 10000))
;(time (pe142 10000))
; real   0.142
; user   0.125
; sys    0.000
1006193

ロジックの解説では、数式の変形や考え方について端折らず丁寧に説明したつもりです*1。そんなに難しい内容ではなく中学校の数学レベルで理解できると思います。あんまり数学が得意じゃないんだけど・・・という人でも数学の問題をプログラムで解くって面白そう!と興味を持ってもらえれば何よりです。

projecteuler.net

*1:数学が得意な人はすごく端折るので、ついていけない事がよくありますw