cp's OEIS Frontend

A189991 Numbers with prime factorization p^4*q^4.

%I A189991 #36 Mar 28 2025 03:27:34
%S A189991 1296,10000,38416,50625,194481,234256,456976,1185921,1336336,1500625,
%T A189991 2085136,2313441,4477456,6765201,9150625,10556001,11316496,14776336,
%U A189991 17850625,22667121,29986576,35153041,45212176,52200625,54700816,57289761,68574961,74805201
%N A189991 Numbers with prime factorization p^4*q^4.
%C A189991 The primes p and q must be distinct, or else the product has factorization p^8 (or q^8, for that matter).
%H A189991 T. D. Noe, <a href="/A189991/b189991.txt">Table of n, a(n) for n = 1..1000</a>
%H A189991 Will Nicholes, <a href="https://willnicholes.com/2010/06/06/list-of-prime-signatures">List of Prime Signatures</a>
%H A189991 <a href="/index/Pri#prime_signature">Index to sequences related to prime signature</a>
%F A189991 Sum_{n>=1} 1/a(n) = (P(4)^2 - P(8))/2 = (A085964^2 - A085968)/2 = 0.000933..., where P is the prime zeta function. - _Amiram Eldar_, Jul 06 2020
%F A189991 a(n) = A006881(n)^4 = A085986(n)^2. - _Chai Wah Wu_, Mar 27 2025
%t A189991 lst = {}; Do[If[Sort[Last/@FactorInteger[n]] == {4, 4}, Print[n]; AppendTo[lst, n]], {n,55000000}]; lst (* Orlovsky *)
%t A189991 lim = 10^8; pMax = PrimePi[(lim/16)^(1/4)]; Select[Union[Flatten[Table[Prime[i]^4 Prime[j]^4, {i, 2, pMax}, {j, i - 1}]]], # <= lim &] (* _Alonso del Arte_, May 03 2011 *)
%t A189991 With[{nn=30},Take[Union[Times@@@(Subsets[Prime[Range[nn]],{2}]^4)],nn]] (* _Harvey P. Dale_, Mar 05 2015 *)
%o A189991 (PARI) list(lim)=my(v=List(),t);forprime(p=2, lim^(1/8), t=p^4;forprime(q=p+1, (lim\t)^(1/4), listput(v,t*q^4))); vecsort(Vec(v)) \\ _Charles R Greathouse IV_, Jul 24 2011
%o A189991 (Python)
%o A189991 from math import isqrt
%o A189991 from sympy import primepi, integer_nthroot, primerange
%o A189991 def A189991(n):
%o A189991     def bisection(f,kmin=0,kmax=1):
%o A189991         while f(kmax) > kmax: kmax <<= 1
%o A189991         kmin = kmax >> 1
%o A189991         while kmax-kmin > 1:
%o A189991             kmid = kmax+kmin>>1
%o A189991             if f(kmid) <= kmid:
%o A189991                 kmax = kmid
%o A189991             else:
%o A189991                 kmin = kmid
%o A189991         return kmax
%o A189991     def f(x): return int(n+x+(t:=primepi(s:=isqrt(y:=integer_nthroot(x,4)[0])))+(t*(t-1)>>1)-sum(primepi(y//k) for k in primerange(1, s+1)))
%o A189991     return bisection(f,n,n) # _Chai Wah Wu_, Feb 22 2025
%Y A189991 Cf. A006881, A137488, A179671, A189990.
%Y A189991 Cf. A085964, A085968, A085986.
%K A189991 nonn
%O A189991 1,1
%A A189991 _Vladimir Joseph Stephan Orlovsky_, May 03 2011