generate.sf

#!/usr/bin/ruby

# a(n) is the least n-gonal number that is the product of n distinct primes, or 0 if there are none.
# https://oeis.org/A359854

# Known terms:
#   6, 66, 0, 11310, 303810, 28962934, 557221665, 15529888374, 1219300152070

# a(n) >= A358862(n), for n >= 3. - ~~~~

# PARI/GP program:
#`(

# General version
squarefree_omega_polygonals(A, B, n, k) = A=max(A, vecprod(primes(n))); (f(m, p, j) = my(list=List()); my(s=sqrtnint(B\m, j)); if(j==1, forprime(q=max(p, ceil(A/m)), s, if(ispolygonal(m*q, k), listput(list, m*q))), forprime(q=p, s, my(t=m*q); if(ceil(A/t) <= B\t, list=concat(list, f(t, q+1, j-1))))); list); vecsort(Vec(f(1, 2, n)));
a(n, k=n) = if(n < 2, return()); if(n==2, return(6)); if(n==4, return(0)); my(x=vecprod(primes(n)), y=2*x); while(1, my(v=squarefree_omega_polygonals(x, y, n, k)); if(#v >= 1, return(v[1])); x=y+1; y=2*x); \\ ~~~~

# Optimized using residue primes:
squarefree_omega_polygonals(A, B, n, k) = A=max(A, vecprod(primes(n))); (f(m, p, j) = my(list=List()); my(s=sqrtnint(B\m, j)); if(j==1, forprime(q=max(p, ceil(A/m)), s, if(ispolygonal(m*q, k), listput(list, m*q))), forprime(q=p, s, my(r=q%k); if(r == 1 || r == 2 || r == 3 || r == 5 || r == 7 || r == 9 || r == 11 || r == 13 || r == 15, my(t=m*q); if(ceil(A/t) <= B\t, list=concat(list, f(t, q+1, j-1)))))); list); vecsort(Vec(f(1, 2, n)));
a(n, k=n) = if(n < 2, return()); if(n==2, return(6)); if(n==4, return(0)); my(x=vecprod(primes(n)), y=2*x); while(1, my(v=squarefree_omega_polygonals(x, y, n, k)); if(#v >= 1, return(v[1])); x=y+1; y=2*x); \\ ~~~~

)

# Find the residues of the prime factors of n-gonal numbers:
#   for n in (3..20) { say [n, 1000.of {|k| polygonal(k, n) }.map{.factor.map{.mod(n)}}.flat.uniq.sort] }

# Resides for each n = 3..n:
#`(
    [3, [0, 1, 2]]
    [4, [1, 2, 3]]
    [5, [0, 1, 2, 3, 4]]
    [6, [1, 2, 3, 5]]
    [7, [0, 1, 2, 3, 4, 5, 6]]
    [8, [1, 2, 3, 5, 7]]
    [9, [1, 2, 3, 4, 5, 7, 8]]
    [10, [1, 2, 3, 5, 7, 9]]
    [11, [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]]
    [12, [1, 2, 3, 5, 7, 11]]
    [13, [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]]
    [14, [1, 2, 3, 5, 7, 9, 11, 13]]
    [15, [1, 2, 3, 4, 5, 7, 8, 11, 13, 14]]
    [16, [1, 2, 3, 5, 7, 9, 11, 13, 15]]
    [17, [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]]
    [18, [1, 2, 3, 5, 7, 11, 13, 17]]
    [19, [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]]
    [20, [1, 2, 3, 5, 7, 9, 11, 13, 17, 19]]
)

func upper_bound(n, from = 2, upto = 2*from) {

    say "\n:: Searching an upper-bound for a(#{n})\n"

    loop {

        var count = n.squarefree_almost_prime_count(from, upto)

        if (count > 0) {

            say "Sieving range: [#{from}, #{upto}]"
            say "This range contains: #{count.commify} elements\n"

            n.squarefree_almost_primes_each(from, upto, {|v|
                if (v.is_polygonal(n)) {
                    say "a(#{n}) = #{v}"
                    return v
                }
            })
        }

        from = upto+1
        upto *= 2
    }
}

upper_bound(13)

__END__
a(2) = 6
a(3) = 66
a(4) = 0
a(5) = 11310
a(6) = 303810
a(7) = 28962934
a(8) = 557221665
a(9) = 15529888374
a(10) = 1219300152070
a(11) = 23900058257790
a(12) = 1231931106828345
a(13) = 500402553453949510
a(14) = 14990069451769732194
a(15) = 610385355391371697410

6, 66, 0, 11310, 303810, 28962934, 557221665, 15529888374, 1219300152070, 23900058257790, 1231931106828345, 500402553453949510, 14990069451769732194, 610385355391371697410