Source code in the prolog programming language

  prime(P) :-
    pascal([1,P|Xs]),
    append(Xs, [1], Rest),
    forall( member(X,Xs), 0 is X mod P).


% To generate the n-th row of a Pascal triangle
% pascal(+N, Row)
pascal(0, [1]).
pascal(N, Row) :-
  N > 0, optpascal( [1, N|Xs] ),
  !,
  pascalize( [1, N|Xs], Row ).

pascalize( Opt, Row ) :-
  % if Opt ends in a pair, then peel off the pair:
  ( append(X, [R,R], Opt) -> true ; append(X, [R], Opt) ),
  reverse(X, Rs), 
  append( Opt, Rs, Row ).

% optpascal(-X) generates optpascal lines:
optpascal(X) :-
  optpascal_successor( [], X).

% optpascal_successor(+P, -Q) is true if Q is an optpascal list beneath the optpascal list P:
optpascal_successor(P, Q) :-
  optpascal(P, NextP),
  (Q = NextP ; optpascal_successor(NextP, Q)).

% optpascal(+Row, NextRow) is true if Row and NextRow are adjacent rows in the Pascal triangle.
% optpascal(+Row, NextRow) where the optpascal representation is used
optpascal(X, [1|Y]) :-
  add_pairs(X, Y).

% add_pairs(+OptPascal, NextOptPascal) is a helper function for optpascal/2.
% Given one OptPascal list, it generates the next by adding adjacent
% items, but if the last two items are unequal, then their sum is
% repeated.  This is intended to be a deterministic predicate, and to
% avoid a probable compiler limitation, we therefore use one cut.
add_pairs([], []).
add_pairs([X], [X]).
add_pairs([X,Y], Ans) :- 
  S is X + Y,
  (X = Y -> Ans=[S] ; Ans=[S,S]),
  !.  % To overcome potential limitation of compiler

add_pairs( [X1, X2, X3|Xs], [S|Ys]) :-
  S is X1 + X2,
  add_pairs( [X2, X3|Xs], Ys).


%%% Task 1: "A method to generate the coefficients of (1-X)^p"

coefficients(N, Coefficients) :-
  pascal(N, X),
  alternate_signs(X, Coefficients).

alternate_signs( [], [] ).
alternate_signs( [A], [A] ).
alternate_signs( [A,B | X], [A, MB | Y] ) :-
  MB is -B,
  alternate_signs(X,Y).

%%% Task 2. "Show here the polynomial expansions of (x − 1)p for p in the range 0 to at least 7, inclusive."

coefficients(Coefficients) :-
  optpascal( Opt), 
  pascalize( Opt, Row ),
  alternate_signs(Row, Coefficients).


%  As required by the problem statement, but necessarily very inefficient:
:- between(0, 7, N), coefficients(N, Coefficients), writeln(Coefficients), fail ; true.


%%% Task 3. Use the previous function in creating [sic]
%%% another function that when given p returns whether p is prime
%%% using the AKS test.

% Even for testing whether a given number, N, is prime, 
% this approach is inefficient, but here is a Prolog implementation:
   
   prime_test_per_requirements(N) :-
     coefficients(N, [1|Coefficients]),
     append(Cs, [_], Coefficients),
     forall( member(C, Cs), 0 is C mod N).


   prime(N) :- optpascal([1,N|Xs]), forall( member(X,Xs), 0 is X mod N).


%%% Task 4. Use your AKS test to generate a list of all primes under 35.

:-   prime(N), (N < 35 -> write(N), write(' '), fail ; nl).

% Output: 1 2 3 5 7 11 13 17 19 23 29 31 

%%% Task 5. As a stretch goal, generate all primes under 50.

:-  prime(N), (N < 50 -> write(N), write(' '), fail ; nl).

% Output: 1 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47


main :- task1(8), nl, task2(50), halt.

task1(N) :-
    pascal(Z),
    length(Rows, N),
    prefix(Rows, Z),
    forall(member(Row, Rows),
        (length(Row, K), succ(DecK, K),
        binomial(x, -1, Row, Expr),
        format("(x-1)**~w = ~w~n", [DecK, Expr]))).

task2(Upto) :-
    primes_upto(Upto, Ps),
    format("The primes upto ~w (via AKS) are: ~p~n", [Upto, Ps]).

pascal(Lz) :-
    lazy_list(pascal_row, [], Lz).

pascal_row([], R1, R1) :- R1 = [1], !.
pascal_row(R0, R1, R1) :-
    sum_adj(R0, Next), R1 = [1|Next].

sum_adj(L, L) :- L = [_], !.
sum_adj([A|As], [C|Cs]) :-
    As = [B|_], C is A + B,
    sum_adj(As, Cs).

% First part of task -- create textual representation of (x-1)^n
%  here we generate expression trees
%
binomial(A, B, Coefs, Expr) :-
    length(Coefs, N), succ(DecN, N),
    binomial(B, DecN, A, 0, Coefs, Exp0),
    reduce(Exp0, Exp1),
    addition_to_subtraction(Exp1, Expr).

binomial(_, _, _, _, [], 0) :- !.
binomial(A, PowA, B, PowB, [N|Ns], Ts + T) :-
    T = N * A**PowA * B**PowB,
    IncPow is PowB + 1,
    DecPow is PowA - 1,
    binomial(A, DecPow, B, IncPow, Ns, Ts).

addition_to_subtraction(A + B, X) :-
    addition_to_subtraction(A, C),
    (make_positive(B, D) -> X = C - D; X = C + B), !.
addition_to_subtraction(X, X).

make_positive(N, Term) :- integer(N), N < 0, !, Term is -N.
make_positive(A*B, Term) :-
    make_positive(A, PosA),
    (PosA = 1 -> Term = B, !; Term = PosA*B).

reduce(A, C) :-
    simplify(A, B),
    (B = A -> C = A; reduce(B, C)).

simplify(_**0, 1) :- !.
simplify(1**_, 1) :- !.
simplify(-1**N, Z) :- integer(N), (0 is N /\ 1 -> Z = 1; Z = -1), !.
simplify(X**1, X) :- !.

simplify(0 + A, A) :- !.
simplify(A + 0, A) :- !.
simplify(A + B, C) :-
    integer(A),
    integer(B), !,
    C is A + B.
simplify(A + B, C + D) :- !,
    simplify(A, C),
    simplify(B, D).

simplify(0 * _, 0) :- !.
simplify(_ * 0, 0) :- !.
simplify(1 * A, A) :- !.
simplify(A * 1, A) :- !.
simplify(A * B, C) :-
    integer(A),
    integer(B), !,
    C is A * B.
simplify(A * B, C * D) :- !,
    simplify(A, C),
    simplify(B, D).

simplify(X, X).

% Second part of task -- Use the coefficients of Pascal's Triangle to check primality.
%

primerow([1, N| Rest]) :- primerow(N, Rest).

primerow(_, End) :- (End = []; End = [1]), !.
primerow(_, [A,A|_]) :- !.  % end when we've seen half the list.
primerow(N, [A|As]) :- A mod N =:= 0, primerow(N, As).

second([_,N|_], N).

primes_upto(N, Ps) :-
    pascal(Z),
    Z = [_, _ | Rows], % we only care about 2nd row on up. ([1,2,1])
    succ(DecN, N), length(CheckRows, DecN), prefix(CheckRows, Rows),
    include(primerow, CheckRows, PrimeRows),
    maplist(second, PrimeRows, Ps).

?- main.