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