123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278 |
- #!/usr/bin/ruby
- # Daniel "Trizen" Șuteu
- # Date: 20 April 2019
- # https://github.com/trizen
- # Continued-fraction factorization method, combined with the mediant inequality to approximate sqrt(n).
- # See also:
- # https://en.wikipedia.org/wiki/Continued_fraction_factorization
- # https://en.wikipedia.org/wiki/Mediant_(mathematics)#Properties
- func gauss_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 farey_factorization (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 a1 = x
- var b1 = 1
- var a2 = x+1
- var b2 = 1
- var y = x
- var z = 1
- var w = x+x
- var r = w
- var (e1, e2) = (1, 0)
- var (f1, f2) = (0, 1)
- #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))
- #var factor_base = B.primes.grep {|p| (p <= 2) || (kronecker(N, p) >= 0) }
- var factor_base = (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
- }
- while (A.len < L) {
- # Continued fraction expansion of sqrt(N)
- y = (r*z - y)
- z = ((N - y*y) / z)
- r = round((x + y) / z)
- var a = ((x*f2 + e2) % n)
- var b = (a*a % n)
- var c = (b > w ? n-b : b)
- resolve_factor(a, c) if c.is_square
- if (c.is_smooth_over_prod(factor_prod)) {
- var c_factors = c.factor_exp
- if (c_factors) {
- A << exponent_signature(c_factors)
- Q << [a, c]
- }
- }
- (f1, f2) = (f2, (r*f2 + f1) % n)
- (e1, e2) = (e2, (r*e2 + e1) % n)
- if (z == 1) {
- return __FUNC__(n, next_multiplier(multiplier))
- }
- # Mediant inequality
- var a3 = a1+a2
- var b3 = b1+b2
- if (a3*a3 < N*b3*b3) {
- (a1, b1) = (a3, b3)
- }
- else {
- (a2, b2) = (a3, b3)
- }
- # There is a good chance that `numerator(m)^2 (mod n)` is B-smooth
- var t = (a3 % n)
- var v = (t*t % n)
- v = (n-v) if (n-v < v)
- # If `v` is a square, then `gcd(t - sqrt(v), n)` is divisor of `n`.
- resolve_factor(t, v) if v.is_square
- if (v.is_smooth_over_prod(factor_prod)) {
- var v_factors = v.factor_exp
- if (v_factors) {
- A << exponent_signature(v_factors)
- Q << [t, v]
- }
- }
- }
- if (A.len < L) {
- A += (L - A.end + 1 -> of(0))
- }
- var I = gauss_elimination(A, L-1)
- var LR = (A.end - A.rindex_by { !.is_zero })
- var rem = n
- var factors = []
- for solution in (I.last(LR)) {
- var X = 1
- var Y = 1
- var done = false
- 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)
- if (rem == 1) { done = true; break }
- }
- }
- }
- break if done
- }
- 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..40) )
- for n in (composites) {
- var factors = farey_factorization(n)
- assert_eq(factors.prod, n)
- say "#{n} = #{factors.map {|f| f.is_prime ? f : \"#{f} (composite)\"}.join(' * ')}"
- }
|