エルガマル暗号?

elgamal.py

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

"""

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

エルガマル暗号系は、g^k ≡ a(mod p)の形の合同式において、p,g,kからaを
計算するのは簡単だが、逆にp,gaから、kを決定するのは難しいこと(離散対数
問題)を利用した、公開鍵暗号系である。

まず、ある人が、大きな素数pと、それから計算した原始根gおよび、秘密鍵k 
を使って、g^k ≡ a(mod p)を計算し、これら法p,原始根g,及び公開鍵aを公表
する。他の人は、その公開鍵を使って、送りたいメッセージを暗号化する。こ
の暗号化は、公開鍵を計算した指数kがない限り解けないようになっている。
これによって、秘密鍵を知っている人だけが、暗号化されたメッセージを見ら
れることが保証される。

I(a) = k(g^k ≡ a(mod p))

底gにpの原始根というある数を選ぶと、I(a)が全単射になる。これの性質を暗
号に利用する。

----

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

秘密鍵をkとし、公開鍵をaとする。まず、暗号化／復号化の際に法として用い
る大きな素数p(150桁以上?)、および、pの原始根のひとつgを選ぶ。そして、
秘密鍵として適当な数kを選び、それらを用いて、公開鍵a≡g^k(mod p)を決定
する。法pと原始根gは、暗号化計算に必要なので、aといっしょに公開鍵とし
て公表しておく。

暗号化(ciphering):

01. 暗号化するメッセージを整数mとする。
02. 指数rをランダムに選ぶ。
03. e_1≡g^r(mod p),e_2≡m*a^r(mod p)を計算する。
04. 暗号化されたメッセージ(e_1,e_2)を得る。

復号化(deciphering):

01. 暗号化されたメッセージを整数(e_1,e_2)とする。
02. m ≡ (e_1^k)^(-1)*e_2(mod p)を計算する。
03. 複合化されたメッセージmを得る。

gが原始根であることにより、a ≡ g^k(mod p)において、各a,kは、1≦n≦p-1
の範囲で、一対一対応する。

式からわかるように、暗号化のための公開鍵aを計算するのは、繰返し二乗法
を一度実行すれば良いだけなので簡単に行えるが、逆に、aに対応する指数kを
計算するには、法pと原始根gから、aに対応するkが何なのかを決定しなければ
ならない。例えば、すべてのk(1≦k≦p-1)について、aが見つかるまで総当り
で計算するなどの方法が考えられるが、pが大きな素数なので現実的には不可
能となる。

----

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 get_prime(n):
    """
    
    nビット長の素数を生成する。

    ほんとは自分で素数テストをして、生成しないと駄目なんだけど、素数テ
    ストについてはまだ知らないので、とりあえず擬似コード。

    """
    return random.choice((
        2074722246773485207821695222107608587480996474721117292752992589912196684750549658310084416732550077,
        2367495770217142995264827948666809233066409497699870112003149352380375124855230068487109373226251983,
        1814159566819970307982681716822107016038920170504391457462563485198126916735167260215619523429714031,
        5371393606024775251256550436773565977406724269152942136415762782810562554131599074907426010737503501,
        6513516734600035718300327211250928237178281758494417357560086828416863929270451437126021949850746381,
        5628290459057877291809182450381238927697314822133923421169378062922140081498734424133112032854812293,
        2908511952812557872434704820397229928450530253990158990550731991011846571635621025786879881561814989,
        2193992993218604310884461864618001945131790925282531768679169054389241527895222169476723691605898517,
        5202642720986189087034837832337828472969800910926501361967872059486045713145450116712488685004691423,
        7212610147295474909544523785043492409969382148186765460082500085393519556525921455588705423020751421))
    
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 solve_congruence(a,b,m):
    """
    
    ax ≡ b(mod m)となる、xを求める。

    説明、、、書くのか？
    
    """
    assert( is_integer(a) and is_integer(b) and is_integer(m))
    assert( m != 0 )
    
    g = gcd(a,m)
    if( b%g != 0 ):
        return None
    else:
        (u,v) = solve_linear_equation(a,m)
        return (b*u/g)%m

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 <<= 1
    return b%m
    
def primitive_root(p):
    """
    
    素数pの最小の原始根を求める

    原始根定理により、素数pについて原始根の存在は保証されている。
    (が、pが素数でないと、、、)

    また、原始根の個数はφ(p-1)個なので、pが大きければ大きいほど、原始
    根をたくさんもつことがわかる。
    
    """
    def e(a,p):
        "a^e≡1(mod p)となる最小のeを求める"
        # 指数の整除性定理があるので、
        # 調べるのはp/2まででOK。
        i = 1
        ph = p/2
        while( i <= ph ):
            if exponentiation_by_squaring(a,i,p) == 1:
                return i
            i += 1
        #else:
        return p-1

    pr = None
    i = 2
    while( i < p ):
        if e(i,p) == p-1:
            pr = i
            break
        i += 1
    return pr

def generate_keys():
    """
    
    秘密鍵k、公開鍵(a,g,p)を生成する。
    
    """
    p = get_prime(None)
    p = 48112959837082048697L
    p = 5915587277L
    p = 1500450271L
    p = 380803L
    g = primitive_root(p)
    assert( g != None )
    #k = random.randint(1,p-1)
    k = 278374
    a = exponentiation_by_squaring(g,k,p)

    return (k,(a,g,p))

def encrypt(m,a,g,p):
    """

    メッセージmを公開鍵(p,g,a)を用いて暗号(e_1,e_2)を得る。

    rは1≦r≦p-1の範囲で、ランダムであれば、なんでもいい。
    
    e_1 = g^r (mod p)
    e_2 = ma^r (mod p).
    
    """
    assert( p != 0 , p)
    assert( 0 <= m <= p-1, m )
    assert( 1 <= g <= p-1, g )
    assert( 1 <= a <= p-1, a )
    
    r = random.randint(1,p-1)
    e_1 = exponentiation_by_squaring(g,r,p)
    e_2 = (m*exponentiation_by_squaring(a,r,p))%p
    return (e_1,e_2)
    
def decrypt(e_1,e_2,k,a,g,p):
    """

    暗号(e_1,e_2)を公開鍵(p,g,a)と秘密鍵kを用いてメッセージmを得る。

    c = e_1^k (mod p)
    u = c^(-1) (mod p)
    u*e_2 ＝ (e_1^k)^(-1)*e_2
          ≡ ((g^r)^k)^(-1)*ma^r
          ≡ ((g^k)^r)^(-1)*ma^r
          ≡ a^(-r)*a^r*m
          ≡ m (mod p).
          
    """
    assert( p != 0, p )
    assert( 1 <= e_1 <= p-1, e_1 )
    assert( 1 <= e_2 <= p-1, e_2 )
    assert( 1 <= g <= p-1, g )
    assert( 1 <= a <= p-1, a )
    assert( 1 <= k <= p-1, k )

    c = exponentiation_by_squaring(e_1,k,p)
    u = solve_congruence(c,1,p)
    assert( (c*u)%p == 1 )
    return (u*e_2)%p

elgamal_test.py

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

"""

elgamal.pyのテスト

"""

import elgamal

def test_gcd():
    print "testing gcd():",
    safe = True
    if elgamal.gcd(1,1) != 1: safe = False
    if elgamal.gcd(1,0) != 1: safe = False
    if elgamal.gcd(0,1) != 1: safe = False
    if elgamal.gcd(0,0) != 0: safe = False
    if elgamal.gcd(15,25) != 5: safe = False
    if elgamal.gcd(24,30) != 6: safe = False
    if elgamal.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 elgamal.solve_linear_equation(0,1) != (0,1):
        safe = False
    if elgamal.solve_linear_equation(1,0) != (1,0):
        safe = False
    if elgamal.solve_linear_equation(1,1) != (1,0) and \
       elgamal.solve_linear_equation(1,1) != (0,1):
        safe = False
    if elgamal.solve_linear_equation(12453,2347) != (304,-1613):
        safe = False
    if elgamal.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 elgamal.exponentiation_by_squaring(0,2,3) != 0:
        safe = False
    if elgamal.exponentiation_by_squaring(2,0,3) != 1:
        safe = False
    if elgamal.exponentiation_by_squaring(7,327,853) != 286:
        safe = False
    if elgamal.exponentiation_by_squaring(2,283976710803262,283976710803263) \
           != 280196559097287:
        safe = False

    if safe:
        print "OK."
    else:
        print "FAIL."
        
def test_solve_congruence():
    print "testing solve_congruence():",

    def t(a,b,m):
        x = elgamal.solve_congruence(a,b,m)
        if (a*x)%m == b:
            return True
        else:
            return False

    safe = True
    safe = t(20,15,19)
    safe = t(18,8,14)
    safe = t(893,266,2432)
        
    if safe:
        print "OK."
    else:
        print "FAIL."

def test_primitive_root():
    print "testing primitive_root():",

    safe = True
    if 2 != elgamal.primitive_root(5):
        safe = False
    if 3 != elgamal.primitive_root(7):
        safe = False
    if 2 != elgamal.primitive_root(11):
        safe = False

    p = 6221
    pr = elgamal.primitive_root(p)
    for i in range(1,p-1):
        if (pr**i)%p == 1:
            safe = False
        
    if safe:
        print "OK."
    else:
        print "FAIL."
        
    
def test_elgamal_crypting():
    print "testing encrypt() and decrypt():"
    
    (k,(a,g,p)) = elgamal.generate_keys()
    message = 123456
    print "k=",k
    print "a=",a
    print "g=",g
    print "p=",p
    print "message=",message

    e_1,e_2 = elgamal.encrypt(message,a,g,p)
    message = elgamal.decrypt(e_1,e_2,k,a,g,p)
    print "encrypted:",e_1,e_2
    print "decrypted:",message

    e_1,e_2 = 61745,206881
    print elgamal.decrypt(e_1,e_2,k,a,g,p)
    e_1,e_2 = 255836,314674
    print elgamal.decrypt(e_1,e_2,k,a,g,p)
    e_1,e_2 = 108147,350768
    print elgamal.decrypt(e_1,e_2,k,a,g,p)


    
# test_gcd()
# test_solve_linear_equation()
# test_exponentiation_by_squaring()
# test_solve_congruence()
# test_primitive_root()

test_elgamal_crypting()