```    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
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)```