prog.sf

#!/usr/bin/ruby

# Smallest squarefree triangular number with exactly n prime factors.
# https://oeis.org/A127637

# Known terms:
#   1, 3, 6, 66, 210, 3570, 207690, 930930, 56812470, 1803571770, 32395433070, 265257422430, 91348974206490, 24630635909489610, 438603767516904990, 14193386885746698630, 2378522762792139793830, 351206814022419685159830, 28791787439593010836313310

# New terms:
#   a(19) = 2402835013540065887743290330
#   a(20) = 120052594044654305809137933570
#   a(21) = 43869525454581224791547259014910

# 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); list=concat(list, f(t, q+1, j-1)))); list); vecsort(Vec(f(1, 2, n)));
a(n, k=3) = if(n==0, return(1)); 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); \\ ~~~~

)

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(3)) {
                    say "a(#{n}) = #{v}"
                    return v
                }
            })
        }

        from = upto+1
        upto *= 2
    }
}

upper_bound(9)