Search code examples
prologpuzzleclpfd

Solving a 4x4 multiplicative puzzle "5040" in Prolog with clp(FD)


Today, I found a puzzle at https://puzzling.stackexchange.com/questions/22064/the-5040-square:

Fill a 4x4 grid with positive integers so that:

Every cell has a different integer

The product of the numbers in each row is 5040, and similarly for the columns

Source: This was an NPR weekly listener challenge, aired on 2005-10-09

Here's my first shot at solving the puzzle using :

:- use_module(library(clpfd)).

m5040_(Mss,Zs) :-
   Mss = [[A1,A2,A3,A4],
          [B1,B2,B3,B4],
          [C1,C2,C3,C4],
          [D1,D2,D3,D4]],
   Zs = [A1,A2,A3,A4,B1,B2,B3,B4,C1,C2,C3,C4,D1,D2,D3,D4],
   Zs ins 1..sup,           % domain: positive integers
   5040 #= A1*A2*A3*A4,     % rows
   5040 #= B1*B2*B3*B4,
   5040 #= C1*C2*C3*C4,
   5040 #= D1*D2*D3*D4,
   5040 #= A1*B1*C1*D1,     % columns
   5040 #= A2*B2*C2*D2,
   5040 #= A3*B3*C3*D3,
   5040 #= A4*B4*C4*D4,
   all_different(Zs).       % pairwise inequality

Sample query:

?- m5040_(Mss,Zs), time(labeling([],Zs)).
% 416,719,535 inferences, 55.470 CPU in 55.441 seconds (100% CPU, 7512588 Lips)
Mss = [[1,3,16,105],[10,14,9,4],[21,8,5,6],[24,15,7,2]], Zs = [1,3,16,105|...] ;
...

My actual question is twofold:

  1. How can I speed up the backtracking search process for one / for all solutions?
  2. Which symmetries / redundancies could I exploit?

Solution

  • my inference counters (using your code) don't match yours... not sure why... the first solution i get (with your code)

    ?- puzzle_5040.
    % 464,043,891 inferences, 158.437 CPU in 160.191 seconds (99% CPU, 2928894 Lips)
    [[1,3,16,105],[10,14,9,4],[21,8,5,6],[24,15,7,2]]
    true 
    

    I thought that reducing the domain could speedup the result

    :- use_module(library(clpfd)).
    :- use_module(library(ordsets)).
    :- use_module(library(apply)).
    
    m5040_(Mss,Zs) :-
        matrix(Mss),
        flatten(Mss, Zs),
        all_factors(Fs),
        make_domain(Fs, Dom),
        Zs ins Dom,
        all_distinct(Zs),
        maplist(m5040, Mss),
        transpose(Mss, Tss), maplist(m5040, Tss).
    
    m5040([A,B,C,D]) :- 5040 #= A * B * C * D.
    length_(L, Xs) :- length(Xs, L).
    
    matrix(Mss) :-
        length_(4, Mss),
        maplist(length_(4), Mss).
    
    factors(L) :-
        L = [A,B,C,D],
        5040 #= 1 * 2 * 3 * U,
        L ins 1..U,
        all_distinct(L),
        A #< B, B #< C, C #< D,
        5040 #= A * B * C * D.
    
    all_factors(AllFs) :-
        findall(L, (factors(L),label(L)), Fs),
        foldl(ord_union, Fs, [], AllFs).
    

    but I was wrong, it was slower actually...

    Since some time ago I tried CLP(FD) solving some Project Euler, and in some cases I found it was slower than raw arithmetic, I arranged a program that doesn't use CLP(FD), but reduces the domain to make it manageable:

    puzzle_5040_no_clp :- time(puzzle_5040_no_clp(S)), writeln(S).
    
    puzzle_5040_no_clp(S) :-
        findall(F, factors(F), Fs),
        factors_group(Fs, G),
        once(solution(G, S)).
    
    disjoint(A, B) :-
        forall(member(X, A), \+ memberchk(X, B)).
    
    eq5040([A,B,C,D]) :-
        5040 =:= A * B * C * D.
    
    factors([A, B, C, D]) :-
        5040 #= 1 * 2 * 3 * U,
        [A, B, C, D] ins 1..U,
        A #< B, B #< C, C #< D,
        5040 #= A * B * C * D,
        label([A, B, C, D]).
    
    all_factors(AllFs) :-   % no more used
        findall(L, factors(L), Fs),
        foldl(ord_union, Fs, [], AllFs).
    
    factors_group(Fs, [A, B, C, D]) :-
        nth1(Ap, Fs, A),
        nth1(Bp, Fs, B), Ap < Bp, disjoint(A, B),
        nth1(Cp, Fs, C), Bp < Cp, disjoint(A, C), disjoint(B, C),
        nth1(Dp, Fs, D), Cp < Dp, disjoint(A, D), disjoint(B, D), disjoint(C, D).
    
    /*
    solution([A,B,C,D], S) :-
        maplist(permutation, [B,C,D], [U,V,Z]),
        transpose([A,U,V,Z], S),
        maplist(eq5040, S).
    */
    solution(T0, [U,V,X,Y]) :-
        peek5040(T0, U, T1),
        peek5040(T1, V, T2),
        peek5040(T2, X, T3),
        peek5040(T3, Y, [[],[],[],[]]).
    
    peek5040([A,B,C,D], [M,N,P,Q], [Ar,Br,Cr,Dr]) :-
        select(M,A,Ar),
        select(N,B,Br), M*N < 5040,
        select(P,C,Cr), M*N*P < 5040,
        select(Q,D,Dr), M*N*P*Q =:= 5040.
    
    % only test
    validate(G) :- maplist(eq5040, G), transpose(G, T), maplist(eq5040, T).
    

    with this approach, getting all solutions

    ?- time(aggregate(count,puzzle_5040_no_clp,N)).
    % 6,067,939 inferences, 1.992 CPU in 1.994 seconds (100% CPU, 3046002 Lips)
    [[1,24,14,15],[3,21,10,8],[16,5,9,7],[105,2,4,6]]
    % 111,942 inferences, 0.041 CPU in 0.052 seconds (79% CPU, 2758953 Lips)
    [[1,24,10,21],[3,15,14,8],[16,7,9,5],[105,2,4,6]]
    ...
    % 62,564 inferences, 0.033 CPU in 0.047 seconds (70% CPU, 1894080 Lips)
    [[1,10,12,42],[15,28,3,4],[16,9,7,5],[21,2,20,6]]
    % 37,323 inferences, 0.017 CPU in 0.027 seconds (65% CPU, 2164774 Lips)
    [[1,14,12,30],[15,2,28,6],[16,9,5,7],[21,20,3,4]]
    % 2,281,755 inferences, 0.710 CPU in 0.720 seconds (99% CPU, 3211625 Lips)
    % 48,329,065 inferences, 18.072 CPU in 27.535 seconds (66% CPU, 2674185 Lips)
    N = 354.