Legendre prime counting function: Difference between revisions

From Rosetta Code
Content added Content deleted
(Added Java solution)
Line 507: Line 507:
π(10^8) = 5761455
π(10^8) = 5761455
π(10^9) = 50847534</pre>
π(10^9) = 50847534</pre>

=={{header|Perl}}==
<lang 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**$_);
}</lang>
{{out}}
<pre>
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
</pre>


=={{header|Phix}}==
=={{header|Phix}}==

Revision as of 19:25, 22 August 2021

Task
Legendre prime counting function
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++

<lang cpp>#include <cmath>

  1. include <iostream>
  2. include <unordered_map>
  3. 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';

}</lang>

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

<lang CoffeeScript> sorenson = require('sieve').primes # Sorenson's extensible sieve from task: Extensible Prime Generator

  1. put in outer scope to avoid recomputing the cache

memoPhi = {} primes = []

isqrt = (x) -> Math.floor Math.sqrt x

pi = (n) ->

   phi = (x, a) ->
       y = memoPhix,a
       return y unless y is undefined
       memoPhix,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

</lang>

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

<lang 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.

</lang>

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

Translation of: Wren
Library: Go-rcu

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. <lang go>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))
   }

}</lang>

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

<lang 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;
   }

}</lang>

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

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. <lang julia>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

</lang>

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

Translation of: CoffeeScript

Run twice to show the benefits of precomputing the library prime functions. <lang julia>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

</lang>

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

<lang Mathematica>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]]</lang>

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. <lang Nim>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
  1. 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</lang>
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

<lang 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**$_);
 }</lang>
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. <lang Picat> 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.

</lang>

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? <lang perl6>use Math::Primesieve;

my $sieve = Math::Primesieve.new;

say "10^$_\t" ~ $sieve.count: exp($_,10) for ^10;

say (now - INIT now) ~ ' elapsed seconds';</lang>

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

Wren

Library: Wren-math

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. <lang ecmascript>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

}</lang>

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