12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196 |
- #!/usr/bin/ruby
- # This script factorizes a natural number given as a command line
- # parameter into its prime factors. It first attempts to use trial
- # division to find very small factors, then uses Brent's version of the
- # Pollard rho algorithm [1] to find slightly larger factors. If any large
- # factors remain, it uses the Self-Initializing Quadratic Sieve (SIQS) [2]
- # to factorize those.
- # [1] Brent, Richard P. 'An improved Monte Carlo factorization algorithm.'
- # BIT Numerical Mathematics 20.2 (1980): 176-184.
- # [2] Contini, Scott Patrick. 'Factoring integers with the self-
- # initializing quadratic sieve.' (1997).
- STDOUT.autoflush(true)
- class SIQSPolynomial (coeff, a=nil, b=nil) {
- method eval (x) {
- var res = 0
- coeff.each {|k|
- res *= x
- res += k
- }
- return res
- }
- }
- class FactorBasePrime (p, t, lp) {
- has soln1 = nil
- has soln2 = nil
- has ainv = nil
- has p_log = p.log
- }
- func fibonacci_factorization (n, upper_bound) {
- var bound = (5 * n.ilog2**2)
- bound = upper_bound if (bound > upper_bound)
- loop {
- return nil if (bound <= 1)
- var B = bound.consecutive_lcm
- var F = B.fibmod(n)
- var g = gcd(F, n)
- if (g == n) {
- bound >>= 1
- next
- }
- if (g > 1) {
- return [g]
- }
- return nil
- }
- }
- func fast_fibonacci_factor (n, upto) {
- var (P, Q) = (3, 1)
- for k in (2 .. upto) {
- var (U, V) = lucasUVmod(P, Q, k, n)
- for t in (U, V, V-1, V+1, V-P) {
- var g = gcd(t, n)
- if (g.is_ntf(n)) {
- return [g]
- }
- }
- }
- return nil
- }
- # Some tuning parameters
- class SIQS (
- MASK_LIMIT = 200, # show Cn if n > MASK_LIMIT, where n ~ log_10(N)
- LOOK_FOR_SMALL_FACTORS = 1,
- TRIAL_DIVISION_LIMIT = 1_000_000,
- PHI_FINDER_ITERATIONS = 100_000,
- FERMAT_ITERATIONS = 100_000,
- NEAR_POWER_ITERATIONS = 1_000,
- PELL_ITERATIONS = 100_000,
- FLT_ITERATIONS = 200_000,
- HOLF_ITERATIONS = 100_000,
- MBE_ITERATIONS = 100,
- DOP_FACTOR_LIMIT = 1000,
- COP_FACTOR_LIMIT = 500,
- MILLER_RABIN_ITERATIONS = 100,
- LUCAS_MILLER_ITERATIONS = 50,
- SIQS_TRIAL_DIVISION_EPS = 25,
- SIQS_MIN_PRIME_POLYNOMIAL = 400,
- SIQS_MAX_PRIME_POLYNOMIAL = 4000,
- ) {
- method factor_base_primes (n, nf) {
- var factor_base = []
- 1e7.each_prime {|p|
- var t = n.sqrtmod(p) || next
- var lp = p.log2.round.int
- factor_base << FactorBasePrime(p, t, lp)
- break if (factor_base.len >= nf)
- }
- return factor_base
- }
- method create_poly (a, b, n, factor_base, first) {
- var b_orig = b
- if (2*b > a) {
- b = (a - b)
- }
- var g = SIQSPolynomial([a*a, 2*a*b, b*b - n], a, b_orig)
- var h = SIQSPolynomial([a, b])
- for fb in (factor_base) {
- next if (fb.p `divides` a)
- fb.ainv = a.invmod(fb.p) if first
- fb.soln1 = ((fb.ainv * ( fb.t - b)) % fb.p)
- fb.soln2 = ((fb.ainv * (-fb.t - b)) % fb.p)
- }
- return (g, h)
- }
- method find_first_poly (n, m, factor_base) {
- var p_min_i = nil
- var p_max_i = nil
- factor_base.each_kv { |i,fb|
- if (!defined(p_min_i) && (fb.p >= SIQS_MIN_PRIME_POLYNOMIAL)) {
- p_min_i = i
- }
- if (!defined(p_max_i) && (fb.p > SIQS_MAX_PRIME_POLYNOMIAL)) {
- p_max_i = i-1
- break
- }
- }
- # The following may happen if the factor base is small
- if (!defined(p_max_i)) {
- p_max_i = factor_base.end
- }
- if (!defined(p_min_i)) {
- p_min_i = 5
- }
- if ((p_max_i - p_min_i) < 20) {
- p_min_i = (p_min_i `min` 5)
- }
- var target0 = ((n.log + 2.log)/2 - m.log)
- var target1 = (target0 - ((log(factor_base[p_min_i].p + factor_base[p_max_i].p) - 2.log) / 2))
- # find q such that the product of factor_base[q_i] is approximately
- # sqrt(2 * n) / m; try a few different sets to find a good one
- var (best_q, best_a, best_ratio)
- 30.times {
- var A = 1
- var log_A = 0
- var Q = Hash()
- while (log_A < target1) {
- var p_i = 0
- while ((p_i == 0) || Q.exists(p_i)) {
- p_i = p_min_i.irand(p_max_i) # inclusive on both sides
- }
- var fb = factor_base[p_i]
- A *= fb.p
- log_A += fb.p_log
- Q{p_i} = fb
- }
- var ratio = exp(log_A - target0)
- # ratio too small seems to be not good
- if (!defined(best_ratio) ||
- (( ratio >= 0.9) && (ratio < best_ratio)) ||
- ((best_ratio < 0.9) && (ratio > best_ratio))
- ) {
- best_q = Q
- best_a = A
- best_ratio = ratio
- }
- }
- var a = best_a
- var b = 0
- var arr = []
- best_q.each_v { |fb|
- var p = fb.p
- var r = a/p
- var gamma = ((fb.t * r.invmod(p)) % p)
- if (gamma > (p >> 1)) {
- gamma = (p - gamma)
- }
- var t = r*gamma
- b += t
- arr << t
- }
- var (g, h) = self.create_poly(a, b, n, factor_base, true)
- return (g, h, arr)
- }
- method find_next_poly (n, factor_base, i, g, arr) {
- # Compute the (i+1)-th polynomials for the Self-Initializing
- # Quadratic Sieve, given that g is the i-th polynomial.
- var v = i.valuation(2)
- var z = ((((i >> (v+1)) & 1) == 0) ? -1 : 1)
- var a = g.a
- var b = ((g.b + 2*z*arr[v]) % a)
- return self.create_poly(a, b, n, factor_base, false)
- }
- method siqs_sieve (factor_base, m) {
- # Perform the sieving step of the SIQS. Return the sieve array.
- var sieve_array = (2*m + 1 -> of(0))
- factor_base.each {|fb|
- fb.p > 100 || next
- fb.soln1 \\ next
- var p = fb.p
- var lp = fb.lp
- var end = 2*m
- var i_start_1 = -((m + fb.soln1) // p) #/
- var a_start_1 = (fb.soln1 + i_start_1*p)
- for i in (RangeNum(a_start_1 + m, end, p)) {
- sieve_array[i] += lp
- }
- var i_start_2 = -((m + fb.soln2) // p) #/
- var a_start_2 = (fb.soln2 + i_start_2*p)
- for i in (RangeNum(a_start_2 + m, end, p)) {
- sieve_array[i] += lp
- }
- }
- return sieve_array
- }
- method siqs_trial_divide (n, factor_base_info) {
- # Determine whether the given number a can be fully factorized into
- # primes from the factors base. If so, return the indices of the
- # factors from the factor base. If not, return `nil`.
- if (n.is_smooth_over_prod(factor_base_info{:prod})) {
- var factor_index = factor_base_info{:index}
- return n.factor_exp.map {|p|
- [factor_index{p[0]}, p[1]]
- }
- }
- return nil
- }
- method siqs_trial_division (n, sieve_array, factor_base_info, smooth_relations, g, h, m, req_relations) {
- # Perform the trial division step of the Self-Initializing Quadratic Sieve.
- var limit = ((m.log + n.log/2)/2.log - SIQS_TRIAL_DIVISION_EPS)
- sieve_array.each_kv { |i,s|
- next if (s < limit)
- var x = (i - m)
- var gx = g.eval(x).abs
- var divisors_idx = self.siqs_trial_divide(gx, factor_base_info) \\ next
- var u = h.eval(x)
- var v = gx
- smooth_relations << [u, v, divisors_idx]
- if (smooth_relations.len >= req_relations) {
- return true
- }
- }
- return false
- }
- method build_matrix (factor_base, smooth_relations) {
- # Build the matrix for the linear algebra step of the Quadratic Sieve.
- var fb = factor_base.len
- var matrix = []
- smooth_relations.each {|sr|
- var row = fb.of(0)
- for u,v in (sr[2]) {
- row[u] = v&1
- }
- matrix << row
- }
- return matrix
- }
- method build_matrix_opt(M) {
- # Convert the given matrix M of 0s and 1s into a list of numbers m
- # that correspond to the columns of the matrix.
- # The j-th number encodes the j-th column of matrix M in binary:
- # The i-th bit of m[i] is equal to M[i][j].
- var m = M[0].len
- var cols_binary = m.of('')
- M.each {|mi|
- mi.each_kv { |j, mij|
- cols_binary[j] += mij
- }
- }
- return (cols_binary.map {|s| Number(s.flip, 2) }, M.len, m)
- }
- method find_pivot_column_opt (M, j) {
- # For a matrix produced by `build_matrix_opt`, return the row of
- # the first non-zero entry in column j, or None if no such row exists.
- var v = M[j]
- if (v == 0) {
- return nil
- }
- return v.valuation(2)
- }
- method solve_matrix_opt (M, n, m) {
- # Perform the linear algebra step of the SIQS. Perform fast
- # Gaussian elimination to determine pairs of perfect squares mod n.
- # Use the optimizations described in [1].
- # [1] Koç, Çetin K., and Sarath N. Arachchige. 'A Fast Algorithm for
- # Gaussian Elimination over GF (2) and its Implementation on the
- # GAPP.' Journal of Parallel and Distributed Computing 13.1
- # (1991): 118-122.
- var row_is_marked = n.of(false)
- var pivots = m.of(-1)
- for j in (m.range) {
- var i = self.find_pivot_column_opt(M, j) \\ next
- pivots[j] = i
- row_is_marked[i] = true
- for k in (m.range) {
- if ((k != j) && M[k].bit(i)) {
- M[k] ^= M[j]
- }
- }
- }
- var perf_squares = []
- for i in (n.range) {
- next if row_is_marked[i]
- var perfect_sq_indices = [i]
- for j in (m.range) {
- if (M[j].bit(i)) {
- perfect_sq_indices << pivots[j]
- }
- }
- perf_squares << perfect_sq_indices
- }
- return perf_squares
- }
- method calc_sqrts (n, square_indices, smooth_relations) {
- # Given on of the solutions returned by `solve_matrix_opt` and
- # the corresponding smooth relations, calculate the pair [a, b], such
- # that a^2 = b^2 (mod n).
- var r1 = 1
- var r2 = 1
- square_indices.each { |i|
- r1 *= smooth_relations[i][0] %= n
- r2 *= smooth_relations[i][1]
- }
- return (r1, r2.isqrt)
- }
- method factor_from_square (n, square_indices, smooth_relations) {
- # Given one of the solutions returned by `solve_matrix_opt`,
- # return the factor f determined by f = gcd(a - b, n), where
- # a, b are calculated from the solution such that a*a = b*b (mod n).
- # Return f, a factor of n (possibly a trivial one).
- var (sqrt1, sqrt2) = self.calc_sqrts(n, square_indices, smooth_relations)
- return gcd(sqrt1 - sqrt2, n)
- }
- method siqs_find_more_factors_gcd(numbers) {
- var res = Hash()
- numbers.each_kv { |i,n|
- res{n} = n
- for k in (i+1 .. numbers.end) {
- var m = numbers[k]
- var f = gcd(n, m)
- if ((f != 1) && (f != n) && (f != m)) {
- if (!res.exists(f)) {
- say "SIQS: GCD found non-trivial factor: #{f}"
- res{f} = f
- }
- var t1 = n//f #/
- var t2 = m//f #/
- res{t1} = t1
- res{t2} = t2
- }
- }
- }
- return res.values
- }
- method siqs_find_factors (n, perfect_squares, smooth_relations) {
- # Perform the last step of the Self-Initializing Quadratic Field.
- # Given the solutions returned by `solve_matrix_opt`, attempt to
- # identify a number of (not necessarily prime) factors of n, and
- # return them.
- var factors = []
- var rem = n
- var non_prime_factors = Hash()
- var prime_factors = Hash()
- perfect_squares.each {|square_indices|
- var fact = self.factor_from_square(n, square_indices, smooth_relations)
- next if (fact <= 1)
- next if (fact > rem)
- if (fact.is_prime) {
- if (!prime_factors.exists(fact)) {
- say "SIQS: Prime factor found: #{fact}"
- prime_factors{fact} = fact
- }
- rem = self.check_factor(rem, fact, factors)
- break if rem.is_one
- if (rem.is_prime) {
- factors << rem
- rem = 1
- break
- }
- if (defined(var root = self.check_perfect_power(rem))) {
- say "SIQS: Perfect power detected with root: #{root}"
- factors << root
- rem = 1
- break
- }
- }
- else {
- if (!non_prime_factors.exists(fact)) {
- say "SIQS: Composite factor found: #{fact}"
- non_prime_factors{fact} = fact
- }
- }
- }
- if ((rem != 1) && non_prime_factors) {
- non_prime_factors{rem} = rem
- var primes = []
- var composites = []
- self.siqs_find_more_factors_gcd(non_prime_factors.values).each {|fact|
- if (fact.is_prime) {
- primes << fact
- }
- elsif (fact > 1) {
- composites << fact
- }
- }
- for fact in (primes + composites) {
- if ((fact != rem) && (fact `divides` rem)) {
- say "SIQS: Using non-trivial factor from GCD: #{fact}"
- rem = self.check_factor(rem, fact, factors)
- }
- break if rem.is_one
- break if rem.is_prime
- }
- }
- if (rem != 1) {
- factors << rem
- }
- return factors
- }
- method siqs_choose_range (n) {
- # Choose m for sieving in [-m, m].
- exp(sqrt(n.log * n.log.log) / 2) -> round.int
- }
- method siqs_choose_nf (n) {
- # Choose parameters nf (sieve of factor base)
- exp(sqrt(n.log * n.log.log))**(2.sqrt / 4) -> round.int
- }
- method siqs_factorize (n, nf) {
- # Use the Self-Initializing Quadratic Sieve algorithm to identify
- # one or more non-trivial factors of the given number n. Return the
- # factors as a list.
- var m = self.siqs_choose_range(n)
- var factors = []
- var factor_base = self.factor_base_primes(n, nf)
- var factor_prod = factor_base.prod_by { .p }
- var factor_index = Hash(
- factor_base.map { .p } ~Z ^factor_base -> flat...
- )
- var factor_base_info = Hash(
- base => factor_base,
- prod => factor_prod,
- index => factor_index,
- )
- var smooth_relations = []
- var required_relations_ratio = 1
- var success = false
- var prev_cnt = 0
- var i_poly = 0
- var (g, h, arr)
- while (!success) {
- say "*** Step 1/2: Finding smooth relations ***"
- say "SIQS sieving range: [-#{m}, #{m}]"
- var required_relations = ((factor_base.len + 1) * required_relations_ratio).round.int
- say "Target: #{required_relations} relations."
- var enough_relations = false
- while (!enough_relations) {
- if (i_poly == 0) {
- (g, h, arr) = self.find_first_poly(n, m, factor_base);
- }
- else {
- (g, h) = self.find_next_poly(n, factor_base, i_poly, g, arr);
- }
- if (++i_poly >= (1 << arr.end)) {
- i_poly = 0
- }
- var sieve_array = self.siqs_sieve(factor_base, m)
- enough_relations = self.siqs_trial_division(
- n, sieve_array, factor_base_info, smooth_relations,
- g, h, m, required_relations
- )
- if ((smooth_relations.len >= required_relations) ||
- (smooth_relations.len > prev_cnt)
- ) {
- printf("Progress: %d/%d relations.\r", smooth_relations.len, required_relations)
- prev_cnt = smooth_relations.len
- }
- }
- say "\n\n*** Step 2/2: Linear Algebra ***"
- say "Building matrix for linear algebra step..."
- var M = self.build_matrix(factor_base, smooth_relations)
- var (M_opt, M_n, M_m) = self.build_matrix_opt(M)
- say "Finding perfect squares using Gaussian elimination...";
- var perfect_squares = self.solve_matrix_opt(M_opt, M_n, M_m)
- say "Finding factors from perfect squares...\n";
- factors = self.siqs_find_factors(n, perfect_squares, smooth_relations)
- if (factors.len >= 2) {
- success = true
- }
- else {
- say "Failed to find a solution. Finding more relations..."
- required_relations_ratio += 0.05
- }
- }
- return factors
- }
- method trial_division_small_primes (n) {
- # Perform trial division on the given number n using all primes up
- # to upper_bound. Initialize the global variable small_primes with a
- # list of all primes <= upper_bound. Return (factors, rem), where
- # factors is the list of identified prime factors of n, and rem is the
- # remaining factor. If rem = 1, the function terminates early, without
- # fully initializing small_primes.
- say "[*] Trial division..."
- var factors = n.trial_factor(TRIAL_DIVISION_LIMIT)
- if (factors.last.is_prime) {
- return (factors, 1)
- }
- var rem = factors.pop
- return (factors, rem)
- }
- method check_perfect_power (n) {
- # Check whether n is a perfect and return its perfect root.
- # Returns undef otherwise.
- if (n.is_power) {
- var power = n.perfect_power
- var root = n.perfect_root
- say "`-> #{n} is #{root}^#{power}"
- return root
- }
- return nil
- }
- method check_factor (Number n, Number i, Array factors) {
- while (i `divides` n) {
- n //= i #/
- factors << i
- if (n.is_prime) {
- factors << n
- return 1
- }
- }
- return n
- }
- method store_factor (Ref rem, Number f, Array factors) {
- f || return false
- f < *rem || return false
- f `divides` *rem || die "error: #{f} does not divide #{*rem}"
- if (f.is_prime) {
- say "`-> prime factor: #{f}"
- *rem = self.check_factor(*rem, f, factors)
- }
- else {
- say "`-> composite factor: #{f}"
- *rem = idiv(*rem, f)
- # Try to find a small factor of f
- var f_factor = self.find_small_factors(f, factors)
- if (f_factor < f) {
- *rem *= f_factor
- }
- else {
- # Use SIQS to factorize f
- var F = self.find_prime_factors(f)
- F.each {|p|
- if (p `divides` *rem) {
- *rem = self.check_factor(*rem, p, factors)
- break if (*rem == 1)
- }
- }
- }
- }
- return true
- }
- method find_small_factors (rem, factors) {
- # Perform up to max_iter iterations of Brent's variant of the
- # Pollard rho factorization algorithm to attempt to find small
- # prime factors. Restart the algorithm each time a factor was found.
- # Add all identified prime factors to factors, and return 1 if all
- # prime factors were found, or otherwise the remaining factor.
- var state = Hash(
- cyclotomic_check => true,
- fast_power_check => true,
- fast_fibonacci_check => true,
- )
- var len = rem.len
- var factorization_methods = [
- func {
- say "=> Miller-Rabin method..."
- miller_factor(rem, (len > 1000) ? 15 : MILLER_RABIN_ITERATIONS)
- },
- func {
- if (len < 3000) {
- say "=> Lucas-Miller method (n+1)..."
- lucas_factor(rem, +1, (len > 1000) ? 10 : LUCAS_MILLER_ITERATIONS)
- }
- },
- func {
- if (len < 3000) {
- say "=> Lucas-Miller method (n-1)..."
- lucas_factor(rem, -1, (len > 1000) ? 10 : LUCAS_MILLER_ITERATIONS)
- }
- },
- func {
- say "=> Phi finder method...";
- phi_finder_factor(rem, PHI_FINDER_ITERATIONS)
- },
- func {
- say "=> Fermat's method..."
- fermat_factor(rem, FERMAT_ITERATIONS)
- },
- func {
- say "=> HOLF method..."
- holf_factor(rem, 2*HOLF_ITERATIONS)
- },
- func {
- say "=> Pell factor..."
- pell_factor(rem, PELL_ITERATIONS)
- },
- func {
- say "=> Pollard p-1 (20K)..."
- pm1_factor(rem, 20_000)
- },
- func {
- say "=> Fermat's little theorem (base 2)..."
- flt_factor(rem, 2, (len > 1000) ? 1e4 : FLT_ITERATIONS)
- },
- func {
- var len_2 = (len * (log(10) / log(2)))
- var iter = (((len_2 * MBE_ITERATIONS) > 1_000) ? int(1_000/len_2) : MBE_ITERATIONS)
- if (iter > 0) {
- say "=> MBE method (#{iter} iter)...";
- mbe_factor(rem, iter)
- }
- },
- func {
- say "=> Fermat's little theorem (base 3)..."
- flt_factor(rem, 3, (len > 1000) ? 1e4 : FLT_ITERATIONS)
- },
- func {
- state{:fast_fibonacci_check} || return nil
- say "=> Fast Fibonacci check..."
- var f = fast_fibonacci_factor(rem, 5000)
- f || do { state{:fast_fibonacci_check} = false }
- f
- },
- func {
- state{:cyclotomic_check} || return nil
- say "=> Fast cyclotomic check..."
- var f = cyclotomic_factor(rem)
- f.len >= 2 || do { state{:cyclotomic_check} = false }
- f
- },
- func {
- say "=> Pollard rho (10M)..."
- prho_factor(rem, 1e10.isqrt)
- },
- func {
- say "=> Pollard p-1 (500K)..."
- pm1_factor(rem, 500_000)
- },
- func {
- say "=> ECM (600)..."
- ecm_factor(rem, 600, 20)
- },
- func {
- say "=> Williams p±1 (500K)..."
- pp1_factor(rem, 500_000)
- },
- func {
- if (len < 1000) {
- say "=> Chebyshev p±1 (500K)..."
- chebyshev_factor(rem, 500_000, 1e6.irand+2)
- }
- },
- func {
- say "=> Williams p±1 (1M)..."
- pp1_factor(rem, 1_000_000)
- },
- func {
- if (len < 1000) {
- say "=> Chebyshev p±1 (1M)..."
- chebyshev_factor(rem, 1_000_000, 1e6.irand+2)
- }
- },
- func {
- say "=> ECM (2K)..."
- ecm_factor(rem, 2000, 10)
- },
- func {
- say "=> Difference of powers..."
- dop_factor(rem, DOP_FACTOR_LIMIT)
- },
- func {
- if (len < 500) {
- say "=> Fibonacci p±1 (500K)..."
- fibonacci_factorization(rem, 500_000)
- }
- },
- func {
- say "=> Pollard-Brent rho (12M)..."
- pbrent_factor(rem, 1e12.isqrt)
- },
- func {
- say "=> Congruence of powers..."
- cop_factor(rem, COP_FACTOR_LIMIT)
- },
- func {
- say "=> Pollard p-1 (5M)..."
- pm1_factor(rem, 5_000_000)
- },
- func {
- say "=> Williams p±1 (3M)..."
- pp1_factor(rem, 3_000_000)
- },
- func {
- say "=> Pollard rho (13M)..."
- prho_factor(rem, 1e13.isqrt)
- },
- func {
- say "=> ECM (160K)..."
- ecm_factor(rem, 160_000, 80)
- },
- ]
- loop {
- break if (rem <= 1)
- if (rem.is_prime) {
- factors << rem
- rem = 1
- break
- }
- len = rem.len
- if ((len >= 25) && (len <= 30)) { # factorize with SIQS directly
- return rem
- }
- printf("\n[*] Factoring %s (%s digits)...\n\n", (len > MASK_LIMIT ? "C#{len}" : rem), len)
- say "=> Perfect power check..."
- var f = self.check_perfect_power(rem)
- if (defined(f)) {
- var pow = rem.perfect_power
- var r = self.factorize(f)
- factors << (r*pow)...
- return 1
- }
- var found_factor = false
- var end = factorization_methods.end
- for (var i = 0; i <= end; ++i) {
- var f = factorization_methods[i].call || next
- f.each {|p|
- if (self.store_factor(\rem, p, factors)) {
- found_factor = true
- }
- }
- if (found_factor) {
- factorization_methods.unshift(factorization_methods.pop_at(i))
- break
- }
- else {
- factorization_methods.push(factorization_methods.pop_at(i))
- --i
- --end
- }
- }
- found_factor && next
- break
- }
- return rem
- }
- method find_prime_factors(n) {
- # Return one or more prime factors of the given number n. Assume
- # that n is not a prime and does not have very small factors, and that
- # the global small_primes has already been initialized. Do not return
- # duplicate factors.
- var factors = Hash()
- if (defined(var root = self.check_perfect_power(n))) {
- factors{root} = root
- }
- else {
- var digits = n.len
- say "\n[*] Using SIQS to factorize #{n} (#{digits} digits)...\n"
- var nf = self.siqs_choose_nf(n)
- var sf = self.siqs_factorize(n, nf)
- factors{sf...} = sf...
- }
- var prime_factors = []
- factors.each_v { |f|
- prime_factors += self.find_all_prime_factors(f)
- }
- return prime_factors
- }
- method verify_prime_factors (Number n, Array factors) {
- factors.prod == n || die "product of factors != n"
- factors.each {|p|
- p.is_prime || die "not prime detected: #{p}"
- }
- factors.sort
- }
- method find_all_prime_factors (n) {
- # Return all prime factors of the given number n. Assume that n
- # does not have very small factors and that the global small_primes
- # has already been initialized.
- var rem = n
- var factors = []
- while (rem > 1) {
- if (rem.is_prime) {
- factors << rem
- break
- }
- self.find_prime_factors(rem).each {|f|
- #~ rem != f || die 'error'
- #~ rem % f == 0 || die 'error'
- #~ is_prime(f) || die 'error'
- rem = self.check_factor(rem, f, factors);
- break if rem.is_one
- }
- }
- return factors
- }
- method factorize (n) {
- # Factorize the given integer n >= 1 into its prime factors.
- var orig = n
- if (n < 1) {
- die "Number needs to be an integer >= 1"
- }
- var len = n.len
- printf("\n[*] Factoring %s (%d digits)...\n", (len > MASK_LIMIT ? "C#{len}" : n), len)
- return [] if n.is_one
- return [n] if n.is_prime
- with (self.check_perfect_power(n)) {|root|
- var pow = n.perfect_power
- var factors = self.factorize(root)
- return self.verify_prime_factors(orig, factors*pow)
- }
- var divisors = (len>=30 ? Math.gcd_factors(
- n, n.dop_factor(DOP_FACTOR_LIMIT) + n.cop_factor(COP_FACTOR_LIMIT)
- ).first(-1) : [])
- if (divisors) {
- say "[*] Divisors found so far: #{divisors.join(', ')}";
- var composites = []
- var factors = []
- divisors.each {|d|
- d > 1 || next
- if (d.is_prime) {
- factors << d
- }
- else {
- composites << d
- }
- }
- factors << composites.reverse.map { self.factorize(_) }.flat...
- var rem = (orig / factors.prod)
- if (rem.is_prime) {
- factors << rem
- }
- elsif (rem > 1) {
- factors << self.factorize(rem)...
- }
- return self.verify_prime_factors(orig, factors)
- }
- var (factors, rem) = self.trial_division_small_primes(n)
- if (factors) {
- say "[*] Prime factors found so far: #{factors.join(', ')}"
- }
- else {
- say "[*] No small factors found..."
- }
- if (rem != 1) {
- if (LOOK_FOR_SMALL_FACTORS) {
- say "[*] Trying to find more small factors..."
- rem = self.find_small_factors(rem, factors)
- }
- else {
- say "[*] Skipping the search for more small factors..."
- }
- if (rem > 1) {
- factors += self.find_all_prime_factors(rem)
- }
- }
- factors.sort!
- assert_eq(factors.prod, n)
- factors.each {|p|
- assert(p.is_prime, "non-prime detected: #{p}")
- }
- return self.verify_prime_factors(orig, factors)
- }
- }
- if (ARGV) {
- var n = Number(ARGV[0])
- if (n.is_nan) {
- # evaluate the expression using PARI/GP
- n = Num(`echo #{ARGV[0].escape} | gp -q -f`)
- }
- if (n.is_nan) {
- die "Invalid expression: #{ARGV[0].dump}"
- }
- var siqs = SIQS()
- say "\nPrime factors: #{siqs.factorize(n)}"
- }
- else {
- say "usage: #{__MAIN__} <N>"
- }
|