# Python program for OEIS A346499
# Michael S. Branicky, Sep 26 2021

# A346499 Prime numbers k where the d(j)-th digit is j for d(j) and j > 0 and d(j) = 0 if and only if j is not a digit of k.
data = [103, 2143, 4001, 4201, 4231, 12503, 21503, 50341, 52301, 53201, 106453, 106543, 146203, 426103, 600401, 600451, 603401, 603541, 620051, 620401, 623401, 1000507, 1004567, 1004657, 1005467, 1006547, 1030067, 1034567, 1035467, 1043657, 1050307, 1050367, 1056347, 1070063, 1075463, 1200007, 1204507, 1205467, 1230067, 1234067, 1234657, 1254367, 1256347, 1264037, 1264537, 1270063, 1276543, 1324007, 1324567, 1326047, 1402567, 1462037, 1504267, 1534207, 1543207, 1564237, 1576243, 1605427, 1634027, 1643027, 1674523, 2104507, 2104567, 2104657, 2105407, 2134007, 2134507, 2135407, 2143567, 2154307, 2164037, 2174063, 2176543, 3010067]

from sympy import isprime
from itertools import combinations
from sympy.utilities.iterables import generate_involutions

def afind(maxdigits):
    alst = []
    for N in range(3, min(9, maxdigits)+1):
        testlst = []
        for inv in generate_involutions(N):
            digs = [d + 1 for d in inv]
            if int(max(digs)) > len(digs): continue
            for j, dj in enumerate(digs, start=1):
                if dj and digs[dj-1] != j:
                    break
            else:
                testlst.append(digs)
                fixed = [i for i in range(2, len(digs)+1) if digs[i-1] == i]
                for zs in range(1, len(fixed)+1):
                    for zlocs in combinations(fixed, zs):
                        t = [di*(i+1 not in zlocs) for i, di in enumerate(digs)]
                        testlst.append(t)
        testlst = sorted([int("".join(map(str, t))) for t in testlst])
        alst.extend([t for t in testlst if isprime(t)])
    return alst

# TEST ON DATA
print(afind(7))
print(data)
print()

# TIME PROGRAM
for N in range(3, 10):
    %time ans = afind(N)
    print(N, len(ans))

# WRITE BFILE FOR ENTIRE SEQUENCE
with open('b346499.txt', 'w') as bfile:
    for n, an in enumerate(ans, start=1):
        bfile.write(f"{n} {an}\n")