%I #36 Feb 22 2025 18:11:30
%S 1296,10000,38416,50625,194481,234256,456976,1185921,1336336,1500625,
%T 2085136,2313441,4477456,6765201,9150625,10556001,11316496,14776336,
%U 16777216,17850625,22667121,29986576,35153041,45212176,52200625
%N Numbers with 25 divisors.
%C Maple implementation: see A030513.
%C Numbers of the form p^24 (24th powers of A000040, subset of A010812) or p^4*q^4 (A189991), where p and q are distinct primes. - _R. J. Mathar_, Mar 01 2010
%H T. D. Noe, <a href="/A137488/b137488.txt">Table of n, a(n) for n = 1..1000</a>
%F A000005(a(n)) = 25.
%F Sum_{n>=1} 1/a(n) = (P(4)^2 - P(8))/2 + P(24) = 0.000933328..., where P is the prime zeta function. - _Amiram Eldar_, Jul 03 2022
%t lst = {}; Do[If[DivisorSigma[0, n] == 25, Print[n]; AppendTo[lst, n]], {n, 55000000}]; lst (* _Vladimir Joseph Stephan Orlovsky_, May 03 2011 *)
%t Select[Range[5221*10^4],DivisorSigma[0,#]==25&] (* _Harvey P. Dale_, Mar 11 2019 *)
%o (Haskell)
%o a137488 n = a137488_list !! (n-1)
%o a137488_list = m (map (^ 24) a000040_list) (map (^ 4) a006881_list) where
%o m xs'@(x:xs) ys'@(y:ys) | x < y = x : m xs ys'
%o | otherwise = y : m xs' ys
%o -- _Reinhard Zumkeller_, Nov 29 2011
%o (PARI) is(n)=numdiv(n)==25 \\ _Charles R Greathouse IV_, Jun 19 2016
%o (Python)
%o from math import isqrt
%o from sympy import primepi, integer_nthroot, primerange
%o def A137488(n):
%o def bisection(f,kmin=0,kmax=1):
%o while f(kmax) > kmax: kmax <<= 1
%o kmin = kmax >> 1
%o while kmax-kmin > 1:
%o kmid = kmax+kmin>>1
%o if f(kmid) <= kmid:
%o kmax = kmid
%o else:
%o kmin = kmid
%o return kmax
%o def f(x): return int(n+x+(t:=primepi(s:=isqrt(y:=integer_nthroot(x,4)[0])))+(t*(t-1)>>1)-sum(primepi(y//k) for k in primerange(1, s+1)))-primepi(integer_nthroot(x,24)[0])
%o return bisection(f,n,n) # _Chai Wah Wu_, Feb 22 2025
%Y Cf. A000005, A010812, A030513-A030516, A030626, A030627, A030634-A030638, A005179, A003680, A096932, A061286, A061283, A135581, A175755.
%K nonn
%O 1,1
%A _R. J. Mathar_, Apr 22 2008