1:- module(poly, [bits/3, poly_in/2, poly_out/2
2 ]). 3
5:- use_module(pac(basic)). 6:- use_module(util(meta2)). 7:- use_module(pac('expand-pac')). 9term_expansion --> pac:expand_pac.
10:- use_module(pac(op)). 11:- style_check(-singleton). 12
22
23gcd(X, Y, Z):- Z is gcd(X, Y).
24
29
30lcm(X, Y, Z):- Z is X*Y // gcd(X, Y).
31
36lcm(X, Y, Z, U, V):- Z is X*Y // gcd(X, Y), U is Z//X, V is Z//Y.
37
40
43
53
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
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).
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
88mono_tail(M, X, T, D):- mono_tail(M, X, T, 0, D).
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
96list_max(L, M):- list_max(L, 0, M).
97
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
106
107int_div(X, D, Q, R):- divmod(X, D, Q, R).
108
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
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
139
147
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
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
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
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
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
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
211
216
221
222poly_in --> flatten(;),
223 maplist(exp_poly),
224 remove([]).
225
227poly_in([], []) --> poly_in.
228
229poly_in(Order, Assoc) -->
230 pred([Order, Assoc], [X, Y]:- preprocess(X, Y, Order, Assoc)),
231 poly_in.
232
234poly_out --> poly_out([]).
235
236poly_out(Assoc) -->
237 maplist(pred(Assoc, [X, Y]:- postprocess(X, Y, Assoc))),
238 maplist(poly_exp).
239
241print_polyset(X):- maplist(writeln, X).
242
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
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
266canonical_form([]).
267canonical_form([1*_|_]).
268
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
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
288zero_poly([]).
289unit_poly([1*[]]).
290constant_poly([_*[]]).
291
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
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
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
323rat_rat(r(P, Q), r(P0, Q0)):- poly(P, P0, lexical),
324 poly(Q, Q0, lexical).
325
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
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
353
354mono_exp(C*M, E):- maplist(pred([X^1, X] & [X, X]), M, M0),
355 mono_exp_fold(C*M0, E).
356
358mono_exp_fold(C*[], C).
359mono_exp_fold(C*[A|R], C*T) :- foldl(pred([U, X, X*U]), R, A, T).
360
362slim_mono_exp(1*M, M).
363slim_mono_exp(M, M).
365rdiv_slash(rdiv(X,Y), X/Y).
366rdiv_slash(X, X).
367
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
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
382exp_exp --> exp_rat, rat_exp.
383
384
386
387
388basic_poly(0,[]):-! .
389basic_poly(X,[X*[]]):-rational(X),! .
390basic_poly(X,[1*[X^1]]):-! .
391
392
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
416
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
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
442make_var(X, 1*[X^1]).
443
444make_constant(I, I*[]).
445
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
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
456apply_poly(-, A, C, _):- prefix_minus_poly(A, C).
457apply_poly(pow(N), A, C, SW):- power_poly(N, A, C, SW).
458
471
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
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
518
519poly_calc(X, Y, Ord):- non_modular_poly(X, Y0, Ord), poly_exp(Y0, Y).
520
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
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
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
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
602
603preprocess(X, Y, Assoc):- once(preprocess(X, Y, [], Assoc)).
604
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
616assoc_numbering(A):- length(A, N),
617 numlist(1, N, L),
618 zip_hyphen(_, L, A).
619
624
625number_vars(X, Assoc, Y):-
626 indeterminates(X, Y, B, []),
627 subtract_assoc(B, Assoc, B0),
628 maplist(pred([X-X]), B0). 629
633subtract_assoc(X, Y, Z):-
634 foldl(pred(Y, ( [A, U, U] :- memberchk(A, Y))
635 & ( [A, [A|U], U]) ),
636 X, Z, []).
637
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
648merge_assoc([], []).
649merge_assoc([X-N|R], [X-N|A]):- merge_assoc(R, R0, X, N),
650 merge_assoc(R0, A).
652merge_assoc([X-N|R], R0, X, N):- !, merge_assoc(R, R0, X, N).
653merge_assoc(R, R, _, _).
654
656zip_hyphen([], [], []).
657zip_hyphen([A|B], [C|D], [A-C|R]):- zip_hyphen(B, D, R).
658
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
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
678
679sort_mono(L, R) :- length(L, N),
680 sort_mono(N, L, _, R).
681
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
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
700
701dim_mono_term([], 0).
702dim_mono_term([_^I|R], N):- dim_mono_term(R, N0), N is I+N0.
703
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
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
733compare_mono(lexical, C, _*M, _*N):- !, compare_mono_term(C, M, N).
734compare_mono(total, C, _*M, _*N):- !, compare_total_order(C, M, N).
738
740merge_mono([X^I|R], [X^J|S]):- merge_mono(R, R0, X, I, J),
741 merge_mono(R0, S).
742merge_mono([], []).
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).
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
758mul_mono_mono(C*M, D*N, E*L):- merge_mono_mono(M, N, L), E is C*D.
759
761mul_mono_poly(M, P, P0):- maplist(mul_mono_mono(M),P, P0).
762
772sort_polyset(F, X, Y):- predsort(compare_poly(F), X, Y).
773
775compare_poly(_, =, [], []).
776compare_poly(_, <, [], _).
777compare_poly(_, >, _, []).
778compare_poly(F, C, [M|_], [N|_]):- compare_mono(F, C, M, N),
779 C\== (=).
780compare_poly(_, <, _, _).
781
788
789sort_poly(F, L, R) :- length(L, N),
790 sort_poly(N, L, _, R, F).
791
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
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
820merge_polyset_polyset(X, Y, Z, F):- merge_polyset_polyset(X, Y, Z, [], F).
821
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
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).
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
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
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
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
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
890
892div_mono_mono_term(X, Y):- div_mono_mono_term(X, Y, _).
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
907div_mono_mono(C*M, C0*M0, C1*M1):- div_mono_mono_term(M, M0, M1),
908 C1 is C rdiv C0.
909
910
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
926
927add_poly_poly(X, Y, Z, F):- merge_poly_poly(X, Y, Z, F).
928
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
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
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 958
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
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
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
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
995add_coeff(C*M, D*_, E*M):- E is C+D, E\==0, !.
996add_coeff(_, _, []).
997
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
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
1023galois_poly(G, P, Q):- rem_poly(P, Q, G).
1024
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
1037
1038power_poly(N, X, Y, SW):- bits(N, [], Bs),
1039 unit_poly(U),
1040 power_poly(Bs, X, U, Y, SW).
1041
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
1058
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
1078mono_lcm_z(X, Y, U, V):- mono_lcm_z(X, Y, _, U, V).
1079
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
1085rev_member(X, [_|R]):- rev_member(X, R).
1086rev_member(X, [X|_]).
1087
1092
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
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
1112mod_exp(P, N, X, Y):- bits(N, [], Bs), !, mod_exp(P, Bs, X, 1, Y).
1113
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)