cp's OEIS Frontend

A113877 Semiprimes to semiprime powers.

%I A113877 #27 Sep 12 2024 12:42:55
%S A113877 256,1296,4096,6561,10000,38416,46656,50625,194481,234256,262144,
%T A113877 390625,456976,531441,1000000,1048576,1185921,1336336,1500625,2085136,
%U A113877 2313441,4477456,5764801,6765201,7529536,9150625,10077696,10556001,11316496,11390625,14776336
%N A113877 Semiprimes to semiprime powers.
%C A113877 This is the semiprime analog of A053810.
%H A113877 Charles R Greathouse IV, <a href="/A113877/b113877.txt">Table of n, a(n) for n = 1..10000</a>
%F A113877 {a(n)} = {a^b where a and b are elements of A001358}.
%F A113877 {a(n)} = {(p*q)^(r*s) = (p^(r*s))*(q^r*s) for distinct primes p, q, r, s} UNION {(p*q)^(p*r) = (p^(p*r))*(q^(p*r)) for distinct primes p, q, r} UNION {(p*q)^(r*r) = (p^(r^2))*(q^(r^2)) for distinct primes p, q, r} UNION {(p*q)^(p*q)= (p^(p*q))*(q^(p*q)) for distinct primes p, q} UNION {(p^2)^(p^2) = p^(2*(p^2)) for prime p}.
%F A113877 a(n) ~ (n log n/log log n)^4. - _Charles R Greathouse IV_, Jun 05 2013
%e A113877 a(1) = 256 = 4^4 = semiprime(1)^semiprime(1).
%e A113877 a(2) = 1296 = 6^4 = semiprime(2)^semiprime(1).
%e A113877 a(3) = 4096 = 4^6 = semiprime(1)^semiprime(2).
%e A113877 a(4) = 6561 = 9^4 = semiprime(3)^semiprime(1).
%e A113877 a(5) = 10000 = 10^4.
%t A113877 lim = 10^8; s = Select[Range[lim^(1/4)], Total[Transpose[FactorInteger[#]][[2]]] == 2 &]; t = {}; j = 1; While[b = s[[j]]; i = 1; While[a = s[[i]]; e = a^b; If[e <= lim, AppendTo[t, e]]; e < lim && i < Length[s], i++]; i > 1, j++]; t = Union[t] (* _T. D. Noe_, Jun 05 2013 *)
%o A113877 (PARI) is(n)=my(b,e=ispower(n,,&b),o); if(e==0,return(0)); o=bigomega(e); (o==2 && bigomega(b)==2) || (e%2==0 && o==3 && isprime(b)) \\ _Charles R Greathouse IV_, Jun 05 2013
%o A113877 (PARI) list(lim)=my(v=List());for(e=4,log(lim\=1+.5)\log(4), if(bigomega(e)!=2, next); for(b=4,(lim+.5)^(1/e), if(bigomega(b)==2, listput(v,b^e)))); vecsort(Vec(v)) \\ _Charles R Greathouse IV_, Jun 05 2013
%o A113877 (Python)
%o A113877 from math import isqrt
%o A113877 from sympy import primepi, primerange, integer_nthroot, factorint
%o A113877 def A113877(n):
%o A113877     def A072000(n): return int(-((t:=primepi(s:=isqrt(n)))*(t-1)>>1)+sum(primepi(n//p) for p in primerange(s+1)))
%o A113877     def f(x): return int(n+x-sum(A072000(integer_nthroot(x, p)[0]) for p in range(4,x.bit_length()) if sum(factorint(p).values())==2))
%o A113877     def bisection(f,kmin=0,kmax=1):
%o A113877         while f(kmax) > kmax: kmax <<= 1
%o A113877         while kmax-kmin > 1:
%o A113877             kmid = kmax+kmin>>1
%o A113877             if f(kmid) <= kmid:
%o A113877                 kmax = kmid
%o A113877             else:
%o A113877                 kmin = kmid
%o A113877         return kmax
%o A113877     return bisection(f,n,n) # _Chai Wah Wu_, Sep 12 2024
%Y A113877 Cf. A001358, A113854, A053810.
%K A113877 easy,nonn
%O A113877 1,1
%A A113877 _Jonathan Vos Post_, Jan 27 2006
%E A113877 Terms corrected by _Charles R Greathouse IV_, Jun 05 2013