Problem 39 - ピタゴラス数

Scheme

オイラーといえば、フェルマーの最終定理に挑んだ数学者のひとり。

フェルマーの最終定理といえば、a^n + b^n = c^n。

というわけで、Problem 39 - PukiWiki

辺の長さが{a,b,c}と整数の3つ組である直角三角形を考え, その周囲の長さをpとする. p = 120のときには3つの解が存在する:

{20,48,52}, {24,45,51}, {30,40,50}

p < 1000 で解の数が最大になる p を求めよ.

何処の数を軸に考えるかというのが問題になりそう。


pを軸に考えると爆発的に計算量が増えてしまうので、a,bを軸にして、a <= bで回してみる。


すると、pをキャッシュしておく必要が出てくるので破壊的にベクタを使ってみる。

(define (problem39 n)
  (letrec ((square (lambda (n) (* n n)))
           (data (make-vector n 0))
           (inc! (lambda (n) (vector-set! data n (+ 1 (vector-ref data n)))))
           (iter (lambda (a b)
                   (let* ((c (sqrt (+ (square a) (square b))))
                          (p (+ a b c)))
                     (cond ((>= a n) data)
                           ((>= p n) (iter (+ a 1) (+ a 1)))
                           (else
                             (if (exact? c)
                                 (inc! p))
                             (iter a (+ b 1))))))))
    (vector-fold
      (lambda (i acc v)
        (if (> v (cdr acc))
            (cons i v)
            acc))
      '(0 . 0)
      (iter 1 1))))

(problem39 121) ; (120 . 3)
(problem39 1000) ; (840 . 8)

いっこづつ増やすのは無駄そうな気がするんだけど、√(a^2 + b^2)が整数になる条件がイマイチわからない。

追記

ピタゴラス数というらしい。

証明もじっくり取り組みたいのでゆっくりやる。

・・・こうやって見てると、フェルマーの最終定理なんて簡単に解けそうな気がしてくる ← 数学ド素人

ピタゴラス数の性質メモ

ピタゴラス数の生成式

a = m^2 - n^2, b = 2mn, c = m^2 + n^2(m>nでm,nは互いに素)を使うと、a,b,cそれぞれが互いに素なピタゴラス数を生成できるが、

今回の場合互いに素ではないピタゴラス数も求めなければならないので、倍数を取る必要がある。

しかぁし・・・重複が出るので、重複を取り除く手段が必要になる・・・・。

できた

やってることはたいしたことないんだけど、alistのユーティリティをいくつか追加したので案外長め。

(define (collect l)
  (letrec ((inc-alist (lambda (x a)
                        (cond ((null? a) (list (cons x 1)))
                              ((eq? (caar a) x) (cons
                                                 (cons x (+ 1 (cdar a)))
                                                 (cdr a)))
                              (else
                                (cons (car a)
                                      (inc-alist x (cdr a))))))))
    (fold-right inc-alist '() l)))

(collect '(1 1 1 2 2 3)) ; ((3 . 1) (2 . 2) (1 . 3))

(define (max-alist l)
  (fold (lambda (x m)
          (if (or (not m) (> (cdr x) (cdr m)))
              x
              m))
        #f
        l))

(max-alist (collect '(1 1 1 2 2 3))) ; (1 . 3)

(define (problem39 x)
  (letrec ((square (lambda (n) (* n n)))
           (calc (lambda (sum l n)
                   (if (>= (* n sum) x)
                       '()
                       (cons (map (lambda (x) (* n x)) l)
                             (calc sum l (+ n 1))))))
           (iter (lambda (m n acc)
                   (let* ((a (- (square m) (square n)))
                          (b (* 2 m n))
                          (c (+ (square m) (square n)))
                          (sum (+ a b c)))
                     (cond ((>= sum x) (if (= n m)
                                           acc
                                           (iter (+ n 1) (+ n 1) acc)))
                           ((not (= (gcd m n) 1)) (iter (+ m 1) n acc))
                           (else
                             (iter (+ m 1) n (lset-union equal?
                                                         acc
                                                         (calc sum
                                                               (list (min a b) (max a b) c)
                                                               1)))))))))
    (max-alist
      (sort (collect (map (lambda (l) (apply + l))
                          (iter 2 1 '())))
            (lambda (a b)
              (< (car a) (car b)))))))

(problem39 121) ; (120 . 3)
(problem39 1000) ; (840 . 8)

メインのループは16^2以下しか回らない。かなり早くなった。後は集計を早くすると・・・。

参考