1:- module(gvec, []).    2:- use_module(gb(vector)).    3term_expansion --> pac:expand_pac.
    4:- use_module(pac(op)).    5
    6% %
    7% init_base(G, true) :- !, new_vector(G, 6).
    8% init_base([], _).
    9
   10%
   11add_to_base(Obj, G0, Index, G1, true):- new_index(Index, Obj, G0, G1).
   12add_to_base(R, G, R, [R|G], _).
   13
   14%
   15spoly(I, J, Spoly, G, Ord, FF, Vec):-
   16	if(FF,	s_poly_z(I, J, Spoly, G, Ord, Vec),
   17		s_poly(I, J, Spoly, G, Ord, Vec)).
   18
   19%
   20gb_BDFM([P|F], G, Opt):-
   21	subset([
   22			order(Ord),
   23			trace(Trace),
   24			assoc(Assoc),
   25			fraction_free(FF),
   26    		prime(Prime),
   27			vector(Vec)
   28		   ], Opt),
   29	if(Vec, SizePred = (vector:v_size), length),
   30	Debug =debug(Trace, Assoc, SizePred),
   31	memberchk(initial_vector_size(Size), Opt),
   32	new_vector(G0, Size),
   33%	init_base(G0, Vec),
   34	add_to_base(P, G0, _, G1, Vec),
   35	poly_pairs(F, [], D, G1, G2, Ord, Vec),
   36	gb_spoly(D, G2, G, Debug, Ord, FF, Prime, Vec).
   37
   38%
   39poly_pairs([], D, D, G, G, _, _).
   40poly_pairs([P|F], D0, D, G0, G, Ord, Vec):-
   41	add_to_base(P, G0, M, G1, Vec),
   42	update_spoly_agenda(P, M, D0, D1, G1, Ord, Vec),
   43	poly_pairs(F, D1, D, G1, G, Ord, Vec).
   44
   45%
   46gb_spoly([], G, G, _, _Ord, _, _, _).
   47gb_spoly([I-J|D], G0, G, Debug, Ord, FF, BL, Vec):-
   48	spoly(I, J, Spoly, G0, Ord, FF, Vec),
   49	once(reduce_head_by_polyset(Spoly, G0, S0, Ord, FF, _Trace, [])),
   50	once(gb:normal_poly(FF, BL, S0, R)),
   51%	gb:log(Debug, D, G0, R),
   52	(  R == []
   53	->	gb_spoly(D, G0, G, Debug, Ord, FF, BL, Vec)
   54	;	add_to_base(R, G0, M, G1, Vec),
   55		update_spoly_agenda(R, M, D, D0, G1, Ord, Vec),
   56		gb_spoly(D0, G1, G, Debug, Ord, FF, BL, Vec)
   57	).
   58
   59s_poly(I, J, Spoly, G, Ord, true):- !, v_ref(I, G, P),
   60	v_ref(J, G, Q),
   61	poly:s_poly(P, Q, Spoly, Ord).
   62s_poly(I, J, Spoly, _, Ord, false):-
   63	poly:s_poly(I, J, Spoly, _,  Ord).
   64
   65%
   66s_poly_z(I, J, Spoly, G, Ord, true):- !, v_ref(I, G, P),
   67	v_ref(J, G, Q),
   68	poly:s_poly_z(P, Q, Spoly, Ord).
   69s_poly_z(I, J, Spoly, _, Ord, _):-
   70	poly:s_poly_z(I, J, Spoly, Ord).
   71
   72%
   73head_term(G, I, X):-  v_ref(I, G, [_*X|_]).
   74
   75%
   76head_term(G, I, X, true):-  head_term(G, I, X).
   77head_term(_, I, I, _).
   78%
   79
   80head_term_lcm(G, I, J, X):- head_term_lcm(G, I, J, X, true).
   81
   82head_term_lcm(G, I, J, X, true):- !, head_term(G, I, Y),
   83	head_term(G, J, Z),
   84	poly:mono_lcm(Y, Z, X).
   85head_term_lcm(_, I, J, X, _):-  poly:mono_lcm(I, J, X).
   86
   87% head_term_lcm(G, I, J, X, true):- !, head_term(G, I, Y),
   88% 	head_term(G, J, Z),
   89% 	poly:mono_lcm(Y, Z, X).
   90% head_term_lcm(_, I, J, X, _):-  poly:mono_lcm(I, J, X).
   91
   92%%
   93update_spoly_agenda(R, M, D, D0, G, Ord, Vec):-
   94	if( (Vec, R=[_*H|_]),
   95		update_spoly_agenda(D, M-H, D0, G, Ord, Vec),
   96		update_spoly_agenda(D, M, D0, G, Ord, Vec)).
   97
   98%
   99update_spoly_agenda(D, M, D0, G, Ord, Vec):-
  100	update_spoly_agenda_B(D, U, [], M, G, Ord, Vec),
  101	update_spoly_agenda_DFM(M, V, G, Ord, Vec),
  102	sort_spoly(V, V0, G),
  103	merge_spoly_agenda(U, V0, D0, G, Vec).
  104
  105%
  106update_spoly_agenda_B([], D, D, _, _, _, _).
  107update_spoly_agenda_B([I-J|R], P, Q, M, G, Ord, Vec):-
  108	b_cond(M, I, J, G, Vec), !,
  109	update_spoly_agenda_B(R, P, Q, M, G, Ord, Vec).
  110update_spoly_agenda_B([A|R], [A|P], Q, M, G, Ord, Vec):-
  111	update_spoly_agenda_B(R, P, Q, M, G, Ord, Vec).
  112
  113%
  114update_spoly_agenda_DFM(M, D, G, Ord, true):-  head_term(G, 1, T1),
  115	update_spoly_agenda_DFM_vec(1-T1, M, D, [], G, Ord).
  116update_spoly_agenda_DFM(M, V, G, Ord, false):-
  117	update_spoly_agenda_DFM_list(G, G, M, V, [], Ord).
  118
  119%
  120update_spoly_agenda_DFM_vec(I-_, J-_, P, P, _, _):- I >= J, !.
  121update_spoly_agenda_DFM_vec(V, W, P, Q, G, Ord):-
  122	V = I-_,
  123	( d_cond(V, W, true);
  124	  f_cond(V, W, G, true);
  125	  m_cond(V, W, G, true) ),
  126	!,
  127	K is I + 1,
  128	head_term(G, K, Tk),
  129	update_spoly_agenda_DFM_vec(K-Tk, W, P, Q, G, Ord).
  130update_spoly_agenda_DFM_vec(I-_, W, [I-J|P], Q, G, Ord):- K is I+1,
  131	W = J-_,
  132	head_term(G, K, Tk),
  133	update_spoly_agenda_DFM_vec(K-Tk, W, P, Q, G, Ord).
  134
  135%
  136update_spoly_agenda_DFM_list([I|Gi], G, J, P, Q, Ord):-
  137	( f_cond(I, J, Gi, false);
  138	  m_cond(I, J, G, false);
  139	  d_cond(I, J, false) ),
  140	!,
  141	update_spoly_agenda_DFM_list(Gi, G, J, P, Q, Ord).
  142update_spoly_agenda_DFM_list([I|Gi], G, J, [I-J|P], Q, Ord):- !,
  143	update_spoly_agenda_DFM_list(Gi, G, J, P, Q, Ord).
  144update_spoly_agenda_DFM_list(_, _,  _, P, P, _).
  145
  146%
  147b_cond(K-Tk, I, J, G, true):- !, K > J,
  148	head_term(G, I, Ti),
  149	head_term(G, J, Tj),
  150	poly:mono_lcm(Ti, Tj, Tij),
  151	poly:div_mono_mono_term(Tij, Tk),
  152	poly:mono_lcm(Ti, Tk, Tik),
  153	Tik \== Tij,
  154	poly:mono_lcm(Tj, Tk, Tjk),
  155	Tjk \== Tij.
  156b_cond([_*Tk|_], [_*Ti|_], [_*Tj|_], _, false):-
  157	poly:mono_lcm(Ti, Tj, Tij),
  158	poly:div_mono_mono_term(Tij, Tk),
  159	poly:mono_lcm(Ti, Tk, Tik),
  160	Tik \== Tij,
  161	poly:mono_lcm(Tj, Tk, Tjk),
  162	Tjk \== Tij.
  163
  164%
  165d_cond(_-Ti, _-Tj, true):- !,
  166	poly: merge_mono_mono(Ti, Tj, T),
  167        poly: mono_lcm(Ti, Tj, T).
  168d_cond([_*Ti|_], [_*Tj|_], false) :-
  169	poly: merge_mono_mono(Ti, Tj, T),
  170        poly: mono_lcm(Ti, Tj, T).
  171
  172%
  173f_cond(I-Ti, _-Tj, G, true) :- !,
  174	poly:mono_lcm(Ti, Tj, Tij),
  175	I0 is I-1,
  176	between(1, I0, K),
  177	head_term(G, K, Tk),
  178	poly:mono_lcm(Tk, Tj, Tij).
  179f_cond([_*Ti|_], [_*Tj|_], G, false) :-
  180	poly:mono_lcm(Ti, Tj, T),
  181%	poly:rev_member([_*Tk|_], G),
  182	v_member([_*Tk|_], G),
  183	poly:mono_lcm(Tk, Tj, T).
  184
  185m_cond(I-Ti, J-Tj, G, true):- !, M is J-1,
  186	poly:mono_lcm(Ti, Tj, Tij),
  187	between(1, M, K),
  188	K\== I,
  189%	call(G, K, Tk),
  190	v_ref(K, G, Tk),
  191	poly:div_mono_mono_term(Tij, Tk),
  192	poly:mono_lcm(Tj, Tk, Tjk),
  193	Tjk \== Tij.
  194m_cond(Pi, Pj, G, false) :-   Pi = [_*Ti|_],
  195	Pj = [_*Tj|_],
  196	poly:mono_lcm(Ti, Tj, Tij),
  197%	poly:rev_member(Pk, G),
  198	v_member(Pk, G),
  199	Pk \== Pi,
  200	Pk = [_*Tk|_],
  201	poly:div_mono_mono_term(Tij, Tk),
  202	poly:mono_lcm(Tj, Tk, Tjk),
  203	Tjk \== Tij.
  204
  205%
  206merge_spoly_agenda([], X, X, _, _).
  207merge_spoly_agenda(X, [], X, _, _).
  208merge_spoly_agenda([P|R], [Q|S], [P|T], G, Vec):-
  209	if(Vec, compare_spoly(G, C, P, Q),
  210		compare_spoly(C, P, Q)),
  211	C = (<),
  212	!,
  213	merge_spoly_agenda(R, [Q|S], T, G, Vec).
  214merge_spoly_agenda(R, [Q|S], [Q|T], G, Vec):-
  215	merge_spoly_agenda(R, S, T, G, Vec).
  216
  217%
  218compare_spoly(G, C, I-J, I0-J0):- head_term_lcm(G, I, J, LCM),
  219	head_term_lcm(G, I0, J0, LCM0),
  220	poly:compare_total_order(C0, LCM, LCM0),
  221	( C0 == (=) ->  C = (<);  C = C0).
  222
  223% ?- bdfm:compare_spoly(C, [1*[a^2]]-[1*[b^2]], [1*[b^1]]-[1*[a^2]]).
  224%@ C =  (>).
  225compare_spoly(C, [_*M|_]-[_*N|_], [_*M0|_]-[_*N0|_]):-
  226	poly:mono_lcm(M,  N,  A),
  227	poly:mono_lcm(M0, N0, A0),
  228	poly:compare_total_order(C0, A, A0),
  229	(C0 == (=) -> C = (<); C = C0 ).
  230%
  231sort_spoly(X, Y, G):- predsort(compare_spoly(G), X, Y).
  232%
  233sort_spoly(X, Y):- predsort(compare_spoly, X, Y).
  234
  235%
  236reduce_head_by_polyset(X, Ps, R, Ord, FF, T0, T):-
  237	reduce_head_by_polyset_one(X, Ps, Z, Ord, FF, T0, T1), !,
  238	reduce_head_by_polyset(Z, Ps, R, Ord, FF, T1, T).
  239reduce_head_by_polyset(X, _, X, _, _, T, T).
  240
  241reduce_head_by_polyset_one(X, Ps, R, Ord, true, [P|T], T):- !,
  242	v_member(P, Ps),
  243	once(reduce_head_by_poly_z(X, P, R, Ord)).
  244reduce_head_by_polyset_one(X, Ps, R, Ord, _, [P|T0], T):-
  245%	poly:rev_member(P, Ps),
  246	v_member(P, Ps),
  247	once(reduce_head_by_poly(X, P, R, Ord, T0, T)).
  248
  249reduce_by_polyset(X, Ps, R, Ord, FF):-
  250	reduce_by_polyset_one(X, Ps, Z, Ord, FF, _Trace, []), !,
  251	reduce_by_polyset(Z, Ps, R, Ord, FF).
  252reduce_by_polyset(X, _, X, _, _).
  253
  254%
  255reduce_by_polyset_one(X, Ps, R, Ord, true, [P|T], T):- !,
  256%	poly:rev_member(P, Ps),
  257	v_member(P, Ps),
  258	once(reduce_by_poly_z(X, P, R, Ord)).
  259reduce_by_polyset_one(X, Ps, R, Ord, _, [P|T], T):-
  260%	poly:rev_member(P, Ps),
  261	v_member(P, Ps),
  262	once(reduce_by_poly(X, P, R, Ord)).
  263
  264
  265% reduce_by_polyset_one(X, Ps, R, Ord, FF, BL):- select(M, X, X0),
  266% 	poly:rev_member([N|Y], Ps),
  267% %	member([N|Y], Ps),
  268% 	poly:div_mono_mono(M, N, D),
  269% 	!,
  270% 	poly:mul_mono_poly(D, Y, Q),
  271% 	poly:subtract_poly_poly(X0, Q, R, Ord).
  272
  273
  274% ?- gb:reduce_by_poly([1*[x^1], 1*[a^1]], [1*[x^1], -1*[]], R, lexical).
  275%@ R = [1*[a^1], 1*[]] .
  276
  277reduce_head_by_poly([M|X], [N|Y], R, Ord):-
  278	poly:div_mono_mono(M, N, D),
  279	!,
  280	poly:mul_mono_poly(D, Y, Q),
  281	poly:subtract_poly_poly(X, Q, R, Ord).
  282
  283% ?- gb:reduce_by_poly_z([1*[x^1], 1*[a^1]], [1*[x^1], -1*[]], R, lexical).
  284%@ R = [1*[a^1], 1*[]] .
  285% ?- gb:reduce_by_poly_z([4*[x^1], 1*[a^1]], [6*[x^1], -1*[]], R, lexical).
  286%@ R = [3*[a^1], 2*[]] .
  287
  288reduce_head_by_poly_z([A*M|X], [B*N|Y], R, Ord):-
  289	poly:div_mono_mono_term(M, N, N0),
  290	!,
  291	poly:lcm(A, B, L),
  292	A0  is L // A,
  293	B0  is L // B,
  294	poly:mul_scalar_poly(A0, X, Q),
  295	poly:mul_mono_poly(B0*N0, Y, Y0),
  296	poly:subtract_poly_poly(Q, Y0, R, Ord).
  297
  298% ?- gb:reduce_by_poly([1*[x^1], 1*[a^1]], [1*[x^1], -1*[]], R, lexical).
  299%@ R = [1*[a^1], 1*[]] .
  300% ?- gb:reduce_by_poly([1*[x^1], 1*[a^1]], [1*[x^2], -1*[]], R, lexical).
  301%@ false.
  302
  303reduce_by_poly(X, [N|Y], R, Ord):- select(M, X, X0),
  304	poly:div_mono_mono(M, N, D),
  305	!,
  306	poly:mul_mono_poly(D, Y, Q),
  307	poly:subtract_poly_poly(X0, Q, R, Ord).
  308
  309% ?- gb:reduce_by_poly_z([1*[x^1], 1*[a^1]], [1*[x^1], -1*[]], R, lexical).
  310%@ R = [1*[a^1], 1*[]] .
  311% ?- gb:reduce_by_poly_z([4*[x^1], 1*[a^1]], [6*[x^1], -1*[]], R, lexical).
  312%@ R = [3*[a^1], 2*[]] .
  313% ?- gb:reduce_by_poly_z([4*[x^1], 1*[a^1]], [6*[x^2], -1*[]], R, lexical).
  314%@ false.
  315
  316reduce_by_poly_z(X, [B*N|Y], R, Ord):- select(A*M, X, X0),
  317	poly:div_mono_mono_term(M, N, N0),
  318	!,
  319	poly:lcm(A, B, L),
  320	A0  is L // A,
  321	B0  is L // B,
  322	poly:mul_scalar_poly(A0, X0, Q),
  323	poly:mul_mono_poly(B0*N0, Y, Y0),
  324	poly:subtract_poly_poly(Q, Y0, R, Ord).
  325
  326normal_poly(true, Prime, Poly, R):- !,
  327	(	Prime == 0
  328	->	once(poly:remove_content_poly(Poly, R))
  329	;	once(poly:galois_poly(Prime, Poly, R))
  330	).
  331normal_poly(_, _, Poly, R):- once(poly:canonical_form(Poly, R))