Toward ZDD library in SWI-Prolog

As always, interesting work… but maybe you alter somewhat in strange ways the module search path ? After loaded library(pac), I’m unable to edit any files with PCEMACS

?- [library(pac)].

?- edit(inv).
% Waiting for editor ... sh: 1: emacsclient: not found

% Editor returned failure; skipped make/0 to reload files

I am not familiar with pack-install, but, I have found that it installs the package at

Kuniaki Mukai

A moment before, I have made an experiment on ZDD in Prolog after Jon Burse.
Although I am not familiar with the problem and don’t know much about
consequences of solving the probolem. The predicate jburse/2 below is a
modificatio of Jan Burse’s original problem. I remember a query as for original
jburse(8, C) took more than 24 hours at that time, while the current query
jburse(11, C) takes, surprisingly, 0.16 seconds. Moreover even jburse(50, C)
takes less than 8 seconds. Unfortunately I don’t know the answers are correct,
but I take these statistics as a sympton that ZDD is powerful for representing
sets of clauses of literals for logical manipulation. I would like to see more
symptons for this, though I have not appropriate samples to try for the moment.

% ?- time(jburse(11, C)).
%@ % 214,911 inferences, 0.016 CPU in 0.016 seconds (100% CPU, 13089962 Lips)
%@ C = 36028797018963968.
% ?- time(jburse(12, C)).
%@ % 271,012 inferences, 0.021 CPU in 0.022 seconds (97% CPU, 12645203 Lips)
%@ C = 73786976294838206464.
% ?- time(jburse(50, C)).
%@ % 65,606,276 inferences, 7.624 CPU in 7.747 seconds (98% CPU, 8604897 Lips)
%@ C = 57775629806269263488289979632051482575456970143627537546788321689581
jburse(N, C):-  N2 is N*N,
		jburse_grounded(N, P),
		zdd(( boole_to_dnf(P, Z), ltr_card(N2, Z, C))).
jburse_grounded(N, P):- problem(N, P),
	term_variables(P, Vs),
	length(Vs, Len),
	Len0 is Len + 1,
	numlist(2, Len0, Ns),
	Vs = Ns.

Kuniaki Mukai

It was a faithful translation of your problem, I thought.
Anyway I thank you, it was more than a couple years ago, and
I will not go into further.

?- listing(problem).

should show listings such like this.

% by Jan Burse. For testing %
problem(N, P):- matrix(N, M), constraint(M, P).
burse(N, P):- matrix(N, M), constraint(M, P).
matrix(N, M):- length(M, N), maplist(col(N), M).
col(N, L):- length(L, N).
inner_product([X], [Y], X*Y):-!.
inner_product([X|Xs], [Y|Ys], #((X*Y),Z)):- inner_product(Xs, Ys, Z).
constraint([V], E=:=true):- !, inner_product(V, V, E).
constraint([V|Vs], (E=:=true)*C*D):- inner_product(V, V, E),
constraint(Vs, D),
inner_constraint_one(V, Vs, C).
inner_constraint_one(V, [U], C=:=false):- !, inner_product(V, U, C).
inner_constraint_one(V, [U|Us], (C=:=false)*D):-
inner_product(V, U, C),
inner_constraint_one(V, Us, D).

I am sorry I am already too busy for old problems to spend time
for new proposed ones.

Kuniaki Mukai

Is it too fast ? I partially agree. Also I don’t remember
answers retuned by my old version. I lost old codes. Are answers below different
from yours for n =< 8 ? There might be bug(s) in my new
codes. I have to take time to review it,

% ?- forall(between(1, 10, N),
% ( jburse(N, C),
% writeln(jburse(N, C))
% )).
%@ jburse(1,1)
%@ jburse(2,2)
%@ jburse(3,8)
%@ jburse(4,64)
%@ jburse(5,1024)
%@ jburse(6,32768)
%@ jburse(7,2097152)
%@ jburse(8,268435456)
%@ jburse(9,68719476736)
%@ jburse(10,35184372088832)
%@ true.

Kuniaki Mukai

I found a serious logic bug in the ltr_card (counting solutions for DNF),
and have fixed it. Then the following query again. The result is still surprising.
I have no idea how to verify this result. As once I had hard time on this problem
to improve time efficiency, I can not believe the drastic current speed. However it is possible
that the codes is correct in some unexpected way.
Anyway I take time for a while for cooling down.

% ?- time(forall(between(1, 20, N),
%			( jburse(N, C),
%			  writeln(jburse(N, C))
%			))).
%@ jburse(1,1)
%@ jburse(2,2)
%@ jburse(3,8)
%@ jburse(4,64)
%@ jburse(5,1024)
%@ jburse(6,32768)
%@ jburse(7,2097152)
%@ jburse(8,268435456)
%@ jburse(9,68719476736)
%@ jburse(10,35184372088832)
%@ jburse(11,36028797018963968)
%@ jburse(12,73786976294838206464)
%@ jburse(13,302231454903657293676544)
%@ jburse(14,2475880078570760549798248448)
%@ jburse(15,40564819207303340847894502572032)
%@ jburse(16,1329227995784915872903807060280344576)
%@ jburse(17,87112285931760246646623899502532662132736)
%@ jburse(18,11417981541647679048466287755595961091061972992)
%@ jburse(19,2993155353253689176481146537402947624255349848014848)
%@ jburse(20,1569275433846670190958947355801916604025588861116008628224)
%@ % 9,061,329 inferences, 0.674 CPU in 0.727 seconds (93% CPU, 13436672 Lips)
%@ true.

Kuniaki Mukai

1 Like

Didn’t you provide a pack before? Should be possible to retrieve that, no?

Thank you for your kind proposal, but
I saw a fragment of log in my local git like this:

?- time(continue((problem(7, P), solc(P, X)))).
%@ % 366,853,945 inferences, 94.389 CPU in 96.556 seconds (98% CPU, 3886607 Lips)
%@ X = 1451520.

I am still not good at GIT. I wish I could grep on all commits in my local repository, for example.

On the other hand, this is a log by the new version with
boole_to_dnf + ltr_card, where “everything” is done simply in ZDD.

% ?- time(jburse(7, C)).
38,378 inferences, 0.004 CPU in 0.004 seconds (84% CPU, 10269735 Lips)
%@ C = 2097152.

The difference suggests that there are bugs in even both of
my two codes. However, as I don’t know the correct answer yet, and
the new version works correctly so far for all test cases
including on propostional logic which I prepared, I will take care of
the new codes only. I will appreciate bug reports.

Kuniaki Mukai

I have found that what the predicate jburse produces is reproduced
by another simple way below. It is interesting that the original problem
by Jan Burse has this property. I will be happy if it is the case because
it supports correctness of my ZDD toward a SWI-Prolog package.

% ?- forall( between(1, 20, I),
%			(	f(I, P),
%				X is 2^P,
%				writeln((2^f(I)=2^P)=X)
%			)).

%@ (2^f(1)=2^1)=2
%@ (2^f(2)=2^3)=8
%@ (2^f(3)=2^6)=64
%@ (2^f(4)=2^10)=1024
%@ (2^f(5)=2^15)=32768
%@ (2^f(6)=2^21)=2097152
%@ (2^f(7)=2^28)=268435456
%@ (2^f(8)=2^36)=68719476736
%@ (2^f(9)=2^45)=35184372088832
%@ (2^f(10)=2^55)=36028797018963968
%@ (2^f(11)=2^66)=73786976294838206464
%@ (2^f(12)=2^78)=302231454903657293676544
%@ (2^f(13)=2^91)=2475880078570760549798248448
%@ (2^f(14)=2^105)=40564819207303340847894502572032
%@ (2^f(15)=2^120)=1329227995784915872903807060280344576
%@ (2^f(16)=2^136)=87112285931760246646623899502532662132736
%@ (2^f(17)=2^153)=11417981541647679048466287755595961091061972992
%@ (2^f(18)=2^171)=2993155353253689176481146537402947624255349848014848
%@ (2^f(19)=2^190)=1569275433846670190958947355801916604025588861116008628224
%@ (2^f(20)=2^210)=1645504557321206042154969182557350504982735865633579863348609024
%@ true.

f(1) = 1.
f(n) = f(n-1) + n.

f(1, 1):-!.
f(N, X):- N>1,
	N0 is N-1,
	f(N0, X0),
	X is X0 + N.

Kuniaki Mukai

Here is an expereiment on ZDD, which might be interesting for the discourse.
It is easy to build such a (big) set F of clauses of literals in ZDD

  1. F has N varialbles. (N: given).
  2. F is valid as DNF.
  3. F is unsatisifiable as CNF.

By 2 the number of solutions of F should be 2^N. By 3 a refutation for F should succeed.

In fact, as log below, predicate ltr_refute below for F succeeds, and ltr_card for F returns 2^N.

Time statisics sounds to say: counting is much faster than refutation. Of course my current codes of ltr_refute is so simple and naive as the source codes listing below. Sadly, I have no idea how to improve codes. Anyway although ZDD implementation of mine is based on my simple-minded background on ZDD, I have already shown up, I hope, that some advertised powerful features of ZDD is possible in SWI-Prolog in small prolog codes.

The current version of the pack pac is 1.5.6. After pack_install the pac, the following log shoud be reproduced:

% ?- use_modul(library(pac)).
% ?- use_modul(zdd(zdd)).
% ?- module(zdd).
% ?- time(zdd((a_big_cnf(100, X), ltr_card(100, X, C)))), C =:= 2^100.
%@ % 2,392,546 inferences, 0.290 CPU in 0.409 seconds (71% CPU, 8246121 Lips)
%@ X = 30001,
%@ C = 1267650600228229401496703205376 .
% ?- call_with_time_limit(100, time(zdd((a_big_cnf(30, X), ltr_refute(X))))).
%@ % 123,551,444 inferences, 14.271 CPU in 14.439 seconds (99% CPU, 8657765 Lips)
%@ X = 2701.
a_big_cnf(X, Y):- shift(a_big_cnf(X, Y)).
a_big_cnf(0, 1, _).
a_big_cnf(N, X, S):- N>0,
	N0 is N-1,
	a_big_cnf(N0, X0, S),
	ltr_insert(N, X0, R, S),
	ltr_insert(-N, X0, L, S),
	zdd_join(L, R, X, S).
% refutation by resolution in ZDD.
ltr_refute(X):- shift(ltr_refute(X)).
ltr_refute(X, S) :- zdd_has_empty(X, S), !.
ltr_refute(X, S) :- X > 1,
	cofact(X, t(A, L, R), S),
	(	L = 0  -> Y = R
	;	cofact(L, t(B, L0, R0), S),
		(	B = -(A)
		-> 	ltr_merge(R0, R, Y0, S),
			ltr_join(L0, Y0, Y, S)
		;   ltr_join(L, R, Y, S)
	ltr_refute(Y, S).

family_card/2 below classifies a given family of sets by cardiality.
The log below is for the powerset of a set U = [1,2,3,…,1000],
and shows the number of subsets of U which has 500 elements.
I think that time 10.355 CPU in 10.540 seconds is reasonable
as for ZDD in SWI-Prolog. The codes below contains reduntant in my sence
because I see multiple time term_hash calll in a redundant way, but
I am almost new for this subttle thing on programming, and busy enough
only for naive version.

% ?- N = 1000, numlist(1, N, Ns),
%	time(zdd((X<<pow(Ns), family_card(X, P),
%	zmemo(family_with_card(500)-L), card(L, C)))).

%@ % 90,174,772 inferences, 10.355 CPU in 10.540 seconds (98% CPU, 8708585 Lips)
%@ C = 270288240945436569515614693625975275496152008446548287007392875106625428705522193898612483924502370165362606085021546104802209750050679917549894219699518475423665484263751733356162464079737887344364574161119497604571044985756287880514600994219426752366915856603136862602484428109296905863799821216320.
% Classify sets in a family by cardinality.

max(A, B, A):- A@>B, !.
max(_, B, B).

family_card(X, Y):- shift(family_card(X, Y)).
family_card(0, 0, _):-!.
family_card(1, 0, S):-!, zmemo(family_with_card(0)-1, S).
family_card(I, P, S):- zmemo(family_card_done(I)-B, S),
	(	nonvar(B) ->  true
	; 	cofact(I, t(A, L, R), S),
		family_card(L, Pl, S),
		family_card(R, Pr, S),
		max(Pl, Pr, P0),
		insert_one_to_family(A, P0, New, S),
		P is P0 + 1,
		zmemo(family_with_card(P)-New, S),
		B = true
insert_one_to_family(A, 0, G, S):-!, zmemo(family_with_card(0)-F, S),
	insert_atom(A, F, G, S).
insert_one_to_family(A, P, G, S):- P0 is P-1,
	insert_one_to_family(A, P0, G0, S),
	zmemo(family_with_card(P)-Fp, S),
	insert_atom(A, Fp, G, S),
	zdd_join(Fp, G0, Gp, S),
	set_memo(family_with_card(P)-Gp, S).

Kuniaki Mukai

Thank you advice for checking the Kochen-Specker paradox paper. It is impressive that
your problem on orthogonal vectors is related to the paper, though I know nothing about the Kochen-Specker paradox. I will check it later.

Unfortunately, I am busy for a couple of programming on ZDD. One is on path counting of a grid graph ( Minato group solved for grid nodes {(x,y) | 1=<x,y=<N } N=< 21 ?? about 10 years ago based on SimPath method by D.Knuth. I know little on this.) The other is related counting a single sorted associative algebra with a given cardinality. On both of them I have only miserable failure history in my life as challage for “Prolog programming” in my sense. Still I can’t expect much about my ZDD implementation, but anyway I would like to post something here on experiment results soon.

Kuniaki Mukai

About the problem of counting associative algebras on a domain D, I have run two versions of codes to get the following results. One, say A, is using ZDD, the other, say B, is not. For |D|=1, 2, 3, A and B returns 1, 8, 113, respectively. For |D|=4, A returns 3492, but B got “timeout.” As the number of algebras on D are 1, 2^4, 3^9, 4^{16}, 5^{25}, \ldots for |D|= 1, 2, 3, 4, 5, \ldots, respectively, the current results should be acceptable, I think. I naively expected on beginning the experiment that there might be some idea on ZDD to use its property of rich structure sharing, but I could not find them. Thus, both codes are naive and brute force. I will appreciate for any pointer to this counting problem. I’m amateur on counting problem in general.

Query log.
% ?-zdd(( count_associative_algebra_with_zdd(1, C))).
%@ C = 1.

% ?-zdd(( count_associative_algebra_with_zdd(2, C))).
%@ C = 8.

% ?-time(zdd(( count_associative_algebra_with_zdd(3, C)))).
%@ % 426,394 inferences, 0.074 CPU in 0.074 seconds (100% CPU, 5789228 Lips)
%@ C = 113.

% ?-time(zdd(( count_associative_algebra_with_zdd(4, C)))).
%@ % 143,325,850 inferences, 25.723 CPU in 25.729 seconds (100% CPU, 5571848 Lips)
%@ C = 3492.

% ?-call_with_time_limit(1800, zdd(( count_associative_algebra_with_zdd(5, C)))).
%@ ERROR: Unhandled exception: Time limit exceeded
Codes using ZDD.
count_associative_algebra_with_zdd(N, C):-
	prepare_count_associative_algebra(N, Track, Etree),
	count_assoc_algebras([h(Etree, [], Track)], 0, C).
prepare_count_associative_algebra(N, Track, Etree):-
	numlist(1, N, Dom),
	raise_list(Dom, 2, Dom2),
	product_list(Dom, Dom2, Dom3),
	maplist(pred([U, a(U, [], [])]), Dom3, Track),
	zdd:map_space(Dom2, Dom, Etree, equational_mapsto).

count_assoc_algebras(X, Y, Z):- shift0(count_assoc_algebras(X, Y, Z)).

count_assoc_algebras([], C, C, _).
count_assoc_algebras([h(0, _, _)|Hs], C, C0, S):-!,
	count_assoc_algebras(Hs, C, C0, S).
count_assoc_algebras([h(_, _, 0)|Hs], C, C0, S):-!,
	count_assoc_algebras(Hs, C, C0, S).
count_assoc_algebras([h(1, _, [])|Hs], C, C0, S):-!, C1 is C+1,
	count_assoc_algebras(Hs, C1, C0, S).
count_assoc_algebras([h(1, _, _)|Hs], C, C0, S):-!,
	count_assoc_algebras(Hs, C, C0, S).
count_assoc_algebras([h(Etree, Es, Tr)|Hs], C, C0, S):-
	cofact(Etree, t(E, L, R), S),
	repeat_law([E|Es], Tr, Tr0),
	(	Tr0 = 0
	->	count_assoc_algebras([h(L, Es, Tr)|Hs], C, C0, S)
	;	count_assoc_algebras([h(L, Es, Tr), h(R, [E|Es], Tr0)|Hs], C, C0, S)

I tried ordinal codes without using ZDD.

Query without using ZDD.
% ?- count_associative_algebra(1, X).
%@ X = 1.

% ?- count_associative_algebra(2, X).
%@ X = 8.

% ?- time(count_associative_algebra(3, X)).
%@ % 87,217,314 inferences, 15.452 CPU in 15.477 seconds (100% CPU, 5644233 Lips)
%@ X = 113.

% ?- call_with_time_limit(180, time(count_associative_algebra(4, X))).
%@ % 984,263,663 inferences, 157.935 CPU in 198.803 seconds (79% CPU, 6232073 Lips)
%@ ERROR: Unhandled exception: Time limit exceeded
Codes without using ZDD.
count_associative_algebra(N, NumOfAssocAlgebras):-
	numlist(1, N, Dom),
	raise_list(Dom, 2, Dom2),
	product_list(Dom, Dom2, Dom3),
	maplist(pred([U, a(U, [], [])]), Dom3, Track0),
	sort(Track0, Track1),
	map_space_in_list(Dom2, Dom, Maps),
	maplist(maplist( pred([[A,B]-C, A+B=C]) ), Maps, FamEs),
	foldl( pred(Track1, (	[Es, C, C0]	:-
							repeat_law(Es, Track1, Track),
							(		Track = [] -> C0 is C+1
							;		C0 = C
		   FamEs,  0,  NumOfAssocAlgebras).

% ?- N is 3, M is 3, time(( numlist(1, N, X), numlist(1, M, Y),
%	map_space_in_list(X, Y, F), length(F, C))).
map_space_in_list([], _, [[]]).
map_space_in_list([X|Xs], Y, Z):- map_space_in_list(Xs, Y, Z0),
	map_space_in_list(Y, X, Z0, Z, []).
map_space_in_list([], _, _, Z, Z).
map_space_in_list([Y|Ys], X, Z, U, V):-
	fold_map_extend(Z, Y, X, U, U0),
	map_space_in_list(Ys, X, Z, U0, V).
fold_map_extend([], _, _, U, U).
fold_map_extend([Z|Zs], Y, X, [[X-Y|Z]|U], V):-
	fold_map_extend(Zs, Y, X, U, V).

Common predicates for checking assoiciative law.
repeat_law(E, [a(_,[A],[A])|X], Y):-!, repeat_law(E, X, Y).
repeat_law(_, [a(_,[_],[_])|_], 0):-!.	% 0 means NON Associativitty.
repeat_law(Eqs, X, X0):- member(E, Eqs),
	apply_law_once(E, X, X1),
	repeat_law(Eqs, X1, X0).
repeat_law(_, X, X).

apply_law_once(E, X, [A0|X0]):-
	select(A, X, X0),
	law(E, A, A0).

% law/3
law( A+B=C,	a([X, A, B], U, []),
			a([X, A, B], U, [X, C])).
law( A+B=C,	a(T, U, [A, B]),
			a(T, U, [C])).
law( A+B=C,	a(T, [A, B], U),
			a(T, [C], U)).
law( A+B=C,	a([A, B, X], [], U),
			a([A, B, X], [C, X], U)).

Kuniaki Mukai

3^9 is the number of algebras on \{1,2,3\}. 113 is the number of associative algebras on {1, 2, 3}.

Kuniaki Mukai

An algebra, here, is taken as a system of equations a+b = c which is viewed as a term rewriting system
The " Common predicates for checking associative law" is just to check confluency of the rewriting system,
In fact, codes is after what I was taught from Fernando Pereira in the middle of 1980’s. It is still one of my favorite programming techniques, among DCG, reverse/3 with stack, etc. Anyway my codes is naive, and brute force except using ZDD. I think if you adopt ZDD as it is, surely your will be able to to go at least beyond my n=4. FOL and model checking are out of my sight for the time being. My interest is whether or not ZDD is good to support for SWI-Prolog package. I need more programming experience on ZDD.

Kuniaki Mukai

Amazing ! Thank you for the pointer.


Let F be the set of algebras on D:

F=\{f\mid f\colon D\times D\rightarrow D\}

Then the set A of associative algebras on D

A=\{f\in F\mid f(x, f(y, z))=f(f(x,y), z)\quad \forall x, y, z\in D\}.

To get |A|, my codes builds faithfully D^3, then construct F as a family of subsets D^3, where tht law predicate in my codes plays a role for the condition part of A. Anyway my codes is written straightforwardly along this line. Thus I said the codes is naive and brute force.
In fact, I am afraid |D|=4 may be the max case which this kind of brute force can reach. Thanks to the pointer you gave I take it a good symptom for ZDD toward SWI-Prolog package, though I have not yet read the amazing article. Thanks, hoping this comment would be useful for you if you read the cryptic codes of mine.


I am on the way of restructuring my codes around ZDD in the pack pac little by little, which still looks like a huge heap of junk codes. However, the recent core of ZDD in that junk codes is short, neat and powerful, I hope, compared with that of previous versions. Unfortunately I’m old enough to do restructuring codes in a fast and modern way. Prolog programming has been always a hobby of mine. With this situation, I am afraid that I am not able to respond quickly to any request.


product_list is for Cartesian product X\times Y. raise_list is for X^n ( n times cartesian product).

Query and codes.
% ?- N=50,  numlist(1, N, Ns), time(raise_list(Ns, 5,X)), length(X, L).
%@ % 318,878,060 inferences, 24.608 CPU in 28.858 seconds (85% CPU, 12958082 Lips)

% ?- N=100,  numlist(1, N, Ns), time(raise_list(Ns, 3,X)), length(X, L).
%@ % 1,010,706 inferences, 0.102 CPU in 0.121 seconds (85% CPU, 9880307 Lips)

raise_list(_, 0, [[]]):-!.
raise_list(L, N, P):- N0 is N-1,
	raise_list(L, N0, P0),
	product_cons(L, P0, P, []).

product_list(X, Y, P) :- product_cons(X, Y, P, []).

product_cons([], _, P, P).
product_cons([A|As], L, P, P0):-
	product_cons_aux(L, A, P, Q),
	product_cons(As, L, Q, P0).
product_cons_aux([], _, P, P).
product_cons_aux([X|Xs], A, [[A|X]|P], Q):-
	product_cons_aux(Xs, A, P, Q).


Why not restrict symmetric algebras f f(x,y)=f(y,x) ? I will check my codes for this later. Anyway, I wish to read the article which you have pointed to us. I love also mathematical solution as well as naive and brute force one. I guess the article must be mathematical one.


Sounds great !