Picat中的分区函数P

问题描述

我得到了Partition函数P的如下实现
在 Prolog 中,从 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),M is K*(3*K+1)//2,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.

如何在 Picat 中实现相同的功能?是吗
真的,aggregate_all+sum 可以替换为 foreach+:= ?
Picat 中的 bignum 怎么样?

解决方法

Bignums 在 Picat 中没有问题。 这是我的 Picat 版本(受 Maple 方法启发):

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,P = S.

你的(整洁的)SWI-Prolog 版本在我的机器上 p(6666) 需要大约 1.9 秒。

?- 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.

我的 Picat 版本大约需要 0.2 秒

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

CPU time 0.206 seconds.

更新 这是在 Picat 中使用 findall 的版本,有点模仿您的方法:

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,M is K*(3*K+1)//2,X is -A-B.

但速度要慢得多(2.6s vs 0.2s):

p(6666) = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863

CPU time 2.619 seconds. Backtracks: 0

我还测试了实现相同的方法,即在 SWI-Prolog 中使用 findall/3

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

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

它比 Picat 的 findall 方法更快,并且与您的版本时间大致相同(稍快但推理更多)。

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

我尝试将自动 Picat 样式转换为一些更高阶的循环结构。然后手动内联高阶循环结构。自动 Picat 风格翻译的输入是:

:- table p/2.
p(0)=1 => !.
p(N)=Z =>
    Z=0,K=1,while(M=<N,(Z:=Z-(-1)^K*p(N-M),K:=K+1,M:=K*(3*K-1)//2)),K:=1,M:=K*(3*K+1)//2,M:=K*(3*K+1)//2)).

在此答案的末尾找到了翻译器代码的链接。自动翻译器给了我:

?- listing(p_a/2).
% example2.pl

:- sys_notrace p_a/2.
p_a(0,1) :-
   !.
p_a(N,A) :-
   Z = 0,while([B,C,D,E]\[F,G,H]\(G is D- -1^E*p(C-B),H is E+1,F is H*(3*H-1)//2),[I,J,_,_]\(I =< J),[M,Z,K],[_,L,_]),O is 1,P is O*(3*O+1)//2,while([Q,R,S,T]\[U,V,W]\(V is S- -1^T*p(R-Q),W is T+1,U is W*(3*W+1)//2),[X,Y,_]\(X =< Y),[P,O],A,_]).

它使用算术函数求值 p(C-B) 和 p(R-Q)。在我的 Prolog 系统算术函数求值使用原生 Java 堆栈,我无法求值 6666:

% ?- p(100,X).
% X = 190569292

% ?- p(6666,X).
% java.lang.StackOverflowError
%   at jekpro.reference.arithmetic.EvaluableElem.moniEvaluate(EvaluableElem.java:207)

同时使用 while/4 元谓词也有点慢。所以我按摩了代码,消除了算术函数评估并内联了 while/4。我也用过穷人表,速度快一点:

:- thread_local p_cache/2.

p_manual(N,X) :- p_cache(N,X),!.
p_manual(0,1) :-
   !,assertz(p_cache(0,1)).
p_manual(N,p21([M,p22([P,assertz(p_cache(N,A)).

p21([B,E],O1) :- B =< C,L is C-B,p_manual(L,M),G is D- -1^E*M,F is H*(3*H-1)//2,p21([F,H],O1).
p21(I1,I1).

p22([Q,T],O2) :- Q =< R,L is R-Q,V is S- -1^T*M,U is W*(3*W+1)//2,p22([U,W],O2).
p22(I2,I2).

现在事情开始看起来不错。 2.743 秒下降到:

/* SWI-Prolog 8.3.21 */
?- time(p_manual(6666,X)).
% 4,155,198 inferences,0.879 CPU in 0.896 seconds (98% CPU,4729254 Lips)
X = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863.

/* Jekejeke Prolog 1.5.0 */
?- time(p_manual(6666,X)).
% Up 736 ms,GC 20 ms,Threads 714 ms (Current 04/14/21 02:16:45)
X = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863

开源:

Picat 风格脚本翻译器 II
https://gist.github.com/jburse/8a24fe5668960c8889770f40c65cdf35#file-picat2-pl

Picat 风格脚本示例 II
https://gist.github.com/jburse/8a24fe5668960c8889770f40c65cdf35#file-example2-pl

Picat 样式脚本内联
https://gist.github.com/jburse/8a24fe5668960c8889770f40c65cdf35#file-tune-pl