```    1:- module(gblex, [gb/2, gb/3, gb/4, rat/2, term_poly/2, poly_term/2]).    2:- use_module(util(misc)).    3% :- expects_dialect(pac).
4term_expansion --> pac:expand_pac.
5:- use_module(pac(op)).    6
7
8% ?- gb([c,f, a,b,d,e, y, x],a*x+b*y-c; d*x+e*y-f, R).
9%@ R = (y*e*a+ -1*y*d*b+d*c+ -1*a*f;x*a+y*b+ -1*c;x*d+y*e+ -1*f).
10% ?- gb([x, c,f, a,b,d,e, y],a*x+b*y-c; d*x+e*y-f, R).
11%@ R = (e*a*x-e*c-d*b*x+b*f;y*b+a*x-c;y*e+d*x-f) .
12% ?- gb([c,f,a,b,d,e,x,y],a*x+b*y-c; d*x+e*y-f, R).
13%@ R = (x*e*a-x*d*b-e*c+b*f;y*b+x*a-c;y*e+x*d-f) .
14% ?- gb(x - t; y - t, A).
15% ?- trace, gb(1, A).
16% ?- gblex:gb([x,y], -x + y + 1;  -y + 2 * x, X).
17% ?- gb(x^2+y^2-1; x+y-z, X).
18%@ X = (-1+y^2+x^2;z-y-x).
19
20% ?- gblex:gb_rev(x^2+y-3;x*y+y^2-1, X).
21% ?- gb(x^2+y^2-1; x+y-k, X).
22%@ X = (x^2-x*k+1/2*k^2-1/2;y+x-k)
23% ?- gblex:gb_rev(x^2+y^2-1; x+y-k, X).
24%@ X = (x^2+y^2-1;k-x-y) .
25% ?- gb(x+y-z;x-y+3*z-2; x^2+y^2+z^2-k^2, X).
26% ?- gblex:gb_ord([a,b,x,y, k], a*b-1; a*x-a*y; x-y-k, X).
27%@ X = (b*a-1;y-x;k) .
28% ?- F1=x^2-x*y-x-y^2+y, F2=x^3-x^2*y-x^2+x*y-y^3, gbtotal:gb([y, x], F1;F2, X).
29% ?- gblex:gb(x*y*z;  y-x; y-2, R).
30%@ R = (-2+x;-2+y;z).
31% ?- gblex:gb([x, i, u], x^2-3; i^2+1; (x*a+ i*b)^2-u, R).
32%@ R = (x^2-3;i^2+1;b^2-2*b*a*i*x-3*a^2+u).
33% ?- gblex:gb([x(j), i, u], x(j)^2-3; i^2+1; (x(j)*a+ i*b)^2-u, R).
34%@ R = (x(j)^2-3;i^2+1;b^2-2*b*a*i*x(j)-3*a^2+u).
35% ?- gblex:gb([x(j), i, u], x(j)^3;x(j)^2-3; i^2+1; (x(j)*a+ i*b)^2-u, R).
36% ?- gblex:gb([x(j), i, u], x(j)^2-3; i^2+1; (x(j)*a+ i*b)^2-u, R).
37%@ R = (x(j)^2-3;i^2+1;b^2-2*b*a*i*x(j)-3*a^2+u)
38% ?- gblex:gb((a+b+3+7)^3, R),write(R).
39%@ b^3+3*b^2*a+30*b^2+3*b*a^2+60*b*a+300*b+a^3+30*a^2+300*a+1000
40% ?- gblex:gb((a+b)^100, R),write(R).
41% ?- gb([i, a], i^2+1; i^3*a+i*b, X).
42% ?- gb(x-y; y-z; z-u; x-1; u,  X).
43% ?- gb([r(3), i, a, b], r(3)^2-1; i^2+1; i^3*a+r(3)^3*i*b, X).
44% ?- gb([r(3), i, b, a], r(3)^2-1; i^2+1; i^3*a+ i*b, X).
45% ?- gb(r(3)^2-1; i^2+1; i^3*a+ i*b, X).
46% ?- gb([ans,r(3),i], r(3)^2+1; i^2+1; (r(3)+i +a)^3; r(3)^4+ i^3-ans, X).
47%@  ans^2-2*ans+2
48% ?- gb([ans,r(3),i], r(3)^2+1; i^2+1; (r(3)+i +a)^3; r(3)^4+ i^4-ans, X).
49%@ ans-2.
50% ?- gb([ans,r(3),i], r(3)^2-3; i^2+1; (r(3)+i +a)^3; r(3)^4+ i^4-ans, X).
51%@  ans-10;
52% ?- trace, gb([], [], r(3)^2-1; i^2+1; i^3*a+ i*b, X).
53%@ X = (b-a;i^2+1;r(3)^2-1) .
54
55gb --> groebner_base.
56
57gb(Ord, X, Y):- math:indeterminates(X, V),
58	subtract(V, Ord, V0),
59	sort(V0, V1),
60	append(Ord, V1, U),
61	groebner_base(U, X, Y).
62
63gb(Pre, Post) --> gb_in, phrase(Pre), buchberger, phrase(Post), gb_out.
64
65%
66gb_rev(X, Y):- math:indeterminates(X, As),
67	reverse(As, Bs),
68	groebner_base(As, Bs, X, Y).
69
70%  ?- gblex:groebner_base(x-y; y-z; z-u; x-1; u,  X).
71groebner_base --> gb_in, !, buchberger, !, gb_out.
72%
73groebner_base(Ord) --> {sort(Ord, Ord0)}, groebner_base(Ord, Ord0).
74%
75groebner_base(Ord, Ord0, F, G):-
76	zip(Ord, Ord0, A),
77	subst(A, F, F0),
78	groebner_base(F0, G0),
79	zip(Ord0, Ord, B),
80	subst(B, G0, G).
81%
82gb_in --> flatten(;),
83	maplist(term_poly),
84	maplist(canonical),
85	remove([]).
86%
87gb_out --> predsort(compare_poly),
88	maplist(poly_term),
89	join_with_semicolon.
90
91% ?- gblex:linear([],a, R).
92% ?- gblex:linear([a],1, R).
93% ?- gblex:linear([a], a*b, R).
94% ?- gblex:linear([x,y], a*x+b*y+c;d*x+e*y+f, R), write(R).
95%@ R = (x*e*a^2-x*d*b*a-f*b*a+e*c*a;y*e*a-y*d*b+f*a-d*c) .
96% ?- E= (a*x+b*y+c;d*x+e*y+f), S = [ (a,1), (b,-1), (c,1), (d,1), (e,1), (f,3)], subst(S, E, E0), gblex:gb(E0, R).
97
98%  Printing polynomial, base, affine form.
99% ?- gblex:print_affine([x*[1*[]], y*[1*[]], [1*[]]]).
100% ?- gblex:print_affine_base([[x*[1*[]], y*[1*[]],[1*[]]], [x*[1*[]], y*[1*[]],[1*[]]]]).
101% ?- gblex:print_poly_base([[1*[x^1], 2*[y^2], 1*[]], [1*[x^1], 2*[y^2], 1*[]]]).
102% ?- gblex:print_poly_base_exp(1+x+2*y^2;1+x+2*y^2).
103
104%
105print_affine([]):- !, writeln(0).
106print_affine(B):- affine_to_poly(B, P),
107	poly_term(P, T),
108	writeln(T).
109
110print_affine_base(B):- maplist(print_affine, B).
111
112%
113print_poly([]):- !, writeln(0).
114print_poly(P):- poly_term(P, P0), writeln(P0).
115
116print_poly_base([]):- !, writeln(0).
117print_poly_base(B):- maplist(print_poly, B).
118
119%
120print_exp(E):- term_poly(E, E0), print_poly(E0).
121
122print_poly_base_exp(B):- flatten((;), B, B1),
123	maplist(print_exp, B1).
124
125%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
126% [2013/09/24] buchberger: Fourth version.
128% clearly correct based on Dickson's theorem on mono ideal.
129% The Dickson theorem implies that no infinite sequence of
131%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
132buchberger(X, Y):-  pairs(X, D),
133	buchberger(D, X, Y0),
134	normal_base(Y0, Y1),
135	maplist(canonical, Y1, Y).
136
137%
138buchberger([], G, G).
139buchberger([(A,B)|R], G, H):- s_poly(A, B, C),
141	canonical(D0, D),
142	buchberger(D, R, G, H).
143
144%
145buchberger(C, R, G, H)	:- zero_poly(C), !, buchberger(R, G, H).
146buchberger(C, _, _, [C]):- unit_poly(C), !.
147buchberger(C, R, G, H)	:- foldr(C^[U, X, [(C,U)|X]], G, R, R0),
148		buchberger(R0, [C|G], H).
149
150%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
151
152% ?- gblex:s_poly([1*[x^1], 1*[a^1]],  [1*[y^1], 1*[b^1]], S).
153%@ S = [1*[y^1,a^1],-1*[x^1,b^1]].
154
155s_poly([C0*V0|P0], [C1*V1|P1], S):-
156      mono_lcm(V0, V1, _, U0, U1),
157      mul_poly_poly([C1*U0], P0, Q0),
158      mul_poly_poly([C0*U1], P1, Q1),
159      subtract_poly_poly(Q0, Q1, S0),
160      canonical(S0, S).
161
162% ?- gblex:normal_base([[1*[x^2], 1*[x^1], 1*[a^1]], [1*[x^1]], [1*[x^2]]], R).
163% ?- gblex:normal_base([[1*[x^2]], [1*[]], [1*[x^2]]], R).
164normal_base(X, Y):- select(P, X, Y0), reduce_poly_one(P, Y0, Q), !,
165	normal_base([Q|Y0], Y).
166normal_base(X, Y):- zero_poly(Z), remove(Z, X, Y).
167
168% ?- gblex:normal_poly([1*[x^1], 1*[a^1]],  [[1*[x^1]], [1*[a^1]]], R).
169normal_poly(P, G, Q):- reduce_poly_one(P, G, P0), !,
170	normal_poly(P0, G, Q).
171normal_poly(P, _, P).
172
173% ?- gblex:normal_poly([1*[x^1], 1*[a^1]],  [[1*[x^2]], [1*[y^1]]]).
174normal_poly(P, G):- reduce_poly_one(P, G, _), !, fail.
175normal_poly(_, _).
176
177% ?- gblex:normal_poly_one([1*[b^1, a^1], 5*[b^1], 3*[a^1]], [1*[b^1]], U).
178%@ U = [1*[b^1],3 rdiv 5*[a^1]].
179%@ U = [1*[b^1],3 rdiv 5*[a^1]].
180
181normal_poly_one([1*P|R], [1*Q|S], U) :- div_mono_mono_v(P, Q, A),
182	mul_mono_poly((-1)*A, S, S0),
184	canonical(U0, U).
185
186% ?-  gblex:reduce_head([1*[x^2], 3*[a^1]], [[1*[x^1], 4*[]]], R).
187%@ R = [3*[a^1],16*[]].
188
189reduce_head([C*M|X], G, P):- member([_C*N|Y], G), % Assuming G is canonical.
190	div_mono_mono_v(M, N, Q),
191	!,
192	mul_mono_poly(C*Q, Y, Y0),
193	subtract_poly_poly(X, Y0, P0),
196
197
198% ?- gblex:reduce_mono(1*[x^1], [1*[x^1],2*[a^1]], P).
199%@ P = [-2*[a^1]].
200% ?- gblex:reduce_mono(1*[x^2], [1*[x^1],2*[a^1]], P).
201%@ P = [-2*[x^1,a^1]].
202%
203reduce_mono(M, [M0|Ms], P):- div_mono_mono(M, M0, D),
204	mul_mono_poly(D, Ms, P0),
205	mul_scalar_poly(-1, P0, P).
206
207% ?- gblex:reduce_poly_one([1*[x^2], 1*[x^1]], [[1*[x^1]]], R).
208%@ R = [1*[x^1]] .
209% ?- gblex:reduce_poly_one([1*[x^2], 1*[x^1]], [[1*[x^1],1*[a^1]]], R).
210%@ R = [1*[x^1,a^1],-1*[x^1]] .
211% ?- gblex:reduce_poly_one([1*[x^2], 1*[x^1]], [[1*[x^1]], [1*[a^1]]], R).
212%@ R = [1*[x^1]] .
213% ?- gblex:reduce_poly_one([1*[a^2]], [[1*[x^1]]], R).
214
215reduce_poly_one(P, Ps, R) :-
216	select(M, P, P0),
217	member(Q, Ps),
218	reduce_mono(M, Q, C),
220
221% ?- gblex:affine([x], [1*[x^1,b^1], 1*[x^1, a^1], [3*[]]], R).
222% ?- gblex:affine([x,y], [1*[y^1], 1*[x^1,b^1], 1*[x^1, a^1], [3*[]]], R).
223% ?- gblex:affine([x,y], [1*[y^1], 1*[x^1,b^1], 1*[x^1, a^1]], R).
224
225affine([], Poly, [Poly]):-!.
226affine([V|Vs], Poly, [V*A|L]):- extract_var(V, Poly, A, R), A\==[], !,
227	affine(Vs, R, L).
228affine([_|Vs], Poly, L):- affine(Vs, Poly, L).
229
230% ?- trace, gblex:affine([x,y], [1*[y^1], 1*[x^1,b^1], 1*[x^1, a^1]], R), gblex:affine_to_poly(R, S).
231%@    Call: (8) gblex:affine([x, y], [1*[y^1], 1*[x^1, b^1], 1*[x^1, a^1]], _G805) ? no debug
232%@ R = [x*[1*[b^1],1*[a^1]],y*[1*[]],[]],
233%@ S = [1*[x^1,b^1],1*[x^1,a^1],1*[y^1]] .
234
235affine_to_poly(L, X):- ditribute_var(L, X0),
236	append(X0, X1),
237	merge_poly(X1, X).
238
239ditribute_var([], []).
240ditribute_var([[]|R], X):- !, ditribute_var(R, X).
241ditribute_var([V*A|R], [B|X]):- mul_mono_poly(1*[V^1], A, B),
242	ditribute_var(R, X).
243ditribute_var([A|R], [A|Ps]):- ditribute_var(R, Ps).
244
245% ?- gblex:extract_var(x, [1*[x^1,b^1], 1*[x^1, a^1], [3*[]]], A, R).
246%@ A = [1*[b^1],1*[a^1]],
247%@ R = [[3*[]]].
248
249extract_var(_, [], [], []):-!.
250extract_var(X, [C*M|P], [C*M0|A], R):- select(X^1, M, M0), !,
251	extract_var(X, P, A, R).
252extract_var(X, [Q|P], A, [Q|R]):- extract_var(X, P, A, R).
253
254
255% Gauss's sweep method.
256linear(Vs) --> gb_in, !,
257	gauss_sweep(Vs),
258	maplist(affine_to_poly),
259	gb_out.
260
261% ?- gblex: elim_linear_var([x*[2*[a^1]], [1*[]]], x, [1*[a^1]], [[1*[]]], R).
262elim_linear_var(F, X, B, G0, H):- select(X*A, F, F0), !,
263	sweep(A, B, F0, G0, H).
264elim_linear_var(F, _, _, _, F).
265
266%
267sweep(A, B, P, Q, R):- div_poly_poly(A, B, J, []), !,
268	mul_scalar_poly(-1, J, J0),
269	mul_poly_affine(J0, Q, Q0),
271sweep(A, B, P, Q, R):- div_poly_poly(B, A, J, []), !,
272	mul_scalar_poly(-1, Q, Q0),
273	mul_poly_affine(J, P, P0),
275sweep(A, B, P, Q, R):-
276	mul_poly_affine(B, P, P0),
277	mul_scalar_poly(-1, A, A0),
278	mul_poly_affine(A0, Q, Q0),
280
281% ?- gblex: mul_scalar_affine(2, [x*[1*[b^1]], y*[1*[]]], R).
282%@ R = [x*[2*[b^1]],y*[2*[]]] .
283mul_scalar_affine(C, F, G):- mul_poly_affine([C*[]], F, G).
284
285% ?- gblex: mul_poly_affine([1*[a^1]],[x*[1*[b^1]], y*[1*[]]], R).
286%@ R = [x*[1*[b^1,a^1]],y*[1*[a^1]]] .
287mul_poly_affine(_, [], []).
288mul_poly_affine(A, [X*B|P], [X*C|Q]):- !, mul_poly_poly(A, B, C),
289	mul_poly_affine(A, P, Q).
290mul_poly_affine(A, [B|P], [C|Q]):- mul_poly_poly(A, B, C),
291	mul_poly_affine(A, P, Q).
292
293% ?- gblex: add_affine_affine([x*[1*[a^1]]],[x*[1*[b^1]], y*[1*[]]], R).
294%@ R = [x*[1*[b^1],1*[a^1]],y*[1*[]]].
295% ?- gblex: add_affine_affine([x*[-1*[a^1]]],[x*[-1*[b^1]], y*[1*[]]], R).
296%@ R = [x*[-1*[b^1],-1*[a^1]],y*[1*[]]].
297% ?- gblex: add_affine_affine([x*[-1*[a^1]]],[x*[1*[a^1]], y*[1*[]]], R).
298%@ R = [x*[],y*[1*[]]].
299
300add_affine_affine(P, Q, [X*C|R]):- select(X*A, P, P0), select(X*B, Q, Q0), !,
303add_affine_affine(P, Q, [X*A|R]):- select(X*A, P, P0), !,
305add_affine_affine(P, Q, [X*A|R]):- select(X*A, Q, Q0), !,
310
311% ?- E = (a*x+b*y+c;d*x+e*y+f), gblex:linear([x,y], E, R).
312%@ R = (x*e*a^2-x*d*b*a-f*b*a+e*c*a;y*e*a-y*d*b+f*a-d*c) .
313
314gauss_sweep(Vs, Ps, S):-  maplist(affine(Vs),  Ps, F),
315	gauss_sweep([], F, Up, Low),
316	append(Up, Low, S0),
317	remove([[]], S0, S).
318
319%
320gauss_sweep(Up, Low, Up0, Low0):- select(P, Low, Low1), select(X*A, P, P0), !,
321	maplist(pred([X, A, P0], [V, W]:- elim_linear_var(V, X, A, P0, W)), Low1, Low2),
322	maplist(pred([X, A, P0], [V, W]:- elim_linear_var(V, X, A, P0, W)), Up, Up1),
323	gauss_sweep([P|Up1], Low2, Up0, Low0).
324gauss_sweep(U, L, U, L).
325
326% [2013/09/23]
327% some tiny
328
329canonical([]).
330canonical([1*_|_]).
331
332canonical([],[]).
333canonical([1*A|X], [1*A|X]):- !.
334canonical([C*M|Ps], [1*M|Qs]):- C0 is 1 rdiv C,
335	maplist(pred(C0, [D*M, D0*M]:- D0 is C0*D), Ps, Qs).
336
337
338% ?- gblex:join_with_semicolon([a,b], R).
339%@ R = (a;b).
340% ?- gblex:join_with_semicolon([0,0], R).
341%@ R = 0.
342
343join_with_semicolon(X, Y):- remove(0, X, X0), join_with_semicolon_(X0, Y).
344
345join_with_semicolon_([], 0):-!.
346join_with_semicolon_([X], X):-!.
347join_with_semicolon_([X, Y], (X;Y)):-!.
348join_with_semicolon_([X, Y|Z], (X;U)):- join_with_semicolon_([Y|Z], U).
349
350% ?- join_term(+, [a,b,c,d], Z).
351% Z = a+b+c+d
352
353join_term(F, [X, Y|Z], U):- !, join_term(F, [Y|Z], U0),
354	U =.. [F,X,U0].
355join_term(_, [X], X).
356
357%
358zero_poly([]).
359unit_poly([1*[]]).
360constant_poly([_*[]]).```
?- gblex:`term_poly((a-b)^3, X)`, `poly_term(X, Y)`.
```  365term_rat(T, R):- rat(T, R0), rat_rat(R0, R).
366
367rat_poly(r(X,[1*[]]), X):-!.
368rat_poly(r(X,[C*[]]), Y):- C\==0, C0 is 1 rdiv C,
369	mul_scalar_poly(C0, X, Y).
370
371term_poly(L=R, P):- !, term_poly(L-R, P).
372term_poly([], [[]]):- !.
373term_poly(T, P):- term_rat(T, R), rat_poly(R, P).
374
375rat_rat(r(P, Q), r(P0, Q0)):- poly(P, P0),
376	poly(Q, Q0).
377
378% internal polynomial -> term
379% ?- gblex:poly_term([1*[a],1*[b],1*[c]], X).
380% X = a+b+c
381
382poly_term([], 0):-!.
383poly_term(P, T) :- maplist(mono_term, P, Ts),
384	join_term(+, Ts, T0),
385	reduce:term_rewrite(gblex:normal_exp_left, [], T0, T1),
386	reduce:term_rewrite(gblex:intro_prefix_minus_right, [], T1, T2),
387	reduce:term_rewrite(gblex:intro_infix_minus, [], T2, T).
388
389%
390mono_term(1*M, T) :- !, mono_term_(M, T).
391mono_term(C*M, C*T) :- mono_term_(M, T).
392
393mono_term_([], 1):-!.
394mono_term_(M, T):- join_term(*,  M, T).
395
396%%%%
397rat_term(R, T):- rat_poly(R, P), !, poly_term(P, T).
398rat_term(r(X, Y), /(X0, Y0)):- poly_term(X, X0), poly_term(Y, Y0).
399
400% ?- gblex: rat_pow( -1 rdiv 3, -3, R).
401% R = -27.
402rat_pow(X, J, Q) :- J >= 0, !, rat_pow(X, J, 1, Q).
403rat_pow(X, J, Q) :- J0 is -J, rat_pow(X, J0, 1, Q0),
404	Q is 1 rdiv Q0.
405
406rat_pow(_, 0, P, P):-!.
407rat_pow(X, J, P, Q):- J0 is J-1, P0 is P*X,
408	rat_pow(X, J0, P0, Q).
409
410%
411term_term --> term_rat, rat_term.
412
413%%%
414%  ?- gblex:poly([1*[x^1]] - [1*[y^1]], P).
415% P = [1*[x^1], -1*[y^1]].
416
417
418basic_poly(0,[]):-! .
419basic_poly(X,[X*[]]):-rational(X),! .
420basic_poly(X,[1*[X^1]]):-! .
421
422
423% ?- gblex:poly(3*x*y, X).
424%@ X = [3*[y^1,x^1]] .
425
426
427poly(0,[]):-! .
428poly([X|Y],[X|Y]):-! .
429poly(X,[X*[]]):-rational(X),! .
430poly(X-Y,A1):-!,(poly(X,A2),poly(Y,A3)),subtract_poly_poly(A2,A3,A1) .
432poly(X*Y,A1):-!,(poly(X,A2),poly(Y,A3)),mul_poly_poly(A2,A3,A1) .
433poly(X^N,A1):-!,poly(X,A2),power_poly(N,A2,A1) .
434poly(X,[1*[X^1]]):-! .
435
436
437% ?-gblex:rat((1+x)/y, R).
438%@ R = r((1*1+x*1)*1, 1*1*y)
439
440
441rat(+X,A1):-!,rat(X,A1) .
442rat(-X,A1):-!,rat(-1*X,A1) .
443rat(X-Y,A1):-!,rat(X+ -1*Y,A1) .
444rat(X+Y,A1):-!,(rat(X,A2),rat(Y,A3)),'pac#431'(A2,A3,A1) .
445rat(X*Y,A1):-!,(rat(X,A2),rat(Y,A3)),'pac#432'(A2,A3,A1) .
446rat(X/Y,A1):-!,(rat(X,A2),rat(Y,A3)),'pac#433'(A2,A3,A1) .
447rat(X^N,A1):-!,rat(X,A2),'pac#434'(A2^N,A1) .
448rat(X,r(X,1)):-! .
449
450'pac#431'(r(A,B),r(C,D),r(A*D+C*B,B*D)):-true .
451'pac#432'(r(A,B),r(C,D),r(A*C,B*D)):-true .
452'pac#433'(r(A,B),r(C,D),r(A*D,B*C)):-true .
453'pac#434'(r(B,C)^N, r(B^N,C^N)):-true .
454
455% :- bekind(rat,[]).
456% +X	=  X.
457% -X	=  (-1)* X.
458% X-Y	=  X + (-1)*Y.
459% X+Y	=  #([r(X0, X1), r(Y0, Y1)]\ `r(X0*Y1+Y0*X1, X1*Y1))@X@Y.
460% X*Y	=  #([r(X0, X1), r(Y0, Y1)]\ `r(X0*Y0, X1*Y1))@X@Y.
461% X/Y	=  #([r(X0, X1), r(Y0, Y1)]\ `r(X0*Y1, X1*Y0))@X@Y.
462% X^N	=  #(N^[r(X0, X1)]\ `r(X0^N, X1^N))@X.
463% X	=  `r(X, 1).
464% :-ekind.
465
466% ?- listing(gblex:rat).
467
468% ?- gblex:mono_merge([c^1, a^1], [d^1, b^2, a^1], X).
469
470mono_merge([], X, X):-!.
471mono_merge(X, [], X):-!.
472mono_merge([X^N|Xs], [X^N0|Xs0], [X^N1|Zs]):- !,
473	N1 is N+N0, mono_merge(Xs, Xs0, Zs).
474mono_merge([X^N|Xs],[X0^N0|Xs0], [X^N|Zs]):-  X @> X0, !,
475	mono_merge(Xs, [X0^N0|Xs0], Zs).
476mono_merge([X^N|Xs], [X0^N0|Xs0], [X0^N0|Zs]):-
477	mono_merge([X^N|Xs], Xs0, Zs).
478
479% ?- gblex:compare_mono(X, [y^1, x^2],[y^2, x^1]).
480%@ X = (<).
481
482compare_mono(C, [], []):-!, C=(=).
483compare_mono(C, [], _):-!, C=(<).
484compare_mono(C, _, []):-!,  C=(>).
485compare_mono(C0, [X^N|Xs], [Y^N0|Xs0]):- X==Y, !, compare(C, N, N0),
486	(  C == (=)  ->  compare_mono(C0, Xs, Xs0) ;  C0 = C  ).
487compare_mono(C, [X^_|_], [Y^_|_]):- X @< Y, !, C= (<).
488compare_mono((>), _, _).
489
490
491% ?- gblex:compare_poly(X, [1*[y^2, x^2]],[1*[y^1, x^2]]).
492%@ X = (>).
493compare_poly(C, [], []):-!, C=(=).
494compare_poly(C, [], _):-!,  C=(<).
495compare_poly(C, _, []):-!,  C=(>).
496compare_poly(C0, [_*M|Ms], [_*N|Ns]):-
497	compare_mono(C, M, N),
498	(	C == (=)
499	->	compare_poly(C0, Ms, Ns)
500	;       C0 = C
501	).
502
503% ?- gblex:sort_poly_list([[1*[a^1]], [1*[a^1]], [1*[b^1]]], X).
504%@ X = [[1*[b^1]],[1*[a^1]],[1*[a^1]]].
505sort_poly_list(X, Y):-
506	predsort( pred(( [A, [_*M|_], [_*M0|_]]:-
507			gblex:compare_mono(A0,  M, M0),
508			(	A0 == (>)
509			->	A  = (<)
510			;	A  = (>))
511		  )),
512		  X, Y).
513
514% ?- gblex:sort_poly([1*[a^1], 1*[b^1], 1*[a^1]], X).
515%@ X = [1*[b^1],1*[a^1],1*[a^1]].
516sort_poly(X, Y):-
517	predsort(pred( ([A, _*M, _*M0]:-
518			gblex:compare_mono(A0,  M, M0),
519			(	A0 == (>)
520			->	A  = (<)
521			;	A  = (>))
522		  )),
523		  X, Y).
524
525merge_poly([], []).
526merge_poly([0*_|Xs], Xs0):- !, merge_poly(Xs, Xs0).
527merge_poly([X], [X]).
528merge_poly([C*M,C0*M|Xs], Xs0):- !, C1 is C+C0,
529	merge_poly([C1*M|Xs], Xs0).
530merge_poly([X|Xs], [X|Xs0]):- merge_poly(Xs, Xs0).
531
532% ?- gblex:mul_mono_mono(1*[c^1, a^1], 2*[c^1, b^2], X).
533%@ X = 2*[c^2,b^2,a^1].
534mul_mono_mono(C*M, C0*M0, C1*M1):- C1 is C0*C, mono_merge(M, M0, M1).
535
536mul_scalar_poly(C, P, Q):- maplist(mul_scalar_mono(C), P, Q).
537
538mul_scalar_mono(C, C0*M, C1*M):- C1 is C0*C.
539
540% ?- gblex:div_mono_mono(-1*[b^2, c^1, a^3], 1*[b^1, a^1], D).
541div_mono_mono(C*M, C0*M0, C1*M1):- C1 is C rdiv C0,
542	div_mono_mono_v(M, M0, M1).
543
544div_mono_mono_v(M, [], M):- !.
545div_mono_mono_v([X^I|Ys], [X^J|Ys0], Zs):- !,
546	( I == J -> div_mono_mono_v(Ys, Ys0, Zs)
547	; I> J, K is I-J, Zs=[X^K|Zs0], div_mono_mono_v(Ys, Ys0, Zs0)
548	).
549div_mono_mono_v([X^I|Ys], [X0^J|Ys0], [X^I|Zs]):- compare((>), X, X0),
550	div_mono_mono_v(Ys, [X0^J|Ys0], Zs).
551
552% div_poly_mono_v(A, B, Q, R)
553
554% % ?- gblex:div_poly_mono([1*[x^2], 2*[x^1], 1*[]], [2*[x^1], 2*[]], X, []).
555
556% div_poly_mono(A, C*B, Q, R):- div_poly_mono_v(A, B, Q0, R),
557% 	C0 is 1 rdiv C,
558% 	mul_scalar_poly(C0, Q0, Q).
559
560% %
561% div_poly_mono_v([], _, [], []).
562% div_poly_mono_v([C*M|A], N, [C*N0|Q], R):- div_mono_mono_v(M, N, N0, []), !,
563% 	div_poly_mono_v(A, N, Q, R).
564% div_poly_mono_v([W|A], N, Q, [W|R]):-
565% 	div_poly_mono_v(A, N, Q, R).
566
567% ?- gblex:div_poly_poly([1*[x^2], 2*[x^1], 1*[]], [2*[x^1], 2*[]], X, []).
568%@ X = [1 rdiv 2*[x^1],1 rdiv 2*[]].
569div_poly_poly([], _, [], []):-!.
570div_poly_poly([M|A], [N|B], [N0|Q], R):- div_mono_mono(M, N, N0), !,
571	mul_mono_poly(N0, B, B0),
572	subtract_poly_poly(A, B0, P),
573	div_poly_poly(P, [N|B], Q, R).
574div_poly_poly(P, _, [], P).
575
576% ?- gblex:div_poly_poly([1*[x^2], 2*[x^1], 1*[]], [2*[x^1], 1*[]], X).
577div_poly_poly(P, Q, R):- div_poly_poly(P, Q, R, []).
578
579% ?- gblex:div_poly_poly([1*[x^2], 2*[x^1], 1*[]], [1*[x^1], 1*[]]).
580%@ true.
581div_poly_poly(P, Q):- div_poly_poly(P, Q, _, []).
582
583% ?- gblex:add_poly_poly([1*[b^1], 2*[a^1]], [1*[a^1]], X).
584%@ X = [1*[b^1],3*[a^1]].
588	(   B==(=)
589	->  C1 is C+C0,
590		(	C1 == 0
591		->	Zs=Zs0
592		;	Zs=[C1*M|Zs0]
593		),
595
596	;   B==(>)
597	->  Zs=[C*M|Zs0],
599
600	;   Zs=[C0*M0|Zs0],
602	).
603
604% ?- gblex:mul_mono_poly(1*[c^1], [1*[b^1], 1*[a^1]], X).
605%@ X = [1*[c^1,b^1],1*[c^1,a^1]].
606mul_mono_poly(M, P, P0):- maplist(mul_mono_mono(M),P, P0).
607
608% ?-numlist(1, 10000, L), maplist([I, I*[x^I]], L, P), gblex:mul_scalar_poly(-1, P, Q), time(100, gblex:subtract_poly_poly(P, Q, S), T).
609%@ T = 1.155068.
610
611subtract_poly_poly([], X,  Y):-	!, mul_scalar_poly(-1, X, Y).
612subtract_poly_poly(X, [],  X):-	!.
613subtract_poly_poly([C*M|Xs],[C0*M0|Xs0], Zs):-  compare_mono(B, M, M0),
614	(   B==(=)
615	->  C1 is C-C0,
616		(	C1 == 0
617		->	Zs=Zs0
618		;	Zs=[C1*M|Zs0]
619		),
620	    subtract_poly_poly(Xs, Xs0, Zs0)
621
622	;   B==(>)
623	->  Zs=[C*M|Zs0],
624	    subtract_poly_poly(Xs, [C0*M0|Xs0], Zs0)
625
626	;   C1 is  -C0,
627	    Zs=[C1*M0|Zs0],
628	    subtract_poly_poly([C*M|Xs], Xs0, Zs0)
629	).
630
631% ?- gblex: (mul_poly_poly([1*[b^1], 1*[a^1]],[1*[b^1], 1*[a^1]], X), mul_poly_poly(X, X, Y)).
632mul_poly_poly([], _, []):- !.
633mul_poly_poly(_, [], []):- !.
634mul_poly_poly([1*[]], P, P):- !.
635mul_poly_poly(P, [1*[]], P):- !.
636mul_poly_poly(P, Q, R):- mul_poly_poly(P, Q, [], R).
637
638mul_poly_poly([], _, P, P):- !.
639mul_poly_poly([M|Ms], P, Q, R):- mul_mono_poly(M, P, P0),
640	add_poly_poly(P0, Q, Q0), mul_poly_poly(Ms, P, Q0, R).
641
642% ?- gblex:power_poly(3, [2*[a^2]], U).
643power_poly(0, _, U):- !, unit_poly(U).
644power_poly(N, X, U):- N0 is N>>1, M is N mod 2,	power_poly(N0, X, U0),
645	mul_poly_poly(U0, U0, U1),
646	(M==1 -> mul_poly_poly(X, U1, U); U=U1).
647
648% ?-gblex:mono_lcm([b^2, a^1], [b^1, a^2], LCM, Y, Z).
649mono_lcm([], Y, Y, Y, []):-!.
650mono_lcm(X, [], X, [], X):-!.
651mono_lcm([V^I|X],[V^J|Y], [V^M|Z], P, Q):- !, compare(C, I, J),
652	( C == (=) -> M=I, mono_lcm(X, Y, Z, P, Q)
653	; C == (>) -> M=I, K is I-J, Q=[V^K|Q0], mono_lcm(X, Y, Z, P, Q0)
654	; M =J, K is J-I, P=[V^K|P0], mono_lcm(X, Y, Z, P0, Q)
655	).
656mono_lcm([U^I|X],[V^J|Y], [U^I|Z], P, [U^I|Q]):- V @< U,!,
657	mono_lcm(X, [V^J|Y], Z, P, Q).
658mono_lcm(X,[B|Y], [B|Z], [B|P], Q):- mono_lcm(X, Y, Z, P, Q).
659
660
661% basic laws  for simplifying normal_exp_left internal expressions
662% ?- term_rewrite(gblex:normal_exp_left, [], T0, T).
663% ?- term_rewrite(gblex:normal_exp_left, [], (1*x+0*x)* (1*x), T).
664%@ T = x*x.
665% ?- term_rewrite(gblex:normal_exp_left, [], (1*x - 0*x)* (-1*x), T).
666%@ T = -1*x*x.
667
668% ?- poly_term([1*[y^1], -1*[x^1]], X).
669%@ X = y-x.
670
671% Stage 1
672normal_exp_left([],	0).
673normal_exp_left(0+X,	X).
674normal_exp_left(+X,	X).
675normal_exp_left(1*X,	X).
676normal_exp_left(0*_,	0).
677normal_exp_left(X^1,	X).
678normal_exp_left(_^0,	1).
679normal_exp_left(A*B,	C):- rational(A), rational(B), C is A*B.
680normal_exp_left(A+B,	C):- rational(A), rational(B), C is A+B.
681normal_exp_left(A^B,	C):- rational(A), rational(B), C is A^B.
682normal_exp_left(-A,	C) :- rational(A), C is -A.
683normal_exp_left(-A,	(-1)*A).
684normal_exp_left(A*(B*C),	(A*B)*C).
685normal_exp_left(A+(B+C),	(A+B)+C).
686normal_exp_left(A*B,	B*A)	:- rational(B).
687normal_exp_left(A+B,	B+A)	:- rational(B).
688normal_exp_left(X*(Y+Z),	X*Y + X*Z).
689normal_exp_left((X+Y)*Z,	X*Z + Y*Z).
690normal_exp_left(X - Y,		X+ (-1)*Y).
691normal_exp_left(rdiv(A,B),	A/B).
692
693% Stage 2
694intro_prefix_minus_right(1*A, A).
695intro_prefix_minus_right(A+B, A+C):- extend_prefix_minus_scope(B, C).
696%
697extend_prefix_minus_scope(A,  -(B)):- rational(A), A<0, B is -A.
698extend_prefix_minus_scope(A*B, C*B):- extend_prefix_minus_scope(A, C).
699extend_prefix_minus_scope(-(A)*B, -(A*B)).
700
701% Stage 3 (final)
702intro_infix_minus(A + -(B), A - B).
703
704% ?- gblex:plain_geo_test(X).
705plain_geo_test(X):- F =
706(a*a1-1; b*b1-1; c*c1-1; d*d1-1; e*e1-1;
707 f*f1-1; g*g1-1; h*h1-1; i*i1-1; j*j1-1;
708 k*k1-1; l*l1-1; m*m1-1;
709e*(a+b)-f*(c+d);
710b*l-m*(c+d);
711b*(f+e)-a*(c+d);
712i*e-k*(c+d);
713g*d-h*(f+e);
714h*(a+b)-g*c;
715d*(a+b)-c*(f+e);
716k*(a+b)-i*f;
717m*(f+e)-l*a;
718m+l-(i+k);
719(a+b)-(f+e)-x),
720gb([    x,
721	a,  b,  c,  d,  e,  f,  g,  h,  j,  m,  l,  i,  k,
722   	a1, b1, c1, d1, e1, f1, g1, h1, j1, m1, l1, i1, k1
723       ], F, X),
724writeln(X)```