A358173 - OEIS

%I #14 Sep 11 2024 00:34:35

%S 36,28,8,36,52,4,16,9,63,36,68,8,32,9,43,16,76,72,27,1,108,16,64,36,

%T 68,4,28,89,36,27,4,69,71,27,29,20,72,77,47,32,128,36,36,136,8,56,25,

%U 91,188,8,188,92,9,99,4,40,144,28,109,62,49,64,49,18,97,11,81

%N First differences of A286708.

%C Consider the sequence of powerful numbers A001694, superset of A246547, the sequence of composite prime powers. Let s = A001694(k) such that omega(s) > 1 be followed by t = A001694(k+1) such that omega(t) = 1.

%C Since A286708 = A001694 \ A246547, prime powers t are missing in A286708. We consider s = A286708(j) and note that the difference A286708(j+1) - A286708(j) > A001694(k+1) - A001694(k).

%C Therefore we see a subset S containing s in A286708 that plots "out of place" with respect to the complementary subset R = A286708 \ S; some of this subset S exceeds the maxima of R in the scatterplot of this sequence. The plot of the R resembles the scatterplot of A001694.

%H Michael De Vlieger, <a href="/A358173/b358173.txt">Table of n, a(n) for n = 1..10000</a>

%H Michael De Vlieger, <a href="/A358173/a358173.png">Scatterplot of a(n)</a>, n = 1..2^16, highlighting in red terms s that in A001694 are succeeded by a prime power t.

%e The number 36 is the smallest powerful number that is not a prime power; the next powerful number that is not a prime power is 72, and their difference is 36, hence a(1) = 36.

%t With[{nn = 2^25}, Differences@ Select[Rest@ Union@ Flatten@ Table[a^2*b^3, {b, nn^(1/3)}, {a, Sqrt[nn/b^3]}], ! PrimePowerQ[#] &]]

%o (Python)

%o from math import isqrt

%o from sympy import integer_nthroot, primepi, mobius

%o def A358173(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     def f(x):

%o         c, l = n+x+1+sum(primepi(integer_nthroot(x, k)[0]) for k in range(2, x.bit_length())), 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 -(a:=bisection(f,n,n))+bisection(lambda x:f(x)+1,a,a) # _Chai Wah Wu_, Sep 10 2024

%Y Cf. A001694, A053707, A076446, A246547, A286708.

%K nonn

%O 1,1

%A _Michael De Vlieger_, Nov 01 2022