A216426 - OEIS

%I #27 Sep 14 2024 12:31:27

%S 72,108,128,200,256,288,392,432,500,512,576,648,675,800,864,968,972,

%T 1125,1152,1323,1352,1372,1568,1600,1728,1800,1944,2000,2048,2187,

%U 2304,2312,2592,2700,2888,2916,3087,3136,3200,3267,3456,3528,3872,3888,4000

%N Numbers of the form a^2*b^3, where a != b and a, b > 1.

%C Terms of A216427 that are not 5th powers of squarefree numbers (A113850) and not 10th powers of primes (A030629). - _Amiram Eldar_, Feb 07 2023

%H Charles R Greathouse IV, <a href="/A216426/b216426.txt">Table of n, a(n) for n = 1..10000</a>

%F Sum_{n>=1} 1/a(n) = 2 + ((zeta(2)-1)*(zeta(3)-1)-1)/zeta(6) - zeta(5)/zeta(10) - P(6) - P(10) = 0.09117811499514578262..., where P(s) is the prime zeta function. - _Amiram Eldar_, Feb 07 2023

%t With[{upto=4000},Select[Union[Flatten[{#[[1]]^2 #[[2]]^3,#[[2]]^2 #[[1]]^3}& /@ Subsets[Range[2,Surd[upto,2]],{2}]]],#<=upto&]](* _Harvey P. Dale_, Jan 04 2014 *)

%t pMx = 25; mx = 2^3 pMx^2; t = Flatten[Table[If[a != b, a^2 b^3, 0], {a, 2, mx^(1/2)}, {b, 2, mx^(1/3)}]]; Union[Select[t, 0 < # <= mx &]] (* _T. D. Noe_, Jan 02 2014 *)

%o (PARI) list(lim)=my(v=List()); for(b=2, sqrtnint(lim\4,3), for(a=2, sqrtint(lim\b^3), if(a!=b, listput(v, a^2*b^3)))); Set(v) \\ _Charles R Greathouse IV_, Jan 02 2014

%o (Python)

%o from math import isqrt

%o from sympy import integer_nthroot, mobius, primepi

%o def A216426(n):

%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         j, b, a, d = isqrt(x), integer_nthroot(x,6)[0], integer_nthroot(x,5)[0], integer_nthroot(x,10)[0]

%o         l, c = 0, n+x-2+primepi(b)+sum(mobius(k)*(j//k**3) for k in range(d+1, b+1))+primepi(d)+sum(mobius(k)*(a//k**2+j//k**3) for k in range(1, d+1))

%o         while j>1:

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

%o             w = sum(mobius(k)*((k2-1)//k**2) for k in range(1, isqrt(k2-1)+1))

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

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

%o         return c+l

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

%Y Cf. A143610.

%Y Subsequence of A216427. - _Zak Seidov_, Jan 03 2014

%Y Cf. A002117, A013661, A013663, A013664, A013668, A030629, A085966, A113850.

%K nonn

%O 1,1

%A _V. Raman_, Sep 07 2012

%E Name corrected by _Charles R Greathouse IV_, Jan 02 2014