筋の良い全探索とは（数学夏祭り問1）

数学夏祭り

出典

攻略法

結論から言うと、この問題の攻略法は、「コンピュータプログラムを組む」が適切であると考える。
ただし、ただ闇雲に組むだけでは駄目である

  1. なぜコンピュータプログラムを組むことが適切なのか
  2. どのようなプログラムを組めば良いのか

について解説していく
注意：途中で登場するソースコードは全てpython(3.x)で書かれたものです。

なぜコンピュータプログラムを組むことが適切なのか

それは普通に問題を解いてみればわかる

通分すると
$$pqr=79(p+q)$$
両辺に{r}を掛けて整理すると
$$(pr-79)(q r-79)=79^2$$
ここで、{p\leq q,r\gt 0}であるため、{pr-79\leq qr-79}となる。つまり以下の可能性だけが残る。（79が素数であることを利用）
$$(pr-79,q r-79)=(1,79^2),(79,79),(-79,-79),(-79^2,-1)$$
（追記：両方とも負の場合を列挙するのを忘れてました。しかしその場合は{pr\leq 0}となるため不適となります。以降の議論に影響は出ません。）
後半2つは不適であるため、以降は前半2つだけを考える。

{pr-79=qr-79=79}のときは、{pr=qr=158}より、{(p,q,r)=(79,79,2),(158,158,1)}である。
（注意：{r\lt 79}であるため、{r=79,158}の可能性は消える）

{pr-79=1,qr-79=79^2}のときは、{pr=80,qr=79\times 80}より{q=79p}となる。
{p}は80の約数を全列挙してみればよい。よって{p=1,80,2,40,4,10,5,16,8,10}となりそうだが、{r\lt 79}であるため、{(p,r)=(1,80)}を除外しなくてはいけない。ここを見落としてミスした人も少なくないと考える。よって{p}は80の正の約数のうち、1を除いたもの全てを見ればよい。そうすると、
$$\frac{1}{p}+\frac{1}{q}=\frac{80}{79p}=\frac{r}{79}$$
より、{p}が「80の約数のうち、1以外のもの」ならば、解の十分性も満たしている。

つまり、手計算で求めようとすると、{p=1}を除き忘れるというミスが起こりやすくなるのである。

プログラムを組めばいいという理由はもう一つある。このようにして解を全列挙した時、{p}について並べ替えなければいけない。すると
$$p=2,4,5,8,10,16,20,40,79,80,158$$
となり、{(pr-79,qr-79)}が{(79,79)}であることから生成できた解と{1,79^2}であることから生成できた解が混じっている。これらの解をミスなく並べ替えるためにはコンピュータプログラムに頼るのが適切である。

ちなみに答えは25280である。ここで、
「前から3番目の{r}」は16であり、（{(p,q,r)=(5,395,16)}）
「後ろから5番目の{q}」は1580となっている。（{(p,q,r)=(20,1580,4)}）

結局の所、このような先着5名が表彰されるようなタイムアタック形式で、このようなケアレスミスを誘発しやすい問題が出題された場合、コンピュータプログラムを組んで突破するのが適切である。電卓の使用はルールの範囲内であるため堂々と使って大丈夫である。

どのようなプログラムを組めば良いのか

コンピュータの使用がルール上問題ないと言っても、正しく使いこなせるようにならなくてはいけない。
解を探す上で、ぱっと思いつく方法は全探索である。しかし、ただ闇雲に全探索してはいけない。

まず、解かどうかを判定するプログラムに注意が必要である。

#ダメな方法
def isanswer(p,q,r):
    if 1/p+1/q==r/79 and p<=q and r<79:
        return True
    else:
        return False

#よい方法
def isanswer(p,q,r):
    if (p*q*r)==79*(p+q) and p<=q and r<79:
        return True
    else:
        return False

なぜ前者がダメで後者がOKなのかというと、浮動小数点演算の誤差を考慮してないからである。ここでは詳しい説明は省くが、要は「プログラミングで小数を扱う場合は誤差が生まれる」のである。そして==という演算子で判定するため、誤差が生まれると正しく判定できなくなる。
これを避けるためには、整数での演算に帰着する必要がある。（ここで正しく式変形できるかが重要になる）

次に全探索する範囲について考える必要がある。
前節での式変形を流用すると
$$(pr-79)(q r-79)=79^2$$
となる。ここで、{p,q,r}は整数であるため、以下の不等式が成り立つ。
$$(pr-79)^2\leq (pr-79)(q r-79)=79^2$$
$$ p-79\leq (pr-79)\leq 79$$
よって{1\leq p\leq 158}で探索すればよいことになる。

{q}については{p}とは別の方法を使う。
$$ q-79\leq (q r-79)\leq 79^2$$
より、
$$ q\leq 79^2+79$$
となる。よって{1\leq q\leq 79^2+79=79\times 80}で探索すればよい。

{r}については、{1\leq r\leq 78}で探索すればよい…というわけではない。
{r}は等式として候補が決まっているため、十分性を満たしているかの判定だけでよい。
具体的にはこうなる。

def isanswer(p,q):
    if (79*(p+q))%(p*q)==0:
        r=(79*(p+q))//(p*q)
        if p<=q and r<79:
            return True
    return False

つまり、{p,q}の2変数だけ決まれば{r}の値は自動的に決まるため、{r}について78個の自然数を探索する必要はない。

結局の所以下のようなプログラムを組めば答えを出力してくれる。

for p in range(1,159):
    for q in range(1,79*80+1):
        if (79*(p+q))%(p*q)==0:
            r=(79*(p+q))//(p*q)
            if p<=q and r<79:
                print(p,q,r)

出力結果は以下のようになる。

2 158 40
4 316 20
5 395 16
8 632 10
10 790 8
16 1264 5
20 1580 4
40 3160 2
79 79 2
80 6320 1
158 158 1

これでめでたくすべての解を出力してくれる。
探索範囲は{158\times 79\times 80=998560}である。この数字が{10^6}よりも小さいオーダーであるため、プログラムの実行時間は現実的な時間となる。これを{r}についても全探索すると{8\times 10^7}くらいの計算量になるため、プログラムの実行に少し時間がかかってしまう。

効率的なアルゴリズム

一般化してみる。

{N}を奇素数として、
$$\frac{1}{p}+\frac{1}{q}=\frac{r}{N}$$
を満たす自然数{p,q,r(p\leq q,r\lt N)}を全列挙せよ

という問題について考える。前述のアルゴリズムだと、計算量は{O(N^3)}となっている。今回出題された問題では{N=79}なので現実的な時間で終わったが、{N}が大きくなっても解けるようなアルゴリズムを考えてみる。

ここで、同様の手法によって
$$(pr,q r)=(N+1,N^2+N),(2N,2N)$$
となっている。{(pr,q r)=(2N,2N)}の場合には、{r=1,2}のみであるため、
{(p,q,r)=(2N,2N,1),(N,N,2)}となる。（{2N}の約数が{1,2,N,2N}の4つのみであることから従う。）
{N+1}については約数全列挙を使えばよい。
ここで、{N}の約数の全列挙は{O(\sqrt{N})}で行うことができる。
これは、{N}という正の約数{d}は{N/d}という約数とセットになっていて、
{d}と{N/d}のどちらかは{\sqrt{N}}より小さいため、{d\leq \sqrt{N}}だけチェックすればOKだからである。すると以下のようになる。

N=int(input())
i=1
while(i*i<=N):
    if N%i==0:
        print(i)
        if i*i!=N:
            print(N//i)
    i+=1

これで{N}の正の約数すべてを出力してくれる。

つまり、この問題を解くためには
①{N+1}の約数を全列挙（{O(\sqrt{N})}）
②{N+1}の約数{d}についてループを走らせる
③　　　{d}について、{(p,q,r)=(d,Nd,(N+1)/d)}が解の候補である。
④　　　{r\lt N}かどうかを確認する
⑤生成した解全てについて{p}でソートする

計算量は{O(\sqrt{N})}であるため、{N\leq 10^{12}}くらいなら現実的な時間で解いてくれる。*1
結局プログラムを書くと以下のようになる。

N=int(input())
D=[]
i=1
while(i*i<=(N+1)):
    if (N+1)%i==0:
        D.append(i)
        if i*i!=(N+1):
            D.append((N+1)//i)
    i+=1
ans=[(2*N,2*N,1),(N,N,2)]
for d in D:
    p,q,r=d,N*d,(N+1)//d
    if r<N:
        ans.append((p,q,r))
ans.sort()
for p,q,r in ans:
    print(p,q,r)

このプログラムは、例えば{N=998244353}を入力に入れても、プログラムの実行は現実的な時間で完了する。
一方、{O(N^3)}の方のプログラムに{N=998244353}を突っ込んだ場合、恐らく宇宙の年齢に匹敵する実行時間がかかるだろう。

*1:約数の個数は元の数よりも十分に小さいため、⑤はボトルネックにはならない。