A336178 - OEIS

%I #12 Sep 15 2024 22:01:48

%S 31,36,67,93,132,140,145,161,166,189,192,220,223,265,280,290,296,311,

%T 316,322,364,384,407,468,537,576,592,602,623,639,644,656,659,661,670,

%U 690,722,769,771,793,828,883,888,890,896,950,961,981,984,987,992,995,1018

%N Numbers k such that there are exactly three powerful numbers between k^2 and (k+1)^2.

%C Positions of 3's in A119241.

%C Shiu (1980) proved that this sequence has an asymptotic density = 0.0770... A more accurate calculation using his formula gives 0.0770742722233...

%D József Sándor, Dragoslav S. Mitrinovic and Borislav Crstici, Handbook of Number Theory I, Springer Science & Business Media, 2005, chapter VI, p. 226.

%H Amiram Eldar, <a href="/A336178/b336178.txt">Table of n, a(n) for n = 1..10000</a>

%H P. Shiu, <a href="https://doi.org/10.1112/S0025579300010056">On the number of square-full integers between successive squares</a>, Mathematika, Vol. 27, No. 2 (1980), pp. 171-178.

%e 31 is a term since there are exactly three powerful numbers, 968 = 2^3 * 11^2, 972 = 2^2 * 3^5 and 1000 = 2^3 * 5^3 between 31^2 = 961 and (31+1)^2 = 1024.

%t powQ[n_] := (n == 1) || Min @@ FactorInteger[n][[;; , 2]] > 1; Select[Range[1000], Count[Range[#^2 + 1, (# + 1)^2 - 1], _?powQ] == 3 &]

%o (Python)

%o from functools import lru_cache

%o from math import isqrt

%o from sympy import mobius, integer_nthroot

%o def A336178(n):

%o     def squarefreepi(n): 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     @lru_cache(maxsize=None)

%o     def g(x):

%o         c, l = 0, 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     def f(x):

%o         c, a = n+x, 1

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

%o             b = g((k+1)**2)

%o             if b == a+4:

%o                 c -= 1

%o             a = b

%o         return c

%o     return bisection(f,n,n) # _Chai Wah Wu_, Sep 14 2024

%Y Cf. A001694, A119241, A119242, A336175, A336176, A336177.

%K nonn

%O 1,1

%A _Amiram Eldar_, Jul 10 2020