123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253 |
- #!/usr/bin/ruby
- # Daniel "Trizen" Șuteu
- # Date: 25 October 2018
- # https://github.com/trizen
- # Simple implementation of the continued fraction factorization method (CFRAC),
- # combined with modular arithmetic (variation of the Brillhart-Morrison algorithm).
- # See also:
- # https://en.wikipedia.org/wiki/Pell%27s_equation
- # https://en.wikipedia.org/wiki/Continued_fraction_factorization
- # https://trizenx.blogspot.com/2018/10/continued-fraction-factorization-method.html
- # "Gaussian elimination" algorithm from:
- # https://github.com/martani/Quadratic-Sieve
- define(
- B_SMOOTH_METHOD = false, # true to use the B-smooth formula for the factor base
- ROUND_DIVISION = true, # true to use round division instead of floor division
- )
- func gaussian_elimination(A, n) {
- var m = A.end
- var I = (m+1).of {|i| 1 << i }
- var nrow = -1
- for col in (0 .. min(m, n)) {
- var npivot = -1
- for row in (nrow+1 .. m) {
- if (A[row].bit(col)) {
- npivot = row
- nrow++
- break
- }
- }
- next if (npivot < 0)
- if (npivot != nrow) {
- A.swap(npivot, nrow)
- I.swap(npivot, nrow)
- }
- for row in (nrow+1 .. m) {
- if (A[row].bit(col)) {
- A[row] ^= A[nrow]
- I[row] ^= I[nrow]
- }
- }
- }
- return I
- }
- func check_factor (n, g, factors) {
- while (g `divides` n) {
- n /= g;
- factors << g
- if (is_prime(n)) {
- factors << n
- return 1
- }
- }
- return n
- }
- func next_multiplier (k) {
- k += 2
- while (!k.is_squarefree) {
- ++k
- }
- return k
- }
- func cffmm (n, multiplier=1) {
- # Check for primes and negative numbers
- return [] if (n <= 1)
- return [n] if n.is_prime
- # Check for perfect powers
- if (n.is_power) {
- var root = n.perfect_root
- var power = n.perfect_power
- var factors = __FUNC__(root)
- return (factors * power -> sort)
- }
- var resolve_factor = {|a,b|
- var g = gcd(a - b.isqrt, n)
- if ((g > 1) && (g < n)) {
- return (
- __FUNC__(g) +
- __FUNC__(n/g) -> sort
- )
- }
- }
- var N = n*multiplier
- var x = N.isqrt
- var y = x
- var z = 1
- var w = x+x
- var r = w
- var B = int(exp(sqrt(log(n) * log(log(n))) / 2)) # B-smooth limit
- var nf = int(exp(sqrt(log(n) * log(log(n))))**(sqrt(2) / 4) / 1.25)
- var factor_base = if (B_SMOOTH_METHOD) {
- B.primes.grep {|p| (p <= 2) || (kronecker(N, p) >= 0) }
- }
- else {
- 1..Inf -> lazy.map { .prime }.grep {|p| (p <= 2) || (kronecker(N, p) >= 0) }.first(nf)
- }
- var factor_prod = factor_base.prod
- var factor_index = Hash(factor_base ~Z ^factor_base -> flat...)
- var L = factor_base.len+1
- var Q = []
- var A = []
- func exponent_signature(factors) {
- var sig = 0
- for p,e in factors {
- sig.setbit!(factor_index{p}) if e.is_odd
- }
- return sig
- }
- var (f1, f2) = (1, x)
- while (A.len < L) {
- y = (r*z - y)
- z = idiv(N - y*y, z) # exact division
- r = (ROUND_DIVISION ? idiv_round(x+y, z) : idiv(x+y, z))
- (f1, f2) = (f2, addmod(mulmod(r, f2, n), f1, n))
- if (z.is_square) {
- resolve_factor(f1, z)
- }
- if (z.abs.is_smooth_over_prod(factor_prod)) {
- var c_factors = z.abs.factor_exp
- if (c_factors) {
- A << exponent_signature(c_factors)
- Q << [f1, z.abs]
- }
- }
- if (z.abs == 1) {
- return __FUNC__(n, next_multiplier(multiplier))
- }
- }
- if (A.len < L) {
- A += (L - A.end + 1 -> of(0))
- }
- var rem = n
- var factors = []
- func cfrac_find_factors(solution) {
- var X = 1
- var Y = 1
- for i in (^Q) {
- if (solution.bit(i)) {
- X *= Q[i][0] %= n #=
- Y *= Q[i][1]
- var g = gcd(X - Y.isqrt, rem)
- if ((g > 1) && (g < rem)) {
- rem = check_factor(rem, g, factors)
- return true if (rem == 1)
- }
- }
- }
- return false
- }
- var I = gaussian_elimination(A, L-1)
- var LR = (A.end - A.rindex_by { !.is_zero })
- for solution in (I.last(LR)) {
- cfrac_find_factors(solution) && break
- }
- var final_factors = []
- for f in (factors) {
- if (f.is_prime) {
- final_factors << f
- }
- else {
- final_factors << __FUNC__(f)...
- }
- }
- if (rem != 1) {
- if (rem != n) {
- final_factors << __FUNC__(rem)...
- }
- else {
- final_factors << rem
- }
- }
- # Failed to factorize n (try again with a multiplier)
- if (rem == n) {
- return __FUNC__(n, next_multiplier(multiplier))
- }
- # Return all prime factors of n
- final_factors.sort
- }
- var composites = (ARGV ? ARGV.map{.to_i} : {|k| irand(2, 1<<k) }.map(2..60) )
- for n in (composites) {
- var factors = cffmm(n)
- assert_eq(factors.prod, n)
- say "#{n} = #{factors.map {|f| f.is_prime ? f : \"#{f} (composite)\"}.join(' * ')}"
- }
|