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