123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289 |
- #!/usr/bin/julia
- # A simple implementation of the continued fraction factorization method.
- # (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
- using Nemo
- const BIG_ONE = ZZ(1)
- const BIG_ZERO = ZZ(0)
- function next_multiplier(n, k)
- k += 2
- while (!is_squarefree(k) || gcd(k,n) != 1)
- k += 1
- end
- return k
- end
- function gaussian_elimination(A, n)
- m = length(A)
- I = [BIG_ONE << k for k in 0:m]
- nrow = 0
- for col in (1 : min(m, n))
- npivot = 0
- for row in (nrow+1 : m)
- if (((A[row] >> (col-1)) & BIG_ONE) == 1)
- npivot = row
- nrow += 1
- break
- end
- end
- npivot == 0 && continue
- if (npivot != nrow)
- A[[npivot,nrow]] = A[[nrow, npivot]]
- I[[npivot,nrow]] = I[[nrow, npivot]]
- end
- for row in (nrow+1 : m)
- if (((A[row] >> (col-1)) & BIG_ONE) == 1)
- A[row] = xor(A[row], A[nrow])
- I[row] = xor(I[row], I[nrow])
- end
- end
- end
- return I
- end
- function check_factor(n, g, factors)
- while (rem(n,g) == 0)
- n = div(n,g)
- push!(factors, g)
- if (is_prime(n))
- push!(factors, n)
- return 1
- end
- end
- return n
- end
- function is_smooth_over_prod(n, k)
- g = gcd(n, k)
- while (g > 1)
- n = div(n, g)
- while (rem(n, g) == 0)
- n = div(n, g)
- end
- n == 1 && return true
- g = gcd(n, g)
- end
- n == 1
- end
- function jacobi(a, n)
- kronecker_symbol(ZZ(a), ZZ(n))
- end
- function cffm(n, multiplier = 1)
- n <= 1 && return []
- is_prime(n) && return [n]
- if (iseven(n))
- v = 0
- while (iseven(n))
- v += 1
- n >>= 1
- end
- arr1 = [ZZ(2) for i in 1:v]
- arr2 = cffm(n)
- return append!(arr1, arr2)
- end
- if (is_square(n))
- f = cffm(isqrt(n))
- append!(f, f)
- return sort(f)
- end
- N = n*multiplier
- x = isqrt(N)
- y = x
- z = 1
- w = x+x
- r = w
- nf = round(Int64, exp(sqrt(log(n) * log(log(n))))^(sqrt(2) / 4) / 1.25)
- factor_base = [2]
- begin
- p = 3
- while (length(factor_base) < nf)
- if (jacobi(N, p) >= 0)
- push!(factor_base, p)
- end
- p = next_prime(p+1)
- end
- end
- factor_prod = prod([ZZ(k) for k in factor_base])
- factor_index = Dict{Int64, Int64}()
- for k in (1:length(factor_base))
- factor_index[factor_base[k]] = k-1
- end
- function exponent_signature(factors)
- sig = BIG_ZERO
- for (p,e) in factors
- if (isodd(e))
- sig |= (BIG_ONE << factor_index[p])
- end
- end
- return sig
- end
- L = length(factor_base)+1
- Q = []
- A = []
- (f1, f2) = (BIG_ONE, x)
- while (length(A) < L)
- y = (r*z - y)
- z = div((N - y*y), z)
- r = div((x + y), z)
- (f1, f2) = (f2, rem((r*f2 + f1), n))
- if (is_square(z))
- g = gcd(f1 - isqrt(z), n)
- if (g > 1 && g < n)
- arr1 = cffm(g)
- arr2 = cffm(div(n,g))
- return sort(append!(arr1, arr2))
- end
- end
- if (z > 1 && is_smooth_over_prod(z, factor_prod))
- push!(A, exponent_signature(factor(z)))
- push!(Q, [f1, z])
- end
- if (z == 1)
- println("z == 1 -> trying again with a multiplier...")
- return cffm(n, next_multiplier(n, multiplier))
- end
- end
- while (length(A) < L)
- push!(A, 0)
- end
- I = gaussian_elimination(A, length(A))
- LR = 0
- for k in (length(A):-1:1)
- if (A[k] != 0)
- LR = k+1
- break
- end
- end
- remainder = n
- factors = []
- function cfrac_find_factors(solution)
- X = BIG_ONE
- Y = BIG_ONE
- for i in 0:length(Q)-1
- if (((solution >> i) & BIG_ONE) == 1)
- X *= Q[i+1][1]
- Y *= Q[i+1][2]
- X %= remainder
- g = gcd(X - isqrt(Y), remainder)
- if (g > 1 && g < remainder)
- remainder = check_factor(remainder, g, factors)
- if (remainder == 1)
- return true
- end
- end
- end
- end
- return false
- end
- for k in LR:length(I)
- cfrac_find_factors(I[k]) && break
- end
- final_factors = []
- for f in (factors)
- if (is_prime(f))
- push!(final_factors, f)
- else
- append!(final_factors, cffm(f))
- end
- end
- if (remainder != 1)
- if (remainder != n)
- append!(final_factors, cffm(remainder))
- else
- push!(final_factors, remainder)
- end
- end
- # Failed to factorize n (try again with a multiplier)
- if (remainder == n)
- println("Trying again with a multiplier...")
- return cffm(n, next_multiplier(n, multiplier))
- end
- # Return all prime factors of n
- return sort(final_factors)
- end
- if (length(ARGS) >= 1)
- for n in (ARGS)
- println(n, " = ", cffm(ZZ(n)))
- end
- else
- @time println("2^128 + 1 = ", cffm(ZZ(2)^128 + 1))
- end
|