複数多項式二次ふるい法（MPQS）

mpqs.rb

#複数多項式二次ふるい法（MPQS）
#    の作りかけ
# 
# 参考：
# https://en.wikipedia.org/wiki/Quadratic_sieve#Example_of_basic_sieve
# http://www.cs.t-kougei.ac.jp/nsim/lecture/2008/ws/QS.pdf
# http://www.asahi-net.or.jp/~KC2H-MSM/mathland/math12/math1207.htm
# http://inaz2.hatenablog.com/entry/2016/01/09/032521
# http://d.hatena.ne.jp/lemniscus/20130226/1361874593#special
# 
# 
# 二次ふるい法は巨大な因数を見つけるアルゴリズム 10^60 程度の数に有効
# 
# 最初からこれ一本で行くのではなくて、
# 先にポラード・ロー法をかけて、小さめの因数をおとしておく
# 
# 

require 'set'
require 'prime'
require './superprime'
require './gauss_mod'



# 素因数分解したいターゲット
#n = 15347   # = 103 * 149        # 4 - 3
#n = 1091 * 1093                   # 10 - 10
#n = (2**19 - 1) * (2**31 - 1)    # 20 - 20
#n = 2**2**6 + 1                   # 100 - 100
#n = 9880973 * 9880991              # 40 - 40
n = 182521213001 * (2**31 - 1)   # 140 - 140
#n = 1325815267337711173 * 549797184491917     # ? - ? 


puts "ターゲット"
puts "n = #{n}"
puts "(#{ Math.log10(n).to_i+1 } digits)"
puts ""


# 使用する基底素数の数
size_of_bases = 140

# 使用する多項式の数
size_of_polynominals = 140



# べき剰余を計算する
def mod_pow(base, exp, mod)
	r = 1
	while exp > 0
		r = (r * base) % mod if exp & 1 == 1
		base = (base * base) % mod
		exp >>= 1;
	end
	r
end



# オイラーの基準を使って，平方剰余記号を計算する
def isQuadratic a, p
	if mod_pow(a, (p-1)/2, p) == 1 
		true
	else
		false
	end
end



# 基底 (n と平方剰余な素数を選ぶ)
def calculate_bases n, size_of_bases
	bases = Array.new

	bases << (-1)

	primes = Prime.each 
	while bases.size < size_of_bases
		p = primes.next
		if isQuadratic(n, p) 
			bases << p
		end
	end

	puts "基底"
	print "B = "
	p bases
	puts ""

	bases
end


bases = calculate_bases(n, size_of_bases)




# sqrt(n) からスタート
x = Math::sqrt(n).to_i

# 多項式の保管用のリスト変数
polynominals = Array.new

# 指数の行列計算用
matrix = []


puts "基底の素因数によって表せる多項式 (B-smooth) を列挙"


def addPolynominal n, x, a, bases, polynominals, matrix
	# mod n の平方数を作る
	y = (x+a)*(x+a) - n

	# y を 基底 で素因数分解 (基底で割り切れない素因数を持つ場合は除外)
	factors = []
	if y < 0 then
		factors << [-1, 1]
	else 
		factors << [-1, 0]
	end
	y_copy = y.abs
	bases.each do |prime|
		if prime > 0 then
			exp = 0
			while y_copy % prime == 0
				exp += 1
				y_copy /= prime
			end
			factors << [prime, exp]
		end
	end


	# factors の指数が mod 2 で全部ゼロの場合は除外
	sum_of_exp = 0
	factors.each do |e|
		sum_of_exp += (e[1] % 2)
	end


	if y.abs != 1 && y_copy.abs == 1 && sum_of_exp > 0
		# 多項式として採用
		polynominals << [(x+a), factors]
		matrix << factors.map{|e| e[1]%2 }

		id = polynominals.size 

		puts "x[#{id}] = #{(x+a)}"
		print "y[#{id}] = #{y} (= #{(x+a)}*#{(x+a)} - #{n}) = "
		if factors.size < 40
			p factors
		else
			puts "[too many factors]"
		end
		puts ""

	end

end


# 多項式を列挙する
a = 1
b = -1
while polynominals.size < size_of_polynominals

	addPolynominal(n, x, a, bases, polynominals, matrix)
	addPolynominal(n, x, b, bases, polynominals, matrix)
	a += 1
	b -= 1
end
puts ""




if size_of_bases < 100
	# 
	puts "行列"

	print "    "
	matrix[0].size.times do |k|
		printf ",%4d", bases[k]
	end
	puts ""


	def print_matrix matrix

		matrix.size.times do |i|
			printf "%3d:", i
			matrix[i].size.times do |k|
				if matrix[i][k] == 1 
					print "    *"
				else
					print "     "
				end
			end
			puts ""
		end
		puts ""

	end

	print_matrix matrix
	#p matrix
end


def trans matrix
	n = matrix.size
	m = matrix[0].size

	result_matrix = Array.new(m).map{ Array.new(n) }

	n.times do |i|
		m.times do |j|
			result_matrix[j][i] = matrix[i][j]
		end
	end

	result_matrix
end



modulo = 2

a_matrix = trans matrix
b_vector = Array.new(size_of_bases, 0)

augmentedMatrix = AugmentedMatrix.new a_matrix, b_vector, modulo

augmentedMatrix.gauss
x_vector = augmentedMatrix.solve

puts ""
p x_vector

dimension = x_vector.size - 1
puts "dimension: #{dimension}" # ベクトル空間の次元
puts ""
puts ""




def to_bin_array number, digit
	number += 2**digit
	array = number.to_s(2).split("").map{|n| n.to_i }.reverse
	array.pop

	array
end




selection = [] # x_vector[dimension] ]
(1..(2**dimension-1)).each do |pattern|

	pattern_array = to_bin_array( pattern, dimension )
	select_vector = Array.new(x_vector[0].size, 0)

	pattern_array.size.times do |i|
		select_vector.size.times do |k|
			select_vector[k] = (select_vector[k] + pattern_array[i]*x_vector[i+1][k]) % modulo
		end
	end

	selection.push select_vector
	#p select_vector
end



selection.each do |s|
	puts "############"
	puts "使用例："
	p s
	puts ""

	
	# x_large
	x_large = 1

	num_of_factors = Array.new(size_of_bases, 0)
	list_of_factors = Array.new(size_of_bases, nil)

	size_of_polynominals.times do |i|
		if s[i] == 1 then
			x_large *= polynominals[i][0]

			factors = polynominals[i][1]
			factors.size.times do |i|
				base = factors[i][0]
				num  = factors[i][1]

				num_of_factors[i] += num
				list_of_factors[i] = base
			end
		end
	end

	# sqrt
	num_of_factors.size.times do |i|
		num_of_factors[i] = num_of_factors[i] / 2
	end

	# y_large
	y_large = 1
	size_of_bases.times do |i|
		y_large *= (list_of_factors[i] ** num_of_factors[i])
	end

	# 約数の計算
	divisor = (x_large-y_large).gcd(n).abs


	# 約数が見つかったとき
	if divisor != 1 && divisor != n 
		puts "X = #{x_large}, Y = #{y_large}"

		puts ""
		puts "gcd(X+Y, n) = #{divisor}"


		puts ""
		puts ""
		puts "* * * RESULT * * *"
		puts "#{n} = #{divisor} * #{n / divisor}"

		puts ""
		puts ""

		break
	end

end