CodeIQ MAGAZINE

問題のおさらい

2つの自然数 a, b（a ≦ b）に対し、直角をはさむ2辺の長さが a と b である直角三角形を考えます。
この直角三角形の外接円を描きます。
さらに、その円に外接し、もとの直角三角形に相似となる直角三角形を描きます。
この直角三角形の直角をはさむ2辺の長さを a’, b’（a’ ≦ b’）とおきます。

例えば、(a, b) = (6, 8) のとき、(a’, b’) = (15, 20) となることが確かめられます。

整数 L に対し、a ≦ b ≦ L の範囲で、a’ と b’ がいずれも整数となるような整数の組 (a, b) の個数を F(L) と定義します。
例えば F(10) = 1、F(50) = 7、F(1000) = 173 となることが確かめられます。

標準入力から、自然数 L（1 ≦ L ≦ 10⁵）が与えられます。
標準出力に F(L) の値を出力するプログラムを書いてください。

方針

外接円と内接円に関する問題です。まずは a’ と b’ をそれぞれ a, b を用いて表すことを考えましょう。一見ややこしそうですが、三角形の外接円の性質と、内接円の性質とを順番に見ていけば大丈夫です。

始めに注目するのは、内側の三角形と外接円の関係です。下記の Wikipedia の記事にもあるように、直角三角形の外心（外接円の中心）は、三角形の斜辺の中点となります。つまり、斜辺の長さを c、円の半径を r とおくと、r = c / 2 の関係があります。

参考：外接円（Wikipedia）

次に注目するのは、外側の三角形と内接円の関係です。三角形の各点を A, B, C とおきます。三角形ABCの面積を二通りに表します。つまり下図において、(底辺)×(高さ)÷2 というよく知られた公式による表し方と、三角形OABと三角形OBCと三角形OCAの和とする表し方です。

c’ を斜辺の長さとおき、二つの表し方をイコールでつないで整理します。すると、a’, b’, c’ と r の関係が見えてきます。

最後の行では、三角形と外接円の関係のところで導いた r = c / 2 を代入しています。

次に、二つの三角形が相似であることを利用します。二つの相似比を k とおくと、各辺について a’ ＝ k a, b’ ＝ k b, c’ ＝ k c という関係が成り立ちます。これらを(※)式に代入します。

相似比 k が求められました。これを a’ ＝ k a と b’ ＝ k b に代入すれば、a’, b’ を a, b, c で表すことができます。

ピタゴラス数を列挙しよう

さて、a’ と b’ が整数となる条件を考えていきましょう。c は内側の直角三角形の斜辺の長さですから、有名なピタゴラスの定理により、c ＝ √(a²＋b²) と表されます。a’, b’ が整数となるためには、少なくとも式中にルート記号が残っていてはいけません。したがって、c が整数（有理数）であることが、a’, b’ が整数となるための必要条件であると分かります。

そのような a, b, c の組はピタゴラス数として知られています。ピタゴラス数を列挙するためには、単純には、全ての (a, b) の組み合わせに対して a²＋b² が平方数となるかどうかを確かめるという方法があります。しかし、本問の L ≦ 10⁵ という条件では、こうした総当たりの方法は非効率です。より高速にピタゴラス数を列挙する方法を用いましょう。

次の方法が有用です。ある互いに素な自然数 m, n が、m＞n であり、かつ m–n が奇数であれば、(m²–n², 2mn, m²+n²) がピタゴラス数となることを利用します。下記の Wikipedia の記事を参照して下さい。

参考：ピタゴラスの定理（Wikipedia）

こうして得られたピタゴラス数 a, b, c を a’, b’ の式に代入します。a’, b’ それぞれの分母が分子で割り切れるかどうかを確認し、どちらも割り切れれば、a’, b’ は整数となると言えます。なお実際には (a, b, c) ＝ (m²–n², 2mn, m²+n²) では、ピタゴラス数のうち最大公約数が 1 のもの（原始ピタゴラス数といいます）のみが得られます。(a, b, c) のそれぞれを 2 倍、3 倍、4 倍、…としていったものも忘れずチェックしましょう。

以上の処理をコードにしたものが以下です（C++）。

#include <iostream>
using namespace std;

// 最大公約数
long long GCD(long long x, long long y) {
    return x ? GCD(y % x, x) : y;
}

int main(int argc, char *argv[])
{
    long long L = 100;
    cin >> L;
    long long answer = 0, n, m, a, b, c, p;
    for (n = 1; 2 * n * n <= L; n++) {
        for (m = n + 1;; m++) {
            if ((n + m) % 2 == 0) continue;
            if (GCD(n, m) != 1) continue;
            a = m * m - n * n, b = 2 * m * n, c = m * m + n * n;
            for (p = 1;; p++) {
                // (a, b) をそれぞれ 1 倍, 2 倍, 3 倍, ... したものについて、
                // 2 辺の長さが L 以下、かつ a', b' の公式の分母が分子を割り切るかチェック。
                if (p * a > L || p * b > L) break;
                if ((p * c * (a + b + c)) % (2 * b) != 0) continue;
                if ((p * c * (a + b + c)) % (2 * a) != 0) continue;
                answer++;
            }
        }
    }
    cout << answer << endl;
}

黒魔術！？ コードをさらに高速化する不思議な行列

さて、本問の L = 10⁵ というリミットでは上記のコードで十分クリア可能です。が、さらなる工夫によって、もっと大きな L に対しても答えを出すことができます。ポイントは 2 つです。

一つは、原始ピタゴラス数の求め方です。上記の方法よりもずっと高速に原始ピタゴラス数を求める方法が存在します。次の 3 つの行列のそれぞれに、ピタゴラス数 (3, 4, 5) をかけると、新たな原始ピタゴラス数が得られるのです。

計算すると、それぞれ (5, 12, 13), (21, 20, 29), (15, 8, 17) となります。これらは全て原始ピタゴラス数です。さらにこれらのそれぞれに対し、同様に 3 つの行列をかけると、新しい原始ピタゴラス数が得られます。このように次々と行列をかけていくことで、原始ピタゴラス数をもれなくダブリなく生み出せるのです。

参考：ピタゴラスの数を生む行列（外部サイト）
参考：Tree of primitive Pythagorean triples（Wikipedia：英語）

もう一つの工夫は、原始ピタゴラス数を 2 倍、3 倍、4 倍、…していったものの扱い方です。上記のコードでは、a' と b' の公式の分母が分子を割り切るかどうかを、原始ピタゴラス数を 2 倍、3 倍、4 倍、…したもの一つ一つについてチェックしていましたが、もっと効率よくできます。原始ピタゴラス数を a' と b' の公式に代入した後で、それぞれ分母と分子を約分し、さらに二つの分母の最小公倍数を計算します。原始ピタゴラス数にこの値をかけたものの整数倍は確実に条件を満たしますので、一気にカウントを行えます。

コードを以下に示します。L = 10¹⁰ ぐらいであっても、十分短い実行時間で答えが得られます。

#include <iostream>
#define MAX(a,b)        ((a)>(b)?(a):(b))
using namespace std;

// 最大公約数
long long GCD(long long x, long long y) {
    return x ? GCD(y % x, x) : y;
}

// 最小公倍数
long long LCM(long long x, long long y) {
    return x * y / GCD(x, y);
}

// 幅優先探索
long long Next(long long a, long long b, long long c, long long L) {
    if (MAX(a, b) > L) return 0;
    long long ga = GCD(c * (a + b + c), 2 * b);
    long long gb = GCD(c * (a + b + c), 2 * a);
    // l: 約分後のa'とb'の分母の最小公倍数
    long long l = LCM(2 * b / ga, 2 * a / gb);
    long long r = L / MAX(a, b) / l;
    if (r == 0) return 0;
    r += Next( a - 2 * b + 2 * c,  2 * a - b + 2 * c,  2 * a - 2 * b + 3 * c, L);
    r += Next( a + 2 * b + 2 * c,  2 * a + b + 2 * c,  2 * a + 2 * b + 3 * c, L);
    r += Next(-a + 2 * b + 2 * c, -2 * a + b + 2 * c, -2 * a + 2 * b + 3 * c, L);
    return r;
}

int main(int argc, char *argv[])
{
    long long L = 100;
    cin >> L;
    long long answer = Next(3, 4, 5, L);
    cout << answer << endl;
}

みんなのコード

本問の挑戦者で、コードを公開して頂いた方のコードを Togetter でまとめました。（随時追加していきます。）こちらもぜひチェックして下さいね。

CodeIQ「ツイン・トライアングル」問題 みんなのコード

CodeIQ運営事務局より

Kawazoeさん、ありがとうございました！
現在、Kawazoeさんの最新問題が出題中です。
ぜひ挑戦してみてくださいね！