arun@bach.ces.cwru.edu (Arun Lakhotia) (12/21/88)
I just derived a program segment that checks if a number is prime from
some pieces code posted by Thomas Sj|land. The code was posted in
reference to finding perfect numbers, where a perfect number can be
stated in terms of (merseme) prime numbers.
The final program (Program 4) was derived in three steps from Program
1. The latter program (except for isprime/1) has been extracted from
Thomas Sj|land's code. The motivation for the transformation was to
improve efficiency. The actual transformation steps applied were
guided by intuition.
%%
%% program: 1
%%
isprime(X) :- divisors(X, [1]).
% divisors(I, D) -- D is a list of divisors of I (excluding I itself).
divisors(I, D) :- J is I//2+1, upto(J, L), divisors(I, L, D).
divisors(_I, [], []).
divisors(I, [H|T], L) :-
(divides(I, H) -> L=[H|R] ; L=R),
divisors(I, T, R).
% upto(N, L) -- L is a list of integers from 1..N-1
upto(N, L) :- upto(N, 1, L).
upto(N, N, []).
upto(N, I, [I|T]) :- I<N, I1 is I+1 , upto(N, I1, T).
divides(I, D) :- I is I//D*D.
%%
%% end of program 1
%%
DISCUSSION
Program 1 first constructs a list of 1..I//2 integers and then extracts
from this list the numbers that are divisors of I. The process of
creating the list of 1..I//2 is extraneous and may be gotten rid of.
This is done by first unfolding `upto(J, L)', introducing a definition
clause
divisors(M, N, I, L, D) :- upto(N, I, L), divisors(I, L, D).
and creating a divisors/5 using Tamaki-Sato's fold/unfold procedure.
Next `recognizing' that the fourth argument 'L' is of no use and
mapping all predicates divisors/5 to divisors/4 as follows
definition(M, N, I, L, D) --> definition(M, N, I, D).
This could have been done in the previous step itself, but then the
variable 'L' would have become local and that sometimes creates
problems when folding.
%%
%% program: 2
%%
isprime(X) :- divisors(X, [1]).
divisors(M, D) :- J is M//2+1, divisors(M, J, 1, D).
divisors(M, N, N, []).
divisors(M, N, I, L) :-
I < N,
I1 is I+1,
(divides(M, I) -> L=[I|R] ; L=R), divisors(M, N, I1, R).
divides(I, D) :- I is I//D*D.
%%
%% end of program 2
%%
DISCUSSION
Program 2 checks if a number is prime by first getting the list of all
divisors of the number, and then checking if that list contains just a
single element '1'. Clearly the process of finding all the divisors,
is a waste as all we need to make sure is their is no other divisor
but 1.
The following program is modified to check the value of the divisor
immediately after it is found. It better be 1 or else isprime_x/3
fails. Also realizing that we do not have any need for the list of
divisors we simply get rid of that argument.
I do not know how to explain this step in terms of known
transformation. Would appreciate input towards that.
%%
%% program 3
%%
isprime(M) :- J is M//2+1, isprime_x(M, J, 1).
isprime_x(M, N, N).
isprime_x(M, N, I) :-
I < N,
I1 is I+1,
(divides(M, I) ->I=1; true), isprime_x(M, N, I1).
divides(I, D) :- I is I//D*D.
%%
%% end of program 3
%%
DISCUSSION
The above program checks for prime-ness of a number N by testing if
any of the numbers from 1..N//2+1 is its divisor. It does so even if
the number is even. This is clearly redundant and can be avoided by
checking early enough if the given number is even (other than 2).
Also checking whether an even number is a divisor for an odd number is
extraneous. This may be avoided by incrementing the sequence in steps
of 2 instead of 1.
The following program uses these observations.
%%
%% program 4
%%
isprime(2).
isprime(M) :-
\+ divides(M, 2),
J is M//2+1,
isprime_x(M, J, 1).
isprime_x(M, N, I) :-
N - I =< 1.
isprime_x(M, N, I) :-
I =< N,
(divides(M, I) -> I = 1 ; true),
I1 is I+2,
isprime_x(M, N, I1).
divides(I, D) :- I is I//D*D.
%%
%% end of program 4
%%