# Greetings from The On-Line Encyclopedia of Integer Sequences! http://oeis.org/
Search: id:a371190
Showing 1-1 of 1
%I A371190 #16 Sep 15 2024 22:02:12
%S A371190 1,4,8,25,32,288,675,968,1152,1369,2700,9800,12167,39200,48668,70225,
%T A371190 235224,332928,465124,1331712,1825200,5724500,7300800,11309768,
%U A371190 78960996,189750625,263672644,384199200,592192224,912670088,1536796800,2368768896,4931691075,5425069447,8957108164
%N A371190 The smaller of a pair of successive powerful numbers without a nonsquarefree number between them.
%H A371190 Amiram Eldar, Table of n, a(n) for n = 1..64
%H A371190 Index entries for sequences related to powerful numbers.
%e A371190 1 is a term since 1 and 4 are successive powerful numbers and the numbers between them, 2 and 3, are both squarefree.
%t A371190 seq[max_] := Module[{pows = Union[Flatten[Table[i^2*j^3, {j, 1, Surd[max, 3]}, {i, 1, Sqrt[max/j^3]}]]], s = {}}, Do[If[AllTrue[Range[pows[[k]] + 1, pows[[k + 1]] - 1], SquareFreeQ], AppendTo[s, pows[[k]]]], {k, 1, Length[pows] - 1}]; s]; seq[10^10]
%o A371190 (PARI) lista(mx) = {my(s = List(), is); for(j = 1, sqrtnint(mx, 3), for(i = 1, sqrtint(mx\j^3), listput(s, i^2 * j^3))); s = Set(s); for(i = 1, #s - 1, is = 1; for(k = s[i]+1, s[i+1]-1, if(!issquarefree(k), is = 0; break)); if(is, print1(s[i], ", ")));}
%o A371190 (Python)
%o A371190 from math import isqrt
%o A371190 from sympy import mobius, integer_nthroot
%o A371190 def A371190_gen(): # generator of terms
%o A371190 def squarefreepi(n): return int(sum(mobius(k)*(n//k**2) for k in range(1, isqrt(n)+1)))
%o A371190 def bisection(f,kmin=0,kmax=1):
%o A371190 while f(kmax) > kmax: kmax <<= 1
%o A371190 while kmax-kmin > 1:
%o A371190 kmid = kmax+kmin>>1
%o A371190 if f(kmid) <= kmid:
%o A371190 kmax = kmid
%o A371190 else:
%o A371190 kmin = kmid
%o A371190 return kmax
%o A371190 def f(x):
%o A371190 c, l, j = x-squarefreepi(integer_nthroot(x,3)[0]), 0, isqrt(x)
%o A371190 while j>1:
%o A371190 k2 = integer_nthroot(x//j**2,3)[0]+1
%o A371190 w = squarefreepi(k2-1)
%o A371190 c -= j*(w-l)
%o A371190 l, j = w, isqrt(x//k2**3)
%o A371190 return c+l
%o A371190 m, w = 1, 1
%o A371190 for n in count(2):
%o A371190 k = bisection(lambda x:f(x)+n,m,m)
%o A371190 if (a:=squarefreepi(k))-w==k-1-m:
%o A371190 yield m
%o A371190 m, w = k, a # _Chai Wah Wu_, Sep 15 2024
%Y A371190 Cf. A001694, A005117, A013929, A240591.
%K A371190 nonn
%O A371190 1,2
%A A371190 _Amiram Eldar_, Mar 14 2024
# Content is available under The OEIS End-User License Agreement: http://oeis.org/LICENSE