技術メモ

役に立てる技術的な何か、時々自分用の覚書。幅広く色々なことに興味があります。

置換可能素数

以前の「Nが現れる素数」に続いて面白い素数が紹介されていた。*1

このような素数を求めてみた。

ルール

<ルール>
素数Pに対して下の条件を満たすn(1,2,...9)が存在する。
1. Pの各桁の中にnを含まない。
2. Pの各桁をnで置換した数も全て素数になる。

(タイトルでは置換可能素数と呼んでいるけれど、そういう言葉はないです。)

コード

python3.5.3で動きます

from sympy import isprime
from math import log10, floor

def replace(num, replace_num):
    digits = len(str(num))
    replaced_nums = [str(num) for i in range(digits)]
    for i, replaced_num in enumerate(replaced_nums):
        replaced_num = replaced_num[:i] + str(replace_num) + replaced_num[i+1:]
        replaced_nums[i] = replaced_num
    return replaced_nums

replace_list = [1,3,7,9]
for prime in range(10,100000):
    if not isprime(prime):
        continue
    for replace_num in replace_list:
        if str(replace_num) in str(prime):
            continue
        replaced_nums = replace(prime, replace_num)
        for replaced_num in replaced_nums:
            if not isprime(replaced_num):
                break
        else:
            print("%8d :(%d) "%(prime, replace_num))

結果

11 :(3)
11 :(7)
13 :(7)
17 :(3)
17 :(9)
19 :(7)
31 :(7)
37 :(1)
41 :(3)
41 :(7)
43 :(1)
43 :(7)
47 :(1)
47 :(3)
61 :(7)
67 :(1)
71 :(3)
73 :(1)
79 :(1)
107 :(3)
107 :(9)
109 :(7)
137 :(9)
139 :(7)
167 :(3)
179 :(3)
197 :(3)
251 :(7)
337 :(1)
347 :(9)
409 :(1)
439 :(1)
449 :(3)
499 :(1)
607 :(1)
643 :(7)
709 :(1)
859 :(3)
1013 :(9)
1019 :(3)
1061 :(3)
1433 :(9)
1499 :(3)
1613 :(9)
2089 :(3)
3637 :(1)
3709 :(1)
3767 :(9)
4337 :(9)
4877 :(1)
4933 :(7)
6997 :(1)
7487 :(1)
8623 :(9)
9433 :(1)
9433 :(7)
16901 :(3)
25633 :(9)
36493 :(7)
47717 :(3)
56269 :(3)
76963 :(1)
86269 :(3)

nCk (mod m) の求め方 [n,kが凄く大きい場合]

大きいnに対してコンビネーションを求める

凄く大きいnに対してCombinationを求める時、例えば {}_{1000} \mathrm{C}_{400}を求めようと思っても愚直に計算するにはプログラムで扱える桁数を超えてしまう。
そこで、mで割った余り(modular)  {}_n \mathrm{C}_k \pmod{m}を求めさせることがある。

pythonで書いたのでサンプルコードを載せておく。
ただし、これを使うにはmはnを超える素数でないといけない

サンプルコード

再帰を使っている。

def mod_combination(n, k, mod):
    # nCk (mod m)
    def mod_permutation(n, k, mod):
        if n<=k:
            return 1
        else:
            return (n * mod_permutation(n-1,k,mod))%mod

    def mod_inv_permutation(k, mod):
        k, mod = int(k), int(mod)
        if k<=1:
            return 1
        else:
            return (pow(k,mod-2,mod) * mod_inv_permutation(k-1, mod))%mod

    return (mod_permutation(n,n-k,mod) * mod_inv_permutation(k, mod))%mod

理屈

割り算に対してモジュラ(余り)を計算するにはどうするかという問題がある。
それにはモジュラ逆数を考える必要がある。
プログラムでは

pow(k, mod-2, mod)

にあたる部分でモジュラー逆数を使っている。

モジュラ逆数とはa^{-1}\equiv x \pmod{m}となるxのこと。

結論から言うと

素数 mと互いに素な aに対して
 a^{-1}\equiv a^{m-2} \pmod{m}

となることを使っている。

(証明)
オイラー関数 \phi(m)を用いて、mが素数の時
\displaystyle
\left\{
\begin{array}{1}
 a^{\phi(m)}\equiv 1 \pmod{m}\\
 \phi(m) = m-1 
\end{array}
\right .
なので

a^{\phi(m) - 1} \equiv a^{-1} \pmod{m} \\
a^{m-2} \equiv a^{-1} \pmod{m}
(Q.E.D)

したがって {}_n \mathrm{C}_k (n < m)に対して

\begin{eqnarray}
{}_n \mathrm{C}_k & =  & n(n-1)\dots(n-k) \times k^{-1}(k-1)^{-1}\dots 1^{-1} \\
& \equiv &  n(n-1)\dots(n-k)\times k^{m-2}(k-1)^{m-2}\dots 1^{m-2} \pmod{m}
\end{eqnarray}


上のプログラムでは

def mod_permutation(n, k, mod)

 n(n-1)\dots(n-k) \pmod{m}  を

def mod_inv_permutation(n, k, mod)

 k^{m-2}(k-1)^{m-2}\dots 1^{m-2} \pmod{m} を
計算している。

おまけ

素数を確認したい時pythonではsympyを使おう

from sympy import isprime
print(isprime(1000000007)) #True