# # # multivariate polynomial algorithms developed by myself # get_null_space([p,q,r,s,...],[x,y,z,...]); # # p,q,r,s,... are mv polynomials in the variables x,y,z,... # this routine returns a basis for the null space in terms of the # polynomials p,q,r,s,... gpolyutil[get_null_space]:=proc(polys::list,vars::list) local all_exponents,i,j; all_exponents := [ op( { seq( op( gpolyutil[mvexps](polys[i],vars)), i = 1 .. nops(polys) ) } ) ]; seq( [ seq( gpolyutil[mvcoeff]( expand(polys[j]), all_exponents[i], vars ), i=1..nops(all_exponents) ) ], j = 1 .. nops(polys)); linalg[matrix]([%]); linalg[transpose](%); linalg[kernel](%) end: ### 1998-12-9 gpolyutil[simplify_poly_table_proc]:= proc( dummy_poly_table::name, poly_table::name, variable_list::list ) local d, idx_param, p, new_poly_table, new_dummy_poly_table, new_poly_list, new_dummy_poly_list, this_dummy_pol, this_sum, hgs_too_high_list, hgs_too_high_index, loop_param, internal_dpt, solutions_list, calculate_hgs_part, hgs_part, hgs_to_onevar, hgs_degree, debug_print; ### subprocedures debug_print:=proc(data) print(data) end: hgs_degree:=proc(poly::polynom, v_list::list) local w, idx_param; subs(seq(v_list[idx_param]=w*v_list[idx_param], idx_param=1..nops(v_list)),poly); collect(expand(%),w); ### WARNING: degree(0,x) now returns -infinity degree(%,w); end: calculate_hgs_part:=proc(deg::integer, poly::polynom, v_list::list) local w, idx_param; subs(seq(v_list[idx_param]=w*v_list[idx_param], idx_param=1..nops(v_list)),poly); collect(expand(%),w); coeff(%,w,deg); end: ### initialization if nops([indices(poly_table)])=0 then ERROR(`null polynomial table supplied`) fi; if type(dummy_poly_table,`table`) then internal_dpt:=eval(dummy_poly_table); # error check if nops([indices(dummy_poly_table)]) <> 1+nops([indices(poly_table)]) then ERROR(`mismatch in lengths of tables`); fi; else if not type(poly_table,`table`) then ERROR(`poly_table should contain a table of polys`) fi; internal_dpt:=table([0=p,seq(op(idx_param)=p[op(idx_param)], idx_param=indices(poly_table))]) fi; hgs_too_high_list:=[]; hgs_too_high_index:=[]; new_poly_list:=[]; new_dummy_poly_list:=[]; solutions_list:=[]; nops([indices(poly_table)]); debug_print(`polynomials remaining ` . %); ### calculate the largest hgs degree d:=max(seq(hgs_degree(expand(idx_param[1]),variable_list), idx_param=entries(poly_table))); debug_print(`max degree is ` . d); for loop_param in indices(poly_table) do hgs_degree(expand(poly_table[op(loop_param)]),variable_list); if ( % < d ) then # list the polys with degree less than d new_poly_list:=[op(new_poly_list),poly_table[op(loop_param)]]; new_dummy_poly_list:=[op(new_dummy_poly_list), eval(internal_dpt[op(loop_param)])]; else if ( % = d ) then # list the polys with degree = d hgs_part:=calculate_hgs_part(d,expand(poly_table[op(loop_param)]), variable_list); hgs_too_high_list := [ op(hgs_too_high_list), hgs_part ]; hgs_too_high_index:= [ op(hgs_too_high_index), op(loop_param) ]; else ERROR(`polynomial has degree higher than highest degree...`) fi fi; od; nops(hgs_too_high_list); debug_print(`number of polys at that degree ` . %); ### now simplify the hgs_too_high_list gpolyutil[get_null_space](hgs_too_high_list,variable_list); # " is a set of VECTORs for loop_param in % do # create new dummy poly corresponding to basis element this_dummy_pol:=sum('loop_param[idx_param]* internal_dpt[hgs_too_high_index[idx_param]]', 'idx_param'=1..nops(hgs_too_high_list)); # now check whether solution or just simpler: this_sum:=sum('loop_param[idx_param]* poly_table[hgs_too_high_index[idx_param]]', 'idx_param'=1..nops(hgs_too_high_list)); if simplify(this_sum) = 0 then # a solution debug_print(`is a solution`); solutions_list:=[op(solutions_list),this_dummy_pol]; else debug_print(`removes hgs part`); # not a solution - we just knocked the top degree hgs part off new_dummy_poly_list:=[op(new_dummy_poly_list),this_dummy_pol]; new_poly_list:=[op(new_poly_list),this_sum]; fi; od; ### results new_dummy_poly_table:=table([0=internal_dpt[0], seq(idx_param=new_dummy_poly_list[idx_param], idx_param=1..nops(new_dummy_poly_list))]); new_poly_table:=table([seq(idx_param=new_poly_list[idx_param], idx_param=1..nops(new_poly_list))]); # call by reference dummy_poly_table:=new_dummy_poly_table; poly_table:=new_poly_table; # return value solutions_list; end: