単純な素因数分解アルゴリズムを実装してみる

「OpenSSLとPythonでRSA暗号の原理を知る」で書いたように、RSA暗号は公開鍵である整数modulusを二つの素数に素因数分解することができれば破ることができる。 ここでは、いくつかの単純な素因数分解アルゴリズムをPythonで実装し、一般的なPCで解くことができるbit数などを調べてみる。 なお、正確には素因数分解ではなく、因数をひとつ見つけるアルゴリズムについて考えるものとする。

環境

Ubuntu 14.04.3 LTS 64bit版、Intel Core i5-4200U (1.6GHz * 2 * 2)

$ uname -a
Linux vm-ubuntu64 3.19.0-25-generic #26~14.04.1-Ubuntu SMP Fri Jul 24 21:16:20 UTC 2015 x86_64 x86_64 x86_64 GNU/Linux

$ lsb_release -a
No LSB modules are available.
Distributor ID: Ubuntu
Description:    Ubuntu 14.04.3 LTS
Release:        14.04
Codename:       trusty

$ grep "processor\|model name" /proc/cpuinfo
processor       : 0
model name      : Intel(R) Core(TM) i5-4200U CPU @ 1.60GHz
processor       : 1
model name      : Intel(R) Core(TM) i5-4200U CPU @ 1.60GHz
processor       : 2
model name      : Intel(R) Core(TM) i5-4200U CPU @ 1.60GHz
processor       : 3
model name      : Intel(R) Core(TM) i5-4200U CPU @ 1.60GHz

素数の生成

素数とその積の生成は、OpenSSLを用いて行うこととする。 128bitのmodulusを生成する場合は次のようになる。

$ openssl genrsa 128 | openssl rsa -noout -text
Generating RSA private key, 128 bit long modulus
.+++++++++++++++++++++++++++
.+++++++++++++++++++++++++++
e is 65537 (0x10001)
Private-Key: (128 bit)
modulus:
    00:bf:27:4e:d6:73:f6:ca:f4:f6:85:f0:b6:a0:27:
    1a:31
publicExponent: 65537 (0x10001)
privateExponent:
    76:67:8f:d3:10:e3:d7:14:d1:a8:40:83:b3:81:1b:
    ad
prime1: 16144814945699257943 (0xe00de7a77ad25e57)
prime2: 15737972014275192503 (0xda688229e3e966b7)
exponent1: 9570317772029090007 (0x84d0994ed566d8d7)
exponent2: 9046986552112602077 (0x7d8d5a560e685bdd)
coefficient: 9750135051940033389 (0x874f702bdaf8736d)

試し割り

もっとも単純なアルゴリズムであり、整数nに対してsqrt(n)以下の整数で順に割っていくというもの。 計算量はO(n^1/2)である。

# trial_division.py

import sys

def isqrt(n):
    x = n
    y = (x + 1) // 2
    while y < x:
        x = y
        y = (x + n//x) // 2
    return x

def trial_division(n):
    k = isqrt(n)

    if n % 2 == 0:
        return 2
    if n % 3 == 0:
        return 3

    i = 6
    while i < k:
        if n % (i-1) == 0:
            return i-1
        if n % (i+1) == 0:
            return i+1
        i += 6

if __name__ == '__main__':
    n = long(sys.argv[1])
    p = trial_division(n)
    if p:
        print "%d = %d * %d" % (n, p, n/p)
    else:
        print "%d is prime" % n

上のコードでは、2と3の倍数となるものについては飛ばすよう若干の効率化を加えている。 また、isqrt関数は整数の範囲で平方根を求める関数である。

次の実行結果の通り、64bitの素因数分解に5分程度かかる。

$ time python trial_division.py 12814570762777948741
12814570762777948741 = 3318288047 * 3861801803

real    4m47.971s
user    4m47.296s
sys     0m0.596s

$ time python trial_division.py 18366865165381711817
18366865165381711817 is prime

real    6m15.418s
user    6m15.348s
sys     0m0.000s

Pollard-Rhoアルゴリズム + Miller-Rabin素数判定法

Pollard-Rhoアルゴリズムは、乱択アルゴリズムにより高速に因数を発見するアルゴリズムである。 整数nに対しその因数をpとしたとき、計算量はO(p^1/2) <= O(n^1/4)となる。

乱択アルゴリズムという性質上、与えられた数が素数だった場合、停止するまでに非常に時間がかかる。 そのため、同じく確率的アルゴリズムであるMiller-Rabin素数判定法と組み合わせ、あらかじめ素数かどうかを判定しておくことが多い。

# pollard_rho.py

import sys
from random import randint
from fractions import gcd

def miller_rabin(n, k=20):
    s, d = 0, n-1

    while d % 2 == 0:
        s += 1
        d /= 2

    for i in xrange(k):
        a = randint(2, n-1)
        x = pow(a, d, n)
        if x == 1:
            continue
        for r in xrange(s):
            if x == n-1:
                break
            x = (x*x) % n
        else:
            return False

    return True

def pollard_rho(n):
    x, y, d = 2, 2, 1

    while d == 1:
        x = (x*x + 1) % n
        y = (y*y + 1) % n
        y = (y*y + 1) % n
        d = gcd(abs(x-y), n)

    if d != n:
        return d

if __name__ == '__main__':
    n = long(sys.argv[1])
    is_prime = miller_rabin(n)
    if is_prime:
        print "%d is prime" % n
    else:
        p = pollard_rho(n)
        print "%d = %d * %d" % (n, p, n/p)

Miller-Rabin素数判定法において、素数判定を誤る確率は繰り返し回数kに対し4^-kとなることが知られている。 上のコードではk = 20としている。

次の実行結果の通り、96bitの素因数分解に1分程度かかる。 なお、素数の場合についてはMiller-Rabin素数判定法により1秒未満で判定できる。

$ time python pollard_rho.py 60766145992321225002169406923
60766145992321225002169406923 = 250117558771727 * 242950340194949

real    1m12.204s
user    1m12.192s
sys     0m0.000s

$ time python pollard_rho.py 71939287897297826407363026419
71939287897297826407363026419 is prime

real    0m0.029s
user    0m0.024s
sys     0m0.004s

フェルマー法

フェルマー法はRSAのように素数同士の積からなる合成数において、二つの素数の差が小さい場合に有効なアルゴリズムである。 具体的には、合成数nの平方根周辺のある値をxとしたとき、(x-k)*(x+k) == nとなるkを小さいものから求める。

# fermat.py

import sys

def isqrt(n):
    x = n
    y = (x + 1) // 2
    while y < x:
        x = y
        y = (x + n//x) // 2
    return x

def is_square(n):
    if not n % 48 in (0, 1, 4, 9, 16, 25, 33, 36):
        return False

    x = isqrt(n)
    return x*x == n

def fermat(n):
    a = isqrt(n)
    b2 = a*a - n
    while not is_square(b2):
        a += 1
        b2 = a*a - n
    return a - isqrt(b2)

if __name__ == '__main__':
    n = long(sys.argv[1])
    p = fermat(n)
    if p:
        print "%d = %d * %d" % (n, p, n/p)
    else:
        print "%d is prime" % n

上のコードにおいて、is_square関数は値が平方数かどうかを判定する関数である。 平方数の判定については、平方剰余をもとに平方数となりえない値を先に判定することで効率化が図ることができ、上では法48について一度判定することで83.4%の剰余について計算を省いている。

次の実行結果の通り、64bitの素因数分解に30秒程度かかる。 素数の場合については大きな時間がかかるが、Pollard-Rhoアルゴリズムの場合と同様にMiller-Rabin素数判定法を組み合わせれば高速に判定することが可能である。

$ time python fermat.py 12814570762777948741
12814570762777948741 = 3318288047 * 3861801803

real    0m23.931s
user    0m23.924s
sys     0m0.000s

また、kが小さい場合についてはnの大きさによらず高速に素因数分解できる。 ただし、そのような素数の組が選ばれる確率は非常に小さい。

次の実行結果では、nが4096bit、kが15bitの場合について、1秒未満で素因数分解できている（参考）。

$ time python fermat.py 1034186074668005446865661651578208369514158803552244005581133348807218259696136975186768057532805396894757503164834709575083359161158060701923008831149233177215596305395749862300700772298386522636562203611668211601359078939366452933448401372783149865070405032946502692195960337152490393257718681775969357953766124434670452760182942622011519373645653657784274264774773270847089540432369489508909290697992927057404389218509541644042197448296417276686963995927414541514778710278820030649227976349765070908207584414253880096569270321293357683127769728901749940165024782079154839570419235490028108256539125090626571373643257455521820767805808214404175900736109428513289545954241551995058498737472560704736894014465695117165251699935296644992760434412206519036326801985540077086806228464411438175416317958133959586691011057808141234300883692932436022122663781160823472736298914304029803683697389880156222528045623045436226522780569221934366978973936408257173404970887883159180079205465195884955186605310665826226975454171035075016280417683011623124380312403051425940705920625660567584606796218038626946234225864706878239646554381165540673005572649955471689205027195933802181344480474070023069194083134722979170599617466579222550596280381011
1034186074668005446865661651578208369514158803552244005581133348807218259696136975186768057532805396894757503164834709575083359161158060701923008831149233177215596305395749862300700772298386522636562203611668211601359078939366452933448401372783149865070405032946502692195960337152490393257718681775969357953766124434670452760182942622011519373645653657784274264774773270847089540432369489508909290697992927057404389218509541644042197448296417276686963995927414541514778710278820030649227976349765070908207584414253880096569270321293357683127769728901749940165024782079154839570419235490028108256539125090626571373643257455521820767805808214404175900736109428513289545954241551995058498737472560704736894014465695117165251699935296644992760434412206519036326801985540077086806228464411438175416317958133959586691011057808141234300883692932436022122663781160823472736298914304029803683697389880156222528045623045436226522780569221934366978973936408257173404970887883159180079205465195884955186605310665826226975454171035075016280417683011623124380312403051425940705920625660567584606796218038626946234225864706878239646554381165540673005572649955471689205027195933802181344480474070023069194083134722979170599617466579222550596280381011 = 32158763574926282399690427421751598974822750157866002942864427634852437380540017586451493854661729909380518733649186624385206737336324813109500237603304026009112696565510846849987937423619262973393969175056759821652138869783215378169757835283584660846583208812725733839059137580944002686113912792569631796916069732431775599320458346937859589815497525828622622652165709271152246464728489927670696601016248559515951932154686633599100945314921834227324381958751184684979824241375253606863601895383658582486045363570755445629865194046700806542078378801136397577730247660070033517187439537339428288763342861366560261414507 * 32158763574926282399690427421751598974822750157866002942864427634852437380540017586451493854661729909380518733649186624385206737336324813109500237603304026009112696565510846849987937423619262973393969175056759821652138869783215378169757835283584660846583208812725733839059137580944002686113912792569631796916069732431775599320458346937859589815497525828622622652165709271152246464728489927670696601016248559515951932154686633599100945314921834227324381958751184684979824241375253606863601895383658582486045363570755445629865194046700806542078378801136397577730247660070033517187439537339428288763342861366560261446073

real    0m0.040s
user    0m0.028s
sys     0m0.012s

Common Factor Attacks

厳密には素因数分解アルゴリズムではないが、異なる二つの合成数が共通の素因数を持つ場合、合成数のbit数によらず高速に素因数を発見することができる。 具体的には、二つの合成数の最大公約数を求めるだけでよい。 最大公約数はユークリッドの互除法によりO(log n)で計算できる。

# common_factor.py

import sys
from itertools import combinations
from fractions import gcd

def common_factor(nums):
    for x, y in combinations(nums, 2):
        p = gcd(x, y)
        if p > 1:
            yield (x, y, p)

if __name__ == '__main__':
    nums = [long(line) for line in sys.stdin]
    for (x, y, p) in common_factor(nums):
        print "%d = %d * %d" % (x, p, x/p)
        print "%d = %d * %d" % (y, p, y/p)

次の実行結果の通り、共通の素因数を持つ二つの合成数において、4096bitの素因数分解が1秒未満でできている。 ただし、フェルマー法同様、そのような合成数の組が見つかる確率は非常に小さい。

$ cat nums.txt
888441243015294102172292601067855336326287332703328874021611868561802676986077663943061900147565825499428361944793850926519424097012115662952834330500310063991436323100617725235375735298880799081272752074568088718329068470304011938234831484558254069480648619902970082268415966404044075480346440255780228196130123810436581349922045986882478809572371947975635740413437554435429379257619351276688537073428555742828688693292492859848882427167863812003144627707272854531325119129193487170520323833544069156767215691882313314486711034734316444571037196897906716976294880217146223065562518732008190599420275772517007856746981048559878253715189698713288512904998165116853386358827890646734214221734467879782398263142454822004494788965625480385484067919453526130347620614241545718696612964390044265731882751353123553453763779825959232927147708421269420912056959539697543946480513329283728353598683158139025560947346471410872276377801166133810455502905521701767968469458987925204027876348478451236285710225645598004719391738052126832244276934547442298148805426263554624932225480952697394164120418952036473986289321176380982895752273796339856601757641830247264277014250785556784199680489778112198639175574756714183440252052420135756327849404813
882780483105414285207808409377995483701534823624777052872370512205663687323012030665008836246779641378802514783938670119510387540080290753721597154663983983212286475091452225736450381989144354075664524345788834071829669406861239991491714910029461119046332512912612670883516196208270591616798036275990449980681900277329633024676002455118771425041101638785071781244260968665572605189724603728273275682283698062764859089709731324988332263001477198449644990702495668478665515369529879885560385782218986658365704612364761778013317060972940381244288115042449405877429592775144911345025062002681580353283005136150916508175831243673046417028814562130687594141584703791313070099191199276692294970536889058503070974721545833923098157515201071570986923012040954935154602172279935762508063846354593464196649444856619523404920553516291837255896358713883164374689514763790547373493144894084431571218539012930850277919063121821877687816010348113929847662784908586555748228552075893306686302435532676222920629055063379540147073916377160422832146045388449909494920704515440022193230206490791891945947544516000231822642964432381556023595277757817268369138829082830952126257601345645720002969113016584602034934040121873671897644510285530925114503625051

$ time { cat nums.txt | python common_factor.py; }
888441243015294102172292601067855336326287332703328874021611868561802676986077663943061900147565825499428361944793850926519424097012115662952834330500310063991436323100617725235375735298880799081272752074568088718329068470304011938234831484558254069480648619902970082268415966404044075480346440255780228196130123810436581349922045986882478809572371947975635740413437554435429379257619351276688537073428555742828688693292492859848882427167863812003144627707272854531325119129193487170520323833544069156767215691882313314486711034734316444571037196897906716976294880217146223065562518732008190599420275772517007856746981048559878253715189698713288512904998165116853386358827890646734214221734467879782398263142454822004494788965625480385484067919453526130347620614241545718696612964390044265731882751353123553453763779825959232927147708421269420912056959539697543946480513329283728353598683158139025560947346471410872276377801166133810455502905521701767968469458987925204027876348478451236285710225645598004719391738052126832244276934547442298148805426263554624932225480952697394164120418952036473986289321176380982895752273796339856601757641830247264277014250785556784199680489778112198639175574756714183440252052420135756327849404813 = 30094504000188239020151081364836978642515567675718773053464099309910028925560991426570939310582088936607583861708569872351068260656205876930160648487885561172643428908022945411031069520952948240317537081178371267446884744948251284183621948670880897605924632767128833848853213483976135804220279591102249759334548097546952384012455150437958938593286685443258091084425290523032227270109892468085164321266154276823025471164370162127401053735334983910367679071627280601662913673154344586359423856921429216035850978379728954418530655343614680839905243246402183753209499080677780527109087925474741203711791946310699728133647 * 29521710775154724738536772527584123117558710613185797920601145487104166977711662880879179271048876405095288045905017016692076383141051745304194231171302588112692719379484816312610044960068902592215853578316022870490297200043945853963854204074128066849603964094184178301412655771436359730063212919768104154891433003679579996888043953631799144984305859356177785182529634481307702685821400386110992830487859752626831497462372631255954636435957631272568434117771136620170717336511538556050260061771074028879608580486351337906026632036083341002564850216617332234416072504420109009193206508593690307919068066529130502806179
882780483105414285207808409377995483701534823624777052872370512205663687323012030665008836246779641378802514783938670119510387540080290753721597154663983983212286475091452225736450381989144354075664524345788834071829669406861239991491714910029461119046332512912612670883516196208270591616798036275990449980681900277329633024676002455118771425041101638785071781244260968665572605189724603728273275682283698062764859089709731324988332263001477198449644990702495668478665515369529879885560385782218986658365704612364761778013317060972940381244288115042449405877429592775144911345025062002681580353283005136150916508175831243673046417028814562130687594141584703791313070099191199276692294970536889058503070974721545833923098157515201071570986923012040954935154602172279935762508063846354593464196649444856619523404920553516291837255896358713883164374689514763790547373493144894084431571218539012930850277919063121821877687816010348113929847662784908586555748228552075893306686302435532676222920629055063379540147073916377160422832146045388449909494920704515440022193230206490791891945947544516000231822642964432381556023595277757817268369138829082830952126257601345645720002969113016584602034934040121873671897644510285530925114503625051 = 30094504000188239020151081364836978642515567675718773053464099309910028925560991426570939310582088936607583861708569872351068260656205876930160648487885561172643428908022945411031069520952948240317537081178371267446884744948251284183621948670880897605924632767128833848853213483976135804220279591102249759334548097546952384012455150437958938593286685443258091084425290523032227270109892468085164321266154276823025471164370162127401053735334983910367679071627280601662913673154344586359423856921429216035850978379728954418530655343614680839905243246402183753209499080677780527109087925474741203711791946310699728133647 * 29333611316534492642114695914361875827353936602039611693387782601066021747108940770284205684743192567435597448844270242635706940719055360846091860177842785258548515263955446733702873336801213730696067625102999759204582858629776679167503390124589638053007027785917641140606694678963338089301468801125801144536374695406033850449781597771312803137398214729077301657987716032422612267398441299477065463230629907039871180083779210842889144946322439175733816114014531390999811837034265672904377154925113301930277530009488651707667660960383253017461910446471490714077271746064722611111519237800570978700971166429844348550133

real    0m0.029s
user    0m0.028s
sys     0m0.000s

関連リンク

• POJ 2429 - GCD & LCM Inverse - fuqinho@競技プログラミング
• a4lg.com | 作問者から見た SECCON 2014 オンライン予選 Crypto 400 (2014-12-14)
• Understanding Common Factor Attacks: An RSA-Cracking Puzzle
• 平方数かどうかを高速に判定する方法 - hnwの日記
• Quadratic residue - Wikipedia, the free encyclopedia