Partition function P in Picat

Question

I got the following realization of the Partition function P
in Prolog, took it from rosetta here:

/* SWI-Prolog 8.3.21 */
:- table p/2.
p(0, 1) :- !.
p(N, X) :-
    aggregate_all(sum(Z), (between(1,inf,K), M is K*(3*K-1)//2,
           (M>N, !, fail; L is N-M, p(L,Y), Z is (-1)^K*Y)), A),
    aggregate_all(sum(Z), (between(1,inf,K), M is K*(3*K+1)//2,
           (M>N, !, fail; L is N-M, p(L,Y), Z is (-1)^K*Y)), B),
    X is -A-B.
 
?- time(p(6666,X)).
% 13,962,294 inferences, 2.610 CPU in 2.743 seconds (95% CPU, 5350059 Lips)
X = 1936553061617076610800050733944860919984809503384
05932486880600467114423441282418165863.

How would one go about an implementation of the same in Picat? Is it
true that aggregate_all+sum can be replaced by foreach+:= ?
What about bignums in Picat?

Answer 1

Bignums is no problem in Picat. Here's my Picat version (inspired by the Maple approach):

table
partition1(0) = 1.
partition1(N) = P =>
  S = 0,
  K = 1,
  M = (K*(3*K-1)) // 2,  
  while (M <= N)
     M := (K*(3*K-1)) // 2,  
     S := S - ((-1)**K)*partition1(N-M),
     K := K + 1
  end,
  K := 1,
  M := (K*(3*K+1)) // 2,
  while (M <= N)
     M := (K*(3*K+1)) // 2,  
     S := S - ((-1)**K)*partition1(N-M),
     K := K + 1
  end,
  P = S.

Your (neat) SWI-Prolog version takes about 1,9s for p(6666) on my machine.

?- time(p(6666,X)), write(X), nl.
% 13,959,720 inferences, 1.916 CPU in 1.916 seconds (100% CPU, 7285567 Lips)
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
X = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863.

My Picat version takes about 0.2s

Picat> time(println('p(6666)'=partition1(6666))) 
p(6666) = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863

CPU time 0.206 seconds.

Update Here's a version using findall in Picat instead, sort of mimicking your approach:

table
p(0, 1) :- !.
p(N, X) :-
    A = sum(findall(Z, (between(1,N,K), M is K*(3*K-1)//2,
           (M>N, !, fail; p(N-M,Y), Z is (-1)**K*Y)))),
    B = sum(findall(Z, (between(1,N,K), M is K*(3*K+1)//2,
           (M>N, !, fail; p(N-M,Y), Z is (-1)**K*Y)))),
    X is -A-B.

But it's much slower (2.6s vs 0.2s):

p(6666) = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863

CPU time 2.619 seconds. Backtracks: 0

I also tested to implement the same approach, i.e. using findall/3 in SWI-Prolog:

:- table p2/2.
p2(0, 1) :- !.
p2(N, X) :-
        findall(Z, (between(1,N,K), M is K*(3*K-1)//2,
                    (M>N, !, fail; L is N-M, p2(L,Y), Z is (-1)**K*Y)), AA),
        sum(AA,A),
        findall(Z, (between(1,N,K), M is K*(3*K+1)//2,
                    (M>N, !, fail; L is N-M, p2(L,Y), Z is (-1)**K*Y)), BB),
        sum(BB,B),
        X is - A - B.

sum(L,Sum) :-
        sum(L,0,Sum).
sum([],Sum,Sum).
sum([H|T],Sum0,Sum) :-
        Sum1 is Sum0 + H,
        sum(T,Sum1,Sum).

It's faster than Picat's findall approach, and about the same time as your version (slightly faster but with more inferences).

?- time(p2(6666,X)).
% 14,636,851 inferences, 1.814 CPU in 1.814 seconds (100% CPU, 8070412 Lips)
X = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863.