123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990 |
- #!/usr/bin/julia
- # The elliptic-curve factorization method (ECM), due to Hendrik Lenstra.
- # Algorithm presented in the YouTube video:
- # https://www.youtube.com/watch?v=2JlpeQWtGH8
- # See also:
- # https://en.wikipedia.org/wiki/Lenstra_elliptic-curve_factorization
- using Primes
- function ecm(N, arange=100, plimit=10000)
- isprime(N) && return N
- # TODO: make sure that N is not a perfect power.
- P = primes(plimit)
- for a in (-arange : arange)
- x = 0
- y = 1
- for B in P
- t = B^trunc(Int64, log(B, plimit))
- (xn, yn) = (0, 0)
- (sx, sy) = (x, y)
- first = true
- while (t > 0)
- if (isodd(t))
- if (first)
- (xn, yn) = (sx, sy)
- first = false
- else
- g = gcd(sx-xn, N)
- if (g > 1)
- g == N && break
- return g
- end
- u = invmod(sx-xn, N)
- L = ((u*(sy-yn)) % N)
- xs = ((L*L - xn - sx) % N)
- yn = ((L*(xn - xs) - yn) % N)
- xn = xs
- end
- end
- g = gcd(2*sy, N) # TODO: if N is odd, use gcd(sy, N) instead
- if (g > 1)
- g == N && break
- return g
- end
- u = invmod(2*sy, N)
- L = ((u*(3*sx*sx + a)) % N)
- x2 = ((L*L - 2*sx) % N)
- sy = ((L*(sx - x2) - sy) % N)
- sx = x2
- sy == 0 && return 1
- t >>= 1
- end
- (x, y) = (xn, yn)
- end
- end
- return 1
- end
- if (length(ARGS) >= 1)
- for n in (ARGS)
- println("One factor of ", n, " is: ", ecm(parse(BigInt, n)))
- end
- else
- @time println("One factor of 2^128 + 1 is: ", ecm(big"2"^128 + 1))
- end
|