全射（Surjection）の個数の求め方と高速化

今回は全射（Surjection）の個数の求め方を使って問題を解いて行きます。まず例題です。

### 【例題１】{a,b}の2文字を使って長さ３以下の単語を作ったときに{a,b}が両方とも含まれている単語はいくつあるでしょう？

3文字以下ということなので、順番に考えて行けば答えがでます。

|n (単語の長さ) |　単語 |　個数 |
|:---:|:---:|:---:|
| 1 | - | 0 |
| 2 | ab, ba | 2 |
| 3 | aab,aba,baa, baa,bab,bba| 6 |
| 合計| | **8** |

こんな小さな数でも、それなりの手間がかかるので、数が大きくなったらこれを計算する公式とかが欲しくなりますね。

### 全射（Surjection）を使って考える

この問題を文字の位置番号の集合$X$から文字の集合$Y$への写像と考えます。例えば aabなら、1  → a, 2  → a, 3  → bという感じです。

```math
\begin{align}
&すると \\
(元の文字が全部&含まれる) = (写像が全射)
\end{align}
```
と問題を置き換えることができます。

### 全射の個数の公式

そこで全射の個数の公式を探してみると、ありました。
[全射の個数の証明とベル数 (高校数学の美しい物語)](https://manabitimes.jp/math/892) によれば、

```math 
\begin{align}
n個の集合&からk個の集合への全射の数は \\
&\sum_{i=0}^{k}(-1)^{k-i}\binom{k}{i}i^n \\
&で表されます
\end{align}
```
### 全射の個数の公式を使う

これを使って例題１を解き、同じ答えが得られました。

```python 
import sympy as sp
def nSurj(n, k):
  return sum((-1)**(k-i) * sp.binomial(k,i) * i**n for i in range(0,k+1))

answer = 0
for n, a in [(1,2),(2,2),(3,2)]:
  ans = nSurj(n, a)
  print(f"n = {n}, a = {a}, ans = {ans}")
  answer += ans
print(f"Answer: {answer}")
# n = 1, a = 2, ans = 0
# n = 2, a = 2, ans = 2
# n = 3, a = 2, ans = 6
# Answer: 8
```

### 【例題２】 アルファベット26文字を使って、長さ32以下の単語を作ったときにアルファベットをすべて含む単語の数を10^9+7で割った余りを求めよ。

これも上の関数**nSurj**を使って、余りをとれば求まります。

```python
answer = 0
a, N = 26, 32
for n in range(a, N+1):
  ans = nSurj(n, a) % M
  print(f"n = {n}, a = {a}, ans = {ans}")
  answer += ans
  answer %= M
print(f"Answer: {answer}")
# n = 26, a = 26, ans = 459042011
# n = 27, a = 26, ans = 123744734
# n = 28, a = 26, ans = 476945811
# n = 29, a = 26, ans = 61156661
# n = 30, a = 26, ans = 46804531
# n = 31, a = 26, ans = 256387855
# n = 32, a = 26, ans = 993916781
Answer: 417998370
```


### 【例題３】100文字を使って、長さ10000以下の単語を作ったときに100文字をすべて含む単語の数を10^7で割った余りを求めよ。

ここまで数が大きくなると、このプログラムでは苦しくなります。高速化を考えます。
このプログラムを数式で表すと

```math 
\begin{align}
&\sum_{n=a}^{N} \color{blue}{\mathrm{nSurj(n,a)}}\\
&= \sum_{n=a}^{N} \color{blue}{\sum_{i=0}^{k}(-1)^{k-i}\binom{k}{i}i^n} \\
&で表されます\\
&= \sum_{n=a}^{N} \color{red}{\sum_{i=0}^{k}(-1)^{k-i}\binom{k}{i}}i^n \\
&するとこの\color{red}{赤い部分}にはnがないので前に出せます\\
&= \color{red}{\sum_{i=0}^{k}(-1)^{k-i}\binom{k}{i}}\sum_{n=a}^{N} i^n \\
&この黒い部分は等比数列の和になるので公式が使えます\\
&= \color{red}{\sum_{i=0}^{k}(-1)^{k-i}\binom{k}{i}}\frac{i^{N+1}-i^a}{i-1} \cdots (1)\\ \\
&これでnでのループを回す必要がなくなります、\\ \\
&さらに二項係数の部分は漸化的に計算できます \\ \\
& \binom{k}{i+1} = \binom{k}{i} \times \frac{k-i}{i+1}\cdots (2) \\ \\
& また式(1)はi=1のとき分母が0となるので別途以下の式で求めます　\\ \\
& (-1)^{k-1}k(n-k+1) \cdots (3)
\end{align}
```
これをコードにしたのがこちらです。例題２の答えが同じであることを確認して、例題３も無事に解くことができました。

```python 
M = 10**9+7
for k, n in [(26, 32), (10**2, 10**5)]:
  bino = k * (k - 1) // 2  # binomial(k,2)
  answer = (-1)**(k-1) * k * (n - k + 1)  #  (3) for i == 1 
  for i in range(2, k + 1):  # Loop for i > 2
    answer += (-1)**(k - i) * bino * (i**(n + 1) -  i**k) // (i - 1)   # (1)
    bino = bino * (k - i) // (i + 1)  # (2)
  print(f"k = {k}, n = {n}, ans = {answer%M}")
# k = 26, n = 32, ans = 417998370
# k = 100, n = 100000, ans = 481055237
```

### 【例題４】10^4文字を使って、長さ10^10以下の単語を作ったときに100文字をすべて含む単語の数を10^7で割った余りを求めよ

まだやるのかという感じですが、Pythonの剰余計算ではべき乗$i^n$を**pow(i,n,M)** を使うことによって$O(log(n))$で行うことができます。

注意点としては剰余の割り算はそのままでは出来ないので **pow(x,-1,M)** を掛けるようにします。

これで無事、例題４も数秒で解答が得られました。

```python 
M = 10**9+7
for k, n in [(26, 32), (10**2, 10**5),(10**4, 10**10),]:
  bino = (k * (k - 1) // 2) % M  # binomial(k,2)
  answer = ((-1)**(k-1) * k * (n - k + 1)) % M #  (3) for i == 1
  for i in range(2, k + 1): # Loop for i > 2
    answer += (-1)**(k - i) * bino * (pow(i,n + 1,M) -  pow(i,k,M)) * pow(i - 1,-1,M)   # (1)
    answer %= M
    bino = bino * (k - i) * pow(i + 1,-1,M)  # (2)
    bino %= M
  print(f"k = {k}, n = {n}, ans = {answer}")
# k = 26, n = 32, ans = 417998370
# k = 100, n = 100000, ans = 481055237
# k = 10000, n = 10000000000, ans = 173510926
```

（開発環境：Google Colab）

この考え方はProject Euler Problem 657: Imcomplete Wordを解くのに役に立ちます。