素因数分解(3) フェルマー法、2次ふるい法
以前の記事でもちらっと書きましたが、素因数分解のアルゴリズムは大きなカテゴリーが2つあり、そのうちの1つにフェルマー法、2次ふるい法、連分数法、数体ふるい法が属しています。
今回はフェルマー法、2次ふるい法を紹介していきます。
素因数分解(1) Pollardのρ法 - wacchoz’s note
素因数分解(2) p-1法、p+1法、楕円曲線法 - wacchoz’s note
素因数分解(3) フェルマー法、2次ふるい法 - wacchoz’s note(今ここ)
目次
フェルマー法
素因数分解したい数がと分解できたとします。そうすればと分解でき素因数分解できます。そこでこのようなを探すことにします。
xの初期値として、とします。
yの初期値はとします。
- を計算する
- ならが素因数である
- ならとする
- ならとする
- 1.に戻る
の初期値が小さい値からインクリメントしていることかわかる通り、の2個の素因子の値が近い場合に見つかるのが早いです。
2次ふるい法
2次ふるい法では
となるを見つけることを目標とします。これが見つかれば
がの倍数になるということなので、うまくいけばがに分かれるかもしれません。そうなればが素因数として見つかることになります。
今、を動かしてをで割った余りをたくさん考えます。
です。ここで、いくつかのをうまく組合せて平方数を作れれば、目的が達成されます。
これは各を素因数分解して、各素数の指数を組み合わせて偶数を作ることを意味し、行列にGauss消去法を適用することにより計算できます。
もう少し具体的にすると、のときのとなったとし、それぞれのを素因数分解すると
となったとします。
ここで指数をでベクトルとして並べて
を作ります。
そしてを行とする行列を作り、Gauss消去法で0だけの要素からなる行を作ります。Gauss消去法でどのような操作を行ったかがわかるように単位行列も用意し、単位行列にも同じ操作をしていきます。これでどの行を組み合わせることで0だけの要素からなる行ができるかが決定できます。
それらの行に対応するを掛け合わせれば平方数になるので、それをとすると、目標としていたが得られる。大抵の場合、で素因子が得られます。
平方数の構成法
具体例としてからいくつかを選びだして、積が平方数になるものを探します。
これらの数の素因子はであるので、各値をそれらで素因数分解したときの指数のをベクトルとして並べ、右側に単位行列を並べます。
この左半分をGauss消去すると以下のようになります。
行の左半分がになっており、平方数になったことを示しています。右半分を見るとどのような組合せなのかわかります。
例えば行目は列目が1になっているため、番目の数、を掛け合わせると平方数になるとわかります。
平方数の構成のプログラム例
# -*- coding: utf-8 -*- from pprint import pprint from math import sqrt # 「nをplistの素数で素因数分解したときの、各素数のべき乗の偶奇」を並べたベクトルを返す # 素因数分解できなければNoneを返す def factoring_vector(n, plist): ans = [] for p in plist: e=0 while n%p==0: n //= p e += 1 ans.append(e%2) if n!=1: return None return ans # Gauss消去 def gaussian_elimination(mat, row, col): for j in range(col): for i in range(j,row): if mat[i][j] == 1: pivot_row = i break else: continue if pivot_row != j: mat[pivot_row],mat[j] = mat[j],mat[pivot_row] pivot_row = j for i in range(pivot_row+1, row): if mat[i][j]==1: mat[i] = [(mat[i][k] + mat[j][k])%2 for k in range(len(plist)+len(rlist))] if __name__ == '__main__': plist = [2,3,5,7,11] rlist = [2,3,4,5,6,7,8,9,10,11] mat = [] # 行列を生成する for k,r in enumerate(rlist): vec = factoring_vector(r, plist) # 単位行列のk行目 identity = [ (1 if i == k else 0) for i in range(len(rlist)) ] print(r,vec) if vec is None: continue mat.append(vec + identity) pprint(mat) # Gauss消去 row = len(mat) col = len(plist) gaussian_elimination(mat, row, col) pprint(mat) # 行列の左半分が0になっている行を探す for i in range(row): allzero = True for j in range(col): if mat[i][j] == 1: allzero = False break if(allzero): n=1 for j in range(len(rlist)): if mat[i][col+j] == 1: n *= rlist[j] if n!=1: print(n,"は",sqrt(n),"の2乗")
2 [1, 0, 0, 0, 0] 3 [0, 1, 0, 0, 0] 4 [0, 0, 0, 0, 0] 5 [0, 0, 1, 0, 0] 6 [1, 1, 0, 0, 0] 7 [0, 0, 0, 1, 0] 8 [1, 0, 0, 0, 0] 9 [0, 0, 0, 0, 0] 10 [1, 0, 1, 0, 0] 11 [0, 0, 0, 0, 1] [[1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0], [1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0], [1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]] [[1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1], [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0]] 4 は 2.0 の2乗 16 は 4.0 の2乗 9 は 3.0 の2乗 100 は 10.0 の2乗 36 は 6.0 の2乗
さて、このプログラムで少し不満なのが、0か1しか値をとらないのに1変数分のメモリを食っているところです。
そこで各行を2進法のbit列とみて、1つの変数で表すことにします。このようにするとメモリの節約になり、Gauss消去の箇所にbit XORが使えるようになります。
# -*- coding: utf-8 -*- from math import sqrt # 「nをplistの素数で素因数分解したときの、各素数のべき乗の偶奇」をbitとして低位桁より並べた整数を返す # 素因数分解できなければNoneを返す def factoring_vector(n, plist): ans = 0 for i,p in enumerate(plist): while n%p==0: n //= p ans ^= 1<<i if n!=1: return None return ans # nのkビット目(0-indexed) def bit(n, k): return (n>>k) & 1 # bit列を逆順にした文字列(出力用) def bin_reverse(n, k): return format(n, 'b')[::-1].ljust(k,'0') # Gauss消去 def gaussian_elimination(mat, row, col): for j in range(col): for i in range(j,row): if bit(mat[i],j) == 1: pivot_row = i break else: continue if pivot_row != j: mat[pivot_row],mat[j] = mat[j],mat[pivot_row] pivot_row = j for i in range(pivot_row+1, row): if bit(mat[i],j)==1: mat[i] ^= mat[j] if __name__ == '__main__': plist = [2,3,5,7,11] rlist = [2,3,4,5,6,7,8,9,10,11] mat = [] for k,r in enumerate(rlist): v = factoring_vector(r, plist) if v is None: continue # 単位行列のk行目 identity = 1 << k mat.append(v + (identity << len(plist))) for v in mat: print(bin_reverse(v, len(plist)+len(rlist))) # Gauss消去 row = len(mat) col = len(plist) gaussian_elimination(mat, row, col) print("-----------") for v in mat: print(bin_reverse(v, len(plist)+len(rlist))) for i in range(row): if (mat[i] & ((1<<col)-1)) > 0: continue n=1 for j in range(len(rlist)): if bit(mat[i],col+j) == 1: n *= rlist[j] if n!=1: print(n,"は",sqrt(n),"の2乗")
100001000000000 010000100000000 000000010000000 001000001000000 110000000100000 000100000010000 100000000001000 000000000000100 101000000000010 000010000000001 ----------- 100001000000000 010000100000000 001000001000000 000100000010000 000010000000001 000000010000000 000001000001000 000000000000100 000001001000010 000001100100000 4 は 2.0 の2乗 16 は 4.0 の2乗 9 は 3.0 の2乗 100 は 10.0 の2乗 36 は 6.0 の2乗
「ふるい」とは?
さて、問題をもう一度振り返ると、を動かしてをで割った余りをたくさん考え、いくつかのをうまく組合せて平方数を作るのでした。
は小さい方が平方数が作りやすいため、はの近くで動かすこととしますが、やみくもに動かしてもを素因数分解したときに多数の素数が出てきて効率が悪くなります。
上の例でいうとは1度しか出てこないため、平方数を作るのに貢献していません。
そこでまずあるを決めて、個の素数を決め、これらで分解されるのみを抽出することを考えます。
(はまたはが平方剰余となる奇素数)
を動かしていった際に、の中で素数で割り切れるものがどこに現れるかを考えます。
の解をと計算さえすれば、毎にの倍数が現れることがわかります。
この部分が2次「ふるい」法と呼ばれるゆえんで、エラトステネスの篩のようなプロセスで篩にかけて抽出していきます。
以外では割れないを求めたいので、最初にを計算し、で割り切れればを引いていき、になればのみで割り切れるということになります。
小さい素数の場合には、で割り切れることも多いですが、が奇素数の場合には、同様に1つの解を計算すれば毎に割り切れます。
[木田]では素数のべきで割れるケースは無視して、がある閾値よりも低くなるときに実際に割り算して確かめる、という方法をとっています。実装上はこの方が楽だと思います。
# -*- coding: utf-8 -*- from math import sqrt, exp, log 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 # Legendre記号 # p: 素数 def legendre(a, p): # p=2のときは1を返しておく if(p==2): return 1 # 奇素数の場合 ret = pow(a, (p-1)//2, p) # print(a,p,ret) if ret==1: return 1 else: return -1 # x^2 ≡ n (mod p)の解の1つを返す def modsqrt(n, p): for i in range(p//2+1): if((i*i -n)%p==0): return i return None # 「nをplistの素数で素因数分解したときの、各素数のべき乗の偶奇」をbitとして低位桁より並べた整数を返す # 素因数分解できなければNoneを返す def factoring_vector(n, plist): ans = 0 for i,p in enumerate(plist): while n%p==0: n //= p ans ^= 1<<i if n!=1: return None return ans # nのkビット目(0-indexed) def bit(n, k): return (n>>k) & 1 # bit列を逆順にした文字列(出力用) def bin_reverse(n, k): return format(n, 'b')[::-1].ljust(k,'0') # Gauss消去 def gaussian_elimination(mat, row, col): for j in range(col): for i in range(j,row): if bit(mat[i],j) == 1: pivot_row = i break else: continue if pivot_row != j: mat[pivot_row],mat[j] = mat[j],mat[pivot_row] pivot_row = j for i in range(pivot_row+1, row): if bit(mat[i],j)==1: mat[i] ^= mat[j] def gcd(a, b): if a<b: return gcd(b,a) r = b while a%b!=0: r = a%b a,b = b, r return r def factoring(n): ret = [] p=1 while p*p<=n: p+=1 if n%p!=0: continue num = 0 while n%p==0: num+=1 n //= p ret.append((p,num)) if n!=1: ret.append((n,1)) return ret if __name__ == '__main__': N = 8800969069 S = int(sqrt(N))+1 kosuP = max(5, int(1.8*exp(sqrt(log(N)*log(log(N))/12)))) print("kosuP=",kosuP) kosuX = int(kosuP*1.1)+4 print("kosuX=", kosuX) # 素数列挙 plist = [2] p = 3 while True: if isprime(p) and legendre(N,p)==1: if len(plist)==kosuP: break plist.append(p) p+=2 # max_iは決め打ちとしている。実際にはkosuX個見つかるまで続ける必要あり max_i = 1000000 work = [0] * max_i for i in range(max_i): work[i] = log((S+i)**2 - N) for p in plist: logp = log(p) root = modsqrt(N,p) for t in ([root, p-root] if p!=2 else [root]): a = (t-S) % p assert(((S+a)**2 - N)%p == 0) while a < max_i: work[a] -= logp a += p mat = [] rlist = [] for i,w in enumerate(work): if w < 3.0: r = (S+i)**2 - N v = factoring_vector(r, plist) if v is None: continue # 単位行列のk行目 identity = 1 << len(mat) mat.append(v + (identity << len(plist))) rlist.append(i) if len(mat) >= kosuX: break for v,r in zip(mat,rlist): print(bin_reverse(v, len(plist)+len(mat))) # Gauss消去 row = len(mat) col = len(plist) gaussian_elimination(mat, row, col) print("-----------") for v in mat: print(bin_reverse(v, len(plist)+len(mat))) for i in range(row): if (mat[i] & ((1<<col)-1)) > 0: continue y2 = 1 x2 = 1 for j in range(len(rlist)): if bit(mat[i],col+j) == 1: y2 *= (S+rlist[j])**2 - N x2 *= S+rlist[j] if y2!=1: y = int(sqrt(y2)+0.5) #print(y2,"は",y,"の2乗") if gcd(x2+y,N) * (N//gcd(x2+y,N)) == N: print("divisor=",gcd(x2+y,N), N//gcd(x2+y,N))
kosuP= 20 kosuX= 26 0100001010000100000010000000000000000000000000 0001010010000000000101000000000000000000000000 0100000000100000011000100000000000000000000000 0011000100000000101000010000000000000000000000 0100010100100010000000001000000000000000000000 0010100001010001000000000100000000000000000000 0101100110001000000000000010000000000000000000 0010001100100100000000000001000000000000000000 0110100100000000110000000000100000000000000000 0011000001010000000100000000010000000000000000 0110011000100001000000000000001000000000000000 0010100000100100000100000000000100000000000000 0100100000000100001000000000000010000000000000 0110011000001000100000000000000001000000000000 0100100011001100000000000000000000100000000000 0000000001010010010000000000000000010000000000 0011101011000000000000000000000000001000000000 0010100001010001000000000000000000000100000000 0011000110100000100000000000000000000010000000 0011001000000000010100000000000000000001000000 0010111000101000000000000000000000000000100000 0110000010001001000000000000000000000000010000 0100110100001000000100000000000000000000001000 0101101000100000100000000000000000000000000100 0000111010000000001000000000000000000000000010 0011001000011010000000000000000000000000000001 ----------- 0100001010000100000010000000000000000000000000 0100000000100000011000100000000000000000000000 0011000100000000101000010000000000000000000000 0001010010000000000101000000000000000000000000 0000110111010001101101010100000000000000000000 0000010100000010011000101000000000000000000000 0000001110100110110101111001000000000000000000 0000000101110001101000100100100000000000000000 0000000011111001110000110110000000000000000000 0000000001110101000100000100000100000000000000 0000000000100111000000001001001000000000000000 0000000000011000101101010101001000100000000000 0000000000001110010101101010000010000000000000 0000000000000110000100111101111000000000000000 0000000000000011101101111100110110000010000000 0000000000000001010101000111111100000010010000 0000000000000000110101000110111100101000000000 0000000000000000010100001101001100010000000000 0000000000000000000000100101111111101000000000 0000000000000000000101000100111000000011000000 0000000000000000000000010100111010101001100000 0000000000000000000000000100000000000100000000 0000000000000000000001000010000000000000001000 0000000000000000000000010001000010000000000100 0000000000000000000000110111111110101000000010 0000000000000000000001101110000010100001010001 divisor= 1 8800969069 divisor= 1 8800969069 divisor= 94349 93281 divisor= 1 8800969069 divisor= 1 8800969069 divisor= 1 8800969069 divisor= 1 8800969069
Pythonで機械設計を自動化した話
仕事でとある機械を設計する仕事をしていました。
3年前ぐらいに設計をPythonで自動化した経験について書きます。
入社したときから設計に配属され、ずっと設計をしていました。
内製の多数の設計ツールを駆使して、手作業でひたすら計算して、設計クライテリアを満たす解を探したりして、最終的に2次元CADを作成します。
ツールだけ見るとかなり古臭いです。実際、インターフェースだけは変わったものの、10*n年前から使われているツールが当たり前という環境です。
手作業での作業は楽しくはあり、設計しているという実感を得られるものの、非常に非効率です。
入社当時は長時間残業が当たり前。
設計者は設計以外の雑用にも追われ、設計に取られる時間も限られています。
そんな中で手作業で計算をしているため、設計の最適解は当然のように得られず、もどかしさを感じていました。
設計というのは、ある制約条件下で何らかの最適解を求める作業です。
開発品で初めての製品を設計するならいざ知らず、類似品を多数設計するのであればプログラムで自動化して最適化すべきものです。
入社N年目に当時の課長に自動化を専任でやらせて欲しいと直訴し、専任でやらせてもらえることになりました。
(なぜあっさりやらせてもらえたのかという背景があるのですが、それはまた別の話。)
私の設計業務は基本設計だったので、内製の計算ツールを複数使い、形状を決定して、2次元CADで計画図を作成するのが主な業務です。(他にも雑務的な仕事は数知れず)
概算計画でも、計算慣れしているベテラン勢で1週間ぐらい、そこまで慣れていない若い人では2週間程度がかかるものです。
実際にはこのあとデザインレビューを行い、正式に図面や図書を発行する流れとなり、案件にもよりますが1ヶ月半~2カ月ぐらいかけるのが普通です。
受注の引合いがあったときの見積り設計は期限が短く、数日で概算設計を終わらせる必要がありました。入札に関わる重要なプロセスにも関わらず、あまり検討をしていない設計で応札するということも多かったです。これでは受注も遠ざかるばかりで、これは改善しなければ、という気持ちでした。
具体的な話
さて、どういう設計を手作業でしていたのかを説明していきます。
設計ツールはいろいろあって
のような感じ20個近いツールがあります。
さらにExcelマクロから呼び出すDLLもあります。
GUIと書いているものは、GUIはフロントエンドだけで、実際にはサーバのFortranプログラムを呼び出しています。
具体的にはGUIから入力テキストファイルを出力し、それをサーバにアップロードし、Fortranプログラムを走らせ、出力されたテキストファイルをダウンロードしています。
Excelシートでの計算ツールや、Excelのマクロを使ったツールがあるのは、他の設計現場でもよくあるのではないでしょうか?
実際の設計では例えば以下のような作業をします。
- 計算ツールAの出力から数値をぬきだし、計算ツールBに入力する
- 計算ツールCの出力のある数値が制約以下になるような入力を探す
サーバのプログラムを走らせたり、Excelを使ったりと設計ツールが多岐に渡りますが、これらを自動的に走らせるのはPythonが向いています。C++とかに手を出すと死亡します。
既存の設計プログラムには一切手を加えず、Pythonでプログラム間の操作を仲介するような方針にしました。
Excelの操作の自動化
Excelの操作にはPythonのwin32comモジュールを使いました。
他にも色々なモジュールがあるようですが、特に比較検討をしたわけではないので、今はもっと良いものがあるのかもしれません。
ほとんどExcelマクロのVBAを書くのと同じ感覚で書くことができます。
PythonからExcelマクロを呼び出すこともできます。
Excelファイルを開いて、セルに数字を入力して、マクロを走らせて、数値を読みだして、ファイルを閉じる、といった作業がPythonから出来ます。
なお注意点として、PythonプログラムからExcelを操作している間にExcelを触ると、エラーで終了してしまうのが地味に不便でした。
サーバプログラムの呼び出しの自動化
これは完全に職場固有の問題なので、あまり参考にならないと思います。
多くの設計ツールがFortranプログラムとしてサーバ内にあるため、GUIプログラムからサーバのプログラムを呼び出して、出力をダウンロードするという処理になっていました。
ここをPythonから自動で行うには、主に以下の3つの関数を実装することになります。
- 計算の入力ファイルをファイルに書き出す関数
- 計算を走らせる関数
- 出力ファイルを変数に読みだす関数
計算はftpでアップロード、ダウンロードして、社内ネットワークなのでtelnetで入ってコマンドを叩く方法としました。これらは全てPythonのモジュールでできます。
出力ファイルを読み込むのは、欲しい数値が固定位置にあるとは限らないので、正規表現で読みだすことにしました。Pythonには正規表現が扱えるモジュールが最初から使えるので便利です。
最適化
設計計算をしていると、例えば出力されたある値が最小(最大)になるような入力値を探したいということがあります。
地道にループで書いてもいいのですが、scipy.optimizeという最適化に使えるモジュールがあるため、これが使えそうな処理であれば使うと楽です。
ログ出力
自動化プログラムを開発するうえで、ログ出力をするとデバッグに便利です。途中でエラーで止まったときにどこに不具合があるか特定しやすいです。
デバッグだけではなく、計算過程を出力しておくことで、あとでなぜその設計になったのかが確認できます。
loggingモジュールがあるので、これを使います。ちょっと使い方にコツがいるので調べてね。
ほかにもデバッグとしては、エラーで止まったときに、「import pdb; pdb.pm()」をやったりとか、あやしそうな箇所で「pdb.set_trace()」したりとかが有効です。
ログの一種としてグラフで経過をわかるようにもしました。matplotlibというグラフを表示するモジュールがあるため、非常に簡単にグラフを表示することができます。
AutoCADの作図の自動化
AutoCADでの作図にはAcadRemoconというフリーソフトを用いました。これはActiveX DLLで、AutoCADに作図コマンドを送ることができます。このActiveX DLLをPythonから呼び出すことで、AutoCADに作図コマンドを送り作図を行うことができます。
私の環境では問題なく動いているのですが、新しいバージョンのWindowsやAutoCADだとサポートされていないようなので、動くかどうかは保証できません。
もしAcadRemoconが使えないような環境であれば、AutoCADの作図コマンドをファイル出力して、ユーザがAutoCADのコマンドウィンドウに貼り付けることで作図できると思います。
作図コマンドで出来ることに限るため、あらゆる作図処理が出来るわけではないと思いますが、線や円弧を作図したり、レイヤーを変えたり、簡単な寸法を入れるぐらいであれば問題なく作図できます。
外部のdxfファイルにある図形を貼り付けるようなこともできます。これはおそらく線分や円弧などのエンティティ毎に貼り付ける操作を実装するしかないはず。少なくとも私はそれ以外は知りません。
このあたりは複雑な処理を行いたければ、AutoCAD.NET等でDLLを開発して、AutoCADからロードすれば出来ますが、資料が少ないのが難点ですね。和書が1冊出てた記憶です。
何が出来たのか
自動化プログラムを開発することで、1週間ほどかかっていた計算が30分ほどで出来るようになりました。
実際のところ、設計自体が複雑でパラメータも多いため、全体を最適化までは出来なかったのですが、全体をワンフローで流し、一部のプロセスに限っては最適化することはできました。
副次的な効果として、以下のようなこともできました。
- 設計ルールが明文化されていない箇所が少なからずあり、ルールを決めることができた(というか決めなければ進められなかった)
- 自動化すると計算ミスが無くなった(もちろん最初の検証作業は必要。検証作業で大量のバグが見つかります)
- 自動化すると人の技量による差がなくなる(ベテランも新人も入力が同じなら同じ答えが出てくる)
なぜ開発がうまくいったのか
この手のプロジェクトは失敗がよくあると思います。
今回は開発がうまくいった例なのですが、一番の理由は設計をよく理解している人間がプログラムを書いたことです。
プログラムを動かして結果やインターフェースがイマイチだったときには、すぐに仕様を変更することができます。
これが外注に頼むと仕様を作りこむだけで数カ月、開発に数カ月、それだけかけてプログラムを動かしてみると期待していた結果が出てこないということになっていたはずです。こうなると追加で開発するために数百万がかかってきます。
設計にもよりますが、今回のような複雑な設計の場合、最初から仕様を固めるのは非常に難しいため、開発しながら考え、柔軟に方針転換ができます。新たなアイディアが思いついたときにも、すぐにプログラムを書き手元で試せます。
結論:設計者にプログラムを書かせよう
おしまい。
素因数分解(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
素因数分解(1) Pollardのρ法
RSA暗号が素因数分解の難しさを安全性の根拠にしていることからわかる通り、大きな数の素因数分解は非常に難しいです。
今のところ、200桁程度の数の素因数分解が限界のようです。
素因数分解(1) Pollardのρ法 - wacchoz’s note(今ここ)
素因数分解(2) p-1法、p+1法、楕円曲線法 - wacchoz’s note
素因数分解(3) フェルマー法、2次ふるい法 - wacchoz’s note
素因数分解アルゴリズムはいくつかの方法が知られています。
ざっくり分けると
の2つのカテゴリーになりますが、今回紹介するPollardのρ法は、どちらにも属しません。
Pollardのρ法は、1975年にポラードが考案した方法です。古典的な方法で、ある程度小さい素因子をもつ場合にしか適用できないものの、他のアルゴリズムとは一風変わった面白い手法です。
ポラード・ロー素因数分解法 - Wikipedia
まず素数を適当に1つとり、
と数列を定義すると、は疑似乱数のように振る舞います。
の数式自体にあまり意味はなく、他の数式でも疑似乱数のように振る舞うものであれば、何でもOKです。
はからの値をとるため、いつか同じ値をとり、その後、周期的に変化することになります。
例として、としてみると、
となり、12項目と5項目が一致しました。この後は周期7で繰り返すことになります。
100000までの素数で、何項目で初めて一致する値が出現するかをプロットしてみます。
プロットしてみると、はからの値をとりうるにもかかわらず、かなり早めに一致する値が出現することがわかります。
pが10^5ぐらいでも、数100項目で一致する値が出現します。
これに似た話は誕生日のパラドックスとしてよく知られています。
ランダムに23人を集めてきたら、その中に同じ誕生日の組がいる確率は50%を超え、70人いれば99.9%にも達します。
一般的には、通りの値をとる乱数があるとき、個の中に同じ値が出現する確率が50%を超えるのは、
で、の値に比較して、かなり少ない個数で一致します。
さて、今までは素数を法とした計算としていたが、素因数分解したい数を法とした計算をしてみます。
として、を計算してみます。
mod 25283の数列について、mod 131で考えると左側の数列と一致するので、9131と3236はmod 131で合同になるはずです。
実際、9131-3236は131で割りきれますし、も131で割り切れるため、
となっています。そこで、
「合成数が与えられたとき、を計算時に、ととの最大公約数を計算し、1でなければ終了。1であれば次のkに進む」
という方法が思いつきます。この方法は素因数をとすると、のオーダーまでにかなりの確率で素因子が見つかりますが、計算量がに比例するため、結局のところ、に比例する計算量となり、試し割り算と変わりません。
もっと効率的に計算するために、フロイドの循環検出法という手法を用います。
フロイドの循環検出法 - Wikipedia
あるから先はとなったとします。つまり数列のある箇所から先が周期になったということです。
このとき、になるようにを大きくとると
となるため、とおくと
が成り立ちます。
ちなみに下図がρ法の名前の由来。
よってとなる可能性が高く、そうなれば素因子が見つかったことになります。
の値はわからないため、小さい順にを計算し、との最大公約数をとります。
具体的には、2系統の数列とを生成し、順にを計算すればよいです。
計算量は平均的にはに比例となるため、でおさえられます。
def gcd(a, b): if a<b: return gcd(b,a) r = b while a%b!=0: r = a%b a,b = b, r return r def f(x): return x*x+1 N = 2**(2**6)+1 print(N) x=2 y=2 for i in range(10**7): x = f(x) % N y = f(f(y)) % N g = gcd(x-y, N) if 1<g<N: print(g) break else: print("Fail")
APR-CL素数判定法
APR-CL素数判定法(Adleman–Pomerance–Rumely-Cohen-Lenstra primality test, APR-CL primality test)について概説&実装してみる。
この素数判定法はAdleman–Pomerance–Rumelyらによって開発された方法で、のちにCohenとLenstraによって改良されています。
数100桁程度の数に対しても現実的な時間で素数判定ができるすごいやつです。
ただし使う数学の道具も高度で、場合分けも多く、実装は重めです。確率的に「ほぼほぼ素数」と判定すればいいだけなら、Miller–Rabin法などを使う方が実装は楽です。
参考文献
- [APR] Adleman, Leonard M.; Pomerance, Carl; Rumely, Robert S. (1983). "On distinguishing prime numbers from composite numbers". Annals of Mathematics. 117 (1): 173–206.
- [Coh-Len1] Cohen, Henri; Lenstra, Hendrik W., Jr. (1984). "Primality testing and Jacobi sums". Mathematics of Computation. 42 (165): 297–330.
- [Coh-Len2] H. Cohen and A. K. Lenstra, Implementation of a new primality test. Math. Comp. 48 (1987), 103-121.
- [和田] 和田秀男, コンピュータと素因子分解, 遊星社, 1987
- [Cohen] H. Cohen, "A course in computational algebraic number theory". Springer-Verlag, Berlin, 1993
- [Chmie] A Chmielowiec, "Primality proving with Gauss and Jacobi sums", 2004.
[APR]はおそらく読まなくても大丈夫。
[Coh-Len1]が[APR]からの改良です。
[Coh-Len2]も実装の上で有用なことが書いています。
和書でAPR法の解説をしているのは[和田]しか見たことがありません。ただ実装できるレベルまでの解説は無いので、実装するなら[Cohen]の本が良いかと思います。
[Chmie]は[Cohen]を簡潔にまとめた解説論文です。
Dirichlet指標とGauss和
まずDirichlet指標とGauss和の導入をします。
乗法群から、(0でない)複素数全体の乗法群への群準同型
を素数を法とするDirichlet指標と呼びます。
つまりmod での乗算構造をそのままに、複素数の乗法に置き換えたもの。
ディリクレ指標 - Wikipedia
具体的に指標の1つを構成します。
奇素数の原始根の1つをとします。とすると
なるがで唯一に定まる。このをと表す。
イメージ的にはのような対数。
をの素因子の1つとし、でありとします。
1の原始乗根を
とすると、
とすればDirichlet指標となっています。
あるいは以下のように書いた方が直感的でわかりやすいかも。
定義から
が成立するのは明らかでしょう。
さてここで、Gauss和を
と定義します。
が成り立つため、Gauss和は
という形に表せ、この形に表したときには一意に決まることが証明できる。
よって
の意味を
と定めてもいいです。
Gauss和を用いた素数判定の原理
具体的に話を進めるために
とする。がの約数となる素数はたくさんある。例えば
であるので、などがの候補となる。
これらのについての積を
と定義する。ここではの素因子分解の中にが何回現れるかという回数である。
のとき、
となる。素数判定したい数をとすると、までの数が素数判定できる。
実装上ははテーブルにしておき、条件を満たすtを探すことになる。
はあらかじめ試し割り算で小さな素因数は無いことがわかっているものとする。つまり、としてよい。
なるを選び、さらにをの素因子の1つとし、でありとします。(最初の方のGauss和の説明と同じ記号にしています。)
このとき、が素数であれば
となる。がからまで動けば、もからまで動くので
が成り立つ。
最後の
がポイントで、もしこの合同式が成立しなければ、は合成数となる。
素数の組合せは多数あるが、が素数であれば、すべての組合せで成立しなければならず、一度でも成立しなければ合成数とわかる。
Fermatの小定理と考え方は同じ。
この合同式がすべてのについて成立しても、残念ながら素数とはいえず、次のステップに進むことになりますが、ひとまずそれは置いておきます。
この計算量を見積もってみます。は100桁ぐらいを想定しておきます。の組合せはいろいろありますが、一番大きいケースでです。
Gauss和を
と表したとき、メモリ量はのオーダーとなる。
Gauss和同士の積は、愚直に計算すると、その2乗となり程度。
の計算は繰り返し2乗法を使えるので、に比例する計算量であり、300回程度の繰り返しとなる。全体としてとなる。さらにはで計算する必要があり、多倍長計算となる。
これを多数のの組合せについて計算するので、不可能ではないレベルだが、実用上は不満のある計算量となる。
これをJacobi和を用いた形に書き換えていく。
Jacobi和への書き換え
Gauss和を用いると計算量が増えるのは、Gauss和はという大きな環に属していたことが原因といえる。
唐突だが、とをでのDirichlet指標とするとき、 Jacobi和を
と定義する。
このJacobi和はという小さな環に属しているので計算量が少ないです。先の例でいえば、が最大である。
をのとき、それ以外のときとなる自明な指標とする。
Jacobi和の性質として、以下が成り立つ。
(1) をのDirichlet指標とすると、
・
・
(2) とをなるのDirichlet指標とすると
・
例えばのとき
のようにJacobi和の積にできたりします。
新たな記号を導入します。
として、をとします。()
とし、の整数係数の和からなる群環を定義します。
と書くとき、
ととするとき、のへの作用を、
で定義します。
少しわかりにくいですが、ただの表記だけの問題で、
のように計算します。
素数判定式を以下のように言い換えます。
「とするとき、が素数であれば、1の乗根が存在し、
となる。」
ここからかなり技巧的に感じる式変形をします。
をなる整数とし、をかつなる整数の集合とするとき、
と
とすると、
と変形でき、Gauss和がJacobi和に書き換えできました。
判定する際には、をとります。ただしという制約のため、のときには場合分けが必要です。
ここから先の詳細は[Cohen]か[Coh-Len1]を参照ということで。
素数判定アルゴリズム
[Cohen]に従い、実際に実装できるレベルでアルゴリズムをまとめておきます。
事前計算
Jacobi和などを事前計算しておきます。
を素数判定したい数の上限とします。
このステップはのみにより、素数判定したい数にはよりません。
- となるを見つける。
- を割りきる素数に対し以下の(a)から(c)を行う。
(a) での原始根を見つける。
に対してとなるようなを計算する。
(b) を割りきる素数に対し、をかつなる整数とする。
Dirichlet指標を、で定義する。
(c) ・またはかつのとき、以下を計算する。
・かつのとき、を上式で計算し、以下を計算する。
素数判定
素数判定したい数を()とする。
Step1. であれば、は合成数である
Step2. に対して、
・かつであれば、
・それ以外であれば、
Step3. となる全ての素数のペアに対し
・のとき、Step 4aへ
・かつのとき、Step 4bへ
・かつのとき、Step 4cへ
・かつのとき、Step 5へ
Step 4a. を未満でと互いに素な自然数の集合とする。
・
・
・
・
・
・
このとき、となるようなの乗根が存在しないとき、は合成数である。
もしが存在し、原始冪根であるとき、とする。
Step 4b. を未満ででまたはに合同な自然数の集合とする。
・
・
・
・
・
・
ここではがかに合同であるときに、それ以外はとする。
このとき、となるようなの乗根が存在しないとき、は合成数である。
もしが存在し、原始冪根であり、かつ、のとき、とする。
Step 4c.
・
・
・のとき、
のとき、
このとき、となるようなの乗根が存在しないとき、は合成数である。
もしが存在し、原始冪根であり、かつ、のとき、とする。
Step 4d.
・
このとき、であれば、は合成数である。
もしかつのとき、とする。
Step 5. なるに対して以下を行う。
・かつとなるをランダムに選ぶ。
・Step 4a-4dをに対して行う。
・何度か繰り返してもとなるがあれば「判定不能」を返す。(現実的にはこれはまず起こらない)
Step 6.
・に対して、を計算する。
・もしあるでがを割りきれば合成数である。そうでなければ素数と判定される。
以上がAPR-CL素数判定アルゴリズムです。
かなり場合分けが煩雑ですね。
上記アルゴリズムだけでは実装には悩む箇所があるので、実装上は[Coh-Len2]を参照するといいです。
を作用させる箇所や、1のべき根であることの判定法などが具体的に記載されています。
Python3による実装
[Cohen]とはできるだけ記号を合わせているので、比較しながら読めばわかるはず。
上のアルゴリズムで事前計算とされている部分は、事前計算にせずにその場で計算しています。
原始根なども保存せずに毎回計算したりしているので、少し工夫すればもう少しは早くなるかも。
手元の環境での実行時間は、100桁だと2秒、2^521-1(157桁)が18秒、2^1279-1(386桁)が5分55秒といったところで、まずまずの速度でしょうか。
GitHub - wacchoz/APR_CL: Adleman–Pomerance–Rumely-Cohen-Lenstra primality test (APR-CL)
# -*- coding: utf-8 -*- import copy import time from math import gcd # version >= 3.5 # primality test by trial division def isprime_slow(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 # v_q(t): how many time is t divided by q def v(q, t): ans = 0 while(t % q == 0): ans +=1 t//= q return ans def prime_factorize(n): ret = [] p=2 while p*p <= n: if n%p==0: num = 0 while n%p==0: num+=1 n//= p ret.append((p,num)) p+= 1 if n!=1: ret.append((n,1)) return ret # calculate e(t) def e(t): s = 1 q_list = [] for q in range(2, t+2): if t % (q-1) == 0 and isprime_slow(q): s *= q ** (1+v(q,t)) q_list.append(q) return 2*s, q_list # Jacobi sum class JacobiSum(object): def __init__(self, p, k, q): self.p = p self.k = k self.q = q self.m = (p-1)*p**(k-1) self.pk = p**k self.coef = [0]*self.m # 1 def one(self): self.coef[0] = 1 for i in range(1,self.m): self.coef[i] = 0 return self # product of JacobiSum # jac : JacobiSum def mul(self, jac): m = self.m pk = self.pk j_ret=JacobiSum(self.p, self.k, self.q) for i in range(m): for j in range(m): if (i+j)% pk < m: j_ret.coef[(i+j)% pk] += self.coef[i] * jac.coef[j] else: r = (i+j) % pk - self.p ** (self.k-1) while r>=0: j_ret.coef[r] -= self.coef[i] * jac.coef[j] r-= self.p ** (self.k-1) return j_ret def __mul__(self, right): if type(right) is int: # product with integer j_ret=JacobiSum(self.p, self.k, self.q) for i in range(self.m): j_ret.coef[i] = self.coef[i] * right return j_ret else: # product with JacobiSum return self.mul(right) # power of JacobiSum(x-th power mod n) def modpow(self, x, n): j_ret=JacobiSum(self.p, self.k, self.q) j_ret.coef[0]=1 j_a = copy.deepcopy(self) while x>0: if x%2==1: j_ret = (j_ret * j_a).mod(n) j_a = j_a*j_a j_a.mod(n) x //= 2 return j_ret # applying "mod n" to coefficient of self def mod(self, n): for i in range(self.m): self.coef[i] %= n return self # operate sigma_x # verification for sigma_inv def sigma(self, x): m = self.m pk = self.pk j_ret=JacobiSum(self.p, self.k, self.q) for i in range(m): if (i*x) % pk < m: j_ret.coef[(i*x) % pk] += self.coef[i] else: r = (i*x) % pk - self.p ** (self.k-1) while r>=0: j_ret.coef[r] -= self.coef[i] r-= self.p ** (self.k-1) return j_ret # operate sigma_x^(-1) def sigma_inv(self, x): m = self.m pk = self.pk j_ret=JacobiSum(self.p, self.k, self.q) for i in range(pk): if i<m: if (i*x)%pk < m: j_ret.coef[i] += self.coef[(i*x)%pk] else: r = i - self.p ** (self.k-1) while r>=0: if (i*x)%pk < m: j_ret.coef[r] -= self.coef[(i*x)%pk] r-= self.p ** (self.k-1) return j_ret # Is self p^k-th root of unity (mod N) # if so, return h where self is zeta^h def is_root_of_unity(self, N): m = self.m p = self.p k = self.k # case of zeta^h (h<m) one = 0 for i in range(m): if self.coef[i]==1: one += 1 h = i elif self.coef[i] == 0: continue elif (self.coef[i] - (-1)) %N != 0: return False, None if one == 1: return True, h # case of zeta^h (h>=m) for i in range(m): if self.coef[i]!=0: break r = i % (p**(k-1)) for i in range(m): if i % (p**(k-1)) == r: if (self.coef[i] - (-1))%N != 0: return False, None else: if self.coef[i] !=0: return False, None return True, (p-1)*p**(k-1)+ r # find primitive root def smallest_primitive_root(q): for r in range(2, q): s = set({}) m = 1 for i in range(1, q): m = (m*r) % q s.add(m) if len(s) == q-1: return r return None # error # calculate f_q(x) def calc_f(q): g = smallest_primitive_root(q) m = {} for x in range(1,q-1): m[pow(g,x,q)] = x f = {} for x in range(1,q-1): f[x] = m[ (1-pow(g,x,q))%q ] return f # sum zeta^(a*x+b*f(x)) def calc_J_ab(p, k, q, a, b): j_ret = JacobiSum(p,k,q) f = calc_f(q) for x in range(1,q-1): pk = p**k if (a*x+b*f[x]) % pk < j_ret.m: j_ret.coef[(a*x+b*f[x]) % pk] += 1 else: r = (a*x+b*f[x]) % pk - p**(k-1) while r>=0: j_ret.coef[r] -= 1 r-= p**(k-1) return j_ret # calculate J(p,q)(p>=3 or p,q=2,2) def calc_J(p, k, q): return calc_J_ab(p, k, q, 1, 1) # calculate J_3(q)(p=2 and k>=3) def calc_J3(p, k, q): j2q = calc_J(p, k, q) j21 = calc_J_ab(p, k, q, 2, 1) j_ret = j2q * j21 return j_ret # calculate J_2(q)(p=2 and k>=3) def calc_J2(p, k, q): j31 = calc_J_ab(2, 3, q, 3, 1) j_conv = JacobiSum(p, k, q) for i in range(j31.m): j_conv.coef[i*(p**k)//8] = j31.coef[i] j_ret = j_conv * j_conv return j_ret # in case of p>=3 def APRtest_step4a(p, k, q, N): print("Step 4a. (p^k, q = {0}^{1}, {2})".format(p,k,q)) J = calc_J(p, k, q) # initialize s1=1 s1 = JacobiSum(p,k,q).one() # J^Theta for x in range(p**k): if x % p == 0: continue t = J.sigma_inv(x) t = t.modpow(x, N) s1 = s1 * t s1.mod(N) # r = N mod p^k r = N % (p**k) # s2 = s1 ^ (N/p^k) s2 = s1.modpow(N//(p**k), N) # J^alpha J_alpha = JacobiSum(p,k,q).one() for x in range(p**k): if x % p == 0: continue t = J.sigma_inv(x) t = t.modpow((r*x)//(p**k), N) J_alpha = J_alpha * t J_alpha.mod(N) # S = s2 * J_alpha S = (s2 * J_alpha).mod(N) # Is S root of unity exist, h = S.is_root_of_unity(N) if not exist: # composite! return False, None else: # possible prime if h%p!=0: l_p = 1 else: l_p = 0 return True, l_p # in case of p=2 and k>=3 def APRtest_step4b(p, k, q, N): print("Step 4b. (p^k, q = {0}^{1}, {2})".format(p,k,q)) J = calc_J3(p, k, q) # initialize s1=1 s1 = JacobiSum(p,k,q).one() # J3^Theta for x in range(p**k): if x % 8 not in [1,3]: continue t = J.sigma_inv(x) t = t.modpow(x, N) s1 = s1 * t s1.mod(N) # r = N mod p^k r = N % (p**k) # s2 = s1 ^ (N/p^k) s2 = s1.modpow(N//(p**k), N) # J3^alpha J_alpha = JacobiSum(p,k,q).one() for x in range(p**k): if x % 8 not in [1,3]: continue t = J.sigma_inv(x) t = t.modpow((r*x)//(p**k), N) J_alpha = J_alpha * t J_alpha.mod(N) # S = s2 * J_alpha * J2^delta if N%8 in [1,3]: S = (s2 * J_alpha ).mod(N) else: J2_delta = calc_J2(p,k,q) S = (s2 * J_alpha * J2_delta).mod(N) # Is S root of unity exist, h = S.is_root_of_unity(N) if not exist: # composite return False, None else: # possible prime if h%p!=0 and (pow(q,(N-1)//2,N) + 1)%N==0: l_p = 1 else: l_p = 0 return True, l_p # in case of p=2 and k=2 def APRtest_step4c(p, k, q, N): print("Step 4c. (p^k, q = {0}^{1}, {2})".format(p,k,q)) J2q = calc_J(p, k, q) # s1 = J(2,q)^2 * q (mod N) s1 = (J2q * J2q * q).mod(N) # s2 = s1 ^ (N/4) s2 = s1.modpow(N//4, N) if N%4 == 1: S = s2 elif N%4 == 3: S = (s2 * J2q * J2q).mod(N) else: print("Error") # Is S root of unity exist, h = S.is_root_of_unity(N) if not exist: # composite return False, None else: # possible prime if h%p!=0 and (pow(q,(N-1)//2,N) + 1)%N==0: l_p = 1 else: l_p = 0 return True, l_p # in case of p=2 and k=1 def APRtest_step4d(p, k, q, N): print("Step 4d. (p^k, q = {0}^{1}, {2})".format(p,k,q)) S2q = pow(-q, (N-1)//2, N) if (S2q-1)%N != 0 and (S2q+1)%N != 0: # composite return False, None else: # possible prime if (S2q + 1)%N == 0 and (N-1)%4==0: l_p=1 else: l_p=0 return True, l_p # Step 4 def APRtest_step4(p, k, q, N): if p>=3: result, l_p = APRtest_step4a(p, k, q, N) elif p==2 and k>=3: result, l_p = APRtest_step4b(p, k, q, N) elif p==2 and k==2: result, l_p = APRtest_step4c(p, k, q, N) elif p==2 and k==1: result, l_p = APRtest_step4d(p, k, q, N) else: print("error") if not result: print("Composite") return result, l_p def APRtest(N): t_list = [ 2, 12, 60, 180, 840, 1260, 1680, 2520, 5040, 15120, 55440, 110880, 720720, 1441440, 4324320, 24504480, 73513440 ] print("N=", N) if N<=3: print("input should be greater than 3") return False # Select t for t in t_list: et, q_list = e(t) if N < et*et: break else: print("t not found") return False print("t=", t) print("e(t)=", et, q_list) # Step 1 print("=== Step 1 ===") g = gcd(t*et, N) if g > 1: print("Composite") return False # Step 2 print("=== Step 2 ===") l = {} fac_t = prime_factorize(t) for p, k in fac_t: if p>=3 and pow(N, p-1, p*p)!=1: l[p] = 1 else: l[p] = 0 print("l_p=", l) # Step 3 & Step 4 print("=== Step 3&4 ===") for q in q_list: if q == 2: continue fac = prime_factorize(q-1) for p,k in fac: # Step 4 result, l_p = APRtest_step4(p, k, q, N) if not result: # composite print("Composite") return False elif l_p==1: l[p] = 1 # Step 5 print("=== Step 5 ===") print("l_p=", l) for p, value in l.items(): if value==0: # try other pair of (p,q) print("Try other (p,q). p={}".format(p)) count = 0 i = 1 found = False # try maximum 30 times while count < 30: q = p*i+1 if N%q != 0 and isprime_slow(q) and (q not in q_list): count += 1 k = v(p, q-1) # Step 4 result, l_p = APRtest_step4(p, k, q, N) if not result: # composite print("Composite") return False elif l_p == 1: found = True break i += 1 if not found: print("error in Step 5") return False # Step 6 print("=== Step 6 ===") r = 1 for t in range(t-1): r = (r*N) % et if r!=1 and r!= N and N % r == 0: print("Composite", r) return False # prime!! print("Prime!") return True if __name__ == '__main__': start_time = time.time() APRtest(2**521-1) # 157 digit, 18 sec # APRtest(2**1279-1) # 386 digit, 355 sec end_time = time.time() print(end_time - start_time, "sec")
某女の子が歌ったり踊ったりするゲームの何か(4)
某女の子が歌ったり踊ったりするゲームの何か(3)
歌うよー!踊るよー!
どうせ暇なので昔作ったプログラムの整理の続き。
「マルチバイト文字セットを使用する」でビルドして下さい。
入力ファイルの設定は"settings.ini"を書き換えて下さい。
(2013/5/5差し替え)
(2018/11 GitHubに移動)
前回のモデル表示までだとコードは2000行強でしたが、歌って踊るまで実装すると一気に6000行位。なかなか大変です。
勢いで書いたので、やっつけで書いたコードとか残ってるけど気にしないで。
いろいろバグが残ってます・・・わかっちゃいるけど・・・