RSA実装?

rsa.py

#-*- coding: utf-8 -*-

"""

RSA公開暗号方式の仕組みを理解するために素朴に実装してみる。もとのメッ
セージは、バイト列として受け取り、整数のリストに暗号化する。

RSAは、(大きな数になると)素数の積を計算するのは簡単だが、逆に、合成数
を素因数分解するのは難しいという事実を利用したアルゴリズムである。

まず、ある人が、大きな素数を使って掛け合わせた合成数を自分専用の暗号化
をするための公開鍵として公表する。他の人は、その公開鍵を使って、送りた
いメッセージを暗号化する。この暗号化は、公開鍵の素因数分解(秘密鍵)を使
わない限り解けないようになっている。これによって、秘密鍵を知っている人
だけが、暗号化されたメッセージを見られることが保証される。


----

大雑把な流れは以下の通り。

秘密鍵を(p,q)とし、公開鍵を(k,m)とする。秘密鍵(p,q)には、二つの大きな
素数(100桁以上?)を選ぶ。公開鍵(k,m)は、m=pqとし、kはmと互いに素な数を
取る。

暗号化(ciphering):

01. 暗号化するメッセージを整数aとする。
02. b ≡ a^k(mod m)を計算する。
03. 暗号化されたメッセージbを得る。

復号化(deciphering):

01. 暗号化されたメッセージを整数bとする。
02. 合同式 a^k ≡ b(mod m) を解く。
03. 複合化されたメッセージaを得る。

復号化の02において合同式を解く必要があるが、これは、mの素因数分解が分
かっていなければ(すなわち、秘密鍵(p,q)を知らなければ)解くことが難しい。
そして、p,qが大きくなれば、mを素因数分解することは困難である。このこと
によって、暗号の安全性が保証される。

手順からわかるように、メッセージは、法mより小さい整数でなければならな
い。もしaがm以上ならば、もとのメッセージは復元できなくなるだろう。

----

aとbがmを法として合同である(a≡b(mod m),b<m)とは、aをmで割った余りがb 
に等しいこと、すなわち、適当な整数qに対して、a-b=qmが成り立つことであ
る。

法の同じ合同式においては、積が成り立つ。すなわち、
a≡b(mod m),c≡d(mod m) ⇒ ac≡bd(mod m)、
このとき、a=c,b=dと置けば、a^2≡b^2(mod m)。

また、gcd(c,m)=1ならば、ac≡bc(mod m)⇔a≡b(mod m)。

つまり、積は常に成り立つが、割り算は、割る数cが法mと互いに素な場合にし
か成り立たない。

これらが成り立つことは、合同式を普通の方程式の形に直せば、簡単に確認で
きる。

----

オイラーのφ関数

φ(n)を1からnまでで、nと互いに素な数(共通因数を持たない数)の個数とする
と、nが素数ならばφ(n)=n-1である。また、m,nを互いに素な数とすれば、
φ(mn)=φ(m)φ(n)が成り立つ。

よって、p,qが共に素数ならば、φ(pq)=φ(p)φ(q)=(p-1)(q-1)。

----

オイラーの公式

gcd(a,m) = 1のとき、a^φ(m) ≡ 1(mod m) が成り立つ。

b_1,b_2,...,b_φ(m)を1からmまでにあるmと互いに素な数のすべてとする。こ
のとき、aもmと互いに素であることから、gcd(ab_i,m)=1である。

ab_1,ab_2,...,ab_φ(m) (mod m)は、それぞれ、b_1,b_2,...,b_φ(m)のうち
のいずれかひとつと合同である({b_i}によって、mと互いに素な数は網羅され
るから)。

また、ab_i≡ab_j(mod m)とすると、ab_i-ab_j=a(b_i-b_j)がmで割り切れる。
aはmと互いに素なので、b_i-b_jは、mで割り切れて、かつ、|b_i-b_j|<mを満
たす数である。よって、b_i-b_j=0、すなわちb_i=b_j。よって、
ab_1,ab_2,...,ab_φ(m)(mod m)は、互いに異なるφ(m)個の数である。

以上により、ab_1,ab_2,...,ab_φ(m)(mod m)は、b_1,b_2,...,b_φ(m)と順序
を除いて同じである。

これらを掛け合わせれば、
ab_1*ab_2*...*ab_φ(m) ＝ a^φ(m)(b_1*b_2*...*b_φ(m))
                       ≡ b_1*b_2*...*b_φ(m)(mod m)

となるが、gcd(m,b_1*b_2*...*b_φ(m))=1なので、a^φ(m)≡1(mod m)が成り
立つ。

"""

import types, random, struct

def is_integer(i):
    return isinstance(i,types.IntType) or \
           isinstance(i,types.LongType)
    
def gcd(a,b):
    """
    ユークリッドの互助法によって、aとbに最大公約数を求める。

    ユークリッドの互助法とは、任意の整数aが必ず、a=bq+r(b,q,rは整数で、
    a≧b>r≧0)の形に表せること(剰余の定理)を利用して、

    a       = b       * q_0      + r_0
    b       = r_0     * q_1      + r_1
    r_0     = r_1     * q_2      + r_2
    r_1     = r_2     * q_3      + r_3
    ...
    r_(n-2) = r_(n-1) * q_n      + r_n
    r_(n-1) = r_n     * q_(n+1)  + 0

    というふうに、次々に割り算していくことによって、r_iは次第に小さく
    なっていき、これをr_(n+1)が0になるまで繰り返すことができる。このと
    き、r_nが、a,bの最大公約数であることが保証されている。

    最後の結果、r_(n-1)=r_n*q_(n_1)を、代入しながら、計算過程を逆にた
    どっていけば、r_nがa,bの公約数であることが確認できる。また、a,bの
    任意の公約数はr_nを割り切るので、r_nは最大公約数であることが確認
    できる。
    
    """
    assert( is_integer(a) and is_integer(b) )

    if a == 0 or b == 0:
        return max(a,b)

    r = a%b
    while( r != 0 ):
        a = b
        b = r
        r = a%b
    return b

def solve_linear_equation(a,b):
    """
    一次方程式 ax + by = gcd(a,b)を満たす(x,y)をひとつ見つける。
    
    一次方程式 ax + by = gcd(a,b)は、必ず解を持つことが保証される。な
    ぜなら、ユークリッド互助法の計算過程によって、

    a       = b       * q_0      + r_0
    b       = r_0     * q_1      + r_1
    r_0     = r_1     * q_2      + r_2
    r_1     = r_2     * q_3      + r_3
    ...    
    r_(n-2) = r_(n-1) * q_n      + r_n
    r_(n-1) = r_n     * q_(n+1)  + 0

    となり、各式は r_(i+1) = r(i-1) - r_i * q_(i+1)と変形できる。これ
    を上から順番に 各式のr_(i+2)に代入していけば、最終的に、
    gcd(a,b)=r_nをax+byの形に表せる。

    """
    assert( is_integer(a) and is_integer(b) )
    assert( a != 0 or b != 0 )

    if (a == 0 and b != 0):
        return (0,1)
    if (a != 0 and b == 0):
        return (1,0)

    q,r = divmod(a,b)
    (x0,y0) = (1,0)
    (x1,y1) = (0,1)
    (x2,y2) = (x0-q*x1,y0-q*y1)
    while( r != 0 ):
        a = b
        b = r
        q,r = divmod(a,b)
        (x0,y0) = (x1,y1)
        (x1,y1) = (x2,y2)
        (x2,y2) = (x0-q*x1,y0-q*y1)
    return (x1,y1)

def exponentiation_by_squaring(a,k,m):
    """

    繰り返し自乗法によって、a^k≡b(mod m)となるb(<m)を求める。

    繰り返し自乗法は、以下のような数の性質を利用して、べき乗の形の合同
    計算を簡単に行うアルゴリズムである。

    - 二つ以上の合同式の間で積が成立すること。

    - a^kを指数法則によって、2のべき乗の指数を持つ、より小さい数の積に
      分割できること。指数法則を使えば、べきを積に、積を和に直せるので
      計算が簡単になる。
    
    - 合同式の二乗の表を作っていく作業は、高々m^2の大きさの数しか扱わ
    　ないで済むので計算が楽(繰り返し3乗法とかだと、表を作るときの計算
    　にm^3が入ってくるので2の方がいい？)。

    まず、kを二進数展開して、

    k = u_n*2^n + ... + u_1*2^1 + u_0*2^0 (u_i∈{0,1})

    を得る。次に、

    a^(2^0) ＝ (a^1)   ≡ A_0 (mod m)
    a^(2^1) ＝ (a^1)^2 ≡ (A_0)^2 ≡ A_1 (mod m)
    a^(2^2) ＝ (a^2)^2 ≡ (A_1)^2 ≡ A_2 (mod m)
    ...
    a^(2^n) ≡ (A_(n-1))^2 ≡ A_n (mod m)

    のように、a(mod m)を次々に2乗していき、a^(2^n)までのテーブルを作る。
    この方法では、毎回mで剰余を取るので、a^(2^i) ≡ A_iは、m^2以上の数
    にはならない。

    以上により、
    
    a^k ＝ a^(u_n*2^n + ... + u_1*2^1 + u_0*2^0)
        ＝ a^(u_n*2^n) * ... * a^(u_1*2^1) * a^(u_0*2^0)
        ≡ A_0 * A_1 * ... * A_n (mod m).

    これが成り立つことは、合同式の積が成り立つことを考えればわかる。
    
    """
    assert( is_integer(a) and is_integer(k) and is_integer(m) )
    assert( k >= 0 )
    assert( m != 0 )
    
    if a == 0:
        return 0

    b = 1
    A = a%m
    bit = 1L<<0
    while( bit < k ):
        if k&bit:
            b = (A*b)%m
        A = (A**2)%m
        bit = bit<<1
    return b%m

def solve_congruence(b,k,m,phi_m):
    """
    
    x^k ≡ b(mod m)となる、xを求める。

    まず、1次方程式 kx + φ(m)y = 1を解く。このとき、gcd(k,φ(m))=1ならば、
    この方程式は解を持つ。この解を(u,v)とすれば、ku + φ(m)v = 1より、

    (x^k)^u = x^(ku) = x^(1-φ(m)v) = x*(x^φ(m))^(-v)。

    オイラーの公式より、x^φ(m) ≡ 1(mod m)なので、

    (x^k)^u ≡ x ≡ b^u(mod m)。

    よって、b^u(mod m)を繰り返し自乗法によって計算すれば、xが求まる。

    uが負の解を得る場合もありうるので、そのときは、φ(m)を加えて調整する(方程式の
    一般解は、(u+tφ(m),v-tk)である)。
    
    """
    assert( is_integer(b) and is_integer(k) \
            and is_integer(m) and is_integer(phi_m) )
    assert( gcd(k,phi_m) == 1 )

    (u,v) = solve_linear_equation(k,phi_m)
    while( u < 0 ):
        u += phi_m
        v -= k
    assert( k*u + phi_m*v == 1 )
    return exponentiation_by_squaring(b,u,m)

def rabin_miller_test(p,n):
    """
    
    ラビン・ミラー判定法を用いて、pの素数判定をn個の異なる数を用いて行
    う。結果が真であるとき、pが合成数である確率は、0.25^nパーセントで
    ある。
    
    """
    assert( (p-1)%2 == 0 )
    assert( 2 < n+2 <= p-1 )
    k = 0
    q = p-1
    while( q%2 == 0 ):
        q /= 2
        k += 1
    for a in range(2,n+2):
        if exponentiation_by_squaring(a,q,p) == 1:
            continue
        e = q
        for i in range(k):
            if exponentiation_by_squaring(a,e,p) == p-1:
                break
            e *= 2
        else:
            return False
    return True

def get_prime(n):
    """
    
    nバイト長の素数を生成する。

    *memo

    いちおうランダムな素数を生成してはいるんだけど、分散してなさすぎな
    気がする。もうちょっと細かく乱数生成してから結合するとか、工夫した
    方がいいかな。

    """
    m = 1L<<(n*8)
    d = 1L<<(n*8-1)
    n = random.randint(m,m+d)
    if n%2 == 0:
        n += 1
    while(not(rabin_miller_test(n,100))):
        n += 2
    return n

def generate_keys():
    """
    
    秘密鍵(p,q)、公開鍵(k,m)を生成する。

    指数kの選び方は、gcd(k,φ(m))=1となる数ならなんでもいいが、3かフェ
    ルマー素数の65537がよく使われるらしい。3はφ(m)と互いに素になる最
    小の数だから(p,q)ともに奇数であることに注意すれば、φ(m)は偶数であ
    ることがわかる)で、65537は、2^16+1と変形でき、べき乗計算が楽だから。

    http://www.cybersyndrome.net/rsa/rsa3.html
    
    """
    p = get_prime(16)
    q = get_prime(16)
    m = p*q
    k = 65537
    return ((p,q),(k,m))

def encrypt(a,k,m):
    """

    メッセージaを公開鍵(k,m)を用いて暗号化する。

    """
    assert( a < m )
    return exponentiation_by_squaring(a,k,m)
    
def decrypt(b,k,m,p,q):
    """

    暗号bを公開鍵(k,m)と秘密鍵(p,q)を用いて復号化する。
    
    """
    phi_m = (p-1)*(q-1)
    assert( m == p*q )
    assert( gcd(k,phi_m) == 1 )
    return solve_congruence(b,k,m,phi_m)

def split(l_int,size):
    """
    
    整数l_intをsizeバイトごとに分割して、そのリストを得る。

    size=4なら、
    0x102030405060708090
     ⇒ (0x00000010,0x20304050,0x60708090)
     
    """
    assert( size != 0 )
    
    mask = 0L
    for i in range(size):
        mask |= 0xFFL<<(i*8)
    i_list = []
    i_list.append( l_int&mask )
    i = 1
    while( l_int > mask ):
        mask <<= size*8
        i_list.append((l_int&mask)>>(i*size*8))
        i += 1
    return i_list

def join(i_list,size):
    """
    
    整数のリストを結合して、ひとつの整数を得る
    
    """
    assert( size != 0 )

    l_int = 0L
    for i in range(len(i_list)):
        l_int |= i_list[i]<<(i*size*8)
    return l_int

def split_size(m):
    """

    法mから、データの分割サイズを決定する。

    mを越えない範囲で、なるべく大きくなるような1<<(8*k)を与えるkを選ぶ。
    (cf. split(),join())

    """
    k = 0
    while( 1L<<(8*(k+1)) < m ):
        k += 1
    return k

rsa_test.py

#! /bin/env python
#-*- coding: utf-8 -*-

"""

rsa.pyのテスト

"""

import rsa

def test_gcd():
    print "testing gcd():",
    safe = True
    if rsa.gcd(1,1) != 1: safe = False
    if rsa.gcd(1,0) != 1: safe = False
    if rsa.gcd(0,1) != 1: safe = False
    if rsa.gcd(0,0) != 0: safe = False
    if rsa.gcd(15,25) != 5: safe = False
    if rsa.gcd(24,30) != 6: safe = False
    if rsa.gcd(7,6673) != 1: safe = False

    if safe:
        print "OK."
    else:
        print "FAIL."

def test_solve_linear_equation():
    print "testing solve_linear_equation():",
    safe = True
    if rsa.solve_linear_equation(0,1) != (0,1):
        safe = False
    if rsa.solve_linear_equation(1,0) != (1,0):
        safe = False
    if rsa.solve_linear_equation(1,1) != (1,0) and \
       rsa.solve_linear_equation(1,1) != (0,1):
        safe = False
    if rsa.solve_linear_equation(12453,2347) != (304,-1613):
        safe = False
    if rsa.solve_linear_equation(22,60) != (11,-4):
        safe = False

    if safe:
        print "OK."
    else:
        print "FAIL."

def test_exponentiation_by_squaring():
    print "testing exponentiation_by_squaring():",

    safe = True
    if rsa.exponentiation_by_squaring(0,2,3) != 0:
        safe = False
    if rsa.exponentiation_by_squaring(2,0,3) != 1:
        safe = False
    if rsa.exponentiation_by_squaring(7,327,853) != 286:
        safe = False
    if rsa.exponentiation_by_squaring(2,283976710803262,283976710803263) \
           != 280196559097287:
        safe = False

    if safe:
        print "OK."
    else:
        print "FAIL."

def test_solve_congruence():
    print "testing solve_congruence():",

    safe = True
    if rsa.solve_congruence(758,131,1073,1008) != 905:
        safe = False
    if rsa.solve_congruence(34781,3968039,27040397,27029856) != 22929826:
        safe = False

    p,q = 2473,3673
    m,phi_m = p*q, (p-1)*(q-1)
    a,k = 324, 6367
    b = rsa.exponentiation_by_squaring(a,k,m)
    if rsa.solve_congruence(b,k,m,phi_m) != a:
        safe = False
        
    if safe:
        print "OK."
    else:
        print "FAIL."

def test_join_and_split():
    print "testing join() and split():"

    i_list = rsa.split(0x102030405060708090L,4)
    for i in i_list:
        print "%x,"%i,
    l_int = rsa.join(i_list,4)
    print "| %x"%l_int

    i_list = rsa.split(0xFFFFFFFFL,4)
    for i in i_list:
        print "%x,"%i,
    l_int = rsa.join(i_list,4)
    print "| %x"%l_int
    
    i_list = rsa.split(0x01FFFFFFFFL,4)
    for i in i_list:
        print "%x,"%i,
    l_int = rsa.join(i_list,4)
    print "| %x"%l_int

def test_rabin_miller_test():
    print "testing rabin_miller_test():"
    safe = True
    if rsa.rabin_miller_test(294409,100) != False:
        safe = False
    if rsa.rabin_miller_test(172947529,100) != False:
        safe = False
    if rsa.rabin_miller_test(7212610147295474909544523785043492409969382148186765460082500085393519556525921455588705423020751421,100) != True:
        safe = False
    if rsa.rabin_miller_test(2074722246773485207821695222107608587480996474721117292752992589912196684750549658310084416732550077,100) != True:
        safe = False
    p = 2425967623052370772757633156976982469681*1451730470513778492236629598992166035067
    if rsa.rabin_miller_test(p,100) != False:
        safe = False
        
    if safe:
        print "OK."
    else:
        print "FAIL."
    
def test_rsa_coding():
    print "testing encrypt() and decrypt():"
    
    ((p,q),(k,m)) = rsa.generate_keys()
    phi_m = (p-1)*(q-1)
    message = 0x123456789
    print "k=",k
    print "m=%x"%m
    print "phi_m=%x"%phi_m
    print "p=%x"%p
    print "q=%x"%q
    print "message=%x"%message

    code = rsa.encrypt(message,k,m)
    message = rsa.decrypt(code,k,m,p,q)
    print "encrypted:%x"%code
    print "decrypted:%x"%message
    
# test_gcd()
# test_solve_linear_equation()
# test_exponentiation_by_squaring()
# test_solve_congruence()
# test_rabin_miller_test()
# test_join_and_split()
test_rsa_coding()