123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105 |
- #!/usr/bin/ruby
- # A simple version of the continued fraction factorization method, without the linear algebra stage.
- # Similar to solving the Pell equation:
- # x^2 - d*y^2 = 1, where `d` is known.
- # See also:
- # https://en.wikipedia.org/wiki/Pell%27s_equation
- func cffm (n) {
- # 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)
- }
- # Check for divisibility by 2
- if (n.is_even) {
- var v = n.valuation(2)
- var t = n>>v
- var factors = v.of(2)
- if (t > 1) {
- factors += __FUNC__(t)
- }
- return factors
- }
- 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 x = n.isqrt
- var y = x
- var z = 1
- var r = x+x
- var (e1, e2) = (1, 0)
- var (f1, f2) = (0, 1)
- var S = Hash()
- do {
- y = (r*z - y)
- z = idiv(n - y*y, z)
- r = idiv(x + y, z)
- var a = addmod(x*f2, e2, n)
- var b = mulmod(a, a, n)
- if (b.is_square) {
- resolve_factor(a, b)
- }
- if (S.exists(b)) {
- resolve_factor(a * S{b}, b*b)
- }
- else {
- S{b} = a
- }
- (f1, f2) = (f2, addmod(r*f2, f1, n))
- (e1, e2) = (e2, addmod(r*e2, e1, n))
- } while (z != 1)
- return [n]
- }
- var composites = (ARGV ? ARGV.map{.to_i} : {|k| 2.of { random_prime(1<<k) }.prod }.map(2..25) )
- for n in (composites) {
- var factors = cffm(n)
- assert_eq(factors.prod, n)
- if (factors.all { .is_prime }) {
- say "#{n} = #{factors.join(' * ')}"
- }
- else {
- say "#{n} = #{factors.join(' * ')} (incomplete factorization)"
- }
- }
|