prog.sf

#!/usr/bin/ruby

# a(n) is the least integer m such that M(m) is divisible by prime(n)^2 or -1 if no such m exists.
# https://oeis.org/A354293

# See also:
#   https://oeis.org/A354292

# Known terms:
#   3, 4, -1, 23, 21, -1, 188, 65, 1010, 2231, -1, -1, 1326, 389, 1092, 13196, 1450, -1, 40466, 85553, 665, -1, 5139193, 333, -1, 408241, -1, 3072, 6702, 1393, 5832, 935, 1071, 77421, 292187, 775383, 493135, 4185, 1784560, 10632, 7935, 743003, 13418, 64499

var known_terms = Hash(
    1 => 3
2=> 4
5=> 21
4=> 23
8=> 65
7=> 188
24=> 333
14=> 389
21=> 665
9=> 1010
15=> 1092
13=> 1326
17=> 1450
10=> 2231
16=> 13196
19=> 40466
20=> 85553
26=> 408241
31 => 5832
32 => 935
33 => 1071
34 => 77421
35 => 292187
30 => 1393,
36 => 775383,
37 => 493135,
28 => 3072,
29 => 6702,
38 => 4185,

41 => 7935,
40 => 10632,
46 => 12176,
43 => 13418,
44 => 64499,

1 => 3
2 => 4
3 => -1
4 => 23
5 => 21
6 => -1
7 => 188
8 => 65
9 => 1010
10 => 2231
11 => -1
12 => -1
13 => 1326
14 => 389
15 => 1092

39 => 1784560,
42 => 743003,

# Term for p = 83 is given in the paper:
#   https://arxiv.org/pdf/2205.09902.pdf

83.pi => 5139193,

5.pi => -1,
13.pi => -1,
31.pi => -1,
37.pi => -1,
61.pi => -1,
79.pi => -1,
97.pi => -1,
103.pi => -1,
)

say "Known terms:\n"
var prev_k = 0
known_terms.keys.map{.to_i}.sort.each {|k|
    k - prev_k == 1 || break

    if (known_terms{k} != -1) {
        if (known_terms{k} < 1e6) {
            assert_eq(motzkin(known_terms{k}) % k.prime**2, 0)
        }
    }

    print(known_terms{k}, ", ")
    prev_k = k
}

say "\n"

#primes -= [  5, 13, 31, 37, 61, 79, 97, 103 ]
#var primes = [39.prime]
var primes = 200.primes

primes = primes.grep {|p|
    !known_terms.has(p.pi)
}

say primes

func catalanmod(n,m) {
    if (is_coprime(n+1, m)) {
        divmod(binomialmod(2*n, n, m), n+1, m)
    }
    else {
        idiv(binomial(2*n, n), (n+1)) % m
    }
}

func isok(n, m) {
    sum(0..n, {|k|
        (-1)**(n-k) * mulmod(binomialmod(n,k,m), catalanmod(k+1, m), m)
    }) % m == 0
}

#say isok(2231, 10.prime**2)
#say isok(5139193, 83**2)

var (a, b, n) = (0, 1, 1)

Inf.times {

    var M = idiv(b,n)
    n += 1
    (a, b) = (b, idiv((3*(n-1)*n*a + (2*n - 1)*n*b), ((n+1)*(n-1))))

    var found = []
    for p in (primes) {
        if (p.sqr `divides` M) {
            say "#{p.pi} => #{n-2}"
            found << p
        }
    }

    if (found) {
        primes -= found
        primes || break
    }
}