A355317 - OEIS

%I #48 Jul 07 2022 08:03:42

%S 5701,27091,50417,55049,57089,60601,61051,63607,66301,72019,132607,

%T 270913,284057,574031,654301,936007,936907,950647,1250609,1461001,

%U 1504417,1580921,1682069,1703287,1750631,1810553,1810573,1837601,2050241,2145089,2485001,2641109,2709169,2800333,2805703

%N Prime numbers that have the same base-10 digits as their prime index (with multiplicity), disregarding zero digits.

%e 15097067 is a term since its index primepi(15097067) = 976571 shares the same multiset of nonzero digits {1,5,6,7,7,9}.

%t a[max_] := Module[{l}, Select[{#, Prime[#]} & /@ Range[max], (l = IntegerDigits[#[[2]]]; SortBy[Tally[l], First] === SortBy[Tally[PadLeft[IntegerDigits[#[[1]]], Length[l]]], First]) &]]; a[10^6][[All, 2]] (* Gives the first 108 terms *)

%o (Python)

%o from sympy import nextprime

%o from itertools import islice

%o def b10s(n): return "".join(sorted(str(n))).lstrip("0")

%o def agen():

%o     k, pk = 1, 2

%o     while True:

%o         if b10s(k) == b10s(pk): yield pk

%o         k, pk = k+1, nextprime(pk)

%o print(list(islice(agen(), 32))) # _Michael S. Branicky_, Jun 28 2022

%o (PARI) digs(k) = vecsort(select(x->(x>0), digits(k)));

%o isok(p) = if (isprime(p), digs(p) == digs(primepi(p))); \\ _Michel Marcus_, Jul 05 2022

%Y Cf. A355318 (the prime indices), A355539.

%Y Cf. A355418.

%K base,easy,nonn

%O 1,1

%A _Xiaofeng Wang_, Jun 28 2022