A376092 - OEIS

%I #22 Sep 10 2024 20:19:03

%S 1,49,3136,253472,23002083,2200079025,215523459072,21348015504200,

%T 2125390162618116,212104218976916644,21190268970925690248,

%U 2118092209873957381248,211765852717674823741924,21174572668805230623003225,2117363857447354911021280900

%N 10^n-th powerful number.

%H Chai Wah Wu, <a href="/A376092/b376092.txt">Table of n, a(n) for n = 0..16</a>

%F a(n) = A001694(10^n).

%F Limit_{n->oo} a(n)/10^(2n) = (zeta(3)/zeta(3/2))^2 = 0.21172829478335...

%o (Python)

%o from math import isqrt

%o from sympy import mobius, integer_nthroot

%o def A376092(n):

%o     def squarefreepi(n):

%o         return int(sum(mobius(k)*(n//k**2) for k in range(1, isqrt(n)+1)))

%o     def bisection(f, kmin=0, kmax=1):

%o         while f(kmax) > kmax: 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     m = 10**n

%o     def f(x):

%o         c, l = m+x, 0

%o         j = isqrt(x)

%o         while j>1:

%o             k2 = integer_nthroot(x//j**2,3)[0]+1

%o             w = squarefreepi(k2-1)

%o             c -= j*(w-l)

%o             l, j = w, isqrt(x//k2**3)

%o         c -= squarefreepi(integer_nthroot(x,3)[0])-l

%o         return c

%o     return bisection(f,m,m)

%Y Cf. A001694, A118896.

%K nonn,more

%O 0,2

%A _Chai Wah Wu_, Sep 09 2024