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

%%