Project Euler - Problem 510

Tangent Circles

Problem 510 (日本語訳)

円 A と B がお互いに, そして線分 L と異なる3点で接している.
円 C が A, B, L の内部空間にあり, それぞれ3つすべてに接している.
A, B, C の半径をそれぞれ r_A, r_B, r_C としよう.



0 < r_A ≤ r_B ≤ n に対し S(n) = Σ r_A + r_B + r_C としよう, ここで r_A, r_B, r_C は整数とする.
0 < r_A ≤ r_B ≤ 5 の時の唯一の解は r_A = 4, r_B = 4, そして r_C = 1, したがって S(5) = 4 + 4 + 1 = 9 となる. また S(100) = 3072 が与えられている.

S(10⁹) を求めよ.

(日本語訳より)

$ H_A H_B = P O_C = \sqrt{ {O_A O_C}^{2} - {O_A P}^{2} } = \sqrt{ {(r_A + r_C)}^{2} - {(r_A - r_C)}^{2} } = 2 \sqrt{ r_C r_A }$
$ H_C H_B = Q O_C = \sqrt{ {O_B O_C}^{2} - {O_B Q}^{2} } = \sqrt{ {(r_B + r_C)}^{2} - {(r_B - r_C)}^{2} } = 2 \sqrt{ r_B r_C }$
$ H_C H_A = R O_A = \sqrt{ {O_A O_B}^{2} - {O_B R}^{2} } = \sqrt{ {(r_B + r_A)}^{2} - {(r_B - r_A)}^{2} } = 2 \sqrt{ r_A r_B }$

$ H_A H_C = H_A H_B + H_C H_B $ なので，

$ \sqrt{ r_A r_B } = \sqrt{ r_B r_C } + \sqrt{ r_C r_A } $

が成立する．添字がだるいので $ r_A = a, \ r_B = b, \ r_C = c $ とする．

$ \sqrt{ ab } = \sqrt{ bc } + \sqrt{ ca } $

S(100) まで適当に調べる

(a, b, c) = 
(4, 4, 1), (8, 8, 2), (9, 36, 4), (12, 12, 3), (16, 16, 4), 
(18, 72, 8), (20, 20, 5), (24, 24, 6), (28, 28, 7), (32, 32, 8), 
(36, 36, 9), (40, 40, 10), (44, 44, 11), (48, 48, 12), (52, 52, 13), 
(56, 56, 14), (60, 60, 15), (64, 64, 16), (68, 68, 17), (72, 72, 18), 
(76, 76, 19), (80, 80, 20), (84, 84, 21), (88, 88, 22), (92, 92, 23), 
(96, 96, 24), (100, 100, 25)

よくよく考えれば GCD(a, b, c) = 1 の場合のみ調べて後でなんとかすればいい．
GCD(a, b, c) = 1 であるもののみ S(1000) まで調べると，

(a, b, c) =
(4, 4, 1), (9, 36, 4), (16, 144, 9), (25, 400, 16), (36, 900, 25),
(100, 225, 36), (441, 784, 144)

どうも a, b, c は平方数になるらしい．示してみる．

1. $ \mathrm{GCD}(a, b, c) = 1 $ かつ $ \sqrt{ ab } = \sqrt{ bc } + \sqrt{ ca } $ ならば $ a, \ b, \ c $ は平方数

$ \sqrt{ ab } = \sqrt{ bc } + \sqrt{ ca } \qquad \cdots (1)$

の両辺 2 乗して整理すると

$ 2c \sqrt{ ab } = ab - bc - ca $

右辺は整数だから左辺も整数で，このとき $ ab $ は平方数．

$ab$ が平方数となるのは，$a_0, b_0$ を整数として
　ⅰ. $a = a_0^{2}, \ b = b_0^{2} $
　ⅱ. ある素因数 $p \ (p > 1)$ が存在して，$ a = p a_0^{2}, \ b = p b_0^{2}$

ここでⅱを仮定して式$(1)$に代入すると，

$ a_0 b_0 p = b_0 \sqrt{ pc } + a_0 \sqrt{ pc } = (a_0 + b_0) \sqrt{ pc } $

左辺整数より右辺も整数で，このとき $ pc $ は平方数．
$ pc $ が平方数となるのは $ c = p c_0^{2} $ と表される場合に限るが，これは $ \mathrm{GCD}(a, b, c) = 1 $ に矛盾．

よってⅰが成立するのでこれを式$(1)$に代入すると，

$ a_0 b_0 = (b_0 + a_0) \sqrt{ c } $

左辺整数より右辺も整数だから，ある整数 $c_0$ を用いて $ c = c_0^{2} $．つまり $c$ は平方数．

以上より，$a, b, c$ は全て平方数．

よって $ a = a_0^{2}, \ b = b_0^{2}, \ c = c_0^{2} $ として，

$ a_0 b_0 = b_0 c_0 + c_0 a_0 $

なる組を拾ってくればいい．適当に拾ってくると，

(a0, b0, c0) =
(2, 2, 1), (3, 6, 2), (4, 12, 3), (5, 20, 4), (6, 30, 5),
(7, 42, 6), (8, 56, 7), (9, 72, 8), (10, 15, 6), (10, 90, 9),
(14, 35, 10), (18, 63, 14), (21, 28, 12), (22, 99, 18), (24, 40, 15),
(30, 70, 21), (33, 88, 24), (36, 45, 20), (44, 77, 28), (55, 66, 30),
(60, 84, 35), (78, 91, 42)

どうも $ a_0 - c_0 , \ b_0 - c_0 $ はともに平方数になるらしい．示してみる．

1. $ a_0 - c_0, \ b_0 - c_0 $ は平方数

$ a_0 b_0 = b_0 c_0 + c_0 a_0 $

を変形すると

$ (a_0 - c_0)(b_0 - c_0) = c_0^{2} \qquad \cdots (1) $

左辺が平方数となるのは，$ m, n $ を整数として
　ⅰ. $a_0 - c_0 = m^{2}, \ b_0 - c_0 = n^{2} $
　ⅱ. ある素因数 $p \ (p > 1)$ が存在して，$ a_0 - c_0 = p m^{2}, \ b_0 - c_0 = p n^{2} $

ここでⅱを仮定して式$(1)$に代入すると，$ c_0 = p m n $ となり，
$ \mathrm{GCD}(a, b, c) = \mathrm{GCD}(a_0^{2}, b_0^{2}, c_0^{2}) = 1 $ に矛盾．

さっきと似たような感じ．

で，$ a_0 - c_0 = m^{2}, \ b_0 - c_0 = n^{2}, \ (a_0 - c_0)(b_0 - c_0) = c_0^{2} $ を解くと，

$ a_0 = m ( m + n ), \ b_0 = n ( m + n ), \ c_0 = m n $

となり，$a_0, \ b_0, \ c_0 $ を生成する式が得られる．なんか原始ピタゴラス数の生成式を思い出す．

$m, \ n$に関する条件は，$ c_0 < a_0 \leq b_0 $ と $ \mathrm{GCD}(a_0, b_0, c_0) = 1 $ から，

$ 1 \leq m \leq n $
$\mathrm{GCD}(m, n) = 1$

で，これと $ b = b_0^{2} = m^{2} (m + n)^{2} \leq 10^{9} $ を満たすように $m, \ n$ のループを適当に回せばいい．

typedef long long ll;

ll GCD(ll m, ll n)
{
    ll k;
    do {
        k = m % n;
        m = n; n = k;
    } while (k != 0);
    return m;
}

int main()
{
    clock_t start, end;
    start = clock();

    const ll MAX = 1000000000;
    ll ans = 0;
    for (ll m = 1; 4 * m * m * m * m <= MAX; m++) {
        for (ll n = m; n * n * (m + n) * (m + n) <= MAX; n++) {
            if (GCD(m, n) != 1) continue;
            ll a0 = m * (m + n);
            ll b0 = n * (m + n);
            ll c0 = m * n;

            ll a = a0 * a0;
            ll b = b0 * b0;
            ll c = c0 * c0;

            ll k = MAX / b;

            ll s = (a + b + c) * k * (k + 1) / 2;
            ans += s;
        }
    }
    cout << ans << endl;

    end = clock();
    printf("%d msec.\n", end - start);
}

3 msec とか