numpyで2次元ベクトルの外積

一般的には2次元の外積の存在は賛否が分かれる所ですが、numpyでは2次元ベクトルの外積が用意されています。

import numpy

numpy.cross([0,1],[ 1,0]) #-1
numpy.cross([0,1],[-1,0]) # 1
numpy.cross([1,0],[ 0,1]) # 0

どうやら内積の90度回転版みたいですね。
キャラクタの方向から見て、点が右にあるのか左にあるのかの判断に使えそうです。

ローカル変数で「L」とか「D」の変数名を使うやつ!

ローカル変数だからって1文字の変数名を付ける奴!
本当にそれでいいと思ってんのか?!


特に大文字の1文字変数を名づけた上に、値を上書きするやつ!!
Rubyだったら怒られてるぞっ!!


このダメプログラマがっ!!
プログラミングの神様に謝れ!!!




ごめんなさい。



ローカル変数で「L」や「D」を使うのはダメプログラマ - kなんとかの日記

高速数値計算ライブラリ「Numpy」覚書き

Pythonで一番有名で普及しているライブラリと言っても過言ではない「Numpy」の覚書きです。かなり多機能な数値計算ライブラリで、内部はC言語で記述されているため超高速に動作します。

ベクトル

ベクトルの長さ&正規化
import numpy
a = numpy.array([[2,2]])

#ベクトルの長さ
length = numpy.linalg.norm(a)
#length=>2.8284271247461903

#ベクトルの正規化
a / numpy.linalg.norm(a)
#=>array([[ 0.70710678,  0.70710678]])
内積外積
import numpy
v1 = numpy.array((1,0,0))
v2 = numpy.array((0,1,0))

#内積
numpy.dot(v1,v2)    #=> 0
#外積
numpy.cross(v1,v2)  #=>[0 0 1] 

一般的には2次元の外積の存在は賛否が分かれる所ですが、numpyでは2次元ベクトルの外積が用意されています。

import numpy

numpy.cross([0,1],[ 1,0]) #-1
numpy.cross([0,1],[-1,0]) # 1
numpy.cross([1,0],[ 0,1]) # 0

内積を90度回転させた感じですね、左右の正と負の領域判定に使えそうです。

行列

単位行列の定義
import numpy
I = numpy.matrix(numpy.identity(3))
#[[ 1.  0.  0.]
# [ 0.  1.  0.]
# [ 0.  0.  1.]]
転置
import numpy
a = numpy.matrix([[1,2],[3,4]])
# [[1, 2],
#  [3, 4]]

a.T
# [[1, 3],
#  [2, 4]]
逆行列行列式
import numpy
a = numpy.matrix([[3,1], [1,2]])

#逆行列
a.I
#[[ 0.4, -0.2],
# [-0.2,  0.6]]

#行列式
numpy.linalg.det(a)
# 5.0
回転行列
import python
from numpy import sin,cos
r = numpy.pi/2.

#2次元回転
numpy.matrix( (
    ( cos(r), sin(r)),
    (-sin(r), cos(r))
) )

#3次元z軸回転
numpy.matrix( (
    ( cos(r), sin(r), 0.),
    (-sin(r), cos(r), 0.),
    (     0.,     0., 1.)
) )

#3次元x軸回転
numpy.matrix( (
    ( 1.,     0.,     0.),
    ( 0., cos(r), sin(r)),
    ( 0.,-sin(r), cos(r))
) )

#3次元y軸回転
numpy.matrix( (
    ( cos(r), 0.,-sin(r)),
    (     0., 1.,     0.),
    ( sin(r), 0., cos(r))
) )
行列とベクトルの積
import numpy
#行列
mat = numpy.matrix( (
    (2.0, 0.0, 10),
    (0.0, 3.0, 10),
) )
#ベクトル
vec = numpy.array( (2.0,2.0, 1.0) )

#行列とベクトルの積
numpy.dot(mat,vec)
#[ 14.  13.]
固有値固有ベクトル&対角化

固有値固有ベクトルは以下の式が成り立つベクトルとスカラ。

A*v = λ*v(λ=Aの固有値,v=Aの固有ベクトルとする)

行列の対角化(diagonal)とは、正方行列Aに対して対角化された行列をPとすると、

P^-1 * A * P = 対角要素が固有値の行列
(P^-1はPの逆行列を表す)

となる行列である。
なお、numpy.diagonal()は対角化する関数ではなく、単に行列の対角要素をベクトルとして返す関数なので注意。

import numpy
A = numpy.array([[5,-2],[-4,7]])
(l,v) = numpy.linalg.eig( A )
#l=>[ 3.,  9.]
#v=>[[-0.70710678,  0.4472136 ],
#    [-0.70710678, -0.89442719]]

その他

平均値&標準偏差
import numpy
pop = [79,81,77,78,83,80,82,78,80,82]

#平均
numpy.mean(pop)  #=>80.0
#標準偏差
numpy.std(pop)   #=>1.8973665961010275
相関係数相関係数行列
import numpy
a = [1,2,3,4,5]
b = [1,1,3,5,5]
c = [5,4,3,1,1]
d = [-1,3,-4,5,-5]

#aとbの相関係数
numpy.corrcoef(a,b)[0][1] #=> 0.948683298051
#aとcの相関係数
numpy.corrcoef(a,c)[0][1] #=> -0.972271824132
#aとdの相関係数
numpy.corrcoef(a,d)[0][1] #=> -0.218797487247

#a,b,c,dの相関係数行列
numpy.corrcoef([a,b,c,d])
#[[ 1.          0.9486833  -0.97227182 -0.21879749]
# [ 0.9486833   1.         -0.97827974 -0.1153164 ]
# [-0.97227182 -0.97827974  1.          0.01933915]
# [-0.21879749 -0.1153164   0.01933915  1.        ]]
高速フーリエ変換(FFT)
import numpy
t = numpy.arange(0.0, 2*numpy.pi, 2*numpy.pi/32)
c = numpy.cos(t)
C = numpy.fft.fft(c)

PythonでA*(A-Star)アルゴリズム

今回はA*アルゴリズムPythonでやってみます。
ゲームプログラマの間では、もはや常識となりつつある最短経路問題解決アルゴリズムです。


A*は、古典的手法であるダイクストラ法」を改良したものです。
スタート地点からノードnを通ってゴールに辿り付くとき、最短距離をf(n)とすると、


f(n) = g(n) + h(n)


とすることができます、g(n)は「スタートからノードnまでの最短距離」、h(n)は「ノードnからゴールまでの最短距離」です。
でも、最初から適切なg(n)とh(n)が判ってるなら苦労しませんよね。


だから、テキトーな予測値を使って、最短経路をある程度予測して効率的に経路探索をしてみようという事です。
テキトーな予測値を使った最短経路距離をf*(n)とすると


f*(n) = g*(n) + h*(n)


となります、f*(n)を求めるためにテキトーなg*(n)とh*(n)を求める訳です。
g*(n)は「スタートからノードnまでの最短距離予測」です、今まで辿ってきた経路を覚えていれば何とかなりそうです。
問題はh*(n)で「ノードnからゴールまでの最短距離予測」です。


A*が最短距離を保証するためにはh*(n)は、「0<= h*(n) <= h(n)」の範囲である必要があります。
つまり予測ではない本物のh(n)よりも、小さい値である必要があります。
ということですので、h*(n)には「直線距離(ユークリッド距離)」を使う場合が多いです。


余談ですがh*(n)が常に0の場合、「ダイクストラ法」と等しくなります、つまり、A*はh*(n)によってテキトーだがそこそこ近いヒント(これをヒューリスティック関数と呼ぶ)を与える事により、効率的に経路検索をしているわけです。


小難しい説明よりも、実際に動くものを見た方が早そうですね。
以下、PythonでA*を実装してみました。

#!/usr/local/bin/python
# -*- coding: utf-8 -*-
map_data = [
'OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO',
'OS  O     O     O         O          O',
'O   O  O  O  O  O         O    OOOO GO',
'O      O     O  O   OOOO  O    O  OOOO',
'OO OOOOOOOOOOOOOOO  O     O    O     O',
'O                O  O     O          O',
'O        OOO     O  O     OOOOOOOOO  O',
'O  OO    O    OOOO  O     O      OO  O',
'O   O    O          O     O  O   O   O',
'O   OOO  O          O        O   O   O',
'O        O          O        O       O',
'OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO',
]
map_width  = max([len(x) for x in map_data])
map_height = len(map_data)

class Node(object):
    """
    f(n) ... startからgoalまでの最短距離
    g(n) ... startからnノードまでの最短距離
    h(n) ... nノードからgoalまでの最短距離
    f(n) = g(n) + h(n)

    関数を推定値にすることにより最短距離を予測する
    h*(n)をnからgoalまでの直線距離と仮定する。

    f*(n) = g*(n) + h*(n)
    """
    start = None    #start位置(x,y)
    goal = None     #goal位置(x,y)

    def __init__(self,x,y):
        self.pos    = (x,y)
        self.hs     = (x-self.goal[0])**2 + (y-self.goal[1])**2
        self.fs     = 0
        self.owner_list  = None
        self.parent_node = None

    def isGoal(self):
        return self.goal == self.pos

class NodeList(list):
    def find(self, x,y):
        l = [t for t in self if t.pos==(x,y)]
        return l[0] if l != [] else None
    def remove(self,node):
        del self[self.index(node)]


#スタート位置とゴール位置を設定
for (i,x) in enumerate(map_data):
    if 'S' in x:
        Node.start = (x.index('S'),i)
    elif 'G' in x:
        Node.goal = (x.index('G'),i)

#OpenリストとCloseリストを設定
open_list     = NodeList()
close_list    = NodeList()
start_node    = Node(*Node.start)
start_node.fs = start_node.hs
open_list.append(start_node)

while True:
    #Openリストが空になったら解なし
    if open_list == []:
        print "There is no route until reaching a goal."
        exit(1);

    #Openリストからf*が最少のノードnを取得
    n = min(open_list, key=lambda x:x.fs)
    open_list.remove(n)
    close_list.append(n)

    #最小ノードがゴールだったら終了
    if n.isGoal():
        end_node = n
        break

    #f*() = g*() + h*() -> g*() = f*() - h*()
    n_gs = n.fs - n.hs

    #ノードnの移動可能方向のノードを調べる
    for v in ((1,0),(-1,0),(0,1),(0,-1)):
        x = n.pos[0] + v[0]
        y = n.pos[1] + v[1]

        #マップが範囲外または壁(O)の場合はcontinue
        if not (0 < y < map_height and
                0 < x < map_width and
                map_data[y][x] != 'O'):
            continue

        #移動先のノードがOpen,Closeのどちらのリストに
        #格納されているか、または新規ノードなのかを調べる
        m = open_list.find(x,y)
        dist = (n.pos[0]-x)**2 + (n.pos[1]-y)**2
        if m:
            #移動先のノードがOpenリストに格納されていた場合、
            #より小さいf*ならばノードmのf*を更新し、親を書き換え
            if m.fs > n_gs + m.hs + dist:
                m.fs = n_gs + m.hs + dist
                m.parent_node = n
        else:
            m = close_list.find(x,y)
            if m:
                #移動先のノードがCloseリストに格納されていた場合、
                #より小さいf*ならばノードmのf*を更新し、親を書き換え
                #かつ、Openリストに移動する
                if m.fs > n_gs + m.hs + dist:
                    m.fs = n_gs + m.hs + dist
                    m.parent_node = n
                    open_list.append(m)
                    close_list.remove(m)
            else:
                #新規ノードならばOpenリストにノードに追加
                m = Node(x,y)
                m.fs = n_gs + m.hs + dist
                m.parent_node = n
                open_list.append(m)


#endノードから親を辿っていくと、最短ルートを示す
n = end_node.parent_node
m = [[x for x in line] for line in map_data]

while True:
    if n.parent_node == None:
        break
    m[n.pos[1]][n.pos[0]] = '+'
    n = n.parent_node
print "\n".join(["".join(x) for x in m])

もうちょっと簡潔に書けるかな・・・?
(2009.7.15 NodeListクラスを作って、少しだけ簡潔にしました)
実行結果は以下の通りです、「S」がスタート位置、「G」がゴール位置、「O」が壁となっています。

OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
OS+ O     O     O         O   +++++++O
O + O  O  O  O  O  +++++++O   +OOOO GO
O +    O     O  O  +OOOO +O   +O  OOOO
OO+OOOOOOOOOOOOOOO +O    +O   +O     O
O ++++++++++++   O +O    +O   ++++++ O
O        OOO +   O +O    +OOOOOOOOO+ O
O  OO    O   +OOOO +O    +O +++++OO+ O
O   O    O   +++++++O    +O +O  +O++ O
O   OOO  O          O    ++++O  +O+  O
O        O          O        O  +++  O
OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO

結構難しい迷路ですけど、ちゃんとやってくれますね。
移動方向は上下左右の4方向のみですが、以下のコードを変更すると斜めを入れた8方向で探索してくれます。
(行が長ぇよ!という方は、移動方向を変数にて定義してください)

#ノードnの移動可能方向のノードを調べる
for v in ((1,0),(-1,0),(0,1),(0,-1)):
    ↓
#ノードnの移動可能方向のノードを調べる
for v in ((1,0),(-1,0),(0,1),(0,-1),(1,1),(-1,1),(1,-1),(-1,-1)):

実行結果は以下の通りです、斜めに進んでくれるので、前よりも効率的に進んでいるのがわかります。

OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
OS  O     O     O         O    +++++ O
O + O  O  O  O  O   +++++ O   +OOOO GO
O +    O     O  O  +OOOO +O   +O  OOOO
OO+OOOOOOOOOOOOOOO +O    +O   +O     O
O  ++++++++++    O +O    +O    ++++  O
O        OOO +   O +O    +OOOOOOOOO+ O
O  OO    O   +OOOO+ O    +O ++++ OO+ O
O   O    O    ++++  O    +O+ O  +O + O
O   OOO  O          O     +  O  +O+  O
O        O          O        O   +   O
OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO

実際にゲームに組み込む時は、もうちょっと工夫が必要かもしれませんね。

RSA暗号で「ふっかつのじゅもん」を作る(3)-完結

Pythonを使って高速素数判定をしてみる - Pashango’s Blog
RSA暗号で「ふっかつのじゅもん」を作る(1) - Pashango’s Blog
RSA暗号で「ふっかつのじゅもん」を作る(2) -RSA暗号鍵の生成


当初はサクッと終わらせるつもりだったのに、意外と長くなってしまいました。

ドラクエ1の「ふっかつのじゅもん」の仕様

まず、ドラクエ1の「ふっかつのじゅもん」の仕様を調べてみました。

  • 使われる文字は64種類のひらがなと数字
  • じゅもんの長さは20文字まで
  • 以下の情報を含む
    • なまえ(0〜63までの数値×4)
    • 魔法の鍵(0〜7までの数値)
    • 薬草(0〜7までの数値)
    • ぶき(0〜7までの数値)
    • ぼうぐ(0〜7までの数値)
    • たて(0〜3までの数値)
    • アイテム(0〜15までの数値×8)
    • フラグ(0〜1までの数値×5)
    • 経験値(0〜65535までの数値)
    • ゴールド(0〜65535までの数値)

呪文の長さは20文字ですので、情報量は6bit×20文字=120bitとなります。
そのうち、最上位bitはRSA暗号化時に失われる可能性がありますので、データ保持として使えるbit数は119bitとなります。
また、上記の情報を保存するのに必要なbit数は107bitですので、119-107=12bitの余裕bitがあります。

暗号鍵を求める(eとdを求める)

前回、暗号鍵の元となる素数p,qを求めましたね。
私もすっかり忘れていましたが、乗数であるeとdを求める関数を作りましょう。
前回のgenerate_pq()と一緒のファイルに書き足してください。

def generate_ed(p,q):
    """与えられた素数ペアからe,dを生成"""
    if p*q < 5: raise ValueError,"p*q must be large than 4."
    l = lcm(p-1,q-1)
    e = random.randint(5,l-1) | 1
    while True:
        e = e+2 if e+2 <= l-1 else 5
        if gcd(e,l) == 1:
            break
    
    d = gcd2(e,l)[0]
    if d < 0: d += l
    return (e,d)

さっそく、鍵を生成してみましょう。
ふっかつのじゅもんに必要なbit数は120bitですので、以下のように求めます。

p,q = generate_pq(120)
e,d = generate_ed(p,q)
n   = p*q

# p = 96568316136941074427
# q = 7372545182284901
# n = 711954273896770180320412909759326727
# e = 26578325672224134549282406558845001
# d = 34518827246730708105140436249657601

pとqはn,e,dが求まれば不要ですので、廃棄してください。
「公開鍵」がeとn、「秘密鍵」がdとnの数値ペアとなります。

勇者の情報を生成する

暗号化する勇者の情報を生成するクラスを作成しました。
このクラスはパラメータの乱数生成に加え、パラメータ<->107bit数値の相互変換機能を持っています。
以下、コードですが特に面白味もないコードなので、サラッと読み飛ばしてください。

#coding:utf8
import random

class DQParams(object):
    NAMES_CHARS = u"*あいうえおかきくけこさしすせそたちつてと"\
                  u"なにぬねのはひふへほまみむめもやゆよ"\
                  u"らりるれろわをんっゃゅょ ー"\
                  u"0123456789S"

    PARAM_INFO  = [
        {'name':u"なまえ",'max':len(NAMES_CHARS),'len':4},
        {'name':u"魔法の鍵",'max':8,'len':1},
        {'name':u"薬草",'max':8,'len':1},
        {'name':u"ぶき",'max':8,'len':1},
        {'name':u"ぼうぐ",'max':8,'len':1},
        {'name':u"たて",'max':4,'len':1},
        {'name':u"アイテム",'max':16,'len':8},
        {'name':u"フラグ1",'max':2,'len':5},
        {'name':u"経験値",'max':65535+1,'len':1},
        {'name':u"ゴールド",'max':65535+1,'len':1}
    ]

    def __init__(self, val=None):
        self.params = {}
        for x in self.PARAM_INFO:
            self.params[x['name']] = [0] * x['len']

        if val == None:
            val = random.randint(0,1<<120)

        dec = lambda val,max:(val/max,val%max)
        for x in reversed(self.PARAM_INFO):
            for y in range(x['len']):
                val,self.params[x['name']][y] = dec(val, x['max'])
            self.params[x['name']].reverse()

    def name2val(self, name):
        return [self.NAMES_CHARS.index(x) for x in name]

    def val2name(self, val):
        return ''.join([self.NAMES_CHARS[x] for x in val])

    def generate_value(self):
        val = 0
        for x in self.PARAM_INFO:
            for y in self.params[x['name']]:
                val = val*x['max'] + y
        return val

    def output(self):
        f = self.PARAM_INFO[0]
        print f['name'],u":",self.val2name(self.params[f['name']])
        for x in self.PARAM_INFO[1:]:
            print x['name'],u":",self.params[x['name']]


それでは、さっそく勇者を生成してみましょう。

obj = DQParams()
obj.output()
#なまえ : かんてつ
#魔法の鍵 : [3L]
#薬草 : [0L]
#ぶき : [2L]
#ぼうぐ : [1L]
#たて : [3L]
#アイテム : [4L, 12L, 12L, 1L, 4L, 0L, 0L, 14L]
#フラグ1 : [1L, 1L, 0L, 1L, 1L]
#経験値 : [46633L]
#ゴールド : [6040L]

val = obj.generate_value()
#val => 17045993125123514293179659130776

おお、よくきたな、ゆうしゃかんてつよ、
今後は、生成された勇者「かんてつ」のパラメータを使います。

チェックサムとランダム値

今回の仕様では、必要情報量107bitー最大情報量119bitで差し引き12bitの余裕bitがあります。
この12bitに「チェックサム」と「ランダム値」を埋め込みます。


チェックサム」は、この数値が正当あるかのチェックに使われ、適当に入力された呪文をはじく効果があります。
「ランダム値」は、データ値に乱数を埋め込むことにより、同じパラメータを暗号化しても毎回呪文が変わるようになり、暗号解析がしにくくなります。


今回の例では12bitの余裕bitのうち、チェックサムに9bit、ランダム値に3bitを割り当てます。
以下、出力された数値にチェックサムとランダム値を埋め込むコードです。

def check_sum(val,base=1<<9):
    csum = 0
    while val > 0:
        csum += val%(base-1)
        csum %= base
        val /= base
    csum = (base-1)-csum
    return csum

csum = check_sum(val)
#csum => 114

val = (val << 9) + csum
val = (val << 3) + random.randint(0,7)
#val => 69820387840505914544863883799659408


なお、チェックサムは((1<<9)-1)から引いた値になっていますが、これはvalが0にならないようにするためです。
valが0ですと暗号化した結果も0になってしまい、暗号になっていません。
チェックサムのビットを反転する事により、valに0が入る事を防いでいます。

ふっかつのじゅもん」を生成

ここまで長い道のりでしたが、ようやく「ふっかつのじゅもん」を生成します。

val = 69820387840505914544863883799659408
n   = 711954273896770180320412909759326727
e   = 26578325672224134549282406558845001
d   = 34518827246730708105140436249657601


val = pow(val,e,n)
#val => 692554535384224777977608304993455189

PASS_CHARS = u"あいうえおかきくけこさしすせそたちつてと"\
             u"なにぬねのはひふへほまみむめもやゆよわ"\
             u"らりるれろ"\
             u"がぎぐげござじずぜぞだぢづでどばびぶべぼ"
s = []
while True:
    s.append(PASS_CHARS[val%64])
    val /= 64
    if val <= 0: break
print "".join(s)
# にざぜくねだひかづそぎけあぬへさときぬめ


勇者「かんてつ」のふっかつのじゅもんは「にざぜくねだひかづそぎけあぬへさときぬめ」となりました。
さっそく、復活できるかチェックしてみましょう。

n   = 711954273896770180320412909759326727
d   = 34518827246730708105140436249657601
val = 692554535384224777977608304993455189

val = pow(val,d,n)
csum = (val>>3) & ((1<<9)-1)
val >>= 12
csum2 = check_sum(val,1<<9)
if csum != csum2:
    print u"ふっかつのじゅもんがちがいます"
else:
    obj = DQParams(val)
    obj.output()

#なまえ : かんてつ
#魔法の鍵 : [3L]
#薬草 : [0L]
#ぶき : [2L]
#ぼうぐ : [1L]
#たて : [3L]
#アイテム : [4L, 12L, 12L, 1L, 4L, 0L, 0L, 14L]
#フラグ1 : [1L, 1L, 0L, 1L, 1L]
#経験値 : [46633L]
#ゴールド : [6040L]


ゆうしゃかんてつよ、よくぞもどってきた

RSA暗号化 vs 旧来の暗号化

RSA暗号化による暗号化と、旧来の暗号化の違いはなんでしょうか?


まずは、保存可能な情報量の違いから。
旧来の暗号化はまるまる120bitに情報を詰め込めます、RSA暗号化は最上位ビットが不定になりますので119bitしか使えません。
余った情報はチェックサムに使われますので、RSA暗号化の方がチェックサムが1bit少なくなります、これはどういう事かと言いますと、適当に打ち込んだ呪文が通ってしまう確率が旧来暗号化に比べて2倍高いという事です。


一方、RSA暗号のメリットは、やはり暗号化と複合化に使われる鍵が2つにわかれている点でしょう。


旧来暗号化の弱点は、暗号化と複合化が必ずセット提供されてしまう点です。
「暗号化の手順=復号化の逆手順」であり、そのまた逆も然りです。
リバースエンジニアリングにより、暗号化の手順が解析されてしまうと、自動的に復号化の手順も解析されてしまうのです。

一方、RSA暗号では、暗号化と復号化のどちらか一方のみの提供が出来ます。
暗号化はサーバーで行い、クライアントは復号化しかできないという扱いも可能です。
RSA暗号でも、クライアントに暗号鍵と復号鍵を埋め込んだらリバースエンジニアリングには無力ですよ、念のため)


うんうん、なんか使えそうなんだけど、意外と使いどころが難しい技術ですね。
でも、工夫次第ですごい面白い事もできそうです。

RSA暗号で「ふっかつのじゅもん」を作る(2) -RSA暗号鍵の生成

Pythonを使って高速素数判定をしてみる - Pashango’s Blog
RSA暗号で「ふっかつのじゅもん」を作る(1) - Pashango’s Blog


さて前回からの続きです、RSA暗号で「ふっかつのじゅもん」を作ってみましょう。


RSA暗号では、暗号化するデータのビット数よりも、1ビット多いn(素数p,qの積)が必要です。
(nでmodするため、nより小さな数は復号ができないため)
ですので、任意ビットの大きさをもつnを生成する素数p,q(RSA暗号鍵)を生成してくれる関数があると非常に便利ですね。


WebではRSA暗号鍵の生成アルゴリズムの資料が無かったので、自分で考えました。もしかしたら鍵生成が甘いかもしれません、その場合は容赦なくツッコミをください。

素数判定「ミラー・ラビンテスト」

暗号鍵生成関数を作るには、「Pythonで高速素数判定」の回で作成した高速素数判定関数が必要です。
ミラー・ラビンテストを使った素数判定で、第1引数が素数ならばTrueを返します。

#coding:utf-8
import random

def is_prime(q,k=50):
    """ミラー・ラビンテストによる素数判定"""
    q = abs(q)
    if q == 2: return True
    if q < 2 or q&1 == 0: return False

    d = (q-1)>>1
    while d&1 == 0:
        d >>= 1
    
    for i in xrange(k):
        a = random.randint(1,q-1)
        t = d
        y = pow(a,t,q)
        while t != q-1 and y != 1 and y != q-1: 
            y = pow(y,2,q)
            t <<= 1
        if y != q-1 and t&1 == 0:
            return False
    return True

任意ビット数の素数生成

この素数判定関数を使って、任意ビットの素数を生成する関数を作りましょう。
ただし、2と3は素数として弱いのではじきますので、任意ビット数は2以上でなくてはなりません。

def generate_prime(bit):
    """任意ビット数の素数を生成する"""
    if bit < 3: raise ValueError,"bit must be large than 2."

    mini = 1 << (bit-1)
    maxi = (mini << 1)-1
    x = random.randint(mini,maxi)
    while True:
        if is_prime(x):
            break
        x = x+1 if x+1 <= maxi else mini
    return x

任意ビット数のnを生成する素数p,q

この素数生成関数を使って、「任意ビット数のn」を生成する素数p,qを返す関数を作ります。

def generate_pq(bit):
    """任意ビット数のnを生成する素数p,qを返す"""
    mini = 1 << (bit-1)
    maxi = (mini << 1)-1
    while True:
        t = random.randint(3,bit-3)
        p = generate_prime(t)
        q = generate_prime(bit-t)
        if mini < (p*q) < maxi:
            break
    return (p,q)


ふぅふぅ、結構バラつきのある素数の組み合わせを生成するのに苦労しました。
さっそく、この関数で遊んでみましょう。
ちなみにPython2.6から「format(x,'b')」で、整数型を2進数表記にできますので、これでビット数を確認します。

p,q = generate_pq(32)
print p*q
print "ビット数:",len(format(p*q,'b'))
# 2477487329
# ビット数: 32

p,q = generate_pq(128)
print p*q
print "ビット数:",len(format(p*q,'b'))
# 317586999734752741757043691807400509529
# ビット数: 128

p,q = generate_pq(257)
print p*q
print "ビット数:",len(format(p*q,'b'))
# 120029623918517010269581862828609021278006347146688949661886393751717303005779
# ビット数: 257


どうやら合ってるようですね、めでたしめでたし。
うう・・・、鍵の生成だけで時間が経ってしまいました、また続きます・・・・


RSA暗号で「ふっかつのじゅもん」を作る(3)-完結 - Pashango’s Blog

RSA暗号で「ふっかつのじゅもん」を作る(1)

オッス、オラ、トンヌラ
前回は、高速素数判定を作りましたが、今回はRSA暗号を使って、昔懐かしの「ふっかつのじゅもん」を作ってみましょう。


Pythonを使って高速素数判定をしてみる - Pashango’s Blog


あ、「今さらRSAかよ」と思いました?
自分でRSAを実装してみると、色々知らない事が出てきて面白いですよ。
あとRSAは、暗号化以外にも応用が利くんで覚えておいて損はしませんよ。

RSA暗号とはなにか?

ゲームプログラマは基本的にゲームばっかやってるんで、一般的な情報処理知識に欠けている場合がほとんどです。
まずはRSA暗号の説明から始めましょう。


RSA暗号とは、2つの鍵「公開鍵」と「秘密鍵」を使う暗号方式です。
「公開鍵」は暗号化キーです、みんなに公開してかまいません。
秘密鍵」は復号化キーです、みんなにバレてはいけません厳重に保管してください、間違ってもネット上には流してはいけません。


データを送ってもらう友達に「公開鍵」だけ渡し、公開鍵でデータを暗号化して送ってもらいます。
暗号化されたデータを復号化できるのは「秘密鍵」を持っているあなただけで、暗号化した友達も、データを横取りしようとしてるハッカーにも復号はできないのです。
このように暗号化キーだけを外部に出し、復号化キーは外部に出さない事で、強固なセキュリティを実現しています。

鍵を作ってみよう

それではRSA暗号の鍵を作ってみましょう。
まずは2つの素数p,qを用意しますが、計算が分かりやすいように小さい素数にしてみます。
p=193,q=11としましょう、次に2つの素数p,qの乗算した値をnとします。

#coding:utf8

p = 193
q = 11
n = p * q


次に(q-1,p-1)の最小公倍数をlとして求めます。
p=193,q=11の場合だと最小公倍数は960になります。

def gcd(a, b):
    """最大公約数を求める"""
    while b:
        a, b = b, a%b
    return a

def lcm(a, b):
    """最小公倍数を求める"""
    return a * b / gcd(a,b)

l = lcm(p - 1,q - 1)  #l = 960


lより小さく、かつlと「素の関係」の適当な数eを選びます。
但し、3は強度が弱いので3以外の数を選んでください、ここでは「77」を選んでみます。

e = 77
print gcd(e,l) #=> 1


次に「ed = 1 (mod l)」を満たすdを求めます、これはeとlが求まっていれば計算で求める事ができます。
eとlは素の関係ですので「gcd(e, l) = 1」ですよね、ここで「拡張ユークリッドの互除法」を使って「ex + ly = 1」のxとyを求めます。
(「ex + ly = 1」を変形すると「ex = ly-1」になるよ!)
求まったxがdとなります、もし「x < 0」になった場合はdにlを足します。

def gcd2(a, b):
    """拡張ユークリッド互除法"""
    if b == 0:
        u = 1
        v = 0
    else:
        q = a / b
        r = a % b
        (u0, v0) = gcd2(b, r)
        u = v0
        v = u0 - q * v0
    return (u, v)

d = gcd2(e,l)[0]
if d < 0:
    d += l


これで暗号化に必要なパラメータe,d,nが揃いました。

「公開鍵」はeとn、「秘密鍵」はdとnの値を使います。

それぞれの値は、e=77,d=773,n=2123となっているはずです。
暗号化は「x^e mod n」とするだけです、さっそく暗号化してみましょう。

#暗号化
org = [13,16,246,21,65,32]

crp = [pow(x,e,n) for x in org]
#crp => [1965, 586, 1589, 1902, 890, 967]


元データと関係ない値になりましたね、それでは復号してみましょう。
復号は「x^d mod n」するだけです。

#復号化
print [pow(x,d,n) for x in crp]
#[13, 16, 246, 21, 65, 32]


おー、復号化できました!!

長くなってしまったので、「ふっかつのじゅもん」の作成は次回に続きます。


RSA暗号で「ふっかつのじゅもん」を作る(2) -RSA暗号鍵の生成 - Pashango’s Blogへ続きます。