素因数分解(2) p-1法、p+1法、楕円曲線法
前回の記事(素因数分解(1) Pollardのρ法 - wacchoz’s note)でちらっと書きましたが、素因数分解のアルゴリズムは大きなカテゴリーが2つあり、そのうちの1つに法、法、楕円曲線法が属しています。
今回はそれらについて書いていきます。
素因数分解(1) Pollardのρ法 - wacchoz’s note
素因数分解(2) p-1法、p+1法、楕円曲線法 - wacchoz’s note(今ここ)
素因数分解(3) フェルマー法、2次ふるい法 - wacchoz’s note
概要
簡単に概要だけ説明すると、法は、のとき、
となることから、約数の多いをとり(具体的には小さな素因数の積とし)、がの倍数となっていれば、
はの倍数となり、たいていの場合はそのものが見つかります。
この方法はが小さな素因数の積になっている場合に、素因数が見つかります。
法も類似した方法で、が小さな素因数の積になっている場合に、素因数が見つかる方法です。
ある数列を計算すると
となるため、約数の多いをとり、がの倍数となっていれば、
はの倍数となり、たいていの場合はそのものが見つかります。
これらはやが小さな素因数の積になっていないとうまくいかず、そうでなければ、をひたすら大きくしていかなければ素因数が見つかりません。
楕円曲線法ではその欠点を克服しています。
で定まる楕円曲線の有理点の集合は、ある演算に対してアーベル群となります。
ここで、群の位数はを変化させると、様々な値をとります。
Hasseによってが証明されていますし、を変化させると、この範囲の整数値をすべて取りうることがDeuringにより証明されています。
ここで、と書くと、上の点に対して、が成り立ちます。
したがってを小さな素数の積として選び、が成り立てば、素因数が見つかります。
法や法はやが小さな素因数の積になっていないとうまくいかず、そうでなければ、をひたすら大きくしていかなければ素因数が見つかりませんでした。楕円曲線法ではを振ることで群の位数を変化させることができ、それらの中には小さな素因数の積になっているものもある可能性があり、法や法よりも素因数が見つかる可能性が高くなっています。
p-1法
上の概要でほぼ語りつくしてしまった感じですが、実装までやってみたいと思います。
整数の素因子の1つがだったとします。
のとき、
となることから、をの倍数とすると
が成り立ちます。
つまりがの倍数となるので、はの倍数となり、たいていの場合、素因子が見つかります。
問題はをどのようにとるかですが、が不明なので、の倍数も当然不明です。そこでとしては約数のできるだけ多い数であれば、素因数分解できる可能性が高くなると考えます。大きさの割に約数が多いものとなると、たくさんの小さな素数の積ということになります。
具体的にはある限界を決めて、以下の素数のべきをすべて掛け合わせたものとします。
は計算時間との兼ね合いになりますが、できるだけ大きめにとります。
Pythonによる実装
以下の実装は『UBASICによるコンピュータ整数論』(日本評論社)を参考にしています。
が非常に大きくなってしまうので、を一気に計算せずに、ごとに計算しています。
from math import gcd # version>=3.5 def isprime(n): if n<2: return False elif n==2 or n==3: return True elif n%2==0: return False else: i = 3 while i*i <= n: if n%i == 0: return False i+=2 return True if __name__ == '__main__': N = 2885500349 L=1000 a=2 aM = 2 for p in range(2,L+1): if not isprime(p): continue Mp = p while Mp*p <= L: Mp*=p #a^m mod Nの計算(素数pごとに計算) aM = pow(aM, Mp, N) g = gcd(aM-1, N) print(g)
実行するとわかりますが、が小さいと答えが1になり素因数が見つかりません。
実装はしませんが、1になったらを大きくしていく必要があります。
また大きすぎて全ての素因子に対してが成り立つと、になってしまうので、その場合は答えが0になってしまいますので、を小さくする必要があるはず。
法ではが小さな素数の積になっていないときはうまくいきませんが、(小さな素数の積)×(ちょっと大きな素数)ぐらいならなんとかなります。(ちょっと大きな素数)をとすると、同様にが素因数となる可能性が高い。も適当な上限を決め、が小さい方から順々にを計算していきます。を増やす時には、べき乗の計算は差分だけにすれば良いです。
前半部を1st step、後半部を2nd stepと呼びます。これは法でも同じようなことが出来ます。
2nd stepまでやっても、に大きな素因子が複数あったり、1つの素因子であってもそれが巨大すぎる場合にはお手上げです。
p+1法
この方法はが小さな素数の積になっているときにうまくいく方法です。
以下の説明は『コンピュータと素因子分解』(和田秀男, 遊星社)によります。
なるを選び、以下の数列を定義します。
の判別式を, 2根を, とおくと、
が成り立ちます。
同様に
と定めると
となる。
単純な計算から
が示せる。
二項定理を用いると
となり、を素数とすると
となります。 まず、です。
はがで平方剰余かどうかによって、となります。
仮にとすると、となるため、
となり
となります。
同様に、仮にとすると、となるため、
となり
となります。
まとめると、がで平方剰余かどうかによってまたはが成り立ちます。
ここでp-1法と同様にを小さな素数の積とします。
例えば、が成立しているとすると、がの約数であればとなり、が素因数となっている可能性が高いというわけです。
事前にが不明なため、がで平方剰余かどうかはわからないので、もしが成り立っていれば、p-1法と同じになってしまいます。(しかもp-1法よりも遅い)
したがって、の組合せを何通りか試す必要があります。
さて大きなに対してを計算するにはどうすればいいでしょうか。
が成り立つため、繰り返し2乗法と同様の考え方で、で計算できます。のペアからまたはのペアに遷移できるためです。
楕円曲線法
楕円曲線の有理点が群をなすことについては、数論好きには有名事実なのですが、説明をするのがなかなか大変なので、tsujimotterさんの記事が非常によくまとまっているので、こちらを参照下さい。
tsujimotter.hatenablog.com
tsujimotter.hatenablog.com
Wikipediaにも結果だけならある程度書いてあります。
楕円曲線 - Wikipedia
となるに対し
の解の全体に無限遠点を追加した集合に対して、以下のように演算を定義します。
2点に対して、の座標を、
の場合、
での場合、
とおき、
と定義します。
での場合、
と定義します。
これらの計算はすべてで計算します。(除算は逆元を掛けます)
これが群の定義を満たします。(結合法則の証明が面倒ですが)
この群をと書くことにします。
と書くと、上の点に対して、が成り立ちます。
はを回加算することを意味します。
をの倍数とするとが成り立ちます。
ここからは法と同じですが、ある限界を決めて
と定め、が成り立つかどうかを調べていきます。
を動かすと、の範囲の整数値をすべてとります。
それらの中には小さな素数の積になっているものをありえます。これが法にまさる理由です。
法ではうまくいかないときにを大きくするしかありませんでした。
楕円曲線法ではうまくいかなければ、を変えればが小さな素数の積になる可能性に期待できます。
具体例を考えていきます。
例としてとして、上の点を次々に倍していくと次のようになり、11倍したところでになります。
1*P = (1,3) 2*P = (2,16) 3*P = (13,11) 4*P = (11,13) 5*P = (6,9) 6*P = (6,8) 7*P = (11,4) 8*P = (13,6) 9*P = (2,1) 10*P = (1,14) 11*P = O
さて、で話を進めてきましたが、実際にはで計算することになります。もちろんが素因数分解したい数です。
環は体ではないので、必ずしも逆元を持つわけではありません。
上の例で考えていきます。をmodとして計算すると以下のようになります。
1*P = (1,3) 2*P = (223,135) 3*P = (13,147) 4*P = (130,64) 5*P = (193,94) 6*P = (57,280) 7*P = (215,242) 8*P = (200,57) 9*P = (291,290) 10*P = (171,252)
を計算しようとするとを計算することになるが、ももmod 17では合同であるためでは逆元が存在しません。しかしながらにより素因数が見つかることになります。
つまりを計算する際に、であれば逆元をもち計算できますが、
であれば逆元が存在しません。しかしながらはの非自明な約数であり、素因数分解できたことになります。
のときは不運で、逆元の計算ができないうえに何も言えません。を小さくして再計算するか、違う曲線から再スタートするかのどちらかとなります。
アルゴリズム
- の大きさに合わせてを定め、とします。
- をランダムにとります。
- となるをとる。
- 点に対しその倍を計算します。もし途中で割り算の計算ができなくなれば、その分母からの約数が求まります。
- もし倍が計算できれば、別のから繰り返します。
なお実際には、ステップ2,3は、をランダムにとってから、を定めるように決めます。また計算にはは使われないので求める必要もありません。
Pythonによる実装
例外を投げて答えを返すクソ実装なので、参考程度に。
from copy import copy import math from math import gcd # version >= 3.5 import random def isprime(n): if n<2: return False elif n==2 or n==3: return True elif n%2==0: return False else: i = 3 while i*i <= n: if n%i == 0: return False i+=2 return True def egcd(a, b): if a == 0: return (b, 0, 1) else: g, y, x = egcd(b % a, a) return (g, x - (b // a) * y, y) # 例外送出用 class ModInvException(Exception): def __init__(self, x, string): self.x = x self.string = string def modinv(a, m): a = a % m g, x, y = egcd(a, m) if g != 1: raise ModInvException(a, 'modular inverse does not exist. (%s,%s)' % (a,m)) else: return x % m class Point(object): a=0 N=0 def __init__(self, x=0, y=0): self.x = x self.y = y def set_curve(a, N): Point.a = a Point.N = N def __add__(self, right): if self.x != right.x: slope = (right.y-self.y) * modinv(right.x-self.x, self.N) elif self.x == right.x and self.y==right.y and self.y!=0: slope = (3*self.x**2 + Point.a) * modinv(2*self.y, self.N) else: raise Exception('result is O!') x3 = (slope**2-self.x - right.x) % self.N y3 = (slope * (self.x - x3) - self.y) % self.N return Point(x3, y3) def __mul__(self, right): if right<1: raise Exception("") r = right - 1 result = copy(self) P = copy(self) while r > 0: if r%2==1: result = (result + P) P = P + P r = r//2 return result def __repr__(self): return "(%s,%s)" % (self.x, self.y) def elliptic_factorization(N): L=min(20000, int(1000*math.exp((math.log(math.sqrt(N))-21.8)*0.13))+100) print("L=",L) while True: a = random.randint(0,N) x = random.randint(0,N) y = random.randint(0,N) Point.set_curve(a,N) P = Point(x,y) for p in range(2,L+1): if not isprime(p): continue Mp = p while Mp*p <= L: Mp*=p try: # M*P mod Nの計算(素数pごとに計算) P = P * Mp except ModInvException as err: print(err.string) divisor = gcd(N, err.x) return divisor except: continue if __name__ == '__main__': N = 18**16+1 divisor = elliptic_factorization(N) print("divisor=", divisor) if N == divisor * (N//divisor): print("%d = %d * %d" % (N, divisor, N//divisor)) else: print("エラー")
出力例
L= 1287 modular inverse does not exist. (54602958734510002506,121439531096594251777) divisor= 3930785153 121439531096594251777 = 3930785153 * 30894471809