文系卒が数学オリンピックをPythonで解く（2021年予選編）

#これまでのあらすじ

• 文系卒の私。でも数学になじみたい。
• ふむふむ、数学オリンピックというのがあるらしい。やってみますか。
• 紙とペンを使ってマジメに解く文系卒。えらい、えらいぞ！
• ５分後にふと衝撃走る。「これPythonでやったほうが早くない？」
• そんな、いつか 誰かに 本気で怒られそうな気づきを 実行したのであった──。

#ルール

• 数学オリンピックの問題をPythonで解く。
• ライブラリは可能な限り使わない。例外的にitertoolsは使う。
• 図形の問題は挑戦しないかも。許してね。
• 式変形すれば解けるようなものも面白味がないので対象外。

#引用元
問題はこちらから引用しております。

第31回(2021年)JMO予選の問題 - 数学オリンピック

#Q1
##問題

互いに素な正の整数m,nが m + n = 90を満たすとき、積mnとしてありうる最大の値を求めよ。
第31回(2021年)JMO予選の問題 - 数学オリンピック

####考える
数学に疎いので恥ずかしながら「互いに素」をググりました。

二つの整数 a, b が互いに素（たがいにそ、英: coprime, relatively prime, prime to）であるとは、a, b を共に割り切る正の整数が 1 のみであることをいう。
互いに素 - Wikipedia

なるほど！
「最大公約数が１である」ということですね。

####仕様

• 以下の２つを満たす整数m,nを網羅する
• 90 - m = nである。
• m,nそれぞれの約数を取得して、その積集合が１のみになる。
• 網羅した整数m,nの積mnを求め、最大値を算出する。

####コーディング
まずは、約数を求める関数を作りましょう。

#約数を取得する関数。set型で返す。
#その数の２番目に大きな約数は、その数の半分以下になる。

def get_divisors(num):
    divisors = {num}
    for _ in range(1, (int(num/2)+1)):
        if num % _ == 0:
            divisors.add(_)
    return divisors

candidatesというリストに、m,nの候補をタプルで格納していきます。
mの約数の集合とnの約数の集合の積集合（共通の要素）が1つならば「互いに素」としています。
(積集合の最大値(=最大公約数)が1の方が適切ですが、見逃してください。)

candidates = []
for m in range(1,91):
    n = 90 - m
    m_divisors = get_divisors(m)
    n_divisors = get_divisors(n)
    if len(m_divisors & n_divisors) == 1:
        candidates.append((n,m))
print(candidates)

# 実行結果 -> [(89, 1), (83, 7), (79, 11), (77, 13), (73, 17), (71, 19), (67, 23), (61, 29), (59, 31), (53, 37), (49, 41), (47, 43), (43, 47), (41, 49), (37, 53), (31, 59), (29, 61), (23, 67), (19, 71), (17, 73), (13, 77), (11, 79), (7, 83), (1, 89)]

整数m,nの候補が網羅できました。
いよいよmnの最大値を求めましょう。

nm_products = [n*m for n,m in candidates]
print(max(nm_products))

# 実行結果 -> 2021

答え：2021

#Q4
##問題

• 黒板に3つの相異なる正の整数が書かれている。

• 黒板に実数a,b,cが書かれているとき、それぞれを(b+c)/2,(a+b)/2,(b+c)/2に同時に書き換えるという操作を考える。
• この操作を2021回行ったところ、最後に黒板に書かれた３つの数はすべて正の整数だった。
• このとき、最初に書かれていた正の整数の和としてありうる最小の値を求めよ。

第31回(2021年)JMO予選の問題 - 数学オリンピック

####考える
問題文が難しい。
３つの数字を要素に持つリストを、他の数の平均にそれぞれ置き換える、という関数を──ん、待てよ？　他の数の平均？　1,1,1を入れると1,1,1が戻り値になるってこと？　
引数と戻り値が同じになるってことじゃない？　1,1,1が最小じゃないか？
と思いきや相異なる正の整数の記載があるのでそんなに簡単ではない。
いや、最初の相異なる正の整数がなんであれ、何回操作しても１未満にはならないから、それって・・・。

わかっちゃったけど書きますね。フフン、私は意外とマジメなんです。フフン。

####仕様

• 相異なる正の整数a,b,cを用意する。
• a,b,c = (b+c)/2,(a+b)/2,(b+c)/2の操作を2021回繰り返す。
• 操作後のa,b,cがすべて整数かどうかを判定する。条件に該当しなければ、別の相異なる正の整数a,b,cを用意する。

####コーディング

3つの相異なる正の整数を作るために組み合わせを使いたいのでitertoolsをimportします。

#イテレータツールちゃんを解禁
import itertools as itr

操作を関数化します。
（全体のことを考えると、再帰的な関数にしたほうがいいのかもしれない。）

def rewrite_blackboard(a,b,c):
    return (b+c)/2, (a+c)/2, (a+b)/2

さて、やりますか。

#どこまでやるか？を決めておく。
n = 8

#3つの相異なる正の整数a,b,cを作るため、rangeから組み合わせを取得。
for a,b,c in itr.combinations(range(1,n), 3):
    #操作前のa,b,cを保存。
    nums_before = [a,b,c]
    print('操作前:',a,b,c)
    i = 0
    while True:
        a,b,c = rewrite_blackboard(a,b,c)
        i += 1
        if i == 2021:
            break
    print('操作後:',a,b,c)
    #整数かどうかを判定
    #要素の数とTrueの数が一致するか（すべての要素がTrueか）
    if sum([(_ % 1 == 0) for _ in [a,b,c]]) == len([a,b,c]):
        break
        
print('最初に書かれていた正の整数の和としてありうる最小の値:',sum(nums_before))
# 実行結果 -> 操作前: 1 2 3
#            操作後: 2.0 2.0 2.0
#            最初に書かれていた正の整数の和としてありうる最小の値: 6

1,2,3が条件を満たす最小の組み合わせでしたので、答えは6。

答え：6

#Q6
##問題

正の整数nに対して、正の整数mであってmとnが互いに素であり、m+1とn+1も互いに素となるようなもののうち最小のものをf(n)で表す。このとき、f(1),f(2),...f(10**10)のうちに現れる正の整数は何種類あるか。
第31回(2021年)JMO予選の問題 - 数学オリンピック

###初見
####考える
問題文が難しいパート２。
m = f(n)ってことでいいのかな。
そして、nを1から10**10まで試して、通算でmのユニーク数がわかればよい、のかな。

ところで、10**10って100億？　100億回計算するのか？
何かもっと上手い解き方があるのでしょうが、ひとまずPythonで殴りましょう・・・。
　
####仕様

• 引数nに対して、以下の条件を満たす戻り値mを返す関数を定義する。
  • 以下の条件を満たす最小のm。
    • mは正の整数
    • mとnが互いに素
    • m+1とn+1が互いに素
• 空のsetを作成し、関数に1から10**10までの引数を与え、戻り値をsetに格納していく。
• setの要素数を数える。

####コーディング（失敗）
「互いに素」というキーワードが出てきましたので、Q１で使った約数を求める関数を流用します。

#約数を取得する関数。set型で返す。
#その数の２番目に大きな約数は、その数の半分以下になる。
def get_divisors(num):
    divisors = {num}
    for _ in range(1, (int(num/2)+1)):
        if num % _ == 0:
            divisors.add(_)
    return divisors

次に、関数fを定義。
もっと上手い書き方があるはずだ。。。絶対に。。。

def f(n):
    m = 0
    while True:
        m += 1
        if (len(get_divisors(n) & get_divisors(m)) == 1):
            if (len(get_divisors(n+1) & get_divisors(m+1)) == 1):
                break
    return m

100億回、スタート。
地球誕生は約46億年前──。
冥王星と太陽の平均距離、約59億km──。
世界総人口、約78億人──。

m_uniques = set()
for _ in range(1,10**10+1):
    m_uniques.add(f(_))
print(m_uniques)
print(len(m_uniques))

いつまでも終わらない悠久の時──。

###再挑戦
####考え直す
私の環境ではいつまでも終わらなかったので、作戦を変更します。
もちろん私の自作関数や実装がイケていないのもあるのですが、根本的に解き方を間違えている気がします。
Pythonなんぞ使わず紙とペンで解いている人たちがいるのですから、何か法則性があるはずです。
少ない回数で実行して、結果を見てみることにしました。

for n in range(1,10+1):
    print(n,f(n))
# 実行結果
1 2
2 1
3 2
4 1
5 4
6 1
7 2
8 1
9 2
10 1

なるほど、nが偶数ならばf(n)は１になるようです。
考えてみればnが偶数なら、n+1は奇数です。
m+1が２(最小値)だとしても、n+1とm+1は互いに素になります。
nが偶数のときの処理を省略すれば、処理回数を半分にできますね。

奇数のときだけ計算するようにしてみました。

for n in range(1,10**3+1):
    if n % 2 != 0:
        print(n,f(n))
# 実行結果
1 2
3 2
5 4
7 2
9 2
11 4
13 2
15 2
17 4
19 2
21 2
23 4
25 2
27 2
29 6

また法則がありますね。2,2,4,2,2,4,...を繰り返しています。
４になるのは、nが6k-1のとき。
おや、n=29のときは6になっていますね。

f(n)が新しい値をとるのは・・・

• nが2kのとき（厳密には2k-1のとき？）
• nが(2*3)k-1のとき
• nが(235)k-1のとき

法則が見つかりました！！！！！
この次はnが(2*3*5*7)k-1のとき、つまり209のときだ！

# 実行結果
205 2
207 2
209 10
211 2
213 2
215 4

当たってるぞ！！！！

つまり、nが100億になるまでに、(2*3*...*X)k-1がいくつ登場するかを調べればいいっぽい。
※本当にそうか証明しなくてはならないとは思うのですが、数学オリンピック予選は途中式を求められないので、レギュレーション準拠でやります。

####仕様ver2

• 列挙した素数すべての積が百億を超えるまで、素数を列挙する。
• 列挙した素数の個数を数える。

####コーディング（成功）
素数列挙の関数を定義します。(おそらくもっとイケてる書き方がある。)

def get_primes(N):
    primes = []
    for i in range(2, N + 1):
        primes.append(i)
        for p in range(2, i):
            if i % p == 0:
                primes.remove(i)
                break
    return primes

列挙した素数すべての積が百億を超えるまで、素数を列挙します。

N = 2
while True:
    N += 1
    primes = get_primes(N)
    prod = 1
    for _ in primes:
        prod *= _
    if prod > 10**10:
        break

列挙した素数の個数を数えます。
正確にはユニーク数（種類数）ですが、素数は重複しないのでlenでOK。

print(primes)
print(len(primes),'種類')

# 実行結果
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31]
11 種類

答え：11種類

#続く──。