# 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, 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 ... Nov 26, 2014 at 17:17
• @belisarius yes, I was rather hoping I could go that high ... Nov 26, 2014 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. Nov 26, 2014 at 18:06
• @Artes I'm not sure I follow - how would this help me find the nth prime power? Nov 26, 2014 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

9999991

PrimePowerPi[10^12]

37607992088

PrimePower

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 :) Nov 26, 2014 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]] Nov 26, 2014 at 21:28
• @DanielLichtblau just seen thas - thank you! :) Nov 28, 2014 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? Nov 26, 2014 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. Nov 26, 2014 at 18:04
• @DanielLichtblau I don't suppose you could elaborate in an answer, could you? Nov 26, 2014 at 18:05