1:- module(gb, []).    2
    3:-use_module(pac('expand-pac')).    4:-use_module(gb('gb-vector')).    5:-use_module(gb('bdfm-generic')).    6
    7term_expansion --> pac:expand_pac.
    8:- use_module(pac(op)).    9
   10% ?- trace, gb:test(X).
   11%
   12default_option([
   13		order(lexical),
   14		prime(0),				% 2147483647,	17
   15		remove_content(24),		% timing bit length
   16		fraction_free(true),
   17		vars([]),
   18		homogeneous(false),
   19		modular(false),
   20		minimal(true),
   21		trace(false),
   22		vector(false),
   23		initial_vector_size(8),
   24		gc(false),
   25		trim(false),
   26		exp(true),
   27		verbose(false)
   28		]).
   29
   30%
   31generic_normal_poly(true, Prime, Poly, R):- !,
   32	(	Prime == 0
   33	->	once(poly:remove_content_poly(Poly, R))
   34	;	once(poly:galois_poly(Prime, Poly, R))
   35	).
   36generic_normal_poly(_, _, Poly, R):- once(poly:canonical_form(Poly, R)).
   37
   38% ?- completing_option([], Y).
   39completing_option(X,  Y):- default_option(Ds),
   40 		completing_option(Ds, X, Y).
   41
   42%
   43completing_option([], X, X).
   44completing_option([A| As], X,  Y):-
   45	(	check_option(A, X)
   46	->	X0 = X
   47	;	X0 = [A|X]
   48	),
   49	completing_option(As, X0, Y).
   50
   51check_option(A, X):- functor(A, F, 1),
   52	functor(B, F, 1),
   53	memberchk(B, X).
   54
   55% ?- gb:help.
   56help :-  default_control(X), writeln(X).
   57
   58% ?- gb:option_to_zip([a(1), b(2)], Z),
   59option_to_zip(X, Y):- maplist(pred([A, K-V]:- A =..[K,V]), X, Y).
   60
   61%
   62% ?- gb:gb((a;b), Y, []).
   63% ?- gb:gb(a, Y, []).
   64% ?- gb:gb((a+b; a-b), Y, [vars([a,b])]).
   65% ?- gb:gb((a+b; a-b), Y, [vars([b,a])]).
   66% ?- gb:gb((a^2+b^2-1; a-b), Y, []).
   67% ?- gb:gb((a^2+b^2-1; a-b), Y, [vars([b,a])]).
   68% ?- gb:gb((a+b)^100, Y, []).
   69% ?- gb:gb((a+b-c)^20, Y, []).
   70% ?- gb:gb(x^2+y^2-4; x*y-1,  Y, []).
   71% ?- gb:gb((a^2+b^2-1; a-b), Y, [homogeneous(true)]).
   72% ?- gb:gb((a^2+b^2-1; a-b), Y, [homogeneous(true), vector(true)]).
   73% ?- gb:gb(a, Y, [vector(true)]).
   74% ?- gb:gb(a;b, Y, [vector(true)]).
   75% ?- gb:gb(a;b;a, Y, [vector(true)]).
   76% ?- gb:gb(a-b;b+a, Y, [vector(true)]).
   77%@ Y = [b, a] .
   78
   79% ?- gb:gb((a+b)^100, Y, [vector(true)]).
   80% ?- gb:gb((a+b-c)^20, Y, [vector(true)]).
   81% ?- gb:gb((a^2+b^2-1; a-b), Y, [vector(true)]).
   82% ?- gb:gb((a^2+a-1; a), Y, [vector(true)]).
   83% ?- gb:gb((a^2+a-1; a-1), Y, [vector(true)]).
   84
   85% ?- gb:gb(r(3)^2+1; i^2+1; (r(3)+i +a)^3; r(3)^4+ i^3-ans, X, [vars([ans,r(3),i])] ).
   86% ?- gb:gb(r(3)^2+1; i^2+1; (r(3)+i +a)^3; r(3)^4+ i^3-ans, X, [vars([ans,r(3),i]), vector(true)] ).
   87
   88% ?- F1 = x^5+y^4+z^3-1, F2=x^3+y^3+z^2-1, time(gb:gb(F1;F2, X, [vars([z,y,x])])).
   89%@ % 123,369,311 inferences, 14.531 CPU in 14.596 seconds (100% CPU, 8490174 Lips)
   90% ?- F1 = x^5+y^4+z^3-1, F2=x^3+y^3+z^2-1, time(gb:gb(F1;F2, X, [vars([z,y,x]),vector(true)])).
   91%@ % 255,189,436 inferences, 44.758 CPU in 44.867 seconds (100% CPU, 5701514 Lips)
   92%@ F1 = x^5+y^4+z^3-1,
   93%@ F2 = x^3+y^3+z^2-1,
   94%@ X = [... + ... + ... * ... - 1617* ... ^ ... - 3228*z^7-1614*z^6+378*z^5+872*z^4+432*z^3+80*z^2, ... + ... + ... * ... - 27268337* ... ^ ... - 14641598*z^4+1214234*z^3+2619854*z^2+139968*z-69984, ... + ... - ... * ... - 26735354* ... ^ ... - 7526542*z^3+650492*z^2+34992*z-34992, ... - ... - ... * ... - 9673136* ... ^ ... - 95162*z^2+31104*z-31104, x^3+y^3+z^2-1, ... + ... + ... ^ ... + z^3-1, ... - ... + ... * ... - 1, ... + ... - ... * ...] .
   95
   96% ?- F1 = x^5+y^4+z^3-1, F2=x^3+y^3+z^2-1, time(gb:gb(F1;F2, X, [order(total), vars([z,y,x]),vector(true)])).
   97% ?- F1 = x^5+y^4+z^3-1, F2=x^3+y^3+z^2-1, time(gb:gb(F1;F2, X, [vars([x,y, z])])).
   98%@ % 406,849 inferences, 0.039 CPU in 0.040 seconds (98% CPU, 10311461 Lips)
   99% ?- F1 = x^5+y^4+z^3-1, F2=x^3+y^3+z^2-1, time(gb:gb(F1;F2, X, [vars([x,y, z]), vector(true)])).
  100%@ % 1,938,943 inferences, 0.225 CPU in 0.228 seconds (99% CPU, 8600858 Lips)
  101
  102
  103%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  104%  Do not try the following two, which seem to take     %
  105%  unexpectedly long time. bug ?                        %
  106%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  107
  108% ?- F1 = x^5+y^4+z^3-1, F2=x^3+y^3+z^2-1, gb:gb(F1;F2, X, [order(total), vars([x,y, z]), vector(true)]).
  109% ?- F1 = x^5+y^4+z^3-1, F2=x^3+y^3+z^2-1, gb:gb(F1;F2, X, [order(total), vars([x, y, z])]).
  110
  111gb(X, Y, Opt):- gb_setup(X, Opt, U, Opt0), !,
  112				gb_BDFM(U, V, Opt0), !,
  113				gb_clean(V, Y, Opt0).
  114
  115gb_setup(X, Opt, Y, [assoc(Assoc)|Opt0]):-
  116		completing_option(Opt, Opt0),
  117
  118		if( memberchk(verbose(true), Opt),
  119			(	writeln("Options:"),
  120				maplist(pred([X]:-format("\t~w\n", [X])), Opt)
  121			),
  122			true),
  123		memberchk(vars(Vars),  Opt0),
  124		poly:poly_in(Vars, Assoc, X, X0),
  125
  126		if( memberchk(homogeneous(true), Opt0),
  127			maplist(poly:make_poly_homogeneous(0), X0, X1),
  128			X1 = X0),
  129
  130		if( memberchk(fraction_free(true), Opt0),
  131			Y = X1,
  132			maplist(poly:canonical_form, X1, Y)).
  133
  134%
  135gb_clean(X, Y, Opt):-
  136	if( memberchk(vector(true), Opt),
  137	    vector:list_to_vector(X0, X),
  138	    X0 = X),
  139
  140	if( memberchk(gc(true),	Opt), garbage_collect),
  141
  142	if( memberchk(trim(true), Opt), trim_stacks),
  143
  144	if( memberchk(minimal(true), Opt),
  145	    minimal_base(X0, X1, Opt),
  146	    X1 = X0),
  147
  148	if( memberchk(homogeneous(true), Opt),
  149		(	memberchk(order(Ord), Opt),
  150			maplist(poly:make_poly_non_homogeneous(Ord), X1, X2)
  151		),
  152	    X2 = X1),
  153
  154	if( memberchk(exp(true), Opt),
  155		(	memberchk(assoc(Assoc), Opt),
  156			poly:poly_out(Assoc, X2, Y)
  157		),
  158		Y = X2).
  159
  160gb_BDFM(X, Y, Opt):-
  161	memberchk(modular(Mod), Opt),
  162	memberchk(vector(Vec),  Opt),
  163	if(Mod,
  164	   gbtrace:gb_BDFM(X,  Y,  Opt),    % <=  Not yet implemented.
  165	   if(Vec,
  166			gvec:gb_BDFM(X, Y, Opt),    % <= using vectors (experimental)
  167			bdfm:gb_BDFM(X, Y, Opt))).	% <= usual list programming
  168
  169%  only for debugging
  170compare_gb_methods(C, Opt1, Opt2, X):-
  171	gb([exp(false)|Opt1], X, Y1, CK1),
  172	gb([exp(false)|Opt2], X, Y2, CK2),
  173	memberchk(order-Ord1, CK1),
  174	memberchk(order-Ord2, CK2),
  175	maplist(pred(Ord1, ([A, B]:- poly:sort_poly(Ord1, A, A0),
  176				poly:canonical_form(A0, B))), Y2, Z2),
  177	maplist(pred(Ord2, ([A, B]:- poly:sort_poly(Ord2, A, A0),
  178				poly:canonical_form(A0, B))), Y1, Z1),
  179	maplist(pred([Y1,Ord1],
  180		     ([U, V]:- reduce_by_polyset(U, Y1, V, Ord1))), Z2, R2),
  181	maplist(pred([Y2,Ord2],
  182		     ([U, V]:- reduce_by_polyset(U, Y2, V, Ord2))), Z1, R1),
  183	compare_gb(C, R1, R2).
  184
  185%
  186compare_gb(C, R1, R2):- maplist(=([]), R1), !,
  187	(	maplist(=([]), R2)
  188	->	C = (=)
  189	;	C = (<)
  190	).
  191compare_gb(C, _, R2):-
  192	(	maplist(=([]), R2)
  193	->	C = (>)
  194	;	C = (incomparable)
  195	).
  196
  197% ?- gb: ideal_membership(a, a, []).
  198% ?- gb: ideal_membership(a, a;b, []).
  199
  200ideal_membership(P, Q, Opt0):-
  201	gb([exp(false)|Opt0], Q, GB, Opt),
  202	memberchk(order(Ord), Opt),
  203	memberchk(assoc(Assoc), Opt),
  204	poly:number_vars(P, Assoc, P0),
  205	poly:poly_in(P0, Ps),
  206	maplist(poly:sort_poly(Ord), Ps, Ps0),
  207	maplist(pred([GB,Ord],
  208		     ([U, R]:- once(reduce_by_polyset(U, GB, R, Ord)))),
  209		Ps0, Rs),
  210	!,
  211	maplist(=([]), Rs).
  212
  213% ?- gb:minimal_base([[1*[x^2], 1*[x^1], 1*[a^1]], [1*[x^1]], [1*[x^2]]], R,
  214
  215minimal_base(X, Y, Opt):-
  216	subset([order(Ord), fraction_free(FF), prime(BL)], Opt),
  217	reduce_poly_head(X, Y0, Ord, FF, BL, _Trace, []), % <== indispensable  ([2014/07/20])
  218	normalize_base(Y0, Y, Ord, FF, BL).
  219%
  220reduce_poly_head(X, Y, Ord, FF, BL, [P|T0], T):-
  221	select(P, X, Y0),
  222	reduce_head_by_polyset_one(P, Y0, Q, Ord, FF, T0, T1),
  223	!,			% <== indispensable ([2014/07/21])
  224	(   poly:zero_poly(Q)
  225	->	reduce_poly_head(Y0, Y, Ord, FF, BL, T1, T)
  226	;	once(generic_normal_poly(FF, BL, Q, Q0)),
  227		reduce_poly_head([Q0|Y0], Y, Ord, FF, BL, T1, T)
  228	).
  229reduce_poly_head(X, X, _, _, _, T, T).
  230
  231%
  232normalize_base(X, Y, Ord, FF, BL):- select(P, X, Y0),
  233	reduce_by_polyset_one(P, Y0, Q, Ord, FF, _Trace, []),
  234	!,
  235	(   poly:zero_poly(Q)
  236	->	normalize_base(Y0, Y, Ord, FF, BL)
  237	;	once(generic_normal_poly(FF, BL, Q, Q0)),
  238		normalize_base([Q0|Y0], Y, Ord, FF, BL)
  239	).
  240normalize_base(X, X, _, _,_).
  241
  242%
  243reduce_head_by_polyset(X, Ps, R, Ord, FF, T0, T):-
  244	reduce_head_by_polyset_one(X, Ps, Z, Ord, FF, T0, T1), !,
  245	reduce_head_by_polyset(Z, Ps, R, Ord, FF, T1, T).
  246reduce_head_by_polyset(X, _, X, _, _, T, T).
  247
  248%
  249% reduce_head_by_polyset_one([M|X0], Ps, R, Ord, FF):-
  250% 	poly:rev_member([N|Y], Ps),
  251% %	member([N|Y], Ps),
  252% 	poly:div_mono_mono(M, N, D),
  253% 	!,
  254% 	poly:mul_mono_poly(D, Y, Q),
  255% 	poly:subtract_poly_poly(X0, Q, R, Ord).
  256
  257reduce_head_by_polyset_one(X, Ps, R, Ord, true, [P|T], T):- !,
  258	poly:rev_member(P, Ps),
  259%	member(P, Ps),
  260	once(reduce_head_by_poly_z(X, P, R, Ord)).
  261reduce_head_by_polyset_one(X, Ps, R, Ord, _, [P|T0], T):-
  262	poly:rev_member(P, Ps),
  263%	member(P, Ps),
  264	once(reduce_head_by_poly(X, P, R, Ord, T0, T)).
  265
  266reduce_by_polyset(X, Ps, R, Ord, FF):-
  267	reduce_by_polyset_one(X, Ps, Z, Ord, FF, _Trace, []), !,
  268	reduce_by_polyset(Z, Ps, R, Ord, FF).
  269reduce_by_polyset(X, _, X, _, _).
  270
  271%
  272reduce_by_polyset_one(X, Ps, R, Ord, true, [P|T], T):- !,
  273	poly:rev_member(P, Ps),
  274%	member(P, Ps),
  275	once(reduce_by_poly_z(X, P, R, Ord)).
  276reduce_by_polyset_one(X, Ps, R, Ord, _, [P|T], T):-
  277	poly:rev_member(P, Ps),
  278%	member(P, Ps),
  279	once(reduce_by_poly(X, P, R, Ord)).
  280
  281% reduce_by_polyset_one(X, Ps, R, Ord, FF, BL):- select(M, X, X0),
  282% 	poly:rev_member([N|Y], Ps),
  283% %	member([N|Y], Ps),
  284% 	poly:div_mono_mono(M, N, D),
  285% 	!,
  286% 	poly:mul_mono_poly(D, Y, Q),
  287% 	poly:subtract_poly_poly(X0, Q, R, Ord).
  288
  289% ?- gb:reduce_by_poly([1*[x^1], 1*[a^1]], [1*[x^1], -1*[]], R, lexical).
  290%@ R = [1*[a^1], 1*[]] .
  291
  292reduce_head_by_poly([M|X], [N|Y], R, Ord):-
  293	poly:div_mono_mono(M, N, D),
  294	!,
  295	poly:mul_mono_poly(D, Y, Q),
  296	poly:subtract_poly_poly(X, Q, R, Ord).
  297
  298% ?- gb:reduce_by_poly_z([1*[x^1], 1*[a^1]], [1*[x^1], -1*[]], R, lexical).
  299%@ R = [1*[a^1], 1*[]] .
  300% ?- gb:reduce_by_poly_z([4*[x^1], 1*[a^1]], [6*[x^1], -1*[]], R, lexical).
  301%@ R = [3*[a^1], 2*[]] .
  302
  303reduce_head_by_poly_z([A*M|X], [B*N|Y], R, Ord):-
  304	poly:div_mono_mono_term(M, N, N0),
  305	!,
  306	poly:lcm(A, B, L),
  307	A0  is L // A,
  308	B0  is L // B,
  309	poly:mul_scalar_poly(A0, X, Q),
  310	poly:mul_mono_poly(B0*N0, Y, Y0),
  311	poly:subtract_poly_poly(Q, Y0, R, Ord).
  312
  313% ?- gb:reduce_by_poly([1*[x^1], 1*[a^1]], [1*[x^1], -1*[]], R, lexical).
  314%@ R = [1*[a^1], 1*[]] .
  315% ?- gb:reduce_by_poly([1*[x^1], 1*[a^1]], [1*[x^2], -1*[]], R, lexical).
  316%@ false.
  317
  318reduce_by_poly(X, [N|Y], R, Ord):- select(M, X, X0),
  319	poly:div_mono_mono(M, N, D),
  320	!,
  321	poly:mul_mono_poly(D, Y, Q),
  322	poly:subtract_poly_poly(X0, Q, R, Ord).
  323
  324% ?- gb:reduce_by_poly_z([1*[x^1], 1*[a^1]], [1*[x^1], -1*[]], R, lexical).
  325%@ R = [1*[a^1], 1*[]] .
  326% ?- gb:reduce_by_poly_z([4*[x^1], 1*[a^1]], [6*[x^1], -1*[]], R, lexical).
  327%@ R = [3*[a^1], 2*[]] .
  328% ?- gb:reduce_by_poly_z([4*[x^1], 1*[a^1]], [6*[x^2], -1*[]], R, lexical).
  329%@ false.
  330
  331reduce_by_poly_z(X, [B*N|Y], R, Ord):- select(A*M, X, X0),
  332	poly:div_mono_mono_term(M, N, N0),
  333	!,
  334	poly:lcm(A, B, L),
  335	A0  is L // A,
  336	B0  is L // B,
  337	poly:mul_scalar_poly(A0, X0, Q),
  338	poly:mul_mono_poly(B0*N0, Y, Y0),
  339	poly:subtract_poly_poly(Q, Y0, R, Ord).
  340
  341normal_poly(true, Prime, Poly, R):- !,
  342	(	Prime == 0
  343	->	once(poly:remove_content_poly(Poly, R))
  344	;	once(poly:galois_poly(Prime, Poly, R))
  345	).
  346normal_poly(_, _, Poly, R):- once(poly:canonical_form(Poly, R))