prog.sf

#!/usr/bin/ruby

# Smallest Niven (or Harshad) number (A005349) with exactly n distinct prime factors.
# https://oeis.org/A359960

# Known terms:
#   1, 2, 6, 30, 210, 2310, 30030, 690690, 14804790, 223092870

# New terms:
#   1, 2, 6, 30, 210, 2310, 30030, 690690, 14804790, 223092870, 8254436190, 200560490130, 8222980095330, 304250263527210, 13082761331670030, 614889782588491410, 32589158477190044730, 1987938667108592728530, 117288381359406970983270, 7858321551080267055879090, 573657473228859495079173570

#`( PARI/GP program:

omega_niven(A, B, n) = A=max(A, vecprod(primes(n))); (f(m, p, j) = my(list=List()); forprime(q=p, sqrtnint(B\m, j), my(v=m*q, r=nextprime(q+1)); while(v <= B, if(j==1, if(v>=A && v%sumdigits(v) == 0, listput(list, v)), if(v*r <= B, list=concat(list, f(v, r, j-1)))); v *= q)); list); vecsort(Vec(f(1, 2, n)));
a(n) = if(n==0, return(1)); my(x=vecprod(primes(n)), y=2*x); while(1, my(v=omega_niven(x, y, n)); 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.omega_prime_count(from, upto)

        if (count > 0) {

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

            n.omega_primes_each(from, upto, {|v|
                if (v.is_div(v.sumdigits)) {
                    say "a(#{n}) = #{v}"
                    return v
                }
            })
        }

        from = upto+1
        upto *= 2
    }
}

upper_bound(15)