###### RATH1 ###### rath1:=module() export main1; global rath1; local extract_function,solve1,fdlst,find_m,coeff1,stdsol,prin,sol1; option package; main1:=proc(eqlist::list) local eq1; global `rath1/eta`,`rath1/m`,`rath1/eq`,lmm,`rath1/lce`,`rath1/lkc`,`rath1/c_par`,`rath1/c_key`, `rath1/settime`; `rath1/settime`:=time(); `rath1/m`:='`rath1/m`'; `rath1/c_par`:=0; `rath1/c_key`:=0; if nargs=2 then `rath1/c_par`:=1; elif nargs=3 then if args[2]=0 then `rath1/c_key`:=1; else `rath1/c_par`:=1; `rath1/c_key`:=1; fi; fi; if nargs<1 then ERROR(` Too few arguments ! `); elif nargs>3 then ERROR(` Too many arguments ! `); elif nops(eqlist)<1 then ERROR(` You must specify a list of equations ! `); elif nops(eqlist)>1 then ERROR(` Too many equations in equation list ! `); fi; if evalb(not(type(op(eqlist),equation))) then ERROR(` The input equation must be in the form of equality ! `); fi; `rath1/eq`:=eqlist; printf(`\n`); printf(`The input equation is:`); print(op(`rath1/eq`)); eq1:=numer(lhs(op(`rath1/eq`)))=rhs(op(`rath1/eq`))*denom(lhs(op(`rath1/eq`))); `rath1/eq`:=[eq1]; solve1(); end: extract_function:=proc(l1::list) local i,tset1,tset2,lw,varp1,opset; global `rath1/name_of_function`,`rath1/varl`,`rath1/func`,`rath1/lce`,`rath1/power_var_list`, `rath1/fraction_list`,`rath1/tset`,`rath1/varnum`; lw:=1; tset1:=l1; tset2:={NULL}; `rath1/lce`:={NULL}; varp1:={NULL}; `rath1/power_var_list`:={NULL}; `rath1/fraction_list`:={NULL}; while evalb(lw>0) do lw:=0; for i from 1 to nops(tset1) do if type(op(i,tset1),`+`) then tset2:=tset2 union {op(op(i,tset1))}; lw:=1; elif type(op(i,tset1),`*`) then tset2:=tset2 union {op(op(i,tset1))}; lw:=1; elif type(op(i,tset1),`^`) then tset2:=tset2 union {op(op(i,tset1))}; if evalb(type(op(1,op(i,tset1)),function) and type(op(2,op(i,tset1)),name)) then `rath1/power_var_list`:=`rath1/power_var_list` union {op(2,op(i,tset1))}; elif evalb(type(op(2,op(i,tset1)),function) and type(op(1,op(i,tset1)),name)) then `rath1/power_var_list`:=`rath1/power_var_list` union {op(1,op(i,tset1))}; elif evalb(type(op(1,op(i,tset1)),function) and type(op(2,op(i,tset1)),fraction)) then `rath1/fraction_list`:=`rath1/fraction_list` union {op(2,op(i,tset1))}; elif evalb(type(op(1,op(i,tset1)),fraction) and type(op(2,op(i,tset1)),function)) then `rath1/fraction_list`:=`rath1/fraction_list` union {op(1,op(i,tset1))}; fi; lw:=1; elif has(op(i,tset1),diff) then tset2:=tset2 union {op(op(i,tset1))}; lw:=1; else tset2:=tset2 union {op(i,tset1)}; fi; od; tset1:=remove(type,tset2,constant); varp1:=select(type,tset1,name); `rath1/lce`:=`rath1/lce` union varp1; tset2:=remove(type,tset1,name); tset1:=tset2; tset2:={NULL}; od; `rath1/func`:=op(tset1); opset:={op(`rath1/func`)}; `rath1/tset`:={op(nops(opset),`rath1/func`)}; `rath1/name_of_function`:=op(0,`rath1/func`); `rath1/varl`:={op(`rath1/func`)}; `rath1/varnum`:=nops(`rath1/varl`); `rath1/lce`:=`rath1/lce` minus `rath1/varl`; `rath1/lce`:=`rath1/lce` minus `rath1/power_var_list`; end: solve1:=proc() local i,j,lrp1,l1,lr,lc,ll1,ll2,ll3,lk1,lrp,eqq,flag,l7,w; global `rath1/eq`,`rath1/tset`,`rath1/lce`,`rath1/lkc`,`rath1/power_var_list`,`rath1/fraction_list`,`rath1/varl`, `rath1/eta`, `rath1/k`,`rath1/varnum`; w:=array(1..20); eqq:=expand(lhs(op(`rath1/eq`))); l1:=[op(eqq)]; eqq:=expand(rhs(op(`rath1/eq`))); lr:={op(eqq)}; extract_function(l1); lrp:=remove(type,lr,constant); `rath1/lkc`:={NULL}; lrp1:=select(hastype,lrp,function); lr:=lr minus lrp1; lrp1:=convert(lrp1,list); lrp1:=expand(-1*lrp1); l1:=[op(l1),op(lrp1)]; flag:=1; lrp1:={NULL}; while evalb(flag>0) do flag:=0; for i from 1 to nops(lrp) do if type(op(i,lrp),`+`) then lrp1:=lrp1 union {op(op(i,lrp))}; flag:=1; elif type(op(i,lrp),`*`) then lrp1:=lrp1 union {op(op(i,lrp))}; flag:=1; elif type(op(i,lrp),`^`) then lrp1:=lrp1 union {op(op(i,lrp))}; flag:=1; elif has(op(i,lrp),diff) then lrp1:=lrp1 union {op(op(i,lrp))}; flag:=1; else lrp1:=lrp1 union {op(i,lrp)}; fi; od; lrp:=lrp1; lrp1:={NULL}; od; lrp:=select(type,lrp,name); `rath1/lce`:=`rath1/lce` union lrp; `rath1/lce`:=`rath1/lce` minus `rath1/varl`; lc:=`rath1/tset` intersect `rath1/varl`; `rath1/varl`:=`rath1/varl` minus lc; l7:=convert(`rath1/varl`,list); l7:=sort(l7); l7:=[op(l7),op(lc)]; if nops(l7)=1 then else ll1:=[`rath1/k`]; for i from 1 to `rath1/varnum`-1 do ll1:=[op(ll1),-`rath1/k`*w[i]]; od; ll3:=[]; for j from 1 to nops(l7) do ll3:=[op(ll3),op(j,ll1)]; od; ll2:=[]; for i from 1 to nops(ll3) do ll2:=[op(ll2),op(i,ll3)*op(i,l7)]; od; `rath1/lkc`:=[op(1,ll1)]; for i from 2 to nops(l7) do lk1:=simplify(op(i,ll1)/op(1,ll1)); `rath1/lkc`:=[op(`rath1/lkc`),-lk1]; od; `rath1/lkc`:=convert(`rath1/lkc`,set); `rath1/eta`:=add(j,j=ll2); fi; fdlst(l1,lr); end: fdlst:=proc(l1::list,lr::set) local j,eqm,cpu_time,lf3,equ1,equ2,equ3,equ,mm,ss,i,lf1,lf2,lf4,lf5,lf6,dnm,lf7,lf8; global `rath1/eta`,`rath1/xi`,`rath1/m`,`rath1/func`,`rath1/power_con`,`rath1/max_rem`,`rath1/ddm`,`rath1/u`; lf1:=[]; if nops(`rath1/func`)>1 then for i from 1 to nops(l1) do equ:=op(i,l1); equ1:=simplify(subs(`rath1/func`=`rath1/u`(`rath1/eta`),equ)); equ2:=simplify(subs(`rath1/eta`=`rath1/xi`,equ1)); equ3:=convert(equ2,diff,`rath1/xi`); lf1:=[equ3,op(lf1)]; od; else lf1:=[]; lf8:=op(`rath1/func`); for i from 1 to nops(l1) do lf7:=simplify(algsubs(lf8=`rath1/xi`,op(i,l1))); lf7:=simplify(subs(`rath1/name_of_function`=`rath1/u`,lf7)); lf1:=[op(lf1),lf7]; od; fi; j:=0; eqm:=add(j,j=lf1); for i from 1 to nops(lr) do eqm:=eqm-op(i,lr); od; mm:=20; ss:=`rath1/u`(`rath1/xi`)=tanh(`rath1/xi`)^mm; if `rath1/fraction_list`<>{NULL} then dnm:=denom(`rath1/fraction_list`); `rath1/ddm`:=ilcm(op(dnm)); lf2:=radsimp(subs(`rath1/u`(`rath1/xi`)=`rath1/v`(`rath1/xi`)^`rath1/ddm`,lf1)); eqm:=radsimp(subs(`rath1/u`(`rath1/xi`)=`rath1/v`(`rath1/xi`)^`rath1/ddm`,eqm)); eqm:=expand(eqm); eqm:=simplify(subs(`rath1/v`(`rath1/xi`)=`rath1/u`(`rath1/xi`),eqm)); lf2:=expand(lf2); lf2:=simplify(subs(`rath1/v`(`rath1/xi`)=`rath1/u`(`rath1/xi`),lf2),power); else lf2:=lf1; fi; `rath1/power_con`:={NULL}; if `rath1/power_var_list`<>{NULL} then lf1:=select(has,lf2,op(`rath1/power_var_list`)); lf2:=convert(lf2,set) minus convert(lf1,set); lf3:={NULL}; for i from 1 to nops(lf1) do if type(op(i,lf1),`*`) then lf3:=lf3 union {op(op(i,lf1))}; else lf3:=lf3 union {op(i,lf1)}; fi; od; if `rath1/fraction_list`<>{NULL} then lf4:=select(hastype,lf3,(function^integer)^algebraic); for i from 1 to nops(lf4) do lf5:={op(op(i,lf4))}; if hastype(op(1,lf5),function^integer) then lf6:=select(type,{op(op(1,lf5))},integer); lf1:=simplify(subs((`rath1/u`(`rath1/xi`)^op(lf6))^op(2,lf5)=`rath1/u`(`rath1/xi`)^(op(lf6) *op(2,lf5)),lf1),power); else lf6:=select(type,{op(op(2,lf5))},integer); lf1:=simplify(subs((`rath1/u`(`rath1/xi`)^op(lf6))^op(1,lf5)=`rath1/u`(`rath1/xi`)^(op(lf6) *op(1,lf5)),lf1),power); fi; od; lf3:=select(has,lf1,op(`rath1/power_var_list`)); lf4:={NULL}; for i from 1 to nops(lf3) do if type(op(i,lf3),`*`) then lf4:=lf4 union {op(op(i,lf3))}; else lf4:=lf4 union {op(i,lf3)}; fi; od; lf5:=select(has,lf4,op(`rath1/power_var_list`)); for i from 1 to nops(lf5) do lf6:={op(op(i,lf5))}; if hastype(op(1,lf6),function) then `rath1/power_con`:=`rath1/power_con` union {op(2,lf6)}; else `rath1/power_con`:=`rath1/power_con` union {op(1,lf6)}; fi; od; else lf4:=select(has,lf3,op(`rath1/power_var_list`)); for i from 1 to nops(lf4) do lf5:={op(op(i,lf4))}; if hastype(op(1,lf5),function) then `rath1/power_con`:=`rath1/power_con` union {op(2,lf5)}; else `rath1/power_con`:=`rath1/power_con` union {op(1,lf5)}; fi; od; fi; lf2:=convert(lf2,set) union convert(lf1,set); lf2:=simplify(subs(op(`rath1/power_var_list`)=1,lf2),power); fi; lf3:=collect(expand(simplify(subs(ss,lf2))),tanh(`rath1/xi`)); lf4:={NULL}; for i from 1 to nops(lf3) do lf4:=lf4 union {degree(op(i,lf3),tanh(`rath1/xi`))}; od; lf4:=lf4 minus {0}; `rath1/m`:='`rath1/m`'; lf5:={NULL}; for i from 1 to nops(lf4) do lf5:=lf5 union {iquo(op(i,lf4),mm)*`rath1/m`+irem(op(i,lf4),mm)}; od; lf3:=convert(convert(lf5,set),list); `rath1/max_rem`:=rem(op(1,lf3),`rath1/m`,`rath1/m`); for i from 2 to nops(lf3) do `rath1/max_rem`:=max(`rath1/max_rem`, rem(op(i,lf3),`rath1/m`,`rath1/m`)); od; find_m(lf3); if evalb(`rath1/m`<>0) then coeff1(eqm); else printf(` No non-trivial solution in terms of tanh-function!\n `); cpu_time:=time()-`rath1/settime`; printf(`The running time is:%a seconds \n`, cpu_time); fi; end: find_m:=proc(lf3::list) local i,hh,m1_max,m2_max,m3_max,m_ineq,m_eq,m_value,max_m_value; global `rath1/m`,`rath1/max_rem`; `rath1/max_rem`:=op(1,lf3)-coeff(op(1,lf3),`rath1/m`)*`rath1/m`; for i from 2 to nops(lf3) do `rath1/max_rem`:=max(`rath1/max_rem`,op(i,lf3)-coeff(op(i,lf3),`rath1/m`)*`rath1/m`); od; m1_max:=max(op(lf3)); m2_max:={convert(m1_max,piecewise)}; m3_max:={op(op(m2_max))}; m_ineq:={NULL}; for i from 1 to nops(m3_max) do if not (type(op(i,m3_max),`+`) or type(op(i,m3_max),`*`)) then m_ineq:=m_ineq union {op(i,m3_max)}; fi; od; m_ineq:=remove(type,m_ineq,constant); m_ineq:=remove(type,m_ineq,name); m_eq:=map(convert,m_ineq, equality); m_value:={NULL}; for i from 1 to nops(m_eq) do if type(lhs(op(i,m_eq)),constant) then hh:=rhs(op(i,m_eq))=lhs(op(i,m_eq)); m_value:=m_value union {hh}; else m_value:=m_value union {op(i,m_eq)}; fi; od; max_m_value:=-100; for i from 1 to nops(m_value) do max_m_value:=max(max_m_value,rhs(op(i,m_value))); od; `rath1/m`:=max_m_value; printf(`The order of the solution is: m=%a`,`rath1/m`); printf(`\n`); end: stdsol:=proc(set1::set,set2::set,sss::set) local i,j,lyy,luy,vsol,lsol,ssol; global `rath1/lll`; `rath1/lll`:=sss; vsol:={NULL}; for i from 1 to nops(`rath1/lll`) do luy:=op(i,`rath1/lll`); if has(lhs(luy),[op(set2)]) then for j from 1 to nops(set1) do if has(rhs(luy),op(j,set1)) then vsol:={solve(luy,{op(j,set1)})}; lyy:=luy; break; fi; od; fi; od; if vsol<>{NULL} then `rath1/lll`:=`rath1/lll` minus {lyy}; lsol:={NULL}; for j from 1 to nops(vsol) do ssol:=simplify(subs(op(j,vsol),`rath1/lll`)); ssol:=ssol union op(j,vsol); lsol:=lsol union {ssol}; od; `rath1/lll`:=lsol; else `rath1/lll`:={`rath1/lll`}; fi; ssol:={NULL}; for i from 1 to nops(`rath1/lll`) do if evalb(SearchText(`_`,convert(op(i,`rath1/lll`),string))<>0) then ssol:=ssol union {op(i,`rath1/lll`)}; fi; od; `rath1/lll`:=`rath1/lll` minus ssol; end: prin:=proc(lae::set,lm::set,a::array) local i,j,zx,lhs1,lhs2,vs,lde1,lde2,le1,le2,luy,xx,lco,ji,ij,luy1, hh,lmm1,lmm2,lmm3,jj,llu,denom2,luy2,dlae,cpu_time,lve,lmm,lhs3, nps,vp,lmm4,xi,lf8,_k; global `rath1/lce`,`rath1/pis1`,`rath1/eta`,`rath1/mn`,`rath1/func`,`rath1/k`,`rath1/u`,`rath1/lkc`; interface(labelling=false); lmm:=lm; xi:=0; _k:=0; #print(nops(lmm),`print_module_start`); lve:=lae; nps:={seq(nops(i),i=lmm)}; if evalb( nops(nps)>1 and nops(lmm)>12) then for i from 1 to nops(nps) do vp[i]:={NULL}; od; for j from 1 to nops(nps) do for i from 1 to nops(lmm) do if nops(op(i,lmm))=op(j,nps) then vp[j]:=vp[j] union {op(i,lmm)}; fi; od; lmm:=lmm minus vp[j]; od; lmm:=[]; for i from 1 to nops(nps) do lmm:=[op(lmm),op(vp[i])] od; fi; lmm1:={NULL}; lmm2:=convert(lmm,set); for i from 1 to nops(lmm) do luy:=op(i,lmm); if not(member(luy,lmm2)) then next; fi; denom2:=1; for xx from 1 to nops(luy) do denom2:=denom2*denom(rhs(op(xx,luy))); od; for j from 1 to nops(lmm2) do luy1:=op(j,lmm2); if evalb(luy<>luy1) then hh:=simplify(subs(luy1,denom2)); if hh <> 0 then lhs1:=simplify(subs(luy1,luy),assume=real); if has(lhs1,{csgn,signum,sign}) then luy2:=simplify(subs({csgn=1,signum=1,sign=1},lhs1)); lhs2:={NULL}; for jj from 1 to nops(luy2) do llu:=simplify(lhs(op(jj,luy2))-rhs(op(jj,luy2)), assume=real); lhs2:=lhs2 union {llu}; od; luy2:=simplify(subs({csgn=-1,signum=-1,sign=-1},lhs1)); lhs3:={NULL}; for jj from 1 to nops(luy2) do llu:=simplify(lhs(op(jj,luy2))-rhs(op(jj,luy2)), assume=real); lhs3:=lhs3 union {llu}; od; if lhs2={0} then lmm1:=lmm1 union {op(j,lmm2)}; elif lhs3={0} then lmm1:=lmm1 union {op(j,lmm2)}; else next; fi; else lhs3:={NULL}; for jj from 1 to nops(lhs1) do llu:=simplify(lhs(op(jj,lhs1))-rhs(op(jj,lhs1)), assume=real); lhs3:=lhs3 union {llu}; od; if lhs3={0} then lmm1:=lmm1 union {op(j,lmm2)}; else next; fi; fi; else next; fi; fi; od; lmm2:=lmm2 minus lmm1; lmm1:={NULL}; od; lmm:=lmm2 minus lmm1; if evalb(`rath1/c_key`=0) then lmm1:=select(has,lmm,I); lmm:=lmm minus lmm1; fi; lmm2:=lmm; dlae:={NULL}; for i from 1 to `rath1/mn` by 2 do dlae:=dlae union {a[i]}; od; dlae:=dlae union {`rath1/k`}; lmm3:={NULL}; for i from 1 to nops(lmm)-1 do luy:=op(i,lmm); lco:={NULL}; for j from 1 to nops(luy) do vs:=lhs(op(j,luy)); lco:=lco union {vs}; od; if evalb(lco intersect {`rath1/k`}<>{NULL}) then luy1:={NULL}; for ij from 1 to nops(luy) do vs:=lhs(op(ij,luy)); if evalb({vs} intersect dlae<>{NULL}) then luy1:=luy1 union {vs=-rhs(op(ij,luy))}; else luy1:=luy1 union {op(ij,luy)}; fi; od; for ji from i+1 to nops(lmm2) do if op(ji,lmm2)=luy1 then lmm3:=lmm3 union {op(ji,lmm2)}; fi; od; else next; fi; od; lmm4:=lmm minus lmm3; xi:='xi'; _k:='_k'; if evalb(nops(lmm4)>0)then printf(`%a required solution(s) is(are) listed as follows `,nops(lmm4)); if nops({op(`rath1/func`)})=1 then lf8:=op(`rath1/func`); #`rath1/pis1`=subs(`rath1/u`=`rath1/name_of_function`,`rath1/pis1`); `rath1/pis1`:=subs(`rath1/xi`=lf8,`rath1/pis1`); `rath1/pis1`:=subs(`rath1/k`=_k,`rath1/pis1`); print(subs(`rath1/xi`=xi,`rath1/pis1`)); else #`rath1/pis1`=subs(`rath1/u`=`rath1/name_of_function`,`rath1/pis1`); `rath1/pis1`:=subs(`rath1/k`=_k,`rath1/pis1`); `rath1/eta`:=subs(`rath1/k`=_k,`rath1/eta`); print(subs(`rath1/xi`=xi,`rath1/pis1`), xi=factor(`rath1/eta`)); fi; print(`where:`); else printf(`No non-trivial solution in terms of tanh-function!`); fi; print(` `); lmm3:=subs(`rath1/k`=_k,lmm4); lmm4:=lmm3; lve:=subs(`rath1/k`=_k,lve); `rath1/lkc`:=subs(`rath1/k`=_k,`rath1/lkc`); `rath1/lce`:=subs(`rath1/k`=_k,`rath1/lce`); for i from 1 to nops(lmm4) do lhs1:=op(i,lmm4); lhs2:={NULL}; for j from 1 to nops(lhs1) do vs:=lhs(op(j,lhs1)); lhs2:=lhs2 union {vs}; od; lde1:=(lve union `rath1/lkc`) minus lhs2; lde2:=`rath1/lce` minus lhs2; le1:=[]; for j from 1 to nops(lhs1) do vs:={lhs(op(j,lhs1))}; if evalb(vs intersect lve <> {NULL})then le1:=[op(le1),op(j,lhs1)]; elif evalb(vs intersect `rath1/lkc` <> {NULL}) then le1:=[op(le1),op(j,lhs1)]; fi; od; le1:=convert(le1,set); le2:=lhs1 minus le1; print(le1); for j from 1 to nops(lde1) do printf(`%a :free `,op(j,lde1)); od; if lde1<>{NULL} then printf(`\n`); fi; if evalb(nops(le2)>0) then printf(`The parameters constraints are: \n`); print(le2); fi; for zx from 1 to nops(lde2) do printf(`%a : free `,op(zx,lde2)); od; print(`\n`); od; cpu_time:=time()-`rath1/settime`; printf(`The running time is:%a seconds \n`, cpu_time); end: sol1:=proc(pp1::list,pp2::list,lae::set,a::array) local tmm1,tmm3,ps3,psr,psr1,psr2,i,j,optm,optm1,luu, jj,sd,pse,pps,pss1,ii,ji,dsb,dss,opt,mmt,lco, luy,spp,pse1,kce,mmt1,ae1,lpt,lkc2, pps1,lmm1,tmm,luy1,rra,rrb,rrc,ps,ps1,lve,ae2,lvc,lvc2,lw, cel,lm,lmm; global `rath1/pis1`,`rath1/lce`,`rath1/c_key`,`rath1/c_par`,`rath1/k`,charsets; interface(labelling=false); # read `d:\\yliu\\charsets2.1\\charsets`; with(charsets); lmm:={NULL}; ps1:=pp1; ps:=pp2; lve:=lae; cel:=sort(convert(`rath1/lce`,list),`lexorder`); ae1:=[]; for i from 0 to `rath1/mn` do for j from 1 to nops(lve) do if i=op(op(j,lve)) then ae1:=[op(ae1),op(j,lve)]; fi; od; od; ae2:=[]; for i from `rath1/mn` to 0 by -1 do for j from 1 to nops(lve) do if i=op(op(j,lve)) then ae2:=[op(ae2),op(j,lve)]; fi; od; od; lkc2:={NULL}; if evalb(`rath1/c_par`=0) then for i from 1 to nops(`rath1/lce`) do lkc2:=lkc2 union {op(i,`rath1/lce`)=0}; od; fi; for i from 1 to nops(`rath1/lkc`) do lkc2:=lkc2 union {op(i,`rath1/lkc`)=0}; od; lkc2:=lkc2 union {a[`rath1/mn`]=0}; kce:={NULL}; for j from 1 to nops(`rath1/lkc`) do kce:=kce union {op(j,`rath1/lkc`)=op(j,`rath1/lkc`)}; od; for j from 1 to nops(`rath1/lce`) do kce:=kce union {op(j,`rath1/lce`)=op(j,`rath1/lce`)}; od; for j from 1 to nops(lve) do kce:=kce union {op(j,lve)=op(j,lve)}; od; if nops(ps1)<=2 then if nops(`rath1/lce`)<=1 then tmm:=simplify(expand(rationalize({csolve(ps1,[op(ae1),op(cel)])}))); else tmm:=simplify(expand(rationalize({csolve(ps1,[op(cel),op(ae1)])}))); fi; else tmm:=simplify(expand(rationalize({csolve(ps1,[op(cel),op(ae2)])}))); fi; lmm1:=simplify(tmm,assume=real); tmm:=simplify(subs({signum=1,csgn=1,sign=1},lmm1)); mmt:=simplify(subs({signum=-1,csgn=-1,sign=1},lmm1)); tmm:=tmm union mmt; mmt:={NULL}; for i from 1 to nops(tmm) do luy:=op(i,tmm); luy1:=radsimp(luy,`ratdenom`); if has(luy1, I) then mmt:=mmt union {luy}; else mmt:=mmt union {luy1}; fi; od; tmm:=mmt; mmt:={NULL}; for i from 1 to nops(tmm) do luy:=op(i,tmm); if luy intersect lkc2<>{NULL} then mmt:=mmt union {luy}; fi; od; tmm:=tmm minus mmt; mmt:={NULL}; for i from 1 to nops(tmm) do if length(op(i,tmm))>1500 then mmt:=mmt union {op(i,tmm)}; fi; od; tmm:=tmm minus mmt; mmt:={NULL}; if evalb(`rath1/c_key`=0) then mmt:=select(has,tmm,I); tmm:=tmm minus mmt; fi; mmt:={NULL}; for i from 1 to nops(tmm) do luy:=op(i,tmm); ps3:=simplify(subs(luy,ps1)); if convert(ps3,set)={0} then mmt:=mmt union {luy}; elif hasfun(ps3,abs) then mmt:=mmt union {luy}; fi; od; tmm:=mmt; for i from 1 to nops(tmm) do optm:=op(i,tmm); lco:={NULL}; for j from 1 to nops(optm) do lco:=lco union {lhs(op(j,optm))}; od; psr:=convert(ps,set); psr1:=simplify(subs(optm,psr)); psr1:=psr1 minus {0}; if evalb(psr1={NULL}) then lmm:=lmm union {optm}; next; fi; psr2:=numer(psr1); pss1:={NULL}; for ii from 1 to nops(psr2)-1 do for j from ii+1 to nops(psr2) do if evalb(ii<>j) then if evalb(op(ii,psr2)+op(j,psr2)=0) then pss1:=pss1 union {op(j,psr2)}; fi; fi; od; od; pps:=psr2 minus pss1; spp:={NULL}; for ii from 1 to nops(pps) do pse:=op(ii,pps); if type(pse,`*`) then pse1:={op(pse)}; else pse1:={pse}; fi; for j from 1 to nops(pse1) do if type(op(j,pse1),constant) then pse:=simplify(pse/op(j,pse1)); fi; od; spp:=spp union {pse}; od; pps:=spp; if hastype(pps,radical) then lpt:=1; elif length(pps)>=4000 then lpt:=1; elif evalb(length(pps)/nops(pps)>=500) then lpt:=1; fi; _EnvExplicit:=true; if evalb(lpt=1) then if nops(pps)=1 then tmm1:={solve(pps,{op(`rath1/lkc`)})}; else tmm1:={solve(pps,{op(`rath1/lce` minus lco),op(`rath1/lkc`),op(lve minus lco)})}; fi; mmt:={NULL}; for j from 1 to nops(tmm1) do mmt1:=op(j,tmm1) minus kce; mmt:=mmt union {mmt1}; od; tmm1:=mmt; mmt:={NULL}; for j from 1 to nops(tmm1) do if op(j,tmm1) intersect lkc2 <> {NULL} then mmt:=mmt union {op(j,tmm1)}; fi; od; tmm1:=tmm1 minus mmt; if evalb(`rath1/c_key`=0) then lmm1:={NULL}; for j from 1 to nops(tmm1) do luy:=op(j,tmm1); if has(luy,I) then lmm1:=lmm1 union {luy}; fi; od; tmm1:=tmm1 minus lmm1; fi; mmt:=rationalize(tmm1); tmm1:=simplify(expand(mmt)); mmt:={NULL}; for j from 1 to nops(tmm1) do luy:=op(j,tmm1); if length(luy)>2500 then mmt:=mmt union {luy}; fi; if evalb(SearchText(`_`,convert(luy,string))<>0) then mmt:=mmt union {luy}; fi; od; tmm1:=tmm1 minus mmt; else pps:=numer(pps); pps1:=factor(pps); if length(pps1)>=2000 then tmm1:={csolve(pps1,{`rath1/k`,op(`rath1/lce` minus lco),op(`rath1/lkc` minus {`rath1/k`}),op(lve minus lco)})}; else lvc:={op(lve minus lco)}; lvc2:=[]; for ii from `rath1/mn` to 0 by -1 do for ji from 1 to nops(lvc) do if ii=op(op(ji,lvc)) then lvc2:=[op(lvc2),op(ji,lvc)]; fi; od; od; tmm1:={csolve(pps1,[`rath1/k`,op(`rath1/lce` minus lco),op(`rath1/lkc` minus {`rath1/k`}),op(lvc2)])}; fi; mmt:=simplify(tmm1); tmm1:=simplify(expand(mmt)); mmt:={NULL}; for j from 1 to nops(tmm1) do luy:=op(j,tmm1); if length(luy)>2500 then mmt:=mmt union {luy}; fi; if evalb(SearchText(`_`,convert(luy,string))<>0) then mmt:=mmt union {luy}; fi; od; tmm1:=tmm1 minus mmt; mmt:={NULL}; for j from 1 to nops(tmm1) do if op(j,tmm1) intersect lkc2 <> {NULL} then mmt:=mmt union {op(j,tmm1)}; fi; od; tmm1:=tmm1 minus mmt; if evalb(`rath1/c_key`=0) then tmm3:=select(has,tmm1,I); tmm1:=tmm1 minus tmm3; fi; fi; for j from 1 to nops(tmm1) do dsb:=op(j,tmm1); optm1:={NULL}; for ji from 1 to nops(optm) do sd:=denom(rhs(op(ji,optm))); optm1:=optm1 union {sd}; od; optm1:=simplify(subs(dsb,optm1)); if evalb(optm1 intersect {0}={NULL}) then opt:=simplify(subs(dsb,optm)); dss:=dsb union opt; lmm:=lmm union {dss}; else next; fi; od; lpt:=0; od; mmt:={NULL}; for i from 1 to nops(lmm) do if length(op(i,lmm))>5000 then mmt:=mmt union {op(i,lmm)}; fi; od; lmm:=lmm minus mmt; if evalb(`rath1/c_key`=0) then mmt:=select(has,lmm,I); lmm:=lmm minus mmt; fi; lw:=1; while(lw>0) do lw:=0; lmm1:={NULL}; for i from 1 to nops(lmm) do stdsol(lae,`rath1/lce`,op(i,lmm)); if {op(i,lmm)}<>`rath1/lll` then lmm1:=lmm1 union `rath1/lll`; lw:=1; else lmm1:=lmm1 union {op(i,lmm)}; fi; od; lmm:=lmm1; od; lw:=1; while(lw>0) do lw:=0; lmm1:={NULL}; for i from 1 to nops(lmm) do stdsol(`rath1/lkc`,`rath1/lce`,op(i,lmm)); if {op(i,lmm)}<>`rath1/lll` then lmm1:=lmm1 union `rath1/lll`; lw:=1; else lmm1:=lmm1 union {op(i,lmm)}; fi; od; lmm:=lmm1; od; lw:=1; while(lw>0) do lw:=0; lmm1:={NULL}; for i from 1 to nops(lmm) do stdsol(lae,`rath1/lkc`,op(i,lmm)); if {op(i,lmm)}<>`rath1/lll` then lmm1:=lmm1 union `rath1/lll`; lw:=1; else lmm1:=lmm1 union {op(i,lmm)}; fi; od; lmm:=lmm1; od; mmt:={NULL}; for i from 1 to nops(lmm) do mmt:=mmt union {op(i,lmm) minus kce }; od; lmm:=mmt; mmt:=expand(lmm); lmm:=simplify(mmt,assume=real); mmt:=simplify(subs({signum=1,csgn=1,sign=1},lmm)); lmm1:=simplify(subs({signum=-1,csgn=-1,sign=1},lmm)); lmm:=mmt union lmm1; lmm1:={NULL}; for j from 1 to nops(lmm) do luy:=op(j,lmm); luy1:=radsimp(luy,`ratdenom`); if has(luy1,I) then lmm1:=lmm1 union {luy}; else lmm1:=lmm1 union {luy1}; fi; od; lmm:=lmm1; mmt:={NULL}; for i from 1 to nops(lmm) do luu:=factor(simplify(subs(op(i,lmm),ps),assume=real)); if convert(luu,set)={0} then mmt:=mmt union {op(i,lmm)}; elif hasfun(luu,abs) then lmm1:=denom(luu); tmm1:=mul(jj,jj=lmm1); if simplify(tmm1,assume=positive)<>0 then if convert(simplify(luu,assume=positive),set)={0} then mmt:=mmt union {op(i,lmm)}; fi; elif simplify(tmm1,assume=negative)<>0 then if convert(simplify(luu,assume=negative),set)={0} then mmt:=mmt union {op(i,lmm)}; fi; elif hasfun(luu,abs(`+`)) then mmt:=mmt union {op(i,lmm)}; fi; fi; od; lmm:=mmt; mmt:={NULL}; for i from 1 to nops(lmm) do if op(i,lmm) intersect lkc2 <>{NULL} then mmt:=mmt union {op(i,lmm)}; fi; od; lmm:=lmm minus mmt; if evalb(`rath1/c_key`=0) then mmt:=select(has,lmm,I); lmm:=lmm minus mmt; fi; lmm1:={NULL}; for i from 1 to nops(lmm)-1 do luy:=op(i,lmm); for j from i+1 to nops(lmm) do luy1:=op(j,lmm); rra:=luy minus luy1; rrb:=luy1 minus luy; rrc:={NULL}; for jj from 1 to nops(rrb) do rrc:=rrc union {rhs(op(jj,rrb))=lhs(op(jj,rrb))}; od; if evalb(rra=rrc) then lmm1:=lmm1 union {luy1}; fi; od; od; lmm:=lmm minus lmm1; lm:=lmm; prin(lae,lm,a); end: coeff1:=proc(eqm::algebraic) local i,opi,da,dq,ds,pse,gcd1,lw,pp1,pp2,pis,ps1,sp,n_1,peqn,cpu_time,a, powvar,g,j,pp,ops,gp,pg,ps,lae,eqmm,eqm1,quos,rems,fterm,fterm1,qn; global `rath1/pis1`,`rath1/mn`,`rath1/max_rem`,`rath1/ddm`,`rath1/k`,`rath1/u`,`rath1/v`; eqmm:=numer(eqm); if evalb(`rath1/power_var_list`<>{NULL}) then powvar:=op(`rath1/power_var_list`); if `rath1/power_con` <>{NULL} then for i from 1 to nops(`rath1/power_con`) do quos:=quo(op(i,`rath1/power_con`),powvar,powvar); rems:=rem(op(i,`rath1/power_con`),powvar,powvar); gp[i]:=`rath1/u`(`rath1/xi`)^op(i,`rath1/power_con`)=(`rath1/v`(`rath1/xi`)^quos)*(`rath1/u`(`rath1/xi`)^rems); od; fi; g[0]:=`rath1/u`(`rath1/xi`)=powvar*`rath1/v`(`rath1/xi`); g[1]:=diff(`rath1/u`(`rath1/xi`),`rath1/xi`)=diff(`rath1/v`(`rath1/xi`),`rath1/xi`); if `rath1/power_con`<>{NULL} then for i from 1 to nops(`rath1/power_con`) do eqmm:=simplify(subs(gp[i],eqmm)); od; fi; for j from 2 to `rath1/max_rem` do pp:=algsubs(`rath1/u`(`rath1/xi`)=`rath1/v`(`rath1/xi`)^(1/powvar),diff(`rath1/u`(`rath1/xi`),`rath1/xi`$j)); g[j]:=0; for i from 1 to nops(pp) do ops:=op(i,pp)*powvar*`rath1/v`(`rath1/xi`)*`rath1/v`(`rath1/xi`)^(-1/powvar); g[j]:=g[j]+ops; od; pg[j]:=diff(`rath1/u`(`rath1/xi`),`rath1/xi`$j)=simplify(g[j]); od; for j from `rath1/max_rem` to 2 by -1 do eqmm:=simplify(subs(pg[j],eqmm)): od; eqmm:=simplify(subs(g[1],expand(eqmm))); eqmm:=simplify(subs(g[0],expand(eqmm))); eqmm:=numer(simplify(algsubs(`rath1/v`(`rath1/xi`)=`rath1/u`(`rath1/xi`),eqmm))); fi; n_1:=sign(`rath1/m`)*denom(`rath1/m`); qn:=1/n_1; `rath1/mn`:=`rath1/m`*n_1; `rath1/v`(`rath1/xi`):='`rath1/v`(`rath1/xi`)'; i:='i'; a:=array(0..`rath1/mn`); pis:=`rath1/v`(`rath1/xi`)=sum(a[i]*tanh(`rath1/xi`)^i,i=0..`rath1/mn`); `rath1/pis1`:=`rath1/func`=sum(a[i]*tanh(`rath1/xi`)^i,i=0..`rath1/mn`)^qn; if evalb(`rath1/power_var_list`<>{NULL}) then powvar:=op(`rath1/power_var_list`); `rath1/pis1`:=`rath1/func`=sum(a[i]*tanh(`rath1/xi`)^i,i=0..`rath1/mn`)^(qn*1/powvar); if evalb(`rath1/fraction_list`<>{NULL}) then `rath1/pis1`:=`rath1/func`=sum(a[i]*tanh(`rath1/xi`)^i,i=0..`rath1/mn`)^(qn*1/(powvar*`rath1/ddm`)); fi; else if evalb(`rath1/fraction_list`<>{NULL}) then `rath1/pis1`:=`rath1/func`=sum(a[i]*tanh(`rath1/xi`)^i,i=0..`rath1/mn`)^(qn*1/`rath1/ddm`); fi; fi; lae:={seq(a[i],i=0..`rath1/mn`)}; eqm1:=simplify(algsubs(`rath1/u`(`rath1/xi`)=`rath1/v`(`rath1/xi`)^qn,eqmm)); eqm1:=numer(eqm1); if hastype(eqm1,function^fraction) then eqm1:=eqm1=0; fterm:={op(lhs(eqm1))}; fterm1:={NULL}; lw:=1; while evalb(lw>0) do lw:=0; for i from 1 to nops(fterm) do if type(op(i,fterm),`+`) then fterm1:=fterm1 union {op(op(i,fterm))}; lw:=1; elif type(op(i,fterm),`*`) then fterm1:=fterm1 union {op(op(i,fterm))}; lw:=1; else fterm1:=fterm1 union {op(i,fterm)}; fi; od; fterm:=fterm1; fterm1:={NULL}; od; fterm1:=select(hastype,fterm,function^fraction); eqm1:=isolate(eqm1,op(1,fterm1)); eqm1:=expand((lhs(eqm1)*denom(rhs(eqm1)))^(1/qn) =numer(rhs(eqm1))^(1/qn)); eqm1:=rhs(eqm1)-lhs(eqm1); else eqm1:=numer(eqm1); fi; peqn:=collect(simplify((subs(pis,eqm1))),tanh(`rath1/xi`)); dq:=0; sp:=[]; for i from 1 to nops(peqn) do opi:=op(i,peqn); ds:=degree(opi,tanh(`rath1/xi`)); if ds<>0 then da:=simplify(opi/tanh(`rath1/xi`)^ds); sp:=[op(sp),da]; else dq:=dq+opi; fi; od: if evalb(dq=0) then else sp:=[op(sp),dq]; fi; sp:=factor(sp); if has(sp,FAIL) then ps1:=[]; ps:=[]; else if nops(sp)>1 then gcd1:=gcd(op(1,sp),op(2,sp)); for i from 3 to nops(sp) do gcd1:=gcd(gcd1,op(i,sp)); od; else gcd1:=1; fi; ps:=[]; for i from 1 to nops(sp) do pse:=op(i,sp)/gcd1; while divide(pse,`rath1/k`) do pse:=simplify(pse/`rath1/k`); od; while divide(pse,a[`rath1/mn`]) do pse:=simplify(pse/a[`rath1/mn`]); od; ps:=[op(ps),pse]; od; ps:=factor(ps); lw:=0; if has(ps,a[0]) then lw:=1 fi; ps1:=[]; if lw=0 then for i from 1 to `rath1/mn` do ps1:=[op(ps1),op(i,ps)]; od; else for i from 1 to `rath1/mn`+1 do ps1:=[op(ps1),op(i,ps)]; od; fi; fi; if ps=[] then printf(`No non-trivial solutions in terms of tanh-function !\n`); cpu_time:=time()-`rath1/settime`; printf(`The running time is:%a seconds \n`, cpu_time); else pp1:=ps1; pp2:=ps; sol1(pp1,pp2,lae,a); fi; end: end module: print("saving rath1 in the library found at:", libname[1]); savelib (rath1); ##### charsets2.1 ##### charsets[charset] := proc() `charsets/charset`(args) end: charsets[mcharset] := proc() `charsets/mcharset`(args) end: charsets[charser] := proc() `charsets/charser`(args) end: charsets[mcs] := proc() `charsets/mcs`(args) end: charsets[ecs] := proc() `charsets/ecs`(args) end: charsets[mecs] := proc() `charsets/mecs`(args) end: charsets[ics] := proc() `charsets/ics`(args) end: charsets[qics] := proc() `charsets/qics`(args) end: charsets[eics] := proc() `charsets/eics`(args) end: charsets[ivd] := proc() `charsets/ivd`(args) end: charsets[pid] := proc() `charsets/pid`(args) end: charsets[remset] := proc() `charsets/remset`(args) end: charsets[cfactor] := proc() `charsets/cfactor`(args) end: charsets[iniset] := proc() `charsets/iniset`(args) end: charsets[csolve] := proc() `charsets/csolve`(args) end: charsets[triser] := proc() `charsets/triser`(args) end: `charsets/remset` := proc(ps,as,ord) local ind,i; if nargs <> 3 then ERROR(`wrong number of arguments`) elif nops(ps) < 1 or nops(as) < 1 then ERROR(`no polynomials specified`) elif nops(ord) < 1 then ERROR(`no indeterminates specified`) elif not type(ord,list) then ERROR(ord,`must be a list`) fi; if member(false,map(type,ord,name)) then ERROR(`bad variable list`) fi; ind := 0; for i to nops(as) do if `charsets/class`(as[i],ord) <= ind then ERROR( `second argument must be a non-contradictory (weak-, quasi-) set`) else ind := `charsets/class`(as[i],ord) fi od; if type(ps,{set,list}) then if member(false,map(type,ps,polynom(polynom(rational),ord))) or member(false,map(type,as,polynom(polynom(rational),ord))) then ERROR(`input must be polynomials over Q in`,ord) fi; `charsets/remseta`(ps,as,ord) else if member(false,map(type,as,polynom(polynom(rational),ord))) then ERROR(`input must be polynomials over Q in`,ord) fi; `charsets/premas`(ps,as,ord) fi end: `charsets/charset`:=proc(ps,lst,medset,y) local mset,ord,qs; if nargs < 2 then ERROR(`too few arguments`) elif nops(ps) < 1 then ERROR(`no polynomials specified`) elif nops(lst) < 1 then ERROR(`no indeterminates specified`) elif 4 < nargs then ERROR(`too many arguments`) fi; if member(false,map(type,lst,name)) then ERROR(`bad variable list`) fi; if member(false,map(type,ps,polynom(polynom(rational),lst))) then ERROR(`input must be polynomials over Q in`,lst) fi; qs := {op(`charsets/expand`(ps))} minus {0}; if type(lst,list) then ord := lst else ord := `charsets/reorder`([op(lst)],`charsets/degord`,qs) fi; if nargs < 3 then mset := 'charsetn' else mset := medset fi; if 3 < nargs then y := ord fi; if member(mset,{'charsetn','wcharsetn','qcharsetn','wbasset', 'qbasset','triset','trisetc','basset','autored'}) then `charsets/charseta`(qs,ord,cat(`charsets/`,mset)) else ERROR(`medial set must be one of basset,wbasset, qbasset,charsetn,wcharsetn,qcharsetn,triset,trisetc`); fi; end: `charsets/expand`:=proc(s) local i; if type(s,list) then ['expand(s[i])' $ ('i' = 1 .. nops(s))] elif type(s,set) then {'expand(s[i])' $ ('i' = 1 .. nops(s))} else expand(s) fi end: `charsets/mcharset`:=proc(ps,lst,medset,y) local mset,ord,qs; if nargs < 2 then ERROR(`too few arguments`) elif nops(ps) < 1 then ERROR(`no polynomials specified`) elif nops(lst) < 1 then ERROR(`no indeterminates specified`) elif 4 < nargs then ERROR(`too many arguments`) fi; if member(false,map(type,lst,name)) then ERROR(`bad variable list`) fi; if member(false,map(type,ps,polynom(polynom(rational),lst))) then ERROR(`input must be polynomials over Q in`,lst) fi; qs := {op(`charsets/expand`(ps))} minus {0}; if type(lst,list) then ord := lst else ord := `charsets/reorder`([op(lst)],`charsets/degord`,qs) fi; if 3 < nargs then y := ord fi; if nargs < 3 then mset := 'charsetn' else mset := medset fi; if member(mset,{'charsetn','wcharsetn','qcharsetn','wbasset','qbasset', 'triset','trisetc','basset','autored'}) then `charsets/fcharseta`(qs,ord,cat(`charsets/`,mset)) else ERROR(`medial set must be one of basset,wbasset,qbasset,charsetn,wcharsetn,qcharsetn,triset,trisetc`) fi end: `charsets/iniset`:=proc(as,ord) local ind,i; ind:= 0; if nargs <> 2 then ERROR(`wrong number of arguments`) elif nops(as) < 1 then ERROR(`no polynomials specified`) elif nops(ord) < 1 then ERROR(`no indeterminates specified`) elif not type(ord,list) then ERROR(ord,`must be a list`) fi; if member(false,map(type,ord,name)) then ERROR(`bad variable list`) fi; if member(false,map(type,as,polynom(polynom(rational),ord))) then ERROR(`input must be polynomials over Q in`,ord) fi; for i to nops(as) do if `charsets/class`(as[i],ord) <= ind then ERROR(`first argument must be a non-contradictory (weak-, quasi-)ascending set`) else ind:=`charsets/class`(as[i],ord) fi od; `charsets/initialset`(`charsets/expand`(as),ord) end: `charsets/charser`:= proc(ps,lst,medset,y) local mset,ord,qs; if nargs < 2 then ERROR(`too few arguments`) elif nops(ps) < 1 then ERROR(`no polynomials specified`) elif nops(lst) < 1 then ERROR(`no indeterminates specified`) elif 4 < nargs then ERROR(`too many arguments`) fi; if member(false,map(type,lst,name)) then ERROR(`bad variable list`) fi; if member(false,map(type,ps,polynom(polynom(rational),lst))) then ERROR(`input must be polynomials over Q in`,lst) fi; qs := {op(`charsets/expand`(ps))} minus {0}; if type(lst,list) then ord := lst else ord := `charsets/reorder`([op(lst)],`charsets/degord`,qs) fi; if 3 < nargs then y := ord fi; if nargs < 3 then mset := 'charsetn' else mset := medset fi; if member(mset,{'charsetn','wcharsetn','wbasset','trisetc','basset', 'autored'}) then `charsets/charseries`(qs,ord,cat(`charsets/`,mset)) else ERROR(`medial set must be one of basset,wbasset,charsetn,wcharsetn,trisetc`) fi end: `charsets/mcs`:= proc(ps,lst,medset,y) local mset,ord,qs; if nargs < 2 then ERROR(`too few arguments`) elif nops(ps) < 1 then ERROR(`no polynomials specified`) elif nops(lst) < 1 then ERROR(`no indeterminates specified`) elif 4 < nargs then ERROR(`too many arguments`) fi; if member(false,map(type,lst,name)) then ERROR(`bad variable list`) fi; if member(false,map(type,ps,polynom(polynom(rational),lst))) then ERROR(`input must be polynomials over Q in`,lst) fi; qs := {op(`charsets/expand`(ps))} minus {0}; if type(lst,list) then ord := lst else ord := `charsets/reorder`([op(lst)],`charsets/degord`,qs) fi; if 3 < nargs then y := ord fi; if nargs < 3 then mset := 'charsetn' else mset := medset fi; if member(mset,{'charsetn','wcharsetn','wbasset','trisetc','basset', 'autored'}) then `charsets/fcharser`(qs,ord,cat(`charsets/`,mset)) else ERROR(`medial set must be one of ``basset``,``wbasset``,`. ```charsetn``,``wcharsetn`` and ``trisetc```) fi end: `charsets/ecs`:= proc(ps,lst,medset,y) local mset,ord,qs; if nargs < 2 then ERROR(`too few arguments`) elif nops(ps) < 1 then ERROR(`no polynomials specified`) elif nops(lst) < 1 then ERROR(`no indeterminates specified`) elif 4 < nargs then ERROR(`too many arguments`) fi; if member(false,map(type,lst,name)) then ERROR(`bad variable list`) fi; if type(ps[1],list) then if member(false,map(type,ps[1],polynom(polynom(rational),lst))) then ERROR(`input must be polynomials over Q in`,lst) fi elif member(false,map(type,ps,polynom(polynom(rational),lst))) then ERROR(`input must be polynomials over Q in`,lst) fi; if type(ps[1],{set,list}) then qs := [{op(`charsets/expand`(ps[1]))} minus {0},ps[2]] else qs := {op(`charsets/expand`(ps))} minus {0} fi; if type(lst,list) then ord := lst else if type(ps[1],{set,list}) then ord := `charsets/reorder`([op(lst)],`charsets/degord`,qs[1]) else ord := `charsets/reorder`([op(lst)],`charsets/degord`,qs) fi fi; if 3 < nargs then y := ord fi; if nargs < 3 then mset := 'charsetn' else mset := medset fi; if member(mset,{'charsetn','wcharsetn','wbasset','trisetc','basset', 'autored'}) then `charsets/excharser`(qs,ord,cat(`charsets/`,mset)) else ERROR(`medial set must be one of ``basset``,``wbasset``,`. ```charsetn``,``wcharsetn`` and ``trisetc```) fi end: `charsets/mecs`:= proc(ps,lst,medset,y) local mset,ord,qs; if nargs < 2 then ERROR(`too few arguments`) elif nops(ps) < 1 then ERROR(`no polynomials specified`) elif nops(lst) < 1 then ERROR(`no indeterminates specified`) elif 4 < nargs then ERROR(`too many arguments`) fi; if member(false,map(type,lst,name)) then ERROR(`bad variable list`) fi; if type(ps[1],list) then if member(false,map(type,ps[1],polynom(polynom(rational),lst))) then ERROR(`input must be polynomials over Q in`,lst) fi elif member(false,map(type,ps,polynom(polynom(rational),lst))) then ERROR(`input must be polynomials over Q in`,lst) fi; if type(ps[1],{set,list}) then qs := [{op(`charsets/expand`(ps[1]))} minus {0},ps[2]] else qs := {op(`charsets/expand`(ps))} minus {0} fi; if type(lst,list) then ord := lst else if type(ps[1],{set,list}) then ord := `charsets/reorder`([op(lst)],`charsets/degord`,qs[1]) else ord := `charsets/reorder`([op(lst)],`charsets/degord`,qs) fi fi; if 3 < nargs then y := ord fi; if nargs < 3 then mset := 'charsetn' else mset := medset fi; if member(mset,{'charsetn','wcharsetn','wbasset','trisetc','basset', 'autored'}) then `charsets/fexcharser`(qs,ord,cat(`charsets/`,mset)) else ERROR(`medial set must be one of ``basset``,``wbasset``,`. ```charsetn``,``wcharsetn`` and ``trisetc```) fi end: `charsets/ics`:=proc(ps,lst,medset,y) local mset,ord,qs; if nargs < 2 then ERROR(`too few arguments`) elif nops(ps) < 1 then ERROR(`no polynomials specified`) elif nops(lst) < 1 then ERROR(`no indeterminates specified`) elif 4 < nargs then ERROR(`too many arguments`) fi; if member(false,map(type,lst,name)) then ERROR(`bad variable list`) fi; if member(false,map(type,ps,polynom(polynom(rational),lst))) then ERROR(`input must be polynomials over Q in`,lst) fi; qs := {op(`charsets/expand`(ps))} minus {0}; if type(lst,list) then ord := lst else ord := `charsets/reorder`([op(lst)],`charsets/degord`,qs) fi; if 3 < nargs then y := ord fi; if nargs < 3 then mset := 'charsetn' else mset := medset fi; if member(mset,{'charsetn','trisetc','basset','autored'}) then `charsets/irrcharser`(qs,ord,cat(`charsets/`,mset)) else ERROR( `medial set must be one of ``basset``,``charsetn```.` and ``trisetc``` ) fi end: `charsets/eics` := proc(ps,lst,medset,y) local mset,ord,qs; if nargs < 2 then ERROR(`too few arguments`) elif nops(ps) < 1 then ERROR(`no polynomials specified`) elif nops(lst) < 1 then ERROR(`no indeterminates specified`) elif 4 < nargs then ERROR(`too many arguments`) fi; if member(false,map(type,lst,name)) then ERROR(`bad variable list`) fi; if type(ps[1],{set,list}) then if member(false,map(type,ps[1],polynom(polynom(rational),lst))) then ERROR(`input must be polynomials over Q in`,lst) fi elif member(false,map(type,ps,polynom(polynom(rational),lst))) then ERROR(`input must be polynomials over Q in`,lst) fi; if type(ps[1],{set,list}) then qs := [{op(`charsets/expand`(ps[1]))} minus {0},ps[2]] else qs := {op(`charsets/expand`(ps))} minus {0} fi; if type(lst,list) then ord := lst else if type(ps[1],{set,list}) then ord := `charsets/reorder`([op(lst)],`charsets/degord`,qs[1]) else ord := `charsets/reorder`([op(lst)],`charsets/degord`,qs) fi fi; if 3 < nargs then y := ord fi; if nargs < 3 then mset := 'charsetn' else mset := medset fi; if member(mset,{'charsetn','trisetc','basset','autored'}) then `charsets/exirrcharser`(`charsets/expand`(qs),ord,cat(`charsets/`,mset)) else ERROR( `medial set must be one of ``basset``,`.```charsetn`` and ``trisetc``` ) fi end: `charsets/qics`:=proc(ps,lst,medset,y) local mset,ord,qs; if nargs < 2 then ERROR(`too few arguments`) elif nops(ps) < 1 then ERROR(`no polynomials specified`) elif nops(lst) < 1 then ERROR(`no indeterminates specified`) elif 4 < nargs then ERROR(`too many arguments`) fi; if member(false,map(type,lst,name)) then ERROR(`bad variable list`) fi; if member(false,map(type,ps,polynom(polynom(rational),lst))) then ERROR(`input must be polynomials over Q in`,lst) fi; qs := {op(`charsets/expand`(ps))} minus {0}; if type(lst,list) then ord := lst else ord := `charsets/reorder`([op(lst)],`charsets/degord`,qs) fi; if 3 < nargs then y := ord fi; if nargs < 3 then mset := 'charsetn' else mset := medset fi; if member(mset,{'charsetn','wcharsetn','wbasset','trisetc','basset', 'autored'}) then `charsets/qirrcharser`(qs,ord,cat(`charsets/`,mset)) else ERROR(`medial set must be one of ``basset``,``wbasset``,`. ```charsetn``,``wcharsetn`` and ``trisetc```) fi end: `charsets/cfactor`:= proc(f,as,ord) local ind,inda,ff,i; global `charsets/das`; if nargs = 1 then RETURN(factor(f)) fi; if nargs = 2 then ERROR(`inproper number of arguments`) elif nops(as) < 1 then ERROR(`no polynomials specified`) elif nops(ord) < 1 then ERROR(`no indeterminates specified`) elif not type(ord,list) then ERROR(ord,`must be a list`) elif 4 < nargs then ERROR(`too many arguments`) fi; if member(false,map(type,ord,name)) then ERROR(`bad variable list`) fi; if member(false,map(type,as,polynom(polynom(rational),ord))) then ERROR(`input must be polynomials over Q in`,ord) fi; ff := numer(f); ind := 0; for i to nops(as) do inda := `charsets/class`(as[i],ord); if inda <= ind then ERROR(`second argument must be a non-contradictory ascending set`) else ind := inda fi od; if 1 < printlevel then lprint(`Warning: be sure the ascending set is irreducible`) fi; if `charsets/class`(ff,ord) <= `charsets/class`(as[nops(as)],ord) then factor(f) else sum('`charsets/degreel`(as[i],ord)','i' = 1 .. nops(as)); if `charsets/degreel`(ff,ord) < % then `charsets/das` := [-1,1,-2,2,-3,false] else `charsets/das` := [1,-1,2,-2,-3,false] fi; `charsets/cfactorsub`(factor(f),as,ord) fi end: `charsets/triser`:=proc(ps,lst,y) local i,ord,qs; if nargs < 1 then ERROR(`too few arguments`) elif nops(ps) < 1 then ERROR(`no polynomials specified`) elif 3 < nargs then ERROR(`too many arguments`) elif nargs = 2 then if nops(lst) < 1 then ERROR(`no indeterminates specified`) fi fi; if nargs = 2 then if member(false,map(type,lst,name)) then ERROR(`bad variable list`) fi fi; if member(false,map(type,ps,polynom(polynom(rational),lst))) then ERROR(`input must be polynomials over Q in`,lst) fi; qs := {op(`charsets/expand`(ps))} minus {0}; if nargs < 2 then ord := `charsets/reorder`( [seq(op(indets(ps[i])),i = 1 .. nops(ps))],`charsets/degord`,qs) elif type(lst,list) then ord := lst else ord := `charsets/reorder`([op(lst)],`charsets/degord`,qs) fi; if 2 < nargs then y := ord fi; `charsets/trisersub`(qs,ord) end: `charsets/csolve`:=proc(ps,lst,y) local i,ord,qsi,qs; if nargs < 1 then ERROR(`too few arguments`) elif nops(ps) < 1 then ERROR(`no polynomials specified`) elif 3 < nargs then ERROR(`too many arguments`) elif nargs = 2 then if nops(lst) < 1 then ERROR(`no indeterminates specified`) fi fi; if nargs = 2 then if member(false,map(type,lst,name)) then ERROR(`bad variable list`) fi fi; if member(false,map(type,ps,polynom(polynom(rational),lst))) then ERROR(`input must be polynomials over Q in`,lst) fi; qs := {op(`charsets/expand`(ps))} minus {0}; if nargs < 2 then ord := `charsets/reorder`( [seq(op(indets(ps[i])),i = 1 .. nops(ps))],`charsets/degord`,qs) elif type(lst,list) then ord := lst else ord := `charsets/reorder`([op(lst)],`charsets/degord`,qs) fi; if 2 < nargs then y := ord fi; qsi := {`charsets/trisersub`(qs,ord)}; if qsi = {{}} then {} else op({seq(`charsets/solveas`(qsi[i],ord),i = 1 .. nops(qsi))}) fi end: `charsets/ivd` := proc(ps,lst,medset,y) local mset,ord,qs; if nargs < 2 then ERROR(`too few arguments`) elif nops(ps) < 1 then ERROR(`no polynomials specified`) elif nops(lst) < 1 then ERROR(`no indeterminates specified`) elif 4 < nargs then ERROR(`too many arguments`) fi; if member(false,map(type,lst,name)) then ERROR(`bad variable list`) fi; if member(false,map(type,ps,polynom(polynom(rational),lst))) then ERROR(`input must be polynomials over Q in`,lst) fi; qs := {op(`charsets/expand`(ps))} minus {0}; if type(lst,list) then ord := lst else ord := `charsets/reorder`([op(lst)],`charsets/degord`,qs) fi; if 3 < nargs then y := ord fi; if nargs < 3 then mset := 'charsetn' else mset := medset fi; if member(mset,{'charsetn','trisetc','basset','autored'}) then `charsets/irrvardec`(qs,ord,cat(`charsets/`,mset)) else ERROR( `medial set must be one of ``basset``,`.```charsetn`` and ``trisetc``` ) fi end: `charsets/pid` :=proc(ps,lst,medset,y) local mset,ord,qs; if nargs < 2 then ERROR(`too few arguments`) elif nops(ps) < 1 then ERROR(`no polynomials specified`) elif nops(lst) < 1 then ERROR(`no indeterminates specified`) elif 4 < nargs then ERROR(`too many arguments`) fi; if member(false,map(type,lst,name)) then ERROR(`bad variable list`) fi; if member(false,map(type,ps,polynom(polynom(rational),lst))) then ERROR(`input must be polynomials over Q in`,lst) fi; qs := {op(`charsets/expand`(ps))} minus {0}; if type(lst,list) then ord := lst else ord := `charsets/reorder`([op(lst)],`charsets/degord`,qs) fi; if 3 < nargs then y := ord fi; if nargs < 3 then mset := 'charsetn' else mset := medset fi; if member(mset,{'charsetn','trisetc','basset','autored'}) then `charsets/pridealdec`(qs,ord,cat(`charsets/`,mset)) else ERROR( `medial set must be one of ``basset``,`.```charsetn`` and ``trisetc``` ) fi end: `charsets/class` := proc(f,ord) local i; options remember,system; for i from nops(ord) by -1 to 1 do if member(ord[i],indets(f)) then RETURN(i) fi od; 0 end: `charsets/lvar`:=proc(f,ord) local i; options remember,system; for i from nops(ord) by -1 to 1 do if member(ord[i],indets(f)) then RETURN(ord[i]) fi od; if 1 < printlevel then lprint(`Warning: lvar is called with constant`) fi; 0 end: `charsets/index`:=proc(f,ord) local i,g,ng; if type(f,list) then ['`charsets/index`(f[i],ord)' $ ('i' = 1 .. nops(f))] elif type(f,set) then {'`charsets/index`(f[i],ord)' $ ('i' = 1 .. nops(f))} else g := expand(f); ng := `charsets/terms`(g); if `charsets/class`(g,ord) = 0 then i := op(1,[op(indets(g)),O]) else i := `charsets/lvar`(g,ord) fi; g := degree(g,i); cat(cat(cat(cat(cat(cat(`[`,ng),` `),i),` `),g),`]`) fi end: `charsets/terms` := proc(g) if type(g,`+`) then nops(g) else 1 fi end: `charsets/initial`:=proc(p,ord) local f; options remember,system; f := expand(p); if `charsets/class`(f,ord) = 0 then 1 else lcoeff(f,`charsets/lvar`(f,ord)); numer(%/lcoeff(%)) fi end: `charsets/mrank`:=proc(f,g,ord) local cf,cg; options remember,system; cf := `charsets/class`(f,ord); cg := `charsets/class`(g,ord); if cf = 0 then true elif cf < cg then true elif cf = cg then cf := `charsets/degreel`(f,ord); cg := `charsets/degreel`(g,ord); if cf < cg then true elif cf = cg then `charsets/mrank`( `charsets/initial`(f,ord),`charsets/initial`(g,ord),ord) else false fi else false fi end: `charsets/rank`:=proc(f,g,ord) local ind,find,cf,cg; options remember,system; find := `charsets/subrank`(f,g,ord,'ind'); if find and ind = 1 then cf := nops(expand(`charsets/initial`(f,ord))); cg := nops(expand(`charsets/initial`(g,ord))); if cf < cg then true elif cf = cg then if nops(expand(f)) < nops(expand(g)) then true else false fi else false fi else find fi end: `charsets/subrank`:=proc(f,g,ord,ind) local cf,cg; options remember,system; cf := `charsets/class`(f,ord); cg := `charsets/class`(g,ord); if cf = 0 then if cg = 0 then ind := 1 fi; true elif cf < cg then true elif cf = cg then cf := `charsets/degreel`(f,ord); cg := `charsets/degreel`(g,ord); if cf < cg then true elif cf = cg then `charsets/subrank`( `charsets/initial`(f,ord),`charsets/initial`(g,ord),ord,'ind') else false fi else false fi end: `charsets/trank`:=proc(f,g,ord) local cf,cg; options remember,system; cf := `charsets/degreel`(f,ord); cg := `charsets/degreel`(g,ord); if cf < cg then true elif cf = cg then `charsets/mrank`( `charsets/initial`(f,ord),`charsets/initial`(g,ord),ord) else false fi end: `charsets/prem`:=proc(uu,vv,x) local r,v,dr,dv,l,t,lu,lv; options remember,system; if type(vv/x,integer) then subs(x = 0,uu) else r := expand(uu); dr := degree(r,x); v := expand(vv); dv := degree(v,x); if dv <= dr then l := coeff(v,x,dv); v := expand(v-l*x^dv) else l := 1 fi; if 1 < printlevel and 500 < nops(r)+nops(v) then lprint(`pseudo-division:`,[nops(r),x,dr],[nops(v),x,dv]) fi; while dv <= dr and r <> 0 do gcd(l,coeff(r,x,dr),'lu','lv'); t := expand(x^(dr-dv)*v*lv); if dr = 0 then r := 0 else r := subs(x^dr = 0,r) fi; r := expand(lu*r)-t; dr := degree(r,x) od; r fi end: `charsets/premas`:=proc(f,as,ord) local remd,i; remd := f; for i from nops(as) by -1 to 1 do remd := `charsets/prem`(remd,as[i],`charsets/lvar`(as[i],ord)) od; if remd <> 0 and type(remd,polynom) then numer(remd/lcoeff(remd)) else 0 fi end: `charsets/remseta`:=proc(ps,as,ord) local i; {'`charsets/premas`(ps[i],as,ord)' $ ('i' = 1 .. nops(ps))} minus {0} end: `charsets/premasb`:=proc(f,as,ord) local remd,i; remd := f; if 1 < nops(as) then for i from nops(as) by -1 to 2 do remd := `charsets/prem`(remd,as[i],`charsets/lvar`(as[i],ord)) od fi; if divide(remd,primpart(as[1],`charsets/lvar`(as[1],ord))) then remd := 0 else remd := `charsets/prem`(remd,as[1],`charsets/lvar`(as[1],ord)) fi; if remd <> 0 then numer(remd/lcoeff(remd)) else 0 fi end: `charsets/remsetb`:=proc(ps,as,ord) local i; {'`charsets/premasb`(ps[i],as,ord)' $ ('i' = 1 .. nops(ps))} minus {0} end: `charsets/reorder`:=proc(ord,p,ps) op(`charsets/reordera`(ord,ps)); [op(`charsets/reorderb`([op({op(ord)} minus {%})],p,ps)),%] end: `charsets/reordera`:=proc(ord,ps) local qs,pp,orb,i; if nops(ps) = 0 then ord else qs := {op(ps)}; orb := {op(ord)}; for i in orb do pp := `charsets/deg0`(ps,i); if nops(pp) = 1 then RETURN( [op(`charsets/reordera`([op(orb minus {i})],qs minus pp)),i] ) fi od; [] fi end: `charsets/reorderb`:=proc(l,p,ps) local n,tn,gap,i,j,temp,v; n := nops(l); tn := p; for i to n do v[i-1] := l[i] od; for gap from 4 while gap <= n do gap := 3*gap+1 od; gap := iquo(gap,3); while 0 < gap do for i from gap to n-1 do temp := v[i]; for j from i-gap by -gap to 0 do if tn(v[j],temp,ps) then break fi; v[j+gap] := v[j] od; v[j+gap] := temp od; gap := iquo(gap,3) od; ['eval(v[i],1)' $ ('i' = 0 .. n-1)] end: `charsets/degord` :=proc(x,y,ps) if op(2,`charsets/degpsmax`(ps,y)) < op(2,`charsets/degpsmax`(ps,x)) then true elif op(2,`charsets/degpsmax`(ps,x)) < op(2,`charsets/degpsmax`(ps,y)) then false elif op(1,`charsets/degpsmax`(ps,y)) < op(1,`charsets/degpsmax`(ps,x)) then true elif op(1,`charsets/degpsmax`(ps,x)) < op(1,`charsets/degpsmax`(ps,y)) then false elif op(2,`charsets/degpsmin`(ps,x)) < op(2,`charsets/degpsmin`(ps,y)) then true elif op(2,`charsets/degpsmin`(ps,y)) < op(2,`charsets/degpsmin`(ps,x)) then false elif op(1,`charsets/degpsmin`(ps,y)) < op(1,`charsets/degpsmin`(ps,x)) then true elif op(1,`charsets/degpsmin`(ps,x)) < op(1,`charsets/degpsmin`(ps,y)) then false elif op(1,`charsets/deg1`(ps,y)) < op(1,`charsets/deg1`(ps,x)) then true elif op(1,`charsets/deg1`(ps,x)) < op(1,`charsets/deg1`(ps,y)) then false elif op(2,`charsets/deg1`(ps,y)) < op(2,`charsets/deg1`(ps,x)) then true else false fi end: `charsets/degpsmax` :=proc(qs,x) local i,m,mm,ps; options remember,system; ps := expand(qs); m := max('degree(ps[i],x)' $ ('i' = 1 .. nops(ps))); mm := 0; for i to nops(ps) do if degree(ps[i],x) = m then mm := mm+m fi od; [mm,m] end: `charsets/degpsmin`:=proc(qs,x) local i,m,mm,ps; options remember,system; ps := expand(qs); {'degree(ps[i],x)' $ ('i' = 1 .. nops(ps))} minus {0}; if % = {} then m := 0 else m := min(op(%)) fi; mm := 0; for i to nops(ps) do if degree(ps[i],x) = m then mm := mm+m fi od; [mm,m] end: `charsets/deg0` := proc(ps,x) local i,ms; ms := {}; for i in ps while nops(ms) < 2 do if member(x,indets(i)) then ms := {i,op(ms)} fi od; ms end: `charsets/deg1`:=proc(PS,x) local i,ps,qs,k; options remember,system; ps := expand(PS); qs := {}; k := op(2,`charsets/degpsmin`(ps,x)); for i to nops(ps) do if degree(ps[i],x) = k then qs := {op(qs),lcoeff(ps[i],x)} fi od; [min('degree(qs[i],indets(qs[i]))' $ ('i' = 1 .. nops(qs))), min('nops(expand(qs[i]))' $ ('i' = 1 .. nops(qs)))] end: `charsets/sort`:=proc(ps,rank,ord,qs) local l,i,qs1; if nops(ps) = 1 then qs := []; ps[1] else l := ps[1]; qs1 := []; for i from 2 to nops(ps) do if rank(ps[i],l,ord) then qs1 := [l,op(qs1)]; l := ps[i] else qs1 := [ps[i],op(qs1)] fi od; qs := qs1; l fi end: `charsets/minus` := proc(ps,qs) [op({op(ps)} minus {op(qs)})] end: `charsets/union` :=proc(ps1,ps2,ps3) [op(({op(ps1)} union {op(ps2)}) union {op(ps3)})] end: `charsets/prod`:=proc(ps) local i; product('ps[i]','i' = 1 .. nops(ps)) end: `charsets/degree` := proc(f,x) degree(collect(f,x,normal),x) end: `charsets/degreel`:=proc(f,ord) expand(f); degree(%,`charsets/lvar`(%,ord)) end: `charsets/basset`:=proc(ps,ord) local qs,qs1,i,b; if nops(ps) < 2 then ps else b := `charsets/sort`(ps,`charsets/rank`,ord,'qs1'); qs := []; if 0 < `charsets/class`(b,ord) then for i in qs1 do if `charsets/brank`(i,b,ord) then qs := [i,op(qs)] fi od else RETURN([b]) fi; [b,op(`charsets/basset`(qs,ord))] fi end: `charsets/wbasset`:= subs(`charsets/brank`=`charsets/wbrank`, `charsets/basset`=`charsets/wbasset`, op(`charsets/basset`)): `charsets/qbasset` := subs(`charsets/brank`=`charsets/qbrank`, `charsets/basset`=`charsets/qbasset`, op(`charsets/basset`)): `charsets/brank`:=proc(a,b,ord) if `charsets/degree`(a,`charsets/lvar`(b,ord)) < `charsets/degreel`(b,ord) then true else false fi end: `charsets/wbrank`:=proc(a,b,ord) if `charsets/class`(b,ord) < `charsets/class`(a,ord) and `charsets/brank`(`charsets/initial`(a,ord),b,ord) then true else false fi end: `charsets/qbrank`:=proc(a,b,ord) if `charsets/class`(b,ord) < `charsets/class`(a,ord) then true else false fi end: `charsets/charseta`:= proc(PS,ord,medset) local ps,cs,rs,l,med; global `charsets/with`; if nops(PS) < 2 then [op(PS)] else ps := `charsets/unigcd`(PS,ord); if medset = `charsets/qcharsetn` then `charsets/with` := {}; med := subs(`charsets/remseta` = `charsets/remsetaA`, `charsets/qcharsetn` = med,op(`charsets/qcharsetn`)); cs := med(ps,ord) else cs := medset(ps,ord) fi; if 0 < `charsets/class`(cs[1],ord) then if member( medset,{`charsets/wbasset`,`charsets/qbasset`,`charsets/basset`} ) then rs := `charsets/remseta`({op(ps)} minus {op(cs)},cs,ord) elif medset = `charsets/qcharsetn` and `charsets/checkwith`( `charsets/with`,`charsets/initialset1`(cs,ord)) then if 1 < printlevel then lprint(`Strategy for 0 remainder succeed`) fi; RETURN(cs) else if medset = `charsets/qcharsetn` and 1 < printlevel then lprint(`Strategy for 0 remainder failed`) fi; rs := `charsets/remsetb`({op(ps)} minus {op(cs)},cs,ord) fi else RETURN([1]) fi; if rs = {} then ['numer(cs[l]/lcoeff(expand(cs[l])))' $ ('l' = 1 .. nops(cs))] elif medset = `charsets/autored` then `charsets/charseta`(`charsets/union`(rs,cs,ps),ord, `charsets/charsetn`) else `charsets/charseta`(`charsets/union`(rs,cs,ps),ord,medset) fi fi end: `charsets/unigcd` := proc(PS,ord) local ps,p,c; ps := {}; for p in PS do c := `charsets/class`(p,ord); if c = 0 then RETURN({1}) elif c = 1 then ps := ps union {p} fi od; if nops(ps) < 2 then PS else {op(PS)} minus ps union {`charsets/GCD`(op(ps))} fi end: `charsets/GCD` := proc() local i; if nargs = 2 then gcd(args) elif nargs > 2 then gcd(args[1],`charsets/GCD`('args'[i]$'i'=2..nargs)) else ERROR(`invalid number of arguments`) fi end: `charsets/fcharseta`:=proc(ps,ord,medset) local csf,fset; csf := `charsets/fcharsetsub`( `charsets/nopower`(ps),ord,medset,[{},indets(ps)],'fset'); csf,`factors removed` = fset[1] end: `charsets/fcharsetsub` :=proc(PS,ord,medset,fset1,fset) local ps,cs,rs,l,fset3,fset2,ts,fmedset,med; global `charsets/with`; if nops(PS) < 2 then fset := fset1; [op(PS)] else ps := `charsets/unigcd`(PS,ord); if member(substring(medset,10 .. length(medset)), {'charsetn','wcharsetn','qcharsetn','triset', 'trisetc','autored'}) then fmedset := cat(`charsets/f`,(substring(medset,10 .. length(medset)))); if fmedset = `charsets/fqcharsetn` then `charsets/with` := {}; med := subs(`charsets/remseta` = `charsets/remsetaA`, `charsets/fqcharsetn` = med,op(`charsets/fqcharsetn`)); cs := med(ps,ord,fset1,'fset2') else cs := fmedset(ps,ord,fset1,'fset2') fi; if 0 < `charsets/class`(cs[1],ord) then if fmedset = `charsets/fqcharsetn` and `charsets/checkwith`( `charsets/with`,`charsets/initialset1`(cs,ord) union fset2[1] ) then if 1 < printlevel then lprint(`Strategy for 0 remainder succeed`) fi; fset := fset2; RETURN(cs) else if medset = `charsets/qcharsetn` and 1 < printlevel then lprint(`Strategy for 0 remainder failed`) fi; rs := `charsets/remsetb`({op(ps)} minus {op(cs)},cs,ord) fi else fset := fset2; RETURN([1]) fi; if rs = {} then fset := fset2; ['numer(cs[l]/lcoeff(expand(cs[l])))' $ ('l' = 1 .. nops(cs))] else `charsets/fcharsetsub`( `charsets/union`(rs,cs,ps),ord,medset,fset2,'fset') fi else cs := medset(ps,ord); fset2 := [fset1[1],fset1[2] union `charsets/initialset1`(cs,ord)]; if 0 < `charsets/class`(cs[1],ord) then rs := `charsets/remseta`({op(ps)} minus {op(cs)},cs,ord); if nops(rs) <> 0 then rs := `charsets/removecont`(rs,ord,'ts'); rs := `charsets/removefactor`(rs,ord,fset2,'fset3'); fset3 := [fset3[1] union ts,fset3[2]] else rs := []; fset3 := fset2 fi else fset := fset2; RETURN([1]) fi; if rs = [] then fset := fset3; ['numer(cs[l]/lcoeff(expand(cs[l])))' $ ('l' = 1 .. nops(cs))] elif medset = `charsets/autored` then `charsets/fcharsetsub`(`charsets/union`(rs,cs,ps), ord,`charsets/charsetn`,fset3,'fset') else `charsets/fcharsetsub`( `charsets/union`(rs,cs,ps),ord,medset,fset3,'fset') fi fi fi end: `charsets/premA`:=proc(uu,vv,x) local r,v,dr,dv,l,t,lu,lv; global `charsets/with`; if type(vv/x,integer) then subs(x = 0,uu) else r := expand(uu); dr := degree(r,x); v := expand(vv); dv := degree(v,x); if dv <= dr then l := coeff(v,x,dv); v := expand(v-l*x^dv) else l := 1 fi; if 1 < printlevel and 500 < nops(r)+nops(v) then lprint(`pseudo-division:`,[nops(r),x,dr],[nops(v),x,dv]) fi; while dv <= dr and r <> 0 do gcd(l,coeff(r,x,dr),'lu','lv'); t := expand(x^(dr-dv)*v*lv); if dr = 0 then r := 0 else r := subs(x^dr = 0,r) fi; r := expand(lu*r)-t; if not type(lu,rational) and type(`charsets/with`,set) then `charsets/with` := `charsets/with` union {lu} fi; dr := degree(r,x) od; r fi end: `charsets/remsetaA`:=proc(ps,as,ord) local i; {'`charsets/premasA`(ps[i],as,ord)' $ ('i' = 1 .. nops(ps))} minus {0} end: `charsets/premasA`:=proc(f,as,ord) local remd,i; remd := f; for i from nops(as) by -1 to 1 do remd := `charsets/premA`(remd,as[i],`charsets/lvar`(as[i],ord)) od; if remd <> 0 then numer(remd/lcoeff(remd)) else 0 fi end: `charsets/checkwith`:=proc(ps1,ps2) local rs,i,j,r; rs := ps1 minus ps2; if rs = {} then true elif ps2 = {} then false else rs := {'`charsets/pfactor`(convert(rs[i],sqrfree))' $ ('i' = 1 .. nops(rs))} ; for i to nops(rs) do r := rs[i]; for j to nops(ps2) do gcd(r,ps2[j],'r'); if type(r,rational) then break fi od; if not type(r,rational) then RETURN(false) fi od; true fi end: `charsets/nopower`:=proc(as) local i; if not type(as,{set,list}) then if type(as,`^`) then op(1,as) elif type(as,`*`) then product('`charsets/nopower`(op(i,as))','i' = 1 .. nops(as)) else as fi else ['`charsets/nopower`(as[i])' $ ('i' = 1 .. nops(as))] fi end: `charsets/charsetn`:=proc(PS,ord) local ps,cs,rs; if nops(PS) < 2 then [op(PS)] else ps := `charsets/unigcd`(PS,ord); cs := `charsets/basset`(ps,ord); if 0 < `charsets/class`(cs[1],ord) then rs := `charsets/remseta`(`charsets/minus`(ps,cs),cs,ord) else RETURN(cs) fi; if rs = {} then cs else `charsets/charsetn`([op(rs),op(cs)],ord) fi fi end: `charsets/wcharsetn`:= subs(`charsets/basset`=`charsets/wbasset`, `charsets/charsetn`=`charsets/wcharsetn`, op(`charsets/charsetn`)): `charsets/qcharsetn`:= subs(`charsets/basset`=`charsets/qbasset`, `charsets/charsetn`=`charsets/qcharsetn`, op(`charsets/charsetn`)): `charsets/fcharsetn`:=proc(PS,ord,fset1,fset) local ps,cs,rs,ts,fset2,fset3; if nops(PS) < 2 then fset := fset1; [op(PS)] else ps := `charsets/unigcd`(PS,ord); cs := `charsets/basset`(ps,ord); fset2 := [fset1[1],fset1[2] union `charsets/initialset1`(cs,ord)]; if 0 < `charsets/class`(cs[1],ord) then rs := `charsets/remseta`(`charsets/minus`(ps,cs),cs,ord); if nops(rs) <> 0 then rs := `charsets/removecont`(rs,ord,'ts'); rs := `charsets/removefactor`(rs,ord,fset2,'fset3'); fset3 := [fset3[1] union ts,fset3[2]] else rs := []; fset3 := fset2 fi else fset := fset2; RETURN(cs) fi; if rs = [] then fset := fset3; cs else `charsets/fcharsetn`([op(rs),op(cs)],ord,fset3,'fset') fi fi end: `charsets/fwcharsetn`:=subs(`charsets/basset`=`charsets/wbasset`, `charsets/fcharsetn`=`charsets/fwcharsetn`, op(`charsets/fcharsetn`)): `charsets/fqcharsetn`:=subs(`charsets/basset`=`charsets/qbasset`, `charsets/fcharsetn`=`charsets/fqcharsetn`, op(`charsets/fcharsetn`)): `charsets/triset`:=proc(PS,ord) local ps,i; global `charsets/@qs`; if nops(PS) < 2 then [op(PS)] else ps := `charsets/unigcd`(PS,ord); for i from 0 to nops(ord) do `charsets/@qs`[i] := [] od; for i in ps do `charsets/@qs`[`charsets/class`(i,ord)] := [op(`charsets/@qs`[`charsets/class`(i,ord)]),i] od; for i from nops(ord) by -1 to 1 do if `charsets/@qs`[0] = [] then `charsets/subtriset`(i,ord) else RETURN([1]) fi od; if `charsets/@qs`[0] <> [] then [1] else ['op(`charsets/@qs`[i])' $ ('i' = 1 .. nops(ord))] fi fi end: `charsets/subtriset`:=proc(i,ord) local ss,ss1,j,p; global `charsets/@qs`; if 1 < nops(`charsets/@qs`[i]) then ss1 := `charsets/sort`(`charsets/@qs`[i],`charsets/trank`,ord,'ss'); `charsets/@qs`[i] := [ss1]; for j in ss do p := `charsets/prem`(j,ss1,ord[i]); if p <> 0 then `charsets/@qs`[`charsets/class`(p,ord)] := [ op(`charsets/@qs`[`charsets/class`(p,ord)]),numer(p/lcoeff(p)) ] fi od; `charsets/subtriset`(i,ord) fi end: `charsets/movefactor`:=proc(f,g,ord,ja) local fg; if not type(f,{set,list}) and not type(g,rational) and divide(f,g,'fg') then if 3 < nargs then if 0 < `charsets/class`(g,ord) then ja := true else ja := false fi fi; `charsets/movefactor`(fg,g,ord) else if 3 < nargs then ja := false fi; f fi end: `charsets/removefactor`:=proc(ps,ord,fset1,fset) local k,rr,ja,fset2,fs,qs,rs,r; if not type(ps,{set,list}) then qs := {ps} else qs := ps fi; rs := {}; fs := fset1; for r in qs do rr := r; fset2 := {}; for k in fs[1] do rr := `charsets/movefactor`(rr,k,ord,'ja') od; for k in fs[2] do if rr <> k then rr := `charsets/movefactor`(rr,k,ord,'ja'); if ja then fset2 := {numer(k/lcoeff(expand(k))),op(fset2)} fi fi od; fs := [fs[1] union fset2,fs[2] minus fset2]; rs := {rr,op(rs)} od; rs := `charsets/nopower`(rs); fset := fs; if not type(ps,{set,list}) then rs[1] else rs fi end: `charsets/removecont`:=proc(ps,ord,ms) local qs,fs,i,cc,pp; if type(ps,{set,list}) then if `charsets/class`(ps[1],ord) = 0 then if 2 < nargs then ms := {} fi; ps else qs := []; fs := {}; for i to nops(ps) do if `charsets/class`(ps[i],ord)>0 then cc := `charsets/nopower`(content(ps[i], `charsets/lvar`(ps[i],ord),'pp')); if 0 < `charsets/class`(cc,ord) then fs := {cc,op(fs)} fi; qs := [op(qs),pp] else qs := [op(qs),ps[i]] fi; od; if 2 < nargs then ms := fs fi; qs fi else qs := `charsets/removecont`([ps],ord,'fs'); ms := fs; qs[1] fi end: `charsets/ftriset`:=proc(ps,ord,fset1,fset) local i,fset2,var; global `charsets/@fact`,`charsets/@qs`; `charsets/@fact` := 1; if nops(ps) < 2 then fset := fset1; ps else fset2 := {}; for i from 0 to nops(ord) do `charsets/@qs`[i] := [] od; for i in ps do `charsets/@qs`[`charsets/class`(i,ord)] := [op(`charsets/@qs`[`charsets/class`(i,ord)]),i] od; var := indets(ps); for i from nops(ord) by -1 to 1 do if `charsets/@qs`[0] = [] then fset2 := fset2 union `charsets/fsubtriset`(i,ord,ps,var) else fset := [fset2,{}]; RETURN([1]) fi od; if `charsets/@qs`[0] <> [] then fset := [fset2,{}]; [1] else fset := [fset2,{}]; ['op(`charsets/@qs`[i])' $ ('i' = 1 .. nops(ord))] fi fi end: `charsets/fsubtriset`:=proc(i,ord,ps,var) local fset2,ss,ts,ss1,j,k,l,p,pp,qq,ja; global `charsets/@qs`,`charsets/@fact`; if 1 < nops(`charsets/@qs`[i]) then ss1 := `charsets/sort`(`charsets/@qs`[i],`charsets/trank`,ord,'ss'); `charsets/@qs`[i] := [ss1]; fset2 := {}; qq := {'numer(ps[l]/lcoeff(expand(ps[l])))' $ ('l' = 1 .. nops(ps))}; for j in ss do pp := `charsets/prem`(j,ss1,ord[i]); if pp <> 0 then if pp <> `charsets/@fact` and `charsets/fsubtrisetsub`(var,`charsets/@fact`) and not member(`charsets/@fact`,qq) then p := `charsets/removecont`(pp,ord,'ts'); fset2 := {op(fset2),op(ts)}; p := `charsets/movefactor`(p,`charsets/@fact`,ord,'ja'); if ja then fset2 := {op(fset2),numer( `charsets/@fact`/lcoeff(expand(`charsets/@fact`)))} fi else p := pp fi; for k in var do if p <> k and not member(k,qq) then p := `charsets/movefactor`(p,k,ord,'ja'); if ja then fset2 := {op(fset2),k} fi fi od; `charsets/@qs`[`charsets/class`(p,ord)] := [ op(`charsets/@qs`[`charsets/class`(p,ord)]), `charsets/nopower`(numer(p/lcoeff(p)))] fi od; `charsets/initial`(ss1,ord); `charsets/@fact` := numer(%/lcoeff(expand(%))); fset2 union `charsets/fsubtriset`(i,ord,ps,var) else {} fi end: `charsets/fsubtrisetsub`:=proc(aa,bb) local i; for i to nops(aa) do if subs(aa[i] = 0,bb) = 0 then RETURN(false) fi od; true end: `charsets/trisetc`:=proc(ps,ord) local ind,cs; global `charsets/@cs`; `charsets/@cs` := `charsets/triset`(ps,ord); cs := `charsets/subtrisetc`(ord,{},'ind'); if ind = 0 then `charsets/charsetn`([op(ps),op(cs)],ord) else cs fi end: `charsets/ftrisetc`:=proc(ps,ord,fset1,fset) local i,ind,cs,fs; global `charsets/@cs`; `charsets/@cs` := `charsets/ftriset`(ps,ord,fset1,'fs'); fset := fs; if fs[1] <> {} then {'op(`charsets/pfactor`(fs[1][i]))' $ ('i' = 1 .. nops(fs[1]))} else {} fi; cs := `charsets/subtrisetc`(ord,%,'ind'); if ind = 0 then `charsets/fcharsetn`([op(ps),op(cs)],ord,fset1,'fset') else cs fi end: `charsets/subtrisetc`:=proc(ord,var,ind) local r,i,j,cs; global `charsets/@cs`; if nops(`charsets/@cs`) = 0 then cs := [] else cs := [`charsets/@cs`[1]] fi; if 1 < nops(`charsets/@cs`) then for i from 2 to nops(`charsets/@cs`) do r := `charsets/premas`(`charsets/@cs`[i],cs,ord); if `charsets/class`(r,ord) <> `charsets/class`(`charsets/@cs`[i],ord) then ind := 0; if r <> 0 then cs := [op(cs),`charsets/nopower`(numer(r/lcoeff(r)))] fi; break else for j in var do r := `charsets/movefactor`(r,j,ord) od; cs := [op(cs),`charsets/nopower`(numer(r/lcoeff(r)))]; if `charsets/class`(r,ord) <> `charsets/class`(`charsets/@cs`[i],ord) then ind := 0; break fi fi od fi; cs end: `charsets/initialset1`:=proc(as,ord) local i,is,iss; is := {}; for i in as do `charsets/initial`(i,ord); if 0 < `charsets/class`(%,ord) then is := {`charsets/pfactor`(%),op(is)} fi od; is := `charsets/compress`(is,ord); iss := {}; for i in is do if 0 < `charsets/class`(i,ord) then iss := {i,op(iss)} fi od; iss end: `charsets/compress`:=proc(ps,ord) local is,is1,i,j,ss; is := ps; if 1 < nops(is) then is1 := []; for i to nops(is)-1 do ss := is[i]; for j from i+1 to nops(is) do ss := `charsets/movefactor`(ss,is[j],ord) od; if 0 < `charsets/class`(ss,ord) then is1 := [`charsets/pfactor`(ss),op(is1)] fi od; is1 := [is[nops(is)],op(is1)]; is := {}; if 1 < nops(is1) then for i to nops(is1)-1 do ss := is1[i]; for j from i+1 to nops(is1) do ss := `charsets/movefactor`(ss,is1[j],ord) od; if 0 < `charsets/class`(ss,ord) then is := {op(is),`charsets/pfactor`(ss)} fi od; {op(is),is1[nops(is1)]} else {op(is1)} fi else {op(is)} fi end: `charsets/pfactor`:=proc(f) local i; if type(f,integer) then op({}) elif type(f,`^`) then op(1,f); numer(%/lcoeff(%)) elif type(f,`*`) then '`charsets/pfactor`(op(i,f))' $ ('i' = 1 .. nops(f)) else numer(f/lcoeff(expand(f))) fi end: `charsets/sfactor`:=proc(f) local i; if type(f,integer) then op({}) elif type(f,`^`) then op(1,f); 'numer(%/lcoeff(%))' $ ('i' = 1 .. op(2,f)) elif type(f,`*`) then '`charsets/sfactor`(op(i,f))' $ ('i' = 1 .. nops(f)) else numer(f/lcoeff(f)) fi end: `charsets/initialset`:=proc(as,ord) local i,is,iss; is := {}; for i in as do `charsets/initial`(i,ord); if 0 < `charsets/class`(%,ord) then is := {%,op(is)} fi od; is := `charsets/factorps`(is); iss := {}; for i in is do if 0 < `charsets/class`(i,ord) then iss := {i,op(iss)} fi od; iss end: `charsets/factorps`:=proc(ps) local qs,i,j,q; qs := {}; for i in ps do q := factor(i); if type(q,`*`) then for j to nops(q) do if not type(op(j,q),integer) then if type(op(j,q),`^`) then qs := {op(qs),numer(op(1,op(j,q))/lcoeff(op(1,op(j,q))))} else qs := {op(qs),numer(op(j,q)/lcoeff(op(j,q)))} fi fi od elif type(q,`^`) then qs := {op(qs),numer(op(1,q)/lcoeff(op(1,q)))} else if not type(q,integer) then qs := {op(qs),numer(q/lcoeff(q))} fi fi od; qs end: `charsets/charseries`:=proc(ps,ord,medset) local qs,cs,iss,n,qsi,qhi,csno,ppi,qqi; if type(ps[1],{set,list}) then qhi := {op(ps)} else qhi := {{op(ps)}} fi; qsi := {}; csno := 0; ppi := {}; qqi := {}; for n from 0 while qhi <> {} do qhi := sort([op(qhi)],`charsets/nopsord`); qs := qhi[1]; ppi := `charsets/select`(ppi,nops(qs)); qqi := {op(ppi[2]),op(qqi)}; if n = 0 then ppi := {} else ppi := {qs,op(ppi[1])} fi; cs := `charsets/charseta`([op(qs)],ord,medset); if 1 < printlevel then csno := csno+1; lprint( `Characteristic set produced`,csno,nops(qhi),nops(qsi),nops(qs) ); print(cs) fi; if 0 < `charsets/class`(cs[1],ord) then iss := `charsets/initialset`(cs,ord); if `charsets/simpa`(iss,cs,ord) <> 0 then qsi := {cs,op(qsi)} fi; iss := `charsets/adjoin`(iss,qs,qqi) else iss := {} fi; if 1 < nops(qhi) then qhi := {op(iss),op(2 .. nops(qhi),qhi)} else qhi := iss fi od; if qsi = {} then [] else op(sort(`charsets/contract`(qsi,ord,0),`charsets/lenord`)) fi end: `charsets/fcharser`:=proc(ps,ord,medset) local qs,cs,iss,n,qhi,qsi,factorset,csno,ppi,qqi; if type(ps[1],{set,list}) then qhi := {op(ps)} else qhi := {{op(ps)}} fi; qsi := {}; csno := 0; ppi := {}; qqi := {}; for n from 0 while qhi <> {} do qhi := sort([op(qhi)],`charsets/nopsord`); qs := qhi[1]; ppi := `charsets/select`(ppi,nops(qs)); qqi := {op(ppi[2]),op(qqi)}; if n = 0 then ppi := {} else ppi := {qs,op(ppi[1])} fi; if nops(qs)-3 < nops(ord) then cs := `charsets/fcharseta`([op(qs)],ord,medset); factorset := op(2,cs[2]); cs := cs[1] else cs := `charsets/charseta`([op(qs)],ord,medset); factorset := {} fi; if 1 < printlevel then csno := csno+1; lprint( `Characteristic set produced`,csno,nops(qhi),nops(qsi),nops(qs) ); print(cs) fi; if 0 < `charsets/class`(cs[1],ord) then iss := `charsets/initialset`(cs,ord); if `charsets/simpa`(iss,cs,ord) <> 0 then qsi := {cs,op(qsi)} fi; iss := iss union `charsets/factorps`(factorset) else iss := `charsets/factorps`(factorset) fi; iss := `charsets/adjoin`(iss,qs,qqi); if 1 < nops(qhi) then qhi := {op(iss),op(2 .. nops(qhi),qhi)} else qhi := iss fi od; if qsi = {} then [] else op(sort(`charsets/contract`(qsi,ord,0),`charsets/lenord`)) fi end: `charsets/excharser`:=proc(ps,ord,medset) local qs,cs,is,iss,i,j,qsi,qhi,r,rr; if type(ps[1],{set,list}) then qhi := {ps} else qhi := {[ps,1]} fi; qsi := {}; while qhi <> {} do qs := qhi[1][1]; cs := `charsets/charseta`([op(qs)],ord,medset); if 1 < printlevel then lprint(`Characteristic set produced`); print(cs) fi; if 0 < `charsets/class`(cs[1],ord) then is := `charsets/initialset`(cs,ord); rr := `charsets/nopower`(`charsets/prod`({op(is),qhi[1][2]})); `charsets/premas`(rr,cs,ord); r := `charsets/simp`(%,cs,ord); if r <> 0 then if r = 1 then qsi := {cs,op(qsi)} else qsi := {op(qsi),[cs,`charsets/simpb`(r,rr)]} fi fi else is := [] fi; iss := {}; if nops(ord) <= nops(ps)+1 then for i in is do iss := {op(iss),[{op(qs),i},qhi[1][2]]} od else for i to nops(is) do if i = 1 then 1 else product('is[j]','j' = 1 .. i-1) fi; iss := {op(iss),[{op(qs),is[i]},%*qhi[1][2]]} od fi; if 1 < nops(qhi) then qhi := {op(iss),op(2 .. nops(qhi),qhi)} else qhi := iss fi od; if qsi = {} then [] else op(sort([op(qsi)],`charsets/lenord`)) fi end: `charsets/simp`:=proc(r,cs,ord) local rr,i,fs,j,ind; if r = 0 then 0 else fs := {`charsets/pfactor`(r)}; rr := 1; for i in fs do if 0 < `charsets/class`(i,ord) then ind := 0; for j in cs do subs(i = 0,j); if % = 0 then if `charsets/class`(`charsets/movefactor`(j,i,ord),ord) = 0 then ind := -1; break fi elif `charsets/class`(%,ord) = 0 then ind := 1; break fi od; if ind = 0 then rr := rr*i elif ind = -1 then break fi fi od; if ind = -1 then 0 else rr fi fi end: `charsets/simpa`:=proc(fs,cs,ord) local i,j,ds; if nops(cs) = 1 then 1 else ds := [op(1 .. nops(cs)-1,cs)]; for i in fs do for j in ds do if subs(i = 0,j) = 0 then if `charsets/class`(`charsets/movefactor`(j,i,ord),ord) = 0 then RETURN(0) fi fi od od; 1 fi end: `charsets/simpb`:=proc(a,b) if `charsets/measure`(a) < `charsets/measure`(b) then a else b fi end: `charsets/measure`:=proc(a) local i; if type(a,`^`) then nops(op(1,a)) elif type(a,`*`) then sum('`charsets/measure`(op(i,a))','i' = 1 .. nops(a)) else nops(a) fi end: `charsets/fexcharser`:=proc(ps,ord,medset) local qs,cs,is,iss,n,i,j,qhi,qsi,r,rr,factorset; if type(ps[1],{set,list}) then qhi := {ps} else qhi := {[ps,1]} fi; qsi := {}; for n from 0 while qhi <> {} do qs := qhi[1][1]; if n = 0 then cs := `charsets/fcharseta`([op(qs)],ord,medset); factorset := op(2,cs[2]); cs := cs[1] else cs := `charsets/charseta`([op(qs)],ord,medset); factorset := {} fi; if 1 < printlevel then lprint(`Characteristic set produced`); print(cs) fi; if 0 < `charsets/class`(cs[1],ord) then is := `charsets/initialset`(cs,ord) union `charsets/factorps`(factorset) ; rr := `charsets/nopower`(`charsets/prod`({op(is),qhi[1][2]})); `charsets/premas`(rr,cs,ord); r := `charsets/simp`(%,cs,ord); if r <> 0 then if r = 1 then qsi := {cs,op(qsi)} else qsi := {op(qsi),[cs,`charsets/simpb`(r,rr)]} fi fi else is := `charsets/factorps`(factorset) fi; iss := {}; if nops(ord) <= nops(ps)+1 then for i in is do iss := {op(iss),[{op(qs),i},qhi[1][2]]} od else for i to nops(is) do if i = 1 then 1 else product('is[j]','j' = 1 .. i-1) fi; iss := {op(iss),[{op(qs),is[i]},%*qhi[1][2]]} od fi; if 1 < nops(qhi) then qhi := {op(iss),op(2 .. nops(qhi),qhi)} else qhi := iss fi od; if qsi = {} then [] else op(sort([op(qsi)],`charsets/lenord`)) fi end: `charsets/irrcharser`:=proc(ps,ord,medset,m) local qs,cs,cst,is,iss,n,ts,qsi,qhi,pi,factorset,ppi,qqi,csno,ind,fset,mind; options remember; if nargs = 3 then ind := 1 else ind := m fi; if type(ps[1],{set,list}) then qhi := {op(ps)} else qhi := {ps} fi; qsi := {}; pi := {}; csno := 0; ppi := {}; qqi := {}; if medset = `charsets/basset` then mind := true else mind := false fi; for n from 0 while qhi <> {} do qhi := sort([op(qhi)],`charsets/nopsord`); qs := qhi[1]; ppi := `charsets/select`(ppi,nops(qs)); qqi := {op(qqi),op(ppi[2])}; if n = 0 then ppi := {} else ppi := {qs,op(ppi[1])} fi; if nops(qs)-3 < nops(ord) then if not mind then cs := cat(`charsets/f`,(substring(medset,10 .. length(medset))))( qs,ord,[{},indets(qs)],'fset'); factorset := fset[1] else cs := `charsets/fcharseta`(qs,ord,medset); factorset := op(2,cs[2]); cs := cs[1]; if 1 < printlevel and ind = 1 then csno := csno+1; lprint(`Characteristic set produced`,csno,nops(qhi), nops(qsi),nops(qs)); print(cs) fi fi; if ind = 0 then cs := [`charsets/fcnormal`(cs,ord,2)]; if 1 < nops(cs) then factorset := factorset union op(2,cs[2]) fi; cs := cs[1] fi elif mind then cs := `charsets/removecont`( `charsets/charseta`([op(qs)],ord,medset),ord,'factorset') else cs := `charsets/removecont`(medset([op(qs)],ord),ord,'factorset') fi; if 0 < `charsets/class`(cs[1],ord) then ts := `charsets/irras`(cs,ord,ind); if ts[2] = 0 then if not mind then if medset = `charsets/autored` then cs := `charsets/charseta`({op(qs),op(cs)}, ord,`charsets/charsetn`) elif not `charsets/subset`(cs,qs) then cs := `charsets/charseta`({op(qs),op(cs)}, ord,medset) fi; if 1 < printlevel and ind = 1 then csno := csno+1; lprint(`Characteristic set produced`,csno,nops(qhi), nops(qsi),nops(qs)); print(cs) fi fi; if not member(cs,pi) then pi := {cs,op(pi)}; if 0 < `charsets/class`(cs[1],ord) then ts := `charsets/irras`(cs,ord,ind); if ts[2] = 0 then qsi := {cs,op(qsi)}; if nops(cs) = nops(ord) then is := `charsets/factorps`(factorset) else is := `charsets/initialset`(cs,ord) union `charsets/factorps`(factorset) fi; iss := `charsets/adjoin`(is,qs,qqi) fi else iss := `charsets/adjoin`( `charsets/factorps`(factorset),qs,qqi) fi else iss := `charsets/adjoin`(`charsets/factorps`(factorset),qs,qqi) fi fi; if ts[2] <> 0 then is := `charsets/factorps`(factorset); if 1 < ts[2] then cst := [op(1 .. ts[2]-1,cs)]; is := is union `charsets/initialset`(cst,ord); iss := `charsets/adjoin`(is,qs,qqi) union `charsets/adjoinb`(ts[1],qs,qqi,cst) else iss := `charsets/adjoin`(is union ts[1],qs,qqi) fi fi else iss := `charsets/adjoin`(`charsets/factorps`(factorset),qs,qqi) fi; if 1 < nops(qhi) then qhi := {op(iss),op(2 .. nops(qhi),qhi)} else qhi := iss fi od; if qsi = {} then [] else op(sort(`charsets/contract`(qsi,ord,1),`charsets/lenord`)) fi end: `charsets/subset`:=proc(ps,qs) local p; for p in ps do if not member(p,qs) then RETURN(false) fi od; true end: `charsets/adjoinb`:=proc(is,qs,qh,cs) local iss,i,j,ind,qhi,itt; iss := {}; qhi := qh minus {qs}; if is <> {} then for i in is do itt := {i,op(qs),op(cs)}; ind := 0; if 0 < nops(qhi) then for j in qhi while ind = 0 do if `charsets/subset`(j,itt) then ind := 1 fi od fi; if ind = 0 then iss := {op(iss),itt} fi od fi; iss end: `charsets/irras`:=proc(as,ord,inda,den) local ind,i,j,p,qs,n,fs,ja,dd; options remember,system; ind := 1; ja := 0; dd := 1; for i to nops(as) do p := factor(as[i]); fs := `charsets/dfactors`(p); qs := {}; for j to nops(fs) do if 0 < `charsets/class`(fs[j],ord) then qs := {op(qs),fs[j]} fi od; `charsets/lvar`(p,ord); if `charsets/degree`(qs[1],%) < `charsets/degree`(p,%) then ja := 1; ind := 0; break fi od; if ind = 1 and 1 < nops(as) then for n while n < nops(as) and `charsets/degreel`(as[n],ord) = 1 do od; if n < nops(as) then qs := `charsets/irrassub`(as,n,ord,inda,'ja','dd') else ja := 0 fi fi; if 3 < nargs then den := dd fi; [qs,ja] end: `charsets/irrassub`:=proc(as,n,ord,ind,ja,den) local m,qs,i,vv,dd; global `charsets/das`; for m from n+1 while m <= nops(as) and `charsets/degreel`(as[m],ord) = 1 do od; dd := 1; if m <= nops(as) then vv := `charsets/lvar`(as[m],ord); if ind = -1 then qs := `charsets/factoras`(as[m],['as[i]' $ ('i' = 1 .. m-1)],ord) else qs := `charsets/cfactor`(as[m],['as[i]' $ ('i' = 1 .. m-1)],ord); dd := denom(qs); qs := {`charsets/qfactor`(numer(qs),ord)} fi; if max('`charsets/degree`(qs[i],vv)' $ ('i' = 1 .. nops(qs))) = `charsets/degree`(as[m],vv) then qs := `charsets/irrassub`(as,m,ord,ind,'ja','dd') else ja := m fi else ja := 0 fi; if 5 < nargs then den := dd fi; qs end: `charsets/dfactors`:=proc(q) local qs,j; if type(q,`*`) then qs := {}; for j to nops(q) do if not type(op(j,q),integer) then if type(op(j,q),`^`) then qs := {op(qs),numer(op(1,op(j,q))/lcoeff(op(1,op(j,q))))} else qs := {op(qs),numer(op(j,q)/lcoeff(op(j,q)))} fi fi od elif type(q,`^`) then qs := {numer(op(1,q)/lcoeff(op(1,q)))} else if not type(q,integer) then qs := {numer(q/lcoeff(q))} fi fi; qs end: `charsets/fcnorm`:=proc(cs,ord,m) local as,i; as := [cs[1]]; for i from 2 to nops(cs) do as := `charsets/fcnormal`([op(as),cs[i]],ord,m); if 1 < nops([as]) then RETURN(cs) fi od; as end: `charsets/fcnormal`:=proc(cs,ord,m) local n,ini,i,j,ggg,gg,ff,ccs,dd,cd,fs,nt; n := nops(cs); if n < 2 then cs else dd := cs[n]; nt := nops(expand(dd)); for i from n-1 by -1 to 1 do ini := `charsets/initial`(dd,ord); if 0 < `charsets/degree`(ini,`charsets/lvar`(cs[i],ord)) then ggg := `charsets/gcdex`(cs[i],ini,`charsets/lvar`(cs[i],ord)); gg := ggg[3]; if 0 < `charsets/degree`(gg,`charsets/lvar`(cs[i],ord)) then ff := cs[i]; gg := {`charsets/pfactor`(gg)}; cd := {}; for j to nops(gg) do if `charsets/class`(gg[j],ord) = `charsets/class`(cs[i],ord) then ff := `charsets/nopower`( `charsets/movefactor`(ff,gg[j],ord)); cd := {op(cd),gg[j]} fi od; if `charsets/class`(ff,ord) = 0 then ccs := [[1]] else ccs := [`charsets/fcnormal`(subs(cs[i] = ff,cs),ord,m)] fi; if nops(ccs) = 1 then RETURN(ccs[1],`common divisors` = cd) else RETURN( ccs[1],`common divisors` = {op(cd),op(op(2,ccs[2]))} ) fi else dd := `charsets/prem`(dd*ggg[2],cs[i],`charsets/lvar`(cs[i],ord)) ; dd := `charsets/movefactor`(dd,`charsets/initial`(cs[i],ord),ord) ; dd := `charsets/nopower`(dd); if 4*m*nt < nops(expand(dd)) then RETURN(cs) fi fi fi od; ccs := [op(1 .. n-1,cs)]; fs := {`charsets/pfactor`(content(dd,`charsets/lvar`(dd,ord),'dd'))}; gg := {}; for i to nops(fs) do if 0 < `charsets/class`(fs[i],ord) then gg := {op(gg),fs[i]} fi od; gg := `charsets/prod`(gg); ini := `charsets/initialset`(ccs,ord); for i to nops(ini) do gg := `charsets/movefactor`(gg,ini[i],ord) od; dd := `charsets/nopower`(gg)*dd; gg := [`charsets/pfactor`(numer(dd/lcoeff(expand(dd))))]; dd := 1; for i to nops(gg) do if 0 < `charsets/class`(gg[i],ord) then dd := dd*gg[i] fi od; if m*nt < nops(expand(dd)) then if 1 < printlevel then lprint(`normalization fails:`,nt,nops(expand(dd))) fi; cs else [op(ccs),dd] fi fi end: `charsets/gcdex`:=proc(A,B,x) local m,pm,cc,cd,c,c1,c2,d,d1,d2,r,r1,r2,q,II; options remember,system; if A = 0 then RETURN([0,1,B]) fi; if B = 0 then RETURN([1,0,A]) fi; cc := content(A,x,c); cd := content(B,x,d); II := readlib(`gcd/degrees`)(expand(c),expand(d),{x},'c','d'); pm := 1; c1 := 1; c2 := 0; d1 := 0; d2 := 1; while d <> 0 do r := prem(c,d,x,'m','q'); divide(r,pm,'r'); divide(m*c1-q*d1,pm,'r1'); divide(c2*m-q*d2,pm,'r2'); c := d; c1 := d1; c2 := d2; d := r; d1 := r1; d2 := r2; pm := m od; subs(II,[c1*cd,c2*cc,c*cc*cd]) end: `charsets/exirrcharser`:=proc(ps,ord,medset) local qs,cs,is,iss,n,i,j,qhi,qsi,r,rr,factorset,mind,fset,ind,ts,den; ind:=0; ind:='ind'; if type(ps[1],{set,list}) then qhi := {ps} else qhi := {[ps,1]} fi; if medset = `charsets/basset` then mind := true else mind := false fi; qsi := {}; for n from 0 while qhi <> {} do qs := qhi[1][1]; if not mind then if n < 20 then cs :=cat(`charsets/f`,(substring(medset,10 .. length(medset))))( qs,ord,[{},indets(qs)],'fset'); factorset := fset[1] else cat(`charsets/`,(substring(medset,10 .. length(medset))))(qs,ord); cs := `charsets/removecont`(%,ord,'factorset') fi else if n < 20 then cs := `charsets/fcharseta`([op(qs)],ord,medset); factorset := op(2,cs[2]); cs := cs[1] else cs := `charsets/charseta`([op(qs)],ord,medset); cs := `charsets/removecont`(cs,ord,'factorset') fi; if 1 < printlevel then lprint(`Characteristic set produced`); print(cs) fi fi; if 0 < `charsets/class`(cs[1],ord) then ts := `charsets/irras`(cs,ord,ind,'den'); if ts[2] = 0 then if not mind then if medset = `charsets/autored` then cs := `charsets/charseta`({op(qs),op(cs)}, ord,`charsets/charsetn`) elif not `charsets/subset`(cs,qs) then cs := `charsets/charseta`({op(qs),op(cs)}, ord,medset) fi; if 1 < printlevel then lprint(`Characteristic set produced`); print(cs) fi fi; if 0 < `charsets/class`(cs[1],ord) then ts := `charsets/irras`(cs,ord,ind,'den'); if ts[2] = 0 then is := `charsets/initialset`(cs,ord) union `charsets/factorps`(factorset); if nops(cs) = nops(ord) then rr := `charsets/nopower`(qhi[1][2]) else rr := `charsets/nopower`( `charsets/prod`({op(is),qhi[1][2]})) fi; `charsets/premas`(rr,cs,ord); r := `charsets/simp`(%,cs,ord); if r <> 0 then if r = 1 then qsi := {cs,op(qsi)} else qsi := {op(qsi),[cs,`charsets/simpb`(r,rr)]} fi fi fi else is := `charsets/factorps`(factorset); ts := [1,0] fi fi; if ts[2] <> 0 then if 1 < ts[2] then is := `charsets/initialset`({op(1 .. ts[2]-1,cs)},ord) union `charsets/factorps`(factorset) else is := `charsets/factorps`(factorset) fi fi else is := `charsets/factorps`(factorset); ts := [1,0] fi; iss := {}; if nops(ord) <= nops(ps)+1 then for i in is do iss := {op(iss),[{op(qs),i},qhi[1][2]]} od else for i to nops(is) do if i = 1 then 1 else product('is[j]','j' = 1 .. i-1) fi; iss := {op(iss),[{op(qs),is[i]},%*qhi[1][2]]} od fi; if ts[2] <> 0 and ts[1] <> {} then if not mind then if medset = `charsets/autored` then `charsets/charseta`({op(qs),op(cs)}, ord,`charsets/charsetn`) elif not `charsets/subset`(cs,qs) then `charsets/charseta`({op(qs),op(cs)},ord,medset) else cs fi; if % <> cs then if ts[2] = 1 then cs := qs else cs := {op(qs),op(1 .. ts[2]-1,cs)} fi fi fi; for i in ts[1] do iss := {op(iss),[{i,op(cs)},`charsets/prod`({op(is),qhi[1][2],den})]} od; if 0 < `charsets/class`(den,ord) and ts[2] < `charsets/class`(ts[1][1],ord) then iss := {op(iss),[{op(cs),den},qhi[1][2]]} fi fi; if 1 < nops(qhi) then qhi := {op(iss),op(2 .. nops(qhi),qhi)} else qhi := iss fi od; if qsi = {} then [] else op(sort([op(qsi)],`charsets/lenord`)) fi end: `charsets/select`:=proc(ppi,n) local i,pp,qq; pp := {}; qq := {}; for i in ppi do if n <= nops(i) then qq := {i,op(qq)} else pp := {i,op(pp)} fi od; [pp,qq] end: `charsets/adjoin`:=proc(is,qs,qh) local iss,i,j,ind,qhi,itt; iss := {}; qhi := qh minus {qs}; if is <> {} then for i in is do itt := {i,op(qs)}; ind := 0; if 0 < nops(qhi) then for j in qhi while ind = 0 do if `charsets/subset`(j,itt) then ind := 1 fi od fi; if ind = 0 then iss := {op(iss),itt} fi od fi; iss end: `charsets/adjoina`:=proc(is,qs,qh) local iss,i,j,ind,qhi,itt; iss := {}; qhi := qh minus {qs}; if is <> {} then for i in is do itt := {i,op(qs)}; ind := 0; if 0 < nops(qhi) then for j in qhi while ind = 0 do if `charsets/subset`(j,itt) then ind := 1 fi od fi; if ind = 0 then iss := {op(iss),[i,op(qs)]} fi od fi; iss end: `charsets/nopsord`:= proc(a,b) options remember,system; if nops(b) < nops(a) then true else false fi end: `charsets/contract`:=proc(cs,ord,irr) local i,j,mem,ts; mem := {}; ts := {}; if nops(cs) < 2 then [op(cs)] else for i to nops(cs)-1 do if not member(i,mem) then for j from i+1 to nops(cs) do if not member(j,mem) then if `charsets/linas`(cs[j],ord,irr) and `charsets/contractsub`(cs[i],cs[j],ord) then ts := {cs[j],op(ts)}; mem := {j,op(mem)} else if `charsets/linas`(cs[i],ord,irr) and `charsets/contractsub`(cs[j],cs[i],ord) then ts := {cs[i],op(ts)} fi fi fi od fi od; [op({op(cs)} minus ts)] fi end: `charsets/contractsub`:=proc(cs1,cs2,ord) local i,is; for i in cs1 do if `charsets/premas`(i,cs2,ord) <> 0 then RETURN(false) fi od; is := `charsets/initialset1`(cs1,ord); for i in is do if `charsets/premas`(i,cs2,ord) = 0 then RETURN(false) fi od; true end: `charsets/linas`:=proc(as,ord,irr) local i,j,n,m; if irr = 1 then true elif irr = 2 then if 1 < nops(as) then for n while n < nops(as) and `charsets/degreel`(as[n],ord) = 1 do od; if n < nops(as) then for m from n+1 while m <= nops(as) and `charsets/degreel`(as[m],ord) = 1 do od; if m <= nops(as) then RETURN(false) fi fi fi; true else for i in as do if 1 < `charsets/degreel`(i,ord) then RETURN(false) fi od; if irr = 0 or nops(as) < 2 then true else for i from 2 to nops(as) do for j to i do if has( `charsets/initial`(as[i],ord),`charsets/lvar`(as[j],ord) ) then RETURN(false) fi od od; true fi fi end: `charsets/qirrcharser`:=proc(ps,ord,medset) local qs,cs,is,iss,n,ts,qsi,qhi,pi,factorset,ppi,qqi,csno,fset,mind; if type(ps[1],{set,list}) then qhi := {op(ps)} else qhi := {ps} fi; qsi := {}; pi := {}; csno := 0; ppi := {}; qqi := {}; if medset <> `charsets/basset` and medset <> `charsets/wbasset` then mind := true else mind := false fi; for n from 0 while qhi <> {} do qhi := sort([op(qhi)],`charsets/nopsord`); qs := qhi[1]; ppi := `charsets/select`(ppi,nops(qs)); qqi := {op(qqi),op(ppi[2])}; if n = 0 then ppi := {} else ppi := {qs,op(ppi[1])} fi; if nops(qs)-3 < nops(ord) then if mind then cs := cat(`charsets/f`,(substring(medset,10 .. length(medset))))( qs,ord,[{},indets(qs)],'fset'); factorset := fset[1] else cs := `charsets/fcharseta`(qs,ord,medset); factorset := op(2,cs[2]); cs := cs[1]; if 1 < printlevel then csno := csno+1; lprint(`Characteristic set produced`,csno,nops(qhi), nops(qsi),nops(qs)); print(cs) fi fi elif not mind and 2 < nops(indets(cs[1])) then cs := `charsets/removecont`( `charsets/charseta`([op(qs)],ord,medset),ord,'factorset') else cs := `charsets/removecont`(medset([op(qs)],ord),ord,'factorset') fi; if 0 < `charsets/class`(cs[1],ord) then ts := `charsets/qirras`(cs,ord); if ts[2] = 0 then if mind then if medset = `charsets/autored` then cs := `charsets/charseta`({op(qs),op(cs)}, ord,`charsets/charsetn`) elif not `charsets/subset`(cs,qs) then cs := `charsets/charseta`({op(qs),op(cs)}, ord,medset) fi; if 1 < printlevel then csno := csno+1; lprint(`Characteristic set produced`,csno,nops(qhi), nops(qsi),nops(qs)); print(cs) fi fi; if not member(cs,pi) then pi := {cs,op(pi)}; if 0 < `charsets/class`(cs[1],ord) then ts := `charsets/qirras`(cs,ord); if ts[2] = 0 then qsi := {cs,op(qsi)}; if nops(cs) = nops(ord) then is := `charsets/factorps`(factorset) else is := `charsets/initialset`(cs,ord) union `charsets/factorps`(factorset) fi; iss := `charsets/adjoin`(is,qs,qqi) fi else iss := `charsets/adjoin`( `charsets/factorps`(factorset),qs,qqi) fi else iss := `charsets/adjoin`(`charsets/factorps`(factorset),qs,qqi) fi fi; if ts[2] <> 0 then is := `charsets/factorps`(factorset) union ts[1]; if 1 < ts[2] then is := is union `charsets/initialset`([op(1 .. ts[2]-1,cs)],ord) fi; iss := `charsets/adjoin`(is,qs,qqi) fi else iss := `charsets/adjoin`(`charsets/factorps`(factorset),qs,qqi) fi; if 1 < nops(qhi) then qhi := {op(iss),op(2 .. nops(qhi),qhi)} else qhi := iss fi od; if qsi = {} then [] else op(sort(`charsets/contract`(qsi,ord,2),`charsets/lenord`)) fi end: `charsets/qirras`:=proc(as,ord) local ind,i,j,p,qs,fs,n,m,ja; options remember,system; qs := []; ind := 1; ja := 0; for i to nops(as) do p := factor(as[i]); fs := `charsets/dfactors`(p); qs := {}; for j to nops(fs) do if 0 < `charsets/class`(fs[j],ord) then qs := {op(qs),fs[j]} fi od; `charsets/lvar`(p,ord); if `charsets/degree`(qs[1],%) < `charsets/degree`(p,%) then ja := 1; ind := 0; break fi od; if ind = 1 and 1 < nops(as) then for n while n < nops(as) and `charsets/degreel`(as[n],ord) = 1 do od; if n < nops(as) then for m from n+1 while m <= nops(as) and `charsets/degreel`(as[m],ord) = 1 do od; if m <= nops(as) then lprint( `Warning: factorization over algebraic field required for ics` ) fi fi fi; [qs,ja] end: `charsets/cfactorsub`:=proc(f,as,ord) local ind,i,ff,fn,ffn,lind; ff := numer(f); if type(ff,`^`) then ind := map(`charsets/newfactoras`,ff,as,ord) elif type(ff,`*`) then fn := {op(ff)}; ind := 1; for i to nops(fn) do if type(fn[i],`^`) then ind := ind*map(`charsets/newfactoras`,fn[i],as,ord) else ind := ind*`charsets/newfactoras`(fn[i],as,ord) fi od else ind := `charsets/newfactoras`(ff,as,ord) fi; lind := lcoeff(ff,`charsets/lvar`(ff,ord))/lcoeff(ind,`charsets/lvar`(ind,ord)) ; for i from nops(as) by -1 to 1 do `charsets/premB`(numer(lind),as[i],`charsets/lvar`(as[i],ord),'fn'); `charsets/premB`(denom(lind),as[i],`charsets/lvar`(as[i],ord),'ffn'); lind := %%*ffn/%/fn od; ind*lind; if type(%,{`*`,`^`}) then simplify(%/denom(f)) else f fi end: `charsets/cfactor_switch` := false: `charsets/newfactoras`:=proc(f,as,ord) local pf,pas,aas,con,vf,va,i,fn; global `charsets/con`,`charsets/cfactor_switch`,terms; vf := `charsets/lvar`(f,ord); if `charsets/class`(vf,ord) <= `charsets/class`(as[nops(as)],ord) then RETURN(f) elif `charsets/degree`(f,vf) = 1 then RETURN(f) fi; aas := []; con := 2; for i to nops(as) do `charsets/degreel`(as[i],ord); if 1 < % then con := con,%; aas := [op(aas),expand(as[i])] fi od; if aas = [] then RETURN(f) fi; va := `charsets/lvar`(aas[1],ord); if nops(aas) = 1 and `charsets/degree`(f,va) = 0 then `charsets/trivial`(f,aas[1],vf,va,'fn'); if type(%,{`*`,`^`}) then RETURN(map(`charsets/newfactoras`,%,as,ord)/fn) fi fi; pas := `charsets/fcnorm`(aas,ord,1); `charsets/fcnormal`([op(pas),f],ord,1); if 1 < nops([%]) then pf := f else pf := %[nops(%)] fi; if `charsets/cfactor_switch` and 1 < nops(pas) then fn := `charsets/cfactor`(pf,[pas[nops(pas)]],ord); if max(op(map(degree,[op(numer(fn))],vf))) < degree(pf,vf) then RETURN(`charsets/cfactorsub`(numer(fn),pas,ord)/denom(fn)) else fn := `charsets/cfactor`(pf,[op(1 .. nops(pas)-1,pas)],ord); if max(op(map(degree,[op(numer(fn))],vf))) < degree(pf,vf) then RETURN(`charsets/cfactorsub`(numer(fn),pas,ord)/denom(fn)) fi fi fi; if {con,`charsets/degreel`(pf,ord)} = {2} and nops(pas) < 3 then `charsets/factoras`(pf,pas,ord) else if 1 < printlevel then lprint(`newfactoras: factorization over algebraic field: degree`, `charsets/degreel`(pf,ord),terms,nops(pf)) fi; `charsets/con` := true; `charsets/newfactorassub`(pf,pas,ord) fi end: `charsets/trivial`:=proc(ff,aa,vf,va,fn) local f,a,da,df,ss,i; f := expand(ff); a := expand(aa); da := degree(a,va); df := degree(f,vf); if numer(f/lcoeff(f)) = subs(va = vf,numer(a/lcoeff(a))) and 2 < nops(f) then fn := 1; ss := [coeff(f,vf,da)]; for i from df-1 by -1 to 0 do ss := [ss[1]*va+coeff(f,vf,i),op(ss)] od; sum('ss[i+1]*vf^(i-1)','i' = 1 .. df); RETURN((vf-va)*%) fi; sum('`charsets/_z`^(da-i)*va^i*coeff(a,va,i)','i' = 0 .. da); sum('`charsets/_z`^(df-i)*vf^i*coeff(f,vf,i)','i' = 0 .. df); `charsets/premB`(%,%%,`charsets/_z`,'fn'); if `charsets/degree`(%,`charsets/_z`) = 0 then factor(%) else f fi end: `charsets/premB`:=proc(uu,vv,x,fn) local r,v,dr,dv,l,t,lu,lv,gn; if type(vv/x,integer) then fn := 1; subs(x = 0,uu) else gn := 1; r := expand(uu); dr := degree(r,x); v := expand(vv); dv := degree(v,x); if dv <= dr then l := coeff(v,x,dv); v := expand(v-l*x^dv) else l := 1 fi; if 1 < printlevel and 500 < nops(r)+nops(v) then lprint(`pseudo-division:`,[nops(r),x,dr],[nops(v),x,dv]) fi; while dv <= dr and r <> 0 do gcd(l,coeff(r,x,dr),'lu','lv'); t := expand(x^(dr-dv)*v*lv); if dr = 0 then r := 0 else r := subs(x^dr = 0,r) fi; r := expand(lu*r)-t; gn := gn*lu; dr := degree(r,x) od; fn := gn; r fi end: `charsets/newfactorassub`:=proc(f,as,ord) local nord,mord,cs,ccs,ccs1,cr,i,con,ff,fff,nas,vf,die,der,ci,fs,m,n,inda,indb, is,CS,fmedset,ncs,cc,bb; global `charsets/with`,`charsets/das`,`charsets/con`; options remember; nas := nops(as); vf := `charsets/lvar`(f,ord); if 1 < nas then for i from 2 to nas do `charsets/newfactorassub`(as[i],[op(1 .. i-1,as)],ord) := as[i] od fi; '`charsets/lvar`(as[i],ord)' $ ('i' = 1 .. nas); nord := [vf,%]; mord := [%%,vf]; con := 0; for i from 2 to nops(nord) do con := con+`charsets/degree`(f,nord[i]) od; if con = 0 then if nas = 1 then m := `charsets/degree`(as[1],nord[2]); n := `charsets/degree`(f,vf); if igcd(m,n) = 1 then RETURN(f) fi fi fi; indets({f,op(as)}) minus {op(nord)}; if nops(%) = 1 and nargs = 3 then RETURN(`charsets/tefactor`(f,as,mord,[op(%)])) fi; indb := true; ci := []; cr := []; fff := f; do if indb then if con <> 0 and `charsets/con` then der := 0; ff := f; con := 0 else if `charsets/das` <> [false] then die := `charsets/das`[1]; `charsets/das` := [op(2 .. nops(`charsets/das`),`charsets/das`)] else die := `charsets/die`() fi; der := sum('(i-die-2)*nord[i]','i' = 2 .. nops(nord)); ff := expand(subs(vf = der+vf,f)) fi; `charsets/con` := false; `charsets/with` := {}; if 1 < printlevel then lprint(`Characteristic set computation:`, `charsets/index`([op(as),ff],nord)) fi; fmedset := subs(`charsets/remseta` = `charsets/remsetaA`, `charsets/fqcharsetn` = fmedset,op(`charsets/fqcharsetn`)); is:=nord[nops(nord)]; if nops(as) = 1 and degree(as[1],is) > 5 and degree(ff,is) > 5 and nops(indets(as) union indets(ff)) = 2 then ccs := factor(resultant(ff,as[1],is)); if type(ccs,`+`) then RETURN(f) elif degree(as[1],is) < degree(ff,is) then ccs := [ccs,as[1]] else ccs := [ccs,ff] fi; fs := [{},{}] else ccs := fmedset([op(as),ff],nord,[{},{}],'fs') fi; is := `charsets/initialset`(ccs,nord); fs := `charsets/factorps`(`charsets/movefactorps`(fs[1],is,nord)); cs := {`charsets/qfactor`( factor(`charsets/movefactorps`(ccs[1],is union fs,ord)),nord)} fi; if not indb or cs = {} then `charsets/with` := {}; if fs = {} then bb := is[1]; is := is minus {is[1]} else bb := fs[1]; fs := fs minus {fs[1]} fi; ccs := subs(`charsets/remseta` = `charsets/remsetaA`,`charsets/charsetn`) ([op(as),ff,bb],nord); is := `charsets/initialset`(ccs,nord); cs := {`charsets/qfactor`( factor(`charsets/movefactorps`(ccs[1],is union fs,ord)),nord)}; indb := true fi; if cs <> {} then if `charsets/checkwith`(`charsets/with`,is union fs) then inda := true else inda := false fi; if `charsets/linearas`(ccs,nord) or 1 < nops(cs) then if nops(cs) = 1 then if inda then ccs1 := [cs[1],op(2 .. nops(ccs),ccs)] else ccs1 := `charsets/charseta`( [op(as),ff,cs[1],op(2 .. nops(ccs),ccs)],nord, `charsets/charsetn`) fi; while `charsets/class`(ccs1[1],nord) = 0 do if fs = {} then bb := is[1]; is := is minus {is[1]} else bb := fs[1]; fs := fs minus {fs[1]} fi; ccs1 := `charsets/charseta`( [op(as),ff,bb],nord,`charsets/charsetn`) od; if ccs1 <> ccs then ccs := ccs1; cs := [`charsets/qfactor`(factor(ccs[1]),nord)] fi fi; if nops(cs) = 1 then if `charsets/linearas`(ccs,nord) then cc := `charsets/algcd`(f,subs(vf = vf-der,ccs[1]),as,mord) ; if 0 < `charsets/degree`(cc,vf) then cs := [`charsets/arrange`(is union fs,vf)]; if 1 < printlevel then lprint(`A non-trivial factor found:`,cc) fi; if nops(cs) = 0 then ci := [fff]; fff := 1 else fff := `charsets/divide`(fff,cc,as,mord); ci := [cc] fi; CS := subs(vf = vf-der,cs); for i to nops(CS) do if 2 < printlevel then lprint(`Algebraic GCD:`,nops(CS),i) fi; cc := `charsets/algcd`(fff,expand(CS[i]),as,mord) ; if 0 < `charsets/degree`(cc,vf) then fff := `charsets/divide`(fff,cc,as,mord); if 1 < printlevel then lprint(`A non-trivial factor found:`,cc ) fi; if `charsets/degree`(cc,vf) = 1 then ci := [op(ci),cc] else cr := [op(cr),cc] fi fi od; break else indb := false fi fi else ncs := nops(cs); cs := [`charsets/arrange`(cs,vf), `charsets/arrange`(is union fs,vf)]; CS := subs(vf = vf-der,cs); for i to nops(CS) do if 2 < printlevel then lprint(`Algebraic GCD:`,nops(CS),i) fi; cc := `charsets/algcd`(fff,expand(CS[i]),as,mord); if 0 < `charsets/degree`(cc,vf) then fff := `charsets/divide`(fff,cc,as,mord); if 1 < printlevel then lprint(`A non-trivial factor found:`,cc) fi; if `charsets/degree`(cc,vf) = 1 then ci := [op(ci),cc] else if i <= ncs then con := op(2 .. nops(ccs),ccs); if inda and not `charsets/vanish`(cs[i],is) then con := [cs[i],con] else con := `charsets/charseta`( {op(as),con,ff,cs[i]},nord, `charsets/charsetn`) fi; if con[1] = cs[i] and `charsets/linearas`(con,nord) then ci := [op(ci),cc] else cr := [op(cr),cc] fi else cr := [op(cr),cc] fi fi fi od; if 0 < nops(ci) or 1 < nops(cr) or nops(cr) = 1 and `charsets/degree`(numer(cr[1]),vf) < `charsets/degree`(f,vf) then break fi fi fi fi od; `charsets/degree`(fff,vf); if 1 < % then cr := [op(cr),fff] elif % = 1 then ci := [op(ci),fff] fi; if nops(ci) = 0 then `charsets/prod`( ['`charsets/newfactorassub`(cr[i],as,ord)' $ ('i' = 1 .. nops(cr))] ) elif nops(cr) = 0 then product('ci[i]','i' = 1 .. nops(ci)) else product('ci[i]','i' = 1 .. nops(ci))*`charsets/prod`( ['`charsets/newfactorassub`(cr[i],as,ord)' $ ('i' = 1 .. nops(cr))] ) fi end: `charsets/algcd`:= proc(f,g,as,mord) local nas,fs; nas := nops(as)+1; if `charsets/degree`(f,mord[nas]) = 0 or `charsets/degree`(g,mord[nas]) = 0 then RETURN(1) fi; if 1 < printlevel then lprint(`GCD computation over algebraic field:`, `charsets/index`([f,g],mord),op(`charsets/index`(as,mord))) fi; op(1,[`charsets/fcnormal`( `charsets/fcharsetnA`(as,[f,g],mord,[{},{}],'fs'),mord,2)]); if nops(%) = nas then %[nas] else 1 fi end: `charsets/fcharsetnA`:=proc(as,ps,ord,fset1,fset) local cs,rs,fset2,fset3,nas; nas := nops(as)+1; cs := [op(as),op(`charsets/basset`(ps,ord))]; fset2 := [fset1[1],fset1[2] union `charsets/initialset1`(cs,ord)]; if 0 < `charsets/degree`(cs[nas],ord[nas]) then rs := `charsets/remseta`(`charsets/minus`(ps,cs),cs,ord); rs := `charsets/removefactor`(rs,ord,fset2,'fset3') else fset := fset2; RETURN([1]) fi; if rs = [] then fset := fset3; cs else `charsets/fcharsetnA`(as,[op(rs),cs[nas]],ord,fset3,'fset') fi end: `charsets/Malgcd`:=proc(f,g,as,mord) local nas,i; nas := nops(as); [f,g]; for i from nas by -1 to 1 do subs(mord[i] = RootOf(as[i],mord[i]),%) od; evala(Gcd(op(%))); for i to nas do subs(RootOf(as[i],mord[i]) = mord[i],%) od; % end: `charsets/divide`:=proc(ff,f,as,ord) local m,q; sprem(ff,f,`charsets/lvar`(ff,ord),'m','q'); `charsets/premas`(q,as,ord) end: `charsets/linearas`:=proc(cs,ord) local i; if nops(cs) = 1 then true else for i from 2 to nops(cs) do if 1 < `charsets/degreel`(cs[i],ord) then RETURN(false) fi od; true fi end: `charsets/arrange`:=proc(ps,x) `charsets/reorderb`([op(ps)],`charsets/arrangesub`,x); op(%) end: `charsets/arrangesub`:=proc(f,g,x) if `charsets/degree`(f,x) < `charsets/degree`(g,x) then true else false fi end: `charsets/die` := rand(3 .. 8): `charsets/movefactorps`:=proc(f,ps,ord) local p,ff,i; if not type(f,{set,list}) then ff := f; for p in ps do ff := `charsets/movefactor`(ff,p,ord) od; ff else {'`charsets/movefactorps`(f[i],ps,ord)' $ ('i' = 1 .. nops(f))} fi end: `charsets/vanish`:=proc(q,ps) local p; for p in ps do if divide(q,p) then RETURN(true) fi od; false end: `charsets/qfactor`:=proc(f,ord) local i; if `charsets/class`(f,ord) = 0 then op({}) elif type(f,`^`) then op(1,f); numer(%/lcoeff(%)) elif type(f,`*`) then '`charsets/qfactor`(op(i,f),ord)' $ ('i' = 1 .. nops(f)) else numer(f/lcoeff(f)) fi end: `charsets/qqfactor`:=proc(f,ord) local i; if `charsets/class`(f,ord) = 0 then op({}) elif type(f,`^`) then op(1,f); 'numer(%/lcoeff(%))' $ ('i' = 1 .. op(2,f)) elif type(f,`*`) then '`charsets/qqfactor`(op(i,f),ord)' $ ('i' = 1 .. nops(f)) else numer(f/lcoeff(f)) fi end: `charsets/tefactor`:=proc(f,as,ord,var) local i,j,k,vf,nv,df,inf,ff,gg,ja,js,fs,ffs,hs,gs,sol,ci,dvar,tvar,mm,tt,yvar,das ; global `charsets/das`,`charsets/dasA`,`charsets/dieA`,`charsets/@m`,FAIR; nv := nops(var); vf := ord[nops(ord)]; df := `charsets/degree`(f,vf); if nv = 1 and nops(as) = 1 then `charsets/prem`(f,as[1],var[1]); if % <> f then factor(%); [`charsets/qqfactor`(%,[vf])]; if type(%%,{`*`,`^`}) and max('`charsets/degree`(%[i],vf)' $ ('i' = 1 .. nops(%))) < df then ['op(2,op(1,[`charsets/fcnormal`([as[1],%[i]],ord,2)]))' $ ('i' = 1 .. nops(%))]; RETURN(`charsets/prod`(map(`charsets/newfactoras`,%,as,ord))) fi fi fi; inf := lcoeff(f,vf); for i to nv do `charsets/@m`[i] := 1 od; js := []; fs := []; `charsets/dasA` := [1,-1,2,-2,3,-3,false]; `charsets/dieA` := rand(-10*nv .. 10*nv); das := proc() global `charsets/dasA`,`charsets/dieA`; if `charsets/dasA` <> [false] then `charsets/dasA`[1]; `charsets/dasA` := [op(2 .. nops(`charsets/dasA`),`charsets/dasA`)] ; %% else `charsets/dieA`() fi end: ; tvar := `charsets/noterms`(nv,`charsets/degree`(f,var)); dvar := 0; ci := {}; gg := _y1; yvar := {_y1}; for mm while ci = {} and dvar < tvar do dvar := `charsets/noterms`(nv,mm); while nops(js) < dvar do sol := 'var[i] = das()' $ ('i' = 1 .. nv); if subs(sol,inf) <> 0 and `charsets/isirr`(subs(sol,as),ord) then js := [op(js),{sol}]; ff := {`charsets/qfactor`( `charsets/cfactor`(subs(sol,f),subs(sol,as),ord),[vf])}; if max('`charsets/degree`(ff[i],vf)' $ ('i' = 1 .. nops(ff))) = df then RETURN(f) else fs := [op(fs),ff] fi fi od; sum('var[i]','i' = 1 .. nops(var)); coeffs(expand(%^mm),var,'tt'); tt := [tt]; nops(gg); gg := gg+sum('cat(_y,(%+i))*tt[i]','i' = 1 .. nops(tt)); yvar := yvar union {'cat(_y,(%%+i))' $ ('i' = 1 .. nops(tt))}; hs := {}; ffs := fs; for j to nops(fs[1]) do gs := [fs[1][j]]; if 1 < nops(fs) then for i from 2 to nops(fs) do `charsets/getclose`(fs[1][j],[op(fs[i])],ord,'ja'); if % = FAIR then if 1 < printlevel then lprint(`Heuristic tefactor failed`) fi; RETURN(`charsets/newfactorassub`(f,as,ord,0)) else gs := [op(gs),%] fi; fs[i] minus {fs[i][ja]}; if i = nops(fs) then fs := [op(1 .. i-1,fs),%] else fs := [op(1 .. i-1,fs),%,op(i+1 .. nops(fs),fs)] fi od fi; sol := {'subs(op(js[k]),gg)-gs[k]' $ ('k' = 1 .. nops(js))}; sol := {solve(sol,yvar)}; if sol <> {} then hs := hs union {expand(subs(op(sol),gg))} fi od; fs := ffs; ff := f; for j in hs do `charsets/divideA`(ff,j,as,ord); if % <> false then ff := %; ci := ci union {j}; if 1 < printlevel then lprint(`A factor found:`,j) fi fi od; if ci = {} and mm <= 1 and _help <> true and nops(ffs[1])^nops(ffs) <= 128 then gs := `charsets/getall`(ffs); for i to nops(gs) while nops(ci) <= nops(fs[1]) do sol := {'subs(op(js[k]),gg)-gs[i][k]' $ ('k' = 1 .. nops(js))}; sol := {solve(sol,yvar)}; if 1 < printlevel and sol = {} then lprint(sol,yvar) fi; if sol <> {} then sol := expand(subs(op(sol),gg)); `charsets/divideA`(ff,sol,as,ord); if % <> false then ff := %; ci := ci union {sol}; if 1 < printlevel then lprint(`A non-trivial factor found:`,j) fi fi fi od fi od; if ci <> {} then ci := ci union {ff}; `charsets/prod`(map(`charsets/newfactoras`,ci,as,ord)) else if 1 < printlevel then lprint(`Heuristic tefactor failed`) fi; `charsets/newfactorassub`(f,as,ord,0) fi end: `charsets/noterms`:=proc(n,d) local i,j; sum('product('n+j-1','j' = 1 .. i)/i!','i' = 1 .. d)+1 end: `charsets/isirr`:=proc(as,ord) local xa,fs,f,as1; if nops(as) = 1 then xa := `charsets/lvar`(as[1],ord); if xa = 0 then false elif `charsets/degree`(as[1],xa) = 1 then true else fs := {`charsets/qfactor`(factor(as[1]),[xa])}; if nops(fs) = 1 then true else false fi fi else as1 := [op(1 .. nops(as)-1,as)]; if not `charsets/isirr`(as1,ord) then false else f := as[nops(as)]; xa := `charsets/lvar`(f,ord); if xa = 0 then false elif `charsets/degree`(f,xa) = 1 then true else fs := {`charsets/qfactor`(`charsets/cfactor`(f,as1,ord),[xa])}; if nops(fs) = 1 then true else false fi fi fi fi end: `charsets/getall`:=proc(fs) local gs,i,j,nf; if nops(fs) = 1 then {'[fs[1][i]]' $ ('i' = 1 .. nops(fs[1]))} else nf := nops(fs); gs := {op(1 .. nf-1,fs)}; gs := `charsets/getall`(gs); {''[op(gs[i]),fs[nf][j]]' $ ('j' = 1 .. nops(fs[nf]))' $ ('i' = 1 .. nops(gs))} fi end: `charsets/getclose`:=proc(g,fs,var,ja) local i,j,gs,hs,dg,nv,vv,jaa,jbb,ts,cv,rr,chs,df; global FAIR; if `charsets/class`(g,var) = 0 then RETURN(FAIR) fi; if nops(fs) = 1 then ja := 1; RETURN(fs[1]) fi; if nops(var) = 1 and _help <> true then gs := []; jaa := []; dg := `charsets/degree`(g,var[1]); for i to nops(fs) do if fs[i] = g or fs[i] = -g then ja := i; RETURN(fs[i]) elif `charsets/degree`(fs[i],var[1]) = dg then jaa := [op(jaa),i]; gs := [op(gs),fs[i]] fi od; if nops(jaa) = 1 then ja := jaa[1]; RETURN(gs[1]) fi; hs := []; jbb := []; cv := [coeffs(expand(g),var[1],'dg')]; for i to nops(jaa) do coeffs(expand(gs[i]),var[1],'df'); if {df} = {dg} then jbb := [op(jbb),jaa[i]]; hs := [op(hs),gs[i]] fi od; if nops(jbb) = 1 then ja := jbb[1]; RETURN(hs[1]) fi; gs := []; jaa := []; dg := ['sign(cv[j])' $ ('j' = 1 .. nops(cv))]; for i to nops(jbb) do ts := [coeffs(expand(hs[i]),var[1])]; ts := ['sign(ts[j])' $ ('j' = 1 .. nops(ts))]; if `charsets/close`(ts,dg) = 0 then jaa := [op(jaa),jbb[i]]; gs := [op(gs),hs[i]] fi od; if nops(gs) = 1 then ja := jaa[1]; RETURN(gs[1]) fi; gs := []; jaa := []; for i to nops(jbb) do ts := [coeffs(expand(hs[i]),var[1])]; ts := ['sign(ts[j])' $ ('j' = 1 .. nops(ts))]; if `charsets/close`(ts,dg) = 1 then jaa := [op(jaa),jbb[i]]; gs := [op(gs),hs[i]] fi od; if nops(gs) = 1 then ja := jaa[1]; RETURN(gs[1]) fi; gs := []; jaa := []; for i to nops(jbb) do ts := [coeffs(expand(hs[i]),var[1])]; ts := ['sign(ts[j])' $ ('j' = 1 .. nops(ts))]; if `charsets/close`(ts,dg) = 2 then jaa := [op(jaa),jbb[i]]; gs := [op(gs),hs[i]] fi od; if nops(gs) = 1 then ja := jaa[1]; RETURN(gs[1]) else RETURN(FAIL) fi else if _help <> true then nv := nops(var); vv := var[nv]; gs := []; jaa := []; dg := `charsets/degree`(g,vv); for i to nops(fs) do if fs[i] = g or fs[i] = -g then ja := i; RETURN(fs[i]) elif `charsets/degree`(fs[i],vv) = dg then jaa := [op(jaa),i]; gs := [op(gs),fs[i]] fi od; if nops(jaa) = 1 then ja := jaa[1]; RETURN(gs[1]) fi; hs := []; jbb := []; hs := []; chs := []; cv := [coeffs(expand(g),vv,'dg')]; for i to nops(gs) do chs := [op(chs),[coeffs(expand(gs[i]),vv,'df')]]; if {df} = {dg} then jbb := [op(jbb),jaa[i]]; hs := [op(hs),gs[i]] fi od; if nops(jbb) = 1 then ja := jbb[1]; RETURN(hs[1]) fi; for i to nops([dg]) do ts := ['chs[j][i]' $ ('j' = 1 .. nops(chs))]; `charsets/getclose`( cv[i],ts,['var[j]' $ ('j' = 1 .. nv-1)],'jbb'); if % = FAIR then RETURN(FAIR) elif % <> FAIL then ja := jaa[jbb]; RETURN(hs[jbb]) fi od else hs := fs; jaa := ['i' $ ('i' = 1 .. nops(hs))] fi; if hs <> [] then # RETURN(FAIR); lprint(` `); lprint(`Please help to choose one polynomial in the list`); print(hs); lprint(`which is closest to the polynomial`); print(g); rr := readstat(`and enter the polynomial number in the list: `); ja := jaa[rr]; RETURN(hs[rr]) else ja := 1; RETURN(fs[1]) fi fi end: `charsets/close`:=proc(ps,qs) local i,m; if ps = qs then 0 else m := 0; for i to nops(ps) do if ps[i] <> qs[i] then m := m+1 fi od; m fi end: `charsets/divideA`:=proc(ff,f,as,ord) local m,q; if `charsets/class`(ff,ord) <> `charsets/class`(f,ord) then RETURN(false) fi; sprem(ff,f,`charsets/lvar`(ff,ord),'m','q'); if `charsets/premas`(%,as,ord) <> 0 then RETURN(false) fi; `charsets/premas`(q,as,ord) end: `charsets/factoras`:=proc(pf,pas,ord) local df,r,s,ff,i,t,fact,gg,hh,fg,z,ind,as,nas,f,vf,sol,con,m,n,mord,nord; global `charsets/das`,`charsets/con`,`charsets/@m`,FAIR; options remember; if not type(pf,polynom) then RETURN(`charsets/factoras`(numer(pf),pas,ord)/denom(pf)) fi; nas := nops(pas); vf := `charsets/lvar`(pf,ord); if 1 < nas then for i from 2 to nas do `charsets/factoras`(pas[i],[op(1 .. i-1,pas)],ord) := pas[i] od fi; '`charsets/lvar`(pas[i],ord)' $ ('i' = 1 .. nas); nord := [vf,%]; mord := [%%,vf]; con := 0; for i from 2 to nops(nord) do con := con+`charsets/degree`(pf,nord[i]) od; if con = 0 then if nas = 1 then m := `charsets/degree`(pas[1],nord[2]); n := `charsets/degree`(pf,vf); if igcd(m,n) = 1 then RETURN(pf) fi fi fi; indets({pf,op(pas)}) minus {op(nord)}; if nops(%) = 1 and nargs = 3 then RETURN(`charsets/tsfactor`(pf,pas,mord,[op(%)])) fi; ind := 0; as := []; r := 0; f := expand(pf); for i to nops(pas) do df := `charsets/degreel`(pas[i],ord); if 1 < df then r := r+1; `charsets/@m`[r] := df; as := [op(as),expand(pas[i])] fi od; z := `charsets/lvar`(f,ord); df := `charsets/degree`(f,z); if df = 1 then f elif r = 0 then f else if 1 < printlevel then lprint(`factoras: factorization over algebraic field -- degree `, `charsets/degreel`(f,ord)) fi; for s to trunc(1/2*df) while ind = 0 do for i to s do `charsets/g`||i := `charsets/summ`( `charsets/@g`[i,'cat(`charsets/@k`,t)' $ ('t' = 1 .. r)]*product ('`charsets/lvar`(as[t],ord)^cat(`charsets/@k`,t)','t' = 1 .. r) ,r) od; for i to df-s do `charsets/h`||i := `charsets/summ`( `charsets/@h`[i,'cat(`charsets/@k`,t)' $ ('t' = 1 .. r)]*product ('`charsets/lvar`(as[t],ord)^cat(`charsets/@k`,t)','t' = 1 .. r) ,r) od; `charsets/g`||0 := 1; `charsets/h`||0 := 1; gg := eval(sum('cat(`charsets/g`,i)*z^(s-i)','i' = 0 .. s)); hh := eval(sum('cat(`charsets/h`,i)*z^(df-s-i)','i' = 0 .. df-s)); ff := numer(f-lcoeff(expand(f),z)*expand(gg*hh)); ff := expand(`charsets/premas`(ff,as,ord)); fact := {}; for i from 0 to df-1 do fact := {op(fact),op(`charsets/coeff`(as,{coeff(ff,z,i)},r,ord))} od; sol := [`charsets/solveps`(fact,`charsets/getvars`(fact))]; if sol = [FAIR] then if 1 < printlevel then lprint(`factoras changes to newfactoras`) fi; `charsets/das` := [-1,1,-2,2,-3,false]; `charsets/con` := true; RETURN(`charsets/newfactorassub`(pf,pas,ord)) fi; fg := f; if sol <> [] then fg := subs(sol[1],gg)* `charsets/factoras`(numer(subs(sol[1],hh)),as,ord); ind := 1 fi od; numer(fg) fi end: `charsets/tsfactor`:=proc(f,as,ord,var) local i,vf,nv,df,inf,ff,sol,das; nv := nops(var); vf := ord[nops(ord)]; df := `charsets/degree`(f,vf); if nv = 1 and nops(as) = 1 then `charsets/prem`(f,as[1],var[1]); if % <> f then factor(%); [`charsets/qqfactor`(%,[vf])]; if type(%%,{`*`,`^`}) and max('`charsets/degree`(%[i],vf)' $ ('i' = 1 .. nops(%))) < df then ['op(2,op(1,[`charsets/fcnormal`([as[1],%[i]],ord,2)]))' $ ('i' = 1 .. nops(%))]; RETURN(`charsets/prod`(map(`charsets/factoras`,%,as,ord))) fi fi fi; inf := lcoeff(expand(f),vf); das := rand(-2*nv .. 3*nv+nops(ord)); sol := 'var[i] = i+1' $ ('i' = 1 .. nv); while subs(sol,inf) = 0 or not `charsets/isirr`(subs(sol,as),ord) do sol := 'var[i] = das()' $ ('i' = 1 .. nv) od; ff := {`charsets/qfactor`(`charsets/cfactor`(subs(sol,f),subs(sol,as),ord),[vf])} ; if max('`charsets/degree`(ff[i],vf)' $ ('i' = 1 .. nops(ff))) = df then RETURN(f) fi; if 1 < printlevel then lprint(`Heuristic tsfactor failed`) fi; `charsets/factoras`(f,as,ord,0) end: `charsets/summ`:= proc(ss,r) global `charsets/@m`; if r = 1 then sum(ss,cat(`charsets/@k`,r) = 0 .. `charsets/@m`[r]-1) else sum(`charsets/summ`(ss,r-1),cat(`charsets/@k`,r) = 0 .. `charsets/@m`[r]-1) fi end: `charsets/coeff` := proc(as,ss,r,ord) global `charsets/@m`; local k,i,j,qs; qs := ss; for j from r by -1 to 1 do qs := {''coeff(qs[i],`charsets/lvar`(as[j],ord),k)' $ ('k' = 0 .. `charsets/@m`[j]-1)' $ ('i' = 1 .. nops(qs))} od; qs minus {0} end: `charsets/getvars`:=proc(as) local ind,ind1,i; if type(as,{set,list}) then {'op(`charsets/getvars`(as[i]))' $ ('i' = 1 .. nops(as))} else ind := {}; ind1 := indets(as); for i in ind1 do if type(i,indexed) then if op(0,i) = `charsets/@g` or op(0,i) = `charsets/@h` then ind := {op(ind),i} fi fi od; ind fi end: `charsets/solveps`:=proc(ps,lst) local cs,ord,sol,j,phi,qs,qs1,n,factorset; global FAIR; options remember; if 1 < printlevel then lprint( `solveps: trying rational solutions of equations:`,nops(ps),nops(lst) ) fi; if 3 < printlevel then lprint(op(`charsets/index`([op(ps)],[op(lst)]))) fi; ord := `charsets/reorder`([op(lst)],`charsets/degord`,ps); sol := {}; cs := `charsets/fqcharsetn`(ps,ord,[{},{}],'factorset'); factorset := factorset[1]; phi := {ps}; for n while phi <> {} do if 6 < n+nops(phi) or sol = {FAIR} then RETURN(FAIR) fi; if sol <> {} then break fi; if 1 < n then cs := `charsets/charseta`(phi[1],ord,`charsets/wcharsetn`); factorset := {} fi; sol := {op(sol),`charsets/solveasr`(cs,ord,'qs1')}; if n = 1 then sol := `charsets/verify`(sol,ps,ord) fi; qs := `charsets/factorps`(qs1) union `charsets/factorps`(factorset); `charsets/qbasset`(qs,ord); qs := [op(%),op(qs minus {op(%)})]; if qs <> {} then if 1 < nops(phi) then phi := {op(2 .. nops(phi),phi), '[op(phi[1]),qs[j]]' $ ('j' = 1 .. nops(qs))} else phi := {'[op(phi[1]),qs[j]]' $ ('j' = 1 .. nops(qs))} fi else if 1 < nops(phi) then phi := {op(2 .. nops(phi),phi)} else phi := {} fi fi od; if sol = {} then op({}) else op(sol) fi end: `charsets/verify`:=proc(sol,ps,ord) local i,j,sss; for i to nops(sol) do if simplify(subs(sol[i],ps)) = {0} then RETURN({sol[i]}) else {'op(1,sol[i][j])-op(2,sol[i][j])' $ ('j' = 1 .. nops(sol[i]))}; sss := {`charsets/solveps`(ps union %,ord)}; if sss <> {} then RETURN(sss) fi fi od; {} end: `charsets/trisersub`:= proc(ps,ord) local qs,cs,iss,n,i,qhi,qsi,factorset,csno,ppi,qqi,ind,mem; options remember; ind := 0; for i to nops(ps) do if nops(expand(ps[i])) < 3 then ind := 1; break fi od; if ind = 1 then cs := `charsets/fcharseta`([op(ps)],ord,`charsets/charsetn`) else cs := `charsets/fcharseta`([op(ps)],ord,`charsets/qcharsetn`) fi; factorset := op(2,cs[2]); cs := cs[1]; qhi := {{op(ps)}}; qsi := {}; csno := 0; ppi := {}; qqi := {}; for n from 0 while qhi <> {} do qhi := sort([op(qhi)],`charsets/nopsord`); qs := qhi[1]; ppi := `charsets/select`(ppi,nops(qs)); qqi := {op(qqi),op(ppi[2])}; if n = 0 then ppi := {} else ppi := {qs,op(ppi[1])}; ind := 0; for i to nops(ps) do if nops(expand(ps[i])) < 3 then ind := 1; break fi od; if ind = 1 then cs := `charsets/nopower`( `charsets/charseta`(qs,ord,`charsets/charsetn`)); factorset := {} elif qs <> mem and 4 < `charsets/degree`(qs[1],ord) then cs := `charsets/fcharseta`(qs,ord,`charsets/qcharsetn`); factorset := op(2,cs[2]); cs := cs[1] elif nops(qs)-3 < nops(ord) then cs := `charsets/fcharseta`(qs,ord,`charsets/wcharsetn`); factorset := op(2,cs[2]); cs := cs[1] else cs := `charsets/nopower`( `charsets/charseta`(qs,ord,`charsets/wcharsetn`)); factorset := {} fi fi; mem := qs; if 1 < printlevel then csno := csno+1; lprint( `Characteristic set produced`,csno,nops(qhi),nops(qsi),nops(qs) ); print(cs) fi; if 0 < `charsets/class`(cs[1],ord) then iss := `charsets/initialset`(cs,ord); if `charsets/simpa`(iss,cs,ord) <> 0 then qsi := {cs,op(qsi)} fi; iss := iss union `charsets/factorps`(factorset) else iss := `charsets/factorps`(factorset) fi; iss := `charsets/adjoina`(iss,qs,qqi); if 1 < nops(qhi) then qhi := {op(iss),op(2 .. nops(qhi),qhi)} else qhi := iss fi od; if qsi = {} then [] else op(sort(`charsets/contract`(qsi,ord,-1),`charsets/lenord`)) fi end: `charsets/solveas`:=proc(cs,ord) local is,ss,sol,solm,i,j,k; sol := {solve({cs[1]},{`charsets/lvar`(cs[1],ord)})}; if 1 < nops(cs) then for i from 2 to nops(cs) do is := `charsets/initial`(cs[i],ord); solm := {}; for j to nops(sol) do ss := {solve({subs(sol[j],cs[i])},{`charsets/lvar`(cs[i],ord)})}; for k to nops(ss) do if simplify(subs(op(sol[j]),ss[k],is)) <> 0 then solm := {op(solm),{op(sol[j]),op(ss[k])}} fi od od; sol := solm od fi; op(sol) end: `charsets/solveasr`:=proc(cs,ord,qs) local is,ss,ts,sol,solm,i,j,k; ts := {}; if 0 < `charsets/class`(cs[1],ord) then sol := {`charsets/solvel`(cs[1],`charsets/lvar`(cs[1],ord))}; if 1 < nops(cs) then for i from 2 to nops(cs) do if 1 <= nops(sol) then is := `charsets/initial`(cs[i],ord); solm := {}; for j to nops(sol) do if simplify(subs(op(sol[j]),is)) = 0 then ts := {is,op(ts)} else ss := {`charsets/solvel`( subs(sol[j],cs[i]),`charsets/lvar`(cs[i],ord))} ; solm := {op(solm), '{op(sol[j]),op(ss[k])}' $ ('k' = 1 .. nops(ss)) } fi od; sol := solm else break fi od fi else sol := {} fi; if 2 < nargs then qs := ts fi; op(sol) end: `charsets/solvel`:= proc(f,x) local g,i,sol; sol := {}; if nops(indets(f)) = 1 then g := `charsets/getfactor`(f,x); for i in g do sol := {op(sol),solve({i},{x})} od else g := `charsets/factorps`({numer(f)}); for i in g do if `charsets/degree`(i,x) = 1 then sol := {op(sol),solve({i},{x})} fi od fi; op(sol) end: `charsets/getfactor`:=proc(f,x) local q,qs,j; q := `charsets/getfact`(f,x); qs := {}; if type(q,`*`) then for j to nops(q) do if not type(op(j,q),integer) then if type(op(j,q),`^`) then qs := {op(qs),numer(op(1,op(j,q))/lcoeff(op(1,op(j,q))))} else qs := {op(qs),numer(op(j,q)/lcoeff(op(j,q)))} fi fi od elif type(q,`^`) then qs := {op(qs),numer(op(1,q)/lcoeff(op(1,q)))} else if not type(q,integer) then qs := {numer(q/lcoeff(q)),op(qs)} fi fi; [op(qs)] end: `charsets/getfact`:=proc(ff,x) local i,f; if `charsets/degree`(ff,x) = 1 then RETURN(ff) fi; f := convert(ff,`sqrfree/sqrfree`,x); if type(f,`^`) then readlib(factor); readlib(`factor/polynom`); readlib(`factor/unifactor`); readlib(`factor/linfacts`)(expand(op(1,f)),x) elif type(f,`*`) then {'`charsets/getfact`(op(i,f),x)' $ ('i' = 1 .. nops(f))}; product(op(i,%),i = 1 .. nops(%)) else readlib(factor); readlib(`factor/polynom`); readlib(`factor/unifactor`); readlib(`factor/linfacts`)(expand(numer(f)),x) fi end: `charsets/irrvardec`:=proc(ps,ord,medset) local phi,psi,i,j,mem,qq,qs; qq := nops(ps); mem := {}; if 1 < printlevel and nargs < 3 then lprint(`Variable order chosen:`,ord) fi; if nargs <> 4 then psi := [`charsets/irrcharser`(ps,ord,medset)] else psi := [`charsets/exirrcharser`([ps,1],ord,medset)]; if psi <> [[]] then phi := psi; psi := op([]); for i in phi do if type(i[1],list) then psi := psi,i[1] else psi := psi,i fi od; psi := [psi] fi fi; phi := []; for i to nops(psi) do if nops(psi[i]) <= qq then phi := [op(phi),psi[i]] fi od; phi := sort(phi,proc(a,b) if nops(a) < nops(b) then true else false fi end) ; psi := []; if phi <> [[]] then if nops(phi) = 1 and nargs <= 4 then RETURN([op(ps)]) fi; if 1 < printlevel then lprint(`ivd: ics finished at`,time(),nops(phi)) fi; for i to nops(phi) do if not member(i,mem) then qs := `charsets/primebasis`(phi[i],ord); if 1 < printlevel then lprint(`ivd: finite basis found:`,nops(phi),i); print(qs) fi; for j from i+1 to nops(phi) do if not member(j,mem) and nops(phi[i]) <> nops(phi[j]) then if `charsets/remseta`(qs,phi[j],ord) = {} then mem := {j,op(mem)} fi fi od; psi := [op(psi),qs] fi od; op(sort([op(psi)],`charsets/lenord`)) else [] fi end: `charsets/primebasis` := proc(cs,ord) local is; if nops(cs) = nops(ord) then is := {} else is := `charsets/initialset`(cs,ord) fi; `charsets/saturbasis`(cs,is,ord) end: `charsets/saturbasis`:=proc(ps,js,ord) local qs,gb,zz,j; if js <> {} then qs := [op(ps),'cat(`charsets/@z`,j)*js[j]-1' $ ('j' = 1 .. nops(js))]; zz := ['cat(`charsets/@z`,(nops(js)-j+1))' $ ('j' = 1 .. nops(js))]; gb := Groebner['gbasis']( qs,plex(op(zz),'ord[nops(ord)-j+1]' $ ('j' = 1 .. nops(ord)))); qs := []; for j to nops(gb) do if {op(zz)} minus indets(gb[j]) = {op(zz)} then qs := [op(qs),gb[j]] fi od else qs := ps fi; qs end: `charsets/pridealdec`:=proc(qs,ord,medset) local ps,phi,psi,fs,pi,ph,e,i,ss,S,j,gb,gb1,urd,vrd,f,gs,k,n,ind,indb,con,sep, wtd,tc,tco; ps := {op(`charsets/expand`(qs))} minus {0}; ps := Groebner['gbasis'](ps,plex(op(`charsets/reverse`(ord)))); phi := [[ps,1,0,{1}]]; psi := {}; con := false; for n while phi <> [] do fs := phi[1]; phi := [op(2 .. nops(phi),phi)]; sep := fs[2]; wtd := fs[3]; tco := fs[4]; tc := tco; fs := fs[1]; ind := true; if 1 < n then fs := Groebner['gbasis'](fs,plex(op(`charsets/reverse`(ord)))); if member(1,fs) then ind := false elif `charsets/contain`(fs,ph,ord,'plex') then ind := false fi fi; if ind then pi := [`charsets/irrvardec`(fs,ord,medset,0,0)]; if pi <> [[]] then e := nops(pi); if 1 < printlevel then lprint(`pid: ivd finished at`,time(),e) fi; for i to e do ss := {}; S := {}; if 1 < e then for j to e do if j <> i then ss := ss union {`charsets/takele`( {op(pi[j])} minus {op(pi[i])},pi[i],ord)} fi od; ss := `charsets/prod`(`charsets/factorps`(ss)); gb := `charsets/saturbasis`(fs,{ss},ord) else ss := 1; gb := Groebner['gbasis'](fs,plex(op(`charsets/reverse`(ord)))) fi; urd := `charsets/maxindset`(pi[i],ord); if urd = [] then f := 1 else vrd := op({op(ord)} minus {op(urd)}); gb1 := Groebner['gbasis'](gb,plex(vrd,op(urd))); f := lcm('op(1,[Groebner['leadmon'](gb1[j],plex(vrd))])' $ ('j' = 1 .. nops(gb1))) fi; f := `charsets/prod`(`charsets/factorps`({f})); gs := `charsets/saturbasis`(gb,{f},ord); k := `charsets/exponent`(gb,f,gs,ord); if 1 < printlevel then lprint(`pid: primary component found:`,nops(psi)+1, nops(phi),e,i); print(gs) fi; if psi = {} then psi := {[gs,pi[i]]}; ph := gs else psi union {[gs,pi[i]]}; if % <> psi then psi := %; ph := `charsets/idealint`(ph,gs,ord); if `charsets/contain`(ps,ph,ord,'plex') then con := true; break fi fi fi; tco := `charsets/idealint`( tco,`charsets/saturbasisR`(gs,{sep*ss},ord),ord); if 0 < `charsets/class`(f,ord) then if `charsets/class`(sep*ss,ord) = 0 then {op(gb)} union {f^k} else `charsets/saturbasis`( {op(gb)} union {f^k},{sep*ss},ord) fi; indb := not `charsets/contain`( %,`charsets/idealintR`(ph,tc,ord),ord,'plex'); if indb then indb := `charsets/satisfy`({op(gb)} union {f^k},sep*ss,ord) ; if indb then phi := [ op(phi),[{op(gb)} union {f^k},sep*ss,wtd+2,tc] ] fi fi fi; S := S union {ss^`charsets/exponent`(fs,ss,gb,ord)} od; if con then break fi; S := {op(fs)} union S; if not member(1,S) then if `charsets/class`(sep,ord) = 0 then S else `charsets/saturbasis`(S,{sep},ord) fi; not `charsets/contain`( %,`charsets/idealintR`(ph,tco,ord),ord,'plex'); if % then phi := [op(phi),[S,sep,wtd+2,tco]] fi fi fi; phi := sort(phi,`charsets/wtdorder`) fi od; if psi = {} or psi = {[1]} then [] else op( sort([`charsets/remred`([op(psi minus {[1]})],ord)],`charsets/lenord`) ) fi end: `charsets/wtdorder` := proc(a,b) if a[3] < b[3] then true else false fi end: `charsets/satisfy`:=proc(ps,g,ord) if `charsets/class`(g,ord) = 0 then true elif member(g,ps) or member(g^2,ps) or member(g^3,ps) then false elif `charsets/contain`(ps,{g},ord,'plex') then false elif `charsets/eicsTest`([ps,g],ord,`charsets/charsetn`) = [] then false else true fi end: `charsets/remred`:=proc(phi,ord) op(`charsets/remrsub`(`charsets/reverse`(`charsets/remrsub`(phi,ord)),ord)) end: `charsets/remrsub`:=proc(phi,ord) local ph,mem,i,j; ph := []; mem := {}; for i to nops(phi) do if not member(i,mem) then for j from i+1 to nops(phi) do if not member(j,mem) then if `charsets/contain`(phi[j][1],phi[i][1],ord,'plex') then mem := {op(mem),j} fi fi od; ph := [op(ph),phi[i]] fi od; ph end: `charsets/contain`:=proc(ps,qs,ord,pt) local X,q; if 1 < printlevel then lprint(`containment test starts at`,time()) fi; X := `charsets/reverse`(ord); for q in qs do if Groebner['normalf'](q,ps,pt(op(X))) <> 0 then RETURN(false) fi od; true end: `charsets/reverse`:=proc(ord) local j; ['ord[nops(ord)-j+1]' $ ('j' = 1 .. nops(ord))] end: `charsets/takele`:=proc(gs,gb,ord) local X,ps,i; X := `charsets/reverse`(ord); ps := `charsets/reverse`(sort([op(gs)],`charsets/nopsord`)); for i to nops(ps) do if Groebner['normalf'](ps[i],gb,plex(op(X))) <> 0 then RETURN(ps[i]) fi od; ERROR(`Check Here!`) end: `charsets/maxindset`:=proc(gb,ord) local x,g,ind,u,urd; urd := {}; for x in ord do ind := true; u := urd union {x}; for g in gb do if indets(op(2,[Groebner['leadmon'](g, plex(op(`charsets/reverse`(ord))))])) minus u = {} then ind := false; break fi od; if ind then urd := u fi od; [op(urd)] end: `charsets/exponent`:=proc(ps,f,qs,ord) local k; for k do if 50 < k then lprint(`Check PID exponent!`,k) fi; if {op(`charsets/idealquo`(ps,f^k,ord))} = {op(qs)} then break fi od; if 2 < printlevel then lprint(`pid: exponent found:`,k) fi; k end: `charsets/idealquo`:=proc(ps,f,ord) local qs,j; if type(f,{set,list}) then qs := `charsets/idealquo`(ps,f[1],ord); if 1 < nops(f) then [op(2 .. nops(f),f)]; qs := `charsets/idealint`(qs,`charsets/idealquo`(ps,%,ord),ord) fi; qs else qs := `charsets/idealint`(ps,f,ord); Groebner['gbasis']({'simplify(qs[j]/f)' $ ('j' = 1 .. nops(qs))}, plex(op(`charsets/reverse`(ord)))) fi end: `charsets/idealint`:=proc(ps,f,ord) local qs,gb,zz,j; zz:=0; zz:='zz'; if 1 < printlevel then lprint(`ideal intersection starts at`,time()) fi; if type(f,{set,list}) then qs := [ 'zz*ps[j]' $ ('j' = 1 .. nops(ps)),'(1-zz)*f[j]' $ ('j' = 1 .. nops(f)) ] else qs := ['zz*ps[j]' $ ('j' = 1 .. nops(ps)),(1-zz)*f] fi; gb := Groebner['gbasis']( qs,plex(zz,'ord[nops(ord)-j+1]' $ ('j' = 1 .. nops(ord)))); qs := []; for j to nops(gb) do if {zz} minus indets(gb[j]) = {zz} then qs := [gb[j],op(qs)] fi od; if 1 < printlevel then lprint(`ideal intersection finished at`,time()) fi; qs end: `charsets/idealintR`:=proc(ps,f,ord) local qs,gb,zz,j; options remember; zz:=0; zz:='zz'; if 1 < printlevel then lprint(`ideal intersection starts at`,time()) fi; if type(f,{set,list}) then qs := [ 'zz*ps[j]' $ ('j' = 1 .. nops(ps)),'(1-zz)*f[j]' $ ('j' = 1 .. nops(f)) ] else qs := ['zz*ps[j]' $ ('j' = 1 .. nops(ps)),(1-zz)*f] fi; gb := Groebner['gbasis']( qs,plex(zz,'ord[nops(ord)-j+1]' $ ('j' = 1 .. nops(ord)))); qs := []; for j to nops(gb) do if {zz} minus indets(gb[j]) = {zz} then qs := [gb[j],op(qs)] fi od; if 1 < printlevel then lprint(`ideal intersection finished at`,time()) fi; qs end: `charsets/saturbasisR`:=proc(ps,js,ord) local qs,gb,zz,j; options remember; if js <> {} then qs := [op(ps),'cat(`charsets/@z`,j)*js[j]-1' $ ('j' = 1 .. nops(js))]; zz := ['cat(`charsets/@z`,(nops(js)-j+1))' $ ('j' = 1 .. nops(js))]; gb := Groebner['gbasis']( qs,plex(op(zz),'ord[nops(ord)-j+1]' $ ('j' = 1 .. nops(ord)))); qs := []; for j to nops(gb) do if {op(zz)} minus indets(gb[j]) = {op(zz)} then qs := [gb[j],op(qs)] fi od else qs := ps fi; qs end: `charsets/eicsTest`:=proc(ps,ord,medset) local qs,cs,is,iss,n,i,j,qhi,qsi,r,rr,factorset,mind,fset,ind,ts,den; ind:=0; ind:='ind'; if type(ps[1],{set,list}) then qhi := {ps} else qhi := {[ps,1]} fi; if medset = `charsets/basset` then mind := true else mind := false fi; qsi := {}; if 1 < printlevel then lprint(`radical membership test starts at`,time()) fi; for n from 0 while qhi <> {} and qsi = {} do qs := qhi[1][1]; if not mind then if n < 20 then cs := cat(`charsets/f`,(substring(medset,10 .. length(medset))))( qs,ord,[{},indets(qs)],'fset'); factorset := fset[1] else cs := cat(`charsets/`,(substring(medset, 10 .. length(medset))))(qs,ord); cs := `charsets/removecont`(cs,ord,'factorset') fi else if n < 20 then cs := `charsets/fcharseta`([op(qs)],ord,medset); factorset := op(2,cs[2]); cs := cs[1] else cs := `charsets/charseta`([op(qs)],ord,medset); cs := `charsets/removecont`(cs,ord,'factorset') fi; if 1 < printlevel then lprint(`characteristic set produced`); print(cs) fi fi; if 0 < `charsets/class`(cs[1],ord) then ts := `charsets/irras`(cs,ord,ind,'den'); if ts[2] = 0 then if not mind then if medset = `charsets/autored` then cs := `charsets/charseta`({op(qs),op(cs)}, ord,`charsets/charsetn`) elif not `charsets/subset`(cs,qs) then cs := `charsets/charseta`({op(qs),op(cs)}, ord,medset) fi; if 1 < printlevel then lprint(`characteristic set produced`); print(cs) fi fi; if 0 < `charsets/class`(cs[1],ord) then ts := `charsets/irras`(cs,ord,ind,'den'); if ts[2] = 0 then is := `charsets/initialset`(cs,ord) union `charsets/factorps`(factorset); if nops(cs) = nops(ord) then rr := `charsets/nopower`(qhi[1][2]) else rr := `charsets/nopower`( `charsets/prod`({qhi[1][2],op(is)})) fi; `charsets/premas`(rr,cs,ord); r := `charsets/simp`(%,cs,ord); if r <> 0 then if r = 1 then qsi := {cs,op(qsi)} else qsi := {[cs,`charsets/simpb`(r,rr)],op(qsi)} fi fi fi else is := `charsets/factorps`(factorset); ts := [1,0] fi fi; if ts[2] <> 0 then if 1 < ts[2] then is := `charsets/initialset`({op(1 .. ts[2]-1,cs)},ord) union `charsets/factorps`(factorset) else is := `charsets/factorps`(factorset) fi fi else is := `charsets/factorps`(factorset); ts := [1,0] fi; iss := {}; if nops(ord) <= nops(ps)+1 then for i in is do iss := {op(iss),[{i,op(qs)},qhi[1][2]]} od else for i to nops(is) do if i = 1 then 1 else product('is[j]','j' = 1 .. i-1) fi; iss := {op(iss),[{op(qs),is[i]},%*qhi[1][2]]} od fi; if ts[2] <> 0 and ts[1] <> {} then if not mind then if medset = `charsets/autored` then `charsets/charseta`({op(qs),op(cs)}, ord,`charsets/charsetn`) elif not `charsets/subset`(cs,qs) then `charsets/charseta`({op(qs),op(cs)},ord,medset) else cs fi; if % <> cs then if ts[2] = 1 then cs := qs else cs := {op(qs),op(1 .. ts[2]-1,cs)} fi fi fi; for i in ts[1] do iss := {op(iss),[{op(cs),i},`charsets/prod`({den,qhi[1][2],op(is)})]} od; if 0 < `charsets/class`(den,ord) and ts[2] < `charsets/class`(ts[1][1],ord) then iss := {op(iss),[{op(cs),den},qhi[1][2]]} fi fi; if 1 < nops(qhi) then qhi := {op(iss),op(2 .. nops(qhi),qhi)} else qhi := iss fi od; if qsi <> {} then op(qsi) else [] fi end: `charsets/lenord`:=proc(a,b) if length(b) < length(a) then true else false fi end: print("saving charset in the library found at:", libname[1]); savelib (charsets);