再びピタゴラス数を求める問題

C Project Euler

同じテーマにばかり執着していて, 木を見て森を見ずの感がなきにしもあらずですが, 新年最初の日記でも, 再び同じ問題を考えてみたいと思います.

Project Euler の Problem 9 (Jan 25, 2002) です. しつこいようですが, 再掲します. 出題内容は次のとおりです.

A Pythagorean triplet is a set of three natural numbers, a < b < c, for which,

a2 + b2 = c2

For example, 32 + 42 = 9 + 16 = 25 = 52.

There exists exactly one Pythagorean triplet for which a + b + c = 1000.

Find the product abc.

準備

この問題の解答が 31875000 (= 200 × 375 × 425) の 1 通りしかないことは, すでに数日前の 日記 (Dec 28, 2011) で解いたばかりでしたので, わかっています. が, a + b + c が 1000 でなければ, 解答は 0 通りの場合もあったり, 2 通り以上の場合もあったりします. ですので, 問題を解く関数を書くにあたっては, 複数のピタゴラス数の組を返すことがありうる点を想定したいと思います.

そこで, まずは, 次のような構造体を定義しておきます.

typedef struct _pythagorean {
    int a;
    int b;
    int c;
    struct _pythagorean *next;
} pythagorean;

すなわち, 条件に合うピタゴラス数が 2 組以上ある場合は, それらを連結リストにすることにします.

また, 条件に合うピタゴラス数の組が見つかるたびに連結リストの末尾に追加していくための関数も用意しておくことにします.

void _add(pythagorean **p0, int a0, int b0, int c0) {
    pythagorean *t, *p;
    t = (pythagorean *)malloc(sizeof(pythagorean));
    t->a = a0;
    t->b = b0;
    t->c = c0;
    t->next = NULL;
    if (*p0 == NULL)
        *p0 = t;
    else {
        pythagorean *p;
        p = *p0;
        while (p->next != NULL)
            p = p->next;
        p->next = t;
    }
    return;
}

ポインタや, ポインタへのポインタの扱い方にまだ慣れていないので, もしかしたらちょっとダサい書き方かもしれませんが, それはまた今後少しずつ改善していくことができたらいいなと思います.

それから, 今回は, ある数とある数が互いに素であるかどうかを判定する必要のある解き方も試してみますので, 2 個の数の最大公約数を返す関数も用意しておきます.

int _gcd(int a0, int b0) {
    int a = a0, b = b0, tmp;
    do {
        tmp = a % b;
        a = b;
        b = tmp;
    } while (tmp != 0);
    return a;
}

ただし, 引数として 0 や負の数が与えられた場合については特に何も想定していませんので, この関数は本当にただ今回の問題のためだけにしか使えないでしょう.

オーソドックスな解法

オーソドックスな解法と呼んでいいのかどうかはよくわかりませんが, まずは a および b を for ループで回して, ひたすら条件に合うものを探し続ける解法について考えてみます. これは, 以前に解いたときにも使用した方法です. ただし, コードはだいぶ書き直しました.

pythagorean *search0(int sum) {
    if (sum % 2 != 0)
        return NULL;
    pythagorean *p;
    p = NULL;
    int a, b, c, b0;
    for (a = 3; a <= (sum - 3) / 3; a++) {
        b0 = a + 1;
        if (a + b0 < sum / 2 + 1)
            b0 = sum / 2 + 1 - a;
        for (b = b0; b < (c = sum - a - b); b++)
            if (a * a + b * b == c * c)
                _add(&p, a, b, c);
    }
    return p;
}

a + b + c が奇数となることはありません. というのも, a + b + c が奇数であるためには, 3 数すべてが奇数である, もしくは 3 数のうち 1 数のみが奇数であることが必要ですが, 以下のように, それはありえないからです.

  • a, b, c がいずれも奇数であるとすると, a2 + b2 = c2 は奇数と奇数を足すと奇数になることを意味する. しかし, これはありえない.
  • a が奇数であり, b と c が偶数であるとすると, a2 + b2 = c2 は奇数と偶数を足すと偶数になることを意味する. しかし, これはありえない.
  • a と c が偶数であり, b が奇数であるとすると, a2 + b2 = c2 は偶数と奇数を足すと偶数になることを意味する. しかし, これはありえない.
  • a と b が偶数であり, c が奇数であるとすると, a2 + b2 = c2 は偶数と偶数を足すと奇数になることを意味する. しかし, これはありえない.

以上から, search0 関数では, 与えられた引数が奇数の場合にはただちに NULL を返して処理を終了することにしています.

a をループで回す際には 3 からスタートしています. 0 < a < b < c であり, a + b > c であることに鑑みますと, a が 0 や 1 になることはありえません. a が 2 になることはありえますが, a + b + c が偶数であることを考慮に加えますと, これもありえなくなります. したがって 3 からスタートしているわけです.

a の最大値については, 次のように考えました. a は最大でも b - 1 であり, 最大でも c - 2 です. ということは, a は, 与えられた引数から 3 を引いたものを 3 で割った数よりも, 大きくはなりません.

もうひとつ, b のループについては, 次のように考えます. 基本は a + 1 からスタートします. ただ, a があまりに小さいと, b のループは a + b が c を越えるまでは何もせずに無駄に繰り返すばかりとなって非効率です. そこで, 常に a + b > c が成り立つように (つまり与えられた引数を 2 で割った数よりも大きくなるように) スタート地点を調整しています.

ちなみに, この search0 関数では, 引数として 1000 が与えられると, 内側の for ループが 20584 回繰り返されます.

別の解法

つづいて, 別の解法をやってみます. すなわち, 昨年の最後の 日記 (Dec 31, 2011) に書いた方法です.

pythagorean *search1(int sum) {
    if (sum % 2 != 0)
        return NULL;
    int s = sum / 2;
    pythagorean *p;
    p = NULL;
    int g, m, n, m0, m1;
    int a, b, c, tmp;
    for (g = 1; g <= s / 6; g++) {
        if (s % g != 0)
            continue;
        m0 = 2;
        m1 = s / g;
        for (m = m0; m * (m + 1) <= m1; m++) {
            if (m1 % m != 0)
                continue;
            n = m1 / m - m;
            if (m <= n || 1 < _gcd(m, n) || (m + n) % 2 == 0)
                continue;
            a = g * (m * m - n * n);
            b = 2 * g * m * n;
            c = g * (m * m + n * n);
            if (a > b) {
                tmp = a;
                a = b;
                b = tmp;
            }
            _add(&p, a, b, c);
        }
    }
    return p;
}

与えられた引数を N としますと, 2 g m (m + n) = N を満たす g, m, n を探すことが主眼となります.

その際, 2 s = N とすれば, g も m も s を割り切れなければなりません. また, n < m であり, m と n は互いに素であり, m と n の偶奇は不一致であることも必要です. それらの条件を判定するため, if 文がやや多くなってしまいました.

g のループは, g が s を 6 で割った数よりも大きくなった時点で終了します. なぜ 6 かというと, n が最小でも 1 で m が最小でも 2 だと考えれば, m (m + n) が最小でも 6 だからです.

m のループは, それが終了する条件はいいとして, スタート地点については改善の余地があります. というのも, 0 < n < m に鑑みれば,

m2 < m (m + n) < 2 m2

が成り立つはずだからです. ただ,

        m0 = sqrt(s / (2 * g)) + 1;

といったように平方根を使うのは, 整数の問題に浮動小数点数の演算を紛れ込ませることによって, 何らかの落とし穴ができるのは嫌だなあと思ってしまいました. そこで, m は 2 からスタートさせることにしました. ((わたしが Project Euler で標準ライブラリの math.h をあまり使わないのは, 本当に単に浮動小数点数の扱いが面倒くさそうという先入観があるからなのですね..))

あと, ささいなところですが, g と m と n を探し出すだけでは a < b が満たされるとはかぎりません. 逆になることがありえます. ですので, その場合は a と b を入れ替えることが必要になります.

ちなみに, この search1 関数では, 引数として 1000 が与えられると, 内側の for ループが 63 回繰り返されます.

両者の比較

ループの繰り返し回数で見ても, これはだいぶ差がつきそうだというのはわかります. 次のようなコードで比較してみました.

int main(void) {
    clock_t t0, t1;
    const int N = 1000;
    pythagorean *(*pf[])(int) = { search0, search1 };
    pythagorean *p;
    int i, j;
    for (i = 0; i < sizeof(pf) / sizeof(pf[0]); i++) {
        printf("[%d]\n", i);
        t0 = clock();
        for (j = 0; j < 10000; j++)
            p = pf[i](N);
        t1 = clock();
        printf("  Elapsed time ... %u\n", t1 - t0);
        if (p == NULL)
            printf("  Not found.\n", i);
        else {
            printf("  Found.\n", i);
            do {
                printf("    %d, %d, %d\n", p->a, p->b, p->c);
                p = p->next;
            } while (p != NULL);
        }
    }
    return EXIT_SUCCESS;
}

それぞれの解法を 10000 回ずつ繰り返してみましたところ, わたしの手元の環境では, search0 が 1300 - 1350 ミリ秒, search1 が 15 - 20 ミリ秒, といった感じでした.

わたしが思うに, おそらく計算量オーダーも違うんじゃないでしょうか. search0 では a のループも b のループもそれぞれ O(N) で, したがって全体としてはおそらく O(N2) でしょう. それに対して, search1 では g のループは O(N) ですが, m のループは O(N1/2) で, したがって全体としてはおそらく O(N3/2) なのではないでしょうか. ちょっと自信はありませんが.

やり残したこと

この問題に関しましては, 少し興味を覚えたまま足踏みをしているテーマがふたつあります.

正解の数

これは大晦日の日記にも書いたことですが, どういう場合に何組の正解があるものでしょうか. a, b, c の合計が 1000 の場合, 正解は (200, 375, 425) の 1 組でした. もし 1200 だとしたら, 大晦日の日記では少なくとも 2 組はありそうだと書きましたが, 実際には 5 組ありました. (48, 575, 577), (75, 560, 565), (200, 480, 520), (240, 450, 510), (300, 400, 500) です.

実際に 200 以下で探してみますと, 次のようになりました.

Sum of a, b, cPattern(s) of a, b, c
12(3, 4, 5)
24(6, 8, 10)
30(5, 12, 13)
36(9, 12, 15)
40(8, 15, 17)
48(12, 16, 20)
56(7, 24, 25)
60(10, 24, 26), (15, 20, 25)
70(20, 21, 29)
72(18, 24, 30)
80(16, 30, 34)
84(12, 35, 37), (21, 28, 35)
90(9, 40, 41), (15, 36, 39)
96(24, 32, 40)
108(27, 36, 45)
112(14, 48, 50)
120(20, 48, 52), (24, 45, 51), (30, 40, 50)
126(28, 45, 53)
132(11, 60, 61), (33, 44, 55)
140(40, 42, 58)
144(16, 63, 65), (36, 48, 60)
150(25, 60, 65)
154(33, 56, 65)
156(39, 52, 65)
160(32, 60, 68)
168(21, 72, 75), (24, 70, 74), (42, 56, 70)
176(48, 55, 73)
180(18, 80, 82), (30, 72, 78), (45, 60, 75)
182(13, 84, 85)
192(48, 64, 80)
198(36, 77, 85)
200(40, 75, 85)

ここには何かの法則があるのでしょうか. まあ, ないわけはないと思うのですが, それが何なのかまだ見えてきません.

abc は 60 の倍数

もうひとつ気になっているテーマは, 積 abc は必ず 60 の倍数になるというものです たまたまインターネット上で何かを検索していて, ピタゴラス数のある性質 という記事で見かけたものです. 今回の問題は abc を答えさせるものでしたので, abc が 60 の倍数であることを使って正解を探し求めるという戦略もありえるでしょうか. もっとも, おもしろそうではあるものの, その性質をどのように活かせばよいのか, わたしにはまださっぱり見えてきません.

まとめ

いったい, いつまで, この問題にこだわっているんだか..