A189991 Numbers with prime factorization p^4*q^4.
1296, 10000, 38416, 50625, 194481, 234256, 456976, 1185921, 1336336, 1500625, 2085136, 2313441, 4477456, 6765201, 9150625, 10556001, 11316496, 14776336, 17850625, 22667121, 29986576, 35153041, 45212176, 52200625, 54700816, 57289761, 68574961, 74805201
Offset: 1
Keywords
Links
- T. D. Noe, Table of n, a(n) for n = 1..1000
- Will Nicholes, List of Prime Signatures
- Index to sequences related to prime signature
Programs
-
Mathematica
lst = {}; Do[If[Sort[Last/@FactorInteger[n]] == {4, 4}, Print[n]; AppendTo[lst, n]], {n,55000000}]; lst (* Orlovsky *) 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 *) With[{nn=30},Take[Union[Times@@@(Subsets[Prime[Range[nn]],{2}]^4)],nn]] (* Harvey P. Dale, Mar 05 2015 *)
-
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
-
Python
from math import isqrt from sympy import primepi, integer_nthroot, primerange def A189991(n): def bisection(f,kmin=0,kmax=1): while f(kmax) > kmax: kmax <<= 1 kmin = kmax >> 1 while kmax-kmin > 1: kmid = kmax+kmin>>1 if f(kmid) <= kmid: kmax = kmid else: kmin = kmid return kmax 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))) return bisection(f,n,n) # Chai Wah Wu, Feb 22 2025
Formula
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
Comments