experimental-projects/oeis-research/Felix Fröhlich/Square array A(n, k) read by antidiagonals downwards: smallest base-n Fermat pseudoprime with k distinct prime factors for k, n >= 2/prog.sf

#!/usr/bin/ruby

# Square array A(n, k) read by antidiagonals downwards: smallest base-n Fermat pseudoprime with k distinct prime factors for k, n >= 2.
# https://oeis.org/A271873

# Previously known terms:
#   341, 561, 91, 11305, 286, 15, 825265, 41041, 435, 124, 45593065, 825265, 11305, 561, 35

# New terms:
#   341, 561, 91, 11305, 286, 15, 825265, 41041, 435, 124, 45593065, 825265, 11305, 561, 35, 370851481, 130027051, 418285, 41041, 1105, 6, 38504389105, 2531091745, 30534805, 2203201, 25585, 561, 21, 7550611589521, 38504389105, 370851481, 68800501, 682465, 62745, 105, 28

#`{

# PARI/GP program:

fermat_psp(A, B, k, base) = A=max(A, vecprod(primes(k))); (f(m, l, lo, k) = my(list=List()); my(hi=sqrtnint(B\m, k)); if(lo > hi, return(list)); if(k==1, forstep(p=lift(1/Mod(m, l)), hi, l, if(isprimepower(p) && gcd(m*base, p) == 1, my(n=m*p); if(n >= A && (n-1) % znorder(Mod(base, p)) == 0, listput(list, n)))), forprime(p=lo, hi, base%p == 0 && next; my(z=znorder(Mod(base, p))); gcd(m,z) == 1 || next; my(q=p, v=m*p); while(v <= B, list=concat(list, f(v, lcm(l, z), p+1, k-1)); q *= p; Mod(base, q)^z == 1 || break; v *= p))); list); vecsort(Set(f(1, 1, 2, k)));
T(n,k) = if(n < 2, return()); my(x=vecprod(primes(k)), y=2*x); while(1, my(v=fermat_psp(x, y, k, n)); if(#v >= 1, return(v[1])); x=y+1; y=2*x);
print_table(n, k) = for(x=2, n, for(y=2, k, print1(T(x, y), ", ")); print(""));
for(k=2, 9, for(n=2, k, print1(T(n, k-n+2)", "))); \\ ~~~~

}

var base2_psp_with_n_prime_factors = Hash(
    2 341
    3 561
    4 11305
    5 825265
    6 45593065
    7 370851481
    8 38504389105
    9 7550611589521
    10 277960972890601
    11 32918038719446881
    12 1730865304568301265
    13 606395069520916762801
    14 59989606772480422038001
    15 6149883077429715389052001
    16 540513705778955131306570201
    17 35237869211718889547310642241
    18 13259400431578770557375638157521
    19 580827911915963785382253469211401
    20 292408776547176479576081475390697801
    21 39498823114155235923831808284152901601
    22 3284710806953570725820888298298841565601
    23 327373784481535488655521620744179013043601
    24 221404014114397213806721960178887462402717201
    25 43691666165877828056799483424028795272585383601
    26 13213974925373194568934435211730355813060799098001
    27 1952204134080476076724242017017925744953021675628161
    28 633922683896008820507659141940205847766668463614120801
    29 208615777921466463779477993429576353802684390620173966001
    30 44091058061027666417635106691235215695970713110442194459201
    31 2884167509593581480205474627684686008624483147814647841436801
    32 2214362930783528605057288166084711828471950070626477770522049201
)

var bfile_fh = File("bfile.txt").open_r
var lookup = Hash()

for k in (2..100) {
    for n in (2..k) {
        bfile_fh.eof && break
        var j = (k - n + 2)
        lookup{[n,j]} = bfile_fh.readline.nums.last
    }
}

func T(n, k) {

    if (n == 2) {
        return base2_psp_with_n_prime_factors{k} if base2_psp_with_n_prime_factors.has(k)
    }

    if (lookup.has([n,k])) {
        return lookup{[n,k]}
    }

    var x = pn_primorial(k)
    var y = 2*x

    loop {
        var arr = k.fermat_psp(n, x, y)

        if (arr) {
            return arr[0]
        }

        x = y+1
        y = 2*x
    }
}

var count = 2

for k in (2..100) {
    for n in (2..k) {
        say (count++, " ", T(n, k - n + 2))
    }
}