OFFSET
1,1
COMMENTS
There are (1, 10, 23, 43, 87, 146, 255, 408, 642) terms below (1, 2, 3, ..., 9)*10^6. Among the 960 terms below 10^7, only {8998991, 8999981, 9899891, 9988991, 9997991, 9999971} end in the digit 1, only 30 end in the digit 3, and 268 end with a digit 7. - M. F. Hasler, Mar 09 2022
LINKS
Michael S. Branicky, Table of n, a(n) for n = 1..10000 (terms 1..1000 from Vincenzo Librandi)
MATHEMATICA
Select[Prime[Range[270000]], Total[IntegerDigits[#]]==53&] (* Harvey P. Dale, Apr 17 2011 *)
PROG
(Magma) [p: p in PrimesUpTo(3800000) | &+Intseq(p) eq 53]; // Vincenzo Librandi, Jul 09 2014
(PARI) select( {is_A106782(p)=sumdigits(p)==53&&isprime(p)}, primes([9e5, 4e6]\1)) \\ M. F. Hasler, Mar 09 2022
(Python)
from itertools import count, islice
from sympy import isprime
from sympy.utilities.iterables import multiset_permutations
def agen(b=10, sod=53): # generator for any base, sum-of-digits
if 0 <= sod < b:
yield sod
nzdigs = [i for i in range(1, b) if i <= sod]
nzmultiset = []
for d in range(1, b):
nzmultiset += [d]*(sod//d)
for d in count(2):
fullmultiset = [0]*(d-1-(sod-1)//(b-1)) + nzmultiset
for firstdig in nzdigs:
target_sum, restmultiset = sod - int(firstdig), fullmultiset[:]
restmultiset.remove(firstdig)
for p in multiset_permutations(restmultiset, d-1):
if sum(p) == target_sum:
t = int("".join(map(str, [firstdig]+p)), b)
if isprime(t):
yield t
if p[0] == target_sum:
break
print(list(islice(agen(), 25))) # Michael S. Branicky, Mar 10 2022
CROSSREFS
KEYWORD
nonn,base
AUTHOR
Zak Seidov, May 16 2005
STATUS
approved