Prime Counting Function

:-dynamic(pi/2).

% Finds the number of prime numbers less than or equal to a given number
%   using Lehmer's formula.
% e.g. lehmer(100000000, X) gives X=5761455
lehmer(1, 0):-!.
lehmer(2, 1):-!.
lehmer(X, Count):-
  S is int(sqrt(X+0.5)),             % x^(1/2)
  lehmer_db(S, B),                   % pi(x^(1/2))
  T is int(aln(ln(X+0.5)/3)),        % x^(1/3)
  lehmer_db(T, C),                   % pi(x^(1/3))
  U is int(sqrt(S+0.5)),             % x^(1/4)
  lehmer_db(U, A),                   % pi(x^(1/4))
  phi(X, A, Phi),                    % phi(x,a)
  Sum is Phi + (B+A-2)*(B-A+1) // 2, % phi(x,a) + (b+a-2)(b-a+1)/2
  primes(B, Ps0),                    % The first B primes
  first_n(A, Ps0, _, Ps),            % Primes starting at the (A+1)-th prime
  lehmer_1(A, B, C, X, Ps, Sum, Count).

lehmer_1(B, B, _, _, _, Sum, Sum):-!.
lehmer_1(I, B, C, X, [P|Ps], Sum0, Sum):-
  I < C, !,
  W is X // P,
  lehmer_db(W, Count),               % pi(x/p_i)       
  Sum1 is Sum0 - Count,
  T is int(sqrt(W+0.5)),
  lehmer_db(T, BI),                  % pi(sqrt(x/p_i))       
  lehmer_2(I, BI, W, [P|Ps], Sum1, Sum2),
  I1 is I + 1,
  lehmer_1(I1, B, C, X, Ps, Sum2, Sum).
lehmer_1(I, B, C, X, [P|Ps], Sum0, Sum):-
  W is X // P,
  lehmer_db(W, Count),               % pi(x/p_i)       
  Sum1 is Sum0 - Count,
  I1 is I + 1,
  lehmer_1(I1, B, C, X, Ps, Sum1, Sum).

lehmer_2(BI, BI, _, _, Sum, Sum):-!.
lehmer_2(J, BI, W, [P|Ps], Sum0, Sum):-
  WP is W // P,
  lehmer_db(WP, Count),              % pi(x/p_i/p_j)    
  Sum1 is Sum0 - Count + J,
  J1 is J + 1,
  lehmer_2(J1, BI, W, Ps, Sum1, Sum).

% Only calculates pi(x) if it is not already in the database.
lehmer_db(X, Count):-
  pi(X, Count), !.      % pi(x) is in the database
lehmer_db(X, Count):-   % pi(x) is not in the database
  lehmer(X, Count),
  assert(pi(X, Count)). % Put pi(x) in the database

% Finds the number of integers less than or equal to x, not divisible by any of
% the primes p1,p2,...,pa:
%   phi(x, a) = phi(x, a-1) - phi(int(x/p_a), a-1) 
% So phi(x, a) = x - sum i=0 to a-1 of phi(int(x/p_a), i)
% e.g. phi(103, 4, X) gives X=24
phi(0, _, 0):-!.
phi(X, _, 1):-X < 3, !.
phi(X, 0, X):-!.
phi(X, 1, Y):-!,
  Y is (X + 1) // 2.
phi(X, A, Y):-
% Enough primes to calculate lehmer(2147483647, X)
  Ps=[5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,
      101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,
      191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,
      281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,
      389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,
      491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,
      607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,
      719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,
      829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,947,
      953,967,971,977,983,991,997,1009,1013,1019,1021,1031,1033,1039,1049,
      1051,1061,1063,1069,1087,1091,1093,1097,1103,1109,1117,1123,1129,1151,
      1153,1163,1171,1181,1187,1193,1201,1213,1217,1223,1229,1231,1237,1249,
      1259,1277,1279,1283,1289],
  X0 is X - (X // 2) - ((X // 3 + 1) // 2),
  phi_1(2, A, Ps, X, X0, Y).

phi_1(A, A, _, _, Y, Y):-!.
phi_1(I, A, [P|Ps], X, Y0, Y):-
  X1 is X // P,
  phi(X1, I, X2),
  Y1 is Y0 - X2,
  I1 is I + 1,
  phi_1(I1, A, Ps, X, Y1, Y).

% Returns the first N and the remaining elements of Xs.
first_n(0, Xs, [], Xs):-!.
first_n(I, [X|Xs], [X|Ys], Zs):-
  I1 is I - 1,
  first_n(I1, Xs, Ys, Zs).

% Returns an ordered list of the first N prime numbers.                                                                
primes(0, []):-!.
primes(N, [2|Ps]):-primes_1(N, 3, Ps).

primes_1(1, _, []):-!.
primes_1(I, P, [P|Ps]):-
  is_prime(P), !,
  P1 is P + 2,
  I1 is I - 1,
  primes_1(I1, P1, Ps).
primes_1(I, P, Ps):-
  P1 is P + 2,
  primes_1(I, P1, Ps).

% Is true if P is a prime number.
is_prime(2):-!.
is_prime(P):-
%  P > 2,
%  P mod 2 =:= 1,
  SqrtP is int(sqrt(P)) + 1,
  is_prime_1(3, SqrtP, P).

is_prime_1(I, SqrtP, _):-
  I >= SqrtP, !.
is_prime_1(I, SqrtP, P):-
  P mod I > 0,
  I1 is I + 2,
  is_prime_1(I1, SqrtP, P).

LPA Index     Home Page

Valid HTML 4.0!