勉強したメモ。
大きい素数でちゃんと動くRSAの実装の記事があんまりなかったので。

目的

256bitRSA暗号を最小構成で実装する。

RSAのやりかた

詳細はwikipediaでも見て頂いて、流れだけ書きます。

1. bit数を決める。(今回は256bit)
2. 大きな素数$p,q$を作る。(128bit。256の半分)
3. $n=pq$
4. $(p-1)(q-1)$より小さく、互いに素な数$e$を決める。
5. $ed\equiv 1 \pmod n$である$d$を求める。

これで準備ができました。暗号化は、

1. 平文$m$を作る($m<n$)
2. 暗号文$c$は、$m^e \pmod n$

復号は、$m=c^d \pmod n$でOKです。

そんなに難しい計算はないですが、ポイントは

• べき乗計算が重い
• 大きな素数$p,q$をどうやって作る？
• $e$をどうやって決める？
• $d$は？

べき乗計算が重い

RSAでは数百桁の数でべき乗するので、素直に指数回掛け算などやってられません。
高速指数演算というのを使うと速くなります。
さらに、求めたいのは冪の剰余であって冪そのものではないので、この点でも最適化できます。
$a$の$b$乗を$n$で割ったあまりを計算するとしましょう。まず$b$を２進数にします。２進数の各桁を$b_i$として、
$$
a^b = a^{\sum 2^{i}b_i}=\prod a^{2^ib_i}
$$
となります。
ここで、
$$
a^{2^{i+1}}=a^{2^i}*a^{2^i}
$$
なので、掛け算一回で次の桁を計算できます。
さらに、計算したいのはべき乗の値を$n$で割ったあまりなので、掛け算の各ステップで随時あまりを取ってあげれば、$n$より小さい値を扱うだけで済みます。
例として、$a=3,b=11,n=21$とすると、まず$b=(1011)_2$つまり、
$$
b=2^3+2^1+2^0
$$

として、
$$
a^b=3^{2^3+2^1+2^0}=3^{2^3}*3^{2^1}*3^{2^0}
$$

$$
3^{2^0}=3
$$

$$
3^{2^1}=3^{2^0}*3^{2^0}=9
$$

$$
3^{2^2}=3^{2^1}*3^{2^1}=81
$$

$$
3^{2^3}=3^{2^2}*3^{2^2}=6561
$$

それぞれ$21$で割ったあまりにして、計算中も21を超えたらあまりを取ります。
$$
3*9*6561\equiv (3*9)*9\equiv 6*9\equiv 12
$$

def modular_exp(a, b, n):
    res = 1
    while b != 0:
        if b & 1 != 0:
            res = (res * a) % n
        a = (a * a) % n
        b = b >> 1

    return res

大きな素数をどうやって作るか

どうやら直接素数を生成するのは結構難しいようです。
実用的には、運良く素数に当たるまで乱数を生成するというのが簡単らしい。

乱数生成

bit数が指定されるので、その数だけ0,1をランダムに並べればOK。

def gen_rand(bit_length):
    bits = [random.randint(0,1) for _ in range(bit_length - 2)]
    ret = 1
    for b in bits:
        ret = ret * 2 + int(b)
    return ret * 2 + 1

最大ビットと最小ビットはなからず１になるようにしました。
最大ビットは桁数をが減らないように、最小ビットは奇数にするためです。
なぜ奇数かと言えば、今は素数がほしくて、偶数は素数ではないからです。

素数判定

巨大な乱数$p$が素数かどうか判定したい。
$p$はとても大きいので、順に割っていくような方法ではおそすぎる。
ミラー–ラビン素数判定法というのがよく使われるらしい。アルゴリズムはwikipediaにも乗ってるのでその通りに実装。
簡単に言えば、素数の必要条件を何回もチェックして、全部通ったらたぶん素数、一つでも通らなかったら合成数だと判定してます。
これは確率的素数判定のアルゴリズムなので、合成数を素数と判定してしまうことがあります。
ただ、その誤検知は１回のチェックで$\frac{1}{4}$程度だそうなので、何十回かやれば実用上ほぼ0になるので問題ないということらしいです。

def mr_primary_test(n, k=100):
    if n == 1:
        return False
    if n == 2:
        return True
    if n % 2 == 0:
        return False
    d = n - 1
    s = 0
    while d % 2 != 0:
        d /= 2
        s += 1

    r = [random.randint(1, n - 1) for _ in range(k)]
    for a in r:
        if modular_exp(a, d, n) != 1:
            pl = [(2 ** rr) * d for rr in range(s)]
            flg = True
            for p in pl:
                if modular_exp(a, p, n) == 1:
                    flg = False
                    break
            if flg:
                return False
    return True

ということで、素数生成はこんな感じです。

def gen_prime(bit):
    while True:
        ret = gen_rand(bit)
        if mr_primary_test(ret):
            break
    return ret

eをどうやって決める？

$(p-1)(q-1)$と互いに素ならなんでもいいので、素数にするのが手早いです。

p = gen_prime(128)
q = gen_prime(128)
e = gen_prime(128)

とかで問題ないのではないかと思います。
wikipediaには$65537$がよく使われると書いてました。決め打ち定数でもいいみたいですね。

dは？

拡張ユークリッドの互除法というやつが使えます。
もともと自然数$a,b$に対して、
$$
ax+by=gcd(a,b)
$$

となるような$x,y$を求めるアルゴリズムです。(gcdは最大公約数)
これを、$e,(p-1)(q-1)$に対して使います。最大公約数は$e$の定義から1なので、
$$
ex+(p-1)(q-1)y=1
$$

両辺を$(p-1)(q-1)$で割ると、
$$
ex\equiv 1
$$
となるので、この$x$が$d$です。
拡張ユークリッドの互除法はここにあったのでコピペしました。$x_0$が$d$です。

def xgcd(b, n):
    x0, x1, y0, y1 = 1, 0, 0, 1
    while n != 0:
        q, b, n = b // n, n, b % n
        x0, x1 = x1, x0 - q * x1
        y0, y1 = y1, y0 - q * y1
    return b, x0, y0


def gen_d(e, l):
    _, x, _ = xgcd(e, l)
    return x % l

実行

あとは平文を用意すれば動きます。

bit_length = 128
p = gen_prime(bit_length)
q = gen_prime(bit_length)
e = gen_prime(bit_length)
d = gen_d(e, (p - 1) * (q - 1))
n = p * q

m = 123456789
c = modular_exp(m, e, n)  # 暗号文
m_ = modular_exp(c, d, n)  # 123456789

1024bitでも一応動きました。(20秒くらいかかった)
実際は、同じ平文を何度も暗号化するのは危険なので、末尾に乱数を追加したり（復号時に削除）、mが小さすぎたらパディングしたりするそうです。