prog.sf

#!/usr/bin/ruby

# Smallest prime p such that p^(2n-1) - 1 is the product of 2n-1 distinct primes.
# https://oeis.org/A359069

# Known terms:
#   3, 59, 47, 79, 347, 6343, 56711, 4523

# Lower-bounds:
#   a(9)  > 1310807
#   a(10) > 467611
#   a(11) > 671219

include("../../../factordb/auto.sf")

Num!VERBOSE = true

func check_partial_factors(f,n) {

    f.uniq.len == f.len || return false
    f.sum {|p| p.is_prime ? 1 : 2 } > n && return false

    if (f.all_prime) {
        if (f.len == n) {
            return true
        }
        return false
    }

    return true
}

func a(n, from=2) {
    assert(from.is_prime)

    for (var k = from; true; k.next_prime!) {

        k.dec.is_squarefree || next

        var v = (k**n - 1)
        v.is_prob_squarefree(1e6) || next

        say "[#{n}] Checking: #{k}"

        var tf = v.trial_factor(1e6)
        check_partial_factors(tf, n) || next
        tf.len.dec + tf.last.ilog(1e6) + 1 >= n || next

        var tf = v.trial_factor(1e7)
        check_partial_factors(tf, n) || next
        tf.len.dec + tf.last.ilog(1e7) + 1 >= n || next

        if (tf.last > 1e60) {
            tf = v.trial_factor(1e8)
            check_partial_factors(tf, n) || next
            tf.len.dec + tf.last.ilog(1e8) + 1 >= n || next
        }

        say "Many factors (at least #{tf.len-1 + (tf.last.is_prime ? 1 : 2)} with C#{tf.last.len}): #{v}"

        var ff = v.special_factor
        check_partial_factors(ff, n) || next

        var f = if (v % (k-1) == 0) {
            factordb(v / (k-1))
        }
        else {
            factordb(v)
        }

        var f3 = gcd_factors(v, tf + ff + f)
        check_partial_factors(f3, n) || next

        if ((f3.last > 1e65) && f3.last.is_composite) {
            tf = v.trial_factor(1e9)
            check_partial_factors(tf, n) || next
            tf.len.dec + tf.last.ilog(1e9) + 1 >= n || next

            say "Strong candidate..."

            f3 = gcd_factors(v, tf + ff + f)
            check_partial_factors(f3, n) || next

            var pf = f3.grep{.is_prime}
            var c = (v / pf.prod)
            pf.len + c.ilog(1e9) + 1 >= n || next

            say "Factoring C#{c.len}: #{c}"
        }

        f3 = f3.map{ .is_prime ? _ : factordb(_) }.flat
        check_partial_factors(f3, n) || next

        f3 = f3.map{ .factor }.flat
        check_partial_factors(f3, n) || next

        if (f3.len == n) {
            return k
        }
    }
}

#var n = 9
#var from = 1310807.next_prime

var n = 11
var from = 671219.prev_prime

say "a(#{n}) = #{a(2*n - 1, from)}"