9

1

Is there a way to generate a counting function for prime powers - i.e. to create a similar function to `PrimePi`

, but including prime powers. The following will, of course, generate a list of these numbers:

```
Select[ Range[1000000], PrimePowerQ]
```

9

1

Is there a way to generate a counting function for prime powers - i.e. to create a similar function to `PrimePi`

, but including prime powers. The following will, of course, generate a list of these numbers:

```
Select[ Range[1000000], PrimePowerQ]
```

11

A naive approach would be this:

```
primePower[n_] := Count[ Range @ n, _?PrimePowerQ]
```

This function works well however it might be very inefficient for large `n`

. It takes a bit to evaluate e.g.

```
primePower[10^6]
```

`78734`

which is only a little bigger than

```
PrimePi[10^6]
```

`78498`

The latter is much more efficient since it uses advanced algorithms for counting primes which exploit sparse caching and sieving techniques ( in documentation pages see Some Notes on Internal Implementation, then J. C. Lagarias, V. S. Miller and A. M. Odlyzko "Computing `π(x)`

: the Meissel-Lehmer method," Math. Comp., 44 (1985) 537-560, another interesting resources The Prime Pages: Prime Number Research, Records and Results).

Therefore we proceed along a different way making use of `PrimePi`

to count prime powers, namely we need to sum `Log[ Prime[k], n]`

(it measure how many times appears given prime with different powers up to `n`

, e.g. for `k == 3`

there are `5, 25, 125, 625, ..., 5^Floor[ Log[ Prime[3], n]]`

) over every prime up to `PrimePi[ Sqrt[n] // Floor]`

:

```
primePowerPi[n_] := PrimePi[n] + Sum[ Floor @ Log[ Prime[k], n] - 1,
{ k, PrimePi[ Sqrt[n] // Floor]}]
```

This might be done as well with:

```
Floor @ Log[ Prime[ Range[ PrimePi[ Sqrt[n] // Floor]]], n]
```

e.g.

```
Total[ Floor @ Log[ Prime[ Range[ PrimePi[ Sqrt[10^5] // Floor]]], 10^5] - 1] ==
primePowerPi[10^5] - PrimePi[10^5]
primePowerPi[10^5] - PrimePi[10^5]
```

`True 108`

Now we can find immediately:

```
primePowerPi /@ { 10^6, 10^7}
```

`{78734, 665134}`

while

```
PrimePi[10^7]
```

`664579`

For some limitations of `Prime`

and `PrimePi`

domains this reference would be helpful What is so special about Prime?.

1an ingenious approach! Thanks for your input on this - it has greatly helped my understanding of generating functions :) – martin – 2013-11-20T08:04:23.720

1this was a really in-depth answer that helped me to understand how *Mathematica* works - again, many thanks :) – martin – 2013-11-20T08:12:03.950

6

Here's a fancy memoized solution:

```
Clear[primePowerCount, primePowerCount`cache, iPrimePowerCount]
iPrimePowerCount[n1_, n2_] := Count[Range[n1, n2], _?PrimePowerQ]
primePowerCount`cache = {1};
primePowerCount[1] = 0;
primePowerCount[n_?Positive] :=
Module[{n0, res},
n0 = First@Nearest[primePowerCount`cache, n];
If[n0 < n,
res = primePowerCount[n0] + iPrimePowerCount[n0 + 1, n],
res = primePowerCount[n0] - iPrimePowerCount[n + 1, n0]
];
primePowerCount`cache = Union[primePowerCount`cache, {n}];
primePowerCount[n] = res
]
```

This is going to be conveniently fast for interactive use with "random" inputs. If you need to evaluate it for a complete `Range`

though, it will not be the best way.

The idea is to store all previously calculated counts, and to compute any newly requested count relative to an already known (cached) one.

You can even pre-initialize this for a range of exponentially increasing values, e.g. `Scan[primePowerCount, Union@Round[1.1^Range[150]]]`

, to speed things up for subsequent computations. These results could be saved and later re-loaded.

EDIT: I just saw Artes's answer, which is the best solution, as it exploits the specific problem. For general slow-to-evaluate brute force counting functions, I believe that my memoized approach is still valuable.

1thank you for your reply - not sure I understand it yet - but I will persevere with it! :) – martin – 2013-11-20T08:06:37.623

5

I can't complete with Artes's mathematical knowledge and approach, but simply as a point of reference, for formulating a brute-force approach it will be more memory efficient to use `Sum`

, though it will still be very slow for large input.

```
Sum[Boole @ PrimePowerQ @ i, {i, 5*^6}] // Timing
MaxMemoryUsed[]
```

`{46.535, 348940} 15083688`

1@MrWizard - thank you for your solution - it is great from a learning point of view :) – martin – 2013-11-20T08:10:13.237

5

Late again, but ...

In addition to`PrimePowerQ`

, the built-in`MangoldtLambda[n]`

also works as in

```
n - Count[MangoldtLambda[Range[n]], 0]
```

where n is the upper limit. However, these functions are both too slow.
The answer by @Artes uses`PrimePi`

effectively for *much* faster results.

A cool method given by Riemann himself is to use`PrimePi`

on fractional powers
of the upper limit n.

```
Sum[PrimePi[n^(1/k)], {k, Floor[Log[2, n]]}]
```

On my machine with `n=10^12`

, the`AbsoluteTiming`

was just 0.015s, whereas`primePowerPi`

required several seconds.

Something like this maybe?

`Module[{count = 0}, Do[If[PrimePowerQ[k], count++], {k, 1, 100000}]; count]`

Do I understand that you want to avoid keeping any long lists in memory? In that case you need to use`Do`

. Though in practice your limiting factor may be processing power and not memory, so`Count[Range[100000], _?PrimePowerQ]`

will do. – Szabolcs – 2013-11-20T02:28:00.790