素因数分解(3) フェルマー法、２次ふるい法

以前の記事でもちらっと書きましたが、素因数分解のアルゴリズムは大きなカテゴリーが２つあり、そのうちの１つにフェルマー法、２次ふるい法、連分数法、数体ふるい法が属しています。
今回はフェルマー法、２次ふるい法を紹介していきます。

素因数分解したい数 $N$ が $N=x^2-y^2$ と分解できたとします。そうすれば $N=(x+y)(x-y)$ と分解でき素因数分解できます。そこでこのような $x, y$ を探すことにします。
xの初期値として、 $x=\lfloor \sqrt{N} \rfloor + 1$ とします。
yの初期値は $y=\lfloor \sqrt{x^2-N} \rfloor$ とします。

$w=x^2-N-y^2$ を計算する
$w=0$ なら $x+y, x-y$ が素因数である
$w>0$ なら $y\leftarrow y+1$ とする
$w<0$ なら $x\leftarrow x+1$ とする
1.に戻る

$y$ の初期値が小さい値からインクリメントしていることかわかる通り、 $N$ の2個の素因子の値が近い場合に見つかるのが早いです。

２次ふるい法

２次ふるい法では
$\displaystyle x^2\equiv y^2 \pmod{N}$
となる $x, y$ を見つけることを目標とします。これが見つかれば
$\displaystyle x^2-y^2=(x+y)(x-y)$
が $N=pq$ の倍数になるということなので、うまくいけば $p, q$ が $x+y, x-y$ に分かれるかもしれません。そうなれば $\text{gcd}(x+y, N)$ が素因数として見つかることになります。

今、 $x$ を動かして $x^2$ を $N$ で割った余り $r$ をたくさん考えます。
$x^2 \equiv r \pmod{N}$
です。ここで、いくつかの $r$ をうまく組合せて平方数を作れれば、目的が達成されます。
これは各 $r$ を素因数分解して、各素数の指数を組み合わせて偶数を作ることを意味し、行列にGauss消去法を適用することにより計算できます。

もう少し具体的にすると、 $x=x_1, x_2, \cdots, x_t$ のときの $r_i=x_i{}^2-N$ となったとし、それぞれの $r_i$ を素因数分解すると
$\displaystyle r_i = p_1^{e_{i1}} p_2^{e_{i2}}\cdots p_m^{e_{im}}$ となったとします。
ここで指数を $\text{mod }2$ でベクトルとして並べて
$\displaystyle v_i=(e_{i1}\text{ mod }2, e_{i2}\text{ mod }2, \cdots, e_{im}\text{ mod }2)$
を作ります。
そして $v_1, v_2, \cdots v_m$ を行とする行列を作り、Gauss消去法で0だけの要素からなる行を作ります。Gauss消去法でどのような操作を行ったかがわかるように単位行列も用意し、単位行列にも同じ操作をしていきます。これでどの行を組み合わせることで0だけの要素からなる行ができるかが決定できます。

それらの行に対応する $r_i$ を掛け合わせれば平方数になるので、それを $y^2$ とすると、目標としていた $x^2\equiv y^2 \text{ mod }N$ が得られる。大抵の場合、 $\text{gcd}(x+y, N)$ で素因子が得られます。

平方数の構成法

具体例として $2,3,4,5,6,7,8,9,10,11$ からいくつかを選びだして、積が平方数になるものを探します。
これらの数の素因子は $2,3,5,7,11$ であるので、各値をそれらで素因数分解したときの指数の $\text{mod }2$ をベクトルとして並べ、右側に単位行列を並べます。

$\displaystyle \left( \begin{array}{ccccc|cccccccccc} 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 \\ \end{array} \right)$

この左半分をGauss消去すると以下のようになります。

$\displaystyle \left( \begin{array}{ccccc|cccccccccc} 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 \\ 0 & 0 & 0 & 0 & 0 & 1 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 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 & 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 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ \end{array} \right)$

$3,5,7,8,9$ 行の左半分が $0$ になっており、平方数になったことを示しています。右半分を見るとどのような組合せなのかわかります。
例えば $5$ 行目は $1,2,5$ 列目が1になっているため、 $1,2,5$ 番目の数、 $2,3,6$ を掛け合わせると平方数になるとわかります。

平方数の構成のプログラム例

# -*- 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乗

「ふるい」とは？

さて、問題をもう一度振り返ると、 $x$ を動かして $x^2$ を $N$ で割った余り $r$ をたくさん考え、いくつかの $r$ をうまく組合せて平方数を作るのでした。
$x^2 \equiv r \pmod{N}$

$r$ は小さい方が平方数が作りやすいため、 $x$ は $\sqrt{N}$ の近くで動かすこととしますが、やみくもに動かしても $r$ を素因数分解したときに多数の素数が出てきて効率が悪くなります。
上の例でいうと $11$ は１度しか出てこないため、平方数を作るのに貢献していません。

そこでまずある $m$ を決めて、 $m$ 個の素数 $p_1, p_2, \cdots, p_m$ を決め、これらで分解される $r$ のみを抽出することを考えます。
（ $p$ は $2$ または $N$ が平方剰余となる奇素数）

$x$ を動かしていった際に、 $x^2-N$ の中で素数 $p$ で割り切れるものがどこに現れるかを考えます。
$x^2\equiv N \pmod{p}$ の解を $x\equiv \pm t_p \pmod{p}$ と計算さえすれば、 $p$ 毎に $p$ の倍数が現れることがわかります。

この部分が2次「ふるい」法と呼ばれるゆえんで、エラトステネスの篩のようなプロセスで篩にかけて抽出していきます。

$p_1, p_2, \cdots, p_m$ 以外では割れない $r$ を求めたいので、最初に $w_i=\log{r_i}$ を計算し、 $p$ で割り切れれば $\log{p}$ を引いていき、 $0$ になれば $p_1, p_2, \cdots, p_m$ のみで割り切れるということになります。

小さい素数 $p$ の場合には、 $p^a$ で割り切れることも多いですが、 $p$ が奇素数の場合には、同様に1つの解を計算すれば $p^a$ 毎に割り切れます。

[木田]では素数のべきで割れるケースは無視して、 $\log{}$ がある閾値よりも低くなるときに実際に割り算して確かめる、という方法をとっています。実装上はこの方が楽だと思います。

# -*- 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