OMC020 プログラミングを用いた解法

# はじめに
数学コンテストサイト「OnlineMathContest」の「OMC020(とある数学のコンテスト)」において、B, C, D, F問題を解く際に、プログラミング (Python) を利用しました。
そのコードと解法・解説を記します。
(コンテスト終了後に、より良さそうな方法やより分かりやすそうな記法などを発見した場合には、そちらを記載しています。)

<b>OMC020(とある数学のコンテスト)</b>
https://onlinemathcontest.com/contests/omc020

<b><font color="Red">注意: 2021年8月より、OnlineMathContest のコンテストでは、プログラミングを含む外部ツールの使用が禁止される予定ですので、ご注意ください。</font></b>
(本記事は、外部ツールの使用が許可されていたコンテストについてのものです。)

# B問題
問題
https://onlinemathcontest.com/contests/omc020/tasks/125

光線の反射回数を $n$ とし、 $m = n + 1$ とします。
解説にあるように、 $m$ と互いに素であるような $1$ 以上 $m$ 以下の整数の個数を数えます。
ただし、二つの整数 $i$ と $m-i$ は片方のみをカウントします。
Pythonでは以下のようにmath.gcdを用いると便利です。

```Python
from math import gcd

n = 2020

m = n + 1
S = set()
for i in range(1, m+1):
    if (not (m-i in S)) and gcd(m, i) == 1:
        S.add(i)

print(len(S))
```

また、「 $m$ と $i$ が互いに素かどうか」と「 $m$ と $m-i$ が互いに素かどうか」の答えは一致します。そのため、以下のようにforループの範囲を半分にすることで、 $i$ と $m-i$ の重複カウントを気にする必要がなくなり、より単純になります。

```Python
from math import gcd

n = 2020

m = n + 1
ans = 0
for i in range(1, m//2 + 1):
    if gcd(m, i) == 1:
        ans += 1

print(ans)
```

なお、問題には $n=4$ の場合に答えが $2$ となることが示されているので、提出前にそのケースで出力が正しくなるかを確認できます。

# C問題
問題
https://onlinemathcontest.com/contests/omc020/tasks/126

最大の数字のインデックス $N$ 、最大の数字 $i$ 、その一つ左の数字 $j$ を固定して考えると、「良い数」は以下のような数列 $b, c$ を合わせたものとして考えることができます。
(図では問題文と異なり、 $a$ の添字を $0$ からにしています。)
![図1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/420990/880ba510-26d9-9421-ea77-235df6ab0db9.png)
この「良い数」が $3$ の倍数であることは、数列 $b$ と数列 $c$ の総和が $3$ の倍数であることと同値です。
（正の整数 $n$ が $3$ の倍数であるかは、$n$ の各位の数の総和が $3$ の倍数であるかを調べることで判定できます。）
そのため、$(数列 b の総和を 3 で割った余り, 数列 c の総和を 3 で割った余り)$ が $(0, 0), (1, 2), (2, 1)$ のどれかである場合は (またそのような場合に限り) 、その「良い数」が $3$ の倍数といえます。
よって、
「 $1$ から $9$ までの数からの $p$ 通りの選び方の中で、最大値が $q$ であり、総和を $3$ で割った余りが $r$ であるものの個数」
「 $0$ から $9$ までの数からの $p$ 通りの選び方の中で、最大値が $q$ であり、総和を $3$ で割った余りが $r$ であるものの個数」
をあらかじめ求めておき、 $(N,i,j)$ の組ごとに $(数列bのパターン数) \times (数列cのパターン数)$ を足し合わせていくと答えが求まります。
Pythonでの組合せの全列挙には、itertools.combinationsが便利です。

```Python
from itertools import combinations

# B[p][q][r]は「1から9までの数からのp通りの選び方の中で、最大値がqであり、総和を3で割った余りがrであるものの個数」
B = [[[0] * 3 for i in range(10)] for j in range(11)]
# C[p][q][r]は「0から9までの数からのp通りの選び方の中で、最大値がqであり、総和を3で割った余りがrであるものの個数」
C = [[[0] * 3 for i in range(10)] for j in range(11)]

# Bの中身をすべて求める
for i in range(1, 10):
    L = list(combinations(range(1, 10), i))
    for j in range(len(L)):
        B[i][max(L[j])][sum(L[j]) % 3] += 1

# Cの中身をすべて求める
for i in range(1, 11):
    L = list(combinations(range(10), i))
    for j in range(len(L)):
        C[i][max(L[j])][sum(L[j]) % 3] += 1

# 答えを求める
ans = 0
R = [(0, 0), (1, 2), (2, 1)]  # 数列b,cの3で割った余りのペア
for N in range(15):
    lb = N  # 数列bの長さ
    lc = 14 - N + 1  # 数列cの長さ
    if 1 <= lb <= 9 and 1 <= lc <= 10:
        for i in range(10):
            for j in range(i):
                for k in range(3):
                    bb, cc = R[k]
                    ans += B[lb][j][bb] * C[lc][i][cc]
print(ans)
```

解説にある解法のほうが簡明なのですが、この方法には問題が複雑であっても解けるという利点があります。
後半の4重ループの回る回数の積が $10^6$ 以下くらいなら問題なく、例えば問題が「 $16$ 進法表記で $25$ 桁の良い数であって、 $15$ の倍数であるものを求めよ」であっても計算できます。（後半のループ回数は高々 $25 \times 16 \times 16 \times 15=96000$ 回）
なお、 $16$ 進数の場合は $3,5,15$ の倍数で、各位の和を見る判定法が使えます。

# D問題
問題
https://onlinemathcontest.com/contests/omc020/tasks/127

解説にあるように、等差数列の公式より $(a+b)(b-a+1) = 2N$ を満たす、 $1≦a<b$ であるような $(a, b)$ の個数が $2021$ 個である数を探せばよく、そのうち $5$ 番目に小さい数が答えとなります。
さらに考察を進められますが、この時点で「 $N$ を決めた場合の $(a, b)$ の個数」を求めることができます。
$xy = 2N$ となるような $(x, y)$ $(x≦y)$ について、$(x, y)=(a+b,\; b-a+1)$ と $(x, y)=(b-a+1,\; a+b)$ をそれぞれ解いたときに適切な解 $(a, b)$ が得られるか確認する、という方法を利用します。
以下のコードで「$N=1, \cdots , n$ に対する $(a, b)$ の個数」を出力できます。

```Python
n = 10000

for i in range(1, n):
    i2 = 2 * i
    ans = 0
    for j in range(2, int(i2**0.5) + 1):
        if i2 % j == 0:
            x = j
            y = i2 // j
            A = [x, y]
            # i2 = x * y となるようなx,yが求められた。
            su = x + y  # x + y = (a+b) + (b-a+1) = 2*b + 1
            if su % 2 == 1:
                b = (su-1) // 2
                # A: x = a+b の場合
                a = x - b
                if 1 <= a < b:
                    ans += 1
                # B: y = a+b の場合(x=yの場合はAと重複になるので行わない)
                if x != y:
                    a = y - b
                    if 1 <= a < b:
                        ans += 1
    print(i, ans)
```

$N$ が増えるにつれて $(a, b)$ の個数が多いものも出てきますが、それでも $N=10^4$ までの最大値で $23$ 個です。
そのため、このコードだけでは答え ( $2021$ 個であるもののうち $5$ 番目に小さい数) は求められなそうとわかります。
( $N=k$ までを調べる時間計算量は $O \( k^{\frac{3}{2}} \)$ であるため、一般のPCでは $N=10^6$ くらいまでの探索が限度となります。)
そのため、解説にあるような考察を進めることになります。

答えが $2^a 3^{336} {p_2}^5$ または $2^a 3^{336} {p_2}^2 {p_3}$ ( $a$ は $0$ 以上の整数、 $p_2, p_3$ は $5$ 以上の素数) と表されるというところまで考察が進むと、再びプログラミングが使えます。
$2^a {p_2}^5$ と $2^a {p_2}^2 {p_3}$ を、 $a≦5, \; p_2, p_3≦23$ くらいまで全列挙し、そのうち $5$ 番目に小さいものに $3^{336}$ をかけたものが答えです。

```Python
P = [5, 7, 11, 13, 17, 19, 23]  # 5以上の素数
ANS = []  # 答え(を3**336で割ったもの)の候補

for i in range(len(P)):
    for a in range(5):
        ANS.append(2**a * P[i]**5)

for i in range(len(P)):
    for j in range(len(P)):
        if i != j:
            for a in range(5):
                ANS.append(2**a * P[i] * P[j]**2)

ANS.sort()

ans = ANS[4]  # 5番目に小さいもの＝最終的な答えを3**336で割ったもの
for i in range(336):
    ans = (ans * 3) % 1000
print(ans)
```

また、ここまでの推論が正しいかについて、最初のコードを使用して確認することができます。
例えば、 $(a, b)$ の個数が $5$ 個である数は、同じように考えると $2^a {p_1}^5$ または $2^a {p_1}^2 {p_2}$ ( $a$ は $0$ 以上の整数、 $p_1, p_2$ は $3$ 以上の素数) と表されることになります。
実際に $45=2^0 \cdot 3^2 \cdot 5,\;\; 63=2^0 \cdot 3^2 \cdot 7,\;\; 75=2^0 \cdot 5^2 \cdot 3,\;\; 90=2^1 \cdot 3^2 \cdot 5,\;\; \cdots, \;\; 243=2^1 \cdot 3^5, \;\; \cdots$ といった数が出力されるため、推論は合っていそうだといえます。

# F問題 (解いてはいません)
問題
https://onlinemathcontest.com/contests/omc020/tasks/129
※F問題はコンテスト終了後に修正があったため、コードと説明もそれに合わせて修正を行ったものとなっています。

最終的な解答は得られないのですが、 $M$ の近似値を求めることができます。
「乱数で $100$ 以上の実数を $4$ 個生成し、『そのうち $2$ 数 $x,y$ に対して計算した、不等式をギリギリ満たす $m$ の値』の<b>最大値</b> ( $m_\mathrm{max}$ とする) を求める」
という試行を何回も行い、 $m_\mathrm{max}$ の<b>最小値</b>をとったものが $M$ の近似値です。
$4$ 数から $2$ 数は自由に選べるため、「」内では最大値をとりますが、 $100$ 以上の実数から $4$ 数が選ばれる際は自由に選べない (最悪の場合＝ $m_\mathrm{max}$ が小さくなってしまう場合、を想定する必要がある) ため最小値をとる、というのがポイントです。

100以上の実数を生成する方法はいろいろありますが、下記のコード内の方法だと、整数部分が $2,3,\cdots,11$ 桁である数がそれぞれ $\frac{1}{10}$の確率で得られます。 (random.random()は $0$ から $1$ の一様乱数を返します。)

```Python
from random import random

ans = 10**9
for i in range(10**5):
    A = [10 ** (random() * 10 + 2) for _ in range(4)]
    max_m = 0
    for i in range(4):
        for j in range(i+1, 4):
            x = A[i]
            y = A[j]
            l = (x*y + 1) ** 2  # 左辺
            r = (x**2 + 1) * (y**2 + 1)  # 右辺
            m = l / r
            max_m = max(max_m, m)
    ans = min(ans, max_m)

print(ans)
```

実際の答えは $M = 0.99998888...$ という値であり、上記のコードを $10$ 回実行した平均値は $0.99999216...$ でした。
なお、解答は小数で提出するわけではないため、この誤差を小さくしていくことで正解をすることはできず、正解するためには解説に書かれた考察を行うことになります。