```    1:- module(math, [boole/2,
2		 pairs/2, cartesian/3,
3		 powerlist/3, powerset/2,
4		 subst/3,
5		 nCr/3, combi/3, combination/4,
6		 fun_to_rel/2, rel_to_fun/2,  principal_filter/3,
7		 filter/3, inverse_image/3, image/3,
8		 bigcap/2,
9		 fact/2, factorial/2, factorial/3]).   10
11% ?- qcompile(misc(math)).
12:- use_module(pac(basic)).   13:- use_module(misc('pac-meta')).   14:- use_module(pac('expand-pac')).   15% :- expects_dialect(pac).
16term_expansion --> pac:expand_pac.
17:- use_module(pac(op)).   18
19
20% ?- qcompile(misc(math)).
21
22% ?- math:list_to_brace([[]], X).
23% ?- math:list_to_brace([[],[]], X).
24% ?- math:list_to_brace([[],[], []], X).
25%@ X = {{}, {{}, {{}}}}.
26
27% ?- listing(list_to_brace).
28
29
30
31:-bekind(list_to_brace, []).   32[] = `({}).
33[X] = `({})@X.
34[X|Y] = `({})@ (`(',')@X@Y).
35:- ekind.   36
37
38% ?- binomial(10, [a,b], E), maplist(val([a=(1/2), b=(1/2)]), E, F).
39		/**********************************
40		*     Eval Boolean expressions    *
41		**********************************/
42
43% ?-eval(math:cond(boolify, (true->true)), R).
44% ?-eval(math:cond(boolify, (1=1->1=2)), R).
45% ?-eval(math:cond(boolify, (1=2->1=2)), R).
46% ?-eval(math:cond(boolify, (1=1; 2=3)), X).
47% ?-eval(math:cond(boolify, (1=2 -> 2=3)), X).
48% ?-eval(math:cond(boolify, (2=2 -> 2=2)), X).
49% ?-eval(math:cond(boolify, (1=2 -> 2=2)), X).
50% ?-eval(math:cond(boole, true), R).
51% ?-eval(math:cond(boole,(-(true*false) +false)), R).
52% ?-eval(math:cond(boole, (-(true+false) +false)), R).
53% ?-eval(math:cond(boole, -(true + false ) +false), R).
54% ?-eval(math:cond(boole, -(true * false ) +false), R).
55
56cond(A1,(t->X),A2):-!,cond(A1,X,A2) .
57cond(_A1,(nil-> _A2),nil):-! .
58cond(A1,(X->Y),A2):-!,cond(A1,X,A3),cond(A1,(A3->Y),A2) .
59cond(_A1,(t;_A2),t):-! .
60cond(A1,(nil;X),A2):-!,cond(A1,X,A2) .
61cond(A1,(X;Y),A2):-!,cond(A1,X,A3),cond(A1,(A3;Y),A2) .
62cond(_A1,(nil,_A2),nil):-! .
63cond(A1,(true,X),A2):-!,cond(A1,X,A2) .
64cond(A1,(X,Y),A2):-!,cond(A1,X,A3),cond(A1,(A3,Y),A2) .
65cond(A1,-t,A2):-!,cond(A1,nil,A2) .
66cond(_A1,-nil,t):-! .
67cond(A1,-X,A2):-!,cond(A1,X,A3),cond(A1,-A3,A2) .
68cond(_A1,t,t):-! .
69cond(_A1,nil,nil):-! .
70cond(_A1,true,t):-! .
71cond(_A1,false,nil):-! .
72cond(S,X,A1):-call(S,X,A1) .
73
74% ?- listing(cond).
75% ?- boolify(a=a, X).
76
77boolify(X,true):-once(call(X)),! .
78boolify(_A1,false).
79
80
81% ?- scm(math).
82% ?- boole(true*false, X).
83
84boole(true*A,A1):-!,boole(A,A1) .
85boole(false* _A1,false):-! .
86boole(X*Y,A1):-!,(boole(X,A2),boole(Y,A3)),canonical_boole(*,A2,A3,A1) .
87boole(true+ _A1,true):-! .
88boole(false+A,A1):-!,boole(A,A1) .
89boole(X+A,A1):-!,(boole(X,A2),boole(A,A3)),canonical_boole(+,A2,A3,A1) .
90boole((X->Y),A1):-!,boole(-X+Y,A1) .
91boole(-true,false):-! .
92boole(-false,true):-! .
93boole(-X,A1):-!,boole(X,A2),canonical_boole(-,A2,A1) .
94boole(not(X),A1):-!,boole(X,A2),canonical_boole(-,A2,A1) .
95boole(true,true):-!.
96boole(false,false).
97
98
99canonical_boole(+, false, X, X).
100canonical_boole(+, X, false, X).
101canonical_boole(+, true, _, true).
102canonical_boole(+, _, true, true).
103canonical_boole(+, X, Y, X+Y).
104canonical_boole(*, true, X, X).
105canonical_boole(*, X, true, X).
106canonical_boole(*, false, _, false).
107canonical_boole(*, _, false, false).
108canonical_boole(*, X, Y, X*Y).
109%
110canonical_boole(-, true, false).
111canonical_boole(-, false, true).
112canonical_boole(-, X, -(X)).
113
114% ?- module(math).
115% ?- set_query_context(math).   % <== this is inserted by handle([sqc]).
116% ?- canonical_boole(A,B,C).
117% ?- eval(math:truth_table(false->true), R).
118% ?- eval(math:truth_table(false + true -> false), R).
119
120truth_table(-true,false):-! .
121truth_table(-false,true):-! .
122truth_table(not(true),false):-! .
123truth_table(not(false),true):-! .
124truth_table(+true,true):-! .
125truth_table(+false,false):-! .
126truth_table(false+false,false):-! .
127truth_table(_A1+ _A2,true):-! .
128truth_table(true*true,true):-! .
129truth_table(_A1* _A2,false):-! .
130truth_table((true->false),false):-! .
131truth_table((_A1-> _A2),true).
132
133% 		/***********************************
134% 		*     invering  power series       *
135% 		***********************************/
136
137% % ?- trace, math:power_series_inverse([1, 1 rdiv 2, 1 rdiv 6], Xs).
138% % Xs = [1, -1 rdiv 2, 1 rdiv 3]
139
140
141% % ?- trace, math:power_series_inverse([0,1], Xs).
142
143% power_series_inverse(Ds, Xs):- power_series(Ds, As, Ts),
144%     reverse([As|Ts], Ws),
145%     Ds =[A1|_],
146%     Rdiv_A is 1 rdiv A1,
147%     find_unknowns(Ws, [Rdiv_A], Ys),
148%     reverse(Ys, Xs).
149
150% %
151% power_series([A1|Ds], As, Ts) :-
152% 	power_series_loop(Ds, [A1], [], As, Ts).
153
154% %
155% power_series_loop([], As, Ts, As, Ts).
156% power_series_loop([D|Ds], As, Ts, Bs, Us):-
157% 	power_series_step(D, As, Ts, Cs, Ws),
158% 	power_series_loop(Ds, Cs, Ws, Bs, Us).
159
160% power_series_step(A, As, Ts, [A|As], Us):-
161%     reverse(As, As_rev),
162% 	innerproduct(As_rev, As, B),
163% 	power_series_step(As_rev, B, Ts, Us).
164% %
165% power_series_step(_, B, [], [[B]]).
166% power_series_step(As, B, [T|Ts], [[B|T]|Us]):-
167% 	innerproduct(As, T, C),
168% 	power_series_step(As, C, Ts, Us).
169
170% :- meta_predicate inverse_binary(3, ?, ?, ?).
171% %
172% inverse_binary(F, X, Y, Z):- call(F, Y, Z, X).
173% %
174% find_unknowns([_], Xs, Xs):-!.
175% find_unknowns([[A]|Ts], Xs, [X|Ys]):-
176% 	maplist(inverse_binary(pred([U,V,[U|V]), Ts, Cs, Us),
177% 	find_unknowns(Us, Xs, Ys),
178% 	innerproduct(Ys,Cs,W),
179% 	X is - (W rdiv A).
180
181% inner product of two "vectors" adjusting to the shorter
182% ?- innerproduct([1,2], [1,2,3], X).
183
184innerproduct(X,Y,Z):- innerproduct(X, Y, 0, Z).
185
186innerproduct([], _, X, X).
187innerproduct(_, [], X, X).
188innerproduct([A|As], [B|Bs], X, Y):-
189	X0 is X+ A * B,
190	innerproduct(As, Bs, X0, Y).
191
192% ?- poly([1,2]+[2,3], X).
193% ?- poly([1,2]+[2,3,4], X).
194% ?- poly(##(poly(1, [1,0,1]), poly(2, [1,1])), X).
195%@ X = poly(2,[1,1,0,0,1,3,3,1]) .
196% ?- poly([0,0,1]*[0,0,1], X).
197%@ X = [0,0,0,0,1] .
198% ?- poly([0,1]^8, X).
199%@ X = [0,0,0,0,0,0,0,0,1] .
200% ?- poly([0,1]^0 mod [1,0,1], X).
201% X = [1]
202% ?- poly([0,1]^1 mod [1,0,1]), X).
203%@ X = [0, 1] .
204% ?- poly([0,1]^2 mod [1,0,1], X).
205% ?- poly([0,1]^3 mod [1,0,1], X).
206% ?- poly([0,1]^4 mod [1,0,1], X).
207% ?- poly(##([0, 0, 0, 1], [1, 1]), X).
208% ?- poly(##([0, 0, 0, 0, 1], [1, 1]), X).
209% ?- poly(##([0, 0, 0, 0, 0, 1], [1, 1]), X).
210
211
212%%% "Little Endian" polynomial calculus
213% ?- poly( poly(2, [1,1]) + poly(3, [1,1]), X).
214% ?- poly( poly(2, [1,1])^2 + poly(3, [1,1])^3, X).
215%@ X = poly(4, [1, 2, 1, 0, 0, 1, 3, 3, 1]) .
216% ?- poly(poly(2, [1,1])*poly(3, [1,1]), X).
217%@ X = poly(5, [1, 2, 1]) .
218% ?- poly(poly(2, [1,1])^5, X).
219% ?- poly((- (1*2*3))^2, X).
220%@ X = [36]
221% ?- poly((- (1*2*3)), X).
222%@ X = [-6].
223% ?- poly((- (rdiv(1,2)*3)), X).
224%@ X = [-3 rdiv 2].
225% ?- poly((- (rdiv(2,4)*3)), X).
226%@ X = [-3 rdiv 2].
227
228% ?- listing(poly).
229
231poly(X-Y,A1):-!,poly(X+[-1]*Y,A1) .
232poly(-X,A1):-!,poly([-1]*X,A1) .
233poly(X*Y,A1):-!,(poly(X,A2),poly(Y,A3)),poly_mul(A2,A3,A1) .
234poly(X^Y,A1):-!,poly(X,X0),poly_exp(X0,Y,A1) .
235poly(##(X,Y),A1):-!,(poly(X,A2),poly(Y,A3)),poly_compose(A2,A3,A1) .
236poly(inverse(P),A1):-!,poly(P,A2),power_series_inverse(A2,A1) .
237poly(P mod Q,A1):-!,(poly(P,A2),poly(Q,A3)),poly_mod_little(A2,A3,A1) .
238poly(rev(P),A1):-!,poly(P,A2),reverse_endian(A2,A1) .
239poly(X,[X]):-rational(X),! .
240poly(X,X).
241
242% ?- poly_mul_bd_poly(poly(0, [1]), poly(0, []), 1, Z).
243% ?- poly_mul_bd_poly(poly(0, [1]), poly(0, [2]), 1, Z).
244% ?- poly_mul_bd_poly(poly(0, [1, 2, 1]), poly(0, [1, 1]), 3, Z).
245% ?- poly_mul_bd_poly(poly(0, [1, 2, 1]), poly(0, [1, 1]), 4, Z).
246% ?- poly_mul_bd_poly(poly(1, [1, 2, 1]), poly(2, [1, 1]), 7, Z).
247
248% Degree of the returned polynomial is less than N.
249
250poly_mul_bd_poly(U, V, N, poly(K, Z)):-
251	(U = poly(I, X); I=0, X=U), !,
252	(V = poly(J, Y); J=0, Y=V), !,
253	K is I + J,
254	N0 is max(N - K, 0),
255	poly_mul_bd(X, Y, N0, Z).
256
257% ?- poly_scalar_bd(1, [1], 2, X).
258% ?- poly_scalar_bd(1, [1, 2, 3], 5, X).
259% ?- poly_scalar_bd(1, [1, 2, 3, 5], 3, X).
260
261poly_scalar_bd(X, Y, Max, Z):-
262	poly_scalar_bd(X, Y, Max, Z, []).
263
264%
265poly_scalar_bd(_, [], _, U, U).
266poly_scalar_bd(_, _, 0, U, U).
267poly_scalar_bd(C, [A|B], J, [D|U], V):- D is A*C,
268	succ(J0, J),
269	poly_scalar_bd(C, B, J0, U, V).
270%
271reverse_endian(poly(X, Y), poly(X, Y0)):- !, reverse(Y, Y0).
272reverse_endian(X, X0):-  reverse(X, X0).
273
274%
275list_plus(A+B,A1):- list_plus(A,A2),
276	list_plus(B,A3),
277	append(A2,A3,A1).
278
279
280% ?- poly_compose(poly(1, [1,0,1]), poly(2,[1,1]), X).
281% ?- poly_compose(poly(1, [1,0,1]), poly(2,[1,1]), X).
282% ?- poly_compose(poly(2, [0,0,0,1]), [1,1], X).
283% ?- poly_compose([1, 1, 1], [1,1], X).
284% ?- poly_compose([1, 1, 1], [0,1], X).
285% ?- poly_compose([0, 1], [2, 3], X).
286
287poly_compose(poly(E0, X0), P, Q):-
288	poly_exp(P, E0, P0),
289	poly_compo0(X0, P, P1),
290	poly_mul(P0, P1, Q).
291
292poly_compose(X0, P, Q):- poly_compose_exact(X0, P, Q).
293
294% ?- scm(math).
295% ?- poly_compose_exact([], poly(1,[1]), X).
296% ?- poly_compose_exact([1], poly(1,[1]), X).
297% ?- poly_compose_exact([1, 1], poly(1,[1,2]), X).
298% ?- poly_compose_exact([1, 1, 1], poly(0,[1,1]), X).
299
300poly_compo0([], _, []).
301poly_compo0([A|As], P, Q):-
302	poly_compo0(As, P, P0),
303	poly_mul(P, P0, P1),
305
306poly_exp(poly(E0,X0), N, poly(E, X)):- bits(N, Bs),
307	poly_exp0(X0, Bs, [1], X),
308	E is E0*N.
309poly_exp(X0, N, X):- bits(N, Bs), poly_exp0(X0, Bs, [1], X).
310
311poly_exp0(_, [], X, X).
312poly_exp0(X, [B|R], Y, Z):-
313	poly_mul_tr(Y, Y, Y0),
314       (B==0 -> Y1=Y0; poly_mul_tr(Y0, X, Y1) ),
315	poly_exp0(X, R, Y1, Z).
316
317poly_scalar(A, poly(E,X), poly(E, Y)):- poly_scalar0(A, X, Y).
318poly_scalar(A, X, Y):- poly_scalar0(A, X, Y).
319
320poly_scalar0(0, _, []).
321poly_scalar0(A, X, Y):- maplist(scalar_mul(A), X, Y).
322
323scalar_mul(A, X, Y) :- Y is A*X.
324
325% ?- math:scalar_exp(1 rdiv 3, 4, R), writeln(R).
326% ?- math:scalar_exp(1 rdiv 3, 3, R), writeln(R).
327% ?- math:scalar_exp(1 rdiv 3, 2, R), writeln(R).
328% ?- math:scalar_exp(1 rdiv 3, 1, R), writeln(R).
329% ?- math:scalar_exp(1 rdiv 3, 0, R), writeln(R).
330
331scalar_exp(X, N, Y):- bits(N, Bs),
332	scalar_exp(Bs, X, 1, Y).
333
334scalar_exp([], _, Y, Y).
335
336scalar_exp([B|R], X, Y, Y0):-  Y1 is Y*Y,
337	( B==0 -> Y2 = Y1;  Y2 is X*Y1	),
338	scalar_exp(R, X, Y2, Y0).
339
340poly_mul(X, Y, Z):- number(X), poly_scalar(X, Y, Z).
341poly_mul(poly(E0,X0),poly(E1,X1), poly(E2,X2)):-
342	poly_mul_tr(X0,X1,X2), E2 is E0+E1.
343poly_mul(poly(E0,X0), X1, poly(E0,X2)):-poly_mul_tr(X0,X1,X2).
344poly_mul(X0, poly(E1,X1), poly(E1,X2)):-poly_mul_tr(X0,X1,X2).
345poly_mul(X0, X1, X2):-poly_mul_tr(X0,X1,X2).
346
347% ?- scm(math).
348% ?- poly_mul_tr([], [], X).
349% ?- poly_mul_tr([1], [2], X).
350% ?- poly_mul_tr([1], [1,1], X).
351% ?- poly_mul_tr([1,1], [1,1], X).
352% ?- poly_mul_tr([1,0,1], [1, 0,1], X).
353% ?- poly_mul_tr([1,1,1], [1, 1, 1], X).
354
355poly_mul_tr(A, B, C):- poly_mul_tr(A, B, [], [], C).
356
357poly_mul_tr([], _, W, S, X):- reverse(W, S, X).
358poly_mul_tr([A|As], B, W, S, X):-
359	poly_scalar(A, B, C),
361	( V == [] ->  U=0, S0=S; V=[U|S0]),
362	poly_mul_tr(As, B, [U|W], S0, X).
363
364% ?- scm(math).
365
366% ?- poly_mul_bd([], [], 3, X).
367% ?- poly_mul_bd([1], [2], 3, X).
368% ?- poly_mul_bd([1], [1,1], 3, X).
369% ?- poly_mul_bd([1,1], [1,1], 2, X).
370% ?- poly_mul_bd([1,1,1], [1, 1, 1], 3, X).
371% ?- poly_mul_bd([1,1,1], [1, 1, 1], 4, X).
372% ?- poly_mul_bd([1,1,1], [1, 1, 1], 5, X).
373% ?- poly_mul_bd([1,1,1], [1, 1, 1], 6, X).
374% ?- poly_mul_bd([1,1,1], [1, 1, 1], 7, X).
375
376poly_mul_bd(A, B, Max, C):- poly_mul_bd(A, B, Max, [], [], C).
377
378%
379poly_mul_bd(_, _, 0, W,  _, R):- !, reverse(W, R).
380poly_mul_bd([], _, N, W, S, X):- !,
382	reverse(W, S0, X).
383poly_mul_bd([A|As], B, N, W, S, X):-
384	poly_scalar_bd(A, B, N, C),
386	succ(N0, N),
387	( V = [] -> U = 0, S0 = S;  V=[U|S0]),
388	poly_mul_bd(As, B, N0, [U|W], S0, X).
389
390%
395
400
402	E is min(E0, E1),
403	D0 is E0-E,
404	D1 is E1-E,
405	shift(D0, X0, Y0),
406	shift(D1, X1, Y1),
408	(maplist(=(0), Y2) -> X=[]; X=Y2).
409
413
414poly_reverse0(poly(E, X), Y) :- repeat(E, cons(0), [], P), reverse(X, P, Y).
415poly_reverse0(X, Y) :- reverse(X, Y).
416
417cut_prefix(_, [], []).
418cut_prefix(X, [X|Xs], Y):- !, cut_prefix(X, Xs, Y).
419cut_prefix(_, Y, Y).
420
421shift(0, X, X):-!.
422shift(_, [], []):-!.
423shift(N, X, [0|Y]):- N1 is N-1, shift(N1, X, Y).
424
425reverse([], X, X).
426reverse([A|As], X, Y):- reverse(As, [A|X], Y).
427
428%  Big endian version
429% ?- poly_mod([1,0,0,-1], [1,1,1], X).
430%@ X = [0,0] .
431% ?- poly_mod([1,0], [1,0], X).
432%@ X = [0].
433% ?- canonical_mpoly([0, 1,0], A, B).
434%@ A = 1,
435%@ B = [0].
436
437poly_mod(P, Q, R):- canonical_mpoly(Q, A, Q0),
438	length(Q0, L0),
439	length(R0, L0),
440	maplist(=(0), R0),
441	poly_mod(P, Q0, R0, R1),
442	poly_scalar0(A, R1, R).
443
444poly_mod([], _, R, R).
445poly_mod([C|P], Q, [R0|Rs], R):-
446	poly_mod(R0, C, Q, Rs, R1),
447	poly_mod(P, Q, R1, R).
448
449poly_mod(A, C, [X|Xs], [Y|Ys], [Z|Zs]):- Z is Y+A*X, poly_mod(A, C, Xs, Ys, Zs).
450poly_mod(A, C, [B], [], [D]):- D is A*B + C.
451
452canonical_mpoly([0|P], A, Q):- !, canonical_mpoly(P, A, Q).
453canonical_mpoly([A|P], A, P0):- A0 is -1 rdiv A, poly_scalar0(A0, P, P0).
454
455% little-endinan poly_mod based on big-endian
456%
458poly_mod_little(P, Q, R):- big_endian(P, P0), big_endian(Q, Q0),
459	poly_mod(P0, Q0, R0),
460	little_endian(R0, R).
461
462big_endian(X, Y):- poly_reverse0(X, X0), cut_prefix(0, X0, Y).
463
464little_endian(X, Y):- cut_prefix(0, X, X0), reverse(X0, Y).
465
466const(F, _, L, V):- call(F, L, V).
467
468arith_curry(F, Y, Z):- build_arith_expr(F, Y, E), !, Z is E.
469
470build_arith_expr(+(X), Y, +(X,Y)).
471build_arith_expr(*(X), Y, *(X,Y)).
472build_arith_expr(^(X), Y, ^(X,Y)).
473build_arith_expr(mod(X), Y, mod(X,Y)).
474build_arith_expr(rdiv(X), Y, rdiv(X,Y)).
475build_arith_expr(E0, Y, E):- E0=..[F|A0], append(A0,[Y], As), E=..[F|As].
476
477%  arith/3
478arith(+, [X, Y|_], Z) :- Z is X+Y.
479arith(*, [X, Y|_], Z) :- Z is X*Y.
480arith(^, [X, Y|_], Z) :- Z is X^Y.
481arith(mod, [X, Y|_], Z) :- Z is X mod Y.
482arith(rdiv, [X, Y|_], Z) :- Z is X rdiv Y.
483arith(unary(OP), [X|_], V):- E=..[OP, X], V is E.
484arith(binary(OP), [X,Y|_], V):- E=..[OP, X, Y], V is E.
485arith(ternary(OP), [X,Y,Z|_], V):- E=..[OP, X, Y, Z], V is E.
486
487% ?- val(arith(1+2), X).
488% ?- val(arith(1 rdiv 2), X).
489% ?- val(arith(1 mod 2), X).
490
491%%%% Endo-map algebra
492%%%% endo map / permutation algebra %%%
493% *  : group operation
494% ^  : integer power
495% ?- endo(([1,3,2] ^ -103)*([3,1,2]^102) , X).
496
497perm(E, V) :- endo(E, V), (_V0= l(V); _V0=p(V); V= _V0).
498
499
500endo_default([],A1):-!,endo_default(l([]),A1) .
501endo_default([X|Y],A1):-!,endo_default(l([X|Y]),A1) .
502
503
504%%% rule of endo-map algebra
505
506
507endo(id(N),A1):-N0 is N,!,endo_id(v(N0),A1) .
508endo(shift(X),A1):-!,endo(X,A2),endo_cycle(A2,A1) .
509endo(shift(J,P),A1):-J0 is J,!,endo(P,A2),endo_cycle_v(J0,A2,A1) .
510endo(cycle(N),A1):-!,endo(shift(1,N),A1) .
511endo(X*Y,A1):-!,(endo(X,A2),endo(Y,A3)),endo_prod(A2,A3,A1) .
512% endo(X^Y,A1):-Y0 is Y,!,endo(X,A2),endo_aux_1(A2,A1) .
513endo(perms(X),A1):-!,list(X,A2),model:permutations(A2,A1) .
514endo(rand(N),A1):-(N0 is N,numlist(1,N0,L),shuffle(L,P)),!,endo(l(P),A1) .
515endo(cast(P),A1):-!,endo(P,A2), endo_aux_2(A2,A1) .
516endo(cast(F,P),A1):- endo(P,A2),cast(F,A2,A1) .
517
518endo_aux_1(A^[B],C):-endo_pow(A,B,C) .
519endo_aux_2(A,B):-'[|]'(v,A,B) .
520
521% ?- cast(v, l([1,2,3]), X).
522% X = v(1, 2, 3).
523cast(v, l(X), Y):- !, Y=..[v|X].
524cast(v, p(L), Y):- !, sort(L, L0), maplist(snd,L0,P), Y=..[v|P].
525cast(v, X, X):-!.
526cast(p, l(X), p(Y)):- !, length(X, N), numlist(1, N, Z), zip(Z,X, Y).
527cast(p, p(X), p(X)):-!.
528cast(p, X, p(Y)):- !, functor(X,_,N), X=..[_|X0], numlist(1, N, Z), zip(Z, X0, Y).
529cast(l, p(L), l(P)):- !, sort(L, L0), maplist(snd, L0, P).
530cast(l, l(X), l(X)):- !.
531cast(l, X, l(Y)):- X=..[_|Y].```
id
```  534endo_id(p, p([])).
535endo_id(l, l([])).
536endo_id(v(N), Id):- functor(Id, v, N), endo_id_v(N, Id).
537
538endo_id_v(0, _):-!.
539endo_id_v(J, Id):- arg(J, Id, J), J0 is J-1, endo_id_v(J0, Id).```
prod
?- `zip([a,b,c],[b,c,a], Z)`, `eval(endo, p(Z)^2, X)`. Z = [ (a, b), (b, c), (c, a) ], X = `p([ (a, c), (b, a), (c, b) ])`
```  546endo_prod(l(X), l(Y), l(Z)):- endo_prod_l(X, Y, Z).
547endo_prod(p(X), p(Y), p(Z)):- endo_prod_p(X, Y, Z).
548endo_prod(X, Y, Z):- endo_prod_v(X, Y, Z).
549
550endo_prod(F, X, Y, Z):- F==l-> endo_prod_l(X, Y, Z); F==p-> endo_prod_p(X, Y, Z); endo_prod_v(X, Y, Z).
551
552endo_prod_l(F, G, H):- endo_prod_l(1, F, G, H).
553
554endo_prod_l(_, _, [], []) :- !.
555endo_prod_l(N, F, [X|Xs], [Y|Ys]) :- (nth1(X, F, Y)->true; Y=X),  N0 is N-1, endo_prod_l(N0, F, Xs, Ys).
556
557endo_prod_p([], Y, Y).
558endo_prod_p(X, [], X).
559endo_prod_p(F, G, H):- maplist(endo_prod_p_(F), G, H).
560
561endo_prod_p_(F, (A,B), (A,C)):- memberchk((B,C), F) -> true; C=B.
562
563endo_prod_v(X, Y, Z):- functor(X, v, N), functor(Z, v, N), endo_prod_v(N, X, Y, Z).
564
565endo_prod_v(0, _, _,_):-!.
566endo_prod_v(J, X, Y, Z):- arg(J, X, K), arg(K, Y, L), arg(J, Z,  L),
567	J0 is J-1, endo_prod_v(J0, X, Y, Z).
568
569%%% inverse
570endo_inverse(l(X), l(Y)):- !, endo_inverse_l(X, Y).
571endo_inverse(p(X), p(Y)):- !, endo_inverse_p(X, Y).
572endo_inverse(X, Y):- endo_inverse_v(X, Y).
573
574endo_inverse(F, X, Y):- F==l -> endo_inverse_l(X, Y); F==v -> endo_inverse_v(X, Y); endo_inverse_p(X, Y).
575
576endo_inverse_l(L, M):- length(L, N), numlist(1, N, L0), zip(L, L0, L1),
577	lists:msort(L1, L2), unzip(L2, _, M).
578
579endo_inverse_v(P, P0):- functor(P, F,  N), functor(P0, F, N), endo_inverse_v(N, P, P0).
580
581endo_inverse_v(0, _, _):- !.
582endo_inverse_v(J, P, P0):- arg(J, P, A), arg(A, P0, J), J0 is J-1, endo_inverse_v(J0, P, P0).
583
584endo_inverse_p([], []).
585endo_inverse_p([(X,Y)|F], [(Y,X)|F0]):- endo_inverse_p(F, F0).```
cyclic
?- `endo_cycle(l([3,2,1]), X)`. X = `p([ (3, 2), (2, 1), (1, 3)])` ?- `endo_cycle(v(3,2,1), X)`. X = `v(1, 3, 2)`.
```  592endo_cycle(l(L), p(L0)):- endo_cycle_l_p(L, L0).
593endo_cycle(X, Y):- endo_cycle_v(1, X, Y).
594
595% ?- endo_cycle_v(1, v(1,2,3), X).
596% X = v(2, 3, 1).
597endo_cycle_v(J, P, Q):- functor(P, F, N),
598						functor(Q, F, N),
599						endo_cycle_v(N, N, J, P, Q).
600
601endo_cycle_v(0, _, _, _, _):- !.
602endo_cycle_v(K, N, J, P, Q):- arg(K, P, V),
603							  V0 is ((V+J-1) mod N)+1,
604							  arg(K, Q, V0), K0 is K-1,
605							  endo_cycle_v(K0, N, J, P, Q).
606% ?- scm(math).
607% ?- endo_cycle_l_p([3, 1, 2], X).
608endo_cycle_l_p([],  []).
609endo_cycle_l_p(Dom, Z):- Dom=[X|Y],
610						 append(Y,[X], Dom0),
611						 maplist(pred([X,Y,(X,Y)]), Dom, Dom0,Z).```
power
?- `endo_pow(l([2,1]), 2, X)`. X = `l([1, 2])` ?- `endo_pow(p([(2,1), (1,2)]), 2, X)`. X = `p([ (2, 2), (1, 1)])` ?- `endo_pow(v(2,1), 2, X)`. X = `v(1, 2)`
```  621endo_pow(l(L), N, l(Y)):- !, endo_pow(l, L, N, [], Y).
622endo_pow(p(L), N, p(Y)):- !, endo_pow(p, L, N, [], Y).
623endo_pow(X, N, Y):- functor(X, _, A), endo_id(v(A), Id), endo_pow(v, X, N, Id, Y).
624
625endo_pow(F, X, N, Id, Y):- N<0, !, N0 is -N, bits(N0, Bits), endo_pow_pos_b(F, X, Bits, Id, X0),
626	endo_inverse(F, X0, Y).
627endo_pow(F, X, N, Id, Y):- bits(N, Bits), endo_pow_pos_b(F, X, Bits, Id, Y).
628
629endo_pow_pos_b(_, _, [], X, X).
630endo_pow_pos_b(F, X, [D|R], Y, Z):- endo_prod(F, Y, Y, Y0),
631      ( D==0 -> Y1=Y0; 	endo_prod(F, Y0, X, Y1) ),
632	endo_pow_pos_b(F, X, R, Y1, Z).
633
634%%%%% conjugate
635% ?- conjugate(v(2,3,1), v(3,1,2)).
636% ?- conjugate(l([2,3,1]), l([1,2,3])).
637
638singleton(X, [X]).
639
640conjugate(X, Y):- endo_type(X, T), endo_type(Y, T).
641
642% ?- conjugate_class([1,2,3], X), maplist(writeln, X).
643%@ X = [[[1,2,3]],[[1,3,2],[2,1,3],[3,2,1]],[[2,3,1],[3,1,2]]] .
644% [[1, 2, 3]]
645% [[1, 3, 2], [2, 1, 3], [3, 2, 1]]
646% [[2, 3, 1], [3, 1, 2]]
647
648conjugate_class(L, C):- eval(endo, perms(L), Ps),
649	maplist(pred([X, l(X)]), Ps, Qs),
650	quotient_set(Qs, endo_type, C).
651
652quotient_set(S, F, Q):-  maplist(F, S, T), zip(T, S, Z),
653	sort(Z, Z0), merge_pairs(Z0, Q0),
654	maplist(maplist(pred([X,Y]:-arg(1, X, Y))), Q0, Q).```
type
?- `endo_type(l, [2,1,3], X)`. X = [1, 2]. ?- `endo_type(v(2,1,3), X)`. X = [1, 2]. ?- `endo_type(p([(1,2), (2,1), (3,3)]), X)`. X = [1, 2]
```  664endo_type(l(X), T):- !, endo_type_l(X, T).
665endo_type(p(X), T):- !, endo_type_p(X, T).
666endo_type(X, T):-  X=..[v|L], endo_type_l(L, T).
667
668endo_type_l(L, T):- maplist(singleton, L, L0),
669	union_find_l(1, L, L0, Ks),
670	maplist(length, Ks, Ls),
671	msort(Ls, T).
672
673endo_type_p(X, T):- field(X, F),
674	maplist(singleton, F, L0),
675	union_find_p(X, L0, Ks),
676	maplist(length, Ks, Ls),
677	msort(Ls, T).
678
679% Union-find
680% ?- union_find_l(1, [2,1], [[1],[2]], X).
681% X = [[1, 2]].
682union_find_l(_, [], Bs, Bs):- !.
683union_find_l(X, [Y|E], Bs0, Bs):- find_cluster(X, Bs0, CX, Bs1),
684   (	memberchk(Y, CX)  ->  Ps = Bs0
685   ;    find_cluster(Y, Bs1, CY, Bs2),
686        append(CX, CY, U),
687	Ps = [U|Bs2]
688   ),
689  X0 is X+1, union_find_l(X0, E, Ps, Bs).
690
691
692% ?- union_find_p([(1,2),(2,1)], [[1],[2]], X).
693% X = [[1, 2]].
694union_find_p([], Bs, Bs):- !.
695union_find_p([(X,Y)|E], Bs0, Bs):- find_cluster(X, Bs0, CX, Bs1),
696   (	memberchk(Y, CX)  ->  Ps = Bs0
697   ;    find_cluster(Y, Bs1, CY, Bs2),
698        append(CX, CY, U),
699	Ps = [U|Bs2]
700   ),
701  union_find_p(E, Ps, Bs).
702
703%%% some tiny
704meet(A, B):- member(X, A), member(X, B).
705
706find_cluster(A, B, C, D):- select(C, B, D), memberchk(A, C), !.
707
708field(L, L):- listp(L), !.
709field(l(L), F):- !, length(L, N), numlist(1, N, F).
710field(p(L),  F):- !, length(L, N), numlist(1, N, F).
711field(X, F):- X=..[v|F].
712
713add_pair((X,Y), L, L0):- sort([X,Y], P), union(P, L, L0).
714
715% merge
716% ?- merge_pairs([(a,b), (a,c)], X).
717% X = [[b, c]].
718
719merge_pairs(X, Y):- merge_pairs_(X, X0), maplist(snd, X0, Y).
720
721merge_pairs_([], []).
722merge_pairs_([(A,B)|R], X):- merge_pairs_(R, Y), merge_pair(A, B, Y, X).
723
724merge_pair(A, B, [(A,X)|R], [(A,[B|X])|R]):-!.
725merge_pair(A, B, R, [(A,[B])|R]).
726
727% apply endo map
728% ?- endo_apply(l([2,1,3]), 2, X).
729% X = 1.
730
731endo_apply(l(F), A, B):- !, nth1(A, F, B).
732endo_apply(p(F), A, B):- !, memberchk((A,B), F).
733endo_apply(F, A, B):- arg(A, F, B).
734
735% group action orbits
736% ?-  group_orbits(l([2,1,3]), X).
737% X = [[1, 2], [3]].
738% ?- group_orbits(l([2,1,3]), [1,2,3], X).
739% X = [[1, 2], [3]].
740
741group_orbits(F, Y):- field(F, D), group_orbits(F, D, Y).
742
743group_orbits(F, X, Y):- group_orbits(F, X, [[]], Y0),
744	maplist(reverse, Y0, Y1), remove([], Y1, Y2),
745	sort(Y2, Y).
746
747group_orbits(_, [], O, O):-!.
748group_orbits(F, [X|Xs],[[]|R], O):- !, group_orbits(F, Xs, [[X]|R], O).
749group_orbits(F, Xs, [H|R], O):- H=[A|_], endo_apply(F, A, A0),
750	(  memberchk(A0, H) -> group_orbits(F, Xs, [[], H|R], O)
751	;  remove(A0, Xs, Ys), group_orbits(F, Ys, [[A0|H]|R], O) ).
752
753%%%
754eval_set_class--> eval([set, subset, ps]).
755set_class  -->  eval(set).
756
757setcalc(X, Y, Z, U, V):- set(X, Y, Z, U, V).  % for compatibility
758
759shorter(<, X, Y):- length(X, L), length(Y, M), L<M, !.
760shorter(>, _, _).
761
762seq(X, Y) :- maplist(eval, X, Y).
763
764seq(S, X, Y) :- maplist(eval(S), X, Y).```
```  768boole_complement(true, false).
769boole_complement(false, true).
770
771% ?- boole(fff, a+b, R).
772% R = true
773% ?- boole(fff, a*b, R).
774% R = false
775
776fff(a, true).
777fff(_, false).
778
779% :- edit(set/2).
780
781
782% ?- time(repeat(1000000, eval([boole(prolog)], 1==1, X))).
783% is faster than this:
784% ?- time(repeat(1000000, eval([boole, prolog], 1==1, X))).
785
786% The index directive :- idnex(boole(0,1,0,0,0,0)) shows
787% no significant effect.
788
789% ?- eval(subset(set), [1,2]+[2,1] < [1,2]+[3], X).
790% ?- eval(subset(set), [1,2]+[2,1]+[3,2] < [1,2]+[3], X).
791% ?- eval(subset(set), [3]+[1]+[2] =< [1,2]+[3], X).
792% ?- eval([boole(subset(set))], ([1]=<[1,2]) * ([2,1]>=[1,2,3]), X).
793% ?- eval([boole(subset(set))], pow([1]+[2]) < pow([1,2]+[3]), X).
794% ?- eval([boole(subset(set))], [1]+[2] < [1,2], X).
795% ?- eval([boole(prolog)], 1==1, X).
796% ?- eval([boole(frprolog)], (between(1, 5, 3)* member(a,[b,a])), X).
797% ?- eval([boole(prolog)], 1==1 -> 2==1, X).
798% ?- listing(boole).
799
800% ?- cond(prolog, (true->fail;fail); fail, R).
801% R = nil
802
803% ?- binomial(10, [a,b], E).
804nil_complement(true, nil).
805nil_complement(nil, true).
806
807%
808proper_subset([], [_|_]):-!.
809proper_subset([X|Y], [X|Z]):- proper_subset(Y,Z).
810
811
812% powerset(X, Y):- once(foldl(powerset, X,  Y, [[]])).
813
814% powerset(_,[], []).
815% powerset(X,[Y,[X|Y]|U],[Y|Z]):- powerset(X, U, Z).
816
817%  powerset(X, Y)  is true if Y is the powerset of X.
818powerset(X, Y):- powerset(X, [[]], Y).
819
820powerset([], X, X).
821powerset([A|R], X, Y):-
822	powerset(X, A, X, X0),
823	powerset(R, X0, Y).
824
825%
826powerset([], _, X, X).
827powerset([X|R], A, S, Y):- powerset(R, A, [[A|X]|S], Y).
828
829%?- ya_powerset([], X).
830%@ X = [[]].
831%?- ya_powerset([1], X).
832%@ X = [[1], []].
833%?- ya_powerset([1,2], X).
834%@ X = [[2], [2, 1], [1], []].
835%?- ya_powerset([1,2,3,4], X), sort(X, X0).
836%@ X = [[4], [4, 1], [4, 2, 1], [4, 2], [4, 3, 2], [4, 3, 2|...], [4, 3|...], [4
837%?- powerset([1,2,3,4], X).
838%@ X = [[], [1], [2], [1, 2], [3], [1, 3], [2, 3], [1|...], [...]|...].
839
840% What is the last argument for ?
841% xreverse(Xs, Ys) :-
842% 	xreverse(Xs, [], Ys, Ys).
843
844% xreverse([], Ys, Ys, []).
845% xreverse([X|Xs], Rs, Ys, [_|Bound]) :-
846% 	xreverse(Xs, [X|Rs], Ys, Bound).
847
848% ?- trace, xreverse([1,2,3|R], X).
849
850% ?-  sublist(X, [a,b,c]).
851sublist([], []).
852sublist(X, [_|As]):- sublist(X, As).
853sublist([A|X], [A|As]):- sublist(X, As).
854
855% ?- permutations([a,b,c],X).
856% X = [[a, b, c], [b, a, c], [b, c, a], [a, c, b], [c, a, b], [c, b, a]]
857permutations([], [[]]).
858permutations([X|Xs], Ps):- permutations(Xs, Qs),
859	permutations(X, Qs, Ps).
860
861permutations(_, [], []).
862permutations(X, [P|Ps], L):- insert_list(X, P, P0),
863	permutations(X, Ps, Qs),
864	append(P0, Qs, L).
865
866% ?- disjoint_pairs([1,2], X).
867disjoint_pairs(X, S) :- bagof((L, R), disjoint_pair(X, L, R), S).
868
869disjoint_pair(X, L, R):-  sublist(L, X),
870						  subtract(X, L, L0),
871						  sublist(R, L0).
872
873% ?- insert_list(a, [1,2], X).
874% X = [[a, 1, 2], [1, a, 2], [1, 2, a]]
875insert_list(X, [], [[X]]).
876insert_list(X, [A|As], [[X,A|As]|L]):-insert_list(X, As, L0),
877	maplist(cons(A), L0, L).
878
879% ?- disjoint_family([[a], [], [b, c]]).
880% ?- disjoint_family([[a], [b], [b, c]]).
881% [2014/09/08]
882disjoint_family([]).
883disjoint_family([A|R]):- disjoint(R, A),
884	disjoint_family(R).
885
886disjoint([], _).
887disjoint([A|_], B):- member(X, A), member(X, B), !, fail.
888disjoint([_|R], B):- disjoint(R, B).
889
890
891% %%% Interesting !!
892% ?- powerset_rev([a,b], X).
893% X = [[], [b], [a], [b, a]]
894
895% powerset_rev(X, Y):- forlist(powerset_rev, X,  [[]], Y).
896
897% powerset_rev(_,[], []).
898% powerset_rev(X,[Y|Z],[Y,[X|Y]|U]):-powerset_rev(X, Z, U).
899
900%
902	foldr(const(inverse(cons([]))), X, P0, [[[]]]),
903	foldl(carp_up, X, P0, P).
904
905graded_powerset(N, X, P):- repeat(N, cons([]), [[[]]], P0),
906	foldl(carp_up, X, P0, P).
907
908%
909carp_up(X, P, Q):- maplist(maplist(cons(X)), P, [_|P0]),
910       zip_append(P, P0, Q).
911
912%
913zip_append([], A, A).
914zip_append(A, [], A).
915zip_append([A|X], [B|Y], [C|Z]):- append(A, B, C), zip_append(X, Y, Z).
916
917% ?- trace, math:cartesian([a,b],[c,d], R).
918
919cartesian(Xs,Ys,Zs) :- prod(basic:ordered_pair, Xs, Ys, Zs).
920
921cartesian(P, Xs, Ys, Zs) :- prod(P, Xs, Ys, Zs).
922
923% ?- product(ordered_pair, [1,2], [a,b,c], X).
924
925product(F, Xs, Ys, Zs):- prod(F, Xs, Ys, Zs, []).
926
927:- meta_predicate prod(:, ?, ?, ?).  928prod(F, X, Y, Z):- prod(F, X, Y, Z, []).
929
930prod(_, [], _, Z, Z).
931prod(F, [A|As], Y, Z, U):- prod_(F, A, Y, Z, V),
932	prod(F, As, Y, V, U).
933
934prod_(_,_,[],Z,Z).
935prod_(F, A, [B|Bs], [C|Cs],U):- call(F, A, B, C),
936	prod_(F, A, Bs, Cs, U).
937
938%  ?- pairs([a, b, c], R).
939%@ R = [(a, b),  (a, c),  (b, c)] .
940
941pairs(A, P):- pairs(A, P, []).
942
943pairs([], X, X).
944pairs([A|X], L, N):- foldl(pred(A, [B, [(A,B)|U], U]), X, L, M),
945	pairs(X, M, N).
946
947% ?- trace, pairs([a,b,c], X).
948% X = [ (a, b), (a, c), (b, c)]
949
950% ?- mapset([a,b,c],[1,2,3,4], X), length(X, L).
951% ?- mapset([],[1], X), length(X, L).
952% ?- mapset([1],[], X), length(X, L).
953
954mapset(Xs, Ys, Fs):- mapset(Xs, Ys, [[]], Fs).
955
956mapset([], _, X, X):-!.
957mapset([X|Xs], Ys, Gs, Fs) :- mapset_extend(Ys, X, Gs, [], Us),
958	mapset(Xs, Ys, Us, Fs).
959
960mapset_extend([], _, _, Fs, Fs).
961mapset_extend([Y|Ys], X, Zs, Gs, Fs):-
962	maplist(cons((X,Y)), Zs, Us),
963	append(Us,Gs,Hs),
964	mapset_extend(Ys, X, Zs, Hs, Fs).
965
966% cons(X, Y, [X|Y]).
```
?- `direct_sum([a,b,c], [a,b,c], X)`.
```  970direct_sum(X, Y, Z):- direct_sum(0, 1, X, Y, Z).
971
972direct_sum(A, B, X, Y, Z):-maplist(ordered_pair(A), X, X1),
973	maplist(ordered_pair(B), Y, Y1),
974	append(X1, Y1, Z).
975
976% ?- nlist([a,b], 3, X).
977%@ X = [[], [a], [a, a], [a, a, a], [a, a, b], [a, b], [a, b, a], [a, b, b], [b], [b, a], [b, a, a], [b, a, b], [b, b], [b, b, a], [b, b, b]] .
978
979nlist(A, N, L):- powerlist_sum(A, N, L0), append(L0, L1), sort(L1, L).
980
981% ?- powerlist_sum([a,b],2,X).
982% X = [[[a, a], [b, a], [a, b], [b, b]], [[a], [b]], [[]]]
983
984powerlist_sum(A, N, L):- repeat(N, powerlist_cons(A), [[[]]], L).
985
986powerlist_cons(A, [X0|Xs], [X, X0|Xs]):- product(cons, A, X0, X).
987
988powerlist(A, N, L):- meta:repeat(N, product(cons, A), [[]], L).
989
990%
991keta(N, R, C) :- keta(N, R, 0, C).
992
993keta(N, R, C, C1):- N mod R =:= 0, !, C2 is C + 1, N1 is N// R,
994	keta(N1, R, C2, C1).
995keta(_, _, C, C).
996
997append(&(X,Y), Z, [X,Y], [U,V], append(U,V,Z)).
998append(X, X, [], [], true) :- listp(X).
999append(X,[X],[], [], true).
1000
1001bits(N, X) :- bits(N, X, []).
1002
1003bits(0, X, X):-!.
1004bits(N, X, Y):- B is N mod 2, N0 is N>>1,
1005	bits(N0, X, [B|Y]).
1006
1007boole_rule_OR_and_over_or(X, Y):-
1008	boole_rule(X, Y)
1009	;
1010	and_over_or(X, Y).
1011
1012% for compatibility
1013
1014conjunctive_normal_form--> boole_cnf, simple_cnf.
1015
1016disjunctive_normal_form--> boole_dnf, simple_dnf.
1017
1018% ?- boole_cnf(or(and(a,b), and(c,d)), X).
1019%@ X = [[a, c], [a, d], [b, c], [b, d]] .
1020
1021% ?- boole_dnf(or(and(a,b), and(c,d)), X).
1022
1023% ?- boole_cnf(not(or(and(a,b), and(c,d))), X).
1024%@ X = [[not(a),not(b)],[not(c),not(d)]] .
1025% ?- boole_cnf(not(or(not(and(a,b)), and(c,d))), X).
1026%@ X = [[a],[b],[not(c),not(d)]] .
1027
1028% ?-  trace, boole_cnf(not(and(a,b)), X).
1029% ?- listing(boole_cnf).
1030
1031
1032
1033boole_cnf(and(A,B),A1):-!,(boole_cnf(A,A2),boole_cnf(B,A3)),union(A2,A3,A1) .
1034boole_cnf(or(A,B),A1):-!,(boole_cnf(A,A2),boole_cnf(B,A3)),union_product(A2,A3,A1) .
1035boole_cnf(not(and(A,B)),A1):-!,boole_cnf(or(not(A),not(B)),A1) .
1036boole_cnf(not(or(A,B)),A1):-!,boole_cnf(and(not(A),not(B)),A1) .
1037boole_cnf(not(not(A)),A1):-!,boole_cnf(A,A1) .
1038boole_cnf(not(false),A1):-!,boole_cnf(true,A1) .
1039boole_cnf(not(true),A1):-!,boole_cnf(false,A1) .
1040boole_cnf(true,[]):-! .
1041boole_cnf(false,[[]]):-! .
1042boole_cnf(A,[[A]]):-! .
1043
1044boole_dnf(X, Y):- boole_cnf(not(X), Z),
1045	maplist(maplist(pred([not(A),A]&[B,not(B)])), Z, Y).
1046
1047% ?- boole_dnf(and(a,b), Y).
1048% ?- boole_cnf(and(a,b), Y).
1049% ?- boole_dnf(true, Y).
1050% ?- boole_dnf(false, Y).
1051
1052% :- bekind(boole_dnf).
1053% or(A,B) = :union@A@B.
1054% and(A,B) =  :union_product@A@B.
1055% not(and(A,B)) = or(not(A), not(B)).
1056% not(or(A,B)) = and(not(A), not(B)).
1057% not(not(A)) = A.
1058% not(false) = true.
1059% not(true) = false.
1060% A  = `[[A]].
1061% :- ekind.
1062
1063union_product([], _, []).
1064union_product([A|R], X, P):- maplist(union(A), X, Q),
1065	union_product(R, X, Q0),
1066	union(Q, Q0, P).
1067
1068% ?- simple_cnf([[a],[b], [a], [not(c),not(d)]], X).
1069%@ X = [[a],[b],[not(c),not(d)]].
1070
1071% ?- simple_cnf([[a],[b], [a], [not(c),not(d), c]], X).
1072%@ X = [[a],[b]].
1073
1074simple_cnf(X, Y):- simple_cnf_one(X, X0), !, simple_cnf(X0, Y).
1075simple_cnf(X, Y):- maplist(sort, X, X0),
1076	sort(X0, Y).
1077
1078simple_cnf_one(X, Y):- select(C, X, X0), cnf_rule(C, X0, Y).
1079
1080cnf_rule([], X, [[]]):- X\==[].
1081cnf_rule(C, X, X):- memberchk(true, C).
1082cnf_rule(C, X, [C0|X]):- select(false, C, C0).
1083cnf_rule(C, X, X):- select(A, C, C0),
1084	complementary(A, B),
1085	memberchk(B, C0).
1086
1087complementary(not(A), A):-!.
1088complementary(A, not(A)).
1089
1090% ?- simple_dnf([[a],[b], [a], [not(c),not(d), c]], X).
1091%@ X = [[a],[b]].
1092
1093simple_dnf(X, Y):- simple_dnf_one(X, X0), !, simple_dnf(X0, Y).
1094simple_dnf(X, Y):- maplist(sort, X, X0),
1095	sort(X0, Y).
1096
1097simple_dnf_one(X, Y):- select(C, X, X0), dnf_rule(C, X0, Y).
1098
1099dnf_rule([], X, X).
1100dnf_rule(C, X, X):- memberchk(false, C).
1101dnf_rule(C, X, [C0|X]):- select(true, C, C0).
1102dnf_rule(C, X, X):- select(A, C, C0),
1103	complementary(A, B),
1104	memberchk(B, C0).
1105
1106
1107flatten_boole(X,A1):-(functor(X,F,2),member(F,[or,and,&]),flatten(F,X,L0),maplist(flatten_boole,L0,L1),Y=..[F,L1]),!,flatten_boole(Y,A1) .
1108flatten_boole(-X,not(A1)):-!,flatten_boole(X,A1) .
1109flatten_boole(not(X),not(A1)):-!,flatten_boole(X,A1) .
1110
1111
1112boole_flat_paramod(and(L), and([X0|L0]), X, X0):- select(X, L, L0).
1113boole_flat_paramod(or(L), or([X0|L0]), X, X0):- select(X, L, L0).
1114
1115% boole_rule/2
1116boole_rule(not(not(L)), L).
1117boole_rule(not(true), false).
1118boole_rule(not(false), true).
1119boole_rule(not(or(L)), and(L0)):- maplist(unary(not), L, L0).
1120boole_rule(not(and(L)), or(L0)):- maplist(unary(not), L, L0).
1121boole_rule(or([]), false).
1122boole_rule(or([X]), X).
1123boole_rule(or(L), true):- member(true, L).
1124boole_rule(or(L), or(L0)):- select(false, L, L0).
1125boole_rule(or(L), or(L0)):- sort(L, L0), L\==L0.
1126boole_rule(or(L), true):-  memq(not(X), L), memq(X, L).
1127boole_rule(or(L), or(L0)):- select(or(X), L, L1), append(X, L1, L0).
1128boole_rule(and([]), true).
1129boole_rule(and([X]), X).
1130boole_rule(and(L), false):- member(false, L).
1131boole_rule(and(L), and(L0)):- select(true, L, L0).
1132boole_rule(and(L), and(L0)):- sort(L, L0), L\==L0.
1133boole_rule(and(L), false):- memq(not(X), L), memq(X,L).
1134boole_rule(and(L), and(L0)):- select(and(X), L, L1), append(X, L1, L0).
1135
1136and_over_or(and(L), or(L1)):- distribute(and, or, L, L1).
1137
1138or_over_and(or(L), and(L1)):- distribute(or, and, L, L1).
1139
1140% ?- term_rewrite:term_rewrite(desugar, ^, (a-b)^(-x), X).
1141% ?- reduce:term_rewrite(desugar, +, (a-b)^(-x), X).
1142% ?- reduce:term_rewrite(desugar, (^,[1]), (a-b)^(-x), X).
1143% ?- reduce:term_rewrite(desugar, (^,[2]), (a-b)^(-x), X).
1144% ?- reduce:subtree((+);(*);(/);(^,[2]), (a*b)^(2/3), X, Y, Z).
1145
1146% rule for reduce
1147% ?- reduce:term_rewrite(ring:desugar, [], -(a-b), X).
1148%@ X = -1* (a+ -1*b).
1149
1150desugar(X-Y, X+(-1)*Y).
1151desugar(+X, X).
1152desugar(-X, (-1)*X).
1153desugar(X/Y, rdiv(X, Y)).
1154
1155% ?- scm(math).
1156% ?- ya_desugar(1, X).
1157% ?- nb_setval(context_module, math).
1158
1159
1160ya_desugar(+X,A1):-!,ya_desugar(X,A1) .
1161ya_desugar(X-Y,A1+(-)*A2):-!,ya_desugar(X,A1),ya_desugar(Y,A2) .
1162ya_desugar(-X,-1*A1):-!,ya_desugar(X,A1) .
1163ya_desugar(X/Y,A1 rdiv A2):-!,ya_desugar(X,A1),ya_desugar(Y,A2) .
1164ya_desugar(X,X):-! .
1165prolog(A,true):-(call(A),!),! .
1166prolog(_A1,false):-! .
1167
1168% ?- math:subst([1-true, 0-false], f(1,0,1), R).
1169%@ R = f(true, false, true).
1170subst_atomic(M, X, Y)	:- atomic(X), !,
1171					   ( key_val_check(M, X, Y) -> true
1172					   ;	Y = X
1173					   ).
1174subst_atomic(M, X, Y)	:-
1175        functor(X, F, N),
1176		functor(Y, F, N),
1177        subst_atomic_args(N, M, X, Y).
1178%
1179subst_atomic_args(0, _, _, _) :- !.
1180subst_atomic_args(J, M, X, Y) :-
1181        arg(J, X, A),
1182        arg(J, Y, B),
1183        subst_atomic(M, A, B),
1184        K is J-1,
1185        subst_atomic_args(K, M, X, Y).
1186%
1187key_val_check([A|As], X, Y) :-
1188	(	(A = (X-Y); A = (X,Y); A = (X:Y) ) -> true
1189	;	key_val_check(As, X, Y)
1190	).```
``` 1193subst(M, X, Y):- subst_atomic(M, X, Y).
1194
1195%  ?- time(fibo(100000, N)).
1196%%%%
1197fibo(0, 1, _).
1198fibo(1, 1, 1).
1199fibo(N, F, F0):- N>1, N0 is N-1, fibo(N0, F0, G), F is F0+G.
1200
1201fibo(N, F) :- fibo(N, F, _).```
?- `time(fibo_tr(100000, N))`.
``` 1206fibonacci(N, F):- fibo_tr(N, F).
1207
1208fibo_tr(0, _, Y, Y).
1209fibo_tr(N, X, Y, Z):- N>0, N0 is N-1, Z0 is X + Y,
1210	fibo_tr(N0, Y, Z0, Z).
1211
1212fibo_tr(0, 1):-!.
1213fibo_tr(1, 1):-!.
1214fibo_tr(N, F):- N > 1, N0 is N-1, fibo_tr(N0, 1, 1, F).
1215
1216% ?-arrange([1,2,3],L).
1217% ?-arrange([1,2,3],[X,Y]).
1218arrange(_, []).
1219arrange(L, [X|R]):- select(X, L, L0),
1220	arrange(L0, R).
1221
1222% %%% shuffling
1223% ?- select_random(A, [taro, jiro, hanako], X).
1224% A = hanako,
1225% X = [taro, jiro]
1226% ?- shuffle([taro,jiro,hanako], X).
1227% X = [hanako, taro, jiro]
1228
1229shuffle(X, Y):- shuffle(X, [], Y).
1230shuffle([], A,  A):-!.
1231shuffle(L, A,  B):- length(L, N),
1232	N0 is random(N),
1233	select_nth(N0, X, L, L0),
1234	shuffle(L0, [X|A], B).
1235
1236select_random(X, L, L0):- length(L, N),
1237	I is random(N),
1238	select_nth(I, X, L, L0).
1239
1240select_nth(0, X, [X|L], L).
1241select_nth(N, X, [Y|L], [Y|L0]):- N>0, N0 is N-1,
1242	select_nth(N0, X, L, L0).
1243
1244random_sublist(L, L0):- length(L, N), N0 is N+1,
1245	R is random(N0),
1246	length(L0, R),
1247	random_sublist_(L0, L).
1248
1249random_sublist_([], _).
1250random_sublist_([X|Xs],  L):- select_random(X, L, L0),
1251	random_sublist_(Xs, L0).
1252
1253%%%
1254modulo(M, X, V):-  V is X mod M.
1255
1256%
1257g_base(lex, X, Y):- !, ensure_loaded('gb-lex'), gblex:gb(X, Y).
1258g_base(total, X, Y):-  ensure_loaded('gb-total'), gbtotal:gbase(X, Y).
1259
1260% split(F, X, Y, [X0, X1], [Y0,  Y1], append(Y0, Y1, Y)):-
1261% 	functor(X, F, 2),
1262% 	arg(1, X, X0),
1263% 	arg(2, X, X1).
1264% split(_, X, [X], [], [], true).
1265
1266%?- rel_to_fun([(a,b),(a,c),(b,d)], X).
1267%@ X = [b-[d], a-[c, b]].
1268
1269rel_to_fun(X, Y):-rel_to_fun(X, [], Y).
1270
1271%
1272rel_to_fun([], X, X).
1273rel_to_fun([P|R], X, Y):- (P=(A-B); P=(A,B)), !,
1274	(	select(A-M, X, X0)
1275	->	rel_to_fun(R, [A-[B|M]|X0], Y)
1276	;	rel_to_fun(R, [A-[B]|X], Y)
1277	).
1278
1279% ?- rel_to_fun([(a,b),(a,c),(b,d)], X), fun_to_rel(X, Y).
1280%@ X = [b-[d], a-[c, b]],
1281%@ Y = [b-d, a-c, a-b].
1282
1283fun_to_rel(X, Y):- fun_to_rel(X, Y, []).
1284%
1285fun_to_rel([], Y, Y).
1286fun_to_rel([A|X], Y, Z):-  (A = (U,L); A = (U-L)), !,
1287	fun_to_rel(L, U, Y, Y0),
1288	fun_to_rel(X, Y0, Z).
1289%
1290fun_to_rel([], _, Z, Z).
1291fun_to_rel([V|Vs], U, [U-V|Y], Z):-
1292	fun_to_rel(Vs, U, Y, Z).
1293
1294%?- math:filter([a,b,c],[c], X).
1295%@ X = [[c],[c,a],[c,b],[c,a,b]].
1296%?- filter([c], [a,b,c],X).
1297
1298filter(X, D, F):- subtract(X, D, Y),
1299	powerset(Y, P),
1300	maplist(pred(D, [X,Y] :- append(D,X,Y)),
1301		P, F).
1302
1303
1304
1305% ?- principal_filter([a,b,c], a, X).
1306principal_filter(D, A, PF):- select(A, D, D0), !,
1307	powerset(D0,  PD),
1308	maplist(pred(A, [X, [A|X]]), PD, PF).
1309principal_filter(_, _, []).
1310
1311% ?- principal_filter_set([a,b,c], [a], X).
1312%@ X = [[[a, c], [a, c, b], [a, b], [a]]].
1313% ?- principal_filter_set([a,b,c], [a,b], X).
1314%@ X = [[[a, c], [a, c, b], [a, b], [a]], [[b, c], [b, c, a], [b, a], [b]]].
1315
1316principal_filter_set(D, As, PFs):-
1317	maplist(pred(D, [A, PF]:- principal_filter(D, A, PF)),
1318		As, PFs).
1319
1320% ?- scm(math).
1321
1322%?- inverse_image([(a,b),(b,c)], [a], M).
1323%?- inverse_image([(a,b),(b,c)], [c], M).
1324
1325inverse_image(F, R, S):- inverse_image(F, R, [], S).
1326
1327inverse_image([(A,B)|C], R, S, S0):- memberchk(B, R), !,
1328	inverse_image(C, R, [A|S], S0).
1329inverse_image([_|C], R, S, S0):- !, inverse_image(C, R, S, S0).
1330inverse_image([], _, S, S).
1331
1332
1333%?-  image([(a,1),(b,2)], [b,a], I).
1334image(F, D, I):- fun_image(F, D, [], I).
1335
1336fun_image([(A,B)|C], D, I, I0):- memberchk(A, D), !,
1337	fun_image(C, D, [B|I], I0).
1338fun_image([_|C], D, I, I0):- !, fun_image(C, D, I, I0).
1339fun_image([], _, I, I).
1340
1341%?- b_subset(X, [a,b]).
1342
1343b_subset([], []).
1344b_subset(S, [A|R]) :- b_subset(S0, R),
1345	(	S = S0
1346	;	S = [A|S0]
1347	).
1348
1349%
1350% ?- time(repeat(10000, eval(list::((set::pow(1..2))+(set::pow(1..2))), X))).
1351% % 890,003 inferences, 0.150 CPU in 0.157 seconds (96% CPU, 5933353 Lips)
1352% true.
1353% ?- time(repeat(10000, eval(ss, list::((set::pow(1..2))+(set::pow(1..2))), X))).
1354% % 760,003 inferences, 0.140 CPU in 0.134 seconds (105% CPU, 5428593 Lips)
1355% true.
1356% ?- time(repeat(100000, eval(ss, list::((set::pow(1..2))+(set::pow(1..2))), X))).
1357% % 7,600,003 inferences, 1.310 CPU in 1.311 seconds (100% CPU, 5801529 Lips)
1358% true.
1359
1360% ?- set(pow([1]+[2]), V).
1361% V = [[],[1],[2],[1,2]].
1362% ?- eval(set::pow([]), X).
1363% ?- set(cup([[1,2],[3,4]]), X).
1364% X = [1,2,3,4]
1365% ?- set(cap([[1,2],[2,3]]), X).
1366% X = [2]
1367% ?-set(cup(pow([1,2,3])), X).
1368% X = [1,2,3]
1369% ?- set(cap(pow([1,2,3])), X).
1370% X = []
1371% ?- set(cup(map([pow([1,2]), pow([a,b])])), R).
1372% R = [[],[1],[1,2],[2],[a],[a,b],[b]]
1373
1374bigcap([A|B], C):- bigcap(B, A, C).
1375
1376bigcap(_, [], []):-!.
1377bigcap([], A, A):-!.
1378bigcap([X|Y],A,B):- intersection(X, A, A0),
1379	bigcap(Y, A0, B).
1380
1381%
1382poly_in(Ps, X, Y, A, B, Q):- member(P, Ps), call(P, X, Y, A, B, Q), !.
1383poly_in(_, X, Y, [], [], true):- poly_in(X, Y).
1384
1385% ?- poly_in((a-b)^(1+3), R).
1386% R = (a+ -1*b)^4
1387% Converting polynomials into internal canonical forms.
1388
1389poly_in(-X,A1):-!,poly_in(-1*X,A1) .
1390poly_in(X-Y,A1):-!,poly_in(X+ -1*Y,A1) .
1391poly_in(X+Y,A1+A2):-!,poly_in(X,A1),poly_in(Y,A2) .
1392poly_in(X*Y,A1*A2):-!,poly_in(X,A1),poly_in(Y,A2) .
1393poly_in(X/Y,A1/A2):-!,poly_in(X,A1),poly_in(Y,A2) .
1394poly_in(X^N,A1^N0):-N0 is N,!,poly_in(X,A1) .
1395poly_in(\$(X),X):-! .
1396poly_in(X mod N,A1 mod N0):-N0 is N,!,poly_in(X,A1) .
1397
1398:-meta_predicate each(1,2,?,?) . 1399
1400each(_A1,_A2,X,X):-(\+compound(X),!),! .
1401each(T,D,X,A1):-(call(T,X),!,X=..[F|Xs]),!,(maplist(each(T,D),Xs,A2),rev(=..,F.A2,A3)),call(D,A3,A1) .
1402each(A1,A2,X,A3):-X=..[F|Xs],!,maplist(each(A1,A2),Xs,A4),rev(=..,F.A4,A3) .
1403ari_val(_A1,A,A):-number(A),! .
1404ari_val(Assoc,A,B):-memberchk(A=B,Assoc),! .
1405ari_val(A1,A*B,A2):-!,(ari_val(A1,A,A3),ari_val(A1,B,A4)),rev(is,A3*A4,A2) .
1406ari_val(A1,A+B,A2):-!,(ari_val(A1,A,A3),ari_val(A1,B,A4)),rev(is,A3+A4,A2) .
1407ari_val(A1,A-B,A2):-!,(ari_val(A1,A,A3),ari_val(A1,B,A4)),rev(is,A3-A4,A2) .
1408ari_val(A1,A^B,A2):-!,(ari_val(A1,A,A3),ari_val(A1,B,A4)),rev(is,A3^A4,A2) .
1409ari_val(A1,-B,A2):-!,ari_val(A1,B,A3),rev(is,-A3,A2) .
1410ari_val(A1,+B,A2):-!,ari_val(A1,B,A3),rev(is,+A3,A2) .
1411si(E,Z):-Z is E .
1412
1413polish_codes(com(F,A),A1):-!,(polish_codes(A,A2),maplist(polish_codes,A2,A3)),phrase(F,A3,A1) .
1414polish_codes([A1|A2],A3):-!,maplist(polish_codes,[A1|A2],A3) .
1415
1416% combinations
1417% ?- scm(ring).
1418% ?- combi(4,3, X).
1419%@ X = 4 .
1420% ?- combi(20,10, X).
1421% ?- fact(20, X).
1422
1423% nCr(N, R, C) <=>
1424% Let X, R be a set  s.t. |X| = N, a number R, respectively,
1425% then C =  |{Y in powerset(X) | |Y| = R}|
1426
1427% ?- time(math:nCr(200,100, R)).
1428%@ % 713 inferences, 0.000 CPU in 0.000 seconds (85% CPU, 4512658 Lips)
1429%@ R = 90548514656103281165404177077484163874504589675413336841320 .
1430% ?- time(math:nCr(10000,5000, R)).
1431
1432combi(N, R, C) :- nCr(N, R, C).
1433
1434nCr(N, R, C):-  N >= 0, R>=0, NR is N - R, NR>=0, !,
1435	(  R < NR -> combination(N, R, NR, C)
1436	;  combination(N, NR, R,C)
1437	).
1438nCr(_,_,0).
1439
1440combination(N, R, S, C):- S1 is S+1,
1441	nat_prod(S1, N, 1, M),
1442	fact(R, FR),
1443	C is M rdiv FR.
1444
1445nat_prod(K, N, X, X):-  K>N, !.
1446nat_prod(K, N, X, Y):-  X0 is K*X, K1 is K+1,
1447	nat_prod(K1, N, X0, Y).
1448
1449%
1450fact(N, F):- N >= 0, factorial(N, 1, F).
1451%
1452factorial(N, F):- N >= 0, factorial(N, 1, F).
1453%
1454factorial(0, X, X).
1455factorial(N, X, Y):- N > 0, X0 is N*X,
1456	succ(N1, N),
1457	factorial(N1, X0, Y).
1458
1459% ?- scm(math).
1460% ?- binomial(3,X).
1461% ?- binomial(10,X).
1462
1463% found incomplete [2019/03/29]
1464% % ?- binomial(3,X).
1465% binomial(N, E) :-
1466% 	do(for(0..N, pred(N, [J, B, [C|B]]:-combi(N, J, C))), [], E).
1467
1468% %  ?- binomial(10, [a,b], E), maplist(val([a=(1/2), b=(1/2)]), E, F).
1469% binomial(N, [P,Q], E) :-
1470%  	for(0..N, pred((P,Q,N), [J, B, [C*(P^J)*(Q^K)|B]]:-(combi(N, J, C),
1471%  					 K is N-J)),
```