1:- module(poly, [bits/3, poly_in/2, poly_out/2
    2				]).    3
    4% Euler - Mersenne prime number  2147483647
    5:- use_module(pac(basic)).    6:- use_module(util(meta2)).    7:- use_module(pac('expand-pac')).    8% :- expects_dialect(pac).
    9term_expansion --> pac:expand_pac.
   10:- use_module(pac(op)).   11:- style_check(-singleton).   12
   13% gcd by Euclidean method of mutual division.
   14% greatest common diviser
   15% gcd(121,33,X) => X = 11
   16% ?- poly:gcd(123, 33, X).
   17% ?- poly:gcd(3, 10, X).
   18% ?- poly:gcd(6, 8, X).
   19% ?- poly:gcd(12, 8, X).
   20% ?- X is  gcd(10, 10).
   21% ?- 10^10 =:= gcd(10^10, 10^10).
   22
   23gcd(X, Y, Z):- Z is gcd(X, Y).
   24
   25% ?- poly:lcm(8, 12, X).
   26% ?- poly:lcm(-3, 4, X).
   27% ?- poly:lcm(6, -4, X).
   28% ?- poly:lcm(6, -4, X).
   29
   30lcm(X, Y, Z):- Z is X*Y // gcd(X, Y).
   31
   32% ?- poly:lcm(6, -4, X, A, B).
   33%@ X = -12,
   34%@ A = -2,
   35%@ B = 3.
   36lcm(X, Y, Z, U, V):- Z is X*Y // gcd(X, Y), U is Z//X, V is Z//Y.
   37
   38% ?- divmod(-10, 3, X, Y).
   39% ?- divmod(-10, -3, X, Y).
   40
   41% TODO
   42% euclid(X, Y, Q, R).
   43
   44% %  given x >= y => (a, b), s.t a*x +b*y= gcd(x, y))
   45% % nat_gcd(x, y, a, b, d) <=> a*x +b*y= d (= gcd(x, y))
   46% % ?- nat_gcd(12, 8, A, B, D).
   47% % A = 1,
   48% % B = -1,
   49% % D = 4.
   50% %
   51% nat_gcd(X, Y, A, B, D):- X < Y, !, nat_gcd_(Y, X, B, A, D).
   52% nat_gcd(X, Y, A, B, D):- nat_gcd_(X, Y, A, B, D).
   53
   54% nat_gcd_(X, Y, A, B, D):- Z is X mod Y,
   55% 	(Z==0 -> A = 0, B = 1, D = Y
   56% 	; nat_gcd_(Y, Z, K, L, D), A0 is X//Y,
   57% 	    A = L, B is K- L*A0
   58% 	).
   59
   60next_prime(X, N, P):-  member(A, X), 0 =:= N rem A, !,
   61	N0 is N + 1,
   62	next_prime(X, N0, P).
   63next_prime(_, N, N).
   64
   65% ?- trace, poly:make_poly_homogeneous(t, [1*[a^1]], PH).
   66% ?- poly:make_poly_homogeneous(t, [1*[a^1], 1*[b^1]], PH).
   67% ?- poly:make_poly_homogeneous(t, [1*[a^2], 1*[b^1]], PH).
   68
   69make_poly_homogeneous(Hvar, P,  PH):-
   70	maplist(pred([C*M, C*M0, T, D]:- mono_tail(M, M0, T, D)),
   71	       P, PH, Ts, Ds),
   72	list_max(Ds, N),
   73	maplist(pred([N, Hvar],
   74				 ([N, []]
   75					&
   76				 ([J, [Hvar^K]]:- K is N-J))),
   77		Ds, Ts).
   78%
   79make_poly_non_homogeneous(SW, X, Y):-
   80	maplist(pred([I*Vs, I*Us]:-
   81		    foldl(pred([0^_, U, U]  & [A, [A|U], U]), Vs, Us, [])),
   82		     X, Y0),
   83	sort_poly(SW, Y0, Y1),
   84	merge_poly(Y1, Y2),
   85	canonical_form(Y2, Y).
   86
   87% ?- poly: mono_tail([a^1, b^1], X, T, D).
   88mono_tail(M, X, T, D):- mono_tail(M, X, T, 0, D).
   89%
   90mono_tail([], U, U, D, D).
   91mono_tail([X^I|R], [X^I|U], T, D0, D):-
   92	D1 is D0 + I,
   93	mono_tail(R, U, T, D1, D).
   94
   95% ?- poly:list_max([1,2,3], X).
   96list_max(L,  M):- list_max(L, 0, M).
   97
   98%
   99list_max([], M, M).
  100list_max([I|R], M, M0):- I> M, !, list_max(R, I, M0).
  101list_max([_|R], M, M0):- list_max(R, M, M0).
  102
  103% ?-  poly:int_div(-13, -3, Q, R).
  104%@ Q = 4,
  105%@ R = -1.
  106
  107int_div(X, D, Q, R):- divmod(X, D, Q, R).
  108
  109%
  110is_sorted_spoly([], _).
  111is_sorted_spoly([_],_).
  112is_sorted_spoly([P, Q|Z], F):-
  113	P = [_*X|_]-[_*Y|_],
  114	Q = [_*U|_]-[_*V|_],
  115	mono_lcm(X,Y,XY),
  116	mono_lcm(U,V,UV),
  117	call(F, C, XY, UV),
  118	C\== (>),
  119	!,
  120	is_sorted_spoly([Q|Z], F).
  121is_sorted_spoly([P, Q|_], _):-
  122	P = [_*X|_]-[_*Y|_],
  123	Q = [_*U|_]-[_*V|_],
  124	mono_lcm(X,Y,XY),
  125	mono_lcm(U,V,UV),
  126	writeln(not_ordered(XY=X-Y, UV=U-V)).
  127
  128% ?- poly:sum_list_abs_two_halves([-1,2,-3,4,-5,6], S1, S2).
  129%@ S1 = 6,
  130%@ S2 = 15.
  131
  132sum_list_abs_two_halves(L, Sum1, Sum2):-
  133	length(L, N),
  134	N0 is N//2,
  135	sum_list_abs_count(N0, L, L0, 0, Sum1),
  136	sum_list_abs(L0, 0, Sum2).
  137
  138% Elegant way for sums  of abs value of integers in the list
  139
  140% sum_list_abs_two_halves(L, Sum1, Sum2):- maplist(abs, L, L0),
  141% 	length(L, N),
  142% 	N0 is N//2,
  143% 	length(Snd_half, N0),
  144% 	append(Fst_half, Snd_half, L0),
  145% 	sum_list(Fst_half, Sum1),
  146% 	sum_list(Snd_half, Sum2).
  147
  148% ?- poly:sum_list_abs_count(3, [-1, -3, -5, -6], L, 0, S).
  149%@ L = [-6],
  150%@ S = 9.
  151
  152sum_list_abs_count(0, L, L, X, X):-!.
  153sum_list_abs_count(N, [A|R], L, X, Y):- X0 is abs(A) + X,
  154	N0 is N-1,
  155	sum_list_abs_count(N0, R, L, X0, Y).
  156
  157% ?- poly:sum_list_abs([-1, -3, -5, -6], 0, S).
  158%@ S = 15.
  159sum_list_abs([], X, X).
  160sum_list_abs([A|R], X, Y):- X0 is abs(A)+ X,
  161	sum_list_abs(R, X0, Y).
  162
  163% ?- poly:remove_content_poly([-6*_, 4*_, -8*_], Cs).
  164% ?- poly:remove_content_poly([], Cs).
  165% ?- poly:remove_content_poly([1*_], Cs).
  166
  167remove_content_poly(_, [], []).
  168remove_content_poly(BL, [C*M|P], R):- 0 =:= (  C>>BL), !,
  169	once(remove_content_poly([C*M|P], R)).
  170remove_content_poly(_, P, P).
  171
  172%
  173remove_content_poly(L,  L0):- maplist(pred([C*M, C]), L, Cs),
  174	once(remove_content(Cs, _, Ds)),
  175	maplist(pred([_*M, C, C*M]), L, Ds, L0).
  176
  177% ?- poly:remove_content([-6, 4, -8], G, Cs).
  178%@ G = 2,
  179%@ Cs = [-3, 2, -4].
  180% ?- poly:remove_content([1, -6, 4, -8], G, Cs).
  181%@ Cs = [1, -6, 4, -8] .
  182% ?- poly:remove_content([-1, -6, -4, -8], G, Cs).
  183%@ G = 1,
  184%@ Cs = [-1, -6, -4, -8] .
  185
  186remove_content([], _, []).
  187remove_content(L, G, Cs):-
  188	sum_list_abs_two_halves(L, Sum1, Sum2),
  189	G1 is gcd(Sum1, Sum2),
  190	maplist(pred(G1, [C, Q, R]:- divmod(C, G1, Q, R)),
  191		L, Qs, Rs),
  192	gcd_list(Rs, G2),
  193	G is gcd(G1, G2),
  194	K is G1//G,
  195	maplist(pred([K,G], [Q, R, C]:- C is Q*K + R//G),
  196		Qs, Rs, Cs).
  197
  198
  199% ?- poly:gcd_list([6, 4, 8], R).
  200gcd_list([A|R], G):- gcd_list(R, A, G).
  201
  202gcd_list([], X, X).
  203gcd_list([A|B], X, Y):-   X0 is gcd(A, X),
  204	gcd_list(B, X0, Y).
  205
  206%
  207% sum_abs_half_half(0, L, B, B, A):- sum_abs_half_half(L, 0, A).
  208% sum_abs_half_half(N, [I|Is], B0, B, A):-  B1 is B0 + abs(I),
  209% 	N0 is N -1,
  210% 	abs_sum_half_half(N0, Is, B1, B, A).
  211
  212%
  213% sum_list_abs([], S, S).
  214% sum_list_abs[A|R], S0, S):- S1 is S0 + abs(A),
  215% 	sum_list_abs(R, S1, S).
  216
  217% ?- poly: ( poly_in(a, Y),  poly_out(Y, Z)).
  218% ?- poly: ( poly_in(a^2;b, Y), sort_polyset(total, Y, Y1),  poly_out([], Y1, Z)).
  219% ?- poly:exp_poly(a^2, P).
  220% ?- poly: ( poly_in(a^2;b, Y), sort_polyset(rev(total), Y, Y1),  poly_out([], Y1, Z)).
  221
  222poly_in --> flatten(;),
  223	maplist(exp_poly),
  224	remove([]).
  225
  226%
  227poly_in([], []) --> poly_in.
  228
  229poly_in(Order, Assoc) -->
  230	pred([Order, Assoc], [X, Y]:- preprocess(X, Y, Order, Assoc)),
  231	poly_in.
  232
  233%
  234poly_out	--> poly_out([]).
  235
  236poly_out(Assoc) -->
  237	maplist(pred(Assoc, [X, Y]:- postprocess(X, Y, Assoc))),
  238	maplist(poly_exp).
  239
  240%
  241print_polyset(X):- maplist(writeln, X).
  242
  243% ?- poly:s_poly([1*[x^1], 1*[a^1]],  [1*[y^1], 1*[b^1]], S, total).
  244%@ S = [1*[y^1, a^1], -1*[x^1, b^1]] .
  245
  246s_poly([C0*V0|P0], [C1*V1|P1], S, F):-
  247      mono_lcm(V0, V1, _, U0, U1),
  248      mul_mono_poly(C1*U0, P0, Q0),
  249      mul_mono_poly(C0*U1, P1, Q1),
  250      subtract_poly_poly(Q0, Q1, S0, F),
  251      canonical_form(S0, S).
  252
  253% ?- poly:s_poly_z([1*[x^1], 1*[a^1]],  [1*[y^1], 1*[b^1]], S, total).
  254%@ S = [1*[y^1, a^1], -1*[x^1, b^1]] .
  255% ?- poly:s_poly_z([4*[x^1], 1*[a^1]],  [6*[y^1], 1*[b^1]], S, total).
  256%@ S = [3*[y^1, a^1], -2*[x^1, b^1]] .
  257
  258s_poly_z([M0|P0], [M1|P1], S, F):-
  259      mono_lcm_z(M0, M1, _, U0, U1),
  260      mul_mono_poly(U0, P0, Q0),
  261      mul_mono_poly(U1, P1, Q1),
  262      subtract_poly_poly(Q0, Q1, S, F).
  263
  264% [2013/09/23]
  265% some tiny
  266canonical_form([]).
  267canonical_form([1*_|_]).
  268
  269%
  270canonical_form([],[]).
  271canonical_form([1*A|X], [1*A|X]):- !.
  272canonical_form([C*M|Ps], [1*M|Qs]):- C0 is 1 rdiv C,
  273	maplist(pred(C0, [D*M, D0*M]:- D0 is C0*D), Ps, Qs).
  274
  275% ?- poly:join_with_semicolon([a,b], R).
  276%@ R = (a;b).
  277% ?- poly:join_with_semicolon([0,0], R).
  278%@ R = 0.
  279
  280join_with_semicolon(X, Y):- remove(0, X, X0), join_with_semicolon_(X0, Y).
  281
  282join_with_semicolon_([], 0):-!.
  283join_with_semicolon_([X], X):-!.
  284join_with_semicolon_([X, Y], (X;Y)):-!.
  285join_with_semicolon_([X, Y|Z], (X;U)):- join_with_semicolon_([Y|Z], U).
  286
  287%
  288zero_poly([]).
  289unit_poly([1*[]]).
  290constant_poly([_*[]]).
  291
  292% ?- poly:exp_poly((a-b)^3, X), poly:poly_exp(X, Y).
  293% ?- poly:(poly_in((a-1)^100, X), poly_out(X, Y)).
  294% ?- poly:(poly_in((a-1)^3, X), sort_poly(total, X, Y0), poly_out(Y0, Y)).
  295% ?- poly:(poly_in((a-b)^3, X), sort_poly(total, X, Y0), poly_out(Y0, Y)).
  296% ?- poly:poly_exp([1*[a],1*[b],1*[c]], X).
  297% ?- poly:poly_exp([1*[a],1*[b],1*[c], -2*[]], X).
  298% ?- poly:poly_exp([(-3)*[x^1, y^3, z^2], (-3)*[x^1, y^3, z^2]], R).
  299% ?- A is -1 rdiv 3, poly:poly_exp([(-3)*[x^1, y^3], A*[x^1, y^2, z^1]], R).
  300
  301exp_rat(T, R):- rat(T, R0), rat_rat(R0, R).
  302
  303rat_poly(r(X,[1*[]]), X):-!.
  304rat_poly(r(X,[C*[]]), Y):- C\==0, C0 is 1 rdiv C,
  305	mul_scalar_poly(C0, X, Y).
  306
  307%
  308exp_poly(L=R, P):- !, exp_poly(L-R, P).
  309exp_poly([], []):- !.
  310exp_poly(T, P):- exp_rat(T, R), rat_poly(R, P).
  311
  312% ?- poly:exp_polyset(a;b;c, X).
  313%@ X = [[1*[a^1]], [1*[b^1]], [1*[c^1]]] .
  314exp_polyset --> exp_polyset(lexical).
  315
  316exp_polyset(Ord)  --> flatten(;),
  317	maplist(exp_poly),
  318	maplist(canonical_form),
  319	remove([]),
  320	sort_polyset(Ord).
  321
  322%
  323rat_rat(r(P, Q), r(P0, Q0)):- poly(P, P0, lexical),
  324	poly(Q, Q0, lexical).
  325
  326%
  327poly_exp([], 0):-!.
  328poly_exp([A|R], T) :-  mono_exp(A, A0),
  329	slim_mono_exp(A0, A1),
  330	foldl(pred(([M, U, V]:- mono_exp(M, M0),
  331			poly_exp(M0, U, V))),
  332	  R, A1, T).
  333
  334%
  335poly_exp(0*_, U, U).
  336poly_exp(1*M, U, U + M).
  337poly_exp(-1*M, U, U - M).
  338poly_exp(C*M, U, U + C0*M):- C>0, rdiv_slash(C, C0).
  339poly_exp(C*M, U, U - C1*M):- C0 is -C,  rdiv_slash(C0, C1).
  340poly_exp(0, U, U).
  341poly_exp(C, U, V):- number(C), !,
  342	(	C> 0
  343	->	rdiv_slash(C, C0),
  344		V =  U + C0
  345	;	C0 is -C,
  346		rdiv_slash(C0, C1),
  347		V = U -C1
  348	).
  349poly_exp(C, U, U + C).
  350
  351% ?- poly:mono_exp((-3)*[x^1, y^3, z^2], R).
  352%@ R = -3* (x*y^3*z^2) .
  353
  354mono_exp(C*M,  E):- maplist(pred([X^1, X] & [X, X]), M, M0),
  355	mono_exp_fold(C*M0, E).
  356
  357%
  358mono_exp_fold(C*[], C).
  359mono_exp_fold(C*[A|R], C*T) :- foldl(pred([U, X, X*U]), R, A, T).
  360
  361%
  362slim_mono_exp(1*M, M).
  363slim_mono_exp(M, M).
  364%
  365rdiv_slash(rdiv(X,Y), X/Y).
  366rdiv_slash(X, X).
  367
  368%%%%
  369rat_exp(R, T):- rat_poly(R, P), !, poly_exp(P, T).
  370rat_exp(r(X, Y), /(X0, Y0)):- poly_exp(X, X0), poly_exp(Y, Y0).
  371
  372% ?- poly: rat_pow( -1 rdiv 3, -3, R).
  373rat_pow(X, J, Q) :- J >= 0, !, rat_pow(X, J, 1, Q).
  374rat_pow(X, J, Q) :- J0 is -J, rat_pow(X, J0, 1, Q0),
  375	Q is 1 rdiv Q0.
  376
  377rat_pow(_, 0, P, P):-!.
  378rat_pow(X, J, P, Q):- J0 is J-1, P0 is P*X,
  379	rat_pow(X, J0, P0, Q).
  380
  381%
  382exp_exp --> exp_rat, rat_exp.
  383
  384
  385%  ?- poly:poly([1*[x^1]] - [1*[y^1]], P, lexical).
  386
  387
  388basic_poly(0,[]):-! .
  389basic_poly(X,[X*[]]):-rational(X),! .
  390basic_poly(X,[1*[X^1]]):-! .
  391
  392
  393% ?- poly:poly(3*x*y, X, total).
  394% ?- poly:poly(3*x*y, X).
  395% ?- poly:poly((3*x*y)^2, X).
  396% ?- poly:poly((x-y+z)^2, X).
  397
  398poly(X, Y) :- poly(X, Y, lexical).
  399
  400poly(0, [], _).
  401poly(X, X, _):- listp(X).
  402poly(X, [X*[]], _):- number(X).
  403poly(X - Y, Z, F):- poly(X, X0, F), poly(Y, Y0, F),
  404	subtract_poly_poly(X0, Y0, Z, F).
  405poly(X + Y, Z, F):- poly(X, X0, F), poly(Y, Y0, F),
  406	add_poly_poly(X0, Y0, Z, F).
  407poly(X * Y, Z, F):- poly(X, X0, F), poly(Y, Y0, F),
  408	mul_poly_poly(X0, Y0, Z, F).
  409poly(X ^ N, Z, F):- poly(X, X0, F), power_poly(N, X0, Z, F).
  410poly('$Var'(X), [1*[X^1]], _).
  411poly(X, [1*[X^1]], _).
  412
  413% ?- poly:poly_stack([1+(2*3)], [], X, lexical).
  414% ?- poly:poly_stack([(1+2)*3], [], X, lexical).
  415% ?- poly:poly_stack([a -b +c ], [], X, lexical).
  416
  417% ?- once(time(poly:poly_stack([(a -b +c)^50 ], [], X, lexical))), fail; true.
  418%@ % 4,350,956 inferences, 0.785 CPU in 0.848 seconds (93% CPU, 5542196 Lips)
  419%@ true.
  420% ?- once(time(poly:non_modular_poly( (a - b + c)^50, R, lexical))), fail; true.
  421% ?- poly:poly_stack([(a -b +c)^10 ], [], X, lexical).
  422
  423poly_stack([],  D, D, _).
  424poly_stack(['$OPR'(F/N)|R],  D,  E, SW):- !, apply_stack(N, F, D, D0, SW),
  425	poly_stack(R, D0, E, SW).
  426poly_stack([X|R],  D, E, SW):-
  427	poly_stack(X, R, D, R0, D0),
  428	poly_stack(R0, D0, E, SW).
  429
  430%
  431poly_stack(X+Y, R, D, [X, Y, '$OPR'((+)/2)|R], D).
  432poly_stack(X-Y, R, D, [X, Y, '$OPR'((-)/2)|R], D).
  433poly_stack(X*Y, R, D, [X, Y, '$OPR'((*)/2)|R], D).
  434poly_stack(+X, R, D, [X|R], D).
  435poly_stack(-X, R, D, [X, '$OPR'((-)/1)|R], D).
  436poly_stack(X^N, R, D, [X, '$OPR'(pow(N)/1)|R], D).
  437poly_stack('$Var'(N), R, D, R, [[V]|D]):- !, make_var(N, V).
  438poly_stack(I, R, D, R, [[C]|D]):- integer(I), !, make_constant(I, C).
  439poly_stack(X, R, D, R, [[V]|D]):- make_var(X, V).
  440
  441%
  442make_var(X, 1*[X^1]).
  443
  444make_constant(I, I*[]).
  445
  446%
  447apply_stack(2, F, [A, B|D], [C|D], SW):- !, apply_poly(F, B, A, C, SW).
  448apply_stack(1, F, [A|D], [C|D], SW):- apply_poly(F, A, C, SW).
  449
  450%
  451apply_poly(+, A, B, C, SW):- add_poly_poly(A, B, C, SW).
  452apply_poly(*, A, B, C, SW):- mul_poly_poly(A, B, C, SW).
  453apply_poly(-, A, B, C, SW):- subtract_poly_poly(A, B, C, SW).
  454
  455%
  456apply_poly(-, A, C, _):- prefix_minus_poly(A, C).
  457apply_poly(pow(N), A, C, SW):- power_poly(N, A, C, SW).
  458
  459% poly(0, [], _).
  460% poly(X, X, _):- listp(X).
  461% poly(X, [X*[]], _):- number(X).
  462% poly(X - Y, Z, F):- poly(X, X0, F), poly(Y, Y0, F),
  463% 	subtract_poly_poly(X0, Y0, Z, F).
  464% poly(X + Y, Z, F):- poly(X, X0, F), poly(Y, Y0, F),
  465% 	add_poly_poly(X0, Y0, Z, F).
  466% poly(X * Y, Z, F):- poly(X, X0, F), poly(Y, Y0, F),
  467% 	mul_poly_poly(X0, Y0, Z, F).
  468% poly(X ^ N, Z, F):- poly(X, X0, F), power_poly(N, X0, Z, F).
  469% poly('$Var'(X), [1*[X^1]], _).
  470% poly(X, [1*[X^1]], _).
  471
  472% ?- poly:modular_poly( 1, R,  lexical, 5).
  473% ?- poly:modular_poly( 2*3*7, R,  lexical, 5).
  474% ?- poly:modular_poly( (x+y)^3, R,  lexical, 5).
  475% ?- poly:modular_poly( (x+y)^100, R,  lexical, 5).
  476% ?- poly:modular_poly( (x+y)^100, R, total, 10).
  477% ?- poly:modular_poly( (x-y)^10, R, total, 3).
  478% ?- poly:poly(x-y+z, R, total), poly:poly(R^2, S, total).
  479% ?- poly:poly(x+y+z, R, total), poly:power_poly(R^100, S, total).
  480% ?- poly:poly(x+y+z, R, total), poly:power_poly(70, R, S, total).
  481% ?- poly:poly(x+y+z, R, total), poly:modular_poly(R^100, S, total, 3).
  482% ?- poly:poly(x+y+z, R, total), poly:modular_poly(R^200, S, total, 11).
  483% ?- poly:modular_poly(x^3, Y, total, 3).
  484% ?- time(poly:modular_poly( (x-y+z)^100, R, total, 3)).
  485% ?- time(poly:modular_poly( (x-y+z)^100, R, lexical, 2)).
  486% ?- time(poly:modular_poly( (x-y+z)^500, R, lexical, 2)).
  487% ?- time(poly:modular_poly( (5*x-7*y)^100, R, total, 3)).
  488% ?- listing(poly:modular_poly).
  489
  490
  491modular_poly(X-Y,C,Ord,G):-!,(modular_poly(X,A,Ord,G),modular_poly(Y,B,Ord,G)),subtract_poly_poly(A,B,AB,Ord),poly_galois(G,AB,C) .
  492modular_poly(X+Y,C,Ord,G):-!,(modular_poly(X,A,Ord,G),modular_poly(Y,B,Ord,G)),add_poly_poly(A,B,AB,Ord),poly_galois(G,AB,C) .
  493modular_poly(X*Y,C,Ord,G):-!,(modular_poly(X,A,Ord,G),modular_poly(Y,B,Ord,G)),mul_poly_poly(A,B,AB,Ord),poly_galois(G,AB,C) .
  494modular_poly(X^N,C,Ord,G):-!,modular_poly(X,A,Ord,G),galois_power_poly(N,A,B,Ord,G),poly_galois(G,B,C) .
  495modular_poly(-X,A1,A2,A3):-!,modular_poly(X,A4,A2,A3),prefix_minus_poly(A4,A1) .
  496modular_poly('$Var'(I),[1*[I^1]],A1,A2):-! .
  497modular_poly(0,[],A1,A2):-! .
  498modular_poly(X,X,A1,A2):-is_list(X),! .
  499modular_poly(X,[Y*[]],A1,G):-(integer(X),galois(G,X,Y)),! .
  500modular_poly(X,[1*[X^1]],A1,A2):-! .
  501
  502
  503% ?- poly:modular_poly_calc(1, Y, lexical, 3).
  504% ?- poly:modular_poly_calc(5^10000, Y, lexical, 3).
  505% ?- poly:modular_poly_calc((a-b)^5, Y, lexical, 2).
  506% ?- poly:modular_poly_calc((a-b-c)^50, Y, lexical, 3).
  507
  508modular_poly_calc(X, Y, Ord, Prim):- exp_poly(X, Poly),
  509	sort_poly(Ord, Poly, Poly0),
  510	poly_galois(Prim, Poly0, Poly1),
  511	modular_poly_list(Poly1, Poly2, Ord, Prim),
  512	poly_exp(Poly2, Y).
  513
  514% ?- poly:poly_calc(1, Y, lexical).
  515% ?- poly:poly_calc(5^10000, Y, lexical).
  516% ?- poly:poly_calc((a-b)^5, Y, lexical).
  517% ?- poly:poly_calc((a-b-c)^50, Y, lexical).
  518
  519poly_calc(X, Y, Ord):- non_modular_poly(X, Y0, Ord), poly_exp(Y0, Y).
  520
  521% ?- poly:modular_poly_list([1,2], X, lexical, 3).
  522% ?- listing(modular_poly_list).
  523
  524modular_poly_list([],[],A1,A2):-! .
  525modular_poly_list([X|Y],[X|Y],A1,A2):-! .
  526modular_poly_list(X-Y,C,Ord,G):-!,(modular_poly_list(X,A,Ord,G),modular_poly_list(Y,B,Ord,G)),subtract_poly_poly(A,B,AB,Ord),poly_galois(G,AB,C) .
  527modular_poly_list(X+Y,C,Ord,G):-!,(modular_poly_list(X,A,Ord,G),modular_poly_list(Y,B,Ord,G)),add_poly_poly(A,B,AB,Ord),poly_galois(G,AB,C) .
  528modular_poly_list(X*Y,C,Ord,G):-!,(modular_poly_list(X,A,Ord,G),modular_poly_list(Y,B,Ord,G)),mul_poly_poly(A,B,AB,Ord),poly_galois(G,AB,C) .
  529modular_poly_list(X^N,C,Ord,G):-!,modular_poly_list(X,A,Ord,G),galois_power_poly(N,A,B,Ord,G),poly_galois(G,B,C) .
  530modular_poly_list(-X,A1,A2,A3):-!,modular_poly_list(X,A4,A2,A3),prefix_minus_poly(A4,A1) .
  531
  532
  533% ?- poly:poly_galois(3, [4*[x^1], 5*[]], R).
  534poly_galois(G, B, C):-
  535	foldl(pred(G, (	[C*X, [D*X|U], U]:- galois(G, C, D), D\==0)
  536				&	[_, U, U]), B,  C, []).
  537
  538% ?- time(poly:non_modular_poly( (5*x-7*y)^100, R, total)).
  539% ?- time(poly:non_modular_poly( (5*x-7*y)^100, R, lexical)).
  540% ?- time(poly:non_modular_poly( (-(x) + -y)^2, R, lexical)).
  541
  542
  543non_modular_poly(X-Y,C,Ord):-!,(non_modular_poly(X,A,Ord),non_modular_poly(Y,B,Ord)),subtract_poly_poly(A,B,C,Ord) .
  544non_modular_poly(X+Y,C,Ord):-!,(non_modular_poly(X,A,Ord),non_modular_poly(Y,B,Ord)),add_poly_poly(A,B,C,Ord) .
  545non_modular_poly(X*Y,C,Ord):-!,(non_modular_poly(X,A,Ord),non_modular_poly(Y,B,Ord)),mul_poly_poly(A,B,C,Ord) .
  546non_modular_poly(X^N,C,Ord):-!,non_modular_poly(X,A,Ord),power_poly(N,A,C,Ord) .
  547non_modular_poly(-X,A1,A2):-!,non_modular_poly(X,A3,A2),prefix_minus_poly(A3,A1) .
  548non_modular_poly(0,[],A1):-! .
  549non_modular_poly('$Var'(N),[1*[N^1]],A1):-! .
  550non_modular_poly(X,[X*[]],A1):-integer(X),! .
  551non_modular_poly(X,X,A1):-listp(X),! .
  552non_modular_poly(X,[1*[X^1]],A1):-! .
  553
  554
  555% ?-gblex:rat((1+x)/y, R).
  556%@ R = r((1*1+x*1)*1, 1*1*y).
  557% ?- poly:rat((a/b)+1, X).
  558%@ X = r(a*1*1+1* (1*b), 1*b*1).
  559
  560
  561rat(+X,A1):-!,rat(X,A1) .
  562rat(-X,A1):-!,rat(-1*X,A1) .
  563rat(X-Y,A1):-!,rat(X+ -1*Y,A1) .
  564rat(X+Y,r(X0*Y1+Y0*X1,X1*Y1)):-!,rat(X,r(X0,X1)),rat(Y,r(Y0,Y1)) .
  565rat(X*Y,r(X0*Y0,X1*Y1)):-!,rat(X,r(X0,X1)),rat(Y,r(Y0,Y1)) .
  566rat(X/Y,r(X0*Y1,X1*Y0)):-!,rat(X,r(X0,X1)),rat(Y,r(Y0,Y1)) .
  567rat(X^N,r(X0^N,X1^N)):-!,rat(X,r(X0,X1)) .
  568rat('$Var'(X),r([1*[X^1]],1)):-! .
  569rat(X,r(X,1)):-! .
  570
  571
  572indeterminates(X-Y, X0-Y0, A, B):-
  573	indeterminates(X, X0, A, A0),
  574	indeterminates(Y, Y0, A0, B).
  575indeterminates(X+Y, X0+Y0, A, B):-
  576	indeterminates(X, X0, A, A0),
  577	indeterminates(Y, Y0, A0, B).
  578indeterminates(X*Y, X0*Y0, A, B):-
  579	indeterminates(X, X0, A, A0),
  580	indeterminates(Y, Y0, A0, B).
  581indeterminates(X/Y, X0/Y0, A, B):-
  582	indeterminates(X, X0, A, A0),
  583	indeterminates(Y, Y0, A0, B).
  584indeterminates(X=Y, X0=Y0, A, B):-
  585	indeterminates(X, X0, A, A0),
  586	indeterminates(Y, Y0, A0, B).
  587indeterminates(X;Y, X0;Y0, A, B):-
  588	indeterminates(X, X0, A, A0),
  589	indeterminates(Y, Y0, A0, B).
  590indeterminates(X^Y, X0^Y, A, B):-
  591	indeterminates(X, X0, A, B).
  592indeterminates(-X, -X0, A, B):-
  593	indeterminates(X, X0, A, B).
  594indeterminates(X, X, A, A):- number(X).
  595indeterminates(X, '$Var'(N), [X-N|A], A).
  596
  597
  598% ?- poly:(preprocess(1+x, X, A), poly(X, Y), postprocess(Y, Z, A)).
  599% ?- poly:(preprocess(1+x, X, A), poly(X, Y), postprocess(Y, Z, A)).
  600% ?- poly:(preprocess((1+x)^100, X, A), poly(X, Y), postprocess(Y, Z, A)).
  601% ?- poly:preprocess(x, X, R).
  602
  603preprocess(X, Y, Assoc):- once(preprocess(X, Y, [], Assoc)).
  604
  605% ?-  poly:preprocess(x+y+a, X, [y,x], R).
  606%@ X = '$Var'(3)+'$Var'(2)+'$Var'(1),
  607%@ R = [a-1, x-3, y-2] .
  608
  609preprocess(X, Y, Order, DefaultAssoc):-
  610	indeterminates(X, Y, B, []),
  611	sort(B, B0),
  612	merge_assoc(B0, DefaultAssoc),
  613	merge_order_to_assoc(Order, DefaultAssoc).
  614
  615%
  616assoc_numbering(A):- length(A, N),
  617	numlist(1, N, L),
  618	zip_hyphen(_, L, A).
  619
  620% ?- poly:number_vars(a+b, [a-1, b-2], Y).
  621%@ Y = '$Var'(1)+'$Var'(2) .
  622% ?- poly:number_vars(a+b+c, [a-1, b-2], Y).
  623%@ Y = '$Var'(1)+'$Var'(2)+'$Var'(c) .
  624
  625number_vars(X, Assoc, Y):-
  626	indeterminates(X, Y, B, []),
  627	subtract_assoc(B, Assoc, B0),
  628	maplist(pred([X-X]), B0).    % pass through unknown new variabls.
  629
  630% ?- poly:subtract_assoc([a-X, b-Y, a-Z, c-U], [a-1], R).
  631%@ X = Z, Z = 1,
  632%@ R = [b-Y, c-U] .
  633subtract_assoc(X, Y, Z):-
  634	foldl(pred(Y,	  (	[A, U, U]  :- memberchk(A, Y))
  635			& (	[A, [A|U], U]) ),
  636		X, Z, []).
  637
  638%
  639postprocess(X, X, []).
  640postprocess(X, Y, Assoc):-
  641	maplist(pred(Assoc, [I*Vs, I*Us]:-
  642		    maplist(pred(Assoc, ( [P^J, Q^J]:- memberchk(Q-P, Assoc), !)
  643					&  [U, U] ),
  644				 Vs, Us)),
  645		     X, Y).
  646
  647%
  648merge_assoc([], []).
  649merge_assoc([X-N|R], [X-N|A]):- merge_assoc(R, R0, X, N),
  650	merge_assoc(R0, A).
  651%
  652merge_assoc([X-N|R], R0, X, N):- !, merge_assoc(R, R0, X, N).
  653merge_assoc(R, R,  _, _).
  654
  655%
  656zip_hyphen([], [], []).
  657zip_hyphen([A|B], [C|D], [A-C|R]):- zip_hyphen(B, D, R).
  658
  659% ?- poly:merge_order_to_assoc([d, c, a], [a-X, b-Y,c-Z]).
  660%@ X = 4,
  661%@ Y = 1,
  662%@ Z = 3.
  663
  664merge_order_to_assoc([], Assoc):- !, assoc_numbering(Assoc).
  665merge_order_to_assoc(L, Assoc):-  zip_hyphen(L, _, Bs),
  666	union(Assoc, Bs, Cs),
  667	assoc_numbering(Cs).
  668
  669%
  670list_size(X, L):- list_size(X, 0, L).
  671
  672list_size([], N, N).
  673list_size([_|R], N0, N):-  N1  is N0 + 1,
  674	list_size(R, N1, N).
  675
  676% ?- poly:sort_mono([a^1, c^1, b^1, c^1], R).
  677%@ R = [c^2, b^1, a^1].
  678
  679sort_mono(L, R) :- length(L, N),
  680	sort_mono(N, L, _, R).
  681
  682%
  683sort_mono(2, [X1, X2|L], L, R) :- !, sort2(X1, X2, R).
  684sort_mono(1, [X|L], L, [X]) :- !.
  685sort_mono(0, L, L, []) :- !.
  686sort_mono(N, L1, L3, R) :-
  687	N1 is N // 2,
  688	plus(N1, N2, N),
  689	sort_mono(N1, L1, L2, R1),
  690	sort_mono(N2, L2, L3, R2),
  691	merge_mono_mono(R1, R2, R).
  692
  693%
  694sort2(X^I, X^J, [X^K]):- !,  K is I+J.
  695sort2(X^I, Y^J, [X^I, Y^J]):- X @> Y, !.
  696sort2(U, V, [V, U]).
  697
  698% ?- poly: dim_mono_term([a^1, b^2], X).
  699%@ X = 3.
  700
  701dim_mono_term([], 0).
  702dim_mono_term([_^I|R], N):- dim_mono_term(R, N0), N is I+N0.
  703
  704% ?- poly: compare_total_order(C, [b^1, a^1], [c^1]).
  705%@ C = (>).
  706% ?- poly: compare_total_order(C, [b^1, a^1], [b^1, a^1]).
  707%@ C = (=) .
  708
  709compare_total_order(C, X, Y):- dim_mono_term(X, D0),
  710	dim_mono_term(Y, D1),
  711	compare(C, D0, D1),
  712	C \== (=),
  713	!.
  714compare_total_order(C, X, Y):- compare_mono_term(C, X, Y).
  715
  716% ?- poly: compare_mono_term(C, [b^1, a^1], [b^1, a^1]).
  717%@ C = (=) .
  718compare_mono_term(C, [X^J|R], [X^K|S]):- !, compare(C0, J, K),
  719	(	C0 == (=)
  720	->	compare_mono_term(C, R, S)
  721	;	C = C0
  722	).
  723compare_mono_term( C , [X^_|_], [Y^_|_]):- !,
  724	(	X @> Y
  725	->	C = (>)
  726	;	C = (<)
  727	).
  728compare_mono_term( = , [], []).
  729compare_mono_term( < , [], _).
  730compare_mono_term( > , _, []).
  731
  732%
  733compare_mono(lexical, C, _*M, _*N):- !, compare_mono_term(C, M, N).
  734compare_mono(total, C, _*M, _*N):- !, compare_total_order(C, M, N).
  735% compare_mono(rev(Order), C, X, Y):- compare_mono(Order, D, X, Y),
  736%  	( D == (<), C = (>);  D == (>), C = (<);  true ),
  737% 	!.
  738
  739%
  740merge_mono([X^I|R], [X^J|S]):-  merge_mono(R, R0, X, I, J),
  741	merge_mono(R0, S).
  742merge_mono([], []).
  743%
  744merge_mono([X^K|R], R0, X, I, J):- !,  I0 is I + K,
  745	merge_mono(R, R0, X, I0, J).
  746merge_mono(R, R, _, I, I).
  747%
  748merge_mono_mono(X, [], X):- !.
  749merge_mono_mono([], X, X):- !.
  750merge_mono_mono([X^I|R], [X^J|S], [X^K|T]):- !, K is I + J,
  751	merge_mono_mono(R, S, T).
  752merge_mono_mono([X^I|R], [Y^J|S], [Y^J|T]):-  X @< Y, !,
  753	merge_mono_mono([X^I|R], S, T).
  754merge_mono_mono([U|R], B, [U|T]):-  merge_mono_mono(R, B, T).
  755
  756% ?- poly:mul_mono_mono(3*[c^1, a^1], 2*[c^1, b^2], X).
  757%@ X = 6*[c^2, b^2, a^1].
  758mul_mono_mono(C*M, D*N, E*L):- merge_mono_mono(M, N, L), E is C*D.
  759
  760% ?- gblex:mul_mono_poly(1*[c^1], [1*[b^1], 1*[a^1]], X).
  761mul_mono_poly(M, P, P0):- maplist(mul_mono_mono(M),P, P0).
  762
  763% ?- poly:compare_poly(lexical, X, [1*[y^2, x^2]], [1*[y^1, x^2]]).
  764%@ X = (>).
  765% ?- poly:compare_poly(rev(lexical), X, [1*[y^2, x^2]], [1*[y^1, x^2]]).
  766%@ X = (<).
  767% ?- poly:compare_poly(total, X, [1*[y^2, x^2]],[1*[y^2, x^2]]).
  768%@ X = (=) .
  769% ?- poly:compare_poly(rev(total), X, [1*[y^2, x^1]],[1*[y^2, x^2]]).
  770%@ X = (>).
  771%
  772sort_polyset(F, X, Y):- predsort(compare_poly(F), X, Y).
  773
  774%
  775compare_poly(_, =, [], []).
  776compare_poly(_, <, [], _).
  777compare_poly(_, >, _, []).
  778compare_poly(F, C, [M|_], [N|_]):- compare_mono(F, C, M, N),
  779	C\== (=).
  780compare_poly(_, <, _, _).
  781
  782% ?- poly:merge_mono([c^1, a^1, a^1], X).
  783% ?- poly:merge_mono_mono([c^1, a^1], [d^1, b^2, a^1], X).
  784% ?- poly:sort_poly(total, [1*[y^1, x^1], 1*[x^3]], X).
  785% ?- poly:sort_poly(lexical, [1*[x^2], 1*[y^1]], X).
  786% ?- poly:sort_poly(total, [1*[y^1, x^1], 1*[x^3]], X).
  787%@ X = [1*[x^3], 1*[y^1, x^1]] .
  788
  789sort_poly(F, L, R) :- length(L, N),
  790	sort_poly(N, L, _, R, F).
  791
  792%
  793sort_poly(2, [X1, X2|L], L, R, F) :- !,
  794	sort2_poly(X1, X2, R, F).
  795sort_poly(1, [X|L], L, [X], _) :- !.
  796sort_poly(0, L, L, [], _) :- !.
  797sort_poly(N, L1, L3, R, F) :-
  798	N1 is N // 2,
  799	plus(N1, N2, N),
  800	sort_poly(N1, L1, L2, R1, F),
  801	sort_poly(N2, L2, L3, R2, F),
  802	merge_poly_poly(R1, R2, R, F).
  803
  804%
  805sort2_poly(U, V, R, F):- compare_mono(F, A, U, V),
  806	(	A == (<),  R = [V, U]
  807	;	A == (>),  R = [U, V]
  808	;	U = C*M,
  809		V = D*_,
  810		E is C + D,
  811		R = [E*M]
  812	).
  813
  814
  815% ?- poly:merge_polyset_polyset([],[], [], total).
  816% ?- poly:merge_polyset_polyset([[2*[]]],[[1*[]]], X, total).
  817% ?- poly:merge_polyset_polyset([[2*[a^2]]],[[1*[b^1]]], X, total).
  818% ?- poly:merge_polyset_polyset([[2*[a^2]]],[[1*[b^1]]], X, lexical).
  819%  is this pred necessary ?
  820merge_polyset_polyset(X, Y, Z, F):- merge_polyset_polyset(X, Y, Z, [], F).
  821
  822%
  823merge_polyset_polyset([], Y, Y, [], _).
  824merge_polyset_polyset(X, [], X, [], _).
  825merge_polyset_polyset([A|X], [A|Y], [A|Z], U, F):- !,
  826	merge_polyset_polyset(X, Y, Z, U, F).
  827merge_polyset_polyset([A|X], [B|Y], [A|Z], U, F):-
  828	compare_poly(F, C, A, B), C == (>), !,
  829	merge_polyset_polyset(X, [B|Y], Z, U, F).
  830merge_polyset_polyset([A|X], [B|Y], [B|Z], U, F):-
  831	merge_polyset_polyset([A|X], Y, Z, U, F).
  832
  833
  834% ?- poly:sort_poly([1*[a^1], 1*[b^1], 1*[a^1]], X, total).
  835%@ X = [1*[b^1], 2*[a^1]] .
  836% ?- poly:sort_poly([1*[a^2], 1*[b^1], 1*[a^1]], X, compare_total_order).
  837%@ X = [1*[a^2], 1*[b^1], 1*[a^1]] .
  838% ?- poly:sort_poly([1*[a^2], 1*[b^1], 1*[a^1]], X, compare_mono_term).
  839%@ X = [1*[b^1], 1*[a^2], 1*[a^1]] .
  840% ?- poly:sort_poly([1*[a^2], 1*[b^1], 1*[a^1]], X, compare_total_order).
  841%@ X = [1*[a^2], 1*[b^1], 1*[a^1]] .
  842% ?- poly:merge_poly([1*[b^1], 1*[b^1], 1*[a^1]], X).
  843% ?- gblex:merge_poly([1*[b^1], 1*[a^1]], X).
  844% ?- gblex:merge_poly([1*[b^1], (-1)*[b^1], 1*[a^1]], X).
  845%@ X = [1*[a^1]].
  846
  847merge_poly([], []).
  848merge_poly([C*M|R], S):- merge_poly(R, R0, M, C, C0),
  849	(	C0\==0
  850	->	S = [C0*M|S0]
  851	;	S0 = S
  852	),
  853	merge_poly(R0, S0).
  854%
  855merge_poly([C*M|R], R0, M, D, E):- !, D0 is D+C,
  856	merge_poly(R, R0, M, D0, E).
  857merge_poly(R, R, _, C, C).
  858
  859% ?- poly:merge_poly_poly([1*[b^1], 1*[a^1]],[1*[b^1], 1*[a^1]], X, lexical).
  860%@ X = [2*[b^1], 2*[a^1]] .
  861
  862merge_poly_poly([], X, X, _).
  863merge_poly_poly(X, [], X, _).
  864merge_poly_poly([M|U], X, P, F):- merge_poly_poly(M, X, X0, P, Q, F),
  865	merge_poly_poly(U, X0, Q, F).
  866
  867%
  868merge_poly_poly(M, [N|V], V0, P0, P, F):- compare_mono(F, A, M, N),
  869	merge_poly_poly(A, M, N, V, V0, P0, P, F).
  870merge_poly_poly(M, [], [], [M|P], P, _).
  871
  872%
  873merge_poly_poly(=, I*X, J*_, V, V, [K*X|P], P, _):-  K is I+J, K\== 0, !.
  874merge_poly_poly(=, _  , _  , V, V, P, P, _):- !.
  875merge_poly_poly(>, M  , N  , V, [N|V], [M|P], P, _).
  876merge_poly_poly(_, M  , N  , V, V0, [N|P], Q, F):-
  877	merge_poly_poly(M, V, V0, P, Q, F).
  878
  879% ?- poly:mul_mono_mono(1*[c^1, a^1], 2*[c^1, b^2], X).
  880%@ X = 2*[c^2, b^2, a^1].
  881prefix_minus_poly(X, Y):- maplist(pred([I*X, J*X]:- J is 0 - I), X, Y).
  882
  883mul_scalar_poly(C, P, Q):- maplist(mul_scalar_mono(C), P, Q).
  884
  885mul_scalar_mono(C, C0*M, C1*M) :- C1 is C0*C.
  886
  887% ?- poly:div_mono_mono(-1*[b^2, c^1, a^3], 1*[b^1, a^1], D).
  888%@ D = -1*[b^1, c^1, a^2].
  889% ?- poly:div_mono_mono_term([b^2, c^1, a^3], [b^1, a^1]).
  890
  891%
  892div_mono_mono_term(X, Y):- div_mono_mono_term(X, Y, _).
  893%
  894div_mono_mono_term(M, [], M):-!.
  895div_mono_mono_term([X^I|Ys], [X^J|Ys0], Zs):- !,
  896	( I == J
  897	->	div_mono_mono_term(Ys, Ys0, Zs)
  898	;	I> J,
  899		K is I-J,
  900		Zs=[X^K|Zs0],
  901		div_mono_mono_term(Ys, Ys0, Zs0)
  902	).
  903div_mono_mono_term([X^I|Ys], [X0^J|Ys0], [X^I|Zs]):- compare((>), X, X0),
  904	div_mono_mono_term(Ys, [X0^J|Ys0], Zs).
  905
  906%
  907div_mono_mono(C*M, C0*M0, C1*M1):- div_mono_mono_term(M, M0, M1),
  908	C1 is C rdiv C0.
  909
  910
  911% ?- poly:div_poly_mono([-1*[b^2, c^1, a^3]], 1*[b^1, a^1], Q, D).
  912%@ Q = [-1*[b^1, c^1, a^2]],
  913%@ D = [].
  914% ?- poly:div_poly_mono([-1*[b^2, c^1, a^3], 1*[b-2]], 1*[b^1, a^1], Q, D).
  915%@ Q = [-1*[b^1, c^1, a^2]],
  916%@ D = [1*[b-2]].
  917
  918div_poly_mono([], _, [], []).
  919div_poly_mono([H|P], M, [Q0|Qs], R):- div_mono_mono(H, M, Q0), !,
  920	div_poly_mono(P, M, Qs, R).
  921div_poly_mono([H|P], M, Q, [H|R]):- div_poly_mono(P, M, Q, R).
  922
  923
  924% ?- poly:merge_poly_poly([2*[a^1]], [1*[a^1]], X, total).
  925% ?- poly:add_poly_poly([2*[a^1]], [1*[a^1]], X, lexical).
  926
  927add_poly_poly(X, Y, Z, F):- merge_poly_poly(X, Y, Z, F).
  928
  929% ?- poly:subtract_poly_poly([3*[x^1]], [1*[y^1], 2*[x^1]], R, total).
  930%@ R = [-1*[y^1], 1*[x^1]] .
  931% ?- poly:subtract_poly_poly([1*[]], [1*[]], R, lexical).
  932%@ R = [] .
  933% more efficient then " minus  then  add ".  [2014/07/13]
  934
  935subtract_poly_poly([], X, Y, _):- !, prefix_minus_poly(X, Y).
  936subtract_poly_poly(X, [], X, _):- !.
  937subtract_poly_poly([M|U], X, P, F):-
  938	subtract_poly_poly(M, X, X0, P, Q, F),
  939	subtract_poly_poly(U, X0, Q, F).
  940
  941%
  942subtract_poly_poly(M, [N|V], V0, P0, P, F):- compare_mono(F, A, M, N),
  943	subtract_poly_poly(A, M, N, V, V0, P0, P, F).
  944subtract_poly_poly(M, [], [], [M|P], P, _).
  945
  946%
  947subtract_poly_poly(=, I*X, J*_, V, V, [K*X|P], P, _):-  K is I-J, K\== 0, !.
  948subtract_poly_poly(=, _,   _, V, V, P, P, _):- !.
  949subtract_poly_poly(>, M,   N, V, [N|V], [M|P], P, _):-!.
  950subtract_poly_poly(_, M,   I*Y, V, V0, [J*Y|P], Q, F):-  J is 0 - I,
  951	subtract_poly_poly(M, V, V0, P, Q, F).
  952
  953
  954		/*****************************************************************
  955		*    Smart multiplication on polynomials with multi-variables    *
  956		*     reducing add (sort/merge) operations.                      *
  957		*****************************************************************/
  958
  959% ?- poly: (mul_poly_poly([1*[b^1], 1*[a^1]],[1*[b^1], 1*[a^1]], X, total), mul_poly_poly(X, X, Y, total)).
  960%@ X = [1*[b^2], 2*[b^1, a^1], 1*[a^2]],
  961%@ Y = [1*[b^4], 4*[b^3, a^1], 6*[b^2, a^2], 4*[b^1, a^3], 1*[a^4]] .
  962
  963mul_poly_poly([], _, [], _).
  964mul_poly_poly([A|As], Bs, X, Ord):-
  965	mul_mono_poly(A, Bs, U),
  966	mul_poly_poly(As, Bs, U, V, X, V, Ord).
  967
  968%
  969mul_poly_poly([], _, U, U, X, X, _).
  970mul_poly_poly([A|As], Bs, U, V, X, Y, Ord):-
  971	mul_mono_poly_head(A, Bs, U, U0, X,  X0, Ord),
  972	mul_poly_poly(As, Bs, U0, V, X0, Y, Ord).
  973
  974%
  975mul_mono_poly_head(_, [], U, U, X, X, _):-!.
  976mul_mono_poly_head(A, Bs, U, V, X, Y, Ord):-
  977	mul_mono_poly(A, Bs, [C|Cs]),
  978	sweep_higher_terms(C, U, U0, X, Y, Ord),
  979	merge_poly_poly(Cs, U0, V, Ord).
  980
  981%
  982sweep_higher_terms(C, [R|Rs], Ps, X, Y, Ord):-
  983	compare_mono(Ord, A, C, R),
  984	(	A == (=)	->	add_coeff(C, R, CR),
  985						(  CR == [] -> Ps = Rs ; Ps = [CR|Rs] ),
  986						X = Y
  987	;	A == (>)	->	Ps = [C, R|Rs],
  988						X = Y
  989	;	X = [R|X0],
  990		sweep_higher_terms(C, Rs, Ps, X0, Y, Ord)
  991	).
  992sweep_higher_terms(C, [], [C], X, X, _).
  993
  994%
  995add_coeff(C*M, D*_, E*M):- E is C+D, E\==0, !.
  996add_coeff(_, _, []).
  997
  998% ?- poly:galois_power_poly(1, [1*[a^1]], R, total, 3).
  999% ?- poly:mul_poly_poly([1*[]], [1*[]], X, total).
 1000% ?- poly:galois_power_poly(3, [1*[a^3]], R, total, 3).
 1001%@ R = [1*[a^9]] .
 1002% ?- poly:galois_power_poly(3, [1*[a^3], 1*[a^2]], R, total, 3).
 1003% ?- poly:galois_power_poly(3, [1*[a^3], 1*[a^2]], R, total, 3).
 1004
 1005galois_power_poly(N, X, Y, SW, G):- bits(N, [], Bs),
 1006	unit_poly(U),
 1007	galois_power_poly(Bs, X, U, Y, SW, G).
 1008
 1009%
 1010galois_power_poly([], _, X, X, _, _):-!.
 1011galois_power_poly([0|Bs], A, X, Y, SW, G):- !,
 1012	mul_poly_poly(X, X, X2, SW),
 1013	galois_poly(G, X2, GX2),
 1014	galois_power_poly(Bs, A, GX2, Y, SW, G).
 1015galois_power_poly([_|Bs], A, X, Y, SW, G):-
 1016	mul_poly_poly(X, X, X2, SW),
 1017	mul_poly_poly(A, X2, AX2, SW),
 1018	galois_poly(G, AX2, Y0),
 1019	galois_power_poly(Bs, A, Y0, Y, SW, G).
 1020
 1021
 1022%   ?- poly:galois_poly(3, [5*_, 6*_, 7*_], R).
 1023galois_poly(G, P, Q):- rem_poly(P, Q, G).
 1024
 1025%  ?- poly:rem_poly([5*_, 7*_], R, 3).
 1026%@ R = [2*_G11, 1*_G12].
 1027
 1028rem_poly([A*M|R], [B*M|R0], P):-  B is  A rem P, B\==0, !,
 1029 	rem_poly(R, R0, P).
 1030rem_poly([_|R], R0, P):- !, rem_poly(R, R0, P).
 1031rem_poly([], [], _).
 1032
 1033% ?- time(poly:power_poly(200, [1*[b^1], 1*[a^1]], R, lexical)).
 1034%@ % 245,396 inferences, 0.064 CPU in 0.074 seconds (86% CPU, 3857821 Lips)
 1035%@ R = [1*[b^200], 200*[b^199, a^1], 19900*[b^198, a^2], 1313400*[b^197, a^3], 64684950*[b^196, a^4], 2535650040*[b^195, ... ^ ...], 82408626300*[... ^ ...|...], 2283896214600*[...|...], ... * ...|...] .
 1036%@ % 311,590 inferences, 0.094 CPU in 0.109 seconds (86% CPU, 3305223 Lips)
 1037
 1038power_poly(N, X, Y, SW):- bits(N, [], Bs),
 1039	unit_poly(U),
 1040	power_poly(Bs, X, U, Y, SW).
 1041
 1042%
 1043power_poly([], _, X, X, _):-!.
 1044power_poly([0|Bs], A, X, Y, SW):- !,
 1045	mul_poly_poly(X, X, X2, SW),
 1046	power_poly(Bs, A, X2, Y, SW).
 1047power_poly([_|Bs], A, X, Y, SW):-
 1048	mul_poly_poly(X, X, X2, SW),
 1049	mul_poly_poly(A, X2, AX2, SW),
 1050	power_poly(Bs, A, AX2, Y, SW).
 1051
 1052% ?- poly: mono_lcm_z(1*[a^1], 1*[a^1], A, B, C).
 1053% ?- poly: mono_lcm_z(6*[a^1], 8*[a^1], A, B, C).
 1054% ?- poly: mono_lcm_z(-6*[b^2, a^1], 8*[b^1, a^2], A, B, C).
 1055% ?- poly: mono_lcm([b^2, a^1], [b^1, a^2], LCM, Y, Z).
 1056% ?- poly: mono_lcm([b^2, a^1], [b^1, a^2], LCM).
 1057% ?- poly: mono_lcm([b^2, a^1], [b^1, a^2], LCM).
 1058
 1059%
 1060mono_lcm(X, Y, Z):- mono_lcm(X, Y, Z, _, _).
 1061
 1062mono_lcm([], Y, Y, Y, []):-!.
 1063mono_lcm(X, [], X, [], X):-!.
 1064mono_lcm([V^I|X],[V^J|Y], [V^M|Z], P, Q):- !, compare(C, I, J),
 1065	( C == (=) -> M=I, mono_lcm(X, Y, Z, P, Q)
 1066	; C == (>) -> M=I, K is I-J, Q=[V^K|Q0], mono_lcm(X, Y, Z, P, Q0)
 1067	; M =J, K is J-I, P=[V^K|P0], mono_lcm(X, Y, Z, P0, Q)
 1068	).
 1069mono_lcm([U^I|X],[V^J|Y], [U^I|Z], P, [U^I|Q]):- V @< U,!,
 1070	mono_lcm(X, [V^J|Y], Z, P, Q).
 1071mono_lcm(X,[B|Y], [B|Z], [B|P], Q):- mono_lcm(X, Y, Z, P, Q).
 1072
 1073
 1074% ?- poly:mono_lcm_z(4*[], 6*[], A, B, C).
 1075%@ A = 12*[],
 1076%@ B = 3*[],
 1077%@ C = 2*[].
 1078mono_lcm_z(X, Y, U, V):- mono_lcm_z(X, Y, _, U, V).
 1079
 1080%
 1081mono_lcm_z(A*X, B*Y, C*Z, D*U, E*V):- mono_lcm(X, Y, Z, U, V),
 1082	lcm(A, B, C, D, E).
 1083
 1084% ?- poly:rev_member(X, [a,b,c]).
 1085rev_member(X, [_|R]):- rev_member(X, R).
 1086rev_member(X, [X|_]).
 1087
 1088% ?- poly:mod((1/2^21000)*2^2100, X, 2147483647).
 1089%@ X = 1024 .
 1090% ?- poly:mod((1 rdiv 2)^10000000, R, 2147483647).
 1091%@ R = 1 rdiv 1048576 .
 1092
 1093% ?- poly:galois(5, -rdiv(-3, -2), Z).
 1094% ?- poly:galois(5, -4, Z).
 1095% P = 2147483647.
 1096galois(P, rdiv(X, Y), Z):- !,
 1097	X0 is (X rem  P),
 1098	Y0 is (Y rem  P),
 1099	Z  is  X0 rdiv Y0.
 1100galois(P, -rdiv(X, Y), Z):- !, galois(P, rdiv(-X, Y), Z).
 1101galois(P, X, X0):-  X0 is (X rem P).
 1102
 1103% ?- poly:bits(3, [], X).
 1104%@ X = [1, 1] .
 1105bits(0, X, X).
 1106bits(1, X, [1|X]).
 1107bits(N, X, Y):- B is  N /\ 1,  N0 is N>>1,
 1108	bits(N0, [B|X], Y).
 1109
 1110% ?- poly:mod_exp(2147483647, 3, 2, X).
 1111% ?- poly:mod_exp(2147483647, 1000000000000000000000000000, 2, X).
 1112mod_exp(P, N, X, Y):- bits(N, [], Bs), !,  mod_exp(P, Bs, X, 1, Y).
 1113
 1114%
 1115mod_exp(P, [B|Bs], X, Y, Z):-!, Ysquare is Y*Y,
 1116	(	B == 0
 1117	->	Y0 is Ysquare
 1118	;	Y0 is X*Ysquare
 1119	),
 1120	galois(P, Y0, Y1),
 1121	mod_exp(P, Bs, X, Y1, Z).
 1122mod_exp(_, [], _,  X, X)