LU Decomposition in Prolog

:- use_module(library(clpr)).

// various obvious bits of matrix machinary.

lu_decompose(M, L*U) :-
 dimensions(M,N*N),
 dimensions(L,N*N),
 dimensions(U,N*N),
 lower(L),
 upper(U),
 matrix_multiply(L,U,M).

% Unprogramming a permutation machine

stack([],[],[]) --> [].
stack([push|Instrs],[I|IS],SS) --> stack(Instrs,IS,[I|SS]).
stack([pop|Instrs],IS,[S|SS]) --> [S], stack(Instrs,IS,SS).

% ?- phrase(stack(Instructions, [1,2,3,4,5], []), [1,2,4,3,5], []).
% Instructions = [push, pop, push, pop, push, push, pop, pop, push, pop] .

Unscrambling without Prolog

Inductive Bit : Set :=
  H  | V
| TL | TR
| BL | BR.

Definition State :=
 (Bit * Bit * Bit *
  Bit *       Bit *
  Bit * Bit * Bit)%type.

Inductive Move : State -> State -> Set :=
| refl : forall S, Move S S

| top_left : forall A B C D E F G H, Move
  (A, B, C,
   D,    E,
   F, G, H)
  (B, C, A,
   D,    E,
   F, G, H)

| left_down : forall A B C D E F G H, Move
  (A, B, C,
   D,    E,
   F, G, H)
  (F, B, C,
   A,    E,
   D, G, H)

| right_down : forall A B C D E F G H, Move
  (A, B, C,
   D,    E,
   F, G, H)
  (A, B, H,
   D,    C,
   F, G, E)

| bottom_left : forall A B C D E F G H, Move
  (A, B, C,
   D,    E,
   F, G, H)
  (A, B, C,
   D,    E,
   G, H, F)

| trans_closure : forall A M Z,
  Move A M -> Move M Z -> Move A Z.


Definition solution state := Move state (TL, H, TR,
                                         V,     V,
                                         BL, H, BR).

Theorem prolog_example :
  solution (H, BR, BL,
            H,     TL,
            TR, V, V).


unfold solution.
eapply trans_closure.
apply top_left.
change (solution (BR, BL, H, H, TL, TR, V, V)).

Ltac whack direction :=
  unfold solution;
  eapply trans_closure;
    [ apply direction
    | match goal with
        | |- Move ?x _ => change (solution x)
      end
    ].

Ltac wtop := whack top_left.
Ltac wleft := whack left_down.
Ltac wright := whack right_down.
Ltac wbottom := whack bottom_left.

wtop.
wtop.

wtop.
wtop.
wbottom.
wbottom.
wleft.
wright.
wbottom.
wbottom.
wleft.
wright.

unfold solution.

apply refl.

Qed.

Theorem automated_example :
  solution (H, BR, BL,
            H,     TL,
            TR, V, V).

Inductive Marker (A : State) : Set := mkMarker.

Ltac save_state :=
  match goal with
    | |- solution ?state =>
      let seen := fresh "seen"
        in pose (seen := mkMarker state)
  end.

Ltac check_progress :=
  match goal with
    | |- solution ?State =>
      match goal with
        | _ : Marker State |- _ =>
          pose (dead_end := True)
        | _ => idtac
      end
  end.

Ltac swtop := save_state; whack top_left; check_progress.
Ltac swleft := save_state; whack left_down; check_progress.
Ltac swright := save_state; whack right_down; check_progress.
Ltac swbottom := save_state; whack bottom_left; check_progress.

Ltac swhack n :=
  (unfold solution; apply refl) ||
    match n with
      | 0 => fail
      | S ?m => swtop; pose (dead_end := False); clear dead_end; swhack m
      | S ?m => swleft; pose (dead_end := False); clear dead_end; swhack m
      | S ?m => swright; pose (dead_end := False); clear dead_end; swhack m
      | S ?m => swbottom; pose (dead_end := False); clear dead_end; swhack m
    end.

swhack 12.

Defined.

Print automated_example.

(*
automated_example = 
let seen := mkMarker (H, BR, BL, H, TL, TR, V, V) in
trans_closure (H, BR, BL, H, TL, TR, V, V) (BR, BL, H, H, TL, TR, V, V)
  (TL, H, TR, V, V, BL, H, BR) (top_left H BR BL H TL TR V V)
  (let dead_end := False in
   let seen0 := mkMarker (BR, BL, H, H, TL, TR, V, V) in
   trans_closure (BR, BL, H, H, TL, TR, V, V) (BL, H, BR, H, TL, TR, V, V)
     (TL, H, TR, V, V, BL, H, BR) (top_left BR BL H H TL TR V V)
     (let dead_end0 := False in
      let seen1 := mkMarker (BL, H, BR, H, TL, TR, V, V) in
      trans_closure (BL, H, BR, H, TL, TR, V, V) (TR, H, BR, BL, TL, H, V, V)
        (TL, H, TR, V, V, BL, H, BR) (left_down BL H BR H TL TR V V)
        (let dead_end1 := False in
         let seen2 := mkMarker (TR, H, BR, BL, TL, H, V, V) in
         ...
*)

Require Import List.
Inductive Move2 : Set := top_left2 | left_down2 | right_down2 | bottom_left2.
Fixpoint display_moves p q (m : Move p q) : list Move2 :=
  match m with
    | refl _ => nil
    | top_left _ _ _ _ _ _ _ _ => top_left2::nil
    | left_down _ _ _ _ _ _ _ _ => left_down2::nil
    | right_down _ _ _ _ _ _ _ _ => right_down2::nil
    | bottom_left _ _ _ _ _ _ _ _ => bottom_left2::nil
    | trans_closure _ _ _ x y => display_moves _ _ x ++ display_moves _ _ y
  end.

Eval compute in (display_moves _ _ automated_example).
(*
     = top_left2
       :: top_left2
       :: left_down2
       :: top_left2
       :: right_down2
       :: top_left2
       :: top_left2
       :: bottom_left2
       :: bottom_left2
       :: left_down2
       :: right_down2
       :: right_down2
       :: nil
     : list Move2
*)

Unscrambling with Prolog

%% I'm playing this Oxyd game (inside Enigma [link])
%% One of the puzzles is this scrambled configuration

starting(([-,⌋,⌊]
         ,[-,  ⌈]
         ,[⌉,|,|])).

%% And you've to hit against the sides to unscramble it into:

solution(([⌈,-,⌉]
         ,[|,  |]
         ,[⌊,-,⌋])).

%% Then it explodes and you win a spring
%% The moves you can do are:

move(top_left,
    ([A,B,C]
    ,[D,  E]
    ,[F,G,H]) --> ([B,C,A]
                  ,[D,  E]
                  ,[F,G,H])).

move(left_down,
    ([A,B,C]
    ,[D,  E]
    ,[F,G,H]) --> ([F,B,C]
                  ,[A,  E]
                  ,[D,G,H])).

move(right_down,
    ([A,B,C]
    ,[D,  E]
    ,[F,G,H]) --> ([A,B,H]
                  ,[D,  C]
                  ,[F,G,E])).

move(bottom_left,
    ([A,B,C]
    ,[D,  E]
    ,[F,G,H]) --> ([A,B,C]
                  ,[D,  E]
                  ,[G,H,F])).

%% I want Prolog to solve this automatically for me though,
%% My first attempt is to use recursion like so,

%: solve(Win,[]) :- solution(Win).
%: solve(State,[Move|Moves]) :- move(Move,State --> Next), solve(Next,Moves).

%! ?- starting(State), solve(State,Moves).
%! ERROR: Out of local stack

%% Since Prolog searches depth first, this just kept trying top_left,
%% again and again, getting nowhere, infact it seeing the same config.
%% so I can pass along a list of seen configurations now, and fail when
%% the same one is seen twice.

%: solve(Win,[],_) :- solution(Win).
%: solve(State,[Move|Moves],Seen) :- move(Move,State --> Next), not(member(Next,Seen)), solve(Next,Moves,[State|Seen]).

%! ?- starting(State), solve(State,Moves,[State]).
%! State = ([-, '⌋', '⌊'], [-, '⌈'], ['⌉', ('|'), ('|')]),
%! Moves = [top_left, top_left, left_down, top_left, top_left, left_down, ...]

%% I got a solution this time... but it's thousands of moves, I think the puzzle can be solved
%% Within at least 10 moves, so Prolog should fail if it gets into a branch in the search tree that deep

solve(Win,[],_,_) :- solution(Win).
solve(State,[Move|Moves],Seen,Limit) :-
 move(Move,State --> Next), not(member(Next,Seen)), succ(LowerLimit,Limit), solve(Next,Moves,[State|Seen],LowerLimit).

solve(Solution) :- starting(State), solve(State,Solution,[State],10).

%! ?- solve(Moves).
%! Moves = [top_left, top_left, bottom_left, bottom_left, left_down, right_down, bottom_left,
%!          bottom_left, left_down, right_down].

Lights out solver using CLP

:- use_module(library(clpfd)).

solve(Puzzle,Solution) :-
 dimensions(Puzzle,Rows,Cols), dimensions(Solution,Rows,Cols),
 flatten(Solution,Variables), Variables ins 0..1,
 cells_of(Puzzle --> configure(Solution)), label(Variables).

configure(Solution,(X,Y,o)) :- around(X,Y,Solution,Group), off(Group).
configure(Solution,(X,Y,*)) :- around(X,Y,Solution,Group), on(Group).

on(Lights) :- sum(Lights, #=, On), 1 #= On mod 2.
off(Lights) :- sum(Lights, #=, Off), 0 #= Off mod 2.

around(X,Y,Grid,[Up,Left,Down,Right,Middle]) :-
 Yu is Y-1, Xl is X-1,
 Yd is Y+1, Xr is X+1,
 ( Grid@([Yu,X]->Up) -> true ; Up = 0 ),
 ( Grid@([Y,Xl]->Left) -> true ; Left = 0 ),
 ( Grid@([Yd,X]->Down) -> true ; Down = 0 ),
 ( Grid@([Y,Xr]->Right) -> true ; Right = 0 ),
 ( Grid@([Y,X]->Middle) -> true ; Middle = 0 ).

game([[*,o,o,o,o,o,*,o],
      [*,*,o,*,*,*,*,*],
      [*,*,o,*,*,*,o,*],
      [*,o,o,o,o,o,o,o],
      [o,*,*,o,o,o,*,*],
      [o,*,o,o,o,o,o,o],
      [*,*,*,*,*,*,o,o],
      [o,*,*,*,*,*,*,*]]).





:- op(500,xfy,@).
:- op(1050,yfx,<-).

flatten([],[]).
flatten([X|Xs],Zs) :- flatten(Xs,Ys), append(X,Ys,Zs).

flip(P,Y,X) :- call(P,X,Y).
dimensions(Grid, Rows, Cols) :- length(Grid, Rows), maplist(flip(length,Cols),Grid).

E@([]->E)         :- !.
[X|_]@([0|Ns]->E) :- !, X@(Ns->E).
[_|X]@([N|Ns]->E) :- !, N > 0, succ(M, N), X@([M|Ns]->E).

_/E@([]<-E)               :- !.
[X|Xs]/[Y|Xs]@([0|Ns]<-E) :- !, X/Y@(Ns<-E).
[X|Xs]/[X|Ys]@([N|Ns]<-E) :- !, N > 0, succ(M, N), Xs/Ys@([M|Ns]<-E).

tuple(Board,Tuples) :- tuple((0,0),Board,Tuples).
tuple(_,[],[]).
tuple((X,Y),[Row|Rows],Tuples) :-
 tuprow((X,Y),Row-Tail,Tuples),
 succ(Y,Y1), tuple((X,Y1),Rows,Tail).
tuprow(_,[]-X,X).
tuprow((I,J),[E|Es]-X,[(I,J,E)|Xs]) :- succ(I,I1),tuprow((I1,J),Es-X,Xs).

cells_of(Grid --> Pred) :- tuple(Grid, Tuples), maplist(Pred, Tuples).

Pure Type Systems type checker in Prolog

They ( http://people.cs.uu.nl/andres/LambdaPi/index.html http://augustss.blogspot.com/2007/10/simpler-easier-in-recent-paper-simply.html ) said it was easy but actually it's very hard, of course these posts are very helpful though. My other reference is the book Lectures on the Curry-Howard Isomorphism (Studies in Logic and Foundations of Mathematics - Vol 149).

Assuming this doesn't have any bugs, does it..?

:- op(1200, xfx, :=).
:- op(150, xfx, ⊢).
:- op(140, xfx, :).
:- op(100, yfx, $).


sort(★).
sort(▢).

axiom(★,▢).

rule(★,★,★). % λ→
rule(★,▢,▢). % λP
rule(▢,★,★). % λ2
rule(▢,▢,▢). % λω


unfold(Term, TermU) :- unfold([],Term,TermU).
unfold(_, Sort, Sort) :- sort(Sort).
unfold(Γ, Ident, Value) :- member(Ident, Γ) -> Ident = Value ; ( Ident := ValueD ), unfold(ValueD, Value).
unfold(Γ, λ(Var:Type,Body), λ(Var:TypeU,BodyU)) :- unfold(Γ, Type, TypeU), unfold([Var|Γ], Body, BodyU).
unfold(Γ, π(Var:Type,Body), π(Var:TypeU,BodyU)) :- unfold(Γ, Type, TypeU), unfold([Var|Γ], Body, BodyU).
unfold(Γ, U $ V, P $ Q) :- unfold(Γ, U, P), unfold(Γ, V, Q).

alpha(X,Y) :- alpha([],X,Y).
alpha(_,Id,Id) :- !.
alpha(Γ,Var1,Var2) :- member(Var1=Var2,Γ).
alpha(Γ,λ(Var1:Type1,Body1), λ(Var2:Type2,Body2)) :- alpha(Γ,Type1,Type2), alpha([Var1=Var2|Γ],Body1,Body2).
alpha(Γ,π(Var1:Type1,Body1), π(Var2:Type2,Body2)) :- alpha(Γ,Type1,Type2), alpha([Var1=Var2|Γ],Body1,Body2).
alpha(Γ,U $ V,P $ Q) :- alpha(Γ,U,P), alpha(Γ,V,Q).

subst(X->Y,X,Y) :- !.
subst(_->_,Var,Var) :- atomic(Var).
subst(X->Y,λ(X:Type1,Body),λ(X:Type2,Body)) :- !, subst(X->Y,Type1,Type2).
subst(X->Y,λ(Z:Type1,Body1),λ(Z:Type2,Body2)) :- subst(X->Y,Type1,Type2), subst(X->Y,Body1,Body2).
subst(X->Y,π(X:Type1,Body),π(X:Type2,Body)) :- !, subst(X->Y,Type1,Type2).
subst(X->Y,π(Z:Type1,Body1),π(Z:Type2,Body2)) :- subst(X->Y,Type1,Type2), subst(X->Y,Body1,Body2).
subst(X->Y,U $ V,P $ Q) :- subst(X->Y,U,P), subst(X->Y,V,Q).

eval(λ(U:_,B) $ V, BetaVal) :- !, subst(U->V,B,Beta), eval(Beta,BetaVal).
eval(X $ Y, XY) :- !, eval(X,Xv), eval(Xv $ Y,XY).
eval(Value, Value).


type(Term, Type) :- unfold(Term,TermU), [] ⊢ TermU : Type.
hasType(Term, Type) :- type(Term, Type1), unfold(Type,TypeU), eval(TypeU, Type2), alpha(Type1,Type2).

_ ⊢ S1 : S2 :- axiom(S1,S2).
Γ ⊢ I : T :- once(member(I:T, Γ)).
Γ ⊢ λ(A:T,B) : π(A:Tv,C) :- Γ ⊢ T : S1, sort(S1), eval(T,Tv), [A:Tv|Γ] ⊢ B : C. %%, /* sanity checks: */ Γ ⊢ π(A:Tv,C) : S2, sort(S2).
Γ ⊢ π(A:T,B) : S :- Γ ⊢ T : S1, [A:T|Γ] ⊢ B : S2, rule(S1,S2,S).
Γ ⊢ P $ Q : T :- Γ ⊢ P : π(X:U,B), Γ ⊢ Q : V, alpha(U,V), eval(λ(X:U,B) $ Q,T).



id := λ(t:★, λ(x:t, x)).
idTy := π(x:★, π(e:x, x)).

const := λ(s:★, λ(t:★, λ(u:s, λ(v:t, u)))).
constTy := π(u:★, π(v:★, π(g0:u, π(g1:v, u)))).

%%bool := π(a:★, a -> a -> a).
bool := π(a:★, π(g1:a, π(g2:a, a))).
true := λ(a:★, λ(x:a, λ(y:a, x))).
false := λ(a:★, λ(x:a, λ(y:a, y))).

%%nat := π(u:★, u -> (u -> u) -> u).
nat := π(u:★, π(g1:u, π(g2:π(g3:u, u), u))).
zero := λ(a:★, λ(x:a, λ(f:π(g:a, a), x))).
%%succ := λ(n:nat, λ(a:★, λ(x:a, λ(f:(a -> a), f $ (n $ a $ x $ f))))).
succ := λ(n:nat, λ(a:★, λ(x:a, λ(f:π(g:a, a), f $ (n $ a $ x $ f))))).

add := λ(m:nat, λ(n:nat, m $ nat $ n $ succ)).

eq := λ(a:★, λ(x:a, λ(y:a,
       λ(p:π(x:a, ★), π(g:p$x, p$y))))).

%%eq(a:★, x:a, y:a) := λ(p:π(x:a, ★), π(g:p$x, p$y)).

Prolog CHR - Prime Seive One Liner

%% http://www.cs.kuleuven.be/~dtai/projects/CHR/

:- use_module(library(chr)).
:- chr_constraint prime/1.

prime(X) \ prime(Y) <=> 0 is Y mod X | true.

%% This simpagation rules states:
%%  given there are clauses prime(X) and prime(Y)
%%  such that 0 is Y mod X, remove the clause prime(Y).


%% To try it out:
primes(1). primes(N) :- prime(N), succ(M,N), primes(M).
%% ?- primes(50) % generate all candidates <= 50
%% % the CHR rule seives out those which are not prime

Executable BNF parser in Prolog

% <digit> ::= 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9
% <sign> ::= + | -
% <number> ::= [ <sign> ] <digit> { <digit> }


:- op(1120, xfx, ::=).
:- op(11, fx, <), op(13, xf, >).


<digit>  ::= 0 ; 1 ; 2 ; 3 ; 4 ; 5 ; 6 ; 7 ; 8 ; 9 .
<sign>   ::= (+) ; (-)                             .
<number> ::= [ <sign> ], <digit>, { <digit> }      .

<expr>   ::= <factor>, { ((+) ; (-)), <factor> }   .
<factor> ::= <term>,   { ((*) ; (/)), <term> }     .
<term>   ::= '(', <expr>, ')' ; <number>           .


parse(Rule)    --> { Rule ::= Body }, parse(Body).
parse(Atom)    --> { atomic(Atom), atom_codes(Atom, Codes) }, Codes.

parse((X , Y)) --> parse(X), parse(Y).

parse((X ; _)) --> parse(X).
parse((_;Y;Z)) --> parse(Y ; Z).
parse((_ ; Z)) --> { Z \= (_ ; _) }, parse(Z).

parse([X])     --> parse(X).
parse([_])     --> {}.

parse({X})     --> parse(X), parse({X}).
parse({_})     --> [].


% phrase(parse(<expr>), "-5*(73+7)").

Compare with a C implementation -- http://cvs.savannah.gnu.org/viewvc/bnf/bnf/src/grio.c?view=markup

An Embedded ALGOL-like language in Prolog

%% Embedded ALGOL-like language in Prolog

change_member(_, N, [], [N]).
change_member(O, N, [O|XS], [N|XS]).
change_member(O, N, [X|XS], [X|YS]) :- change_member(O, N, XS, YS).


:- op(900, xfy, :=).
:- op(100, fx, [if, then, else]).
:- op(125, fx, [while,  return]).
:- op(125, yfx, do).


run({Block},                          EnvI, O, Ret) :-
  run(Block, EnvI, O, Ret).

run((First; Second),                  EnvI, O, Ret) :-
  run(First, EnvI, M, _), run(Second, M, O, Ret).

run(Place := Expr,                    EnvI, O, void) :-
  eval(Expr, Value, EnvI), change_member(Place = _, Place = Value, EnvI, O), !.

run(if(Cond) then {This} else {That}, EnvI, O, Ret) :-
  eval(Cond, Value, EnvI),
  ( Value = true -> Body = This ; Body = That ),
  run(Body, EnvI, O, Ret).

run(while(Cond) do {Body},            EnvI, O, Ret) :-
  eval(Cond, Value, EnvI),
  ( Value = true -> run(Body, EnvI, M, _),
                    run(while(Cond) do {Body}, M, O, Ret)
  ; O = EnvI, Ret = void
  ).

run(return Value,                     EnvI, O, Ret) :-
  eval(Value, Ret, EnvI), O = EnvI.


eval(Var,   Val,   Env) :- member(Var = Val, Env), !.
eval(Num,   Num,   _  ) :- number(Num), !.
eval(true,  true,  _  ).
eval(false, false, _  ).

eval(X + Y, Z, Env) :- eval(X, XV, Env), eval(Y, YV, Env), Z is XV + YV.
eval(X - Y, Z, Env) :- eval(X, XV, Env), eval(Y, YV, Env), Z is XV - YV.
eval(X * Y, Z, Env) :- eval(X, XV, Env), eval(Y, YV, Env), Z is XV * YV.
eval(X / Y, Z, Env) :- eval(X, XV, Env), eval(Y, YV, Env), Z is truncate(XV / YV).

eval(X = Y, B, Env) :- eval(X, XV, Env), eval(Y, YV, Env),
  ( XV = YV -> B = true ; B = false ).
eval(X < Y, B, Env) :- eval(X, XV, Env), eval(Y, YV, Env),
  ( XV < YV -> B = true ; B = false ).
eval(X > Y, B, Env) :- eval(X, XV, Env), eval(Y, YV, Env),
  ( XV > YV -> B = true ; B = false ).
eval(not(P), B, Env) :- eval(P, PV, Env),
  ( PV = true -> B = false ; B = true ).





factorial(N, {i := 1;
              n := N;
              while(n > 0) do {
                i := i * n;
                n := n - 1
              };
              return i}).

%% ?- factorial(5, P), run(P, [], _, X).

%% X = 120



fibonacci(F, {i := 1; a := 1; b := 1;
              while(not(i = F)) do {
                a := a + b;
                b := a - b;
                i := i + 1
              };
              return b}).

%% ?- fibonacci(7, P), run(P, [], _, X).

%% X = 13



pow2(N, {i := N;
         n := 1;
         while(not(i = 0)) do {
           i := i - 1;
           n := n * 2
         };
         return n}).

%% ?- pow2(8, P), run(P, [], _, X).

%% X = 256



gcd(X, Y,
    {x := X;
     y := Y;
     while(not(x = 0)) do {
       if(x < y) then {
         y := y - x
       }
       else {
         x := x - y
       }
     };
     return y}).

%% ?- A is 2*2*5*7*13, B is 5*13*7*7, gcd(A, B, P), run(P, [], _, X).

%% A = 1820,
%% B = 3185,
%% X = 455



isqrt(N, {q := N; p := (q + N / q) / 2;
          while(q > p) do {
            q := p;
            p := (q + N / q) / 2
          };
          return p}).

%% ?- isqrt(81, P), run(P, [], _, X).

%% X = 9

Type inference for The Simply Typed Lambda Calculus

%% Type inference for The Simply Typed Lambda Calculus

:- op(150, xfx, ⊢).
:- op(140, xfx, :).
:- op(100, xfy, ->).
:- op(100, yfx, $).

Γ ⊢    Term : Type          :- atom(Term), member(Term : Type, Γ).
Γ ⊢ λ(A, B) : Alpha -> Beta :- [A : Alpha|Γ] ⊢ B : Beta.
Γ ⊢   A $ B : Beta          :- Γ ⊢ A : Alpha -> Beta, Γ ⊢ B : Alpha.

/** Examples:

% Typing fix
?- Γ ⊢ y $ f : Y, Γ ⊢ f $ (y $ f) : Y.
Γ = [y: (Y->Y)->Y, f:Y->Y|_G330]

% Typing the Y combinator
?- Γ ⊢ λ(f, (λ(f,f $ f)) $ (λ(g,f $ (g $ g)))) : Y.
Y = (_G309 -> _G309) -> _G309 

% Typing flip id
?- Id = λ(i, i), Flip = λ(f, λ(y, λ(x, f $ x $ y))),
   Γ ⊢ Flip $ Id : T.
T = _G395 -> (_G395 -> _G411) -> _G411

**/

Generate, Test and Intertwine

%% Generate, Test and Intertwine
%% -----------------------------

%% The aim here is to generate 4x4 magic squares. First of all, We start
%% with the 3x3 case, since that's easier.
%% This way any improvements can be reused, assuming they are general enough

%% I'll define a magic square as, an nxn table of numbers such that all the
%% rows, columns and diagonals sum to the same value, each number in the range
%% 1 .. n^2 appears once in the grid.

%% In Prolog, We'll represent squares as lists:

%: [4,8,2,
%:  3,5,7,
%:  8,1,6]

%% So to generate this, We'll select numbers from a list of digits and insert
%% them into the table and check that everything adds up.

saturn_digits([1,2,3,4,5,6,7,8,9]).

saturn([A,B,C,
        D,E,F,
        G,H,I]) :-
        saturn_digits(D1),
        select(A, D1, D2),
        select(B, D2, D3),
        select(C, D3, D4),
        select(D, D4, D5),
        select(E, D5, D6),
        select(F, D6, D7),
        select(G, D7, D8),
        select(H, D8, D9),
        select(I, D9, []),
        S is A + B + C,
        S is D + E + F,
        S is G + H + I,
        S is A + D + G,
        S is B + E + H,
        S is C + F + I,
        S is A + E + I,
        S is C + E + G.

%% That first attempt works fine:

%? ?- saturn(Square).
%? Square = [2, 7, 6, 9, 5, 1, 4, 3, 8] ;
%? Square = [2, 9, 4, 7, 5, 3, 6, 1, 8] ;
%? Square = [4, 3, 8, 9, 5, 1, 2, 7, 6] ;
%? Square = [4, 9, 2, 3, 5, 7, 8, 1, 6] ;
%? Square = [6, 1, 8, 7, 5, 3, 2, 9, 4] ;
%? Square = [6, 7, 2, 1, 5, 9, 8, 3, 4] ;
%? Square = [8, 1, 6, 3, 5, 7, 4, 9, 2] ;
%? Square = [8, 3, 4, 1, 5, 9, 6, 7, 2] ;

%% It's quite clear the program is correct,
%% The two important properties are
%%  * Soundness -- Generates only magic squares
%%  * Completeness -- Generates every magic square
%%
%% Since digits are selected from the list, It's possible for the square to
%% be filled in every possible way, and due to:
%:   select(I, D9, []),
%% It's ensured that every digit is used. Every row, column and diagonal is
%% summed, since S is used for all summations, this ensures that all the sums
%% are equal.


%% There is this pattern in the code though, which I don't very much like:

%: select(A, D1, D2),
%: select(B, D2, D3),
%: select(C, D3, D4),
%: select(D, D4, D5),
%: select(E, D5, D6),
%: select(F, D6, D7),
%: select(G, D7, D8),
%: select(H, D8, D9),
%: select(I, D9, []),

%% Prolog has a macro, DCG style, that lets us erase this pattern.
%% If we use --> instead of :- , then this is implicit, everything is sewn
%% together like we did manually, except automatically.
%% {} is used to escape it.
%% So now, It can be rewritten to the equivalent program:

saturn_2([A,B,C,
          D,E,F,
          G,H,I]) -->
        select(A),
        select(B),
        select(C),
        select(D),
        select(E),
        select(F),
        select(G),
        select(H),
        select(I),
        { S is A + B + C },
        { S is D + E + F },
        { S is G + H + I },
        { S is A + D + G },
        { S is B + E + H },
        { S is C + F + I },
        { S is A + E + I },
        { S is C + E + G }.

%? ?- saturn_digits(D), saturn_2(Square, D, []).
%? D = [1, 2, 3, 4, 5, 6, 7, 8, 9],
%? Square = [2, 7, 6, 9, 5, 1, 4, 3, 8] ;
%? ..

%% It produces all the same answers in the same order. One of the fantastic
%% thing about Prolog is you could actually implement DCG in it, if it wasn't
%% already there, Or even your own macros.

%% Anyway, there is an optmisation which is quite blatant: select always
%% succeeds, It's the tests in {} which cause backtracking, So if we were to
%% intertwine the generating (done by select) with the testing, (as done by
%% summation) then backtracking which does occur, will happen sooner.

saturn_3([A,B,C,
          D,E,F,
          G,H,I]) -->
        select(A),
        select(B),
        select(C),
        { S is A + B + C },
        select(D),
        select(G),
        { S is A + D + G },
        select(E),
        { S is C + E + G },
        select(F),
        { S is D + E + F },
        select(H),
        { S is B + E + H },
        select(I),
        { S is G + H + I },
        { S is C + F + I },
        { S is A + E + I }.


%% So now, It's straightforward to use the same process for the 4x4 case.
%% We use DCG first in this case though.

jupiter_digits([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]).

jupiter([A,B,C,D,
         E,F,G,H,
         I,J,K,L,
         M,N,O,P]) -->
        select(A),
        select(B),
        select(C),
        select(D),
        select(E),
        select(F),
        select(G),
        select(H),
        select(I),
        select(J),
        select(K),
        select(L),
        select(M),
        select(N),
        select(O),
        select(P),
        
        { S is A + B + C + D },
        { S is E + F + G + H },
        { S is I + J + K + L },
        { S is M + N + O + P },
        
        { S is A + E + I + M },
        { S is B + F + J + N },
        { S is C + G + K + O },
        { S is D + H + L + P },
        
        { S is A + F + K + P },
        { S is D + G + J + M }.

%% And then rearrange

jupiter_2([A,B,C,D,
           E,F,G,H,
           I,J,K,L,
           M,N,O,P]) -->
        select(A),
        select(B),
        select(C),
        select(D),
        %% * * * *
        %% o o o o
        %% o o o o
        %% o o o o
        { S is A + B + C + D },
        
        select(E),
        select(I),
        select(M),
        %% * * * *
        %% * o o o
        %% * o o o
        %% * o o o
        { S is A + E + I + M },
        
        select(F),
        select(K),
        select(P),
        %% * * * *
        %% * * o o
        %% * o * o
        %% * o o *
        { S is A + F + K + P },
        
        select(G),
        select(H),
        %% * * * *
        %% * * * *
        %% * o * o
        %% * o o *
        { S is E + F + G + H },
        
        select(L),
        %% * * * *
        %% * * * *
        %% * o * *
        %% * o o *
        { S is D + H + L + P },
        
        select(J),
        %% * * * *
        %% * * * *
        %% * * * *
        %% * o o *
        { S is I + J + K + L },
        { S is D + G + J + M },
        
        select(N),
        %% * * * *
        %% * * * *
        %% * * * *
        %% * * o *
        { S is B + F + J + N },
        
        select(O),
        %% * * * *
        %% * * * *
        %% * * * *
        %% * * * *
        { S is M + N + O + P },
        { S is C + G + K + O }.

%? ?- jupiter_digits(D), jupiter_2(Square, D, []).
%? 
%? D = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16],
%? Square = [1, 2, 15, 16, 12, 14, 3, 5, 13, 7, 10, 4, 8, 11, 6, 9] ;
%? ...

%% There is a huge speed difference in the two versions.
%% Here's the final program without comments:

jupiter_3([A,B,C,D,
           E,F,G,H,
           I,J,K,L,
           M,N,O,P]) -->
        select(A),
        select(B),
        select(C),
        select(D),
        { S is A + B + C + D },
        select(E),
        select(I),
        select(M),
        { S is A + E + I + M },
        select(F),
        select(K),
        select(P),
        { S is A + F + K + P },
        select(G),
        select(H),
        { S is E + F + G + H },
        select(L),
        { S is D + H + L + P },
        select(J),
        { S is I + J + K + L },
        { S is D + G + J + M },
        select(N),
        { S is B + F + J + N },
        select(O),
        { S is M + N + O + P },
        { S is C + G + K + O }.