1:- module(path_count4, [
    2	]).    3
    4:- use_module(library(apply)).    5:- use_module(library(apply_macros)).    6:- use_module(library(clpfd)).    7:- use_module(library(statistics)).    8:- use_module(zdd('zdd-array')).    9:- use_module(util(math)).   10:- use_module(util(meta2)).   11:- use_module(pac(basic)).   12:- use_module(pac(meta)).   13:- use_module(util(misc)).   14:- use_module(pac('expand-pac')). % For the kind block.
   15:- use_module(zdd('zdd-misc')).   16:- use_module(pac(op)).   17:- use_module(zdd(zdd)).   18:- use_module(zdd('update-link')).   19
   20:- set_prolog_flag(stack_limit, 10_200_147_483_648).   21
   22 :- op(1060, xfy, ~).		% equivalence
   23 :- op(1060, xfy, #).		% exor
   24 :- op(1060, xfy, <->).		% equivalence
   25 :- op(1050, yfx, <-).   26 :- op(1060, xfy, <=> ).	% equivalence
   27 :- op(1040, xfy, \/).		% OR
   28 :- op(1030, xfy, /\).		% AND
   29 :- op(1020, fy, \).		% NOT
   30 :- op(700, xfx, :=).		% Assignment
   31 :- op(1000, xfy, &).   32
   33% for pac query.
   34 :- pac:op(1000, xfy, &).   35 :- pac:op(700, xfx, :=).   36
   37term_expansion --> pac:expand_pac.
   38
   39
   40% ?- zdd((dfs_dg_path_count(dg([a,b], [a-b]), C, a-b))).
   41% ?- zdd((dfs_dg_path_count(dg([a,b], [a-b, b-a]), C, a-b))).
   42% ?- zdd((dfs_dg_path_count(dg([a,b,c], [a-b]), C, a-b))).
   43% ?- zdd((dfs_dg_path_count(dg([a,b,c], [a-b, b-c]), C, a-c))).
   44% ?- zdd((dfs_dg_path_count(dg([a,b,c,d], [a-b, b-c, b-d, a-c, c-b]), C, a-d))).
   45% ?- R = rect(1,1),  zdd((rect_dg(R, N, E), dfs_dg_path_count(dg(N, E), C, p(0,0)-p(1,1)))).
   46
   47% ?-zdd((dfs_rect_path_count(rect(1,2), C))).
   48%@ C = 4 .
   49% ?-time(zdd((dfs_rect_path_count(rect(15,1), C)))).
   50%@ % 106,749,564 inferences, 11.306 CPU in 12.134 seconds (93% CPU, 9441435 Lips)
   51%@ C = 32768 .
   52% ?-time(zdd((dfs_rect_path_count(rect(2,2), C)))).
   53%@ % 17,867 inferences, 0.031 CPU in 0.057 seconds (54% CPU, 583412 Lips)
   54%@ C = 12 .
   55
   56% ?-time(zdd((dfs_rect_path_count(rect(3,3), C)))).
   57%@ % 484,761 inferences, 0.196 CPU in 0.340 seconds (58% CPU, 2468610 Lips)
   58%@ C = 184 .
   59% ?-time(zdd((dfs_rect_path_count(rect(4,4), C)))).
   60%@ % 33,077,014 inferences, 3.574 CPU in 4.034 seconds (89% CPU, 9254250 Lips)
   61%@ C = 8512 .
   62% ?-time(zdd((dfs_rect_path_count(rect(5,5), C)))).
   63%@ % 7,143,092,797 inferences, 898.272 CPU in 962.381 seconds (93% CPU, 7952034 Lips)
   64%@ C = 1262816 .
   65
   66dfs_rect_path_count(R, C):- @(dfs_rect_path_count(R, C)).
   67%
   68dfs_rect_path_count(R, C, S):-	R = rect(W, H),
   69	rect_dg(R, N, E),
   70	dfs_dg_path_count(dg(N, E), C, p(0,0)-p(W, H), S).
   71
   72%
   73dfs_dg_path_count(DG, C, ST):- @(dfs_dg_path_count(DG, C, ST)).
   74%
   75
   76dfs_dg_path_count(DG, C, ST, S):- DG = dg(N, E),
   77	ST = A-_,
   78	build_suc(dg(N, E), S),
   79	findall(X-X, member(X, N), Id),
   80	linear_zdd(Id, Zid, S),
   81	dfs_build_qp(Zid, A, ST, S),
   82	count_child_path(Zid, C, S).
   83%
   84dfs_build_qp(X, A, ST, S):- memo(visited(X, A)-B, S),
   85	(	nonvar(B) -> true
   86	;	B = visited,
   87		memo(suc(A)-M, S),
   88		dfs_build_qp(M, X, A, ST, S)
   89	).
   90%
   91dfs_build_qp([], _, _, _, _).
   92dfs_build_qp([N|Ns], X, A, ST, S):-
   93	dfs_qp_add_child(A-N, X, ST, S),
   94	dfs_build_qp(Ns, X, A, ST, S).
   95
   96% ?- zdd((X<<{[a-a, b-b]}, dfs_qp_add_child(a-b, X, a-b),
   97%	memo(child(X)-U),
   98%	dfs_qp_add_child(b-a, X, a-b),
   99%	memo(child(X)-V))).
  100
  101dfs_qp_add_child(E, X, ST):- @(dfs_qp_add_child(E, X, ST)).
  102%
  103
  104dfs_qp_add_child(E, X, ST, S):- update_qp(E, X, Y, S),
  105	(	Y = 0 ->  true
  106	;	(	final_qp(Y, ST, S) ->  Y0 = 1
  107		;	Y0 = Y
  108		),
  109		add_child(child(X), E-Y0, S),
  110		E = _- R,
  111		dfs_build_qp(Y0, R, ST, S)		% Depth-first search.
  112	).
  113
  114% ?- zdd((X << {[1-2]}, final_qp(X, 1-2))).
  115% ?- zdd((X << {[1-2,3-4]}, final_qp(X, 1-2))).
  116% ?- zdd((X << {[1-1,2-2]}, final_qp(X, 1-2))).
  117% ?- zdd((X << {[]}, final_qp(X, 1-2))).
  118
  119final_qp(X, ST):- @(final_qp(X, ST)).
  120%
  121final_qp(X, ST, S):- X>1,
  122	cofact0(X, A, Y, S),
  123	(	A = (I-I) -> final_qp(Y, ST, S)
  124	;	A = ST -> qp_id(Y, S)
  125	).
  126%
  127qp_id(1, _):-!.
  128qp_id(X, S):- X>1,
  129	cofact0(X, I-I, Y, S),
  130	qp_id(Y, S).
  131
  132% ?- zdd((build_suc(dg([a,b], [a-b])), memo(suc(a)-U))).
  133% ?- zdd((build_suc(dg([a,b,c], [a-b, a-c, a-b, b-a])), memo(suc(a)-U),
  134%	memo(suc(b)-V), memo(suc(c)-W))).
  135
  136build_suc(R):- @(build_suc(R)).
  137%
  138build_suc(dg(Ns, Rel), S):-
  139	maplist(pred(S, [X]:- memo(suc(X)-[], S)), Ns),
  140	rel_to_coa(suc, Rel, S).
  141
  142% ?- zdd((rel_to_coa(suc, [a-b]), memo(suc(a)-U))).
  143% ?- zdd((rel_to_coa(suc, [a-b, a-c, a-b, b-a]), memo(suc(a)-U),
  144%	memo(suc(b)-V), memo(suc(c)-W))).
  145
  146rel_to_coa(Key, R):- @(rel_to_coa(Key, R)).
  147%
  148rel_to_coa(_, [], _):-!.
  149rel_to_coa(Key, [A-B|R], S):-functor(U, Key, 1),
  150	arg(1, U, A),
  151	add_child(U, B, S),
  152	rel_to_coa(Key, R, S).
  153%
  154count_child_path(X, C):- @(count_child_path(X, C)).
  155%
  156count_child_path(X, X, _):- X<2, !.
  157count_child_path(X, C, S):- memo(count_child_path(X)-C, S),
  158	(	nonvar(C) -> true
  159	;	memo(child(X) - M, S),
  160		count_child_path(M, 0, C, S)
  161	).
  162%
  163count_child_path([], C, C, _).
  164count_child_path([_-X|R], C, D, S):- count_child_path(X, U, S),
  165	C0 is U + C,
  166	count_child_path(R, C0, D, S).
  167
  168		/**************************************
  169		*     breadth first search version    *
  170		**************************************/
  171
  172
  173% ?-time(zdd((bfs_rect_path_count(rect(15,0), C)))).
  174% ?-time(zdd((bfs_rect_path_count(rect(0, 120), C)))).
  175% ?-time(zdd((bfs_rect_path_count(rect(1,1), C)))).
  176% ?-time(zdd((bfs_rect_path_count(rect(2,2), C)))).
  177%@ % 17,728 inferences, 0.002 CPU in 0.002 seconds (97% CPU, 8330827 Lips)
  178%@ C = 12 .
  179% ?-time(zdd((bfs_rect_path_count(rect(3,3), C)))).
  180%@ % 483,059 inferences, 0.044 CPU in 0.045 seconds (97% CPU, 11082639 Lips)
  181%@ C = 184 .
  182% ?-time(zdd((bfs_rect_path_count(rect(4,4), C)))).
  183
  184%@ % 33,025,049 inferences, 3.659 CPU in 3.785 seconds (97% CPU, 9024986 Lips)
  185%@ C = 8512 .
  186
  187% ?-profile(time(zdd((bfs_rect_path_count(rect(4,4), C))))).
  188%@ % 33,025,163 inferences, 4.752 CPU in 4.915 seconds (97% CPU, 6950417 Lips)
  189%@ C = 8512.
  190
  191
  192
  193
  194% ?-time(zdd((bfs_rect_path_count(rect(5,5), C)))).
  195%@ % 5,517,389,364 inferences, 816.623 CPU in 842.814 seconds (97% CPU, 6756348 Lips)
  196%@ C = 1262816 .
  197
  198bfs_rect_path_count(R, C):- @(bfs_rect_path_count(R, C)).
  199%
  200bfs_rect_path_count(R, C, S):-	R = rect(W, H),
  201	rect_dg(R, N, E),
  202	bfs_dg_path_count(dg(N, E), C, p(0,0)-p(W, H), S).
  203%
  204bfs_dg_path_count(DG, C, ST):- @(bfs_dg_path_count(DG, C, ST)).
  205%
  206bfs_dg_path_count(DG, C, ST, S):- DG = dg(N, E),
  207	ST= A-_,
  208	build_suc(dg(N, E), S),
  209	findall(X-X, member(X, N), Id),
  210	linear_zdd(Id, Zid, S),
  211	bfs_build_qp([A-Zid], [], ST, S),
  212	count_child_path(Zid, C, S).
  213%
  214% bfs_build_qp(X, Y, _, _):- length(X, L), length(Y, M), writeln(L-M), fail.
  215bfs_build_qp([], [], _, _):-!.
  216bfs_build_qp([], X, ST, S):-!, bfs_build_qp(X, [], ST, S).
  217bfs_build_qp([A-X|R], R0, ST, S):- memo(visited(X, A)-B, S),
  218	(	nonvar(B)
  219	->	% write(.),  % many hits.
  220		bfs_build_qp(R, R0, ST, S)
  221	;	B = visited,
  222		memo(suc(A)-M, S),
  223		bfs_build_qp(M, X, A, U-R0, ST, S),
  224		bfs_build_qp(R, U, ST, S)
  225	).
  226
  227%
  228bfs_build_qp([], _, _, U-U, _, _).
  229bfs_build_qp([N|Ns], X, A, U-V, ST, S):-
  230	bfs_qp_add_child(A-N, X, Y, ST, S),
  231	(	( Y = 0; Y = 1) -> U0 = U
  232	;	U = [N-Y|U0]
  233	),
  234	bfs_build_qp(Ns, X, A, U0-V, ST, S).
  235
  236% ?- zdd((X<<{[a-a, b-b]}, qp_add_child(a-b, X, a-b),
  237%	memo(child(X)-U),
  238%	qp_add_child(b-a, X, a-b),
  239%	memo(child(X)-V))).
  240
  241bfs_qp_add_child(E, X, Y, ST):- @(bfs_qp_add_child(E, X, Y, ST)).
  242%
  243bfs_qp_add_child(E, X, Y, ST, S):- update_qp(E, X, Y0, S),
  244	(	Y0 = 0 ->  Y = 0
  245	;	(	final_qp(Y0, ST, S) ->  Y = 1
  246		;	Y = Y0
  247		),
  248		add_child(child(X), E-Y, S)
  249	).
  250
  251
  252		/****************
  253		*     bridge    *
  254		****************/
  255
  256% ?- zdd((count_path_line(p(1,1),p(2,1), C))).
  257% ?- zdd((count_path_line(p(1,2), p(3,2), C))).
  258% ?- zdd((count_path_line(p(1,2),p(10,2), C))).
  259% ?- zdd((count_path_line(p(1,2),p(100,2), C))).
  260% ?- time(zdd((count_path_line(p(1,2),p(1000, 2), C)))).
  261%@ % 4,691,349 inferences, 0.982 CPU in 1.154 seconds
  262% ?- time(zdd((count_path_line(p(1,2),p(10000,2), C)))).
  263%@ % 61,745,569 inferences, 59.803 CPU in 60.290 seconds
  264
  265count_path_line(X, Y, C):- @(count_path_line(X, Y, C)).
  266%
  267count_path_line(X, Y, C, S):-
  268	qp_line(X, Y, R, S),
  269	memo_final_line(R, X, Y, S),
  270	qp_list(Q, [X-X], S),
  271	path_tree(Q, Tree, S),
  272	card(Tree, C, S).
  273%
  274memo_final_line([], _, _, _).
  275memo_final_line([A|R], X, Y, S):-
  276	qp_list(A, List, S),
  277	(	List = [X-Y]
  278	->  memo(qp_final(A)-true, S)
  279	;  true
  280	),
  281	memo_final_line(R, X, Y, S).
  282
  283%
  284qp_line(X, Y, R):- @(qp_line(X, Y, R)).
  285%
  286qp_line(X, X, [Q], S):-!,
  287	zdd_insert_atoms([X-X], 1, Q, S),
  288	memo(qp_suc(Q)-List, S),
  289	( var(List)-> List=[]; true).
  290qp_line(X, Y, R, S):- X = p(I, K), Y = p(J, _),
  291	M is (I+J)//2,
  292	qp_line(X, p(M, K), R1, S),
  293	M0 is M+1,
  294	qp_line(p(M0,K), Y, R2, S),
  295	qp_joint([p(M,K)-p(M0,K), p(M0,K)-p(M,K)], R1, R2, R3, S),
  296	prune_qp_line_local(R3, X, Y, R, S).
  297%
  298prune_qp_line_local([], _, _, [], _).
  299prune_qp_line_local([Q|Qs], X, Y, Rs, S):-
  300	X = p(I,K),
  301	Y = p(J,K),
  302	qp_list(Q, List, S),
  303	(	member(p(A,_)-p(B,_), List),
  304		A \== B,
  305		(	I<A, A<J
  306		;	I<B, B<J
  307		)
  308	->  Rs0 = Rs
  309	;   Rs = [Q|Rs0]
  310	),
  311	prune_qp_line_local(Qs, X, Y, Rs0, S).
  312
  313%% test basics.
  314% ?- zdd((
  315%	qp_list(X, [a-a]), qp_list(Y, [b-b]),
  316%	qp_joint([a-b, b-a], [X], [Y], Z),
  317%	maplist(pred(([U]:-qp_list(U, List), writeln(List))), Z))).
  318% ?- zdd((
  319%	qp_list(X, [a-a]), qp_list(Y, [b-b]),
  320%	qp_joint([a-b, b-a], [X], [Y], Z),
  321%	zdd_join(X, Y, A),
  322%	memo(qp_suc(A)-L), writeln(qp_suc(A)-L),
  323%	maplist(pred(([U]:-qp_list(U, List), writeln(List))), Z))).
  324% ?- zdd((
  325%	qp_list(X, [a-a]), qp_list(Y, [b-b]),
  326%	qp_joint([a-b, b-a], [X], [Y], Z),
  327%	zdd_join(X, Y, A),
  328%	memo(qp_suc(A)-L), writeln(qp_suc(A)-L),
  329%	maplist(pred(([U]:-qp_list(U, List), writeln(List))), Z))).
  330
  331qp_joint(X, Y, Z, U):-@(qp_joint(X, Y, Z, U)).
  332%
  333qp_joint(X, Y, Z, U, S):- map_prod_sum(Y, Z, V, [], S),
  334	qp_bridge_links(X, Y, Z, U, V, S).
  335
  336% ?- zdd((qp_list(X, [a-a]), qp_list(Y, [b-b]),
  337%	qp_bridge(a-b, X, Y, Z), qp_list(Z, U))).
  338% ?- zdd((qp_list(X, [a-a]), qp_list(Y, [b-b]),
  339%	qp_bridge(b-a, X, Y, Z), qp_list(Z, U))).
  340
  341qp_bridge(E, X, Y, Z):- @(qp_bridge(E, X, Y, Z)).
  342%
  343qp_bridge(A-B, X, Y, Z, S):- % assuming qp  X < qp Y
  344	qp_list(X, L, S),
  345	qp_list(Y, M, S),
  346	(	select(U-A, L, L0),
  347		select(B-V, M, M0)
  348	;	select(B-V, L, L0),
  349		select(U-A, M, M0)
  350	),
  351	!,
  352	ord_union(L0, M0, N0),
  353	ord_union([U-V], N0, N),
  354	qp_list(Z, N, S).
  355qp_bridge(_, _, _, 0, _).
  356
  357%
  358qp_bridge_links(X, Y, Z, U, V):- @(qp_bridge_links(X, Y, Z, U, V)).
  359%
  360qp_bridge_links([], _, _, Z, Z, _).
  361qp_bridge_links([E|Es], X, Y, Z, U, S):-
  362	qp_bridge_link(X, E, Y, Z, Z0, S),
  363	qp_bridge_links(Es, X, Y, Z0, U, S).
  364%
  365qp_bridge_link([], _, _, Z, Z, _).
  366qp_bridge_link([X|Xs], E, Y, Z, U, S):-
  367	qp_bridge_basic(Y, X, E, Z, Z0, S),
  368	qp_bridge_link(Xs, E, Y, Z0, U, S).
  369%
  370qp_bridge_basic([], _, _, Z, Z, _).
  371qp_bridge_basic([Y|Ys], X, E, Z, U, S):- qp_bridge(E, X, Y, V, S),
  372	(	V = 0
  373	->	Z0 = Z
  374	;	Z  = [V|Z0],
  375		select_goto(E, X, Y, V, S)
  376	),
  377	qp_bridge_basic(Ys, X, E, Z0, U, S).
  378%
  379select_goto(E, X, Y, Z):- @(select_goto(E, X, Y, Z)).
  380%
  381select_goto(A-B, X, Y, Z, S):-
  382	G = (A-B)-Z,
  383	( A@<B -> W = X ; W = Y ),
  384	getmemo(qp_suc(W)-U, V, S),
  385	(	var(V) -> V = [] ; true ),
  386	add_new(G, V, U).
  387
  388		/***********************
  389		*     end of bridge    *
  390		***********************/
  391
  392% ?-  belong(3, 3, 4).
  393% ?-  belong(3, 4, 1).
  394belong(X, I, J):- I @=< J, !, I @=< X, X @=< J.
  395belong(X, I, J):- J @=< X, X @=< I.
  396%
  397sub_interval(I, J, X, Y):- belong(I, X, Y), belong(J, X, Y).
  398%
  399disjoint_interval(I, J, X, Y):- \+ belong(I, X, Y), \+ belong(J, X, Y).
  400
  401% ?- compatible_interval(1, 2, 3, 4).
  402% ?- compatible_interval(1, 3, 2, 4). % false
  403% ?- compatible_interval(2, 4, 5, 7).
  404% ?- compatible_interval(7, 5, 2, 4).
  405
  406compatible_interval(I, J, X, Y):- disjoint_interval(I, J, X, Y), !.
  407compatible_interval(I, J, X, Y):- sub_interval(I, J, X, Y), !.
  408compatible_interval(I, J, X, Y):- sub_interval(X, Y, I, J).
  409
  410% ?- rect_dg(rect(1,1), X, Y).
  411% ?- rect_dg(rect(0,1), X, Y).
  412% ?- rect_dg(rect(0,0), X, Y).
  413
  414rect_dg(rect(W, H), Nodes, Links):-
  415	rect_nodes(W, H, Nodes),
  416	rect_links(W, H, Nodes, Links).
  417
  418% ?- rect_nodes(4,5,Ns), length(Ns, L).
  419rect_nodes(W, H, Nodes):-
  420	findall(p(I, J),
  421			(	between(0, W, I),
  422				between(0, H, J)
  423			),
  424			Nodes0),
  425	sort(Nodes0, Nodes).
  426
  427%
  428rect_links(W, H, Ns, Links):-
  429	findall(P-Q,
  430			(	member(P, Ns),
  431				P = p(I, J),
  432				(	Q = p(I, J1),
  433					(	J1 is J-1,
  434						J1 >= 0
  435					;	J1 is J+1,
  436						J1 =< H
  437					)
  438				;   Q = p(I1, J),
  439					(	I1 is I-1,
  440						I1 >= 0
  441					;	I1 is I+1,
  442						I1 =< W
  443					)
  444				)
  445			),
  446			Links0),
  447	sort(Links0, Links).
  448
  449		/********************
  450		*   frontier based  *
  451		********************/
  452
  453% ?- zdd(path_count(rect(0,0), C)).
  454% ?- zdd(path_count(rect(1,0), C)).
  455% ?- zdd(path_count(rect(0,1), C)).
  456% ?- zdd(path_count(rect(1,1), C)).
  457% ?- zdd(path_count(rect(2, 1), C)).
  458% ?- zdd(path_count(rect(3, 1), C)).
  459% ?- zdd(path_count(rect(1, 3), C)).
  460% ?- zdd(path_count(rect(4, 1), C)).
  461% ?- zdd(path_count(rect(2, 3), C)).
  462% ?- time(zdd(path_count(rect(4,3), C))).
  463% ?- time(zdd(path_count(rect(4,4), C))).
  464%@ C= 8512
  465% ?- call_with_time_limit(360, zdd(path_count(rect(5, 5), C))).
  466
  467path_count(Rect, C):- @(path_count(Rect, C)).
  468
  469%
  470path_count(rect(0,0), 1, _):-!.
  471path_count(Rect, C, S):-
  472	rect_dg(Rect, Ns, Links),
  473	findall(A-A, member(A, Ns), Ids),
  474	Rect = rect(W, H),
  475	set_key(final, p(W, H), S),
  476	qp_list(Q, Ids, S),
  477	memo(qp_suc(Q)-[], S),
  478	frontier(Links, [Q], FinalQs, S),
  479	fold_qp_memo_final(FinalQs, S),
  480	path_tree(Q, Tree, S),
  481	card(Tree, C, S).
  482
  483%
  484frontier([], Qs, Qs, _).
  485frontier([E|Es], Qs, Qs0, S):-
  486	frontier_link(E, Qs, Qs1, Qs, S),
  487	prune_frontier(Qs1, E, Qs2, S),
  488	frontier(Es, Qs2, Qs0, S).
  489
  490%
  491frontier_link(_, [], Qs, Qs, _):-!.
  492frontier_link(E, [X|Qs], Qs0, Qs1, S):- X<2, !,
  493 	frontier_link(E, Qs, Qs0, Qs1, S).
  494frontier_link(E, [Q|Qs], Next, Qs1, S):-
  495	qp_update(E, Q, Q0, S),
  496	(	Q0 < 2
  497	->	Next = Qs0
  498	;	qp_memo_final(Q0, S),
  499		Next = [Q0|Qs0]
  500	),
  501	qp_update_succs(E-Q0, Q, S),
  502	frontier_link(E, Qs, Qs0, Qs1, S).
  503
  504% ?- zdd((qp_list(Q, [a]), qp_list(Q0, [b]), qp_list(Q1, [c]),
  505%	map_prod_sum([Q, Q0], [Q1], Z),
  506%	maplist(pred([A, M]:-qp_list(A, M)), Z, W))).
  507%
  508map_prod_sum(X, Y, Z):- map_prod_sum(X, Y, Z, []).
  509%
  510map_prod_sum(X, Y, Z, U):- @(map_prod_sum(X, Y, Z, U)).
  511
  512map_prod_sum([], _, U, U, _).
  513map_prod_sum([Q|Y], Z, U, V, S):- qp_list(Q, L, S),
  514	map_qp_sum(Z, L, U, U0, S),
  515	map_prod_sum(Y, Z, U0, V, S).
  516%
  517map_qp_sum([], _, U, U, _).
  518map_qp_sum([Q|As], L, [Q0|U], V, S):-qp_list(Q, M, S),
  519	ord_union(L, M, N),
  520	qp_list(Q0, N, S),
  521	map_qp_sum(As, L, U, V, S).
  522
  523%
  524prune_frontier(F, Link, F0):- @(prune_frontier(F, Link, F0)).
  525%
  526prune_frontier([], _, [], _).
  527prune_frontier([Q|F], Link, [Q|F0], S):- admissible_qp(Q, Link, S), !,
  528	prune_frontier(F, Link, F0, S).
  529prune_frontier([_|F], Link, F0, S):- prune_frontier(F, Link, F0, S).
  530
  531%
  532admissible_qp(Q, Link, S):- qp_list(Q, List, S),
  533	admissible_qp_list(List, Link).
  534
  535% ?- admissible_qp_list([p(0,0)-p(3, 1), p(3, 2)-p(3, 0)], p(0,0)-p(0,2)).
  536admissible_qp_list([], _):-!.
  537admissible_qp_list([X-X|List], Link):-!, admissible_qp_list(List, Link).
  538admissible_qp_list([_-X|_], P-_):- X @< P, !, fail.
  539admissible_qp_list([_|List], Link):- admissible_qp_list(List, Link).
  540
  541% ?- zdd((path_count(rect(1,0), C))).
  542% ?- zdd((path_count(rect(0,1), C))).
  543% ?- zdd((path_count(rect(1,1), C))).
  544% ?- zdd((path_count(rect(3,1), C))).
  545% ?- time(zdd((path_count(rect(1,8), C)))).   % take time.
  546%@ % 86,018,490 inferences, 10.604 CPU in 11.323
  547% ?- time(zdd((path_count(rect(8,1), C)))).
  548%@ % 9,254,518 inferences, 3.089 CPU in 4.850
  549% ?- zdd((path_count(rect(2,3), C))).
  550% ?- zdd((path_count(rect(3,2), C))).
  551% ?- time(zdd((path_count(rect(3,3), C)))).
  552%@ % 50,584,360 inferences, 7.362 CPU in 9.480
  553%@ C = 184 .
  554% ?- time(zdd((path_count(rect(4,4), C)))).
  555%@ % 9,474,429,287 inferences, 1256.779 CPU in 1433.734 seconds (88% CPU, 7538660 Lips)
  556%@ C = 8512 .
  557
  558qp_update(E, X, Y):- @(qp_update(E, X, Y)).
  559%
  560qp_update(_, X, 0, _):- X < 2.
  561qp_update(E, X, Y, S):- qp_list(X, U, S),
  562	memo(qp_update(E, X)-Y, S),
  563	(	nonvar(Y)  ->	true
  564	;	select_touch(E, [A, B], U, U0)
  565	->	connect_pairs(E, [A, B], W),
  566		ord_union([W], U0, V),
  567		qp_list(Y, V, S)
  568	;	Y = 0
  569	).
  570qp_update(_, _, 0,  _).
  571
  572% ?- zdd((qp_list(X, [a-a, b-b]), memo(qp_suc(X)-[]))).
  573qp_update_succs(G, Q):- @(qp_update_succs(G, Q)).
  574
  575% ?- zdd(qp_update_succs(a, 1)).
  576qp_update_succs(G, Q, S):- getmemo(qp_suc(Q)-U, V, S),
  577	G = E-Q0,
  578	(	var(V) -> U = [E-Q0]
  579	;	add_new(E-Q0, V, U)
  580	).
  581
  582qp_admissible(L):-
  583	forall((select(A-B, L, L0), select(C-D, L0, _)),
  584		   compatible_interval(A, B, C, D)).
  585
  586% ?- compatible_interval(1, 2, 3, 4).
  587% ?- compatible_interval(1, 3, 2, 4). % false
  588% ?- compatible_interval(2, 4, 5, 7).
  589% ?- compatible_interval(7, 5, 2, 4).
  590
  591% compatible_interval(I, J, X, Y):- disjoint_interval(I, J, X, Y), !.
  592% compatible_interval(I, J, X, Y):- sub_interval(I, J, X, Y), !.
  593% compatible_interval(I, J, X, Y):- sub_interval(X, Y, I, J).
  594%@ true.
  595%@ path-count4.pl qcompiled.
  596
  597% ?- select_touch(a-b, [A, B], [a-a, b-b], X).
  598% ?- select_touch(a-b, [A, B], [a-a, b-b], X),
  599%	connect_pairs(a-b, [A, B], Y).
  600%@ false.
  601
  602% ?- select_touch(a-b, [A, B], [a-a, b-c], X),
  603%	connect_pairs(a-b, [A, B], Y).
  604% ?- select_touch(a-b, [A, B], [a-a, c-b], X).
  605
  606select_touch(_, [], V, V).
  607select_touch(E, [A|As], [A|U], V):- touch(E, A), !,
  608	select_touch(E, As, U, V).
  609select_touch(E, As, [A|U], [A|V]):-select_touch(E, As, U, V).
  610
  611% ?- connect_pairs(a-b, [X-Y, U-V], W).
  612connect_pairs(X-Y, [U-X, Y-V], U-V).
  613connect_pairs(X-Y, [Y-V, U-X], U-V).
  614
  615% connect_pair(X-Y, U-X, U-Y).
  616% connect_pair(X-Y, Y-V, X-V).
  617
  618touch(A-_, _-A).
  619touch(_-A, A-_).
  620
  621% ?- add_new(a,[b,c,a], X).
  622add_new(X, Y, Y):- memberchk(X, Y), !.
  623add_new(X, Y, [X|Y]).
  624
  625% ?- zdd((add_child(a, 1), memo(a-X))).
  626% ?- zdd((add_child(a, 1), add_child(a, 2), memo(a-X))).
  627% ?- zdd((add_child(a, 1), add_child(a, 2), add_child(a, 1), memo(a-X))).
  628% ?- numlist(1, 100000, Ns),
  629%	time(zdd(( maplist(pred([Child]:- add_child(a, Child)), Ns), memo(a-X)))).
  630%@ % 2,600,319 inferences, 81.730 CPU in 81.817
  631
  632% ?- zdd((qp_list(Q, [p(0,0)-p(1,0)]), qp_list(Q, L))).
  633%@ L = [p(0, 0)-p(1, 0)].
  634%
  635qp_list(Q, List):- @(qp_list(Q, List)).
  636%
  637qp_list(0, 0, _):-!.
  638qp_list(1, [], _):-!.   

 ???

  639qp_list(Q, List, S):- var(Q), !,
  640	zdd_insert_atoms(List, 1, Q, S),
  641	memo(qp_suc(Q)-U, S),
  642	(	var(U) -> U = []
  643	;	true			% Already registered.
  644	).
  645qp_list(Q, [A|As], S):- cofact(Q, t(A, 0, R), S),
  646	qp_list(R, As, S).
  647
  648%
  649qp_memo_final(X):- @(qp_memo_final(X)).
  650%
  651qp_memo_final(Q, S):-  memoq(qp_final(Q)-true, S), !.
  652qp_memo_final(Q, S):- qp_list(Q, List, S),
  653	get_key(final, F, S),
  654	(	qp_final(List, F)
  655	->	memo(qp_final(Q)-true, S)
  656	;	true
  657	).
  658%
  659fold_qp_memo_final([], _).
  660fold_qp_memo_final([Q|Qs], S):- qp_memo_final(Q, S),
  661	fold_qp_memo_final(Qs, S).
  662
  663%
  664path_tree(Q, T):- @(path_tree(Q, T)).
  665%
  666
  667path_tree(X, 0,  _):- X<1, !.
  668path_tree(Q, 1, S):- memoq(qp_final(Q)-true, S), !.
  669path_tree(Q, T, S):- memo(path_tree(Q)-T, S),
  670	(	nonvar(T) -> true
  671	;	memo(qp_suc(Q) - L, S),		% writeln(qp_suc(Q)-L),
  672		path_tree_list(L, 0, T, S)
  673	).
  674%
  675path_tree_list([], C, C, _).
  676path_tree_list([E-Q|Qs], C0, C, S):-
  677	path_tree(Q, A, S),
  678	zdd_insert(E, A, A0, S),
  679	zdd_join(C0, A0, C1, S),
  680	path_tree_list(Qs, C1, C, S).
  681
  682% ?- qp_final([p(0,0)-p(2,3),p(0,3)-p(0,3),p(1,3)-p(1,3)], p(2,3)).
  683
  684qp_final([p(0,0)-F|R], F):-!, qp_all_id(R).
  685qp_final([A-A|R], F):- qp_final(R, F).
  686
  687%
  688qp_all_id([]).
  689qp_all_id([A-A|R]):- qp_all_id(R).
  690
  691%
  692dag_card(Q, C):- @(dag_card(Q, C)).
  693%
  694dag_card(Q, C, S):- open_state(M),
  695	dag_card(Q, C, S, M),
  696	close_state(M).
  697%
  698dag_card(Q, 0,  _, _):- Q<2, !.
  699dag_card(Q, 1, S, _):- memoq(qp_final(Q)-true, S), !.
  700dag_card(Q, C, S, M):- memo(dag_card(Q)-C, M),
  701	(	nonvar(C) -> true
  702	;	memo(qp_suc(Q) - L, S),
  703		dag_card_list(L, 0, C, S, M)
  704	).
  705%
  706dag_card_list([], C, C, _, _).
  707dag_card_list([_-Q|Qs], C0, C, S, M):-
  708	dag_card(Q, A, S, M),
  709	C1 is C0 + A,
  710	dag_card_list(Qs, C1, C, S, M).
  711
  712
  713		/***************
  714		*     debug    *
  715		***************/
  716%
  717snap_qp_update(E, X, Y, S):- qp_update(E, X, Y, S),
  718	qp_list(X, X0, S),
  719	qp_list(Y, Y0, S),
  720	format("~w\n~w\n~w\n\n", [qp_update(E, X, Y), X0, Y0]).
  721%
  722dump_frontier(Qs, S):- maplist(pred(S,
  723	([Q]:- qp_list(Q, U, S), writeln(U))), Qs).
  724%
  725dump_qp(X, _):-  X<2, !, writeln(X).
  726dump_qp(X, S):-  qp_list(X, L, S),
  727	writeln(X=L).
  728
  729%
  730delta_star([], Q, Q, _).
  731delta_star([E|Es], Q,  Q0, S):-  memo(qp_suc(Q)-L, S),
  732	writeln(E),
  733	memberchk(E-Q1, L),
  734	dump_qp(Q1, S),
  735	delta_star(Es, Q1, Q0, S).
  736
  737%
  738mirror(X, Y):- @(mirror(X, Y)).
  739%
  740mirror(X, X, _):- X<2, !.
  741mirror(X, Y, S):- cofact(X, t(p(I,J)-p(K,L), Lx, Rx), S),
  742	mirror(Lx, MLx, S),
  743	mirror(Rx, MRx, S),
  744	zdd_insert(p(J, I)-p(L, K), MRx, MRx0, S),
  745	zdd_join(MLx, MRx0, Y, S).
  746
  747%
  748rect_path_tree(Rect, T):-@(rect_path_tree(Rect, T)).
  749
  750rect_path_tree(Rect, T, S):-
  751	sort_links_by_rank(Rect, Rs, Ids),
  752	Rect = rect(W, H),
  753	set_key(final, p(W, H), S),
  754	qp_list(Q, Ids, S),
  755	repeat_frontier(Rs, [Q], FinalQs, S),
  756	fold_qp_memo_final(FinalQs, S),
  757	path_tree(Q, T, S).
  758
  759% ?- spy(remove_end_links).
  760% ?- remove_end_links([1-[c-a]], R0, a-d).
  761% ?- remove_end_links([1-[c-a, d-c, a-d]], R0, a-d).
  762
  763remove_end_links(R, R0, ST):- maplist(pred(ST,
  764	(	[J-Es, J-Es0]:- foldl(
  765		pred([ST, Es, Es0],
  766			 ( [P-Q, U, V]:- ST = S-T,
  767					(	( Q==S; P==T)
  768					->  V = U
  769					;	U = [P-Q|V]
  770					))),  Es,  Es0, []))),  R, R0).
  771
  772%
  773rest(Q, X, Y, S):- memo(qp_suc(Q)-L, S),
  774	(	select(Q, X, X0)-> true
  775	;	X0 = X
  776	),
  777	rest_list(L, X0, Y, S).
  778%
  779rest_list([], X, X, _).
  780rest_list([_-Q|Qs], X, Y, S):-
  781	rest(Q, X, X0, S),
  782	rest_list(Qs, X0, Y, S).
  783
  784
  785% ?- zdd((dbg_qp_update(p(0,0)-p(1,0), [p(0,0)-p(0,0),p(1,0)-p(1,0)], L))).
  786% ?- spy(dbg_qp_update).
  787% ?- zdd((dbg_qp_update(p(2,2)-p(2,3), [p(0,0)-p(2,2),p(0,3)-p(0,3), p(1,3)-p(1,3), p(2,3)-p(2,3)], L))).
  788
  789dbg_qp_update(E, X, Y):- @(dbg_qp_update(E, X, Y)).
  790%
  791dbg_qp_update(E, X, Y, S):- qp_list(Q, X, S),
  792	qp_update(E, Q, Q0, S),
  793	qp_list(Q0, Y, S).
  794
  795
  796% ?- zdd((rect_path_tree(rect(3,2), T), card(T, C),
  797%   zdd_insert_atoms([p(0,0)-p(2,3),p(0,3)-p(0,3),p(1,3)-p(1,3)], 1, W),
  798%	mirror(T, MT), card(MT, CMT),
  799%	rect_path_tree(rect(2,3), T0), card(T0, C0), mirror(T0, MT0),
  800%	card(MT0, CM0),
  801%	zdd_subtr(T, MT0, D), card(D, CD), psa(D)
  802%)).
  803
  804
  805% Test to show getmemo/2,3 works well.
  806% ?- time((numlist(1, 10000, Ns),
  807%	zdd((maplist( pred([A]:- getmemo(a(A)-A, _)), Ns),
  808%		 maplist( pred(([A]:- getmemo(a(A)-B, Old), B is 2*Old)),  Ns),
  809%		 foldl( pred(([I, B, C]:- getmemo(a(I)-J, J),
  810%							 C is B + J)), Ns, 0, Sum))))).
  811%@ % 938,719 inferences, 0.074 CPU in 0.075 seconds (99% CPU, 12701698 Lips)
  812%@ Ns = [1, 2, 3, 4, 5, 6, 7, 8, 9|...],
  813%@ Sum = 100010000.