is_prob_squarefree.sf

#!/usr/bin/ruby

# Daniel "Trizen" Șuteu
# Date: 09 January 2019
# https://github.com/trizen

# Test if a large integer (>50 digits) is probably squarefree.

# See also:
#   https://en.wikipedia.org/wiki/Square-free_integer

func is_squarefree_over_prod (n, k) {

    var g = gcd(n, k)

    if (g > 1) {
        var r = (n / g)
        return true if (r == 1)
        return false if (gcd(r, k) > 1)
        return false if r.is_power
    }

    return true
}

func store_random_factor(a, f) {
    *a = f()
    (*a).len > 1
}

func find_random_factor(n, small=true, trial=false) {

    var e = []

    (trial && store_random_factor(\e, { n.trial_factor(1e6) })) ||
    (small && store_random_factor(\e, { n.pbrent_factor(1, 15000) })) ||
    (small && store_random_factor(\e, { n.pminus1_factor(50000, 500000) })) ||

    store_random_factor(\e, { n.pminus1_factor(1e6) })  ||
    store_random_factor(\e, { n.pplus1_factor(1e6) })   ||
    store_random_factor(\e, { n.ecm_factor(800, 10) })  ||
    store_random_factor(\e, { n.ecm_factor(8000, 20) }) ||

    #store_random_factor(\e, { n.ecm_factor(80000, 40) }) ||
    #store_random_factor(\e, { n.ecm_factor(320000, 40) }) ||
    #store_random_factor(\e, { n.ecm_factor(1000000, 80) }) ||

    store_random_factor(\e, { n.holf_factor(1e6) }) ||
    store_random_factor(\e, { n.squfof_factor(1e8) });

    return e
}

func is_prob_squarefree(r, tries=50) {

    if (r.len <= 50) {
        return r.is_squarefree
    }

    return false if r.is_power

    static FP = 5.of {|k|
        10**(k+2) -> primorial
    }

    FP.all { is_squarefree_over_prod(r, _) } || return false

    loop {
        var f = r.trial_factor(1e6)

        break if (f.len == 1)

        r /= f[0]
        return false if (f[0] `divides` r)
    }

    tries.times { |k|

        var len = r.len

        if (len <= 50) {
            return r.is_squarefree
        }

        if (len <= 10_000) {
            return true if r.is_prob_prime
        }

        var e = find_random_factor(r, k <= 3)

        e.len > 1 || break
        e.pop

        e.each {|p|

            if (p*p `divides` r) {
                return false
            }

            p.is_squarefree || return false

            p.factor.each { |z|
                if (z*z `divides` r) {
                    return false
                }
            }
        }

        e.each {|p| r /= p }

        if (r.is_power) {
            return false
        }
    }

    return true
}

# Example for testing several k such that 2^k + 1 is divisible by a square > 1.
# See: https://oeis.org/A049096

var a = [231, 234, 237, 243, 410, 411, 417, 891, 897, 903, 2370, 2373, 41390]

a.each {|k|
    var t = is_prob_squarefree(2**k + 1)
    say "2^#{k} + 1 is #{t ? 'prob squarefree' : 'divisible by a square'}"
}

say "\nLarge non-squarefree tests:"
say is_prob_squarefree(249860117627119815706149728420771785501032147277376706113)              # false
say is_prob_squarefree(1444134648503819781988817515640867565630266453262667463566673)          # false
say is_prob_squarefree(16240879415462221539256806220654226179619821855946633706672633593)      # false
say is_prob_squarefree(300973095636437173473855290756399701276460223130678583012549023809113)  # false