結果
| 問題 | No.3127 Multiple of Twin Prime |
| ユーザー |
 Nanashi.
|
| 提出日時 | 2025-04-25 23:09:16 |
| 言語 | Ruby (3.4.1) |
| 結果 |
AC
|
| 実行時間 | 1,472 ms / 2,500 ms |
| コード長 | 58,687 bytes |
| コンパイル時間 | 410 ms |
| コンパイル使用メモリ | 9,088 KB |
| 実行使用メモリ | 28,416 KB |
| 最終ジャッジ日時 | 2025-04-25 23:09:37 |
| 合計ジャッジ時間 | 20,181 ms |
|
ジャッジサーバーID （参考情報） |
judge5 / judge2 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 1 |
| other | AC * 12 |
コンパイルメッセージ
Main.rb:2247: warning: 'frozen_string_literal' is ignored after any tokens Syntax OK
ソースコード
# Kept libraries:
# (none)
# Expanded libraries:
# - faster_prime
# - nanacl/bsearch_right
# - faster_prime/base
# - faster_prime/core_ext
# - faster_prime/version
# - faster_prime/utils
# - faster_prime/primality_test
# - faster_prime/prime_factorization
# - faster_prime/sieve
# Errored libraries:
# (none)
# Removed libraries:
# (none)
#
# ------------------------------------------------------------------------------
# frozen_string_literal: true
# This file is expanded by nanacl.
main = -> do # =================================================================
# require "faster_prime" # (expanded: L2222)
# require "nanacl/bsearch_right" # (expanded: L2227)
in_t = gets.chomp.to_i
twin_prime_multiples = []
primes = FasterPrime.each(10**7).to_a
primes.each_cons(2) { |p1, p2| twin_prime_multiples << p1 * p2 if p2 - p1 == 2 }
in_t.times do
in_n = gets.chomp.to_i
puts (twin_prime_multiples.bsearch_right { _1 <= in_n } || -1)
end
end # --------------------------------------------------------------------------
# === dependencies -------------------------------------------------------------
# == faster_prime/utils from faster_prime/primality_test -----------------------
module FasterPrime
SMALL_PRIMES = []
SMALL_PRIME_TABLE = []
def self.setup_table(n)
n = (n + 59) / 60 * 60
SMALL_PRIMES.clear << 2
n /= 2
SMALL_PRIME_TABLE.replace([true] * n)
i = 1
while i < n
if SMALL_PRIME_TABLE[i]
SMALL_PRIMES << j = i * 2 + 1
k = i + j
while k < n
SMALL_PRIME_TABLE[k] = false
k += j
end
end
i += 1
end
end
setup_table(100000)
module Utils
module_function
FLOAT_BIGNUM = Float::RADIX ** (Float::MANT_DIG - 1)
# Find the largest integer m such that m <= sqrt(n)
def integer_square_root(n)
if n < FLOAT_BIGNUM
Math.sqrt(n).floor
else
# newton method
a, b = n, 1
a, b = a / 2, b * 2 until a <= b
a = b + 1
a, b = b, (b + n / b) / 2 until a <= b
a
end
end
# Find b such that a * b = 1 (mod n)
def mod_inv(a, n)
raise ArgumentError, "not coprime for mod_inv" if a == 0 || a.gcd(n) != 1
# Extended GCD Algorithm:
# Find $$(x, y)$$ such that $$ax + by = gcd(a, b), x > 0 and y > 0$$
x0, x1 = 1, 0
y0, y1 = 0, 1
b = n
while b != 0
q = a / b
a, b = b, a % b
x0, x1 = x1, x0 - x1 * q
y0, y1 = y1, y0 - y1 * q
end
x0 % n
end
# Calculate a^b mod n
def mod_pow(a, b, n)
r = 1
return 1 if b == 0
len = b.bit_length
len.times do |i|
if b[i] == 1
r = (a * r) % n
return r if i == len - 1
end
a = (a * a) % n
end
raise "assert not reached"
end
# Find x such that x^2 = a (mod prime)
# See ``Complexity of finding square roots'' section of wikipedia:
# http://en.wikipedia.org/wiki/Quadratic_residue
def mod_sqrt(a, prime)
return 0 if a == 0
case
when prime == 2
a.odd? ? 1 : 0
when prime % 4 == 3
mod_pow(a, (prime + 1) / 4, prime)
when prime % 8 == 5
# Let v = (2a)^((prime-5)/8) and r = (2av^2 - 1)av.
#
# 2^((prime-1)/2) = -1 because (2|prime) = -1.
# a^((prime-1)/2) = 1 because (a|prime) = 1.
# So, (2a)^((prime-1)/2) = -1.
#
# r = ((2a)^((prime-1)/4) - 1) * a * (2a)^((prime-5)/8)
# r^2 = ((2a)^((prime-1)/2) - 2(2a)^((prime-1)/4) + 1)
# * a^2 * (2a)^((prime-5)/4)
# = -2(2a)^((prime-1)/4) * a^2 * (2a)^((prime-5)/4)
# = -a * (2a)^((prime-1)/2)
# = a
#
# So r is the solution of x^2 = a (mod prime).
v = mod_pow(a * 2, (prime - 5) / 8, prime)
(2 * a * v * v - 1) * a * v % prime
else # prime % 8 == 1
# Use Tonelli–Shanks algorithm, see:
# http://en.wikipedia.org/wiki/Tonelli%E2%80%93Shanks_algorithm
# Calculate (q, s) such that prime - 1 = 2^s * q and q is odd
s, q = 3, prime - 1
s += 1 until q[s] == 1
q >>= s
# Find quadratic nonresidue z.
z = (2...prime).find {|z_| kronecker(z_, prime) != 1 }
# The main part of the algorithm
c = mod_pow(z, q, prime)
x = mod_pow(a, (q - 1) / 2, prime)
r = a * x % prime # a^((q+1)/2)
t = r * x % prime # a^q
until t == 1
i, d = 0, t
i, d = i + 1, d * d % prime while d != 1
b = mod_pow(c, 1 << (s - i - 1), prime)
c = b * b % prime
r = r * b % prime
t = t * c % prime
s = i
end
r
end
end
# p-adic valuation: Find the highest exponent v such that p^v divides n
def valuation(n, p)
e = 0
while n % p == 0
n /= p
e += 1
end
e
end
# Enumerate divisors of n, including n itself
def each_divisor(n)
return enum_for(:each_divisor, n) unless block_given?
t = 1
while t * t < n
yield t if n % t == 0
t += 1
end
yield t if n == t * t
while t > 1
t -= 1
yield n / t if n % t == 0
end
end
# Calculate the number of positive integers less than or equal to n that
# are coprime to n
def totient(n)
n = n.abs
raise "n must not be zero" if n == 0
r = n
PrimeFactorization.prime_factorization(n) do |prime, exp|
r = r * (prime - 1) / prime
end
r
end
# Finds the least primitive root modulo n.
# Algorithm 1.4.4 (Primitive Root).
def primitive_root(n)
t = Utils.totient(n)
ps = PrimeFactorization.prime_factorization(t).to_a
2.upto(t) do |a|
found = true
ps.each do |q,|
if mod_pow(a, t / q, n) == 1
found = false
break
end
end
return a if found
end
nil
end
# Calculate Kronecker symbol (a|b) in O((log a)(log b)), see pp. 29--30,
# Cohen, Henri (1993), A Course in Computational Algebraic Number Theory
def kronecker(a, b)
return a.abs != 1 ? 0 : 1 if b == 0
return 0 if a.even? && b.even?
k = 1
case a & 7
when 1, 7 then b = b / 2 while b.even?
when 3, 5 then b, k = b / 2, -k while b.even?
end
if b < 0
b = -b
k = -k if a < 0
end
while a != 0
case b & 7
when 1, 7 then a = a / 2 while a.even?
when 3, 5 then a, k = a / 2, -k while a.even?
end
k = -k if 2 & a & b == 2
a = a.abs
a, b = b % a, a
a -= b if a > b / 2
end
return b > 1 ? 0 : k
end
end
end
# == faster_prime/utils from faster_prime/prime_factorization ------------------
module FasterPrime
SMALL_PRIMES = []
SMALL_PRIME_TABLE = []
def self.setup_table(n)
n = (n + 59) / 60 * 60
SMALL_PRIMES.clear << 2
n /= 2
SMALL_PRIME_TABLE.replace([true] * n)
i = 1
while i < n
if SMALL_PRIME_TABLE[i]
SMALL_PRIMES << j = i * 2 + 1
k = i + j
while k < n
SMALL_PRIME_TABLE[k] = false
k += j
end
end
i += 1
end
end
setup_table(100000)
module Utils
module_function
FLOAT_BIGNUM = Float::RADIX ** (Float::MANT_DIG - 1)
# Find the largest integer m such that m <= sqrt(n)
def integer_square_root(n)
if n < FLOAT_BIGNUM
Math.sqrt(n).floor
else
# newton method
a, b = n, 1
a, b = a / 2, b * 2 until a <= b
a = b + 1
a, b = b, (b + n / b) / 2 until a <= b
a
end
end
# Find b such that a * b = 1 (mod n)
def mod_inv(a, n)
raise ArgumentError, "not coprime for mod_inv" if a == 0 || a.gcd(n) != 1
# Extended GCD Algorithm:
# Find $$(x, y)$$ such that $$ax + by = gcd(a, b), x > 0 and y > 0$$
x0, x1 = 1, 0
y0, y1 = 0, 1
b = n
while b != 0
q = a / b
a, b = b, a % b
x0, x1 = x1, x0 - x1 * q
y0, y1 = y1, y0 - y1 * q
end
x0 % n
end
# Calculate a^b mod n
def mod_pow(a, b, n)
r = 1
return 1 if b == 0
len = b.bit_length
len.times do |i|
if b[i] == 1
r = (a * r) % n
return r if i == len - 1
end
a = (a * a) % n
end
raise "assert not reached"
end
# Find x such that x^2 = a (mod prime)
# See ``Complexity of finding square roots'' section of wikipedia:
# http://en.wikipedia.org/wiki/Quadratic_residue
def mod_sqrt(a, prime)
return 0 if a == 0
case
when prime == 2
a.odd? ? 1 : 0
when prime % 4 == 3
mod_pow(a, (prime + 1) / 4, prime)
when prime % 8 == 5
# Let v = (2a)^((prime-5)/8) and r = (2av^2 - 1)av.
#
# 2^((prime-1)/2) = -1 because (2|prime) = -1.
# a^((prime-1)/2) = 1 because (a|prime) = 1.
# So, (2a)^((prime-1)/2) = -1.
#
# r = ((2a)^((prime-1)/4) - 1) * a * (2a)^((prime-5)/8)
# r^2 = ((2a)^((prime-1)/2) - 2(2a)^((prime-1)/4) + 1)
# * a^2 * (2a)^((prime-5)/4)
# = -2(2a)^((prime-1)/4) * a^2 * (2a)^((prime-5)/4)
# = -a * (2a)^((prime-1)/2)
# = a
#
# So r is the solution of x^2 = a (mod prime).
v = mod_pow(a * 2, (prime - 5) / 8, prime)
(2 * a * v * v - 1) * a * v % prime
else # prime % 8 == 1
# Use Tonelli–Shanks algorithm, see:
# http://en.wikipedia.org/wiki/Tonelli%E2%80%93Shanks_algorithm
# Calculate (q, s) such that prime - 1 = 2^s * q and q is odd
s, q = 3, prime - 1
s += 1 until q[s] == 1
q >>= s
# Find quadratic nonresidue z.
z = (2...prime).find {|z_| kronecker(z_, prime) != 1 }
# The main part of the algorithm
c = mod_pow(z, q, prime)
x = mod_pow(a, (q - 1) / 2, prime)
r = a * x % prime # a^((q+1)/2)
t = r * x % prime # a^q
until t == 1
i, d = 0, t
i, d = i + 1, d * d % prime while d != 1
b = mod_pow(c, 1 << (s - i - 1), prime)
c = b * b % prime
r = r * b % prime
t = t * c % prime
s = i
end
r
end
end
# p-adic valuation: Find the highest exponent v such that p^v divides n
def valuation(n, p)
e = 0
while n % p == 0
n /= p
e += 1
end
e
end
# Enumerate divisors of n, including n itself
def each_divisor(n)
return enum_for(:each_divisor, n) unless block_given?
t = 1
while t * t < n
yield t if n % t == 0
t += 1
end
yield t if n == t * t
while t > 1
t -= 1
yield n / t if n % t == 0
end
end
# Calculate the number of positive integers less than or equal to n that
# are coprime to n
def totient(n)
n = n.abs
raise "n must not be zero" if n == 0
r = n
PrimeFactorization.prime_factorization(n) do |prime, exp|
r = r * (prime - 1) / prime
end
r
end
# Finds the least primitive root modulo n.
# Algorithm 1.4.4 (Primitive Root).
def primitive_root(n)
t = Utils.totient(n)
ps = PrimeFactorization.prime_factorization(t).to_a
2.upto(t) do |a|
found = true
ps.each do |q,|
if mod_pow(a, t / q, n) == 1
found = false
break
end
end
return a if found
end
nil
end
# Calculate Kronecker symbol (a|b) in O((log a)(log b)), see pp. 29--30,
# Cohen, Henri (1993), A Course in Computational Algebraic Number Theory
def kronecker(a, b)
return a.abs != 1 ? 0 : 1 if b == 0
return 0 if a.even? && b.even?
k = 1
case a & 7
when 1, 7 then b = b / 2 while b.even?
when 3, 5 then b, k = b / 2, -k while b.even?
end
if b < 0
b = -b
k = -k if a < 0
end
while a != 0
case b & 7
when 1, 7 then a = a / 2 while a.even?
when 3, 5 then a, k = a / 2, -k while a.even?
end
k = -k if 2 & a & b == 2
a = a.abs
a, b = b % a, a
a -= b if a > b / 2
end
return b > 1 ? 0 : k
end
end
end
# == faster_prime/utils from faster_prime/sieve --------------------------------
module FasterPrime
SMALL_PRIMES = []
SMALL_PRIME_TABLE = []
def self.setup_table(n)
n = (n + 59) / 60 * 60
SMALL_PRIMES.clear << 2
n /= 2
SMALL_PRIME_TABLE.replace([true] * n)
i = 1
while i < n
if SMALL_PRIME_TABLE[i]
SMALL_PRIMES << j = i * 2 + 1
k = i + j
while k < n
SMALL_PRIME_TABLE[k] = false
k += j
end
end
i += 1
end
end
setup_table(100000)
module Utils
module_function
FLOAT_BIGNUM = Float::RADIX ** (Float::MANT_DIG - 1)
# Find the largest integer m such that m <= sqrt(n)
def integer_square_root(n)
if n < FLOAT_BIGNUM
Math.sqrt(n).floor
else
# newton method
a, b = n, 1
a, b = a / 2, b * 2 until a <= b
a = b + 1
a, b = b, (b + n / b) / 2 until a <= b
a
end
end
# Find b such that a * b = 1 (mod n)
def mod_inv(a, n)
raise ArgumentError, "not coprime for mod_inv" if a == 0 || a.gcd(n) != 1
# Extended GCD Algorithm:
# Find $$(x, y)$$ such that $$ax + by = gcd(a, b), x > 0 and y > 0$$
x0, x1 = 1, 0
y0, y1 = 0, 1
b = n
while b != 0
q = a / b
a, b = b, a % b
x0, x1 = x1, x0 - x1 * q
y0, y1 = y1, y0 - y1 * q
end
x0 % n
end
# Calculate a^b mod n
def mod_pow(a, b, n)
r = 1
return 1 if b == 0
len = b.bit_length
len.times do |i|
if b[i] == 1
r = (a * r) % n
return r if i == len - 1
end
a = (a * a) % n
end
raise "assert not reached"
end
# Find x such that x^2 = a (mod prime)
# See ``Complexity of finding square roots'' section of wikipedia:
# http://en.wikipedia.org/wiki/Quadratic_residue
def mod_sqrt(a, prime)
return 0 if a == 0
case
when prime == 2
a.odd? ? 1 : 0
when prime % 4 == 3
mod_pow(a, (prime + 1) / 4, prime)
when prime % 8 == 5
# Let v = (2a)^((prime-5)/8) and r = (2av^2 - 1)av.
#
# 2^((prime-1)/2) = -1 because (2|prime) = -1.
# a^((prime-1)/2) = 1 because (a|prime) = 1.
# So, (2a)^((prime-1)/2) = -1.
#
# r = ((2a)^((prime-1)/4) - 1) * a * (2a)^((prime-5)/8)
# r^2 = ((2a)^((prime-1)/2) - 2(2a)^((prime-1)/4) + 1)
# * a^2 * (2a)^((prime-5)/4)
# = -2(2a)^((prime-1)/4) * a^2 * (2a)^((prime-5)/4)
# = -a * (2a)^((prime-1)/2)
# = a
#
# So r is the solution of x^2 = a (mod prime).
v = mod_pow(a * 2, (prime - 5) / 8, prime)
(2 * a * v * v - 1) * a * v % prime
else # prime % 8 == 1
# Use Tonelli–Shanks algorithm, see:
# http://en.wikipedia.org/wiki/Tonelli%E2%80%93Shanks_algorithm
# Calculate (q, s) such that prime - 1 = 2^s * q and q is odd
s, q = 3, prime - 1
s += 1 until q[s] == 1
q >>= s
# Find quadratic nonresidue z.
z = (2...prime).find {|z_| kronecker(z_, prime) != 1 }
# The main part of the algorithm
c = mod_pow(z, q, prime)
x = mod_pow(a, (q - 1) / 2, prime)
r = a * x % prime # a^((q+1)/2)
t = r * x % prime # a^q
until t == 1
i, d = 0, t
i, d = i + 1, d * d % prime while d != 1
b = mod_pow(c, 1 << (s - i - 1), prime)
c = b * b % prime
r = r * b % prime
t = t * c % prime
s = i
end
r
end
end
# p-adic valuation: Find the highest exponent v such that p^v divides n
def valuation(n, p)
e = 0
while n % p == 0
n /= p
e += 1
end
e
end
# Enumerate divisors of n, including n itself
def each_divisor(n)
return enum_for(:each_divisor, n) unless block_given?
t = 1
while t * t < n
yield t if n % t == 0
t += 1
end
yield t if n == t * t
while t > 1
t -= 1
yield n / t if n % t == 0
end
end
# Calculate the number of positive integers less than or equal to n that
# are coprime to n
def totient(n)
n = n.abs
raise "n must not be zero" if n == 0
r = n
PrimeFactorization.prime_factorization(n) do |prime, exp|
r = r * (prime - 1) / prime
end
r
end
# Finds the least primitive root modulo n.
# Algorithm 1.4.4 (Primitive Root).
def primitive_root(n)
t = Utils.totient(n)
ps = PrimeFactorization.prime_factorization(t).to_a
2.upto(t) do |a|
found = true
ps.each do |q,|
if mod_pow(a, t / q, n) == 1
found = false
break
end
end
return a if found
end
nil
end
# Calculate Kronecker symbol (a|b) in O((log a)(log b)), see pp. 29--30,
# Cohen, Henri (1993), A Course in Computational Algebraic Number Theory
def kronecker(a, b)
return a.abs != 1 ? 0 : 1 if b == 0
return 0 if a.even? && b.even?
k = 1
case a & 7
when 1, 7 then b = b / 2 while b.even?
when 3, 5 then b, k = b / 2, -k while b.even?
end
if b < 0
b = -b
k = -k if a < 0
end
while a != 0
case b & 7
when 1, 7 then a = a / 2 while a.even?
when 3, 5 then a, k = a / 2, -k while a.even?
end
k = -k if 2 & a & b == 2
a = a.abs
a, b = b % a, a
a -= b if a > b / 2
end
return b > 1 ? 0 : k
end
end
end
# == faster_prime/version from faster_prime/base -------------------------------
module FasterPrime
VERSION = "1.0.1"
end
# == faster_prime/utils from faster_prime/base ---------------------------------
module FasterPrime
SMALL_PRIMES = []
SMALL_PRIME_TABLE = []
def self.setup_table(n)
n = (n + 59) / 60 * 60
SMALL_PRIMES.clear << 2
n /= 2
SMALL_PRIME_TABLE.replace([true] * n)
i = 1
while i < n
if SMALL_PRIME_TABLE[i]
SMALL_PRIMES << j = i * 2 + 1
k = i + j
while k < n
SMALL_PRIME_TABLE[k] = false
k += j
end
end
i += 1
end
end
setup_table(100000)
module Utils
module_function
FLOAT_BIGNUM = Float::RADIX ** (Float::MANT_DIG - 1)
# Find the largest integer m such that m <= sqrt(n)
def integer_square_root(n)
if n < FLOAT_BIGNUM
Math.sqrt(n).floor
else
# newton method
a, b = n, 1
a, b = a / 2, b * 2 until a <= b
a = b + 1
a, b = b, (b + n / b) / 2 until a <= b
a
end
end
# Find b such that a * b = 1 (mod n)
def mod_inv(a, n)
raise ArgumentError, "not coprime for mod_inv" if a == 0 || a.gcd(n) != 1
# Extended GCD Algorithm:
# Find $$(x, y)$$ such that $$ax + by = gcd(a, b), x > 0 and y > 0$$
x0, x1 = 1, 0
y0, y1 = 0, 1
b = n
while b != 0
q = a / b
a, b = b, a % b
x0, x1 = x1, x0 - x1 * q
y0, y1 = y1, y0 - y1 * q
end
x0 % n
end
# Calculate a^b mod n
def mod_pow(a, b, n)
r = 1
return 1 if b == 0
len = b.bit_length
len.times do |i|
if b[i] == 1
r = (a * r) % n
return r if i == len - 1
end
a = (a * a) % n
end
raise "assert not reached"
end
# Find x such that x^2 = a (mod prime)
# See ``Complexity of finding square roots'' section of wikipedia:
# http://en.wikipedia.org/wiki/Quadratic_residue
def mod_sqrt(a, prime)
return 0 if a == 0
case
when prime == 2
a.odd? ? 1 : 0
when prime % 4 == 3
mod_pow(a, (prime + 1) / 4, prime)
when prime % 8 == 5
# Let v = (2a)^((prime-5)/8) and r = (2av^2 - 1)av.
#
# 2^((prime-1)/2) = -1 because (2|prime) = -1.
# a^((prime-1)/2) = 1 because (a|prime) = 1.
# So, (2a)^((prime-1)/2) = -1.
#
# r = ((2a)^((prime-1)/4) - 1) * a * (2a)^((prime-5)/8)
# r^2 = ((2a)^((prime-1)/2) - 2(2a)^((prime-1)/4) + 1)
# * a^2 * (2a)^((prime-5)/4)
# = -2(2a)^((prime-1)/4) * a^2 * (2a)^((prime-5)/4)
# = -a * (2a)^((prime-1)/2)
# = a
#
# So r is the solution of x^2 = a (mod prime).
v = mod_pow(a * 2, (prime - 5) / 8, prime)
(2 * a * v * v - 1) * a * v % prime
else # prime % 8 == 1
# Use Tonelli–Shanks algorithm, see:
# http://en.wikipedia.org/wiki/Tonelli%E2%80%93Shanks_algorithm
# Calculate (q, s) such that prime - 1 = 2^s * q and q is odd
s, q = 3, prime - 1
s += 1 until q[s] == 1
q >>= s
# Find quadratic nonresidue z.
z = (2...prime).find {|z_| kronecker(z_, prime) != 1 }
# The main part of the algorithm
c = mod_pow(z, q, prime)
x = mod_pow(a, (q - 1) / 2, prime)
r = a * x % prime # a^((q+1)/2)
t = r * x % prime # a^q
until t == 1
i, d = 0, t
i, d = i + 1, d * d % prime while d != 1
b = mod_pow(c, 1 << (s - i - 1), prime)
c = b * b % prime
r = r * b % prime
t = t * c % prime
s = i
end
r
end
end
# p-adic valuation: Find the highest exponent v such that p^v divides n
def valuation(n, p)
e = 0
while n % p == 0
n /= p
e += 1
end
e
end
# Enumerate divisors of n, including n itself
def each_divisor(n)
return enum_for(:each_divisor, n) unless block_given?
t = 1
while t * t < n
yield t if n % t == 0
t += 1
end
yield t if n == t * t
while t > 1
t -= 1
yield n / t if n % t == 0
end
end
# Calculate the number of positive integers less than or equal to n that
# are coprime to n
def totient(n)
n = n.abs
raise "n must not be zero" if n == 0
r = n
PrimeFactorization.prime_factorization(n) do |prime, exp|
r = r * (prime - 1) / prime
end
r
end
# Finds the least primitive root modulo n.
# Algorithm 1.4.4 (Primitive Root).
def primitive_root(n)
t = Utils.totient(n)
ps = PrimeFactorization.prime_factorization(t).to_a
2.upto(t) do |a|
found = true
ps.each do |q,|
if mod_pow(a, t / q, n) == 1
found = false
break
end
end
return a if found
end
nil
end
# Calculate Kronecker symbol (a|b) in O((log a)(log b)), see pp. 29--30,
# Cohen, Henri (1993), A Course in Computational Algebraic Number Theory
def kronecker(a, b)
return a.abs != 1 ? 0 : 1 if b == 0
return 0 if a.even? && b.even?
k = 1
case a & 7
when 1, 7 then b = b / 2 while b.even?
when 3, 5 then b, k = b / 2, -k while b.even?
end
if b < 0
b = -b
k = -k if a < 0
end
while a != 0
case b & 7
when 1, 7 then a = a / 2 while a.even?
when 3, 5 then a, k = a / 2, -k while a.even?
end
k = -k if 2 & a & b == 2
a = a.abs
a, b = b % a, a
a -= b if a > b / 2
end
return b > 1 ? 0 : k
end
end
end
# == faster_prime/primality_test from faster_prime/base ------------------------
# require "faster_prime/utils" # (expanded: L22)
module FasterPrime
module PrimalityTest
module_function
DEFAULT_RANDOM = Random.new
def prime?(n, random = DEFAULT_RANDOM)
r = probable_prime?(n, random)
return r if r != nil
return APRCL.prime?(n)
end
# Miller-Rabin primality test
# Returns true for a prime, false for a composite, nil if unknown.
def probable_prime?(n, random = Random.new)
return n >= 2 if n <= 3
return true if n == 2 || n == 3
return true if n == 5
return false unless 30.gcd(n) == 1
return SMALL_PRIME_TABLE[n / 2] if n / 2 < SMALL_PRIME_TABLE.size
m = n % 6
return false if m != 1 && m != 5
list = DETERMINISTIC_MR_TEST_TABLE.find {|t,| n < t }
if list
ret = true
list = list.last
else
ret = nil
list = (1..20).map { random.rand(n - 2) + 1 }
end
d, s = n - 1, 0
while d.even?
d >>= 1
s += 1
end
list.each do |a|
y = Utils.mod_pow(a, d, n)
if y != 1
f = s.times do
break if y == n - 1
y = (y * y) % n
end
return false if f
end
end
ret
end
# These numbers are taken from:
# http://code.google.com/p/sympy/
# http://primes.utm.edu/prove/prove2_3.html
DETERMINISTIC_MR_TEST_TABLE = [
[ 1_373_653, [2, 3]],
[ 170_584_961, [350, 3_958_281_543]],
[ 4_759_123_141, [2, 7, 61]],
[ 75_792_980_677, [2, 379_215, 457_083_754]],
[ 1_000_000_000_000, [2, 13, 23, 1_662_803]],
[ 2_152_302_898_747, [2, 3, 5, 7, 11]],
[ 3_474_749_660_383, [2, 3, 5, 7, 11, 13]],
[341_550_071_728_321, [2, 3, 5, 7, 11, 13, 17]],
]
# APRCL Primarity Test
#
# Henry Cohen. A course in computational algebraic number theory.
# 9.1 The Jacobi Sum Test
class APRCL
# usage:
#
# (1) simplest
#
# APRCL.prime?(N) #=> true or false
#
# (2) reuse-table
#
# aprcl = APRCL.new(B)
# aprcl.prime?(N1) #=> true or false
# aprcl.prime?(N2) #=> true or false
# aprcl.bound #=> upper bound (>= B)
#
# (3) manual-t (for test or debug)
#
# aprcl = APRCP.new(B)
# aprcl.set_t(t)
# aprcl.prime?(N) #=> true or false
# An internal exception for terminating the algorithm because n is composite
class Composite < StandardError
end
# An exception representing that the test has failed (highly improbable)
class Failed < StandardError
end
# Returns true if n is prime, false otherwise.
def self.prime?(n)
new(n).prime?
end
# Creates an APRCL instance to test an integer that is less than bound
def initialize(n, jacobi_sum_table: DefaultJacobiSumTable, t: nil)
# an integer in question
@n = n
# precompute n-1 and (n-1)/2 because n is normally big
@n1 = n - 1
@n1_2 = @n1 / 2
# @eta[pk]: a set of p^k-th roots of unity
@eta = {}
# a cache of Jacobi sums
@js = jacobi_sum_table
# compute e(t)
if t
raise "t must be even" if t.odd?
@t = t
@et = compute_e(@t)
raise "t is too small to test n" if @n >= @et ** 2
else
@t = find_t(@n)
@et = compute_e(@t)
end
end
# Algorithm 9.1.28 (Jacobi Sum Primality test)
def prime?
return false if @n <= 1
# 1. [Check GCD]
g = @n.gcd(@t * @et)
return g == @n && PrimalityTest.prime?(g) if g != 1
# 2. [Initialize]
lp = {}
PrimeFactorization.prime_factorization(@t) do |p, |
lp[p] = false if p == 2 || Utils.mod_pow(@n, p - 1, p ** 2) == 1
end
# 3. [Loop on characters]
Utils.each_divisor(@t) do |d|
q = d + 1
next unless PrimalityTest.prime?(q)
PrimeFactorization.prime_factorization(d) do |p, k|
# 4. [Check (*beta)]
lp.delete(p) if step4(q, p, k)
end
end
# 5. [Check conditions Lp]
lp.keys.each {|p| step5(p) }
# 6. [Final trial division]
r = 1
1.upto(@t - 1) do
r = r * @n % @et
return false if @n % r == 0 && 1 < r && r < @n
end
return true
rescue Composite
return false
end
private
def step4(q, p, k)
case
when p >= 3 then step4a(q, p, k)
when k >= 3 then step4b(q, k)
when k == 2 then step4c(q)
when k == 1 then step4d(q)
end
end
# p is odd
def step4a(q, p, k)
pk = p ** k
# s1 = J(p,q)^theta, s3 = J(p,q)^alpha
s1 = s3 = ZZeta.one(p, k)
r = @n % pk
jpq = @js.get_j(p, q)
1.upto(pk - 1) do |x|
next if x % p == 0
j = jpq.substitute(x).reduce_mod!(@n)
s1 = s1 * j.mod_pow(x, @n)
s3 = s3 * j.mod_pow(r * x / pk, @n)
end
s2 = s1.mod_pow(@n / pk, @n)
# S(p,q) = s2 J(p,q)^alpha
s = (s2 * s3).reduce_mod!(@n)
i = check_root_of_unity(p, k, pk, s)
i % p != 0 # it is a primitive root of unity
end
# p = 2, k >= 3
def step4b(q, k)
pk = 2 ** k
# s1 = J_3(q)^theta, s3 = J_3(q)^alpha
s1 = s3 = ZZeta.one(2, k)
r = @n % pk
j3, j2 = @js.get_j3_j2(q)
diff = 1
0.step(pk - 1, 4) do |x|
x, diff = x + diff, -diff # x = 8m+1 or 8m+3, i.e., 1, 3, 9, 11, ...
j = j3.substitute(x).reduce_mod!(@n)
s1 = (s1 * j.mod_pow(x , @n)).reduce_mod!(@n)
s3 = (s3 * j.mod_pow(r * x / pk, @n)).reduce_mod!(@n)
end
# S(2,q) = s2 J_3(q)^alpha if r = 8m+1 or 8m+3
# s2 J_3(q)^alpha J_2(q) if r = 8m+5 or 8m+7
s = s3 * s1.mod_pow(@n / pk, @n)
s = s * j2 if r[2] == 1 # r = 8m+5 or 8m+7
s.reduce_mod!(@n)
i = check_root_of_unity(2, k, pk, s)
i[0] == 1 && Utils.mod_pow(q, @n1_2, @n) == @n1
end
# p = 2, k = 2
def step4c(q)
s0 = @js.get_j(2, q).sqr
s1 = s0.scalar(q)
s2 = s1.mod_pow(@n / 4, @n)
# S(2,q) = s2 if n = 4m+1
# s2 J(2,q)^2 if n = 4m+3
s = (@n[1] == 0 ? s2 : s2 * s0).reduce_mod!(@n)
i = check_root_of_unity(2, 2, 4, s)
i[0] != 1 && Utils.mod_pow(q, @n1_2, @n) == @n1
end
# p = 2, k = 1
def step4d(q)
# S(2,q) = (-q)^((n-1)/2) mod N
s = Utils.mod_pow(-q, @n1_2, @n)
raise Composite if s != 1 && s != @n1
s == @n && @n[1] == 0
end
# check whether there exists a p^k-th root of unity e such that s = e mod n
def check_root_of_unity(p, k, pk, s)
unless @eta[pk]
table = {}
pk.times do |i|
table[ZZeta.monomial(p, k, i).reduce_mod!(@n)] = i
end
@eta[pk] = table
end
@eta[pk][s] || raise(Composite)
end
TRIAL = 2000
def step5(p)
bound_k = 4
# Find q and k by brute force.
# We first try to find a small k because step4 might take long time
# when k is big (>= 4). (An example case: prime?(19541) with TRIAL = 128)
# If we cannot do so, we use any k.
2.times do
q = p + 1
TRIAL.times do
if @et % q != 0 && PrimalityTest.prime?(q)
if @n % q == 0
return if @n == q
raise Composite
end
k = Utils.valuation(q - 1, p)
return if k <= bound_k && step4(q, p, k)
end
q += p
end
bound_k = Float::INFINITY
end
raise Failed, "APRCL primality test failed (highly improbable)"
end
# The pairs of t candidate and floor(log_2(e^2(t))).
# generated by tool/build-aprcl-t-table.rb
T_TABLE = [
[12, 31], [24, 33], [30, 34], [36, 54], [60, 65], [72, 68], [108, 70],
[120, 78], [180, 102], [240, 104], [360, 127], [420, 136], [540, 153],
[840, 165], [1008, 169], [1080, 178], [1200, 192], [1260, 206],
[1620, 211], [1680, 222], [2016, 225], [2160, 244], [2520, 270],
[3360, 279], [3780, 293], [5040, 346], [6480, 348], [7560, 383],
[8400, 396], [10080, 426], [12600, 458], [15120, 527], [25200, 595],
[30240, 636], [42840, 672], [45360, 684], [55440, 708], [60480, 771],
[75600, 775], [85680, 859], [100800, 893], [110880, 912], [128520, 966],
[131040, 1009], [166320, 1042], [196560, 1124], [257040, 1251],
[332640, 1375], [393120, 1431], [514080, 1483], [655200, 1546],
[665280, 1585], [786240, 1661], [831600, 1667], [917280, 1677],
[982800, 1728], [1081080, 1747], [1179360, 1773], [1285200, 1810],
[1310400, 1924], [1441440, 2001], [1663200, 2096], [1965600, 2166],
[2162160, 2321], [2751840, 2368], [2827440, 2377], [3326400, 2514],
[3341520, 2588], [3603600, 2636], [3931200, 2667], [4324320, 3028],
[5654880, 3045], [6652800, 3080], [6683040, 3121], [7207200, 3283],
[8648640, 3514], [10810800, 3725], [12972960, 3817], [14414400, 3976],
[18378360, 3980], [21621600, 4761], [36756720, 5067], [43243200, 5657],
[64864800, 5959], [73513440, 6423], [86486400, 6497],
]
# Find t such that e(t) > n^(1/2).
def find_t(n)
n_log2 = n.bit_length
# e^2(t) >= 2^(et2_log2) >= 2^n_log2 > n
t, = T_TABLE.find {|_t, et2_log2| et2_log2 >= n_log2 }
raise "too large" unless t
t
end
# Compute e(t) (the definition comes from Theorem 9.1.12)
def compute_e(t)
r = 2
Utils.each_divisor(t) do |d|
q = d + 1
r *= q ** (Utils.valuation(t, q) + 1) if PrimalityTest.prime?(q)
end
r
end
class JacobiSumTableClass
# Algorithm 9.1.27 (Precomputations)
# Note that this computes the table lazily, i.e., on demand
def initialize
@f = {} # a table of the function f(x) in 2. (1)
@j = {} # J(p,q)
@j3_j2 = {} # J3(p,q) and J2(p,q)
end
def clear
@f.clear
@j.clear
@j3_j2.clear
end
# compute a table of the function f(x) such that 1 - g^x = g^f(x),
# where g is a primitive root modulo q
def calc_f(q)
return @f[q] if @f[q]
g = Utils.primitive_root(q)
f, h = [], {}
1.upto(q - 2) {|x| h[(1 - Utils.mod_pow(g, x, q)) % q] = x }
1.upto(q - 2) {|fx| f[h[Utils.mod_pow(g, fx, q)]] = fx }
@f[q] = f
end
# compute J(p,q)
def get_j(p, q)
# assumes that p >= 3 or (p == 2 && k >= 2)
return @j[[p, q]] if @j[[p, q]]
f = calc_f(q)
k = Utils.valuation(q - 1, p)
pk = p ** k
# J(p,q) = \Sum_{1<=x<=q-2} (zeta_{p^k}^{x+f(x)})
coefs = [0] * pk
1.upto(q - 2) {|x| coefs[(x + f[x]) % pk] += 1 }
@j[[p, q]] = ZZeta.new(coefs, p, k).reduce!
end
# compute J_3(q) and J_2(q)
def get_j3_j2(q)
# assumes that p == 2 && k >= 3
return @j3_j2[q] if @j3_j2[q]
f = calc_f(q)
k = Utils.valuation(q - 1, 2)
pk = 2 ** k
# coefs3: \Sum_{1<=x<=q-2} (zeta_{2^k}^{2x+f(x)})
# coefs2: \Sum_{1<=x<=q-2} (zeta_{2^3}^{3x+f(x)})
coefs3 = [0] * pk
coefs2 = [0] * pk
1.upto(q - 2) do |x|
coefs3[(2 * x + f[x]) % pk] += 1
coefs2[(3 * x + f[x]) % 8 * (pk / 8)] += 1
end
# J_3(q) = J(2,q) coefs3, J_2(q) = ( coefs2 )^2
j3 = (get_j(2, q) * ZZeta.new(coefs3, 2, k)).reduce!
j2 = ZZeta.new(coefs2, 2, k).reduce!.sqr.reduce!
@j3_j2[q] = [j3, j2]
end
end
DefaultJacobiSumTable = JacobiSumTableClass.new
# A group algebra Z[zeta_{p^k}]
# (linear sum of p^k-th roots of unity)
class ZZeta
def initialize(coefs, p, k)
@coefs = coefs
@p, @k = p, k
@pk1 = p ** (k - 1) # p^(k-1)
@pk = @pk1 * p # p^k
@p1pk1 = @pk - @pk1 # (p-1) p^(k-1)
if @pk != @coefs.size
raise "coefs.size (#{ coefs.size }) must be equal to #{ @p }**#{ @k }"
end
end
attr_reader :coefs
protected :coefs
# generates 1
def self.one(p, k)
@@one[p ** k] ||= monomial(p, k, 0)
end
@@one = {}
# generates zeta_{p^k}^i
def self.monomial(p, k, i)
coefs = [0] * (p ** k)
coefs[i] = 1
ZZeta.new(coefs, p, k)
end
# for hashing
def ==(other); @coefs == other.coefs; end
def hash; @coefs.hash; end
alias eql? ==
# scalar multiplication
def scalar(n)
ZZeta.new(@coefs.map {|x| x * n }, @p, @k)
end
# polynomial multiplication
def *(other)
coefs = [0] * @pk
other.coefs.each_with_index do |oc, j|
next if oc == 0
@pk.times do |i|
coefs[(i + j) % @pk] += @coefs[i] * oc
end
end
ZZeta.new(coefs, @p, @k)
end
# polynomial power modulo n
def mod_pow(exp, n)
return ZZeta.one(@p, @k) if exp == 0
r = nil
z = reduce_mod!(n)
len = exp.bit_length
len.times do |i|
if exp[i] == 1
r = r ? (z * r).reduce_mod!(n) : z
return r if i == len - 1
end
z = z.sqr.reduce_mod!(n)
end
raise "assert not reached"
end
# polynomial square (self * self)
# assumes self is already reduced
def sqr
coefs = [0] * @pk
@p1pk1.times do |j|
oc = @coefs[j]
next if oc == 0
@p1pk1.times do |i|
coefs[(i + j) % @pk] += @coefs[i] * oc
end
end
ZZeta.new(coefs, @p, @k)
end
# generate a reminder divided by CyclotomicPolynomial(p^k)
def reduce!
# CyclotomicPolynomial(n) = \Prod_{d|n} (x^{n/d} - 1)^\moeb{d}
#
# CyclotomicPolynomial(p^k) =
# (x^{p^ k } - 1)^\moeb{p^0} (case d = p^0)
# * (x^{p^{k-1}} - 1)^\moeb{p^1} (case d = p^1)
# ...
# * (x^{p^ 0 } - 1)^\moeb{p^k} (case d = p^k)
# = (x^{p^k} - 1) / (x^{p^{k-1}} - 1)
# = \Sum_{0 <= i < p} (x^{i p^{k-1}})
@p1pk1.times do |i|
@coefs[i] -= @coefs[@p1pk1 + (i % @pk1)]
end
@pk1.times do |i|
@coefs[@p1pk1 + i] = 0
end
self
end
# reduce modulo n
def reduce_mod!(n)
@pk.times do |i|
# I would use centermod if it was built-in...
@coefs[i] = (@coefs[i] - @coefs[@p1pk1 + (i % @pk1)]) % n
end
self
end
# replace zeta with zeta^x
def substitute(x)
coefs = []
@pk.times {|i| coefs << @coefs[(x * i) % @pk] }
ZZeta.new(coefs, @p, @k)
end
def inspect
s = []
(@pk - 1).downto(0) do |i|
c = @coefs[i]
next if c == 0
s << (c > 0 ? ?+ : ?-) +
case
when i == 0 then c.abs.to_s
when c.abs == 1 then i == 1 ? "z" : "z^#{ i }"
else "#{ c.abs }#{ i == 1 ? "z" : "z^#{ i }" }"
end
end
return "0" if s.empty?
s.first[0, 1] = "" if s.first[0] == ?+
s.join
end
end
end
end
end
# == faster_prime/prime_factorization from faster_prime/base -------------------
# require "faster_prime/utils" # (expanded: L22)
module FasterPrime
module PrimeFactorization
module_function
# Factorize an integer
def prime_factorization(n)
return enum_for(:prime_factorization, n) unless block_given?
return if n == 0
if n < 0
yield [-1, 1]
n = -n
end
SMALL_PRIMES.each do |prime|
if n % prime == 0
c = 0
begin
n /= prime
c += 1
end while n % prime == 0
yield [prime, c]
end
if prime * prime > n
yield [n, 1] if n > 1
return
end
return if n == 1
end
if PrimalityTest.prime?(n)
yield [n, 1]
else
d = nil
until d
[PollardRho, MPQS].each do |algo|
begin
d = algo.try_find_factor(n)
rescue Failed
else
break
end
end
end
pe = Hash.new(0)
prime_factorization(n / d) {|p, e| pe[p] += e }
prime_factorization(d) {|p, e| pe[p] += e }
pe.keys.sort.each do |p|
yield [p, pe[p]]
end
end
end
end
class Failed < StandardError
end
# An implementation of Pollard's rho algorithm (Brent's variant)
#
# https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm#Variants
module PollardRho
module_function
def self.try_find_factor(n)
x0 = rand(n)
y = x0
m = Math.log(n, 2).ceil
c = 2
r = q = 1
while true
x = y
r.times { y = (y * y + c) % n }
k = 0
while true
ys = y
[m, r - k].min.times do
y = (y * y + c) % n
q = (q * (x - y)) % n
end
g = q.gcd(n)
k += m
break if k >= r || g > 1
end
r *= 2
break if g > 1
end
if g == n
while true
ys = (ys * ys + c) % n
g = (x - ys).gcd(n)
break if g > 1
end
end
if g == n
raise Failed, "failed to find a factor: #{ n }"
else
g
end
end
end
# An implementation of MPQS factoring algorithm.
#
# R. D. Silverman,
# "The Multiple Polynomial Quadratic Sieve." (1987)
# http://www.ams.org/journals/mcom/1987-48-177/S0025-5718-1987-0866119-8/S0025-5718-1987-0866119-8.pdf
class MPQS
# MPQS parameters
# (N digits, factor base size, sieve interval, and large prime tolerance)
PARAMETERS = [
[24, 100, 250, 1.5 ],
[30, 200, 1_250, 1.5 ],
[36, 400, 1_250, 1.75],
[42, 900, 2_500, 2.0 ],
[48, 1200, 7_000, 2.0 ],
[54, 2000, 12_500, 2.2 ],
[60, 3000, 17_500, 2.4 ],
[66, 4500, 25_000, 2.6 ],
]
SMALL_PRIMES = [2, 3, 5, 7, 11]
# Generates a table for finding k such that kN = 1 mod 4 and that kN is
# rich in small quadratic residues, i.e., (kN | p) = 1 for all p in
# SMALL_PRIMES.
def self.k_table
return @@k_table if defined?(@@k_table)
k_table = { 1 => {}, 3 => {} }
Sieve.each do |p|
next if p <= SMALL_PRIMES.last
state = SMALL_PRIMES.map {|q| Utils.kronecker(p, q) }
k_table[p % 4][state] ||= p
break if k_table.all? {|r, t| t.size == 2 ** SMALL_PRIMES.size }
end
@@k_table = k_table
end
def self.try_find_factor(n)
new(n).try_find_factor
end
S = 65536
def initialize(n)
@n = n
if @n <= 1 || PrimalityTest.prime?(@n)
@trivial_factor = 1
return
end
SMALL_PRIMES.each do |p|
if @n % p == 0
@trivial_factor = p
return
end
end
m = Utils.integer_square_root(n)
if m * m == n
@trivial_factor = m
return
end
# select parameters
# @k: a multiplier
@k = MPQS.k_table[@n % 4][SMALL_PRIMES.map {|q| Utils.kronecker(@n, q) }]
@kn = @k * @n
# @f: the size of the factor base
# @m: the length of the sieve interval (2M+1)
# @t: the large prime tolerance
_digit, @f, @m, @t =
PARAMETERS.find {|d,| Math.log(@n, 10) < d } || PARAMETERS.last
# compute factor base
# @ps: the factor base
# @log_ps: log of each factor base
# @sqrts: sqrt(kN) mod p
@ps, @log_ps, @sqrts = [], [], []
Sieve.each do |p|
if @n % p == 0
@trivial_factor = p
return
end
if p >= 3 && Utils.kronecker(@kn, p) == 1
@ps << p
@log_ps << (Math.log(p) * S).floor
@sqrts << Utils.mod_sqrt(@kn, p)
break if @ps.size >= @f
end
end
@pmax = @ps.last # max factor base
# for enumerating coefficients
@base_d = Math.sqrt(Math.sqrt(@kn) / (Math.sqrt(2) * @m)).floor | 3
@diff_d = 0 # 0, 4, -4, 8, -8, 12, -12, ...
# sieve table
@ws = [0] * (2 * @m + 1)
# a test value for sieve
@test_value = (Math.log(@m * Math.sqrt(@kn) / @pmax ** @t) * S).round
# relations found
@incomplete_relations = {}
# for gaussian elimination
@elim = GaussianElimination.new(@ps.size + 2)
end
# try to find any factor (not necessarily a prime)
def try_find_factor
return @trivial_factor if defined?(@trivial_factor)
catch(:factor_found) do
while true
p :foo
select_q
solve_roots
sieve
@ws.each_with_index do |w, i|
check_relation(i - @m) if w > @test_value
end
end
end
end
# Selects Q(x) = Ax^2 + Bx + C for quadratic sieve
def select_q
# find a prime D such that (kN | D) = 1 and D = 3 mod 4
begin
@d = @base_d + @diff_d
@diff_d = @diff_d > 0 ? -@diff_d : 4 - @diff_d
end until Utils.kronecker(@kn, @d) == 1 && PrimalityTest.probable_prime?(@d) != false && !@ps.include?(@d)
# select coefficients (see the paper)
@a = @d * @d
h0 = Utils.mod_pow(@kn, @d / 4, @d)
h1 = (@kn * h0) % @d
# assert: h1*h1 = kN mod D
# assert: h0*h1 = 1 mod D
h2 = (Utils.mod_inv(2, @d) * h0 * ((@kn - h1 * h1) / @d)) % @d
# assert: kN - h1^2 = 0 mod D
@b = (h1 + h2 * @d) % @a
@b -= @a if @b[0] == 0
# assert: B^2 = kN mod 4A
#@c = (@b * @b - @kn) / (4 * @a) # not needed
# compute real roots of Q(x) = 0
s = Utils.integer_square_root(@kn)
@real_root_1 = (-@b - s) / (2 * @a)
@real_root_2 = (-@b + s) / (2 * @a)
# x is between real roots if @real_root_1 <= x && x <= @real_root_2
end
# Computes the roots of Q(x) = 0 mod p (p in the factor base)
def solve_roots
@roots = []
@ps.zip(@sqrts) do |p, s|
a_inv = Utils.mod_inv(2 * @a, p)
x1 = (-@b + s) * a_inv % p
x2 = (-@b - s) * a_inv % p
@roots << [x1, x2]
end
end
# Performs the quadratic sieve
def sieve
@ws.fill(0)
@ps.zip(@log_ps, @roots) do |p, lp, (s1, s2)|
((@m + s1) % p).step(2 * @m, p) {|i| @ws[i] += lp }
((@m + s2) % p).step(2 * @m, p) {|j| @ws[j] += lp }
end
end
# Tests whether Q(x) is actually a product of factor bases
def check_relation(x)
h, qx = compute_q(x)
return if qx == 0
es, l = exponent_bitvector(qx)
# discard this x if the residue L is too big
return if l > @pmax ** @t
if l == 1
# complete relation found:
# Q(x) = p0^e0 * p1^e1 * ... * pk^ek (pi in the factor base)
qx_vec = Hash.new(0)
PrimeFactorization.prime_factorization(qx) {|p, e| qx_vec[p] += e }
collect_relation(es, h, qx_vec)
elsif @incomplete_relations[l]
# large prime procedure:
# make a complete relation by multiplying two incomplete relations
es2, h2, qx2 = @incomplete_relations[l]
# XXX: use FactoredInteger
qx_vec = Hash.new(0)
PrimeFactorization.prime_factorization(qx) {|p, e| qx_vec[p] += e }
PrimeFactorization.prime_factorization(qx2) {|p, e| qx_vec[p] += e }
collect_relation(es ^ es2, h * h2 % @kn, qx_vec)
else
@incomplete_relations[l] = [es, h, qx]
end
end
# Adds a complete relation found to gaussian elimination list
def collect_relation(es, h, q)
#p "%0#{ @ps.size + 1 }b" % es
@elim[es] = [h, q]
if @elim.size > 0.9 * @f && @elim.size % [1, @f / 100].max == 0
@elim.eliminate do |relations|
# factor candidate found
check_factor_candidate(relations)
end
end
if @elim.size > @f * 1.5
raise Failed, "failed to find a factor: #{ @n }"
end
end
# Computes and checks a factor candidate
def check_factor_candidate(relations)
# computes the factor candidate
p1 = 1
p2_vec = Hash.new(0) # XXX: use FactoredInteger
relations.each do |h, qx_vec|
p1 = p1 * h % @kn
qx_vec.each {|p, e| p2_vec[p] += e }
end
p2 = 1
p2_vec.each {|p, e| p2 = (p2 * Utils.mod_pow(p, e / 2, @kn)) % @kn }
return if p1 == p2 || (p1 + p2) % @kn == 0
factor = (p1 + p2).gcd(@n)
# check the factor candidate
if factor != 1 && factor != @n
throw :factor_found, factor
end
end
# Computes Q(x)
def compute_q(x)
h = (2 * @a * x + @b) * Utils.mod_inv(2 * @d, @kn) % @kn
q = h * h % @kn
q -= @kn if @real_root_1 <= x && x <= @real_root_2
# assert: q == @a*x*x + @b*x + @c
return h, q
end
# Factorizes n in factor base and returns an exponent bitvector
def exponent_bitvector(n)
vec = 0
if n < 0
vec = 1
n = -n
end
bit = 2
([2] + @ps).each do |p|
e = false
while n % p == 0
n /= p
e = !e
end
vec |= bit if e
bit <<= 1
end
return vec, n
end
# Incremental gaussian elimination.
#
# A. K. Lenstra, and M. S. Manasse,
# "Compact incremental Gaussian Elimination over Z/2Z." (1988)
# http://dl.acm.org/citation.cfm?id=896266
class GaussianElimination
def initialize(m)
@n = 0
@m = m
@d = []
@e = []
@u = []
@v = []
end
def []=(es, val)
@d << es
@e << nil
@v << val
end
def size
@d.size
end
def eliminate
n2 = @d.size
@n.times do |i|
if @u[i]
@e[@u[i]] = i
@n.upto(n2 - 1) do |j|
@d[j] ^= @d[i] if @d[j][@u[i]] == 1
end
end
end
@n.upto(n2 - 1) do |j|
found = false
@m.times do |i|
if @d[j][i] == 1 && !@e[i]
@u[j] = i
@e[@u[j]] = j
@d[j] &= ~(1 << @u[j])
(j + 1).upto(n2 - 1) do |k|
@d[k] ^= @d[j] if @d[k][@u[j]] == 1
end
found = true
break
end
end
if !found
# linear dependency found
@u[j] = nil
c = {}
c[j] = true
@m.times do |k|
c[@e[k]] = true if @d[j][k] == 1
end
yield c.keys.map {|k| @v[k] }
end
end
@n = n2
end
end
end
end
# == faster_prime/sieve from faster_prime/base ---------------------------------
# Sieve of Atkin
# http://cr.yp.to/papers/primesieves.pdf
# require "faster_prime/utils" # (expanded: L22)
module FasterPrime
module Sieve
module_function
TBL = []
[
[:mark_type_1, ->(f, g) { 4*f*f + g*g }, 15, 30, 1, 4],
[:mark_type_2, ->(f, g) { 3*f*f + g*g }, 10, 30, 7, 12],
[:mark_type_3, ->(f, g) { 3*f*f - g*g }, 10, 30, 11, 12],
].each do |type, expr, fmax, gmax, n, mod|
gmax.times do |g|
fmax.times do |f|
d = expr[f, g] % 60
TBL << [type, f, g, d] if d % mod == n && d.gcd(60) == 1
end
end
end
MODULO60 = [1, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 49, 53, 59]
def each(ubound = nil, &blk)
if ubound && SMALL_PRIMES.last >= ubound
SMALL_PRIMES.each(&blk)
return
end
u = SMALL_PRIMES.last / 60 * 60
SMALL_PRIMES.pop until SMALL_PRIMES.last < u
SMALL_PRIMES.each(&blk)
b = Utils.integer_square_root(ubound * 3000) if ubound
b = 100000 if !b || b > 100000
l = u / 60
# Algorithm 3.1, 3.2, 3.3
flags = [0] * b
while true
# step
flags.fill(0)
# sieve
TBL.each do |name, f, g, d|
send(name, flags, f, g, d, l, b)
end
# squarefree
squarefree(flags, l, b)
# enumerate primes
if ubound && ubound < 60 * (l + b)
enumerate(flags, l, b) do |p|
return if p >= ubound
yield p
end
else
enumerate(flags, l, b, &blk)
end
l += b
end
end
def enumerate(flags, l, b)
b.times do |x|
base = 60 * (l + x)
flag = flags[x]
MODULO60.each do |d|
next if flag[d] == 0
p = base + d
SMALL_PRIMES << p
yield p
end
end
end
def squarefree(flags, l, b)
SMALL_PRIMES.each do |q|
q2 = q * q
break if q2 > 60 * (l + b)
next if q < 7
minus_inv_60 = -Utils.mod_inv(60, q2)
MODULO60.each do |d|
# 60(l+x)+d == 0 (mod q2)
# x = -d * 60^{-1} - l (mod q2)
x = (d * minus_inv_60 - l) % q2
m = ~(1 << d)
while x < b
flags[x] &= m
x += q2
end
end
end
end
# Algorithm 4.1
def mark_type_1(flags, f, g, d, l, b)
x, y0 = f, g
k0 = (4 * f * f + g * g - d) / 60 - l
m = 1 << d
while k0 < l + b
k0 += 2*x + 15
x += 15
end
while true
x -= 15
k0 -= 2*x + 15
return if x <= 0
while k0 < 0
k0 += y0 + 15
y0 += 30
end
k, y = k0, y0
while k < b
flags[k] ^= m # mark
k += y + 15
y += 30
end
end
end
# Algorithm 4.2
def mark_type_2(flags, f, g, d, l, b)
x, y0 = f, g
k0 = (3 * f * f + g * g - d) / 60 - l
m = 1 << d
while k0 < b
k0 += x + 5
x += 10
end
while true
x -= 10
k0 -= x + 5
return if x <= 0
while k0 < 0
k0 += y0 + 15
y0 += 30
end
k, y = k0, y0
while k < b
flags[k] ^= m # mark
k += y + 15
y += 30
end
end
end
# Algorithm 4.3
def mark_type_3(flags, f, g, d, l, b)
x, y0 = f, g
k0 = (3 * f * f - g * g - d) / 60 - l
m = 1 << d
while true
while k0 >= b
return if x <= y0
k0 -= y0 + 15
y0 += 30
end
k, y = k0, y0
while k >= 0 && y < x
flags[k] ^= m # mark
k -= y + 15
y += 30
end
k0 += x + 5
x += 10
end
end
end
end
# == faster_prime/core_ext from faster_prime/base ------------------------------
class Integer
def Integer.from_prime_division(pd)
FasterPrime.int_from_prime_division(pd)
end
def prime_division(generator = nil)
FasterPrime.prime_division(self)
end
def prime?
FasterPrime.prime?(self)
end
def Integer.each_prime(ubound, &block) # :yields: prime
Faster.Prime.each(ubound, &block)
end
end
# == faster_prime/base from faster_prime ---------------------------------------
# require "faster_prime/version" # (expanded: L705)
# require "faster_prime/utils" # (expanded: L22)
# require "faster_prime/primality_test" # (expanded: L939)
# require "faster_prime/prime_factorization" # (expanded: L1499)
# require "faster_prime/sieve" # (expanded: L1960)
# require "faster_prime/core_ext" # (expanded: L2137)
module FasterPrime
module Core
include Enumerable
class SuccEnumerator < Enumerator
alias succ next
end
def each(ubound = nil, &block)
return SuccEnumerator.new(Float::INFINITY) do |y|
Sieve.each(ubound) {|p| y << p }
end unless block_given?
Sieve.each(ubound, &block)
end
def prime?(n)
raise ArgumentError, "Expected an integer, got #{ n }" unless n.respond_to?(:integer?) && n.integer?
PrimalityTest.prime?(n)
end
def prime_division(n)
PrimeFactorization.prime_factorization(n).to_a
end
def int_from_prime_division(pd)
pd.inject(1) {|v, (p, e)| v * (p ** e) }
end
end
def self.instance
self
end
include Core
extend Core
end
# == faster_prime/core_ext from faster_prime -----------------------------------
class Integer
def Integer.from_prime_division(pd)
FasterPrime.int_from_prime_division(pd)
end
def prime_division(generator = nil)
FasterPrime.prime_division(self)
end
def prime?
FasterPrime.prime?(self)
end
def Integer.each_prime(ubound, &block) # :yields: prime
Faster.Prime.each(ubound, &block)
end
end
# == faster_prime from main ----------------------------------------------------
# require "faster_prime/base" # (expanded: L2156)
# require "faster_prime/core_ext" # (expanded: L2137)
# == nanacl/bsearch_right from main --------------------------------------------
# frozen_string_literal: true
class Array
def bsearch_right(&)
index = bsearch_index_right(&)
index && self[index]
end
def bsearch_index_right(&block)
right = bsearch_index { |elem| !block.call(elem) }
if right.nil?
size - 1
elsif right == 0
nil
else
right - 1
end
end
end
class Range
def bsearch_right(&block)
right = bsearch { |elem| !block.call(elem) }
if right.nil?
last
elsif right == first
nil
else
right - 1
end
end
end
# ==============================================================================
main.call