sidef-scripts/Math/continued_fraction_factorization_method.sf

#!/usr/bin/ruby

# Daniel "Trizen" Șuteu
# Date: 25 October 2018
# https://github.com/trizen

# Simple implementation of the continued fraction factorization method (CFRAC),
# combined with modular arithmetic (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

define(
    B_SMOOTH_METHOD = false,    # true to use the B-smooth formula for the factor base
    ROUND_DIVISION  = true,     # true to use round division instead of floor division
)

func gaussian_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 cffmm (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 y = x
    var z = 1
    var w = x+x
    var r = w

    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) / 1.25)

    var factor_base = if (B_SMOOTH_METHOD) {
        B.primes.grep {|p| (p <= 2) || (kronecker(N, p) >= 0) }
    }
    else {
        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
    }

    var (f1, f2) = (1, x)

    while (A.len < L) {

        y = (r*z - y)
        z = idiv(N - y*y, z)         # exact division
        r = (ROUND_DIVISION ? idiv_round(x+y, z) : idiv(x+y, z))

        (f1, f2) = (f2, addmod(mulmod(r, f2, n), f1, n))

        if (z.is_square) {
            resolve_factor(f1, z)
        }

        if (z.abs.is_smooth_over_prod(factor_prod)) {
            var c_factors = z.abs.factor_exp

            if (c_factors) {
                A << exponent_signature(c_factors)
                Q << [f1, z.abs]
            }
        }

        if (z.abs == 1) {
            return __FUNC__(n, next_multiplier(multiplier))
        }
    }

    if (A.len < L) {
        A += (L - A.end + 1 -> of(0))
    }

    var rem     = n
    var factors = []

    func cfrac_find_factors(solution) {

        var X = 1
        var Y = 1

        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)
                    return true if (rem == 1)
                }
            }
        }

        return false
    }

    var I  = gaussian_elimination(A, L-1)
    var LR = (A.end - A.rindex_by { !.is_zero })

    for solution in (I.last(LR)) {
        cfrac_find_factors(solution) && break
    }

    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..60) )

for n in (composites) {

    var factors = cffmm(n)
    assert_eq(factors.prod, n)

    say "#{n} = #{factors.map {|f| f.is_prime ? f : \"#{f} (composite)\"}.join(' * ')}"
}