# Update

Is there any better way of generating the nth prime power?

Chip Hurst gave a great solution for a list below, so

PrimePowersUpTo[x_] := Union @@ Table[Prime[Range[PrimePi[x^(1/n)]]]^n, {n, Log2[x]}]
pp=PrimePowersUpTo[10^7];
pp[[n]]


would work, but I was wondering whether there was a more direct way?

# Original question

I am trying to get a list of prime powers, but I can only seem to get one that goesx to 10^7.

xmax = 10000000;
rr = DeleteCases[Table[If[PrimePowerQ[n] == True, n, 0], {n, 1, xmax}], 0];


or

Select[Range[10], PrimePowerQ]


I tried using Artes' idea of taking log base prime as per here, but didn't get very far.

Is there any way I can go to 10^12?

• PrimePi[10^12] is a very large number ... – Dr. belisarius Nov 26 '14 at 17:17
• @belisarius yes, I was rather hoping I could go that high ... – martin Nov 26 '14 at 17:19
• You can estimate reasonably PrimePi with RiemannR even if Mathematica implementation doesn't help, e.g. RiemannR[10^12] // N yields 3.76079*10^10. See also Approximation to the prime counting function. – Artes Nov 26 '14 at 18:06
• @Artes I'm not sure I follow - how would this help me find the nth prime power? – martin Nov 26 '14 at 18:13

Here's my way of finding the nth prime power. I do a binary search. Unfortunately it's slow for large inputs.

lower[_] = 2;
upper[x_] := Prime[x]

PrimePower[n_] := Module[{lo = lower[n], hi = upper[n], mid, pp},
mid = Quotient[lo + hi, 2];

While[lo < mid < hi && (pp = PrimePowerPi[mid]) != n,
If[pp > n,
hi = mid,
lo = mid
];
mid = Quotient[lo + hi, 2];
];

PreviousPrimePower[mid]
]

PreviousPrimePower[x_] := Max[Table[NextPrime[x^(1/n), -1]^n, {n, Log2[x]}]]


Examples:

PrimePowerPi[10^7]

665134

PrimePower[665134]

9999991

PrimePowerPi[10^12]

37607992088

PrimePower[37607992088]

999999999989


From Daniel Lichtblau:

Feel free to use this variant (divides more finely).

PrimePowerb[n_] := Module[
{lo = lower[n], hi = upper[n], mid, ppmid, pphi, pplo},
mid = Round[.95*hi];
ppmid = PrimePowerPi[mid];
pplo = 1;
pphi = PrimePowerPi[hi];
While[lo < mid < hi && ppmid != n,
If[ppmid > n,
hi = mid;
pphi = ppmid;
mid = mid + Round[(lo - mid)*(n - ppmid)/(pplo - ppmid)];
,
lo = mid;
pplo = ppmid;
mid = mid + Round[(hi - mid)*(n - ppmid)/(pphi - ppmid)];
];
If[mid == lo, mid++];
If[mid == hi, mid--];
ppmid = PrimePowerPi[mid];
];
PreviousPrimePower[mid]]

• thank you! this works really well :) – martin Nov 26 '14 at 18:45
• Feel free to use this variant (divides more finely). PrimePowerb[n_] := Module[ {lo = lower[n], hi = upper[n], mid, ppmid, pphi, pplo}, mid = Round[.95*hi]; ppmid = PrimePowerPi[mid]; pplo = 1; pphi = PrimePowerPi[hi]; While[lo < mid < hi && ppmid != n, If[ppmid > n, hi = mid; pphi = ppmid; mid = mid + Round[(lo - mid)*(n - ppmid)/(pplo - ppmid)]; , lo = mid; pplo = ppmid; mid = mid + Round[(hi - mid)*(n - ppmid)/(pphi - ppmid)]; ]; If[mid == lo, mid++]; If[mid == hi, mid--]; ppmid = PrimePowerPi[mid]; ]; PreviousPrimePower[mid]] – Daniel Lichtblau Nov 26 '14 at 21:28
• @DanielLichtblau just seen thas - thank you! :) – martin Nov 28 '14 at 18:03

PrimePowersUpTo[x_] := Union @@ Table[Prime[Range[PrimePi[x^(1/n)]]]^n, {n, Log2[x]}]

PrimePowersUpTo[10^7] // Length

665134


Now about getting the first 10^12 prime powers. You would want to get them in batches (segmented sieve) because there are about 38 billion of them:

PrimePowerPi[x_] := Sum[PrimePi[x^(1/n)], {n, Log2[x]}]

PrimePowerPi[10^12]

37607992088

• is there any way of generating the nth prime power? – martin Nov 26 '14 at 17:35
• You might generate the nth prime and then subtract your way down to account for the relatively fewer powers that also near that value. Possibly you can use a tighter estimate since you have ways to approximate how many of those there are, but then you also have to count carefully so as not to include powers that are too large. – Daniel Lichtblau Nov 26 '14 at 18:04
• @DanielLichtblau I don't suppose you could elaborate in an answer, could you? – martin Nov 26 '14 at 18:05