#### Affine Planes and Projective Planes

This is based on a description given in "A Short Course in Combinatorial Designs" by Ian Anderson and Iiro Honkala.

/* affine_plane(N, Cover) is true if N is prime, and Cover is an order N     */
/*   affine plane, that is, a (N^2, N, 1) design, that is, a                 */
/*   c(N^2, N, 2, 2, N^2+N) cover.                                           */
/* e.g. affine_plane(3, X)                                                   */
/*        gives X=[[1,2,3],[1,4,7],[1,5,9],[1,6,8],[2,4,9],[2,5,8],          */
/*                 [2,6,7],[3,4,8],[3,5,7],[3,6,9],[4,5,6],[7,8,9]]          */
affine_plane(N, Cover):-
is_prime(N),
N1 is N - 1,
integers(0, N1, Is),
findall(Vector, vector(Is, Vector), Vectors),
Is=[_,_|Is1],
equivalence_classes(Vectors, Is1, N, EquivalenceClasses),
aff_cover(EquivalenceClasses, EquivalenceClasses, N, UnsortedCover),
sort(UnsortedCover, Cover).

/* aff_cover(Xs, AllXs, N, Iss) is true if each element of Iss is a list of  */
/*   the N indices (adjusted to start at 1) of those equivalence classes     */
/*   (a,b,c) in AllXs, for which a>0 and ad+be+cf=0 (using multiplication    */
/*   modulo N), where (d,e,f) is the corresponding equivalence class in Xs.  */
aff_cover([], _, _, []).
aff_cover([X|Xs], AllXs, N, [Is|Iss]):-
Index is -N,
aff_cover_1(N, AllXs, X, Index, N, Is),
Is=[_|_], !,
aff_cover(Xs, AllXs, N, Iss).
aff_cover([_|Xs], AllXs, N, Iss):-
aff_cover(Xs, AllXs, N, Iss).

aff_cover_1(0, _, _, _, _, []):-!.
aff_cover_1(J, [v(A,B,C)|Xs], X0, Index, N, [Index|Is]):-
A > 0,
X0=v(D,E,F),
(A*D + B*E + C*F) mod N =:= 0, !,
Index1 is Index + 1,
J1 is J - 1,
aff_cover_1(J1, Xs, X0, Index1, N, Is).
aff_cover_1(J, [_|Xs], X0, Index, N, Is):-
Index1 is Index + 1,
aff_cover_1(J, Xs, X0, Index1, N, Is).

/* projective_plane(N, Cover) is true if N is prime, and Cover is an order   */
/*   N projective plane, that is, a (N^2+N+1, N+1, 1) design, that is, a     */
/*   c(N^2+N+1, N+1, 2, 2, N^2+N+1) cover.                                   */
/* e.g. projective_plane(2, X)                                               */
/*        gives X=[[1,2,3],[1,4,5],[1,6,7],[2,4,6],[2,5,7],[3,4,7],[3,5,6]]  */
projective_plane(N, Cover):-
is_prime(N),
N1 is N - 1,
integers(0, N1, Is),
findall(Vector, vector(Is, Vector), Vectors),
Is=[_,_|Is1],
equivalence_classes(Vectors, Is1, N, EquivalenceClasses),
proj_cover(EquivalenceClasses, EquivalenceClasses, N, UnsortedCover),
sort(UnsortedCover, Cover).

/* proj_cover(Xs, AllXs, N, Iss) is true if each element of Iss is a list of */
/*   the N+1 indices of those equivalence classes (a,b,c) in AllXs, for      */
/*   which ad+be+cf=0 (using multiplication modulo N), where (d,e,f) is the  */
/*   corresponding equivalence class in Xs.                                  */
proj_cover([], _, _, []).
proj_cover([X|Xs], AllXs, N, [Is|Iss]):-
proj_cover_1(N, AllXs, X, 1, N, Is),
proj_cover(Xs, AllXs, N, Iss).

proj_cover_1(-1, _, _, _, _, []):-!.
proj_cover_1(J, [v(A,B,C)|Xs], X0, Index, N, [Index|Is]):-
X0=v(D,E,F),
(A*D + B*E + C*F) mod N =:= 0, !,
Index1 is Index + 1,
J1 is J - 1,
proj_cover_1(J1, Xs, X0, Index1, N, Is).
proj_cover_1(J, [_|Xs], X0, Index, N, Is):-
Index1 is Index + 1,
proj_cover_1(J, Xs, X0, Index1, N, Is).

/* vector(Is, Vector) is true if Vector is a 3-element non-zero vector       */
/*   containing elements from the list Is.                                   */
vector(Is, v(X,Y,Z)):-
member(X, Is),
member(Y, Is),
member(Z, Is),
X + Y + Z > 0.

/* equivalence_classes(Vs, Is, N, Ws) is true if Ws is the subset of Vs      */
/*   consisting of the first encountered vector of each equivalence class.   */
/*   Two vectors are equivalent if one of them can be obtained from the      */
/*   other by multiplying (modulo N) by an element of Is. (An equivalence    */
/*   class contains N-1 elements.)                                           */
equivalence_classes([], _, _, []).
equivalence_classes([V|Vs], Is, N, [V|Ws]):-
N2 is N - 2,
equivalence_classes_1(N2, Vs, V, Is, N, Vs1),
equivalence_classes(Vs1, Is, N, Ws).

equivalence_classes_1(0, Zs, _, _, _, Zs):-!.
equivalence_classes_1(I, [V|Vs], V0, Is, N, Zs):-
are_equivalent(Is, V0, V, N), !,
I1 is I - 1,
equivalence_classes_1(I1, Vs, V0, Is, N, Zs).
equivalence_classes_1(I, [V|Vs], V0, Is, N, [V|Zs]):-
equivalence_classes_1(I, Vs, V0, Is, N, Zs).

are_equivalent([I|_], v(A,B,C), v(D,E,F), N):-
D is (I * A) mod N,
E is (I * B) mod N,
F is (I * C) mod N, !.
are_equivalent([_|Is], V0, V, N):-
are_equivalent(Is, V0, V, N).

/* integers(M, N, Is) is true if Is is the list of integers from M to N      */
/*   inclusive.                                                              */
integers(N, N, [N]):-!.
integers(I, N, [I|Is]):-I < N, I1 is I + 1, integers(I1, N, Is).

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

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

/* member(X, Xs) is true if the element X is contained in the list Xs.       */
%member(X, [X|_]).
%member(X, [_|Xs]):-member(X, Xs).