数学の力

京大生が数学の定理・公式の証明や入試問題の解説をするブログ.

プロジェクトオイラー157. ~ Solving the Diophantine equation 1/a+1/b=p/10^n ~


スポンサードリンク

No.157 Solving the Diophantine equation 1/a+1/b=p/10^n

前回の記事で  1/a+1/b=p/10 という問題について書きました(→方程式 1/a+1/b=p/10).

今回の問題はさらに右辺の分母が  10^n になったことで少し難しくなっています.

 

問題.

 a, b, p, n を正の整数, a\leq b として,ディオファントス方程式  \dfrac{1}{a}+\dfrac{1}{b}=\dfrac{p}{10^n} を考える. n=1 のとき,この方程式は20個の解をもつ.(具体的な解は前回の記事を御覧ください)
 1\leq n\leq 9 におけるこの方程式の解は全部でいくつあるか.



基本的な解き方(考え方)は前の記事に書いた  n=1 の場合と同様ですが,少し工夫しないと計算時間が足りなくなります.



考え方とプログラム例(Python3).

まず,方程式の両辺に  10^npab をかけて,変形します.

\begin{align*}
10^npa+10^npb &=  p^2ab\\
p^2ab-10^npa-10^npb+10^{2n} &=  10^{2n}\\
(pa-10^n)(pb-10^n) &=  10^{2n}
\end{align*}


右辺の  10^{2n} 10^{2n}=PQ, と自然数の積の形に書けるとすると,

\begin{align*}
pa &=  10^n+P\\
pb &=  10^n+Q
\end{align*}

となります. p 10^n+P, 10^n+Q の公約数です.



そこで, 10^n+P, 10^n+Q の公約数の個数を数え上げればいいことが分かります.



 n を固定したときの解の個数を  f(n) とすると,

\begin{align}
f(n) = \sum_{\substack{PQ=10^{2n}\\P\leq Q}} \mathrm{d}(\gcd(10^n+P, 10^n+Q))
\end{align}

と書けます.


 nが大きい場合について普通に公約数を数えようとすると大変なので,もう少し解析を進めます.



 10^{2n}の素因数は2と5のみなので, P=2^{x_1}5^{x_2}, Q=2^{2n-x_1}5^{2n-x_2} とすると,

\begin{align*}
10^n+P &=  2^{\min(x_1, n)}5^{\min(x_2, n)}\left\{2^{\max(n-x_1, 0)}5^{\max(n-x2, 0)}+2^{\max(x_1-n, 0)}5^{\max(x_2-n, 0)}\right\}\\
10^n+Q &=  2^{\min(2n-x_1, n)}5^{\min(2n-x_2, n)}\left\{2^{\max(x_1-n, 0)}5^{\max(x2-n, 0)}+2^{\max(n-x_1, 0)}5^{\max(n-x_2, 0)}\right\}
\end{align*}

となるので,


\begin{align*}
\gcd(10^n+P, 10^n+Q) &=  2^{\min(x_1, 2n-x_1)}5^{\min(x_2, 2n-x_2)}\\&\quad\times\left\{2^{\max(n-x_1, 0)}5^{\max(n-x2, 0)}+2^{\max(x_1-n, 0)}5^{\max(x_2-n, 0)}\right\}
\end{align*}

これで  \gcd(10^n+P, 10^n+Q)素因数分解が途中までできている形になっているので,約数の個数が数えやすくなります.

import math
def number_of_divisors(A, B, N):
    ret = 1
    n = 2
    while N != 1 or n <= 5:
        m = 0
        while N % n == 0:
            m = m + 1
            N = N / n
        if n == 2:
            ret = ret * (1 + A + m)
        elif n == 5:
            ret = ret * (1 + B + m)
        else:
            ret = ret * (1 + m)
        n = n + 1
    return ret

def f(x1, x2):
    return number_of_divisors(min(x1, 2*n-x1), min(x2, 2*n-x2), 2**max(x1-n, 0)*5**max(x2-n,0)+2**max(n-x1,0)*5**max(n-x2,0))

ans = 0
for n in range(1, 10):
    for x1 in range(n):
        for x2 in range(2*n+1):
            ans = ans + f(x1, x2)
    for x2 in range(n+1):
        ans = ans + f(n, x2)
print(ans)