1:- module(zdd_euler, []).    2:- use_module(zdd('zdd-array')).    3:- use_module(zdd(math)).    4
    5% ?- [misc(zdd)], module(zdd_euler).
    6
    7% ?- open_array(A), new_elem(hello, A, I), get_elem(I, A, V), set_elem(I, A, world), get_elem(I, A, W).
    8% ?- open_array(A), hash_memo(f(a), A, V), V=a(hello), setarg(1, V, world), hash_memo(f(a), A, W).
    9% ?- derive_orbits([[p(0,0)]], Os, 1).
   10% ?- derive_orbits([[p(0,0)]], Os, 2).
   11% ?- derive_orbits([[p(0,0)]], Os, 3).
   12% ?- derive_orbits([[p(0,0)], [p(0,0), p(0, 1)]], Os, 3).
   13% ?- derive_orbits([[p(0,0), p(1,0)], [p(0,0), p(0, 1)]], Os, 3).
   14% ?- derive_orbits([[p(1, 1), p(2, 1)]], R, 2).
   15
   16derive_orbits(Ps, Orbits, W):-
   17	maplist(orbit_with_bound, Ps, Bs),
   18	append(Bs, Bs0),
   19	shift_variants_list(Bs0, Cs, h(W)),
   20	append(Cs, Ds),
   21	sort(Ds, Orbits).
   22%
   23shift_variants_list([], [], _).
   24shift_variants_list([X|Xs], [Y|Ys], W):-
   25	all_shift_variant(X, W, Y),
   26	shift_variants_list(Xs, Ys, W).
   27
   28%
   29prepare_euler(W, Ps):-
   30		derive_orbits(Ps, Orbits, W),
   31		W0 is W - 1,
   32		findall(p(I,0), between(0, W0, I), Floor),
   33		b_setval(euler, e(Orbits, W, Floor)).
   34
   35% e(Orbits, Width, Floor, DG).
   36euler(E)	:- b_getval(euler, E).
   37orbits(X)	:- b_getval(euler, E), arg(1, E, X).
   38width(X)	:- b_getval(euler, E), arg(2, E, X).
   39floor(X)	:- b_getval(euler, E), arg(3, E, X).
   40
   41% basic access to the array.
   42euler_memo(X, I, E):- zdd_memo(X-I, zdd_euler),
   43	  ( var(I)
   44	  ->	b_getval(zdd_euler, R),
   45			new_elem(E, R, I)
   46	  ;		euler_get(I, E)
   47	  ).
   48%
   49euler_memo(X, I):- zdd_memo(X-I, zdd_euler).
   50%
   51euler_get(I, E):- b_getval(zdd_euler, R),
   52	  arg(3, R, Array),
   53	  arg(I, Array, E).
   54
   55		/**************
   56		*     tiny    *
   57		**************/
   58%
   59cons(X, Y, [X|Y]).
   60%
   61zip([],[],[]).
   62zip([A|As],[B|Bs],[A-B|Cs]):- zip(As, Bs, Cs).
   63%
   64disjoint([], _):-!.
   65disjoint([A|As], X):- \+ memberchk(A, X), disjoint(As, X).
   66
   67% ?- max(1, 2, X).
   68max(X, Y, M):- (X > Y -> M = X; M = Y).
   69min(X, Y, M):- (X < Y -> M = X; M = Y).
   70
   71% ?- euler_tile_count(p(2, 2), [[p(0,0),p(1,0)]], H).
   72
   73% ?-time(zdd_euler:euler_tile_count(p(14, 14), [[p(0,0),p(1,0)]], H)).
   74%@ % 13,781,361 inferences, 48.576 CPU in 48.716 seconds (100% CPU, 283706 Lips)
   75%@ H = 112202208776036178000000 .
   76
   77%@ % 13,723,536 inferences, 47.405 CPU in 47.534 seconds (100% CPU, 289495 Lips)
   78%@ H = 112202208776036178000000 .
   79
   80euler_tile_count(p(U, V), Ps, Count):-
   81		open_zdd(zdd_euler, 2),
   82%		b_setval(zdd_euler, #(0, #(2, #([],[])), #(_,_))),
   83		(	U > V -> X = V, Y = U
   84		;	X = U, Y = V
   85		),
   86		prepare_euler(X, Ps),
   87		euler_memo([], I, E),
   88		E = []-_,
   89		euler_dg([I]),
   90		euler_count(Y, I, Count).
   91
   92% ?- euler_tile_count(p(1, 1), [[p(0,0)]], H).
   93euler_count(Y, I, Count):-
   94		unfold_dg(H),
   95		b_getval(zdd_euler, #(N,_,_)),
   96		functor(Cs, #, N),
   97		initial_args(Cs, 0),
   98		setarg(I, Cs, 1),
   99		repeat_matrix(H^Y, Cs, Ds),
  100		arg(I, Ds, Count).
  101
  102% ?- repeat_matrix(h([1-2,2-2],[1-3,2-3])^3, g(1,0), R).
  103% ?- repeat_matrix(h([1-2,2-2],[1-3,2-3]), g(1,0), R).
  104% ?- repeat_matrix(h([1-2,2-2],[1-3,2-3])^10, g(1,0), R).
  105repeat_matrix(_^0, Cs, Cs):-!.
  106repeat_matrix(H^1, Cs, Ds):-!,
  107	repeat_matrix(H, Cs, Ds).
  108repeat_matrix(H^I, Cs, Ds):-!,
  109	repeat_matrix(H, Cs, Cs0),
  110	I0 is I - 1,
  111	repeat_matrix(H^I0, Cs0, Ds).
  112repeat_matrix(H, Cs, Ds):-functor(Cs, F, N),
  113	functor(Ds, F, N),
  114	repeat_matrix(N, H, Cs, Ds).
  115%
  116repeat_matrix(0, _, _, _):-!.
  117repeat_matrix(I, H, Cs, Ds):-
  118		arg(I, H, Row),
  119		linear_sum(Row, Cs, 0, P),
  120		arg(I, Ds, P),
  121		I0 is I - 1,
  122		repeat_matrix(I0, H, Cs, Ds).
  123%
  124linear_sum([], _, P, P).
  125linear_sum([I-A|M], Cs, P, Q):-
  126		arg(I, Cs, V),
  127		P0 is P + A*V,
  128		linear_sum(M, Cs, P0, Q).
  129
  130% euler_dg/1.
  131euler_dg([]).
  132euler_dg([I|R]):- euler_get(I, E), !,
  133		E = A-S,
  134		(	nonvar(S)	->	R0 = R
  135		;	orbits(Orbits),
  136			euler_dg(Orbits, A, S0, []),
  137			euler_dg(S0, S),
  138			union(S, R, R0)
  139		),
  140		euler_dg(R0).
  141% euler_dg/2.
  142euler_dg([], []).
  143euler_dg([X|Xs], [I|Is]):-
  144		euler_memo(X, I),
  145		(	nonvar(I) -> true
  146		;	euler_memo(X, I, X-_)
  147		),
  148		euler_dg(Xs, Is).
  149% euler_dg/4.
  150euler_dg(Orbits, A, R, R0):- width(W),
  151		W0 is W - 1,
  152		between(0, W0, K),
  153		\+ memberchk(p(K,0), A),
  154		!,
  155		euler_dg(Orbits, A, p(K,0), R, R0).
  156euler_dg(_, A, [B|R], R):-
  157		cursor_up(A, B).
  158% euler_dg/5.
  159euler_dg([], _,  _, L, L).
  160euler_dg([Orbit|R], A, P, [B0|L], L0):- memberchk(P, Orbit),
  161		disjoint(A, Orbit), !,
  162		union(A, Orbit, B),
  163		sort(B, B0),
  164		euler_dg(R, A, P, L, L0).
  165euler_dg([_|R], U, P, L, L0):-
  166		euler_dg(R, U, P, L, L0).
  167
  168%
  169unfold_dg(H):- b_getval(zdd_euler, #(N,_,G)),
  170	functor(H, #, N),
  171	floor(Floor),
  172	unfold_dg(N, H, G, Floor).
  173%
  174unfold_dg(0, _H, _G, _Floor):- !.
  175unfold_dg(I, H, G, Floor):- unfold_dg(I, _Bag, H, G, Floor),
  176	I0 is I - 1,
  177	unfold_dg(I0, H, G, Floor).
  178%
  179unfold_dg(I, Bag, H, G, Floor):- arg(I, H, U),
  180	(	nonvar(U) -> Bag = U
  181	;	arg(I, G, A-L),
  182		(	subset(Floor, A) ->  maplist(bagfy, L, Bag)
  183		;	unfold_dg(L, [], Bag, H, G, Floor)
  184		),
  185		setarg(I, H, Bag)
  186	).
  187%
  188unfold_dg([], Bag, Bag, _, _, _).
  189unfold_dg([I|Is], Bag, OutBag, H, G, Floor):- unfold_dg(I, B, H, G, Floor),
  190		sum_row_row(B, Bag, Bag1),
  191		unfold_dg(Is, Bag1, OutBag, H, G, Floor).
  192%
  193bagfy(X, X-1).
  194
  195% ?- sum_row_row([1-1,2-2], [1-2,3-1], X).
  196sum_row_row([], X, X).
  197sum_row_row(X, [], X).
  198sum_row_row([A-I|R], [A-J|S], [A-K|T]):- !, K is I + J, sum_row_row(R, S, T).
  199sum_row_row([A-I|R], [B-J|S], [A-I|T]):- A<B, !, sum_row_row(R, [B-J|S], T).
  200sum_row_row(X, [B-J|S], [B-J|T]):- sum_row_row(X, S, T).
  201
  202% ?- sum_row_row(2, 3, [1-1, 2-2], [1-1, 2-2], X).
  203sum_row_row(P, Q, A, B, C):- mul_scalar_row(P, A, A0),
  204	mul_scalar_row(Q, B, B0),
  205	sum_row_row(A0, B0, C).
  206
  207%
  208initial_args(H, A):- functor(H, _, N),
  209				   initial_args(N, H, A).
  210%
  211initial_args(0, _, _).
  212initial_args(I, H, A):- arg(I, H, A),
  213					I0 is I - 1,
  214					initial_args(I0, H, A).
  215
  216%
  217cursor_up(X, Y):- shift_down(X, X0), sort(X0, Y).
  218
  219% ?- orbit_with_bound([p(1,1)], F).
  220% ?- orbit_with_bound([p(0,0), p(1,0)], F).
  221% ?- orbit_with_bound([p(0,0), p(1,0), p(0, 1)], F).
  222% ?- orbit_with_bound([p(-1,-1), p(-1,0), p(0, -1)], F), length(F, N).
  223% ?- orbit_with_bound([p(1,1), p(2,2)], F).
  224% ?- orbit_with_bound([p(1,1)], F).
  225
  226orbit_with_bound(X, Y):- apply_mirror_rotate(X, L),
  227	maplist(normalize_shape, L, L0),
  228	sort(L0, Y).
  229
  230% ?- all_shift_variant([p(0,0)]-p(1,1), p(2,2), R0).
  231% ?- all_shift_variant([p(0,0)]-p(1,1), v(2), R0).
  232% ?- all_shift_variant([p(0,0)]-p(1,1), h(2), R0).
  233% ?- shift_variant([p(0,0)]-p(1,1), p(2,2), R0).
  234% ?- shift_variant([p(0,0)]-p(1,1), v(2), R0).
  235% ?- shift_variant([p(0,0)]-p(1,1), h(2), R0).
  236all_shift_variant(U, B, Rs):-
  237	findall(R, shift_variant(U, B, R), Rs).
  238%
  239shift_variant(B, h(I), R0):-!, shift_variant_hori(B, I, R0).
  240shift_variant(B, v(I), R0):-!, shift_variant_vert(B, I, R0).
  241shift_variant(R-p(X,Y), p(I,J), R0):- H is I-X,
  242	V is J-Y,
  243	between(0, H, A),
  244	between(0, V, B),
  245	move_region(p(A, B), R, R0).
  246%
  247shift_variant_hori(R-p(X,_), I, R0):- H is I-X,
  248	between(0, H, A),
  249	move_region(p(A, 0), R, R0).
  250%
  251shift_variant_vert(R-p(_,Y), I, R0):- V is I-Y,
  252	between(0, V, A),
  253	move_region(p(0, A), R, R0).
  254%
  255apply_mirror_rotate(X, [X, M, R1, R2, R3]):- mirror(X, M),
  256	rotate(X, R1),
  257	rotate(R1, R2),
  258	rotate(R2, R3).
  259
  260% ?- normalize_shape([p(1,1)], R).
  261% ?- normalize_shape([p(1,1),p(1,2)], R).
  262% ?- normalize_shape([p(2,2),p(3,3)], R).
  263% ?- normalize_shape([p(-2,-2),p(0,0)], R).
  264normalize_shape(X, Y-p(J,K)):- boundary_rect(X, p(A,B)-p(C,D)),
  265		H is -A,
  266		V is -B,
  267		move_region(p(H,V), X, Y0),
  268		sort(Y0, Y),
  269		J is C + H + 1,
  270		K is D + V + 1.
  271
  272% ?- move_region(p(2,3),[p(1,1)], R).
  273move_region(_, [], []).
  274move_region(p(H,V), [p(X,Y)|R], [p(X0, Y0)|R0]):-
  275	X0 is X + H,
  276	Y0 is Y + V,
  277	move_region(p(H,V), R, R0).
  278
  279% Generic subset.
  280% ?- sub_region([p(1,1)], p(1,0)-p(2,3)).
  281% ?- sub_region([p(1,1)], p(2,3)).
  282% ?- sub_region([p(1,1)], v(0, 2)).
  283% ?- sub_region([p(1,1)], h(0, 2)).
  284% ?- sub_region([p(1,1)], [p(2,2), p(1,1)]).
  285sub_region(R,  p(I,J)-p(I0, J0)):- !, sub_rect(R, I, I0, J, J0).
  286sub_region(R,  p(I, J)):- !, sub_rect(R, 0, I, 0, J).
  287sub_region(R,  v(I, J)):- !, sub_vert_band(R, I, J).
  288sub_region(R,  v(J)):- !, sub_vert_band(R, 0, J).
  289sub_region(R,  h(I, J)):- !, sub_hor_band(R, I, J).
  290sub_region(R,  h(J)):- !, sub_hor_band(R, 0, J).
  291sub_region(R,  L):- subset(R, L).
  292%
  293sub_rect([], _, _, _, _).
  294sub_rect([p(X,Y)|R], I, I0, J, J0):- I=<X, X=<I0, J=<Y, Y=<J0,
  295	sub_rect(R, I, I0, J, J0).
  296%
  297sub_vert_band([], _, _).
  298sub_vert_band([p(X,_)|R], I, I0):- I=<X, X=<I0,
  299	sub_vert_band(R, I, I0).
  300%
  301sub_hor_band([], _, _).
  302sub_hor_band([p(_,Y)|R], J, J0):- J=<Y, Y=<J0,
  303	sub_hor_band(R, J, J0).
  304
  305% ?- boundary_rect([p(1,1)], R).
  306% ?- boundary_rect([p(1,1),p(1,2), p(2,1)], R).
  307boundary_rect(S, R):- S = [P|Ps], boundary_rect(Ps, P-P, R).
  308%
  309boundary_rect([], X, X).
  310boundary_rect([P|Ps], X, Y):-  boundary_rect_one(P, X, X0),
  311	boundary_rect(Ps, X0, Y).
  312%
  313boundary_rect_one(P, X-Y, U-V):- choose_low(P, X, U), choose_high(P, Y, V).
  314%
  315choose_low(p(I,J), p(A,B), p(U,V)):- min(I, A, U), min(J, B, V).
  316%
  317choose_high(p(I,J), p(A,B), p(U,V)):- max(I, A, U), max(J, B, V).
  318%
  319mirror([], []).
  320mirror([p(I,J)|R], [p(J, I)|R0]):- mirror(R, R0).
  321%
  322rotate([], []).
  323rotate([p(I,J)|R], [p(J0,I)|R0]):- J0 is -J, rotate(R, R0).
  324%
  325boundary_chk([], _).
  326boundary_chk([p(X,_)|R], Width):- X>=0, X<Width,
  327		boundary_chk(R, Width).
  328%
  329shift_horizontal(X, 0, X):-!.
  330shift_horizontal([], _, []).
  331shift_horizontal([p(X, Y)|R], I, [p(X0, Y)|R0]):-
  332		X0 is I + X,
  333		shift_horizontal(R, I, R0).
  334%
  335shift_down([], []).
  336shift_down([p(_,0)|R], S):-!, shift_down(R, S).
  337shift_down([p(X,I)|R], [p(X, I0)|S]):- I0 is I - 1,
  338	   shift_down(R, S).
  339
  340		/***************************
  341		*     Matrix operations    *
  342		***************************/
  343
  344pred_rows(U, X, Pred):-
  345	( integer(U) -> N = U, F = (#)
  346	; functor(U, F, N)
  347	),
  348	functor(X, F, N),
  349	pred_arg_row(N, Pred, X).
  350%
  351pred_arg_row(0, _, _):-!.
  352pred_arg_row(I, Pred, X):-	call(Pred, I, V), !,
  353	arg(I, X, V),
  354	I0 is I - 1,
  355	pred_arg_row(I0, Pred, X).
?- matrix(unit(10), B), writeln(B). ?- matrix(unit(10)^100000000000000000000, B), writeln(B). ?- matrix((2*unit(10))^100, B), writeln(B). ?- matrix(dim(10, unit), B), writeln(B). ?- matrix(dim(10, step), B), writeln(B). ?- matrix(dim(10, zero), B), writeln(B). ?- matrix(dim(10, nil_potent(10)), R), writeln(R).
  365unit(I, [I-1]).
  366zero(_, []).
  367step(I, [I-I]).
  368nil_potent(N, N, []):-!.
  369nil_potent(_, I, [J-1]):- J is I + 1.
  370%
  371unit_matrix(X, Y):- pred_rows(X, Y, unit).
  372
  373% ?- matrix(1 + (a([1-1],[2-1])-a([1-1],[2-1]))^1, R), writeln(R).
  374% ?- matrix(a([1-1],[2-1])^0, R), writeln(R).
  375% ?- matrix(-1 + a([1-1],[2-1]), R), writeln(R).
  376% ?- N=2*100+1, matrix((-1 * a([1-1],[2-1]))^N, R), writeln(R).
  377% ?- matrix(dim(10, step)^3, R), writeln(R).
  378% ?- matrix(dim(10, nil_potent(10))^3, R), writeln(R).
  379% ?- matrix(unit(10)^3, R), writeln(R).
  380% ?- matrix(tr(unit(10))^3, R), writeln(R).
  381% ?- matrix(tr(dim(10, nil_potent(10)))^10, R), writeln(R).
  382% ?- matrix(tr(dim(10, nil_potent(10)))^5, R), writeln(R).
  383
  384matrix(A + B, C):-!, matrix(A, A0),
  385	matrix(B, B0),
  386	matrix(+, A0, B0, C).
  387matrix(A - B, C):-!, matrix(A + (-1)*B, C).
  388matrix(-(A), C):-!, matrix((-1) * A, C).
  389matrix(A * B, C):-!, matrix(A, A0),
  390	matrix(B, B0),
  391	matrix(*, A0, B0, C).
  392matrix(A ^ B, C):-!, B0 is B,
  393	matrix(A, A0),
  394	matrix(^, A0, B0, C).
  395matrix(unit(A), B):-!, pred_rows(A, B, unit).
  396matrix(dim(A, Pred), B):-!, pred_rows(A, B, Pred).
  397matrix(tr(A), B):-!, matrix(A, A0), transpose_matrix(A0, B).
  398matrix(A, A).
  399
  400%
  401matrix(F, A, B, C):- atomic(A), atomic(B), !,
  402	(	F == (+) -> C is A + B
  403	;   F == (*) -> C is A * B
  404	;	C is A ^ B
  405	).
  406matrix(*, A, B, C):-!,
  407	(	number(A) -> mul_scalar_matrix(A, B, C)
  408	;	number(B) -> mul_scalar_matrix(B, A, C)
  409	;	mul_matrix(A, B, C)
  410	).
  411matrix(+, A, B, C):-!,
  412	(	number(A) -> functor(B,_,N),
  413					 unit_matrix(N, U),
  414					 mul_scalar_matrix(A, U, V),
  415					 sum_matrix(V, B, C)
  416	;	number(B) -> functor(A,_,N),
  417					 unit_matrix(N, U),
  418					 mul_scalar_matrix(B, U, V),
  419					 sum_matrix(V, A, C)
  420	;	sum_matrix(A, B, C)
  421	).
  422matrix(^, A, N, C):-
  423	(	number(A) -> C is A**N
  424	;	pow_matrix(A, N, C)
  425	).
  426
  427%
  428mul_scalar_matrix(1, X, X):-!.
  429mul_scalar_matrix(0, X, Y):-!, functor(X,F,N),functor(Y,F,N),
  430	initial_args(N, Y, []).
  431mul_scalar_matrix(A, X, Y):- functor(X, F, N), functor(Y, F, N),
  432	mul_scalar_matrix(N, A, X, Y).
  433%
  434mul_scalar_matrix(0, _, _, _):-!.
  435mul_scalar_matrix(I, A, X, Y):- arg(I, X, R),
  436	 mul_scalar_row(A, R, S),
  437	 arg(I, Y, S),
  438	 I0 is I - 1,
  439	 mul_scalar_matrix(I0, A, X, Y).
  440%
  441sum_matrix(X, Y, Z):- functor(X, F, N), functor(Z, F, N),
  442	 sum_matrix(N, X, Y, Z).
  443%
  444sum_matrix(0, _, _, _):-!.
  445sum_matrix(I, X, Y, Z):- arg(I, X, R),
  446	 arg(I, Y, S),
  447	 sum_row_row(R, S, T),
  448	 arg(I, Z, T),
  449	 I0 is I - 1,
  450	 sum_matrix(I0, X, Y, Z).
  451
  452% ?- mul_matrix(#([2- (-1)], [1-1]), #([2- (-1)], [1-1]), R), writeln(R).
  453mul_matrix(X, Y, Z):- functor(X, F, N), functor(Z, F, N),
  454	 mul_matrix(N, X, Y, Z).
  455%
  456mul_matrix(0, _, _, _):-!.
  457mul_matrix(I, X, Y, Z):- arg(I, X, R),
  458	  mul_row_matrix(R, Y, [], U),
  459	  arg(I, Z, U),
  460	  I0 is I - 1,
  461	  mul_matrix(I0, X, Y, Z).
  462%
  463mul_row_matrix([], _, U, U).
  464mul_row_matrix([J-A|R], X, U, V):- arg(J, X, B),
  465	mul_scalar_row(A, B, B0),
  466	sum_row_row(B0, U, U0),
  467	mul_row_matrix(R, X, U0, V).
  468%
  469mul_scalar_row(_, [], []):-!.
  470mul_scalar_row(A, [I-V|R], [I-U|S]):- U is A*V,
  471	mul_scalar_row(A, R, S).
  472%
  473pow_matrix(A, 0, C):-!, functor(A,_,N), unit_matrix(N, C).
  474pow_matrix(A, 1, A):-!.
  475pow_matrix(A, N, B):- divmod(N, 2, K, M),
  476	pow_matrix(A, K, B0),
  477	matrix(*, B0, B0, B1),
  478	(	M =	0 -> B = B1
  479	;  	matrix(*, A, B1, B)
  480	).
  481
  482% ?- transpose_matrix(a([2-1], [1-2]), B), transpose_matrix(B, C),
  483%	transpose_matrix(C, D).
  484transpose_matrix(X, Y):- number(X), !, Y=X.
  485transpose_matrix(X, Y):- functor(X, F, N), functor(Y, F, N),
  486	initial_args(N, Y, []),
  487	transpose_matrix(N, X, Y).
  488%
  489transpose_matrix(0, _, _):-!.
  490transpose_matrix(I, X, Y):- arg(I, X, Row),
  491	transpose_row(Row, I, Y),
  492	I0 is I - 1,
  493	transpose_matrix(I0, X, Y).
  494%
  495transpose_row([], _, _).
  496transpose_row([J-V|R], I, X):-
  497	arg(J, X, U),
  498	setarg(J, X, [I-V|U]),
  499	transpose_row(R, I, X).
  500
  501% ?- column(0, a([0-1],[1-1],[0-2,3-1]), C).
  502column(I, M, C):- functor(M, F, N),
  503	functor(C, F, N),
  504	column(N, I, M, C).
  505
  506column(0, _I, _M, _C):-!.
  507column(J, I, M, C):- arg(J, M, Row),
  508	(	ord_memberchk(I, V, Row) -> true
  509	;	V = 0
  510	),
  511	arg(J, C, V),
  512	J0 is J - 1,
  513	column(J0, I, M, C).
  514%
  515ord_memberchk(I, V, [I-V|_]):-!.
  516ord_memberchk(I, V, [J-_|R]):- J<I,
  517	ord_memberchk(I, V, R).
  518
  519%
  520swap_row(I, J, A):- arg(I, A, U),
  521					arg(J, A, V),
  522					setarg(I, A, V),
  523					setarg(J, A, U).
  524
  525% ?- A = a([2-1], [1-2]), swap_column(1,2, A).
  526swap_column(I, J, A):- functor(A,_,N),
  527	(	I=<J -> swap_column(N, I, J, A)
  528	;	swap_column(N, J, I, A)
  529	).
  530swap_column(0, _, _, _):-!.
  531swap_column(K, I, J, A):- arg(K, A, L),
  532	swap_key(L, L0, I, J),
  533	setarg(K, A, L0),
  534	K0 is K - 1,
  535	swap_column(K0, I, J, A).
  536
  537% ?- swap_key(1,2,[1-a,2-b], Y).
  538% ?- swap_key(1,2,[], Y).
  539% ?- swap_key(1,3,[1-a, 2-b, 3-c], Y).
  540% ?- swap_key(1,3,[2-b, 3-c], Y).
  541% ?- swap_key(1,3,[2-b], Y).
  542% ?- swap_key(2,3,[2-b], Y).
  543% Swap keys in an assoc list.
  544swap_key(X, Y, I, J):- swap_key(I-_, J-_, X, Y, []).
  545%
  546swap_key(_, _, [], X, X):-!.
  547swap_key(I-A, J-B, [I-A|S], [I-B|X], Y):-!,
  548	swap_key_sec(A, J-B, S, X, Y).
  549swap_key(I-A, J-B, [K-C|S], [K-C|X], Y):- K<I, !,
  550	swap_key(I-A, J-B, S, X, Y).
  551swap_key(I-A, J-B, S, [I-B|X], Y):-!, A = 0,
  552	swap_key_sec(A, J-B, S, X, Y).
  553%
  554swap_key_sec(A, J-0, [], [J-A|X], X):-!.
  555swap_key_sec(A, J-B, [J-B|S], [J-A|X], Y):-!, append(S, Y, X).
  556swap_key_sec(A, J-B, [K-C|S], [K-C|X], Y):- K<J, !,
  557	swap_key_sec(A, J-B, S, X, Y).
  558swap_key_sec(A, J-0, S, [J-A|X], Y):- append(S, Y, X).
  559
  560% Find a pivot.
  561% ?- find_pivot(a([1-1], [2-1], [3-1]), 3, 1, 1, I).
  562% ?- find_pivot(a([1-1,2-2], [2-1], [3-1]), 3, 2, 1, I).
  563% ?- find_pivot(a([1-1], [2-1], [3-1]), 3, 2, 2, I).
  564% ?- find_pivot(a([1-1], [3-1], [2-1]), 3, 2, 2, I).
  565find_pivot(Mat, Dim, I, J, K):-
  566	between(I, Dim, J),
  567	arg(J, Mat, Row),
  568	between(I, Dim, K),
  569	memberchk(K-V, Row),
  570	V \== 0.
  571
  572% ?- M=a([1-1]), sum_row_to_arg(1, M, [1-1]).
  573% ?- M=a([1-1]), sum_row_to_arg(1, M, [1-1, 2-2]).
  574sum_row_to_arg(I, M, A):- arg(I, M, B),
  575	sum_row_row(A, B, C),
  576	setarg(I, M, C).
  577
  578% ?- M=a([1-1]), sum_row_to_arg(1, M, [1-1, 2-2], -1).
  579sum_row_to_arg(I, M, A, U):- mul_scalar_row(U, A, B),
  580	sum_row_to_arg(I, M, B).
  581
  582% ?- M = a([0-16, 1-2, 2-4], [0-5, 1 - 1, 2-1]),  gauss(M).
  583gauss(M):- functor(M,_,Dim),
  584		   gauss(1, M, Dim).
  585
  586% ?- M = a([1-1]),  gauss(1, M, 1).
  587gauss(I, _, Dim):- I > Dim, !.
  588gauss(I, M, Dim):- find_pivot(M, Dim, I, J, K), !,
  589		(	K \== I	->	swap_row(I, K, M)	; true),
  590		(	J \== I	->	swap_column(I, J, M); true),
  591		arg(I, M, Row),
  592		memberchk(I-A, Row),
  593		A0 is 1 rdiv A,
  594		mul_scalar_row(A0, Row, Row0),
  595		setarg(I, M, Row0),
  596		sweep(Dim, I, Row0, M),
  597		I0 is I + 1,
  598		gauss(I0, M, Dim).
  599gauss(_,_,_).
  600
  601%
  602sweep(0, _I, _Row, _M):-!.
  603sweep(I, I, Row, M):-!, K is I - 1,
  604	sweep(K, I, Row, M).
  605sweep(K, I, Row, M):- arg(K, M, KRow),
  606	( memberchk(I-V, KRow), V \== 0
  607	->	V0 is -V,
  608		sum_row_to_arg(K, M, Row, V0)
  609	;	true
  610	),
  611	K0 is K - 1,
  612	sweep(K0, I, Row, M).
  613
  614% Gauss method with history.
  615% ?- M = a([1-1]),  gauss(1, M, 1, P, []).
  616% ?- M = a([1-1], [2-2]),  gauss(1, M, 2, P, []).
  617% ?- M = a([], [2-2]),  gauss(1, M, 2, P, []).
  618% ?- M = a([2- -1], [1-1]),  gauss(1, M, 2, P, []).
  619% ?- M = a([2- -1], [1-1]),  gauss(1, M, 2, P, []).
  620% ?- M = a([0-1, 1-1], [0-2,2-1]),  gauss(2, M, 2, P, []).
  621% ?- M = a([0-1, 2 - -1], [0-2, 1 - 1]),  gauss(M, P).
  622% ?- M = a([0-16, 1-2, 2-4], [0-5, 1 - 1, 2-1]),  gauss(M).
  623
  624gauss(Matrix, History):- functor(Matrix, _, Dim),
  625			  gauss(1, Matrix, Dim, History, []).
  626%
  627gauss(I, _, Dim, P, P):- I > Dim, !.
  628gauss(I, M, Dim, P, Q):- find_pivot(M, Dim, I, J, K), !,
  629		(	K \== I
  630		->	swap_row(I, K, M),
  631			P = [swap_row(I,K)|P1]
  632		;	P1 = P
  633		),
  634		(	J \== I
  635		->	swap_column(I, J, M),
  636			P1 = [swap_col(I,J)|P2]
  637		;	P2 = P1
  638		),
  639		arg(I, M, Row),
  640		memberchk(I-A, Row),
  641		A0 is 1 rdiv A,
  642		mul_scalar_row(A0, Row, Row0),
  643		setarg(I, M, Row0),
  644		P2 = [mul_scalar(I, A0)|P3],
  645		sweep(Dim, I, Row0, M, P3, P4),
  646		I0 is I + 1,
  647		gauss(I0, M, Dim, P4, Q).
  648gauss(_,_,_, P, P).
  649
  650%
  651sweep(0, _I, _Row, _M, P, P):-!.
  652sweep(I, I, Row, M, P, Q):-!, K is I - 1,
  653	sweep(K, I, Row, M, P, Q).
  654sweep(K, I, Row, M, P, Q):- arg(K, M, KRow),
  655	( memberchk(I-V, KRow), V \== 0
  656	->	V0 is -V,
  657		sum_row_to_arg(K, M, Row, V0),
  658		P = [sweep(K, I, V0)|P1]
  659	;	P1 = P
  660	),
  661	K0 is K - 1,
  662	sweep(K0, I, Row, M, P1, Q).
  663
  664% ?- inverse_matrix(a([1-2, 2-4],[1-1,2-1]), Inv).
  665% ?- inverse_matrix(a([3- (-1)], [2-1], [1-1]), Inv).
  666inverse_matrix(Matrix, OutMatrix):-
  667	duplicate_term(Matrix, DupX),
  668	gauss(DupX, History),
  669	unit_matrix(DupX, OutMatrix),
  670	act_matrix(History, OutMatrix),
  671	mul_matrix(OutMatrix, Matrix, TobeUnit),
  672	writeln(Matrix * OutMatrix = TobeUnit).
  673
  674% act_matrix/2
  675act_matrix([], _).
  676act_matrix([A|As], M):-act_basic(A, M), !,
  677	act_matrix(As, M).
  678%
  679act_basic(swap_row(I,J), M):- swap_row(I, J, M).
  680act_basic(swap_column(I,J), M):- swap_column(I, J, M).
  681act_basic(mul_scalar(I, A), M):- arg(I, M, B),
  682	mul_scalar_row(A, B, B0),
  683	setarg(I, M, B0).
  684act_basic(sweep(K, I, V), M):- arg(I, M, IRow),
  685	mul_scalar_row(V, IRow, IRow0),
  686	arg(K, M, KRow),
  687	sum_row_row(IRow0, KRow, KRow0),
  688	setarg(K, M, KRow0)