sidef-scripts/Math/chebyshev_factorization_method.sf

#!/usr/bin/ruby

# Daniel "Trizen" Șuteu
# Date: 11 December 2019
# https://github.com/trizen

# A simple factorization method, using the Chebyshev T_n(x) polynomials, based on the identity:
#   T_{m n}(x) = T_m(T_n(x))

# where:
#   T_n(x) = (1/2) * V_n(2x, 1)

# where V_n(P, Q) is the Lucas V sequence.

# See also:
#   https://oeis.org/A001075
#   https://en.wikipedia.org/wiki/Lucas_sequence
#   https://en.wikipedia.org/wiki/Iterated_function
#   https://en.wikipedia.org/wiki/Chebyshev_polynomials

func lucasVmod(P, n, m) {    # assumes Q = 1

    var(V1 = 2, V2 = P)

    for bit in (n.bits) {
        if (bit) {
            V1 = submod(mulmod(V2, V1, m), P, m)
            V2 = mulmod(V2, V2, m)-2
        }
        else {
            V2 = submod(mulmod(V2, V1, m), P, m)
            V1 = mulmod(V1, V1, m)-2
        }
    }

    return V1
}

func chebyshev_factorization(n, B1 = n.ilog(6)**3, a = 3) {

    var x = a
    var G = B1*B1
    var i = invmod(2, n)

    func chebyshevTmod(a, x) {
        mulmod(lucasVmod(2*x, a, n), i, n)
    }

    primes_each(2, B1.isqrt, {|p|
        G.ilog(p).times {
            x = chebyshevTmod(p, x)      # T_k(x) (mod n)
        }
    })

    primes_each(B1.isqrt+1, B1, {|p|
        x = chebyshevTmod(p, x)      # T_k(x) (mod n)
        is_coprime(x-1, n) || return gcd(x-1, n)
    })

    return gcd(x-1, n)
}

say chebyshev_factorization(257221 * 470783,               1000);     #=> 257221           (p+1 is 1000-smooth)
say chebyshev_factorization(1124075136413 * 3556516507813, 4000);     #=> 1124075136413    (p+1 is 4000-smooth)
say chebyshev_factorization(7553377229 * 588103349,         800);     #=> 7553377229       (p+1 is  800-smooth)
say chebyshev_factorization(333732865481 * 1632480277613,  3000);     #=> 333732865481     (p-1 is 3000-smooth)

say ''

say chebyshev_factorization(15597344393 * 12388291753,   3000)        #=> 15597344393       (p-1 is 3000-smooth)
say chebyshev_factorization(43759958467 * 59037829639,   3200)        #=> 43759958467       (p+1 is 3200-smooth)
say chebyshev_factorization(112601635303 * 83979783007,   700)        #=> 112601635303      (p-1 is  700-smooth)
say chebyshev_factorization(228640480273 * 224774973299, 2000)        #=> 228640480273      (p-1 is 2000-smooth)

say ''

say chebyshev_factorization(5140059121 * 8382882743,     2500, 2)     #=> 5140059121      (p-1 is 2500-smooth)
say chebyshev_factorization(18114813019 * 17402508649,   6000, 2)     #=> 18114813019     (p+1 is 6000-smooth)
say chebyshev_factorization(533091092393 * 440050095029,  300, 2)     #=> 533091092393    (p+1 is  300-smooth)

say ''

for k in (3..20) {

    var n = 2.of { random_nbit_prime(k) }.prod
    var p = chebyshev_factorization(n, 2 * n.ilog2**2)

    if (p.is_prime) {
        say "#{n} = #{p} * #{n/p}"
    }
}

__END__
143 = 11 * 13
323 = 17 * 19
2773 = 59 * 47
12091 = 107 * 113
42781 = 179 * 239
181427 = 433 * 419
685603 = 967 * 709
3240379 = 1999 * 1621
12254063 = 3089 * 3967
38271817 = 6763 * 5659
119391583 = 9511 * 12553
367073731 = 20929 * 17539
3578349739 = 60961 * 58699
15888478703 = 121571 * 130693
36182994919 = 143593 * 251983
545782820497 = 668713 * 816169
3656558375893 = 1912639 * 1911787
12346497466267 = 3737821 * 3303127
31369209076123 = 5655031 * 5547133
148520973663353 = 13394431 * 11088263
412638129424819 = 17334017 * 23805107