A353126 - OEIS

%I #25 May 07 2022 10:31:55

%S 2,2700,40071,100363

%N Numbers k such that k = PrimePi(k * digsum(k)).

%C Numbers k such that k = A000720(A057147(k)).

%C No further terms < 3600000. - _Michael S. Branicky_, Apr 25 2022

%C a(5) > 10^10, if it exists. - _David A. Corneth_, Apr 25 2022

%e Consider number 2: digsum(2) equal 2; 2*2 = 4; The number of primes not exceeding 4 is 2, the number itself. Thus, 2 is in this sequence.

%e Consider number 11: 11*digsum(11) = 22; PrimePi(22) = 8. Thus, 11 is not in this sequence.

%t Select[Range[1000000], # == PrimePi[# Total[IntegerDigits[#]]] &]

%o (PARI) isok(k) = k==primepi(k*sumdigits(k)); \\ _Michel Marcus_, Apr 25 2022

%o (PARI) upto(n) = {q = 2; t = 0; res = List(); forprime(p = 3, n, t++; s = nextmultiple(q, t); forstep(i = s, p - 1, t, if(i % t == 0, c = i/t; if(sumdigits(t) == c, listput(res, t)); ) ); q = p; ); res }

%o nextmultiple(n, m) = my(d = ((n-m)%m)); n + !!d*m - d \\ _David A. Corneth_, Apr 25 2022

%o (Python)

%o from sympy import primepi, sieve

%o sieve.extend(12*10**6)

%o def sod(n): return sum(map(int, str(n)))

%o def afind(limit):

%o     for k in range(1, limit+1):

%o         if k == primepi(k*sod(k)):

%o             print(k, end=", ")

%o afind(100363) # _Michael S. Branicky_, Apr 25 2022

%Y Cf. A000720, A007953, A057147, A353134.

%K nonn,base,more

%O 1,1

%A _Tanya Khovanova_, Apr 24 2022