Legendre prime counting function: Difference between revisions
PicoLisp version |
Thundergnat (talk | contribs) m syntax highlighting fixup automation |
||
Line 23: | Line 23: | ||
=={{header|C++}}== |
=={{header|C++}}== |
||
< |
<syntaxhighlight lang="cpp">#include <cmath> |
||
#include <iostream> |
#include <iostream> |
||
#include <unordered_map> |
#include <unordered_map> |
||
Line 83: | Line 83: | ||
for (int i = 0, n = 1; i < 10; ++i, n *= 10) |
for (int i = 0, n = 1; i < 10; ++i, n *= 10) |
||
std::cout << "10^" << i << "\t" << counter.prime_count(n) << '\n'; |
std::cout << "10^" << i << "\t" << counter.prime_count(n) << '\n'; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 100: | Line 100: | ||
=={{header|CoffeeScript}}== |
=={{header|CoffeeScript}}== |
||
<syntaxhighlight lang="coffeescript"> |
|||
<lang CoffeeScript> |
|||
sorenson = require('sieve').primes # Sorenson's extensible sieve from task: Extensible Prime Generator |
sorenson = require('sieve').primes # Sorenson's extensible sieve from task: Extensible Prime Generator |
||
Line 135: | Line 135: | ||
console.log "10^#{i}\t#{pi(n)}" |
console.log "10^#{i}\t#{pi(n)}" |
||
n *= 10 |
n *= 10 |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{Out}} |
{{Out}} |
||
<pre> |
<pre> |
||
Line 150: | Line 150: | ||
</pre> |
</pre> |
||
=={{header|Erlang}}== |
=={{header|Erlang}}== |
||
<syntaxhighlight lang="erlang"> |
|||
<lang Erlang> |
|||
-mode(native). |
-mode(native). |
||
Line 203: | Line 203: | ||
true -> M + 1 |
true -> M + 1 |
||
end. |
end. |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{Out}} |
{{Out}} |
||
<pre> |
<pre> |
||
Line 225: | Line 225: | ||
The alternative of sieving a billion numbers and then counting them takes about 4.8 seconds using the present library function though this could be cut to 1.2 seconds using a third party library such as primegen. |
The alternative of sieving a billion numbers and then counting them takes about 4.8 seconds using the present library function though this could be cut to 1.2 seconds using a third party library such as primegen. |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 270: | Line 270: | ||
fmt.Printf("10^%d %d\n", i, pi(n)) |
fmt.Printf("10^%d %d\n", i, pi(n)) |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 287: | Line 287: | ||
=={{header|Java}}== |
=={{header|Java}}== |
||
< |
<syntaxhighlight lang="java">import java.util.*; |
||
public class LegendrePrimeCounter { |
public class LegendrePrimeCounter { |
||
Line 341: | Line 341: | ||
return primes; |
return primes; |
||
} |
} |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 360: | Line 360: | ||
Memoization utilities: |
Memoization utilities: |
||
< |
<syntaxhighlight lang="haskell">data Memo a = Node a (Memo a) (Memo a) |
||
deriving Functor |
deriving Functor |
||
Line 385: | Line 385: | ||
split [] = ([],[]) |
split [] = ([],[]) |
||
split [x] = ([x],[]) |
split [x] = ([x],[]) |
||
split (x:y:xs) = let (l,r) = split xs in (x:l, y:r)</ |
split (x:y:xs) = let (l,r) = split xs in (x:l, y:r)</syntaxhighlight> |
||
Computation of Legendre function: |
Computation of Legendre function: |
||
< |
<syntaxhighlight lang="haskell">isqrt :: Integer -> Integer |
||
isqrt n = go n 0 (q `shiftR` 2) |
isqrt n = go n 0 (q `shiftR` 2) |
||
where |
where |
||
Line 411: | Line 411: | ||
where a = legendrePi (floor (sqrt (fromInteger n))) |
where a = legendrePi (floor (sqrt (fromInteger n))) |
||
main = mapM_ (\n -> putStrLn $ show n ++ "\t" ++ show (legendrePi (10^n))) [1..7]</ |
main = mapM_ (\n -> putStrLn $ show n ++ "\t" ++ show (legendrePi (10^n))) [1..7]</syntaxhighlight> |
||
<pre>λ> main |
<pre>λ> main |
||
Line 427: | Line 427: | ||
This is "almost a primitive" in J: |
This is "almost a primitive" in J: |
||
< |
<syntaxhighlight lang="j"> require'format/printf' |
||
{{echo '10^%d: %d' sprintf y;p:inv 10^y}}"0 i.10 |
{{echo '10^%d: %d' sprintf y;p:inv 10^y}}"0 i.10 |
||
10^0: 0 |
10^0: 0 |
||
Line 438: | Line 438: | ||
10^7: 664579 |
10^7: 664579 |
||
10^8: 5761455 |
10^8: 5761455 |
||
10^9: 50847534</ |
10^9: 50847534</syntaxhighlight> |
||
A complete implementation of the legendre prime counting function would be <tt>(1&p: + p:inv)</tt>: |
A complete implementation of the legendre prime counting function would be <tt>(1&p: + p:inv)</tt>: |
||
< |
<syntaxhighlight lang="j"> (1&p: + p:inv) 10 |
||
4 |
4 |
||
(1&p: + p:inv) 11 |
(1&p: + p:inv) 11 |
||
5</ |
5</syntaxhighlight> |
||
But we can simplify that to p:inv when we know that primes are not being tested. |
But we can simplify that to p:inv when we know that primes are not being tested. |
||
Line 465: | Line 465: | ||
'''Preliminaries''' |
'''Preliminaries''' |
||
< |
<syntaxhighlight lang="jq"># For the sake of the accuracy of integer arithmetic when using gojq: |
||
def power($b): . as $in | reduce range(0;$b) as $i (1; . * $in); |
def power($b): . as $in | reduce range(0;$b) as $i (1; . * $in); |
||
Line 483: | Line 483: | ||
| reduce (2, 1 + (2 * range(1; $s))) as $i (.; erase($i)) |
| reduce (2, 1 + (2 * range(1; $s))) as $i (.; erase($i)) |
||
| map(select(.)) ; |
| map(select(.)) ; |
||
</syntaxhighlight> |
|||
</lang> |
|||
'''Phi and the Legendre Counting Function''' |
'''Phi and the Legendre Counting Function''' |
||
< |
<syntaxhighlight lang="jq">def legendre: |
||
(sqrt | floor + 1 | eratosthenes) as $primes |
(sqrt | floor + 1 | eratosthenes) as $primes |
||
Line 522: | Line 522: | ||
; |
; |
||
task</ |
task</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 540: | Line 540: | ||
=== Using the Julia Primes library === |
=== Using the Julia Primes library === |
||
This would be faster with a custom sieve that only counted instead of making a bitset, instead of the more versatile library function. |
This would be faster with a custom sieve that only counted instead of making a bitset, instead of the more versatile library function. |
||
< |
<syntaxhighlight lang="julia">using Primes |
||
function primepi(N) |
function primepi(N) |
||
Line 550: | Line 550: | ||
println("10^", rpad(power, 5), primepi(10^power)) |
println("10^", rpad(power, 5), primepi(10^power)) |
||
end |
end |
||
</ |
</syntaxhighlight>{{out}} |
||
<pre> |
<pre> |
||
10^0 0 |
10^0 0 |
||
Line 567: | Line 567: | ||
{{trans|CoffeeScript}} |
{{trans|CoffeeScript}} |
||
Run twice to show the benefits of precomputing the library prime functions. |
Run twice to show the benefits of precomputing the library prime functions. |
||
< |
<syntaxhighlight lang="julia">using Primes |
||
const maxpi = 1_000_000_000 |
const maxpi = 1_000_000_000 |
||
Line 594: | Line 594: | ||
end |
end |
||
</ |
</syntaxhighlight>{{out}} |
||
<pre> |
<pre> |
||
10^0 1 |
10^0 1 |
||
Line 621: | Line 621: | ||
=={{header|Mathematica}}/{{header|Wolfram Language}}== |
=={{header|Mathematica}}/{{header|Wolfram Language}}== |
||
< |
<syntaxhighlight lang="mathematica">ClearAll[Phi,pi] |
||
$RecursionLimit = 10^6; |
$RecursionLimit = 10^6; |
||
Phi[x_, 0] := x |
Phi[x_, 0] := x |
||
Phi[x_, a_] := Phi[x, a] = Phi[x, a - 1] - Phi[Floor[x/Prime[a]], a - 1] |
Phi[x_, a_] := Phi[x, a] = Phi[x, a - 1] - Phi[Floor[x/Prime[a]], a - 1] |
||
pi[n_] := Module[{a}, If[n < 2, 0, a = pi[Floor[Sqrt[n]]]; Phi[n, a] + a - 1]] |
pi[n_] := Module[{a}, If[n < 2, 0, a = pi[Floor[Sqrt[n]]]; Phi[n, a] + a - 1]] |
||
Scan[Print[pi[10^#]] &, Range[0,9]]</ |
Scan[Print[pi[10^#]] &, Range[0,9]]</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>0 |
<pre>0 |
||
Line 641: | Line 641: | ||
=={{header|Nim}}== |
=={{header|Nim}}== |
||
The program runs in 2.2 seconds on our laptop. Using int32 instead of naturals (64 bits on our 64 bits computer) saves 0.3 second. Using an integer rather than a tuple as key (for instance <code>x shl 32 or a</code>) saves 0.4 second. |
The program runs in 2.2 seconds on our laptop. Using int32 instead of naturals (64 bits on our 64 bits computer) saves 0.3 second. Using an integer rather than a tuple as key (for instance <code>x shl 32 or a</code>) saves 0.4 second. |
||
< |
<syntaxhighlight lang="nim">import math, strutils, sugar, tables |
||
const |
const |
||
Line 674: | Line 674: | ||
for i in 0..9: |
for i in 0..9: |
||
echo "π(10^$1) = $2".format(i, π(n)) |
echo "π(10^$1) = $2".format(i, π(n)) |
||
n *= 10</ |
n *= 10</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 689: | Line 689: | ||
=={{header|Perl}}== |
=={{header|Perl}}== |
||
< |
<syntaxhighlight lang="perl">#!/usr/bin/perl |
||
use strict; # https://rosettacode.org/wiki/Legendre_prime_counting_function |
use strict; # https://rosettacode.org/wiki/Legendre_prime_counting_function |
||
Line 717: | Line 717: | ||
{ |
{ |
||
printf "%d %12d %10d %10d\n", $_, 10**$_, pi(10**$_), prime_count(10**$_); |
printf "%d %12d %10d %10d\n", $_, 10**$_, pi(10**$_), prime_count(10**$_); |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 734: | Line 734: | ||
=={{header|Phix}}== |
=={{header|Phix}}== |
||
<!--< |
<!--<syntaxhighlight lang="phix">(phixonline)--> |
||
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
||
<span style="color: #000080;font-style:italic;">-- |
<span style="color: #000080;font-style:italic;">-- |
||
Line 782: | Line 782: | ||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
||
<span style="color: #0000FF;">?</span><span style="color: #7060A8;">elapsed</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">time</span><span style="color: #0000FF;">()-</span><span style="color: #000000;">t0</span><span style="color: #0000FF;">)</span> |
<span style="color: #0000FF;">?</span><span style="color: #7060A8;">elapsed</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">time</span><span style="color: #0000FF;">()-</span><span style="color: #000000;">t0</span><span style="color: #0000FF;">)</span> |
||
<!--</ |
<!--</syntaxhighlight>--> |
||
The commented-out length(get_primes_le()) gives the same results, in about the same time, |
The commented-out length(get_primes_le()) gives the same results, in about the same time, |
||
and in both cases re-running the main loop a second time finishes near-instantly. |
and in both cases re-running the main loop a second time finishes near-instantly. |
||
Line 803: | Line 803: | ||
=={{header|Picat}}== |
=={{header|Picat}}== |
||
Performance is surprisingly good, considering that the Picat's library SoE is very basic and that I recompute the base primes each time pi() is called. This version takes 23 seconds on my laptop, outperforming the CoffeeScript version which takes 30 seconds. |
Performance is surprisingly good, considering that the Picat's library SoE is very basic and that I recompute the base primes each time pi() is called. This version takes 23 seconds on my laptop, outperforming the CoffeeScript version which takes 30 seconds. |
||
<syntaxhighlight lang="picat"> |
|||
<lang Picat> |
|||
table(+, +, nt) |
table(+, +, nt) |
||
phi(X, 0, _) = X. |
phi(X, 0, _) = X. |
||
Line 819: | Line 819: | ||
writef("10^%w %w%n", K, pi(N)), |
writef("10^%w %w%n", K, pi(N)), |
||
N := N * 10. |
N := N * 10. |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{Out}} |
{{Out}} |
||
<pre> |
<pre> |
||
Line 836: | Line 836: | ||
=={{header|PicoLisp}}== |
=={{header|PicoLisp}}== |
||
Uses code from the Sieve of Eratosthenes to generate primes. |
Uses code from the Sieve of Eratosthenes to generate primes. |
||
<syntaxhighlight lang="picolisp"> |
|||
<lang PicoLisp> |
|||
(setq Legendre-Max (** 10 9)) |
(setq Legendre-Max (** 10 9)) |
||
Line 862: | Line 862: | ||
(bye) |
(bye) |
||
</syntaxhighlight> |
|||
</lang> |
|||
{{Out}} |
{{Out}} |
||
<pre> |
<pre> |
||
Line 879: | Line 879: | ||
=={{header|Python}}== |
=={{header|Python}}== |
||
Using <code>primesieve</code> module (<code>pip install primesieve</code>). |
Using <code>primesieve</code> module (<code>pip install primesieve</code>). |
||
< |
<syntaxhighlight lang="python">from primesieve import primes |
||
from math import isqrt |
from math import isqrt |
||
from functools import cache |
from functools import cache |
||
Line 902: | Line 902: | ||
for e in range(10): |
for e in range(10): |
||
print(f'10^{e}', legpi(10**e))</ |
print(f'10^{e}', legpi(10**e))</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>10^0 0 |
<pre>10^0 0 |
||
Line 917: | Line 917: | ||
=={{header|Raku}}== |
=={{header|Raku}}== |
||
Seems like an awful lot of pointless work. Using prime sieve anyway, why not just use it? |
Seems like an awful lot of pointless work. Using prime sieve anyway, why not just use it? |
||
<lang |
<syntaxhighlight lang="raku" line>use Math::Primesieve; |
||
my $sieve = Math::Primesieve.new; |
my $sieve = Math::Primesieve.new; |
||
Line 923: | Line 923: | ||
say "10^$_\t" ~ $sieve.count: exp($_,10) for ^10; |
say "10^$_\t" ~ $sieve.count: exp($_,10) for ^10; |
||
say (now - INIT now) ~ ' elapsed seconds';</ |
say (now - INIT now) ~ ' elapsed seconds';</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>10^0 0 |
<pre>10^0 0 |
||
Line 939: | Line 939: | ||
=={{header|Sidef}}== |
=={{header|Sidef}}== |
||
< |
<syntaxhighlight lang="ruby">func legendre_phi(x, a) is cached { |
||
return x if (a <= 0) |
return x if (a <= 0) |
||
__FUNC__(x, a-1) - __FUNC__(idiv(x, a.prime), a-1) |
__FUNC__(x, a-1) - __FUNC__(idiv(x, a.prime), a-1) |
||
Line 957: | Line 957: | ||
legendre_prime_count(10**n), prime_count(10**n)) |
legendre_prime_count(10**n), prime_count(10**n)) |
||
assert_eq(legendre_prime_count(10**n), prime_count(10**n)) |
assert_eq(legendre_prime_count(10**n), prime_count(10**n)) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 973: | Line 973: | ||
Alternative implementation of the Legendre phi function, by counting k-rough numbers <= n. |
Alternative implementation of the Legendre phi function, by counting k-rough numbers <= n. |
||
< |
<syntaxhighlight lang="ruby">func rough_count (n,k) { |
||
# Count of k-rough numbers <= n. |
# Count of k-rough numbers <= n. |
||
Line 1,028: | Line 1,028: | ||
legendre_prime_count(10**n), prime_count(10**n)) |
legendre_prime_count(10**n), prime_count(10**n)) |
||
assert_eq(legendre_prime_count(10**n), prime_count(10**n)) |
assert_eq(legendre_prime_count(10**n), prime_count(10**n)) |
||
}</ |
}</syntaxhighlight> |
||
=={{header|Swift}}== |
=={{header|Swift}}== |
||
< |
<syntaxhighlight lang="swift">import Foundation |
||
extension Numeric where Self: Strideable { |
extension Numeric where Self: Strideable { |
||
Line 1,109: | Line 1,109: | ||
print("π(10^\(i)) = \(π(n: n))") |
print("π(10^\(i)) = \(π(n: n))") |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 1,129: | Line 1,129: | ||
To sieve a billion numbers and then count the primes up to 10^k would take about 53 seconds in Wren so, as expected, the Legendre method represents a considerable speed up. |
To sieve a billion numbers and then count the primes up to 10^k would take about 53 seconds in Wren so, as expected, the Legendre method represents a considerable speed up. |
||
< |
<syntaxhighlight lang="ecmascript">import "/math" for Int |
||
var limit = 1e9.sqrt.floor |
var limit = 1e9.sqrt.floor |
||
Line 1,155: | Line 1,155: | ||
System.print("10^%(i) %(pi.call(n))") |
System.print("10^%(i) %(pi.call(n))") |
||
n = n * 10 |
n = n * 10 |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
Revision as of 17:45, 27 August 2022
![Task](http://static.miraheze.org/rosettacodewiki/thumb/b/ba/Rcode-button-task-crushed.png/64px-Rcode-button-task-crushed.png)
You are encouraged to solve this task according to the task description, using any language you may know.
The prime-counting function π(n) computes the number of primes not greater than n. Legendre was the first mathematician to create a formula to compute π(n) based on the inclusion/exclusion principle.
To calculate:
Define
φ(x, 0) = x φ(x, a) = φ(x, a−1) − φ(⌊x/pa⌋, a−1), where pa is the ath prime number.
then
π(n) = 0 when n < 2 π(n) = φ(n, a) + a - 1, where a = π(√n), n ≥ 2
The Legendre formula still requires the use of a sieve to enumerate primes; however it's only required to sieve up to the √n, and for counting primes, the Legendre method is generally much faster than sieving up to n.
- Task
Calculate π(n) for values up to 1 billion. Show π(n) for n = 1, 10, 100, ... 109.
For this task, you may refer to a prime number sieve (such as the Sieve of Eratosthenes or the extensible sieve) in an external library to enumerate the primes required by the formula. Also note that it will be necessary to memoize the results of φ(x, a) in order to have reasonable performance, since the recurrence relation would otherwise take exponential time.
C++
#include <cmath>
#include <iostream>
#include <unordered_map>
#include <vector>
std::vector<int> generate_primes(int limit) {
std::vector<bool> sieve(limit >> 1, true);
for (int p = 3, s = 9; s < limit; p += 2) {
if (sieve[p >> 1]) {
for (int q = s; q < limit; q += p << 1)
sieve[q >> 1] = false;
}
s += (p + 1) << 2;
}
std::vector<int> primes;
if (limit > 2)
primes.push_back(2);
for (int i = 1; i < sieve.size(); ++i) {
if (sieve[i])
primes.push_back((i << 1) + 1);
}
return primes;
}
class legendre_prime_counter {
public:
explicit legendre_prime_counter(int limit);
int prime_count(int n);
private:
int phi(int x, int a);
std::vector<int> primes;
std::unordered_map<int, std::unordered_map<int, int>> phi_cache;
};
legendre_prime_counter::legendre_prime_counter(int limit) :
primes(generate_primes(static_cast<int>(std::sqrt(limit)))) {}
int legendre_prime_counter::prime_count(int n) {
if (n < 2)
return 0;
int a = prime_count(static_cast<int>(std::sqrt(n)));
return phi(n, a) + a - 1;
}
int legendre_prime_counter::phi(int x, int a) {
if (a == 0)
return x;
auto& map = phi_cache[x];
auto i = map.find(a);
if (i != map.end())
return i->second;
int result = phi(x, a - 1) - phi(x / primes[a - 1], a - 1);
map[a] = result;
return result;
}
int main() {
legendre_prime_counter counter(1000000000);
for (int i = 0, n = 1; i < 10; ++i, n *= 10)
std::cout << "10^" << i << "\t" << counter.prime_count(n) << '\n';
}
- Output:
10^0 0 10^1 4 10^2 25 10^3 168 10^4 1229 10^5 9592 10^6 78498 10^7 664579 10^8 5761455 10^9 50847534
CoffeeScript
sorenson = require('sieve').primes # Sorenson's extensible sieve from task: Extensible Prime Generator
# put in outer scope to avoid recomputing the cache
memoPhi = {}
primes = []
isqrt = (x) -> Math.floor Math.sqrt x
pi = (n) ->
phi = (x, a) ->
y = memoPhi[[x,a]]
return y unless y is undefined
memoPhi[[x,a]] =
if a is 0 then x
else
p = primes[a - 1]
throw "You need to generate at least #{a} primes." if p is undefined
phi(x, a - 1) - phi(x // p, a - 1)
if n < 2
0
else
a = pi isqrt n
phi(n, a) + a - 1
maxPi = 1e9
gen = sorenson()
primes = while (p = gen.next().value) < isqrt maxPi then p
n = 1
for i in [0..9]
console.log "10^#{i}\t#{pi(n)}"
n *= 10
- Output:
10^0 0 10^1 4 10^2 25 10^3 168 10^4 1229 10^5 9592 10^6 78498 10^7 664579 10^8 5761455 10^9 50847534
Erlang
-mode(native).
-define(LIMIT, 1000000000).
main(_) ->
put(primes, array:from_list(primality:sieve(floor(math:sqrt(?LIMIT))))),
ets:new(memphi, [set, named_table, protected]),
output(0, 9).
nthprime(N) -> array:get(N - 1, get(primes)).
output(A, B) -> output(A, B, 1).
output(A, B, _) when A > B -> ok;
output(A, B, N) ->
io:format("10^~b ~b~n", [A, pi(N)]),
output(A + 1, B, N * 10).
pi(N) ->
Primes = get(primes),
Last = array:get(array:size(Primes) - 1, Primes),
if
N =< Last -> small_pi(N);
true ->
A = pi(floor(math:sqrt(N))),
phi(N, A) + A - 1
end.
phi(X, 0) -> X;
phi(X, A) ->
case ets:lookup(memphi, {X, A}) of
[] ->
Phi = phi(X, A-1) - phi(X div nthprime(A), A-1),
ets:insert(memphi, {{X, A}, Phi}),
Phi;
[{{X, A}, Phi}] -> Phi
end.
% Use binary search to count primes that we have already listed.
small_pi(N) ->
Primes = get(primes),
small_pi(N, Primes, 0, array:size(Primes)).
small_pi(_, _, L, H) when L >= (H - 1) -> L + 1;
small_pi(N, Primes, L, H) ->
M = (L + H) div 2,
P = array:get(M, Primes),
if
N > P -> small_pi(N, Primes, M, H);
N < P -> small_pi(N, Primes, 0, M);
true -> M + 1
end.
- Output:
10^0 1 10^1 4 10^2 25 10^3 168 10^4 1229 10^5 9592 10^6 78498 10^7 664579 10^8 5761455 10^9 50847534
Go
This runs in about 2.2 seconds which is surprisingly slow for Go. However, Go has the same problem as Wren in not being able to use lists for map keys. I tried using a 2-element array for the key but this was a second slower than the Cantor function.
The alternative of sieving a billion numbers and then counting them takes about 4.8 seconds using the present library function though this could be cut to 1.2 seconds using a third party library such as primegen.
package main
import (
"fmt"
"log"
"math"
"rcu"
)
var limit = int(math.Sqrt(1e9))
var primes = rcu.Primes(limit)
var memoPhi = make(map[int]int)
func cantorPair(x, y int) int {
if x < 0 || y < 0 {
log.Fatal("Arguments must be non-negative integers.")
}
return (x*x + 3*x + 2*x*y + y + y*y) / 2
}
func phi(x, a int) int {
if a == 0 {
return x
}
key := cantorPair(x, a)
if v, ok := memoPhi[key]; ok {
return v
}
pa := primes[a-1]
memoPhi[key] = phi(x, a-1) - phi(x/pa, a-1)
return memoPhi[key]
}
func pi(n int) int {
if n < 2 {
return 0
}
a := pi(int(math.Sqrt(float64(n))))
return phi(n, a) + a - 1
}
func main() {
for i, n := 0, 1; i <= 9; i, n = i+1, n*10 {
fmt.Printf("10^%d %d\n", i, pi(n))
}
}
- Output:
10^0 0 10^1 4 10^2 25 10^3 168 10^4 1229 10^5 9592 10^6 78498 10^7 664579 10^8 5761455 10^9 50847534
Java
import java.util.*;
public class LegendrePrimeCounter {
public static void main(String[] args) {
LegendrePrimeCounter counter = new LegendrePrimeCounter(1000000000);
for (int i = 0, n = 1; i < 10; ++i, n *= 10)
System.out.printf("10^%d\t%d\n", i, counter.primeCount((n)));
}
private List<Integer> primes;
private Map<Integer, Map<Integer, Integer>> phiCache = new HashMap<>();
public LegendrePrimeCounter(int limit) {
primes = generatePrimes((int)Math.sqrt((double)limit));
}
public int primeCount(int n) {
if (n < 2)
return 0;
int a = primeCount((int)Math.sqrt((double)n));
return phi(n, a) + a - 1;
}
private int phi(int x, int a) {
if (a == 0)
return x;
Map<Integer, Integer> map = phiCache.computeIfAbsent(x, k -> new HashMap<>());
Integer value = map.get(a);
if (value != null)
return value;
int result = phi(x, a - 1) - phi(x / primes.get(a - 1), a - 1);
map.put(a, result);
return result;
}
private static List<Integer> generatePrimes(int limit) {
boolean[] sieve = new boolean[limit >> 1];
Arrays.fill(sieve, true);
for (int p = 3, s = 9; s < limit; p += 2) {
if (sieve[p >> 1]) {
for (int q = s; q < limit; q += p << 1)
sieve[q >> 1] = false;
}
s += (p + 1) << 2;
}
List<Integer> primes = new ArrayList<>();
if (limit > 2)
primes.add(2);
for (int i = 1; i < sieve.length; ++i) {
if (sieve[i])
primes.add((i << 1) + 1);
}
return primes;
}
}
- Output:
10^0 0 10^1 4 10^2 25 10^3 168 10^4 1229 10^5 9592 10^6 78498 10^7 664579 10^8 5761455 10^9 50847534
Haskell
Memoization utilities:
data Memo a = Node a (Memo a) (Memo a)
deriving Functor
memo :: Integral a => Memo p -> a -> p
memo (Node a l r) n
| n == 0 = a
| odd n = memo l (n `div` 2)
| otherwise = memo r (n `div` 2 - 1)
nats :: Integral a => Memo a
nats = Node 0 ((+1).(*2) <$> nats) ((*2).(+1) <$> nats)
memoize :: Integral a => (a -> b) -> a -> b
memoize f = memo (f <$> nats)
memoize2 :: (Integral a, Integral b) => (a -> b -> c) -> a -> b -> c
memoize2 f = memoize (memoize . f)
memoList :: [b] -> Integer -> b
memoList = memo . mkList
where
mkList (x:xs) = Node x (mkList l) (mkList r)
where (l,r) = split xs
split [] = ([],[])
split [x] = ([x],[])
split (x:y:xs) = let (l,r) = split xs in (x:l, y:r)
Computation of Legendre function:
isqrt :: Integer -> Integer
isqrt n = go n 0 (q `shiftR` 2)
where
q = head $ dropWhile (< n) $ iterate (`shiftL` 2) 1
go z r 0 = r
go z r q = let t = z - r - q
in if t >= 0
then go t (r `shiftR` 1 + q) (q `shiftR` 2)
else go z (r `shiftR` 1) (q `shiftR` 2)
phi = memoize2 phiM
where
phiM x 0 = x
phiM x a = phi x (a-1) - phi (x `div` p a) (a - 1)
p = memoList (undefined : primes)
legendrePi :: Integer -> Integer
legendrePi n
| n < 2 = 0
| otherwise = phi n a + a - 1
where a = legendrePi (floor (sqrt (fromInteger n)))
main = mapM_ (\n -> putStrLn $ show n ++ "\t" ++ show (legendrePi (10^n))) [1..7]
λ> main 1 4 2 25 3 168 4 1229 5 9592 6 78498 7 664579 8 5761455 9 50847534
J
This is "almost a primitive" in J:
require'format/printf'
{{echo '10^%d: %d' sprintf y;p:inv 10^y}}"0 i.10
10^0: 0
10^1: 4
10^2: 25
10^3: 168
10^4: 1229
10^5: 9592
10^6: 78498
10^7: 664579
10^8: 5761455
10^9: 50847534
A complete implementation of the legendre prime counting function would be (1&p: + p:inv):
(1&p: + p:inv) 10
4
(1&p: + p:inv) 11
5
But we can simplify that to p:inv when we know that primes are not being tested.
jq
Adapted from Wren
Works with gojq, the Go implementation of jq (*)
The following implementation freely uses the following conveniences of jq syntax:
- the jq expression {$x} is an abbreviation for {"x" : $x}
- the jq expression {x} is an abbreviation for {"x" : .x}
The key-value pairs can be similarly specified in JSON objects with multiple keys.
(*) gojq struggles to get beyond [5,9592] without using huge amounts of memory; jq becomes excessively slow after [7,664579].
Preliminaries
# For the sake of the accuracy of integer arithmetic when using gojq:
def power($b): . as $in | reduce range(0;$b) as $i (1; . * $in);
# Input: $n, which is assumed to be positive integer,
# Output: an array of primes less than or equal to $n (e.g. 10|eratosthenes #=> [2,3,5,7]
def eratosthenes:
# erase(i) sets .[i*j] to false for integral j > 1
def erase(i):
if .[i] then
reduce range(2; (1 + length) / i) as $j (.; .[i * $j] = false)
else .
end;
(. + 1) as $n
| (($n|sqrt) / 2) as $s
| [null, null, range(2; $n)]
| reduce (2, 1 + (2 * range(1; $s))) as $i (.; erase($i))
| map(select(.)) ;
Phi and the Legendre Counting Function
def legendre:
(sqrt | floor + 1 | eratosthenes) as $primes
# Input: {x, a}
# Output: {phi: phi(x,a), memo} where .memo might have been updated
| def phi:
. as {x: $x, a: $a, memo: $memo}
| if $a == 0 then {phi: $x, $memo}
else "\($x),\($a)" as $ix
| $memo[ $ix ] as $m
| if $m then {phi: $m, $memo}
else .a += -1
| phi as {phi: $phi1, memo: $memo}
| (($x / $primes[$a - 1])|floor) as $x
| ({$x, a, $memo} | phi) as {phi: $phi2, memo: $memo}
| ($phi1 - $phi2) as $phi
| {$phi, $memo}
| .memo[$ix] = $phi
end
end;
def l:
. as $n
| if . < 2 then 0
else ($n|sqrt|floor|l) as $a
| ({x: $n, $a, memo: {}} | phi).phi + $a - 1
end;
l;
def task:
range(0;10)
| . as $i
| [., (10|power($i)|legendre)]
;
task
- Output:
[0,0] [1,4] [2,25] [3,168] [4,1229] [5,9592] [6,78498] [7,664579] <terminated>
Julia
Using the Julia Primes library
This would be faster with a custom sieve that only counted instead of making a bitset, instead of the more versatile library function.
using Primes
function primepi(N)
delta = round(Int, N^0.8)
return sum(i -> count(primesmask(i, min(i + delta - 1, N))), 1:delta:N)
end
@time for power in 0:9
println("10^", rpad(power, 5), primepi(10^power))
end
- Output:
10^0 0 10^1 4 10^2 25 10^3 168 10^4 1229 10^5 9592 10^6 78498 10^7 664579 10^8 5761455 10^9 50847534 1.935556 seconds (1.64 k allocations: 415.526 MiB, 2.15% gc time)
Using the method given in the task
Run twice to show the benefits of precomputing the library prime functions.
using Primes
const maxpi = 1_000_000_000
const memoφ = Dict{Vector{Int}, Int}()
const pₐ = primes(isqrt(maxpi))
function π(n)
function φ(x, a)
if !haskey(memoφ, [x, a])
memoφ[[x, a]] = a == 0 ? x : φ(x, a - 1) - φ(x ÷ pₐ[a], a - 1)
end
return memoφ[[x, a]]
end
n < 2 && return 0
a = π(isqrt(n))
return φ(n, a) + a - 1
end
@time for i in 0:9
println("10^", rpad(i, 5), π(10^i))
end
@time for i in 0:9
println("10^", rpad(i, 5), π(10^i))
end
- Output:
10^0 1 10^1 4 10^2 25 10^3 168 10^4 1229 10^5 9592 10^6 78498 10^7 664579 10^8 5761455 10^9 50847534 12.939211 seconds (69.44 M allocations: 3.860 GiB, 27.98% gc time, 1.58% compilation time) 10^0 1 10^1 4 10^2 25 10^3 168 10^4 1229 10^5 9592 10^6 78498 10^7 664579 10^8 5761455 10^9 50847534 0.003234 seconds (421 allocations: 14.547 KiB)
Mathematica /Wolfram Language
ClearAll[Phi,pi]
$RecursionLimit = 10^6;
Phi[x_, 0] := x
Phi[x_, a_] := Phi[x, a] = Phi[x, a - 1] - Phi[Floor[x/Prime[a]], a - 1]
pi[n_] := Module[{a}, If[n < 2, 0, a = pi[Floor[Sqrt[n]]]; Phi[n, a] + a - 1]]
Scan[Print[pi[10^#]] &, Range[0,9]]
- Output:
0 4 25 168 1229 9592 78498 664579 5761455 50847534
Nim
The program runs in 2.2 seconds on our laptop. Using int32 instead of naturals (64 bits on our 64 bits computer) saves 0.3 second. Using an integer rather than a tuple as key (for instance x shl 32 or a
) saves 0.4 second.
import math, strutils, sugar, tables
const
N = 1_000_000_000
S = sqrt(N.toFloat).int
var composite: array[3..S, bool]
for n in countup(3, S, 2):
if n * n > S: break
if not composite[n]:
for k in countup(n * n, S, 2 * n):
composite[k] = true
# Prime list. Add a dummy zero to start at index 1.
let primes = @[0, 2] & collect(newSeq, for n in countup(3, S, 2): (if not composite[n]: n))
var cache: Table[(Natural, Natural), Natural]
proc phi(x, a: Natural): Natural =
if a == 0: return x
let pair = (x, a)
if pair in cache: return cache[pair]
result = phi(x, a - 1) - phi(x div primes[a], a - 1)
cache[pair] = result
proc π(n: Natural): Natural =
if n <= 2: return 0
let a = π(sqrt(n.toFloat).Natural)
result = phi(n, a) + a - 1
var n = 1
for i in 0..9:
echo "π(10^$1) = $2".format(i, π(n))
n *= 10
- Output:
π(10^0) = 0 π(10^1) = 4 π(10^2) = 25 π(10^3) = 168 π(10^4) = 1229 π(10^5) = 9592 π(10^6) = 78498 π(10^7) = 664579 π(10^8) = 5761455 π(10^9) = 50847534
Perl
#!/usr/bin/perl
use strict; # https://rosettacode.org/wiki/Legendre_prime_counting_function
use warnings;
no warnings qw(recursion);
use ntheory qw( nth_prime prime_count );
my (%cachephi, %cachepi);
sub phi
{
return $cachephi{"@_"} //= do {
my ($x, $aa) = @_;
$aa <= 0 ? $x : phi($x, $aa - 1) - phi(int $x / nth_prime($aa), $aa - 1) };
}
sub pi
{
return $cachepi{$_[0]} //= do {
my $n = shift;
$n < 2 ? 0 : do{ my $aa = pi(int sqrt $n); phi($n, $aa) + $aa - 1 } };
}
print "e n Legendre ntheory\n",
"- - -------- -------\n";
for (1 .. 9)
{
printf "%d %12d %10d %10d\n", $_, 10**$_, pi(10**$_), prime_count(10**$_);
}
- Output:
e n Legendre ntheory - - -------- ------- 1 10 4 4 2 100 25 25 3 1000 168 168 4 10000 1229 1229 5 100000 9592 9592 6 1000000 78498 78498 7 10000000 664579 664579 8 100000000 5761455 5761455 9 1000000000 50847534 50847534
Phix
with javascript_semantics -- -- While a phix dictionary can handle keys of {x,a}, for this -- task performance was dreadful (7,612,479 entries, >3mins), -- so instead memophix maps known x to an index to memophia -- which holds the full [1..a] for each x, dropping to a much -- more respectable (albeit not super-fast) 14s -- constant memophix = new_dict() sequence memophia = {} -- 1..a (max 3401) for each x function phi(integer x, a) if a=0 then return x end if integer adx = getd(x,memophix), res if adx=NULL then memophia = append(memophia,repeat(-1,a)) adx = length(memophia) setd(x,adx,memophix) else object ma = memophia[adx] integer l = length(ma) if a>l then memophia[adx] = 0 -- kill refcount memophia[adx] = ma & repeat(-1,a-l) else res = ma[a] if res>=0 then return res end if end if ma = 0 -- kill refcount end if res = phi(x, a-1) - phi(floor(x/get_prime(a)), a-1) memophia[adx][a] = res return res end function function pi(integer n) if n<2 then return 0 end if integer a = pi(floor(sqrt(n))) return phi(n, a) + a - 1 end function atom t0 = time() for i=0 to iff(platform()=JS?8:9) do printf(1,"10^%d %d\n",{i,pi(power(10,i))}) -- printf(1,"10^%d %d\n",{i,length(get_primes_le(power(10,i)))}) end for ?elapsed(time()-t0)
The commented-out length(get_primes_le()) gives the same results, in about the same time, and in both cases re-running the main loop a second time finishes near-instantly.
- Output:
10^0 0 10^1 4 10^2 25 10^3 168 10^4 1229 10^5 9592 10^6 78498 10^7 664579 10^8 5761455 10^9 50847534 "14.0s"
(It is about 4 times slower under pwa/p2js so output is limited to 10^8, unless you like staring at a blank screen for 52s)
Picat
Performance is surprisingly good, considering that the Picat's library SoE is very basic and that I recompute the base primes each time pi() is called. This version takes 23 seconds on my laptop, outperforming the CoffeeScript version which takes 30 seconds.
table(+, +, nt)
phi(X, 0, _) = X.
phi(X, A, Primes) = phi(X, A - 1, Primes) - phi(X // Primes[A], A - 1, Primes).
pi(N) = Count, N < 2 => Count = 0.
pi(N) = Count =>
M = floor(sqrt(N)),
A = pi(M),
Count = phi(N, A, math.primes(M)) + A - 1.
main =>
N = 1,
foreach (K in 0..9)
writef("10^%w %w%n", K, pi(N)),
N := N * 10.
- Output:
10^0 0 10^1 4 10^2 25 10^3 168 10^4 1229 10^5 9592 10^6 78498 10^7 664579 10^8 5761455 10^9 50847534
PicoLisp
Uses code from the Sieve of Eratosthenes to generate primes.
(setq Legendre-Max (** 10 9))
(load "plcommon/eratosthenes.l") # see task "Sieve of Eratosthenes, 2x3x5x7 wheel version.
# Create an index tree of the first N primes up to √Legendre-Max
(setq Sieve (sieve (sqrt Legendre-Max)))
(balance 'Primes (mapcar 'cons (range 1 (length Sieve)) Sieve))
(de prime (N) (cdr (lup Primes N)))
(de ϕ (X A)
(cache '(NIL) (cons X A)
(if (=0 A)
X
(- (ϕ X (dec A)) (ϕ (/ X (prime A)) (dec A))))))
(de π (N)
(if (< N 2)
0
(let A (π (sqrt N))
(+ (ϕ N A) A -1))))
(for N 10
(prinl "10\^" (dec N) "^I" (π (** 10 (dec N)))))
(bye)
- Output:
10^0 0 10^1 4 10^2 25 10^3 168 10^4 1229 10^5 9592 10^6 78498 10^7 664579 10^8 5761455 10^9 50847534
Python
Using primesieve
module (pip install primesieve
).
from primesieve import primes
from math import isqrt
from functools import cache
p = primes(isqrt(1_000_000_000))
@cache
def phi(x, a):
res = 0
while True:
if not a or not x:
return x + res
a -= 1
res -= phi(x//p[a], a) # partial tail recursion
def legpi(n):
if n < 2: return 0
a = legpi(isqrt(n))
return phi(n, a) + a - 1
for e in range(10):
print(f'10^{e}', legpi(10**e))
- Output:
10^0 0 10^1 4 10^2 25 10^3 168 10^4 1229 10^5 9592 10^6 78498 10^7 664579 10^8 5761455 10^9 50847534
Raku
Seems like an awful lot of pointless work. Using prime sieve anyway, why not just use it?
use Math::Primesieve;
my $sieve = Math::Primesieve.new;
say "10^$_\t" ~ $sieve.count: exp($_,10) for ^10;
say (now - INIT now) ~ ' elapsed seconds';
- Output:
10^0 0 10^1 4 10^2 25 10^3 168 10^4 1229 10^5 9592 10^6 78498 10^7 664579 10^8 5761455 10^9 50847534 0.071464489 elapsed seconds
Sidef
func legendre_phi(x, a) is cached {
return x if (a <= 0)
__FUNC__(x, a-1) - __FUNC__(idiv(x, a.prime), a-1)
}
func legendre_prime_count(n) is cached {
return 0 if (n < 2)
var a = __FUNC__(n.isqrt)
legendre_phi(n, a) + a - 1
}
print("e n Legendre builtin\n",
"- - -------- -------\n")
for n in (1 .. 9) {
printf("%d %12d %10d %10d\n", n, 10**n,
legendre_prime_count(10**n), prime_count(10**n))
assert_eq(legendre_prime_count(10**n), prime_count(10**n))
}
- Output:
e n Legendre builtin - - -------- ------- 1 10 4 4 2 100 25 25 3 1000 168 168 4 10000 1229 1229 5 100000 9592 9592 6 1000000 78498 78498 7 10000000 664579 664579
Alternative implementation of the Legendre phi function, by counting k-rough numbers <= n.
func rough_count (n,k) {
# Count of k-rough numbers <= n.
func (n,p) is cached {
if (p > n.isqrt) {
return 1
}
if (p == 2) {
return (n >> 1)
}
if (p == 3) {
var t = idiv(n,3)
return (t - (t >> 1))
}
var u = 0
var t = idiv(n,p)
for (var q = 2; q < p; q.next_prime!) {
var v = __FUNC__(t - (t % q), q)
if (v == 1) {
u += prime_count(q, p-1)
break
}
u += v
}
t - u
}(n*k, k)
}
func legendre_phi(n, a) {
rough_count(n, prime(a+1))
}
func legendre_prime_count(n) is cached {
return 0 if (n < 2)
var a = __FUNC__(n.isqrt)
legendre_phi(n, a) + a - 1
}
print("e n Legendre builtin\n",
"- - -------- -------\n")
for n in (1 .. 9) {
printf("%d %12d %10d %10d\n", n, 10**n,
legendre_prime_count(10**n), prime_count(10**n))
assert_eq(legendre_prime_count(10**n), prime_count(10**n))
}
Swift
import Foundation
extension Numeric where Self: Strideable {
@inlinable
public func power(_ n: Self) -> Self {
return stride(from: 0, to: n, by: 1).lazy.map({_ in self }).reduce(1, *)
}
}
func eratosthenes(limit: Int) -> [Int] {
guard limit >= 3 else {
return limit < 2 ? [] : [2]
}
let ndxLimit = (limit - 3) / 2 + 1
let bufSize = ((limit - 3) / 2) / 32 + 1
let sqrtNdxLimit = (Int(Double(limit).squareRoot()) - 3) / 2 + 1
var cmpsts = Array(repeating: 0, count: bufSize)
for ndx in 0..<sqrtNdxLimit where (cmpsts[ndx >> 5] & (1 << (ndx & 31))) == 0 {
let p = ndx + ndx + 3
var cullPos = (p * p - 3) / 2
while cullPos < ndxLimit {
cmpsts[cullPos >> 5] |= 1 << (cullPos & 31)
cullPos += p
}
}
return (-1..<ndxLimit).compactMap({i -> Int? in
if i < 0 {
return 2
} else {
if cmpsts[i >> 5] & (1 << (i & 31)) == 0 {
return .some(i + i + 3)
} else {
return nil
}
}
})
}
let primes = eratosthenes(limit: 1_000_000_000)
func φ(_ x: Int, _ a: Int) -> Int {
struct Cache {
static var cache = [String: Int]()
}
guard a != 0 else {
return x
}
guard Cache.cache["\(x),\(a)"] == nil else {
return Cache.cache["\(x),\(a)"]!
}
Cache.cache["\(x),\(a)"] = φ(x, a - 1) - φ(x / primes[a - 1], a - 1)
return Cache.cache["\(x),\(a)"]!
}
func π(n: Int) -> Int {
guard n > 2 else {
return 0
}
let a = π(n: Int(Double(n).squareRoot()))
return φ(n, a) + a - 1
}
for i in 0..<10 {
let n = 10.power(i)
print("π(10^\(i)) = \(π(n: n))")
}
- Output:
π(10^0) = 0 π(10^1) = 4 π(10^2) = 25 π(10^3) = 168 π(10^4) = 1229 π(10^5) = 9592 π(10^6) = 78498 π(10^7) = 664579 π(10^8) = 5761455 π(10^9) = 50847534
Wren
This runs in about 5.7 seconds which is not too bad for the Wren interpreter. As map keys cannot be lists, the Cantor pairing function has been used to represent [x, a] which is considerably faster than using a string based approach for memoization.
To sieve a billion numbers and then count the primes up to 10^k would take about 53 seconds in Wren so, as expected, the Legendre method represents a considerable speed up.
import "/math" for Int
var limit = 1e9.sqrt.floor
var primes = Int.primeSieve(limit)
var memoPhi = {}
var phi // recursive function
phi = Fn.new { |x, a|
if (a == 0) return x
var key = Int.cantorPair(x, a)
if (memoPhi.containsKey(key)) return memoPhi[key]
var pa = primes[a-1]
return memoPhi[key] = phi.call(x, a-1) - phi.call((x/pa).floor, a-1)
}
var pi // recursive function
pi = Fn.new { |n|
if (n < 2) return 0
var a = pi.call(n.sqrt.floor)
return phi.call(n, a) + a - 1
}
var n = 1
for (i in 0..9) {
System.print("10^%(i) %(pi.call(n))")
n = n * 10
}
- Output:
10^0 0 10^1 4 10^2 25 10^3 168 10^4 1229 10^5 9592 10^6 78498 10^7 664579 10^8 5761455 10^9 50847534