fermat_hybrid_factorization_method.sf

#!/usr/bin/ruby

# Daniel "Trizen" Șuteu
# Date: 27 February 2018
# https://github.com/trizen

# A hybrid factorization algorithm, using:
#   * Pollard's p-1 algorithm
#   * Pollard's rho algorithm
#   * A simple version of the continued-fraction factorization method
#   * Fermat's factorization method

# See also:
#   https://en.wikipedia.org/wiki/Quadratic_sieve
#   https://en.wikipedia.org/wiki/Dixon%27s_factorization_method
#   https://en.wikipedia.org/wiki/Fermat%27s_factorization_method
#   https://en.wikipedia.org/wiki/Pollard%27s_p_%E2%88%92_1_algorithm

func fermat_hybrid_factorization (n) {

    if ((n <= 1) or n.is_prime) {
        return [n]
    }

    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 p = n.isqrt
    var x = p
    var q = (p*p - n)

    var t = 1
    var u = 1
    var v = 1

    var h = 1
    var k = 4
    var z = n.random_prime

    var g = 1
    var c = q+p

    var a0 = 1
    var a1 = (a0*a0 + c)
    var a2 = (a1*a1 + c)

    var c1 = p
    var c2 = 1

    var r = p+p

    var (e1, e2) = (1, 0)
    var (f1, f2) = (0, 1)

    while (!q.is_square) {

        # Trizen's random factorization method (for small factors)
        u = powmod(q * u, k, n)
        v = powmod(p * v, k, n)

        k += 16

        g = gcd(u - v, n)

        if ((g > 1) && (g < n)) {
            return (
                __FUNC__(n/g) +
                __FUNC__(g) -> sort
            )
        }

        q += ((p++ << 1) + 1)

        # Pollard's rho algorithm
        g = gcd(n, a2-a1)

        if ((g > 1) && (g < n)) {
            return (
                __FUNC__(n/g) +
                __FUNC__(g) -> sort
            )
        }

        a1 = ((a1*a1 + c) % n)
        a2 = ((a2*a2 + c) % n)
        a2 = ((a2*a2 + c) % n)

        # Simple version of the continued-fraction factorization method
        c1 = (r*c2 - c1)
        c2 = ((n - c1*c1) // c2)        #/

        var a = ((x*f2 + e2) % n)
        var b = (a*a % n)

        if (b.is_square) {
            g = gcd(a - b.isqrt, n)

            if ((g > 1) && (g < n)) {
                return (
                    __FUNC__(g) +
                    __FUNC__(n/g) -> sort
                )
            }
        }

        r = round((x + c1) / c2)

        (f1, f2) = (f2, (r*f2 + f1) % n)
        (e1, e2) = (e2, (r*e2 + e1) % n)

        # Pollard's p-a algorithm (random variation)
        t = z
        h = h.next_prime
        z = powmod(z, h, n)
        g = gcd(z * powmod(t, n.irand, n) - 1, n)

        if (g > 1) {

            if (g == n) {
                h = 1
                z = n.random_prime
                next
            }

            return (
                __FUNC__(n/g) +
                __FUNC__(g) -> sort
            )
        }
    }

    # Fermat's method
    var s = q.isqrt

    __FUNC__(p + s) +
    __FUNC__(p - s) -> sort
}

if (ARGV) {
    say fermat_hybrid_factorization(ARGV[0].to_i)
    return 1
}

for k in (1..63) {
    var n = irand(2, 1<<k)
    var prime_factors = fermat_hybrid_factorization(n)
    say "#{n} = #{prime_factors.join(' * ')}"
    assert_eq(prime_factors.prod, n)
}