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
  207poly(X+Y,A1):-!,(poly(X,A2),poly(Y,A3)),poly_add(A2,A3,A1) .
  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),
  281	poly_add([A], P1, Q).
  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),
  337	poly_add(C, S, V),
  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):- !,
  358	head(N, S, S0),
  359	reverse(W, S0, X).
  360poly_mul_bd([A|As], B, N, W, S, X):-
  361	poly_scalar_bd(A, B, N, C),
  362	poly_add(C, S, V),
  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%
  368head(_, [], []):-!.
  369head(0, _, []):-!.
  370head(N, [X|Y], [X|H]):- succ(N0, N),
  371	head(N0, Y, H).
  372
  373poly_add(poly(E0,X0), poly(E1,X1), poly(E, X)):- poly_add(E0,X0,E1,X1,E,X).
  374poly_add(poly(E0,X0), X1, poly(E, X)):- poly_add(E0,X0, 0,X1,E,X).
  375poly_add(X0, poly(E1,X1), poly(E, X)):- poly_add(0,X0,E1,X1,E,X).
  376poly_add(X0, X1, X):- poly_add0(X0,X1,X).
  377
  378poly_add(E0,X0,E1,X1,E,X):-
  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),
  384	poly_add0(Y0, Y1, Y2),
  385	(maplist(=(0), Y2) -> X=[]; X=Y2).
  386
  387poly_add0([],X, X).
  388poly_add0(X,[], X).
  389poly_add0([X|Xs], [Y|Ys], [Z|Zs]):- Z is X+Y, poly_add0(Xs,Ys,Zs).
  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%
  434% ?- poly_add([1,0], [1,0], R).
  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%
  894graded_powerset(X, P):-
  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])