OFFSET
1,1
COMMENTS
Subsequence of A058369.
If k is a term, then digsum(k) = 19, 37 or 73, for k < 10^9.
If k is an integer such that digsum(k) = digsum(k^2) = p, with p prime, then p == 1 (mod 9) (A061237).
This sequence has infinitely many terms of the form 1999...9 (A067272). If p is a prime with p == 1 (mod 9), i.e., p = 9m + 1 for some m, then t = 2*10^m - 1 = 1999...9, i.e., 1 followed by m 9's, is in this sequence since digsum(t) = 9m + 1 = p and t^2 = 39...960...01, where there are (m - 1) 9's and (m - 1) 0's, so digsum(t^2) = 3 + 9*(m - 1) + 6 + 1 = 9m + 1 = p. Dirichlet's theorem guarantees the existence of infinitely many primes of the form 9w + 1 and hence infinitely many terms of this sequence.
2*10^m - 1 is the least number with digit sum 9*m + 1. Since the next prime congruent to 1 (mod 9) after 73 is 109 = 9*12 + 1, the first term with digit sum other than 19, 37 or 73 is 2*10^12 - 1. - Robert Israel, Jul 07 2024
EXAMPLE
199 is a term, because its digital sum is 1 + 9 + 9 = 19 and 199^2 = 39601, whose digital sum is 3 + 9 + 6 + 0 + 1 = 19, which is prime.
MAPLE
ds:= n -> convert(convert(n, base, 10), `+`):
filter:= proc(n) local p;
p:= ds(n);
isprime(p) and ds(n^2) = p
end proc:
select(filter, [seq(i, i=1..1000, 9)]); # Robert Israel, Jul 05 2024
MATHEMATICA
Select[Range[5490], PrimeQ[dg=DigitSum[#]]&&(dg==DigitSum[#^2])&] (* Stefano Spezia, Jul 05 2024 *)
PROG
(PARI) isok(k) = my(s=sumdigits(k)); isprime(s) && (s==sumdigits(k^2)); \\ Michel Marcus, Jul 06 2024
CROSSREFS
KEYWORD
nonn,base
AUTHOR
Gonzalo MartÃnez, Jul 05 2024
STATUS
approved