受験数学をプログラミングで倒してみた

はじめに

しばらく前に大学受験を終えたluna_moonlightです。突然ですが、数学に苦しめられたという経験はありますか?おそらく多くの方は数学、特に受験数学に苦しめられたことがあるでしょう。そんな数学の問題にも、プログラミングを使えば一瞬にして答えがわかってしまうようなものがあります。この記事では、受験数学をプログラミングで倒してみなさんの無念を晴らす(?)とともに、そのプログラムの高速化の方法について解説していきます。

一橋大学 2020年 大問1(1)

$10^{10}$を$2020$で割った余りを求めよ。

数学解法

解法

$$\begin{split} &10^1 \equiv 10 \pmod {2020}\\ &10^2 \equiv 100 \pmod {2020}\\ &10^3 \equiv 1000 \pmod {2020}\\ &10^4 \equiv 10000 \equiv 1920 \pmod {2020}\\ &10^5 \equiv 10^4 \cdot 10 \equiv 1920 \cdot 10 \equiv 19200 \equiv 1020 \pmod {2020}\\ &10^6 \equiv 10^5 \cdot 10 \equiv 1020 \cdot 10 \equiv 10200 \equiv 100 \pmod {2020}\\ \end{split}$$

周期性より、$10^{10} \equiv 100 \pmod {2020}$

プログラミング解法

愚直

愚直に$10$を$10$回かけて、$2020$で割った余りを求めます。

a = 10
b = 10
mod = 2020

ans = 1
for _ in range(b):
    ans *= a
    ans %= mod

print(ans)

この方法だと、指数の肩を$N$として$O(N)$で答えが求まります。

高速化

$x^{2N}=x^N \cdot x^N$なので、この計算を繰り返し用いることで、指数の肩が2の冪乗の場合は高速に計算できます。なので、指数の肩を2の冪乗に分解することで高速に計算できます。分解方法は、2進数表記への変換と対応しています。このような方法を、繰り返し二乗法と呼びます。

実際に$10^{10}$を分解すると、

$$10^{10} = 10^{8 + 2} = 10^8 \cdot 10^2$$

このようになります。

指数の肩が2の冪乗になる値を求めると、

$$\begin{split} &10^2 \equiv 10^1 \cdot 10^1 \equiv 100 \pmod {2020}\\ &10^4 \equiv 10^2 \cdot 10^2 \equiv 100 \cdot 100 \equiv 1920 \pmod {2020}\\ &10^8 \equiv 10^4 \cdot 10^4 \equiv 1920 \cdot 1920 \equiv 1920 \pmod {2020} \end{split}$$

このようになるので、$10^{10} \equiv 10^8 \cdot 10^2 \equiv 1920 \cdot 100 \pmod {2020}$と求まります。

a = 10
b = 10
mod = 2020

ans = 1
while b > 0:
    if b % 2 == 1:
        ans *= a
        ans %= mod
    a *= a
    a %= mod
    b //= 2

print(ans)

この方法だと、指数の肩を$N$として$O(\log N)$で答えが求まります。

実行時間

参考までに、指数の肩の大きさによる実行時間の違いを載せておきます。

	$10^2$	$10^4$	$10^6$	$10^8$
愚直	$2.98 \cdot 10^{-5} \mathrm{s}$	$1.40 \cdot 10^{-3} \mathrm{s}$	$0.110 \mathrm{s}$	$9.88 \mathrm{s}$
繰り返し二乗法	$2.02 \cdot 10^{-5} \mathrm{s}$	$2.38 \cdot 10^{-5} \mathrm{s}$	$2.67 \cdot 10^{-5} \mathrm{s}$	$3.00 \cdot 10^{-5} \mathrm{s}$

一橋大学 2021年 大問1

$1000$以下の素数は$250$個以下であることを示せ。

数学解法

解法

$1000$以下の$2, 3, 5$の倍数の集合をそれぞれ$A, B, C$とする。
$1000$以下の素数$2,3,5$以外の$2,3,5$の倍数の集合$D$の個数は、

$$\begin{split} |D| &=|A \cup B \cup C| - 3\\ &= |A| + |B| + |C| - |A \cap B| - |B \cap C| - |C \cap A| + |A \cap B \cap C| - 3\\ &= \left\lfloor \dfrac{1000}{2} \right\rfloor + \left\lfloor \dfrac{1000}{3} \right\rfloor + \left\lfloor \dfrac{1000}{5} \right\rfloor - \left\lfloor \dfrac{1000}{2\cdot3} \right\rfloor - \left\lfloor \dfrac{1000}{3\cdot5} \right\rfloor - \left\lfloor \dfrac{1000}{5 \cdot 2} \right\rfloor + \left\lfloor \dfrac{1000}{2\cdot3\cdot5} \right\rfloor - 3\\ &= 500 + 333 + 200 -166 -66 -100 + 33 - 3\\ &= 731 \end{split}$$

$7,11,13,17,19,23,29$から$2$つ選んでかけたものは$1000$以下の合成数であり、$2,3,5$の倍数ではない。
ゆえに、${}_7 C_2=21$個は$D$に含まれない$1000$以下の合成数である。
以上より、$1000$以下の素数でない自然数は$752$個以上であるため、$1000$以下の素数は$250$個以下である。

別解

$n$と互いに素な$n$以下の自然数の個数を$\varphi(n)$^[1]とする。
$n$の素因数分解が$n=\prod_{i=1}^{d}p_i^{e_i}$と表されるとき、$\varphi(n)=n\prod_{i=1}^{d}\left(1 - \dfrac{1}{p_i} \right)$である。
ここで、$1050=2\cdot3\cdot5^2\cdot7$について考えると、

$$\varphi(1050)=1050\left(1 - \dfrac{1}{2}\right)\left(1 - \dfrac{1}{3}\right)\left(1 - \dfrac{1}{5}\right)\left(1 - \dfrac{1}{7}\right)=240$$

ゆえに、$1050$以下の素数の個数は、先ほどの結果に$2,3,5,7$を加えて$244$個以下である。
よって、明らかに$1000$以下の素数の個数は、$250$個以下である。

プログラミング解法

プログラミング解法では、$1000$以下の素数の個数を求めます。

愚直

$1000$以下のすべての自然数に対して、素数判定を行います。$N$が素数か否かの判定は、$2$から$N-1$までのすべての自然数で割り切ることができないか否かで判定することができます。

N = 1000

ans = 0
for i in range(1, N + 1):
    if i == 1:
        continue
    for j in range(2, i):
        if i % j == 0:
            break
    else:
        ans += 1

print(ans)

この方法だと、$N$以下の素数の個数は$O(N^2)$で求まります。

高速化1

先ほどの素数判定には無駄があります。$2$から$N-1$までのすべての自然数で割り切ることができないか判定していましたが、実は$2$から$\lfloor\sqrt N\rfloor$までで十分です。

N = 1000

ans = 0
for i in range(1, N + 1):
    if i == 1:
        continue
    for j in range(2, int(i ** 0.5) + 1):
        if i % j == 0:
            break
    else:
        ans += 1

print(ans)

この方法だと、$N$以下の素数の個数は$O(N\sqrt N)$で求まります。

高速化2

素数$p$が見つかったら、$p$を除く$p$の倍数をふるい落とすことで、合成数を消すことができます。これを次々と行うことで、最終的に素数だけを残すことができます。この方法をエラトステネスのふるいと呼びます。

試しに、$30$以下の素数をエラトステネスのふるいを使って求めてみましょう。
まず、$2$に素数ラベルをつけて、$2$を除く$2$の倍数をふるい落とします。

次に、$3$に素数ラベルをつけて、$3$を除く$3$の倍数をふるい落とします。

$4$はすでにふるい落とされているため飛ばします。
$5$に素数ラベルをつけて、$5$を除く$3$の倍数をふるい落とします。

この操作を繰り返し続けていくと、以下のようになります。

これで、$30$以下の素数は$2,3,5,7,11,13,17,19,23,29$であることがわかります。

N = 1000

ans = 0
P = [True] * (N + 1)
P[0] = False
P[1] = False
for i in range(N + 1):
    if P[i]:
        ans += 1
        j = 2
        while j * i <= N:
            P[j * i] = False
            j += 1

print(ans)

この方法だと、$N$以下の素数の個数は$O(N\log{\log N})$^[2]で求まります。
実は、高速化1と同じく、ふるい落とし操作は$\lfloor\sqrt N\rfloor$までで十分です。このことも考慮すると、少しだけ速くすることができます。

エラトステネスのふるいの計算量は$O(N\log{\log N})$でとても高速ですが、代わりに空間計算量が$O(N)$になってしまいます。

実行時間

参考までに、素数判定する範囲による実行時間の違いを載せておきます。

	$10^2$	$10^4$	$10^6$	$10^8$
愚直	$9.08 \cdot 10^{-5} \mathrm{s}$	$0.299\mathrm{s}$	$1.98 \cdot 10^{3} \mathrm{s}$	-
$\sqrt N$素数判定	$5.70 \cdot 10^{-5} \mathrm{s}$	$9.03 \cdot 10^{-3} \mathrm{s}$	$3.31 \mathrm{s}$	$2.37 \cdot 10^{3} \mathrm{s}$
エラトステネスのふるい	$4.20 \cdot 10^{-5} \mathrm{s}$	$3.42 \cdot 10^{-3} \mathrm{s}$	$0.346 \mathrm{s}$	$63.3 \mathrm{s}$

京都大学 2022年 大問1

$5.4<\log_4{2022}<5.5$を示せ。ただし、$0.301<\log_{10}{2}<0.3011$

数学解法

解法

$2000<2022<2048$より、$\log_4{2000}<\log_4{2022}<\log_4{2048}$

$$\log_4{2048} = \log_4{2^{11}} = 11\log_4{2} = 11\cdot\dfrac{\log_2{2}}{\log_2{4}} = 11\cdot\dfrac{1}{2} = 5.5$$

$$\begin{split} \log_4{2000} = \log_4{(2\cdot10^3)} = \log_4{2} + \log_4{10^3} = \log_4{2} + 3\log_4{10}\\ = \dfrac{\log_2{2}}{\log_2{4}} + 3 \cdot \dfrac{\log_{10}{10}}{\log_{10}{4}} = \dfrac{1}{2} + 3 \cdot \dfrac{1}{2\log_{10}{2}}\\ > \dfrac{1}{2} + 3 \cdot \dfrac{1}{2 \cdot 0.3011} > 5.4 \end{split}$$

以上より、$5.4<\log_4{2022}<5.5$

プログラミング解法

多くのプログラミング言語では、$\log$の計算がサポートされているため、それを使えば一発でわかります。しかしこの記事では、高速化の方法について解説したいので、$\log$を使わずに解いていきます。
また、有理数乗がサポートされているという前提で進めていきます。

プログラミング解法では、$4^x=2022$となる$x$の近似値を求めます。なぜなら、$\log_4{2022}=x$とおくと、$4^x=2022$となるためです。
また、$x$の存在範囲はおよそわかっているものとします。

愚直

$x$が$L$から$R$までの範囲に含まれているとします。$L$から$R$まで$0.01$刻みで$4^x$を計算して、$2022$を下回る最大値と上回る最小値を記録します。

from decimal import Decimal

a = 4
b = 2022

L = 0
R = 10
l_ans = 0
r_ans = 10
x = Decimal(L)
while L <= x <= R:
    if 4 ** x < b:
        l_ans = max(l_ans, x)
    else:
        r_ans = min(r_ans, x)
    x += Decimal('0.01')

print(f'{l_ans}<log_{a}({b})<{r_ans}')

この方法だと、$x$の存在範囲の大きさを$N$として$O(N)$で近似値が求まります。

高速化

先ほどの方法には無駄があります。$4^x$は単調増加するので、小さい値から順番に$4^x$が$2022$よりも小さいかどうかを確かめるのは無駄が多いです。探索範囲の半分の場所で$4^x$が$2022$よりも小さいかどうかを確かめれば、探索範囲を半分にすることができます。この方法を二分法と呼びます。

実際に、二分法を行なっている様子を見てみましょう。
まず、$x$の存在範囲の左端を$L=0$、右端を$R=10$とします。

$x$の存在範囲の中心 $\text{mid}=5$について、$4^{\text{mid}}$は$2022$よりも小さいので、$L=\text{mid}=5$とします。

次に、$x$の存在範囲の中心 $\text{mid}=7.5$について、$4^{\text{mid}}$は$2022$よりも大きいので、$R=\text{mid}=7.5$とします。

さらに、$x$の存在範囲の中心 $\text{mid}=6.25$について、$4^{\text{mid}}$は$2022$よりも大きいので、$R=\text{mid}=6.25$とします。

この操作を繰り返していくと、$x$の存在範囲を半分ずつ縮めながら$x$の近似値を求めることができます。

from decimal import Decimal

a = 4
b = 2022

l_ans = Decimal('0')
r_ans = Decimal('10')
while (r_ans - l_ans) > Decimal('0.01'):
    mid = (l_ans + r_ans) / 2
    if 4 ** mid < b:
        l_ans = mid
    else:
        r_ans = mid

print(f'{l_ans}<log_{a}({b})<{r_ans}')

この方法だと、$x$の存在範囲の大きさを$N$として$O(\log N)$で近似値が求まります。

実行時間

参考までに、$x$の存在範囲の大きさによる実行時間の違いを載せておきます。

	$10^1$	$10^2$	$10^3$	$10^4$	$10^5$
愚直	$5.80 \cdot 10^{-2} \mathrm{s}$	$0.562 \mathrm{s}$	$5.65 \mathrm{s}$	$57.2 \mathrm{s}$	$575 \mathrm{s}$
二分法	$7.92 \cdot 10^{-4} \mathrm{s}$	$1.01 \cdot 10^{-3} \mathrm{s}$	$1.30 \cdot 10^{-3} \mathrm{s}$	$1.43 \cdot 10^{-3} \mathrm{s}$	$1.64 \cdot 10^{-3} \mathrm{s}$

まとめ

この記事では、受験数学をプログラミングで倒すとともに、そのプログラムの高速化の方法について解説しました。これらの高速化は、競技プログラミングでもよく使うテクニックなので、ぜひ他の例でも試して習得してみてください。

脚注

1. オイラーのトーシェント関数と呼ばれる。

2. 素数の逆数和$\dfrac{1}{2}+\dfrac{1}{3}+\dfrac{1}{5}+\dfrac{1}{7}+...$は$\log{\log N}$漸近することが知られている。