1:- module(gblex, [gb/2, gb/3, gb/4, rat/2, term_poly/2, poly_term/2]). 2:- use_module(util(misc)). 4term_expansion --> pac:expand_pac.
5:- use_module(pac(op)). 6
7
19
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
66gb_rev(X, Y):- math:indeterminates(X, As),
67 reverse(As, Bs),
68 groebner_base(As, Bs, X, Y).
69
71groebner_base --> gb_in, !, buchberger, !, gb_out.
73groebner_base(Ord) --> {sort(Ord, Ord0)}, groebner_base(Ord, Ord0).
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).
82gb_in --> flatten(;),
83 maplist(term_poly),
84 maplist(canonical),
85 remove([]).
87gb_out --> predsort(compare_poly),
88 maplist(poly_term),
89 join_with_semicolon.
90
97
103
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
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
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
132buchberger(X, Y):- pairs(X, D),
133 buchberger(D, X, Y0),
134 normal_base(Y0, Y1),
135 maplist(canonical, Y1, Y).
136
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
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
151
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
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
169normal_poly(P, G, Q):- reduce_poly_one(P, G, P0), !,
170 normal_poly(P0, G, Q).
171normal_poly(P, _, P).
172
174normal_poly(P, G):- reduce_poly_one(P, G, _), !, fail.
175normal_poly(_, _).
176
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
188
189reduce_head([C*M|X], G, P):- member([_C*N|Y], G), 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
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
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
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
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
248
(_, [], [], []):-!.
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
256linear(Vs) --> gb_in, !,
257 gauss_sweep(Vs),
258 maplist(affine_to_poly),
259 gb_out.
260
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
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
283mul_scalar_affine(C, F, G):- mul_poly_affine([C*[]], F, G).
284
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
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
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
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
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
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
352
353join_term(F, [X, Y|Z], U):- !, join_term(F, [Y|Z], U0),
354 U =.. [F,X,U0].
355join_term(_, [X], X).
356
358zero_poly([]).
359unit_poly([1*[]]).
360constant_poly([_*[]]).
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
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
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
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
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
411term_term --> term_rat, rat_term.
412
416
417
418basic_poly(0,[]):-! .
419basic_poly(X,[X*[]]):-rational(X),! .
420basic_poly(X,[1*[X^1]]):-! .
421
422
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
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
465
467
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
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
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
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
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
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
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
553
555
559
566
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
577div_poly_poly(P, Q, R):- div_poly_poly(P, Q, R, []).
578
581div_poly_poly(P, Q):- div_poly_poly(P, Q, _, []).
582
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
606mul_mono_poly(M, P, P0):- maplist(mul_mono_mono(M),P, P0).
607
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
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
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
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
667
670
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
694intro_prefix_minus_right(1*A, A).
695intro_prefix_minus_right(A+B, A+C):- extend_prefix_minus_scope(B, C).
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
702intro_infix_minus(A + -(B), A - B).
703
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)