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.
  127% policy: keep  head unique.
  128% clearly correct based on Dickson's theorem on mono ideal.
  129% The Dickson theorem implies that no infinite sequence of
  130% incomparables head terms exists.
  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),
  140	reduce_head(C, G, D0),
  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),
  183	add_poly_poly(S0, R, U0),
  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),
  194	reduce_head(P0, G, P).
  195reduce_head(P, _, P).
  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),
  219	add_poly_poly(C, P0, R).
  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),
  270	add_affine_affine(P, Q0, R).
  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),
  274	add_affine_affine(P0, Q0, R).
  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),
  279	add_affine_affine(P0, Q0, R).
  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), !,
  301	add_poly_poly(A, B, C),
  302	add_affine_affine(P0, Q0, R).
  303add_affine_affine(P, Q, [X*A|R]):- select(X*A, P, P0), !,
  304	add_affine_affine(P0, Q, R).
  305add_affine_affine(P, Q, [X*A|R]):- select(X*A, Q, Q0), !,
  306	add_affine_affine(P, Q0, R).
  307add_affine_affine([], P, P):- !.
  308add_affine_affine(P, [], P):- !.
  309add_affine_affine([P], [Q], [R]):- add_poly_poly(P, Q, R).
  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) .
  431poly(X+Y,A1):-!,(poly(X,A2),poly(Y,A3)),add_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]].
  585add_poly_poly([], X, X):-!.
  586add_poly_poly(X, [], X):-!.
  587add_poly_poly([C*M|Xs],[C0*M0|Xs0], Zs):-  compare_mono(B, M, M0),
  588	(   B==(=)
  589	->  C1 is C+C0,
  590		(	C1 == 0
  591		->	Zs=Zs0
  592		;	Zs=[C1*M|Zs0]
  593		),
  594	    add_poly_poly(Xs, Xs0, Zs0)
  595
  596	;   B==(>)
  597	->  Zs=[C*M|Zs0],
  598	    add_poly_poly(Xs, [C0*M0|Xs0], Zs0)
  599
  600	;   Zs=[C0*M0|Zs0],
  601	    add_poly_poly([C*M|Xs], Xs0, 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)