3
$\begingroup$

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?

$\endgroup$
4
  • $\begingroup$ PrimePi[10^12] is a very large number ... $\endgroup$ Commented Nov 26, 2014 at 17:17
  • $\begingroup$ @belisarius yes, I was rather hoping I could go that high ... $\endgroup$ Commented Nov 26, 2014 at 17:19
  • $\begingroup$ 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. $\endgroup$ Commented Nov 26, 2014 at 18:06
  • $\begingroup$ @Artes I'm not sure I follow - how would this help me find the nth prime power? $\endgroup$ Commented Nov 26, 2014 at 18:13

2 Answers 2

4
$\begingroup$

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]] 
$\endgroup$
1
  • 3
    $\begingroup$ 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]] $\endgroup$ Commented Nov 26, 2014 at 21:28
3
$\begingroup$

How about

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 
$\endgroup$
3
  • $\begingroup$ is there any way of generating the nth prime power? $\endgroup$ Commented Nov 26, 2014 at 17:35
  • $\begingroup$ 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. $\endgroup$ Commented Nov 26, 2014 at 18:04
  • $\begingroup$ @DanielLichtblau I don't suppose you could elaborate in an answer, could you? $\endgroup$ Commented Nov 26, 2014 at 18:05

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.