1:- module(frtvec, [udg_path_count/3,
    2				   udg_path_count/4,
    3				   rect_path_count/4,
    4				   rect_links/2,
    5				   prepare_udg/1,
    6				   udg_path/2,
    7				   hamilton_filter/3
    8				  ]).    9
   10:- use_module(zdd('zdd-array')).   11:- use_module(zdd(zdd)).   12:- use_module(pac(op)).   13%
   14udg_path(End, PathSet):- get_key(coa, Coa),
   15	set_key(ends, End),
   16	udg_mate_prune(Coa, 1, PathSet).
   17
   18		/******************************
   19		*     counting path in UDG    *
   20		******************************/
   21% ?- zdd.
   22% ?- udg_path_count(a*d, [a-b, b-c, c-d], C).  % fail.
   23% ?- udg_path_count(a-x, [a-b, b-c, c-d], C).  % fail.
   24% ?- udg_path_count(a-d, [a-b, b-c, c-d], C).
   25
   26% A big circular graph.
   27% ?- N = 10000,
   28%	findall(A,	( between(1, N, I), I<N, J is I + 1, A=(I-J)), Ls),
   29%	time(udg_path_count(1-2,[1-N |Ls], C)).
   30
   31% Counting paths in complete graph.
   32% ?- N = 11, findall(I-J, ( between(1, N, I), between(1, N, J), I < J), Ls),
   33%  time(zdd udg_path_count(1-N, Ls, C)).
   34%@ % 1,835,060,694 inferences, 233.677 CPU in 238.492 seconds (98% CPU, 7852993 Lips)
   35%@ N = 11,
   36%@ Ls = [1-2, 1-3, 1-4, 1-5, 1-6, 1-7, 1-8, 1-9, ... - ...|...],
   37%@ C = 986410.
   38
   39% [2024/01/04] Is global variable efficient ?
   40% ?- N = 11, findall(I-J, ( between(1, N, I), between(1, N, J), I < J), Ls),
   41%  time(zdd udg_path_count(1-N, Ls, C)).
   42%@ % 1,775,162,524 inferences, 211.846 CPU in 216.653 seconds (98% CPU, 8379506 Lips)
   43%@ N = 11,
   44%@ Ls = [1-2, 1-3, 1-4, 1-5, 1-6, 1-7, 1-8, 1-9, ... - ...|...],
   45%@ C = 986410.
   46%
   47udg_path_count(Ends, Links, C):- udg_path_count(Ends, Links, C, _).
   48
   49%
   50udg_path_count(Ends, Links, C, X):-
   51		prepare_udg(Ends, Links),
   52		ends_frontier(EF),
   53		!,
   54		get_key(coa, Coa),
   55		udg_mate_prune(EF, Coa, 1, X0),
   56		zdd_sort_paths(X0, X),
   57		card(X, C).
   58%
   59udg_mate_prune(EF, Ls, X, Y):-
   60	add_links(EF, Ls, X, Y0),
   61	get_key(ends, Pair),
   62	prune_final(Pair, Y0, Y).
   63
   64
   65% ?- zdd.
   66
   67% ?- time((rect_path_count(p(0,0)-p(1,1), rect(1,1), C, _))).
   68% ?- time((rect_path_count(p(0,0)-p(2,2), rect(2,2), C, _))).
   69% ?- time((rect_path_count(p(3,3)-p(4,4), rect(6,6), C, _))).
   70% ?- time((rect_path_count(p(0,0)-p(6,6), rect(6,6), C, _))).
   71% ?- time(rect_path_count(rect(1,0), C)).
   72% ?- time(rect_path_count(rect(1,1), C)).
   73% ?- time(rect_path_count(rect(2,1), C)).
   74% ?- time(rect_path_count(rect(1,2), C)).
   75% ?- time(rect_path_count(rect(1,3), C)).
   76% ?- time(rect_path_count(rect(3,1), C)).
   77% ?- time(rect_path_count(rect(4,1), C)).
   78% ?- time(rect_path_count(rect(5,1), C)).
   79% ?- time(rect_path_count(rect(1,2), C)).
   80% ?- time(rect_path_count(rect(1,3), C)).
   81% ?- time(rect_path_count(rect(1,4), C)).
   82% ?- time(rect_path_count(rect(1,5), C)).
   83% ?- time(rect_path_count(rect(2,2), C)).
   84% ?- time(rect_path_count(rect(3,3), C)).
   85%@ % 552,165 inferences, 0.037 CPU
   86%@ C = 184.
   87% ?- time(rect_path_count(rect(4,4), C)).
   88%@ % 4,588,052 inferences, 0.303 CPU
   89%@ C = 8512.
   90
   91% ?- time(rect_path_count(rect(4,4), C, X)).
   92%@ % 4,624,626 inferences, 0.307 CPU
   93%@ C = 8512,
   94%@ X = 2611.
   95% ?- time(rect_path_count(rect(5,5), C)).
   96%@ % 38,687,913 inferences, 1.799 CPU
   97%@ C = 1262816.
   98% ?- time(rect_path_count(rect(6,6), C)).
   99%@ % 105,853,302 inferences, 10.325 CPU in 12.362 seconds (84% CPU, 10252435 Lips)
  100%@ C = 575780564.
  101%@ % 105,853,302 inferences, 10.255 CPU in 12.287 seconds (83% CPU, 10322399 Lips)
  102%@ C = 575780564.
  103
  104% ?- time(rect_path_count(rect(7,7), C)).
  105%@ % 531,968,509 inferences, 56.027 CPU in 66.247 seconds (85% CPU, 9494888 Lips)
  106%@ C = 789360053252.
  107%@ % 531,968,509 inferences, 57.000 CPU in 67.227 seconds (85% CPU, 9332805 Lips)
  108%@ C = 789360053252.
  109% ?- time(rect_path_count(rect(8,8), C)).
  110%@ % 2,619,632,064 inferences, 294.280 CPU in 345.288 seconds (85% CPU, 8901837 Lips)
  111%@ C = 3266598486981642.
  112
  113% ?- time(rect_path_count(rect(9,9), C)).
  114%@ % 12,584,483,532 inferences, 1528.588 CPU in 1775.464 seconds (86% CPU, 8232752 Lips)
  115%@ C = 41044208702632496804.
  116
  117% ?- time(rect_path_count(rect(10,10), C)).
  118%@ % 59,632,082,303 inferences, 7437.316 CPU in 8649.928 seconds (86% CPU, 8017958 Lips)
  119%@ C = 1568758030464750013214100.
  120
  121% ?- time(rect_path_count(rect(11,11), C)).
  122%@ % 277,342,057,605 inferences, 35622.993 CPU in 41124.669 seconds (87% CPU, 7785479 Lips)
  123%@ C = 182413291514248049241470885236.
  124
  125% [2024/09/03] 13 x 13 grid graph passed also by the simple frontier vector.
  126% ?- time(rect_path_count(rect(12,12), C)).
  127% 1,273,378,663,129 inferences, 176187.026 CPU in 244597.424 seconds (72% CPU, 7227426 Lips)
  128% C = 64528039343270018963357185158482118.
  129
  130% ?- C = 182413291514248049241470885236,
  131%	C1 = 64528039343270018963357185158482118,
  132%	D is C1//C.
  133
  134% ?- forall(between(1, 13, I), (X is 2^(I*I), writeln(I=> X))).
  135
  136
  137rect_path_count(R, C):- rect_path_count(R, C, _).
  138%
  139rect_path_count(R, C, X):-
  140	R = rect(W, H),
  141	rect_path_count( p(0,0) - p(W,H), R, C, X).
  142%
  143rect_path_count(Pair, R, C, X):- rect_links(R, Links),
  144	udg_path_count(Pair, Links, C, X).
  145
  146% ?- rect_links(rect(1,1), Links).
  147% ?- rect_links(rect(10,10), Links),length(Links, L).
  148rect_links(rect(W, H), Links):-
  149	findall( p(I,J) - p(K,L),
  150				 (	between(0, W, I),
  151					between(0, H, J),
  152					(  L = J, K is I + 1, K =< W
  153					;  K = I, L is J + 1, L =< H
  154					)
  155				 ),
  156				 Links).
  157
  158		/********************************
  159		*     Prepare UDG in coalgebra  *
  160		********************************/
  161
  162% ?- udg_path_count(a-b,[a-b], C, X).
  163% ?- udg_path_count(a-c,[a-b, a-d, b-c, d-c], C).
  164% ?- udg_path_count(a-f,[a-b, a-d, b-c, d-c, d-e, e-f, c-f], C).
  165% ?- udg_path_count(a-k,
  166%	[a-b, a-d, b-c,
  167%	d-c, d-e, e-f,
  168%	c-f, b-g, g-h,
  169%	h-k, c-h, f-k], C).
  170
  171% ?- prepare_udg([a-b, b-c, c-d]), memo(node_id(a)-Id, memo_nodes),
  172% get_key(dom, X).
  173
  174% ?- prepare_udg(a-d, [a-b, b-c, c-d]),
  175%	get_key(ends, ST),
  176%	get_key(frontier, F),
  177%	get_key(dom, Dom),
  178%	get_key(coa, U).
  179%@ ST = 1-4,
  180%@ F = #(..),
  181%@ Dom = [1, 2, 3, 4],
  182%@ U = [4-[], 3-[4], 2-[3], 1-[2]].
  183
  184prepare_udg(ST, Links):-
  185	prepare_udg(Links),
  186	prepare_ends(ST, Pair),
  187	set_key(ends, Pair).
  188%
  189prepare_udg(Links):-
  190	open_memo(memo_nodes),
  191	prepare_udg(Links, Succs, D, Vec),
  192	length(D, N),
  193	completing_succs(Succs, Succs0, N),
  194	set_key(coa, Succs0),
  195	set_key(dom, D),
  196	set_key(frontier, Vec).
  197%
  198prepare_udg(Links, Succs, D, Vec):-
  199	intern_links(Links, Links0),
  200	normal_mate_list(Links0, Links1),
  201	sort(Links1, Links2),
  202	domain_of_links(Links2, D),
  203	rel_to_fun(Links2, Succs),
  204	Vec=..[#|D],
  205	setup_frontier(Links1, Vec).
  206%
  207prepare_ends(A-B, R):-!, R = A0-B0,
  208	memo(node_id(A)-I, memo_nodes),
  209	memo(node_id(B)-J, memo_nodes),
  210	(	nonvar(I), nonvar(J) -> sort([I, J], [A0, B0])
  211	;	format("No link at ~w or ~w\n", [A,B]),
  212		fail
  213	).
  214prepare_ends(E, _):-
  215	format("Unexpected form of end nodes ~w \n", [E]),
  216	fail.
  217
  218% ?-completing_succs([], Y, 2).
  219% ?-completing_succs([2-[a]], Y, 3).
  220completing_succs(X, X, 0):-!.
  221completing_succs([I-A|Ls], [I-A|Ms], I):-!, J is I - 1,
  222	completing_succs(Ls, Ms, J).
  223completing_succs(Ls, [I-[]|Ms], I):- J is I - 1,
  224	completing_succs(Ls, Ms, J).
  225
  226% ?- normal_mate_list([1-2], X).
  227% ?- normal_mate_list([2-1, 1-2], X).
  228normal_mate_list([], []).
  229normal_mate_list([P|R], [P0|R0]):- P = I-J,
  230	(	J @< I -> P0 = J - I
  231	;	P0 = P
  232	),
  233	normal_mate_list(R, R0).

 rel_to_fun(+R, -F) is det

convert set R of links to a list F of successor lists. In other words: F is a function derived from the relation R such that F(x) = P (x in dom(R)) if P = { y | R(x,y)} e.g. R=[a-b, a-c, b-d, b-e] => F=[a-[b,c], b-[d,e]]

  241% ?- rel_to_fun([], R).
  242% ?- rel_to_fun([a-b, a-c, b-d, b-e], R).
  243rel_to_fun(L, R):- sort(L, L0), rel_to_fun(L0, [], R).
  244%
  245rel_to_fun([], X, X).
  246rel_to_fun([A-B|L], [A-U|V], R):-!,
  247	rel_to_fun(L, [A-[B|U]|V], R).
  248rel_to_fun([A-B|L], U, R):-!,
  249	rel_to_fun(L, [A-[B]|U], R).
  250
  251% ?- domain_of_links([a-b, b-c, a-c], Y).
  252domain_of_links(X, Y):-
  253	findall(A , (	member(L, X),
  254					( L = (A - _)
  255					; L = (_ - A)
  256					)
  257				),
  258		   Y0),
  259   sort(Y0, Y).
  260
  261% ?- open_memo(memo_nodes), node_id(a, 0, C).
  262node_id(N, C, C0):- node_id(N, _, C, C0).
  263
  264% ?- open_memo(memo_nodes),
  265%   numlist(1, 10000, Ns),
  266%	foldl(pred(( [I, C, C0]:- node_id(st(I), K, C, C0))), Ns, 0, R).
  267node_id(N, I, C, C0):- memo(node_id(N)-I, memo_nodes),
  268	(	nonvar(I) -> C0 = C
  269	;	C0 is C+1,
  270		I = C0
  271	).
  272
  273% ?- open_memo(memo_nodes), intern_links([a-b, b-a], R).
  274intern_links(L, L0):- intern_links(L, L0, 0, _).
  275%
  276intern_links([], [], C, C).
  277intern_links([A-B|L], [I-J|M], C, D):-
  278	node_id(A, I, C, C0),
  279	node_id(B, J, C0, C1),
  280	intern_links(L, M, C1, D).
  281
  282		/******************
  283		*     frontier    *
  284		******************/

 off_frontier(+I, +J, +F) is det

True if node I must not accessible directly by a link from a node less than J, w.r.t. the froniter vector F.

  291% ?- X=f(1,2,3), setup_frontier([1-2,2-3], X), off_frontier(1, 3, X). %false
  292% ?- X=f(1,2,3), setup_frontier([1-2,2-3], X), off_frontier(3, 2, X). %true
  293
  294off_frontier(I, J, F):- arg(I, F, K), J =< K.

 on_frontier(+I, +J, +F) is det

True if node I may be accessible directly by a link from a node less than J, w.r.t. the frontier vector F.

  301% ?- X=f(1,2,3), setup_frontier([1-2,2-3], X), on_frontier(1, 3, X).
  302on_frontier(I, J, F):- arg(I, F, K), K < J.
  303
  304
  305% ?- X=f(1,2,3), setup_frontier([1-2,2-3], X).
  306setup_frontier([], _).
  307setup_frontier([I-J|L], F):-
  308	update_frontier(I, J, F),
  309	!,
  310	setup_frontier(L, F).
  311
  312% ?- X=f(1,2), update_frontier(1,2, X).
  313% ?- X=f(1,2,3), update_frontier(2, 3, X), update_frontier(1, 2, X).
  314update_frontier(I, J, V):-
  315	arg(J, V, A),
  316	(	I < A -> setarg(J, V, I)
  317	;	true
  318	).
  319
  320		/*******************
  321		*	 Helpers       *
  322		*******************/
  323
  324% ?- arrow_symbol(_->_, F).
  325% ?- arrow_symbol(a->b, F, X, Y).
  326arrow_symbol( _ -> _).
  327%
  328arrow_symbol(A, A0):- functor(A, A0, 2).
  329arrow_symbol(A, A0, A1, A2):- functor(A, A0, 2),
  330		arg(1, A, A1),
  331		arg(2, A, A2).
  332
  333% The most basic helpers.
  334composable_pairs(A-B, A-C, B, C).
  335composable_pairs(A-B, C-A, B, C).
  336composable_pairs(B-A, A-C, B, C).
  337composable_pairs(B-A, C-A, B, C).
  338%
  339normal_pair(A-B, B-A):- B < A, !.
  340normal_pair(A->B, B->A):- B < A, !.
  341normal_pair(X, X).
  342%
  343strong_less_than(_-A, B-_):- A<B.
  344%
  345ends_frontier(efr(E, Fr)):-
  346	get_key(ends, E),
  347	get_key(frontier, Fr).
  348
  349		/************************
  350		*     core predicates   *
  351		************************/
  352add_links(_, [], X, X):-!,
  353	writeln(done).
  354add_links(EF, [A-Ns|Ls], X, Y):-!,
  355	writeln(A),
  356	bdd_cons(X0, A-A, X),
  357	add_links_(A, Ns, X0, X1),
  358	prune_by_frontier(EF, X1, X2, A),
  359	slim_gc(X2, X3),
  360	add_links(EF, Ls, X3, Y).
  361
  362%
  363add_links_(_, [], X, X):-!.
  364add_links_(A, [B|Ns], X, Y):- add_link(A-B, X, X0), !,
  365	zdd_join(X, X0, X1),
  366	add_links_(A, Ns, X1, Y).
  367%
  368add_link(_, X, 0):- X<2, !.
  369add_link(U, X, Y):- % memo useless here.
  370	cofact(X, t(A, L, R)),
  371	add_link(U, L, L0),
  372	(	A = (_->_) -> R0 = 0
  373	;	U = A  -> R0 = 0	% cycle found
  374	;   strong_less_than(U, A) -> R0 = 0  %
  375	;  (	composable_pairs(U, A, V, W) ->
  376			U = (Ul-Ur),
  377		    subst_node([Ul->Ur], V, W, R, R0)
  378	   ;	add_link(U, R, R1),
  379		    zdd_insert(A, R1, R0)
  380	   )
  381	),
  382	zdd_join(L0, R0, Y).

  385subst_node(_, _, _, X, 0):- X<2, !.
  386subst_node(Es, A, P, X, Y):- % memo useless here
  387	cofact(X, t(U, L, R)), % replace A with P
  388	subst_node(Es, A, P, L, L0),
  389	arrow_symbol(U, F, Lu, Ru),
  390	(	F = (->) ->  R0 = 0
  391	;	Ru = A	->
  392		normal_pair(Lu-P, V),
  393		mate_ord_insert([V|Es], R, R0)
  394%		zdd_ord_insert([V|Es], R, R0)
  395	;	Lu = A	->
  396		normal_pair(P-Ru, V),
  397		mate_ord_insert([V|Es], R, R0)
  398%		zdd_ord_insert([V|Es], R, R0)
  399	;	subst_node(Es, A, P, R, R1),
  400		zdd_insert(U, R1, R0)
  401	),
  402	zdd_join(L0, R0, Y).
  403
  404% Despite of extra normalizing,
  405% a little bit better than zdd_ord_insert.
  406%
  407mate_ord_insert(_, 0, 0):-!.
  408mate_ord_insert([], X, X):-!.
  409mate_ord_insert(As, 1, Y):-!, bdd_append(As, 1, Y).
  410mate_ord_insert(As, X, Y):-		% memo is useless here.
  411		As = [A|U],
  412		(	A = (_->_) ->
  413			insert_arrows_after_links(As, X, Y)
  414		;   cofact(X, t(B, L, R)),
  415			zdd_compare(C, A, B),
  416			(	C = (<) ->
  417				mate_ord_insert(U, X, Y0),
  418				bdd_cons(Y, A, Y0)
  419			;	C = (=) ->
  420				mate_ord_insert(U, L, L1),
  421				mate_ord_insert(U, R, R1),
  422				zdd_join(R1, L1, Y0),
  423				bdd_cons(Y, A, Y0)
  424			;	mate_ord_insert(As, L, L1),
  425				mate_ord_insert(As, R, R1),
  426				cofact(Y, t(B, L1, R1))
  427			)
  428		).
  429%
  430insert_arrows_after_links(_, 0, 0):-!.
  431insert_arrows_after_links([], X, X):-!.
  432insert_arrows_after_links(As, 1, Y):-!, bdd_append(As, 1, Y).
  433insert_arrows_after_links(As, X, Y):- cofact(X, t(B, L, R)),
  434	(	B = (_->_) -> bdd_append(As, X, Y)
  435	;	insert_arrows_after_links(As, L, L0),
  436		insert_arrows_after_links(As, R, R0),
  437		cofact(Y, t(B, L0, R0))
  438	).
  439
  440		/***************************
  441		*     Prune by frontier    *
  442		***************************/
  443
  444% prune_by_frontier(EF, I, X, Y):-
  445% 	prune_by_frontier(X, Y, I, Pair, V).

 prune_by_frontier(+X, -Y, +I, +M, +V) is det

Y is unified with pruned X. 1) A-A is removed from path that has A-A when A is off_frontier. 2) Path which has A-B with off_frontier A or B is removed.

  452prune_by_frontier(_, X, X, _):- X<2, !.
  453prune_by_frontier(EF, X, Y, I):- cofact(X, t(A, L, R)),
  454	EF = efr(Pair, V),
  455	(	A = (_->_) -> Y = X
  456	; 	prune_by_frontier(EF, L, L0, I),
  457		classify_pair(A, I, Pair, V, C),
  458		(	C = arrow -> bdd_cons(R0, A, R)
  459		; 	C = keep ->
  460			prune_by_frontier(EF, R, R1, I),
  461			zdd_insert(A, R1, R0)
  462		;	C = ignore ->
  463			prune_by_frontier(EF, R, R0, I)
  464		;	R0 = 0
  465		),
  466		zdd_join(L0, R0, Y)
  467	).
  468
  469% helper.
  470on_pair(I, J-K):- K=I; J=I.
  471
  472% works! [2024/01/11]
  473classify_pair((_->_), _, _, _, arrow):-!.
  474classify_pair(J-J, I, _, V, C):-!,
  475	(	on_frontier(J, I, V) -> C = keep
  476	;	C = ignore
  477	).
  478classify_pair(J-K, I, Pair, V, C):-		% J\==K.
  479	(	on_frontier(J, I, V) ->
  480		(	on_frontier(K, I, V) -> C = keep
  481		;	on_pair(K, Pair) -> C = keep
  482		;	C = 0
  483		)
  484	;	on_frontier(K, I, V) ->
  485		(	on_pair(J, Pair) -> C = keep
  486		;	C = 0
  487		)
  488	;	(J-K) = Pair -> C = keep
  489	;   C = 0
  490	).
  491% ?- zdd.
  492% ?- X<< {[1-6,2-2,5-5,(1->3),(3->4),(4->6)]},
  493%	prune_final(1-6, X, Y), psa(Y).
  494
  495% ?- X<< +[*[a-b, a->b]], prune_final(a-b, X, Y), psa(X), psa(Y).
  496prune_final(_, X, 0):- X<2, !.
  497prune_final(Pair, X, Y):- cofact(X, t(A, L, R)),
  498	prune_final(Pair, L, L0),
  499	(	A = (_->_) -> R0 = 0, writeln('unexpected ***arrow***')
  500	;  	A = Pair -> prune_final0(R, R0)
  501	;	A = V-V -> prune_final(Pair, R, R0)
  502	;	R0 = 0
  503	),
  504	zdd_join(L0, R0, Y).
  505%
  506prune_final0(X, X):- X<2, !.
  507prune_final0(X, Y):- cofact(X, t(A, L, R)),
  508	(	A = (_-_) -> Y = X
  509	; 	prune_final0(L, L0),
  510		(	A = (_->_) -> bdd_cons(R0, A, R)
  511		;	A = (V-V) -> prune_final0(R, R0)
  512		;	R0 = 0
  513		),
  514		zdd_join(L0, R0, Y)
  515	).
  516
  517%
  518hamilton_prune_final(P, P, _, 1):-!.
  519hamilton_prune_final(_, _, X, 0):- X<2, !.
  520hamilton_prune_final(P, Q, X, Y):- cofact(X, t(A, L, R)),
  521	hamilton_prune_final(P, Q, L, L0),
  522	(	A = (_->_) -> Y = 0
  523	;  	A = P-Q -> hamilton_prune_final0(R, R0)
  524	;	R0 = 0
  525	),
  526	zdd_join(L0, R0, Y).
  527%
  528hamilton_prune_final0(X, X):- X<2, !.
  529hamilton_prune_final0(X, Y):- cofact(X, t(A, L, _R)),
  530	(	A = (_->_) -> Y = X
  531	; 	hamilton_prune_final0(L, L0),
  532		R0 = 0,
  533		zdd_join(L0, R0, Y)
  534	).
  535
  536% ?- X<<{[1->2, 1->3, 2->3]},
  537%	hamilton_filter([1,2,3], X, Y), card(Y, C).
  538% ?- X<<{[1->2]},
  539%	hamilton_filter([1,2], X, Y), card(Y, C).
  540% ?- X<<{[1->2, 1->3, 2->4, 3->4]},
  541%	hamilton_filter([1,2,3,4], X, Y), card(Y, C).
  542
  543hamilton_filter([I, J], X, Y):-!, hamilton_filter_special(I, J, X, Y).
  544hamilton_filter(D, X, Y):- hamilton_filter_list(D, X, Y).
  545%
  546hamilton_filter_special(_, _, X, 0):- X<2, !.
  547hamilton_filter_special(I, J, X, Y):- cofact(X, t(A, L, R)),
  548	hamilton_filter_special(I, J, L, L0),
  549	(	A=(I->J) ->
  550		(	R = 1 -> R0 = 1
  551		;	R0 = 0
  552		)
  553	;  R0 = 0
  554	),
  555	cofact(Y, t(A, L0, R0)).
  556%
  557hamilton_filter_list([], X, X).
  558hamilton_filter_list([J|Js], X, Y):-
  559	hamilton_filter(2, J, X, X0),
  560	hamilton_filter_list(Js, X0, Y).
  561%
  562hamilton_filter(0, I, X, Y):-!, without_node(I, X, Y).
  563hamilton_filter(_, _, X, 0):- X<2, !.
  564hamilton_filter(K, I, X, Y):- memo(hamilton(K, I, X)-Y),
  565	(	nonvar(Y) -> true	%,	write(.)
  566	;	cofact(X, t(A, L, R)),
  567		hamilton_filter(K, I, L, L0),
  568		(	A=(E-_), I < E -> R0 = 0
  569		;	(	on_arrow(I, A) ->
  570				K0 is K-1,
  571				hamilton_filter(K0, I, R, R0)
  572			;	hamilton_filter(K, I, R, R0)
  573			)
  574		),
  575		cofact(Y, t(A, L0, R0))
  576	).
  577%
  578without_node(_, X, X):- X<2, !.
  579without_node(I, X, Y):- memo(without_node(I,X)-Y),
  580	(	nonvar(Y) -> true	%, write(+)
  581	;	cofact(X, t(A, L, R)),
  582		without_node(I, L, L0),
  583		(	on_arrow(I, A) -> R0 = 0
  584		;	without_node(I, R, R0)
  585		),
  586		cofact(Y, t(A, L0, R0))
  587	).
  588%
  589on_arrow(X, Y->Z):- (X=Y; X=Z).
  590
  591% Sort each path eliminating duplicate occurrence.
  592% [b, c, a, b, c] => [a,b,c]
  593
  594zdd_sort_paths(X, X):- X<2, !.
  595zdd_sort_paths(X, Y):- memo(sort_paths(X)-Y),
  596	(	nonvar(Y)->true %, write(.) % many hits.
  597	;	cofact(X, t(A, L, R)),
  598		zdd_sort_paths(L, L0),
  599		zdd_sort_paths(R, R0),
  600		zdd_join(L0, R0, X0),
  601		zdd_insert(A, X0, Y)
  602	)