# globals: # List_Undets - list of all undeterminate variables # List_Consts - list of all constant variables # List_glConsts - list of all global constant variables, such as Pi... # sp_type - Type of scalarproduct for H-Bases (options: "D" for Differential and "M" for Monomial, standard is "D" without any option) # basis_type - Type of Gamma-Basis (options: "G" for grlex Groebner-Basis and "H" for H-Basis, standard is "H" without any option) # is_Basis - true or false, should be used carefully # # menge_Pspace - Basis for vector space, should be restarted after change of List_Undets # menge_Vspace - Basis for vector space, should be restarted after change of List_Undets or menge_Poly # menge_Vscal - scalarproducts for orthogonal space, should be restarted after change of List_Undets or menge_Poly # menge_Undets - list_coeff[i,j] for solving menge_Vscal # menge_Hpoly - Polynoms with coefficients from Vspace # ---------------------------------------------------------------------- with(LinearAlgebra): ElemIdeal:=proc(menge_Poly::list, menge_U::list)::list; local zaehl_i, menge_T; menge_T:=[]; for zaehl_i from 1 to nops(menge_Poly) do if (indets(menge_Poly[zaehl_i]) subset {op(menge_U)}) then menge_T:=[op(menge_T), menge_Poly[zaehl_i]]; fi; od; return(menge_T); end proc: IndSets:=proc(menge_Poly::list)::list; local menge_S, menge_T, zaehl_i; menge_T:=[]; menge_S:=combinat[powerset](List_Undets)[2..-1]; for zaehl_i from 1 to nops(menge_S) do if (ElemIdeal(menge_Poly, menge_S[zaehl_i])=[]) then menge_T:=[op(menge_T), menge_S[zaehl_i]]; fi; od; return(menge_T); end proc: dimId:=proc(menge_Poly::list)::integer; local menge_G, menge_I; menge_G:=gltSet(GBasis(menge_Poly)); menge_I:=IndSets(menge_G); if (nops(menge_I)>0) then return(max(seq(nops(menge_I[zaehl_i]), zaehl_i=1..nops(menge_I)))); else return(0); fi; end proc: dimIdeal:=proc(menge_Poly::list)::integer; global List_Undets,basis_type; local menge_H, zaehl_i, zaehl_j, max, menge_S, count; if (basis_type="G") then menge_H:=gltSet(GBasis(menge_Poly)); else menge_H:=hltSet(HBasis(menge_Poly)); fi; max:=0; menge_S:={solve({op(menge_H)}, {op(List_Undets)})}; for zaehl_i from 1 to nops(menge_S) do count:=0; for zaehl_j from 1 to nops(menge_S[zaehl_i]) do if (rhs(menge_S[zaehl_i][zaehl_j])=lhs(menge_S[zaehl_i][zaehl_j])) then count:=count+1; fi; od; if (count>max) then max:=count; fi; od; return(max); end proc: mdeg := proc(poly_p::polynom)::integer; description "returns degree resp. List_Undets"; global List_Undets; return(degree(poly_p,convert(List_Undets,set))); end proc: isMonom:=proc(poly_p::polynom)::bool; local monom::bool; monom:=true; if (mdeg(poly_p)>0) then if (convert(op(0,poly_p),string)="+") then monom:=false; fi; fi; return(monom); end proc: # ---------------------------------------------------------------------- setConsts:=proc(menge_C::list) global List_Consts, List_glConsts, constants; if (is(List_Consts,list)) then List_Consts:=menge_C; else List_glConsts:=[constants]; List_Consts:=menge_C; fi; constants:=op(List_glConsts), op(List_Consts): userinfo(3, 'information', cat("List of constants was set to ",convert(List_Consts,string),".")); end proc: setUndets:=proc(menge_U::list) global List_Undets, menge_Pspace; menge_Pspace:='menge_Pspace'; List_Undets:=menge_U; userinfo(3, 'information', cat("List of undets was set to ",convert(List_Undets,string),".")); end proc: setSP:=proc(art::string) global sp_type; sp_type:=art; userinfo(3, 'information', cat("Type of scalar product was set to ",sp_type,".")); end proc: setBasisType:=proc(art::string) global basis_type; basis_type:=art; userinfo(3, 'information', cat("Type of basis was set to ",basis_type,".")); end proc: setBasis:=proc(value) global is_basis; if (value=true or value=false) then is_basis:=value; userinfo(3, 'information', cat("Current ideal as basis was set to ",convert(is_basis,string),".")); else error "Value should be true or false. Nothing changed."; fi; end proc: # ---------------------------------------------------------------------- IdentityVector:=proc(int_len::integer, int_pos::integer)::Vector; description "returns Vector with length int_len and all zeros but position int_pos=1"; local vector_v::Vector; if (int_len<=0) then error "Vector length must be greater than 0."; fi; if (int_pos<=0) then error "Position index must be greater than 0."; fi; if (int_len0) then if not (isMonom(poly_f)) then for zaehl_i from 1 to nops(poly_f) do int_coeff:=gleadcoeff(op(zaehl_i,poly_f)); int_gcdc:=lcm(numer(int_gcdc),denom(int_coeff))/gcd(denom(int_gcdc),numer(int_coeff)); od; int_gcdc:=int_gcdc; else int_gcdc:=1/gleadcoeff(poly_f); fi; fi; if (int_gcdc=0) then int_gcdc:=1; fi; return(mcollect(poly_f*int_gcdc)); end proc: shortSet:=proc(menge_Poly::list)::list; description "shorten all polynoms in set"; local menge_H::list, zaehl_i::integer; menge_H:=[]; for zaehl_i from 1 to nops(menge_Poly) do menge_H:=[op(menge_H), shorten(menge_Poly[zaehl_i])]; od; return(menge_H); end proc: gleadterm := proc(poly_p::polynom)::polynom; description "gets leading term of a polynom"; local poly_g::polynom, poly_f::polynom; if (mdeg(poly_p)<=0) then poly_g:=poly_p; else poly_f:=mcollect(poly_p); if not (isMonom(poly_f)) then poly_g:=op(1,poly_f); else poly_g:=poly_f; fi; fi; return(poly_g); end proc: gleadcoeff := proc(poly_p::polynom)::integer; description "gets leading coefficient of a polynom"; local poly_f::polynom, int_coeff::integer, zaehl_i::integer; global List_Undets; int_coeff:=poly_p; if (mdeg(poly_p)>0) then poly_f:=gleadterm(poly_p); int_coeff:=1; if (convert(op(0,poly_f),string)="*") then for zaehl_i from 1 to nops(poly_f) while not (member(op(1,op(zaehl_i,poly_f)),List_Undets)) do int_coeff:=int_coeff*op(zaehl_i,poly_f); od; fi; fi; if (int_coeff=0) then int_coeff:=1; fi; return(int_coeff); end proc: gleadmon := proc(poly_p::polynom)::polynom; description "gets leading monom of a polynom"; local poly_f::polynom, poly_g::polynom; if (mdeg(poly_p)=0) then poly_g:=1; else poly_f:=gleadterm(poly_p); if (convert(op(0,poly_f),string)="*") then poly_g:=poly_f/gleadcoeff(poly_f); else poly_g:=poly_f; fi; fi; return(mcollect(poly_g)); end proc: gltSet:=proc(menge_Poly::list)::list; description "returns set of leading terms"; local menge_H::list, zaehl_i::integer; menge_H:=[]; for zaehl_i from 1 to nops(menge_Poly) do menge_H:=[op(menge_H),gleadterm(menge_Poly[zaehl_i])]; od; return(menge_H); end proc: SyzProd :=proc(syz_Poly::Vector, menge_Poly::list)::polynom; description "product between two vectors"; if (op(1,syz_Poly)<>nops(menge_Poly)) then error "Size of syzygy vector (%1) and number of polynomials (%2) doesn't match.", op(1,syz_Poly), nops(menge_Poly); fi; return(mcollect(add(syz_Poly[zaehl_i]*menge_Poly[zaehl_i], zaehl_i=1..nops(menge_Poly)))); end proc: gSyzSet :=proc(poly_f::polynom, poly_g::polynom)::Vector; description "returns Vector with two syzygy polynoms"; local poly_u::polynom, poly_v::polynom, poly_temp::polynom, vector_v::Vector; vector_v:=Vector(2); if (poly_f<>0 and poly_g<>0) then poly_u := gleadterm(poly_f); poly_v := gleadterm(poly_g); poly_temp := lcm(numer(poly_u),numer(poly_v))/lcm(denom(poly_u),denom(poly_v)); poly_u := poly_temp/poly_u; poly_v := poly_temp/poly_v; poly_v := -1*poly_v/gleadcoeff(poly_u); poly_u := gleadmon(poly_u); vector_v[1]:=poly_u; vector_v[2]:=poly_v; fi; return(vector_v); end proc: gSyzPoly :=proc(poly_f::polynom, poly_g::polynom)::polynom; description "returns remainder after multiplication with syzygy polynoms"; return(SyzProd(gSyzSet(poly_f,poly_g),[poly_f,poly_g])); end proc: gSyzSets :=proc(menge_Poly::list)::list; description "returns syzygy set"; local zaehl_i::integer, zaehl_j::integer, menge_S::list, vector_v::Vector, vector_temp::Vector; vector_v:=Vector(2); vector_temp:=Vector(nops(menge_Poly)); menge_S:=[]; for zaehl_i from 1 to nops(menge_Poly) do for zaehl_j from zaehl_i+1 to nops(menge_Poly) do vector_v:=gSyzSet(menge_Poly[zaehl_i],menge_Poly[zaehl_j]); vector_temp:=Vector(nops(menge_Poly)); vector_temp[zaehl_i]:=vector_v[1]; vector_temp[zaehl_j]:=vector_v[2]; menge_S:=[op(menge_S),vector_temp]; od; od; return(menge_S); end proc: gDivAlg := proc(menge_Poly::list,poly_g::polynom)::list; description "division algorithm returning coefficients and remainder"; local poly_t::polynom, menge_C::list, zaehl_k::integer, dva::bool; # # print(menge_Poly, poly_g); for zaehl_k from 1 to nops(menge_Poly)+1 do menge_C[zaehl_k]:=0; od; poly_t:=poly_g; while (poly_t <> 0) do zaehl_k:=1; dva:=false; # # print(zaehl_k, gltSet(menge_Poly), gleadmon(poly_t)); while (zaehl_k <= nops(menge_Poly) and dva=false) do # # print("polynom?", gleadmon(poly_t)/gleadmon(menge_Poly[zaehl_k])); if (type(gleadmon(poly_t)/gleadmon(menge_Poly[zaehl_k]),polynom)) then menge_C[zaehl_k]:=menge_C[zaehl_k]+gleadterm(poly_t)/gleadterm(menge_Poly[zaehl_k]); # # print(menge_C); poly_t:=mcollect(simplify(poly_t-(gleadterm(poly_t)/gleadterm(menge_Poly[zaehl_k]))*menge_Poly[zaehl_k])); dva:=true; else zaehl_k:=zaehl_k+1; fi; od; if not (dva) then menge_C[nops(menge_Poly)+1]:=menge_C[nops(menge_Poly)+1]+gleadterm(poly_t); poly_t:=mcollect(poly_t-gleadterm(poly_t)); fi; od; for zaehl_k from 1 to nops(menge_Poly)+1 do menge_C[zaehl_k]:=mcollect(menge_C[zaehl_k]); od; return(convert(menge_C,list)); end proc: gcfs := proc(menge_Poly::list,poly_g::polynom)::list; description "returns coeffs of division algorithm"; local menge_H::list; menge_H:=gDivAlg(menge_Poly,poly_g); return(menge_H[1..nops(menge_H)-1]); end proc: grmd := proc(menge_Poly::list,poly_g::polynom)::polynom; description "returns remainder of division algorithm"; global List_Undets; local menge_H::list; menge_H:=gDivAlg(menge_Poly,poly_g); return(menge_H[nops(menge_H)]); # return(Groebner[normalf](poly_g, menge_Poly, tdeg(op(List_Undets)))); end proc: canReduce:=proc(menge_Poly::list, int_pos::integer)::Vector; description "returns index of polynom that could be reduced"; local menge_H::list, zaehl_i::integer, vector_r::Vector, vector_v::Vector, menge_I::list; if (nops(menge_Poly)>0) then menge_H:=menge_Poly; zaehl_i:=int_pos+1; vector_r:=Vector(nops(menge_H)); while (10) then vector_r[1]:=zaehl_i; fi; od; else vector_r:=Vector(1); fi; return(vector_r); end proc: reduceSet:=proc(menge_Poly::list)::list; description "reduces a set of polynoms rsp. G-Basis"; local menge_H::list, vector_r::Vector, vector_v::Vector, int_pos::integer, poly_r::polynom; menge_H:=menge_Poly; vector_r:=canReduce(menge_H, nops(menge_H)); int_pos:=vector_r[1]; while (int_pos>0) do vector_v:=convert([op(menge_H[1..int_pos-1]),op(menge_H[int_pos+1..nops(menge_H)])],Vector); poly_r:=mcollect(menge_H[int_pos]-Transpose(vector_v).vector_r[2..nops(menge_H)]); if (poly_r<>0) then menge_H:=[op(menge_H[1..int_pos-1]),poly_r,op(menge_H[int_pos+1..nops(menge_H)])]; else menge_H:=[op(menge_H[1..int_pos-1]),op(menge_H[int_pos+1..nops(menge_H)])]; fi; int_pos:=int_pos-1; if (int_pos=0) then int_pos:=nops(menge_H); fi; vector_r:=canReduce(menge_H, int_pos); int_pos:=vector_r[1]; od; return(menge_H); end proc: IGBMatrix:=proc(menge_Poly::list)::Matrix; description "returns Matrix for G = F*M, F set of polynoms and G as reduced G-Basis"; local menge_G::list, poly_r::polynom, zaehl_i::integer, zaehl_j::integer, vector_v::Vector, menge_N::list, vector_r::Vector, menge_H::list, vector_t::Vector, menge_M::list, poly_h::polynom, vector_temp::Vector, int_pos::integer, matrix_M::Matrix, change::bool, temp, menge_R; menge_G:=menge_Poly; menge_M:=[seq(IdentityVector(nops(menge_Poly),zaehl_i), zaehl_i=1..nops(menge_Poly))]; # print ("ich reduziere", nops(menge_G), time()); vector_r:=canReduce(menge_G, nops(menge_G)); int_pos:=vector_r[1]; while (int_pos>0) do vector_v:=convert([op(menge_G[1..int_pos-1]),op(menge_G[int_pos+1..nops(menge_G)])],Vector); matrix_M:=convert([op(menge_M[1..int_pos-1]),op(menge_M[int_pos+1..nops(menge_M)])],Matrix); poly_r:=mcollect(menge_G[int_pos]-Transpose(vector_v).vector_r[2..nops(menge_G)]); # print (int_pos, "reduziertes Polynom"); if (poly_r<>0) then menge_G:=[op(menge_G[1..int_pos-1]),poly_r,op(menge_G[int_pos+1..nops(menge_G)])]; vector_temp:=collectVector(menge_M[int_pos]-matrix_M.vector_r[2..nops(menge_G)]); menge_M:=[op(menge_M[1..int_pos-1]),vector_temp,op(menge_M[int_pos+1..nops(menge_M)])]; else menge_G:=[op(menge_G[1..int_pos-1]),op(menge_G[int_pos+1..nops(menge_G)])]; menge_M:=[op(menge_M[1..int_pos-1]),op(menge_M[int_pos+1..nops(menge_M)])]; fi; int_pos:=int_pos-1; if (int_pos=0) then int_pos:=nops(menge_G); fi; vector_r:=canReduce(menge_G,int_pos); int_pos:=vector_r[1]; od; # print ("reduzieren beendet", nops(menge_G), time()); while (menge_G<>menge_H) do menge_H:=menge_G; menge_N:=[]; for zaehl_i from 1 to nops(menge_H) do for zaehl_j from zaehl_i+1 to nops(menge_H) do poly_h:=gSyzPoly(menge_H[zaehl_i],menge_H[zaehl_j]); # print(zaehl_i, zaehl_j, "Syz-Polynom", poly_h); menge_R:=gDivAlg(menge_H, poly_h); poly_r:=menge_R[-1]; # print(zaehl_i, zaehl_j, "Rest", poly_r); if (poly_r <> 0) then vector_v:=gSyzSet(menge_H[zaehl_i],menge_H[zaehl_j]); vector_t:=convert(menge_R[1..-2],Vector); vector_temp:=Vector(nops(menge_H)); vector_temp[zaehl_i]:=vector_v[1]; vector_temp[zaehl_j]:=vector_v[2]; vector_temp:=vector_temp-vector_t; menge_N:=[op(menge_N),collectVector(convert(menge_M,Matrix).vector_temp)]; menge_G:=[op(menge_G),poly_r]; fi; od; od; menge_M:=[op(menge_M),op(menge_N)]; # print ("ich reduziere", nops(menge_G), time()); vector_r:=canReduce(menge_G, nops(menge_G)); int_pos:=vector_r[1]; while (int_pos>0) do vector_v:=convert([op(menge_G[1..int_pos-1]),op(menge_G[int_pos+1..nops(menge_G)])],Vector); matrix_M:=convert([op(menge_M[1..int_pos-1]),op(menge_M[int_pos+1..nops(menge_M)])],Matrix); poly_r:=mcollect(menge_G[int_pos]-Transpose(vector_v).vector_r[2..nops(menge_G)]); # print (int_pos, "reduziertes Polynom"); if (poly_r<>0) then menge_G:=[op(menge_G[1..int_pos-1]),poly_r,op(menge_G[int_pos+1..nops(menge_G)])]; vector_temp:=collectVector(menge_M[int_pos]-matrix_M.vector_r[2..nops(menge_G)]); menge_M:=[op(menge_M[1..int_pos-1]),vector_temp,op(menge_M[int_pos+1..nops(menge_M)])]; else menge_G:=[op(menge_G[1..int_pos-1]),op(menge_G[int_pos+1..nops(menge_G)])]; menge_M:=[op(menge_M[1..int_pos-1]),op(menge_M[int_pos+1..nops(menge_M)])]; fi; int_pos:=int_pos-1; if (int_pos=0) then int_pos:=nops(menge_G); fi; vector_r:=canReduce(menge_G, int_pos); int_pos:=vector_r[1]; od; # print ("reduzieren beendet", nops(menge_G), time()); # print("vorlaeufige G-Basis", menge_G); od; return(convert(menge_M,Matrix)); end proc: GBasis:=proc(menge_Poly::list)::list; description "returns reduced Groebner-Basis for set of polynoms"; global List_Undets, is_basis; if not (is_basis) then return(shortSet(convert(collectVector(Transpose(convert(menge_Poly,Vector)).(IGBMatrix(menge_Poly))),list))); else return(menge_Poly); fi; end proc: GBMatrix:=proc(menge_Poly::list)::Matrix; description "returns Matrix for F = G*M, F set of polynoms and G as reduced G-Basis"; local basis::list, menge_M::list, zaehl_i::integer; menge_M:=[]; basis:=GBasis(menge_Poly); for zaehl_i from 1 to nops(menge_Poly) do menge_M:=[op(menge_M),convert(gcfs(basis,menge_Poly[zaehl_i]),Vector)]; od; return(convert(menge_M,Matrix)); end proc: gSMatrix:=proc(menge_Poly::list)::Matrix; local basis::list, menge_S::list, menge_R::list, poly_h::polynom, vector_temp::Vector, zaehl_i::integer, zaehl_j::integer, zaehl_k::integer; zaehl_k:=0; basis:=GBasis(menge_Poly); vector_temp:=Vector(nops(basis)); menge_S:=gSyzSets(basis); menge_R:=[]; for zaehl_i from 1 to nops(basis) do for zaehl_j from zaehl_i+1 to nops(basis) do zaehl_k:=zaehl_k+1; poly_h:=gSyzPoly(basis[zaehl_i],basis[zaehl_j]); vector_temp:=menge_S[zaehl_k]-convert(gcfs(basis,poly_h),Vector); menge_R:=[op(menge_R),vector_temp]; od; od; return(Transpose(convert(menge_R, Matrix))); end proc: aSyzMod:=proc(menge_Poly::list)::list; description "returns list G with for F*G = 0 with F set of polynoms (basis of syzygy modul)"; global is_basis, infolevel; local menge_G::list, zaehl_i::integer, menge_M::list, vector_temp::Vector, int_pos::integer, matrix_M::Matrix, matrix_N::Matrix, matrix_Q::Matrix, matrix_R::Matrix, value; matrix_M:=IGBMatrix(menge_Poly); menge_G:=shortSet(convert(collectVector(Transpose(convert(menge_Poly,Vector)).(IGBMatrix(menge_Poly))),list)); # print("G-Basis", menge_G); value:=is_basis; infolevel[information]:=2; setBasis(true); menge_M:=[]; for zaehl_i from 1 to nops(menge_Poly) do menge_M:=[op(menge_M),convert(gcfs(menge_G,menge_Poly[zaehl_i]),Vector)]; od; matrix_N:=convert(menge_M,Matrix); matrix_R:=gSMatrix(menge_G); # print(matrix_M); # print(matrix_N); # print(matrix_R); int_pos:=op(1,[op(1,matrix_R)]); matrix_Q:=IdentityMatrix(nops(menge_Poly))-Transpose(matrix_N).Transpose(matrix_M); matrix_Q:=Matrix(nops(menge_Poly)+int_pos,nops(menge_Poly)); matrix_Q[1..nops(menge_Poly),1..nops(menge_Poly)]:=collectMatrix(simplify(IdentityMatrix(nops(menge_Poly))-Transpose(matrix_N).Transpose(matrix_M))); matrix_Q[nops(menge_Poly)+1..nops(menge_Poly)+int_pos,1..nops(menge_Poly)]:=collectMatrix(matrix_R.Transpose(matrix_M)); # print(matrix_Q); matrix_Q:=Transpose(matrix_Q); int_pos:=op(2,[op(1,matrix_Q)]); menge_M:=[]; for zaehl_i from 1 to int_pos do vector_temp:=Column(matrix_Q,zaehl_i); if (Norm(vector_temp)<>0) then menge_M:=[op(menge_M),vector_temp]; fi; od; # # print(menge_M); setBasis(value); infolevel[information]:=3; return(menge_M); end proc: # ----------------------------------------------------------------------------------- mip := proc (mi_a::list)::polynom; description "returns polynom of a multiindex"; global List_Undets; if (nops(mi_a) <> nops(List_Undets)) then error "Size of array (%1) and number of undeterminates (%2) does not match.",nops(mi_a), nops(List_Undets); fi; return(mcollect(mul(List_Undets[zaehl_i]^mi_a[zaehl_i], zaehl_i=1..nops(List_Undets)))); end proc: coeffmi := proc(poly_p::polynom, mi_a::list)::integer; description "gets the coefficient of a polynom with multiindex a"; local poly_f::polynom, poly_t::polynom, zaehl_i::integer, int_coeff::integer, poly_a::polynom; global List_Undets; if (nops(mi_a) <> nops(List_Undets)) then error "Size of array (%1) and number of undeterminates (%2) does not match.",nops(mi_a), nops(List_Undets); fi; int_coeff:=0; poly_a:=mip(mi_a); if (mdeg(poly_p)<=0) then if (poly_a=1) then int_coeff:=poly_p; fi; else poly_f:=mcollect(poly_p); if not (isMonom(poly_f)) then for zaehl_i from 1 to nops(poly_f) do poly_t:=op(zaehl_i,poly_f); if (gleadmon(poly_t) = poly_a) then int_coeff:=gleadcoeff(poly_t); break; fi; od; else if (gleadmon(poly_f) = poly_a) then int_coeff:=gleadcoeff(poly_f); fi; fi; fi; return(int_coeff); end proc: pmi := proc (poly_p::polynom)::list; description "returns multiindex of monom"; local poly_f::polynom, poly_g::polynom, mi_a::list, zaehl_j::integer, int_pos; global List_Undets; if not (isMonom(poly_p)) then error "Could not calculate multiindex of polynomial %1.",poly_p; fi; mi_a:=[seq(0, zaehl_j=1..nops(List_Undets))]; if (mdeg(poly_p)>0) then poly_f:=gleadmon(poly_p); if (convert(op(0,poly_f),string)="^") then if (member(op(1,poly_f), List_Undets, 'int_pos')) then mi_a[int_pos]:=op(2,poly_f); fi; else for zaehl_j from 1 to nops(poly_f) do poly_g:=op(zaehl_j,poly_f); if (member(op(1,poly_g), List_Undets, 'int_pos')) then if (convert(op(0,poly_g),string)="^") then mi_a[int_pos]:=op(2,poly_g); else mi_a[int_pos]:=1; fi; fi; od; fi; fi; return(mi_a); end proc: facmi := proc(mi_a::list)::integer; description "factorial of a multiindex as list"; local zaehl_i::integer, temp::integer; temp:=mul(factorial(mi_a[zaehl_i]), zaehl_i=1..nops(mi_a)); return (temp); end proc: scalarproduct := proc(poly_p::polynom,poly_q::polynom)::integer; description "scalar product between two polynoms"; local mi_a::list, zaehl_i::integer, int_scalprod::integer, poly_f::polynom, poly_g::polynom, poly_t::polynom; global sp_type; int_scalprod:=0; poly_f:=mcollect(poly_p); poly_g:=mcollect(poly_q); if (mdeg(poly_f)<=0 and mdeg(poly_g)<=0) then int_scalprod:=poly_f*poly_g; elif (mdeg(poly_f)<=0 or mdeg(poly_g)<=0) then mi_a:=pmi(0); # constant term int_scalprod:=coeffmi(poly_f,mi_a)*coeffmi(poly_g,mi_a); else if (isMonom(poly_f) or isMonom(poly_g)) then if (isMonom(poly_f)) then mi_a:=pmi(poly_f); else mi_a:=pmi(poly_g); fi; if (sp_type="M") then int_scalprod:=coeffmi(poly_f,mi_a)*coeffmi(poly_g,mi_a); else int_scalprod:=coeffmi(poly_f,mi_a)*coeffmi(poly_g,mi_a)*facmi(mi_a); fi; else # f and g both polynoms if (nops(poly_f) > nops(poly_g)) then poly_t:=poly_f; poly_f:=poly_g; poly_g:=poly_t; # shorter polynom at first fi; for zaehl_i from 1 to nops(poly_f) do poly_t:=op(zaehl_i,poly_f); mi_a:=pmi(poly_t); if (sp_type="M") then int_scalprod:=int_scalprod+coeffmi(poly_f,mi_a)*coeffmi(poly_g,mi_a); else int_scalprod:=int_scalprod+coeffmi(poly_f,mi_a)*coeffmi(poly_g,mi_a)*facmi(mi_a); fi; od; fi; fi; return(int_scalprod); end proc: hleadterm := proc(poly_p::polynom)::polynom; description "gets the homogeneous leading monom of a polynom"; local poly_f::polynom, poly_g::polynom, zaehl_i::integer; poly_f:=mcollect(poly_p); poly_g:=0; if not (isMonom(poly_f)) then for zaehl_i from 1 to nops(poly_f) while mdeg(poly_f)=mdeg(op(zaehl_i,poly_f)) do poly_g:=poly_g+op(zaehl_i,poly_f); od; else poly_g:=poly_f; fi; return(poly_g); end proc: Pspace := proc(int_d::integer)::list; description "builds a basis for vector space PId"; global menge_Pspace, List_Undets; local menge_G::list, zaehl_j::integer, zaehl_k::integer, poly_p::polynom, menge_H::list; menge_G:=[]; if (int_d>=0) then if (is(menge_Pspace[int_d],list)=false) then if (int_d=0) then menge_G:=[1]; else menge_G:=Pspace(int_d-1); menge_H:=[]; for zaehl_j from 1 to nops(List_Undets) do for zaehl_k from 1 to nops(menge_G) do poly_p:=menge_G[zaehl_k]*List_Undets[zaehl_j]; if not (member(poly_p, menge_H)) then menge_H:=[op(menge_H), poly_p]; fi; od; od; menge_G:=menge_H; menge_Pspace[int_d]:=menge_G; fi; else menge_G:=menge_Pspace[int_d]; fi; fi; return(menge_G); end proc: Vspace := proc(int_d::integer,menge_Poly::list)::list; description "builds a generating set for a subspace of PId from polynoms"; local menge_G::list, zaehl_i::integer, zaehl_j::integer, int_max::integer, int_min::integer, int_temp::integer, poly_temp::polynom, poly_temp2::polynom, menge_H::list, menge_Temp::list, temp_deg; menge_Temp:=[]; for zaehl_i from 1 to nops(menge_Poly) do poly_temp:=shorten(hleadterm(menge_Poly[zaehl_i])); if (poly_temp<>0 and not member(poly_temp, menge_Temp)) then menge_Temp:=[op(menge_Temp), poly_temp]; fi; od; temp_deg:=seq(mdeg(menge_Temp[zaehl_i]), zaehl_i=1..nops(menge_Temp)); int_max:=max(temp_deg); int_min:=min(temp_deg); menge_G:=[]; for zaehl_i from int_min to int_max do menge_H[zaehl_i]:=Pspace(int_d-zaehl_i); od; for zaehl_i from 1 to nops(menge_Temp) do poly_temp:=menge_Temp[zaehl_i]; int_temp:=mdeg(poly_temp); for zaehl_j from 1 to nops(menge_H[int_temp]) do poly_temp2:=mcollect(poly_temp*menge_H[int_temp][zaehl_j]); if not (member(poly_temp2, menge_G)) then menge_G:=[op(menge_G), poly_temp2]; fi; od; od; return(LinSpace(shortSet(menge_G))); end proc: Wspace := proc(int_d::integer,menge_Poly::list)::list; description "builds a generating set for orthogonal subspace of PId from polynoms"; local menge_V::list, menge_P::list, menge_G::list, poly_temp::polynom, list_coeff, menge_S::set, zaehl_i::integer, zaehl_j::integer, menge_Temp::set, menge_U::set, matrix_M::matrix, vector_temp::Vector; menge_P:=Pspace(int_d); menge_V:=Vspace(int_d,menge_Poly); menge_G:=[]; menge_Temp:={}; menge_U:=[]; # print(menge_V, menge_P); for zaehl_i from 1 to nops(menge_P) do menge_U:=[op(menge_U), list_coeff[zaehl_i]]; od; # print(menge_U); poly_temp:=add(list_coeff[zaehl_i]*menge_P[zaehl_i], zaehl_i=1..nops(menge_P)); # print(poly_temp); for zaehl_i from 1 to nops(menge_V) do menge_Temp:=menge_Temp union {scalarproduct(menge_V[zaehl_i],poly_temp)=0}; od; # print(menge_Temp); (matrix_M, vector_temp):=GenerateMatrix(menge_Temp, menge_U); # print(matrix_M, vector_temp); menge_S:=NullSpace(matrix_M); # print(menge_S); for zaehl_j from 1 to nops(menge_S) do menge_G:=[op(menge_G), Transpose(convert(menge_P, Vector)).menge_S[zaehl_j]]; od; return(LinSpace(shortSet(menge_G))); end proc: LinSpace:=proc(menge_Poly::list)::list; local menge_LT, menge_V, int_length, zaehl_i, poly_temp, vector_temp, zaehl_j, poly_lt, int_coeff, menge_B, menge_U; # menge_LT:=[]; # menge_V:=[]; # print("Konvertiere",nops(menge_Poly),"Polynome in Vektoren"); int_length:=0; for zaehl_i from 1 to nops(menge_Poly) do poly_temp:=menge_Poly[zaehl_i]; vector_temp:='vector_temp'; # vector_temp:=[]; while (poly_temp<>0) do poly_lt:=gleadmon(poly_temp); int_coeff:=gleadcoeff(poly_temp); # # print(int_coeff, poly_lt); poly_temp:=mcollect(poly_temp-int_coeff*poly_lt); if (member(poly_lt, convert(menge_LT,list), 'int_pos')) then # # print("member"); if (convert(op(vector_temp),string)="vector_temp") then for zaehl_j from 1 to int_pos do # vector_temp:=[op(vector_temp), 0]; vector_temp[zaehl_j]:=0; od; else for zaehl_j from nops(op(op(vector_temp)))+1 to int_pos do # vector_temp:=[op(vector_temp), 0]; vector_temp[zaehl_j]:=0; od; fi; vector_temp[int_pos]:=int_coeff; else int_length:=int_length+1; # # print(int_length, nops(vector_temp)); if (convert(op(vector_temp),string)="vector_temp") then for zaehl_j from 1 to int_length do # vector_temp:=[op(vector_temp), 0]; vector_temp[zaehl_j]:=0; od; else for zaehl_j from nops(op(op(vector_temp)))+1 to int_length do # vector_temp:=[op(vector_temp), 0]; vector_temp[zaehl_j]:=0; od; fi; vector_temp[int_length]:=int_coeff; # menge_LT:=[op(menge_LT), poly_lt]; menge_LT[int_length]:=poly_lt; fi; od; # menge_V:=[op(menge_V), vector_temp]; menge_V[zaehl_i]:=op(vector_temp); od; # # print(menge_V); for zaehl_i from 1 to nops(menge_Poly) do vector_temp:=menge_V[zaehl_i]; # # print(vector_temp); for zaehl_j from nops(convert(vector_temp,list))+1 to int_length do vector_temp[zaehl_j]:=0; od; # # print(vector_temp); menge_V[zaehl_i]:=convert(convert(vector_temp,list), Vector); od; menge_V:=convert(menge_V,list); menge_LT:=convert(convert(menge_LT,list), Vector); # print(menge_V); # print(menge_LT); # print("Berechne Basis von", nops(menge_V), "Vektoren der Laenge", int_length); if (type(menge_V[1],Vector)) then # # print(time()); # menge_B:=GramSchmidt(menge_V); # # print(time()); menge_B:=Basis(menge_V); # # print(time()); else menge_B:=[]; fi; # print(menge_B); # print("Konvertiere",nops(menge_B),"Vektoren in Polynome"); if (nops(menge_B)<>0) then for zaehl_i from 1 to nops(menge_B) do menge_U[zaehl_i]:=Transpose(menge_B[zaehl_i]).menge_LT; od; menge_U:=convert(menge_U, list); else menge_U:=[]; fi; return(menge_U); end proc: hltSet:=proc(menge_Poly::list)::list; description "returns set of homogeneous leading terms"; local menge_H::list, zaehl_i::integer; menge_H:=[]; for zaehl_i from 1 to nops(menge_Poly) do menge_H:=[op(menge_H),shorten(hleadterm(menge_Poly[zaehl_i]))]; od; return(menge_H); end proc: hcoeffPoly :=proc (menge_Poly::list,poly_g::polynom)::list; description "returns polynoms h_i for lm(g) = sum (f_i*h_i) + r with (f_i*h_i) in V and r in W"; local menge_H::list, int_gamma::integer, zaehl_i::integer, menge_P::list, poly_r::polynom, menge_Temp::set, menge_V::list, menge_S1::set, menge_S2::set, int_temp::integer, zaehl_k::integer, zaehl_j::integer, list_coeff, menge_VS, menge_U, vector_t, matrix_M, vector_v, menge_ind; global List_Undets, menge_Vspace, menge_Vscal, menge_Undets, menge_Hpoly; # # print("hcoeffPoly", poly_g); menge_Temp:=convert(menge_Poly,set); if (poly_g<>0) then menge_S1:={}; menge_S2:={}; menge_U:=[]; menge_H:=[]; int_gamma:=mdeg(poly_g); # # print("Grad", int_gamma); if (int_gamma<0) then int_gamma:=-1; fi; if (op(0,menge_Vspace[int_gamma])<>list) then for zaehl_i from 1 to nops(menge_Poly) do if (menge_Poly[zaehl_i]<>0) then menge_P:=Pspace(int_gamma-mdeg(menge_Poly[zaehl_i])); else menge_P:=Pspace(-1); fi; menge_H:=[op(menge_H), add(list_coeff[zaehl_i,zaehl_j]*menge_P[zaehl_j], zaehl_j=1..nops(menge_P))]; for zaehl_k from 1 to nops(menge_P) do menge_U:=[op(menge_U), list_coeff[zaehl_i,zaehl_k]]; od; od; menge_Hpoly[int_gamma]:=menge_H; menge_Undets[int_gamma]:=menge_U; # # print("Berechne Vspace"); menge_Vspace[int_gamma]:=Vspace(int_gamma,menge_Poly); menge_V:=menge_Vspace[int_gamma]; poly_r:=-add(menge_H[zaehl_j]*hleadterm(menge_Poly[zaehl_j]), zaehl_j=1..nops(menge_Poly)); # # print(poly_r); # print("Bilde Skalarprodukte global" , nops(menge_V)); menge_Vscal[int_gamma]:=[]; for zaehl_i from 1 to nops(menge_V) do menge_Vscal[int_gamma]:=[op(menge_Vscal[int_gamma]), scalarproduct(menge_V[zaehl_i],poly_r)]; if (zaehl_i mod 100=0) then # print("aktuell", zaehl_i); fi; od; menge_VS:=menge_Vscal[int_gamma]; else menge_V:=menge_Vspace[int_gamma]; menge_VS:=menge_Vscal[int_gamma]; menge_U:=menge_Undets[int_gamma]; menge_H:=menge_Hpoly[int_gamma]; fi; # # print("grad", int_gamma); # # print("V", menge_V); # # print("VS", menge_VS); poly_r:=hleadterm(poly_g); # # print("Leitterm", poly_r, poly_g); # print("Bilde Skalarprodukte lokal" , nops(menge_V)); menge_Temp:=[]; for zaehl_i from 1 to nops(menge_V) do menge_Temp:=[op(menge_Temp), scalarproduct(menge_V[zaehl_i],poly_r)+menge_VS[zaehl_i]=0]; if (zaehl_i mod 1000=0) then # print("aktuell", zaehl_i); fi; od; # # print("T", menge_Temp); # # print("Koeffizienten", menge_U); # # print ("Loese das System"); (matrix_M,vector_v):=GenerateMatrix(expand(menge_Temp),menge_U); # # print(matrix_M,vector_v); vector_t:=LinearSolve(matrix_M,vector_v); menge_ind:=convert(indets(vector_t),list); menge_ind:=[seq(menge_ind[zaehl_i]=0, zaehl_i=1..nops(menge_ind))]; # # print(vector_t, subs(menge_ind,vector_t)); if (Norm(subs(menge_ind,vector_t))=0) then for zaehl_i from 1 to nops(menge_ind) do menge_ind[zaehl_i]:=lhs(menge_ind[zaehl_i])=1; od; fi; vector_t:=subs(menge_ind,vector_t); # # print(vector_t); for zaehl_i from 1 to nops(menge_U) do menge_S2:=menge_S2 union {menge_U[zaehl_i]=vector_t[zaehl_i]}; od; menge_H:=subs(menge_S2,menge_H); fi; # # print(menge_H); return(menge_H); end proc: hDivAlg := proc(menge_Poly::list,poly_g::polynom)::list; description "division algorithm returning coefficients and remainder"; local menge_H::list, poly_r::polynom, poly_t::polynom, zaehl_i::integer, menge_C::list, zaehl_j::integer; for zaehl_i from 1 to nops(menge_Poly)+1 do menge_C[zaehl_i]:=0; od; poly_t:=poly_g; while (poly_t <> 0) do menge_H:=hcoeffPoly(menge_Poly,poly_t); for zaehl_i from 1 to nops(menge_H) do menge_C[zaehl_i]:=menge_C[zaehl_i]+menge_H[zaehl_i]; od; poly_r:=hleadterm(poly_t)-add(menge_H[zaehl_j]*hleadterm(menge_Poly[zaehl_j]), zaehl_j=1..nops(menge_Poly)); menge_C[nops(menge_Poly)+1]:=menge_C[nops(menge_Poly)+1]+poly_r; poly_t:=mcollect(poly_t-add(menge_H[zaehl_j]*menge_Poly[zaehl_j], zaehl_j=1..nops(menge_Poly))-poly_r); od; for zaehl_i from 1 to nops(menge_Poly)+1 do menge_C[zaehl_i]:=mcollect(menge_C[zaehl_i]); od; return(convert(menge_C,list)); end proc: hcfs := proc(menge_Poly::list,poly_g::polynom)::list; description "returns coeffs of division algorithm"; local menge_H::list; menge_H:=hDivAlg(menge_Poly,poly_g); return(menge_H[1..nops(menge_H)-1]); end proc: hrmd := proc(menge_Poly::list,poly_g::polynom)::polynom; description "returns remainder of division algorithm"; local menge_H::list; menge_H:=hDivAlg(menge_Poly,poly_g); return(menge_H[nops(menge_H)]); end proc: isGBasis:=proc(menge_Poly::list)::bool; description "checks if set is already a G-Basis <=> s*F -> 0 mod F for each s in Syz(LT(F))"; local zaehl_i, zaehl_j, gb::bool, poly_r, poly_h; gb:=true; # print ("Teste auf G-Basis"); for zaehl_i from 1 to nops(menge_Poly) do for zaehl_j from zaehl_i+1 to nops(menge_Poly) do poly_h:=gSyzPoly(menge_Poly[zaehl_i],menge_Poly[zaehl_j]); # print("Syz-Polynom Nr.", zaehl_i,zaehl_j,poly_h); # print("Div", gDivAlg(menge_Poly,poly_h)); poly_r:=grmd(menge_Poly,poly_h); # print("Rest", poly_r); if (poly_r <> 0) then gb:=false; break; fi; od; if not (gb) then # print ("keine G-Basis"); break; fi; od; return(gb); end proc: isHBasis:=proc(menge_Poly::list, full)::bool; description "checks if set is already a H-Basis <=> s*F -> 0 mod F for each s in Syz(LT(F))"; local poly_h::polynom, menge_S::list, zaehl_i::integer, hb::bool, zaehl_j::integer, allMonomials::bool, menge_H, menge_L, onlyZero::bool; hb:=true; allMonomials:=true; menge_H:=hltSet(menge_Poly); # # if menge_Poly G-Basis and all homogeneous leading terms are monomials then menge_Poly is already a H-Basis # print ("Teste auf Monome"); for zaehl_i from 1 to nops(menge_H) do if not (isMonom(menge_H[zaehl_i])) then # print ("keine Monome"); allMonomials:=false; break; fi; od; if (allMonomials=false or isGBasis(menge_Poly)=false) then onlyZero:=true; # print ("Teste auf Nullstellen"); menge_L:={allvalues({op({solve(convert(menge_H,set), convert(List_Undets,set))})})}; # print(menge_L, nops(menge_L)); while (nops(menge_L)=1 and convert(op(0,op(1,menge_L)),string)="set") do menge_L:=op(menge_L); od; if (nops(menge_L)>1 and convert(op(0,op(1,menge_L)),string)="set") then onlyZero:=false; else # print(menge_L, op(menge_L)); for zaehl_i from 1 to nops(menge_L) do if (rhs(menge_L[zaehl_i])<>0) then onlyZero:=false; break; fi; od; fi; if (not full) then # print ("Keine H-Basis bisher"); fi; # print("only Zero", onlyZero, "full", full, "allMonomials", allMonomials); if (not onlyZero) then if (full) then # print ("Berechne Basis der Syz-Modul"); menge_S:=aSyzMod(menge_H); # print ("Anzahl Basis Syz-Modul:", nops(menge_S)); for zaehl_i from 1 to nops(menge_S) do poly_h:=mcollect(add((menge_S[zaehl_i][zaehl_j])*menge_Poly[zaehl_j],zaehl_j=1..nops(menge_Poly))); # print("Polynom", poly_h); poly_h:=hrmd(menge_Poly, poly_h); # print("Rest", poly_h); if (poly_h<>0) then hb:=false; break; fi; od; else hb:=false; fi; fi; fi; return(hb); end proc: HBasis:=proc(menge_Poly::list)::list; description "calculates H-Basis from a polynom set"; local menge_H::list, menge_G::list, poly_h::polynom, menge_S::list, zaehl_i::integer, zaehl_j::integer, full; global menge_Vspace, menge_Vscal, menge_Undets, menge_Hpoly; if (is_basis) then return(menge_Poly); fi; menge_H:=reduceSet(shortSet(menge_Poly)); menge_G:=[]; # print("H-Basis", menge_H); full:=false; if not (isHBasis(menge_H, full)) then while (menge_H<>menge_G) do menge_G:=menge_H; menge_Vspace:='menge_Vspace'; menge_Vscal:='menge_Vscal'; menge_Hpoly:='menge_Hpoly'; menge_Undets:='menge_Undets'; # print ("Berechne Basis der Syz-Modul"); menge_S:=aSyzMod(hltSet(menge_H)); # print ("Anzahl Basis Syz-Modul:", nops(menge_S)); for zaehl_i from 1 to nops(menge_S) do poly_h:=mcollect(add((menge_S[zaehl_i][zaehl_j])*menge_G[zaehl_j],zaehl_j=1..nops(menge_G))); # print(zaehl_i, "Polynom", poly_h); poly_h:=hrmd(menge_G, poly_h); # print(zaehl_i, "Rest", poly_h); if (poly_h<>0) then menge_H:=[op(menge_H), poly_h]; fi; # print("H-Basis reduzieren"); menge_H:=reduceSet(shortSet(menge_H)); # print("H-Basis reduzieren beendet"); od; # print("H-Basis ist", menge_H); if (isHBasis(menge_H, full)) then menge_G:=menge_H; fi; od; else menge_G:=menge_H; fi; menge_Vspace:='menge_Vspace'; menge_Vscal:='menge_Vscal'; menge_Hpoly:='menge_Hpoly'; menge_Undets:='menge_Undets'; return(shortSet(menge_G)); end proc: # ----------------------------------------------------------------------------------------- isZeroDim:=proc(menge_Poly::list)::bool; description "checks if ideal is zero-dimensional"; global basis_type; local zerodim::bool; zerodim:=true; if (dimIdeal(menge_Poly)>0) then zerodim:=false; fi; return(zerodim); end proc: NForm:=proc(menge_Poly::list)::list; description "returns basis for space of normal forms"; global basis_type; local menge_H::list, menge_B::list, menge_W::list, zaehl_i::integer, zaehl_j::integer, menge_P::list; menge_B:=[]; zaehl_i:=0; menge_W:=[1]; if not (isZeroDim(menge_Poly)) then error("Ideal is not zero-dimensional. Could not calculate normal form space."); fi; if (basis_type="G") then menge_H:=gltSet(GBasis(menge_Poly)); # menge_H:=gltSet(menge_Poly); while (menge_W<>[]) do menge_W:=[]; menge_P:=Pspace(zaehl_i); for zaehl_j from 1 to nops(menge_P) do # if (Groebner[normalf](menge_P[zaehl_j],menge_H,tdeg(op(List_Undets)))<>0) then if (grmd(menge_H, menge_P[zaehl_j])<>0) then menge_W:=[op(menge_W),menge_P[zaehl_j]]; fi; od; # print("NForm-Basis", menge_W); menge_B:=[op(menge_B), op(menge_W)]; zaehl_i:=zaehl_i+1; od; else menge_H:=HBasis(menge_Poly); # menge_H:=menge_Poly; while (menge_W<>[]) do menge_W:=Wspace(zaehl_i, menge_H); # print("NForm-Basis", menge_W); menge_B:=[op(menge_B), op(menge_W)]; zaehl_i:=zaehl_i+1; od; fi; return(shortSet(menge_B)); end proc: MultTable:=proc(menge_Poly::list, poly_h::polynom)::Matrix; description "calculates the multiplication table resp. basis menge_P for NF(menge_poly) and polynom poly_h"; global basis_type; local zaehl_i::integer, zaehl_j::integer, menge_U::list, poly_r::polynom, menge_H::list, menge_C::list, menge_P::list, poly_p::polynom, menge_Temp::list, menge_L::set, list_coeff; if (basis_type="G") then menge_H:=GBasis(menge_Poly); else menge_H:=HBasis(menge_Poly); fi; # print(menge_H); menge_P:=NForm(menge_H); # print(menge_P); menge_Temp:=[]; poly_p:=add(list_coeff[zaehl_j]*menge_P[zaehl_j], zaehl_j=1..nops(menge_P)); menge_U:=[seq(list_coeff[zaehl_j], zaehl_j=1..nops(menge_P))]; # print("Polynom p", poly_p); # print("Unbekannte", menge_U); for zaehl_i from 1 to nops(menge_P) do if (basis_type="G") then poly_r:=grmd(menge_H, menge_P[zaehl_i]*poly_h); else poly_r:=hrmd(menge_H, menge_P[zaehl_i]*poly_h); fi; # print(poly_r); # print("Polynom", poly_p-poly_r); menge_C:=[coeffs(mcollect(poly_p-poly_r), List_Undets)]; menge_L:={seq(menge_C[zaehl_i]=0, zaehl_i=1..nops(menge_C))}; # print(menge_L); # print(GenerateMatrix(menge_L, menge_U)); menge_Temp:=[op(menge_Temp),LinearSolve(GenerateMatrix(menge_L,menge_U, 'augmented'=true))]; # print(menge_Temp); od; return(Transpose(convert(menge_Temp, Matrix))); end proc: MultTableSet:=proc(menge_Poly::list, menge_Undet::list)::list; description "calculates the multiplication table resp. basis menge_P for NF(menge_poly) and polynom poly_h"; global basis_type; local matrix_M::Matrix, zaehl_i::integer, zaehl_j::integer, menge_U::list, poly_r::polynom, menge_H::list, menge_P::list, poly_p::polynom, menge_Temp::list, menge_L::set, list_coeff, menge_C::list, menge_M::list, zaehl_k::integer; if (basis_type="G") then menge_H:=GBasis(menge_Poly); else menge_H:=HBasis(menge_Poly); fi; menge_M:=[]; # print("Basis", menge_H); menge_P:=NForm(menge_H); # print("Normalform", menge_P); for zaehl_k from 1 to nops(menge_Undet) do # print(zaehl_k, "Unbekannte", menge_Undet[zaehl_k]); matrix_M:=Matrix(nops(menge_P), nops(menge_P), symbol=list_coeff); menge_Temp:=[]; for zaehl_i from 1 to nops(menge_P) do if (basis_type="G") then poly_r:=grmd(menge_H, menge_P[zaehl_i]*menge_Undet[zaehl_k]); else poly_r:=hrmd(menge_H, menge_P[zaehl_i]*menge_Undet[zaehl_k]); fi; # # print("Polynom r", poly_r); poly_p:=0; menge_U:=[]; for zaehl_j from 1 to nops(menge_P) do menge_U:=[op(menge_U), list_coeff[zaehl_j]]; poly_p:=poly_p+list_coeff[zaehl_j]*menge_P[zaehl_j]; od; # # print(menge_U); # # print("Polynom p", poly_p); menge_C:=[coeffs(mcollect(poly_r-poly_p), List_Undets)]; menge_L:={seq(menge_C[zaehl_i]=0, zaehl_i=1..nops(menge_C))}; # # print(menge_L); menge_Temp:=[op(menge_Temp),LinearSolve(GenerateMatrix(menge_L,menge_U, 'augmented'=true))]; od; matrix_M:=Transpose(convert(menge_Temp, Matrix)); menge_M:=[op(menge_M), matrix_M]; od; return(menge_M); end proc: traceMatrix:=proc(menge_Poly::list, poly_h::polynom)::Matrix; description "calculates the trace matrix resp. basis menge_P for NF(menge_poly) and polynom poly_h"; global basis_type; local zaehl_i::integer, zaehl_j::integer, menge_H::list, menge_P::list, matrix_T::Matrix, int_p::integer, matrix_M::Matrix; if (basis_type="G") then menge_H:=GBasis(menge_Poly); else menge_H:=HBasis(menge_Poly); fi; menge_P:=NForm(menge_H); int_p:=nops(menge_P); matrix_T:=Matrix(int_p, int_p); matrix_M:=Matrix(int_p, int_p); for zaehl_i from 1 to int_p do for zaehl_j from zaehl_i to int_p do matrix_M:=MultTable(menge_H, menge_P[zaehl_i]*menge_P[zaehl_j]*poly_h); matrix_T[zaehl_i,zaehl_j]:=Trace(matrix_M); matrix_T[zaehl_j,zaehl_i]:=matrix_T[zaehl_i,zaehl_j]; od; od; return(matrix_T); end proc: listEl:=proc(int_i::integer, int_j::integer, int_k::integer, int_n::integer)::integer; local menge_G::set, int_res::integer; menge_G:=sort([int_i,int_j,int_k]); int_res:=(menge_G[1]-1)*int_n^2+(menge_G[2]-1)*int_n+(menge_G[3]-1)+1; return(int_res); end proc: TraceMatrix:=proc(menge_Poly::list)::Matrix; description "returns the trace matrix of a polynom set"; local menge_H::list, menge_P::list, int_p::integer, zaehl_i::integer, zaehl_j::integer, matrix_T::Matrix, menge_C::list, poly_h::polynom, menge_L::set, menge_S::set, poly_p::polynom, menge_U::set, zaehl_k::integer, zaehl_l::integer, menge_P2::list, list_coeff, vector_v; global menge_Vspace, menge_Vscal, menge_Undets, menge_Hpoly, basis_type; if (basis_type="G") then menge_H:=GBasis(menge_Poly); # menge_H:=menge_Poly; else menge_H:=HBasis(menge_Poly); fi; menge_P:=NForm(menge_H); int_p:=nops(menge_P); # print(menge_P); # TraceMatrix berechnen matrix_T:=Matrix(int_p, int_p); menge_Vspace:='menge_Vspace'; menge_Vscal:='menge_Vscal'; menge_Hpoly:='menge_Hpoly'; menge_Undets:='menge_Undets'; for zaehl_i from 1 to int_p do for zaehl_j from zaehl_i to int_p do poly_h:=menge_P[zaehl_i]*menge_P[zaehl_j]; for zaehl_k from zaehl_j to int_p do # print("Rest berechnen:", zaehl_i, zaehl_j, zaehl_k, "von", menge_P[zaehl_k]*poly_h); if (basis_type="G") then menge_P2[listEl(zaehl_i, zaehl_j, zaehl_k, int_p)]:=grmd(menge_H, menge_P[zaehl_k]*poly_h); # menge_P2[listEl(zaehl_i, zaehl_j, zaehl_k, int_p)]:=Groebner[normalf](menge_P[zaehl_k]*poly_h, menge_H, tdeg(op(List_Undets))); else menge_P2[listEl(zaehl_i, zaehl_j, zaehl_k, int_p)]:=hrmd(menge_H, menge_P[zaehl_k]*poly_h); fi; od; od; od; for zaehl_i from 1 to int_p do for zaehl_j from 1 to int_p do poly_h:=menge_P[zaehl_i]*menge_P[zaehl_j]; menge_U:=[]; menge_L:={}; # print(zaehl_i, zaehl_j); for zaehl_k from 1 to int_p do poly_p:=0; for zaehl_l from 1 to int_p do menge_U:=[op(menge_U), list_coeff[zaehl_k,zaehl_l]]; poly_p:=poly_p+list_coeff[zaehl_k,zaehl_l]*menge_P[zaehl_l]; od; menge_C:=[coeffs(mcollect(menge_P2[listEl(zaehl_i, zaehl_j, zaehl_k, int_p)]-poly_p), List_Undets)]; menge_L:=menge_L union {seq(menge_C[zaehl_i]=0, zaehl_i=1..nops(menge_C))}; od; vector_v:=LinearSolve(GenerateMatrix(menge_L,menge_U, 'augmented'=true)); # # print(vector_v); matrix_T[zaehl_i,zaehl_j]:=add(vector_v[1+(int_p+1)*zaehl_l], zaehl_l=0..int_p-1); # Trace of Matrix od; od; # print("Tracematrix", matrix_T); return(matrix_T); end proc: Radical:=proc(menge_Poly::list)::list; description "returns radical of polynom-set (einfache Nullstellen), faster"; global basis_type; local menge_H::list, menge_P::list, int_p::integer, poly_temp::polynom, zaehl_i::integer, matrix_T::Matrix, menge_G::list, menge_N::set, vector_temp::Vector; if (basis_type="G") then menge_H:=GBasis(menge_Poly); else menge_H:=HBasis(menge_Poly); fi; menge_P:=NForm(menge_H); int_p:=nops(menge_P); # print(menge_P); matrix_T:=TraceMatrix(menge_H); menge_N:=NullSpace(matrix_T); menge_G:=menge_Poly; vector_temp:=convert(menge_P, Vector); for zaehl_i from 1 to nops(menge_N) do poly_temp:=Transpose(menge_N[zaehl_i]).vector_temp; menge_G:=[op(menge_G),poly_temp]; od; if (basis_type="G") then return(GBasis(menge_G)); else return(HBasis(menge_G)); fi; end proc: realSol:=proc(menge_Poly::list)::integer; description "returns number of real solutions"; global basis_type; local menge_H::list, menge_P::list, int_p::integer, poly_temp::polynom, exp_x, zaehl_i::integer, matrix_T::Matrix, int_v1, int_v2, poly_p, rsol; if (basis_type="G") then menge_H:=GBasis(menge_Poly); # menge_H:=menge_Poly; else menge_H:=HBasis(menge_Poly); fi; menge_P:=NForm(menge_H); int_p:=nops(menge_P); # print(menge_P); matrix_T:=TraceMatrix(menge_H); poly_p:=CharacteristicPolynomial(matrix_T, exp_x); # print("charpoly", poly_p); int_v1:=0; int_v2:=0; rsol:=0; for zaehl_i from 1 to nops(poly_p) do if (int_v1=0) then int_v1:=signum(coeff(op(zaehl_i,poly_p), exp_x, degree(op(zaehl_i,poly_p)))); else int_v2:=signum(coeff(op(zaehl_i,poly_p), exp_x, degree(op(zaehl_i,poly_p)))); if (int_v1<>int_v2) then rsol:=rsol+1; fi; int_v1:=int_v2; int_v2:=0; fi; od; rsol:=2*rsol-int_p+ldegree(poly_p); return(rsol); end proc: # --------------------------------------------------------------------------------- CartUnion:=proc(menge_M::list)::list; description "builds the union of elements of some listlist"; local menge_N::list, zaehl_i::integer, zaehl_j::integer, menge_Temp::list; # # print("ein"); # # print(menge_M, nops(menge_M)); menge_N:=[]; if (nops(menge_M)>=2) then if (nops(menge_M)=2) then # # print(nops(menge_M),nops(menge_M[1]),nops(menge_M[2])); for zaehl_i from 1 to nops(menge_M[1]) do for zaehl_j from 1 to nops(menge_M[2]) do # # print(zaehl_i, zaehl_j); menge_N:=[op(menge_N), [ [op(menge_M[1][zaehl_i][1]), op(menge_M[2][zaehl_j][1])],[op(menge_M[1][zaehl_i][2]), op(menge_M[2][zaehl_j][2])] ] ]; od; od; else # # print(nops(menge_M)); menge_Temp:=CartUnion(menge_M[2..-1]); for zaehl_i from 1 to nops(menge_M[1]) do for zaehl_j from 1 to nops(menge_Temp) do menge_N:=[op(menge_N), [ [op(menge_M[1][zaehl_i][1]), op(menge_Temp[zaehl_j][1])],[op(menge_M[1][zaehl_i][2]), op(menge_Temp[zaehl_j][2])] ] ]; od; od; fi; elif (nops(menge_M)=1) then # # print(nops(menge_M)); # print(menge_M); for zaehl_j from 1 to nops(menge_M[1]) do menge_N:=[op(menge_N), [ [op(menge_M[1][zaehl_j][1])],[op(menge_M[1][zaehl_j][2])] ] ]; od; fi; # # print(nops(menge_N[1][1])); # # print("aus"); return(menge_N); end proc: Eigenspaces:=proc(menge_Poly::list)::list; global List_Undets; local menge_M, menge_EW, matrix_M, zaehl_i, zaehl_j, menge_NS, menge_S, menge_Temp; menge_S:=[]; menge_M:=MultTableSet(menge_Poly,List_Undets); # print("MultTables", menge_M); for zaehl_i from 1 to nops(menge_M) do menge_Temp:=[]; # print("Matrix", menge_M[zaehl_i]); menge_EW:=convert(Eigenvalues(menge_M[zaehl_i], output=list),set); # print("Eigenwerte", menge_EW); for zaehl_j from 1 to nops(menge_EW) do # print("Eigenwert", menge_EW[zaehl_j]); matrix_M:=menge_M[zaehl_i]-menge_EW[zaehl_j]*IdentityMatrix(op(1,menge_M[zaehl_i])[1]); menge_NS:=collectSet(NullSpace(matrix_M)); # print("Eigenraum", menge_NS); menge_Temp:=[op(menge_Temp),[[menge_EW[zaehl_j]],[menge_NS]]]; od; menge_S:=[op(menge_S), menge_Temp]; od; return(menge_S); end proc: solveH:=proc(menge_Poly::list)::list; description "solves the problem menge_Poly=0"; global List_Undets; local menge_E::list, zaehl_i::integer, zaehl_j::integer, zaehl_k::integer, menge_Vl::list, menge_Wl::list, menge_CU::list, menge_Temp1::list, menge_S::list, menge_NF::list, menge_W::list, menge_Temp2::list, match::bool, menge_B::list; menge_S:=[]; menge_E:=[]; menge_W:=[]; menge_Wl:=[]; menge_Vl:=[]; menge_E:=Eigenspaces(menge_Poly); # print("Eigenraeume", menge_E); for zaehl_i from 1 to nops(menge_E) do if (nops(menge_E[zaehl_i])>1) then # <=> Vielfachheit EW = dim(ER) = dim(NF) menge_W:=[op(menge_W),menge_E[zaehl_i]]; menge_Wl:=[op(menge_Wl),List_Undets[zaehl_i]]; else menge_Vl:=[op(menge_Vl),List_Undets[zaehl_i]=menge_E[zaehl_i][1][1][1]]; fi; od; # print("Unabhaengige Variablen", menge_Vl); # print("Abhaengige Variablen", menge_Wl, menge_W); menge_CU:=CartUnion(menge_W); # print("Vereinigung", menge_CU); # print("Vereinigung abgeschlossen"); for zaehl_i from 1 to nops(menge_CU) do menge_Temp1:=menge_CU[zaehl_i][2]; match:=true; for zaehl_j from 1 to nops(menge_Temp1) do for zaehl_k from zaehl_j+1 to nops(menge_Temp1) do # print("Zaehler:", zaehl_i,zaehl_j,zaehl_k,"von",nops(menge_CU),nops(menge_Temp1),nops(menge_Temp1)); menge_Temp2:=[op(menge_Temp1[zaehl_j]), op(menge_Temp1[zaehl_k])]; # print("Vereinigung einzeln", menge_Temp2, nops(menge_Temp2)); menge_B:=Basis(menge_Temp2); # # print(nops(menge_B)); if (nops(menge_Temp2)>1 and nops(menge_B)>=nops(menge_Temp2)) then match:=false; fi; od; od; if (match) then # print("true, bei #Basis(Vektoren) < #Vektoren"); menge_Temp2:=menge_CU[zaehl_i][1]; for zaehl_j from 1 to nops(menge_Temp2) do menge_Temp2[zaehl_j]:=menge_Wl[zaehl_j]=menge_Temp2[zaehl_j]; od; menge_Temp2:=[op(menge_Temp2), op(menge_Vl)]; menge_S:=[op(menge_S), menge_Temp2]; fi; od; if (nops(menge_S)=0) then menge_S:=menge_Vl; fi; return(menge_S); end proc: checkSol:=proc(menge_Poly::list)::bool; description "checks the solution if it really solves the problem"; local menge_S::list, menge_G::list, int_sum::ineteger, zaehl_i::integer, zaehl_j::integer, check::bool; check:=true; menge_S:=solveH(menge_Poly); for zaehl_i from 1 to nops(menge_S) do menge_G:=simplify(subs(menge_S[zaehl_i], menge_Poly)); int_sum:=add(abs(menge_G[zaehl_j]), zaehl_j=1..nops(menge_G)); if (int_sum<>0) then check:=false; fi; od; return(check); end proc: # ------------------------------------------------------------------------------------------ drawKin :=proc(kin::string, menge_S::set, menge_T::list, value) local menge_C, zaehl_i, zaehl_j, int_l1, int_l2, int_l3, int_x, int_y, int_z, step1, step2, step3, nsteps, int_temp, vector_B1, vector_B2, vector_B3, vector_A1, vector_A2, vector_A3, points_P, vector_X, menge_T2, vector_O, vector_O3, menge_U, menge_C2, menge_Sol, matrix_M; global infolevel, List_Consts, List_Undets; nsteps:=10; points_P:=[]; vector_O:=[0,0]; vector_O3:=[0,0,0]; if (nops(menge_T)=0) then return; fi; if (kin="2D1s") then menge_C:={a}; if (nops(indets(subs(menge_S, menge_C)))>0) then error "All constants %1 must be substituted. %2 is missing.",menge_C, indets(subs(menge_S, menge_C)); fi; if (nops(menge_T) mod 2 <> 0) then error "List of parameters must contain 2*n entries."; fi; if (nops(menge_T)=2) then menge_T2:=[op(menge_T), op(menge_T)]; else menge_T2:=menge_T; fi; vector_A1:=[-a,0]; vector_A2:=[a,0]; if (value=inverse) then for zaehl_i from 2 to nops(menge_T2) by 2 do if (menge_T2[zaehl_i]>=0) then error "List of parameters in y could only contain negative entries."; fi; od; # int_l1:=sqrt(subs(menge_S, x^2-2*a*x+a^2+y^2)); # int_l2:=sqrt(subs(menge_S, x^2+2*a*x+a^2+y^2)); for zaehl_i from 3 to nops(menge_T2)-1 by 2 do step1:=(menge_T2[zaehl_i]-menge_T2[zaehl_i-2])/nsteps; step2:=(menge_T2[zaehl_i+1]-menge_T2[zaehl_i-1])/nsteps; for zaehl_j from 0 to nsteps do vector_X:=[menge_T2[zaehl_i-2]+zaehl_j*step1, menge_T2[zaehl_i-1]+zaehl_j*step2]; points_P:=[op(points_P), plot(evalf(subs(menge_S, [vector_A1, vector_O, vector_A2, vector_X, vector_A1])))]; od; od; else for zaehl_i from 1 to nops(menge_T2) do if (menge_T2[zaehl_i]<0) then error "List of parameters could only contain positive entries."; fi; od; int_x:=-1/4*(-l2^2+l1^2)/a; int_y:=Re(-1/4*(-l2^4+8*l1^2*a^2-16*a^4+8*a^2*l2^2+2*l2^2*l1^2-l1^4)^(1/2)/a); for zaehl_i from 3 to nops(menge_T2)-1 by 2 do step1:=(menge_T2[zaehl_i]-menge_T2[zaehl_i-2])/nsteps; step2:=(menge_T2[zaehl_i+1]-menge_T2[zaehl_i-1])/nsteps; for zaehl_j from 0 to nsteps do vector_X:=subs({l1=menge_T2[zaehl_i-2]+zaehl_j*step1, l2=menge_T2[zaehl_i-1]+zaehl_j*step2}, [int_x, int_y]); points_P:=[op(points_P), plot(evalf(subs(menge_S, [vector_A1, vector_O, vector_A2, vector_X, vector_A1])))]; od; od; fi; elif (kin="2D1a") then menge_C:={x1,y1,x2,y2}; if (nops(indets(subs(menge_S, menge_C)))>0) then error "All constants %1 must be substituted. %2 is missing.",menge_C, indets(subs(menge_S, menge_C)); fi; if (nops(menge_T) mod 2 <> 0) then error "List of parameters must contain 2*n entries."; fi; if (nops(menge_T)=2) then menge_T2:=[op(menge_T), op(menge_T)]; else menge_T2:=menge_T; fi; vector_A1:=[x1,y1]; vector_A2:=[x2,y2]; if (value=inverse) then for zaehl_i from 2 to nops(menge_T2) by 2 do if (menge_T2[zaehl_i]>=0) then error "List of parameters in y could only contain negative entries."; fi; od; # int_l1:=sqrt(subs(menge_S, (x1-x)^2+(y+y1)^2)); # int_l2:=sqrt(subs(menge_S, (x2-x)^2+(y+y2)^2)); for zaehl_i from 3 to nops(menge_T2)-1 by 2 do step1:=(menge_T2[zaehl_i]-menge_T2[zaehl_i-2])/nsteps; step2:=(menge_T2[zaehl_i+1]-menge_T2[zaehl_i-1])/nsteps; for zaehl_j from 0 to nsteps do vector_X:=[menge_T2[zaehl_i-2]+zaehl_j*step1, menge_T2[zaehl_i-1]+zaehl_j*step2]; points_P:=[op(points_P), plot(evalf(subs(menge_S, [vector_A1, vector_O, vector_A2, vector_X, vector_A1])))]; od; od; else for zaehl_i from 1 to nops(menge_T2) do if (menge_T2[zaehl_i]<0) then error "List of parameters could only contain positive entries."; fi; od; int_x:=Re((4*x2^3+4*x1*y1^2+4*x1*y2^2-4*x2*x1^2-4*x1*x2^2+4*x2*y2^2+4*x2*y1^2+4*x1*l2^2-4*x1*l1^2-4*x2*l2^2+4*x2*l1^2-8*x2*y1*y2+4*x1^3-8*x1*y1*y2-4*(24*y1^2*y2^2*x2*x1-16*y1*y2^3*x2*x1-16*y2*y1^3*x2*x1-4*y2^2*l1^2*x2*x1-4*y2^2*l2^2*x2*x1-4*y1^2*l1^2*x2*x1-4*y1^2*l2^2*x2*x1-4*y1*l1^2*l2^2*y2-4*y2*l2^2*x2^2*y1-8*y2*x1*x2^3*y1-4*y2*l1^2*x1^2*y1+12*y1*x2^2*x1^2*y2-4*y2*l1^2*x2^2*y1-8*y2*x1^3*x2*y1-4*y1*l2^2*x1^2*y2-y1^2*l2^4-y1^2*l1^4-y1^2*x2^4+2*y2^4*l2^2-2*x1^2*y2^4+6*y2*y1^5-15*y2^2*y1^4+2*y1^4*l2^2-15*y1^2*y2^4+6*y1*y2^5-2*x1^2*y1^4+20*y2^3*y1^3+2*y1^4*l1^2-2*y2^4*x2^2-2*x2^2*y1^4-y2^2*l2^4-y2^2*l1^4-y2^2*x2^4-y2^2*x1^4+8*x1^2*y2^3*y1+2*y1*x1^4*y2+2*y1^2*l1^2*l2^2+2*y2^2*l2^2*x2^2+4*y2^4*x2*x1+12*y1^2*y2^2*l2^2-8*y1*y2^3*l2^2+4*y2^2*x1*x2^3+4*y1^4*x2*x1+12*y1^2*l1^2*y2^2-8*y1*l1^2*y2^3+8*y2*x1^2*y1^3-12*y2^2*x1^2*y1^2-8*y2*y1^3*l2^2-8*y2*y1^3*l1^2+8*y2*x2^2*y1^3-12*y2^2*x2^2*y1^2+2*y1^2*l2^2*x1^2+4*y2^2*x1^3*x2+2*y2^2*l1^2*x1^2+2*y2^2*l1^2*l2^2-6*y1^2*x2^2*x1^2+2*y2^2*l1^2*x2^2+2*y1^2*l1^2*x2^2+2*y2^2*l2^2*x1^2+2*y1^2*l1^2*x1^2+2*y1^2*l2^2*x2^2+4*y1^2*x1^3*x2+4*y1^2*x1*x2^3-6*y2^2*x2^2*x1^2+8*y2^3*x2^2*y1+2*y1*l2^4*y2+2*y1*l1^4*y2+2*y1*x2^4*y2+2*l1^2*y2^4-y1^2*x1^4+8*y2*l1^2*x2*x1*y1+8*y2*l2^2*x2*x1*y1-y1^6-y2^6)^(1/2))/(2*(-8*y1*y2-8*x2*x1+4*y2^2+4*x1^2+4*y1^2+4*x2^2))); int_y:=Re((8*y1*x1*x2-4*y2*l1^2+4*y1*y2^2-4*y1*x2^2-4*y2*x2^2+4*y2*l2^2+4*y1*l1^2-4*y2*x1^2+4*y2*y1^2-4*y1*l2^2-4*y1^3-4*y1*x1^2-4*y2^3+8*x1*x2*y2-4*(20*x2^3*x1^3-15*x2^4*x1^2+2*x1^4*l2^2+6*x2^5*x1+6*x1^5*x2-15*x1^4*x2^2+2*l1^2*x2^4+2*x1^4*l1^2-x1^2*l2^4-x1^2*l1^4-l1^4*x2^2+12*y1^2*y2^2*x2*x1-8*y1*y2^3*x2*x1-8*y2*y1^3*x2*x1-4*y2^2*l1^2*x2*x1-4*y2^2*l2^2*x2*x1-4*y1^2*l1^2*x2*x1-4*y1^2*l2^2*x2*x1-4*y2*l2^2*x2^2*y1-16*y2*x1*x2^3*y1-4*y2*l1^2*x1^2*y1+24*y1*x2^2*x1^2*y2-4*y2*l1^2*x2^2*y1-16*y2*x1^3*x2*y1-4*y1*l2^2*x1^2*y2-2*y1^2*x2^4-x1^2*y2^4-x1^2*y1^4+12*x1^2*l2^2*x2^2+12*x1^2*l1^2*x2^2+2*x1^2*l1^2*l2^2+2*l1^2*l2^2*x2^2-8*x1^3*l2^2*x2-8*l1^2*x2^3*x1-8*x1^3*l1^2*x2+2*x1*l2^4*x2+2*x1*l1^4*x2-8*l2^2*x2^3*x1-y2^4*x2^2-x2^2*y1^4-2*y2^2*x2^4-2*y2^2*x1^4+4*x1^2*y2^3*y1+4*y1*x1^4*y2+2*y2^2*l2^2*x2^2+2*y2^4*x2*x1+8*y2^2*x1*x2^3+2*y1^4*x2*x1+4*y2*x1^2*y1^3-6*y2^2*x1^2*y1^2+4*y2*x2^2*y1^3-6*y2^2*x2^2*y1^2+2*y1^2*l2^2*x1^2+8*y2^2*x1^3*x2+2*y2^2*l1^2*x1^2-12*y1^2*x2^2*x1^2+2*y2^2*l1^2*x2^2+2*y1^2*l1^2*x2^2+2*y2^2*l2^2*x1^2+2*y1^2*l1^2*x1^2+2*y1^2*l2^2*x2^2+8*y1^2*x1^3*x2+8*y1^2*x1*x2^3-12*y2^2*x2^2*x1^2+4*y2^3*x2^2*y1+4*y1*x2^4*y2+2*l2^2*x2^4-l2^4*x2^2-x1^6-x2^6-4*l1^2*l2^2*x2*x1-2*y1^2*x1^4+8*y2*l1^2*x2*x1*y1+8*y2*l2^2*x2*x1*y1)^(1/2))/(2*(-8*y1*y2-8*x2*x1+4*y2^2+4*x1^2+4*y1^2+4*x2^2))); for zaehl_i from 3 to nops(menge_T2)-1 by 2 do step1:=(menge_T2[zaehl_i]-menge_T2[zaehl_i-2])/nsteps; step2:=(menge_T2[zaehl_i+1]-menge_T2[zaehl_i-1])/nsteps; for zaehl_j from 0 to nsteps do vector_X:=subs({l1=menge_T2[zaehl_i-2]+zaehl_j*step1, l2=menge_T2[zaehl_i-1]+zaehl_j*step2}, [int_x, int_y]); points_P:=[op(points_P), plot(evalf(subs(menge_S, [vector_A1, vector_O, vector_A2, vector_X, vector_A1])))]; od; od; fi; elif (kin="2D2") then menge_C:={a,b}; if (nops(indets(subs(menge_S, menge_C)))>0) then error "All constants %1 must be substituted. %2 is missing.",menge_C, indets(subs(menge_S, menge_C)); fi; if (nops(menge_T) mod 2 <> 0) then error "List of parameters must contain 2*n entries."; fi; if (nops(menge_T)=2) then menge_T2:=[op(menge_T), op(menge_T)]; else menge_T2:=menge_T; fi; vector_A1:=[-a,0]; vector_A2:=[a,0]; if (value=inverse) then for zaehl_i from 2 to nops(menge_T2) by 2 do if (menge_T2[zaehl_i]>=0) then error "List of parameters in y could only contain negative entries."; fi; od; # int_l1:=((-2*x*a*y^2-2*x^3*a+2*a*y*b*(y^2+x^2)^(1/2)+2*x^2*y^2+x^4+y^2*b^2+x^2*b^2+x^2*a^2+y^4+y^2*a^2)/(y^2+x^2))^(1/2); # int_l2:=(( 2*x*a*y^2+2*x^3*a+2*a*y*b*(y^2+x^2)^(1/2)+2*x^2*y^2+x^4+y^2*b^2+x^2*b^2+x^2*a^2+y^4+y^2*a^2)/(y^2+x^2))^(1/2); for zaehl_i from 3 to nops(menge_T2)-1 by 2 do step1:=(menge_T2[zaehl_i]-menge_T2[zaehl_i-2])/nsteps; step2:=(menge_T2[zaehl_i+1]-menge_T2[zaehl_i-1])/nsteps; for zaehl_j from 0 to nsteps do vector_X:=[menge_T2[zaehl_i-2]+zaehl_j*step1, menge_T2[zaehl_i-1]+zaehl_j*step2]; vector_B2:=b*convert(Normalize(convert([1,-vector_X[1]/vector_X[2]], Vector), Euclidean), list); vector_B1:=-vector_B2; points_P:=[op(points_P), plot(evalf(subs(menge_S, [vector_A1, vector_O, vector_A2, vector_X+vector_B2, vector_X+vector_B1, vector_A1, vector_O, vector_X])))]; od; od; else for zaehl_i from 1 to nops(menge_T2) do if (menge_T2[zaehl_i]<0) then error "List of parameters could only contain positive entries."; fi; od; int_x:=-1/4/a*(-l2^2+l1^2); int_y:=Re(-1/12*(24*(66*b^2*l2^2*l1^2+8*a^6+120*l1^2*b^2*a^2+120*l2^2*b^2*a^2-12*l2^2*a^4-264*a^4*b^2-12*l1^2*a^4-264*a^2*b^4+8*b^6-12*b^4*l1^2-12*b^4*l2^2-3*l1^4*l2^2+12*l1^2*l2^2*a^2-21*b^2*l2^4-21*b^2*l1^4-3*l1^2*l2^4+6*l2^4*a^2-l1^6-l2^6+6*a^2*l1^4+3*(-84*l1^8*b^2*a^2-192*l1^2*l2^6*b^2*a^2-84*l2^8*b^2*a^2+456*a^4*l1^6*b^2+456*a^4*l2^6*b^2+816*a^4*b^4*l2^4+816*a^4*b^4*l1^4-1200*a^6*l2^4*b^2-1200*a^6*l1^4*b^2+1080*a^4*l1^2*b^2*l2^4+1080*a^4*l1^4*b^2*l2^2+1872*l1^4*b^4*l2^2*a^2+1872*l1^2*b^4*l2^4*a^2-5472*l1^2*b^6*l2^2*a^2+432*l1^4*b^6*a^2-336*l1^6*b^4*a^2+1536*a^8*l1^2*b^2+1536*a^8*l2^2*b^2+96*l1^2*b^8*l2^2+1536*l1^2*b^8*a^2-72*l1^2*b^6*l2^4-72*l1^4*b^6*l2^2-324*l1^2*b^4*l2^6-324*l1^6*b^4*l2^2+558*l1^4*b^4*l2^4-1536*l1^2*b^6*a^4-12*l1^6*b^2*l2^4+6*l1^8*b^2*l2^2+6*l1^2*b^2*l2^8-12*l1^4*b^2*l2^6-1536*a^6*b^4*l1^2-1536*a^6*b^4*l2^2-336*l2^6*b^4*a^2+1536*l2^2*b^8*a^2+432*l2^4*b^6*a^2-1536*l2^2*b^6*a^4-2208*a^6*l1^2*l2^2*b^2-216*l1^4*l2^4*b^2*a^2-192*l1^6*l2^2*b^2*a^2-4704*a^4*b^4*l2^2*l1^2+3072*a^4*b^8-768*a^2*b^10-48*l1^4*b^8+45*l1^8*b^4+72*l1^6*b^6+6*l1^10*b^2-48*l2^4*b^8+72*l2^6*b^6+45*l2^8*b^4-768*a^10*b^2+3072*a^8*b^4-4608*a^6*b^6+6*l2^10*b^2)^(1/2))^(1/3)-864*(1/9*l1^2*b^2-1/18*l1^2*l2^2+1/9*l2^2*b^2-1/36*l1^4+1/9*l2^2*a^2-14/9*b^2*a^2+1/9*a^2*l1^2-1/36*l2^4-1/9*b^4-1/9*a^4)/(66*b^2*l2^2*l1^2+8*a^6+120*l1^2*b^2*a^2+120*l2^2*b^2*a^2-12*l2^2*a^4-264*a^4*b^2-12*l1^2*a^4-264*a^2*b^4+8*b^6-12*b^4*l1^2-12*b^4*l2^2-3*l1^4*l2^2+12*l1^2*l2^2*a^2-21*b^2*l2^4-21*b^2*l1^4-3*l1^2*l2^4+6*l2^4*a^2-l1^6-l2^6+6*a^2*l1^4+3*(-84*l1^8*b^2*a^2-192*l1^2*l2^6*b^2*a^2-84*l2^8*b^2*a^2+456*a^4*l1^6*b^2+456*a^4*l2^6*b^2+816*a^4*b^4*l2^4+816*a^4*b^4*l1^4-1200*a^6*l2^4*b^2-1200*a^6*l1^4*b^2+1080*a^4*l1^2*b^2*l2^4+1080*a^4*l1^4*b^2*l2^2+1872*l1^4*b^4*l2^2*a^2+1872*l1^2*b^4*l2^4*a^2-5472*l1^2*b^6*l2^2*a^2+432*l1^4*b^6*a^2-336*l1^6*b^4*a^2+1536*a^8*l1^2*b^2+1536*a^8*l2^2*b^2+96*l1^2*b^8*l2^2+1536*l1^2*b^8*a^2-72*l1^2*b^6*l2^4-72*l1^4*b^6*l2^2-324*l1^2*b^4*l2^6-324*l1^6*b^4*l2^2+558*l1^4*b^4*l2^4-1536*l1^2*b^6*a^4-12*l1^6*b^2*l2^4+6*l1^8*b^2*l2^2+6*l1^2*b^2*l2^8-12*l1^4*b^2*l2^6-1536*a^6*b^4*l1^2-1536*a^6*b^4*l2^2-336*l2^6*b^4*a^2+1536*l2^2*b^8*a^2+432*l2^4*b^6*a^2-1536*l2^2*b^6*a^4-2208*a^6*l1^2*l2^2*b^2-216*l1^4*l2^4*b^2*a^2-192*l1^6*l2^2*b^2*a^2-4704*a^4*b^4*l2^2*l1^2+3072*a^4*b^8-768*a^2*b^10-48*l1^4*b^8+45*l1^8*b^4+72*l1^6*b^6+6*l1^10*b^2-48*l2^4*b^8+72*l2^6*b^6+45*l2^8*b^4-768*a^10*b^2+3072*a^8*b^4-4608*a^6*b^6+6*l2^10*b^2)^(1/2))^(1/3)-3*(32*a^4+3*l2^4-6*l1^2*l2^2+3*l1^4+32*b^2*a^2-16*a^2*l1^2-16*l2^2*a^2)/a^2)^(1/2)); for zaehl_i from 3 to nops(menge_T2)-1 by 2 do step1:=(menge_T2[zaehl_i]-menge_T2[zaehl_i-2])/nsteps; step2:=(menge_T2[zaehl_i+1]-menge_T2[zaehl_i-1])/nsteps; for zaehl_j from 0 to nsteps do vector_X:=subs({l1=menge_T2[zaehl_i-2]+zaehl_j*step1, l2=menge_T2[zaehl_i-1]+zaehl_j*step2}, [int_x, int_y]); vector_B2:=b*convert(Normalize(convert([1,-vector_X[1]/vector_X[2]], Vector), Euclidean), list); vector_B1:=-vector_B2; points_P:=[op(points_P), plot(evalf(subs(menge_S, [vector_A1, vector_O, vector_A2, vector_X+vector_B2, vector_X+vector_B1, vector_A1, vector_O, vector_X])))]; od; od; fi; elif (kin="3D1s") then menge_C:={a}; if (nops(indets(subs(menge_S, menge_C)))>0) then error "All constants %1 must be substituted. %2 is missing.",menge_C, indets(subs(menge_S, menge_C)); fi; if (nops(menge_T) mod 3 <> 0) then error "List of parameters must contain 3*n entries."; fi; if (nops(menge_T)=3) then menge_T2:=[op(menge_T), op(menge_T)]; else menge_T2:=menge_T; fi; vector_A1:=[0,a,0]; vector_A2:=[-1*sqrt(3)*a/2, -a/2,0]; vector_A3:=[sqrt(3)*a/2, -a/2,0]; if (value=inverse) then for zaehl_i from 2 to nops(menge_T2) by 3 do if (menge_T2[zaehl_i]>=0) then error "List of parameters in y could only contain negative entries."; fi; int_temp:=menge_T2[zaehl_i]; menge_T2[zaehl_i]:=menge_T2[zaehl_i+1]; menge_T2[zaehl_i+1]:=int_temp; od; # int_l1:=sqrt(x^2+y^2+(a-z)^2); # int_l2:=sqrt((-1/2*3^(1/2)*a-x)^2+y^2+(-1/2*a-z)^2); # int_l3:=sqrt((1/2*3^(1/2)*a-x)^2+y^2+(-1/2*a-z)^2); for zaehl_i from 4 to nops(menge_T2)-2 by 3 do step1:=(menge_T2[zaehl_i]-menge_T2[zaehl_i-3])/nsteps; step2:=(menge_T2[zaehl_i+1]-menge_T2[zaehl_i-2])/nsteps; step3:=(menge_T2[zaehl_i+2]-menge_T2[zaehl_i-1])/nsteps; for zaehl_j from 0 to nsteps do vector_X:=[menge_T2[zaehl_i-3]+zaehl_j*step1, menge_T2[zaehl_i-2]+zaehl_j*step2, menge_T2[zaehl_i-1]+zaehl_j*step3]; points_P:=[op(points_P), plots[surfdata](evalf(subs(menge_S, [[vector_A1,vector_A2,vector_A3,vector_A1],[vector_X,vector_X,vector_X,vector_X]])), style=wireframe)]; od; od; else for zaehl_i from 1 to nops(menge_T2) do if (menge_T2[zaehl_i]<0) then error "List of parameters could only contain positive entries."; fi; od; int_x:=1/6*3^(1/2)*(l2^2-l3^2)/a; int_y:=Re(-1/3/a*(-l3^4-9*a^4-l2^4+l2^2*l3^2+3*l3^2*a^2+3*l1^2*a^2+l2^2*l1^2+3*l2^2*a^2+l3^2*l1^2-l1^4)^(1/2)); int_z:=-1/6*(-l2^2+2*l1^2-l3^2)/a; # print(int_x,int_y, int_z); for zaehl_i from 4 to nops(menge_T2)-2 by 3 do step1:=(menge_T2[zaehl_i]-menge_T2[zaehl_i-3])/nsteps; step2:=(menge_T2[zaehl_i+1]-menge_T2[zaehl_i-2])/nsteps; step3:=(menge_T2[zaehl_i+2]-menge_T2[zaehl_i-1])/nsteps; for zaehl_j from 0 to nsteps do vector_X:=subs({l1=menge_T2[zaehl_i-3]+zaehl_j*step1, l2=menge_T2[zaehl_i-2]+zaehl_j*step2, l3=menge_T2[zaehl_i-1]+zaehl_j*step3}, [int_x, int_z, int_y]); # print(vector_X); points_P:=[op(points_P), plots[surfdata](evalf(subs(menge_S, [[vector_A1,vector_A2,vector_A3,vector_A1],[vector_X,vector_X,vector_X,vector_X]])), style=wireframe)]; od; od; fi; elif (kin="3D1a") then menge_C:={x1,y1,z1,x2,y2,z2,x3,y3,z3}; if (nops(indets(subs(menge_S, menge_C)))>0) then error "All constants %1 must be substituted. %2 is missing.",menge_C, indets(subs(menge_S, menge_C)); fi; if (nops(menge_T) mod 3 <> 0) then error "List of parameters must contain 3*n entries."; fi; if (nops(menge_T)=3) then menge_T2:=[op(menge_T), op(menge_T)]; else menge_T2:=menge_T; fi; vector_A1:=[x1,z1,y1]; vector_A2:=[x2,z2,y2]; vector_A3:=[x3,z3,y3]; if (value=inverse) then for zaehl_i from 2 to nops(menge_T2) by 3 do if (menge_T2[zaehl_i]>=0) then error "List of parameters in y could only contain negative entries."; fi; int_temp:=menge_T2[zaehl_i]; menge_T2[zaehl_i]:=menge_T2[zaehl_i+1]; menge_T2[zaehl_i+1]:=int_temp; od; # int_l1:=sqrt((x1-x)^2+(y1-y)^2+(z1-z)^2); # int_l2:=sqrt((x2-x)^2+(y2-y)^2+(z2-z)^2); # int_l3:=sqrt((x3-x)^2+(y3-y)^2+(z3-z)^2); for zaehl_i from 4 to nops(menge_T2)-2 by 3 do step1:=(menge_T2[zaehl_i]-menge_T2[zaehl_i-3])/nsteps; step2:=(menge_T2[zaehl_i+1]-menge_T2[zaehl_i-2])/nsteps; step3:=(menge_T2[zaehl_i+2]-menge_T2[zaehl_i-1])/nsteps; for zaehl_j from 0 to nsteps do vector_X:=[menge_T2[zaehl_i-3]+zaehl_j*step1, menge_T2[zaehl_i-2]+zaehl_j*step2, menge_T2[zaehl_i-1]+zaehl_j*step3]; points_P:=[op(points_P), plots[surfdata](evalf(subs(menge_S, [[vector_A1,vector_A2,vector_A3,vector_A1],[vector_X,vector_X,vector_X,vector_X]])), style=wireframe)]; od; od; else for zaehl_i from 1 to nops(menge_T2) do if (menge_T2[zaehl_i]<0) then error "List of parameters could only contain positive entries."; fi; od; infolevel[information]:=2; if (type(List_Undets,list)) then menge_U:=List_Undets; fi; if (type(List_Consts,list)) then menge_C2:=List_Consts; fi; setUndets([x,z,y]); setConsts([l1,l2,l3]); menge_Sol:=solveH(subs(menge_S, [(x1-x)^2+(y1-y)^2+(z1-z)^2-(l1)^2, (x2-x)^2+(y2-y)^2+(z2-z)^2-(l2)^2, (x3-x)^2+(y3-y)^2+(z3-z)^2-(l3)^2]))[2]; if (type(menge_U,list)) then setUndets(menge_U); else List_Undets:='List_Undets'; fi; if (type(menge_C2,list)) then setConsts(menge_C2); else List_Consts:='List_Consts'; fi; infolevel[information]:=3; # print(menge_Sol); for zaehl_i from 1 to nops(menge_Sol) do if (lhs(menge_Sol[zaehl_i])='x') then int_x:=rhs(menge_Sol[zaehl_i]); elif (lhs(menge_Sol[zaehl_i])='y') then int_y:=Re(rhs(menge_Sol[zaehl_i])); elif (lhs(menge_Sol[zaehl_i])='z') then int_z:=rhs(menge_Sol[zaehl_i]); fi; od; # print(int_x,int_y, int_z); for zaehl_i from 4 to nops(menge_T2)-2 by 3 do step1:=(menge_T2[zaehl_i]-menge_T2[zaehl_i-3])/nsteps; step2:=(menge_T2[zaehl_i+1]-menge_T2[zaehl_i-2])/nsteps; step3:=(menge_T2[zaehl_i+2]-menge_T2[zaehl_i-1])/nsteps; for zaehl_j from 0 to nsteps do vector_X:=subs({l1=menge_T2[zaehl_i-3]+zaehl_j*step1, l2=menge_T2[zaehl_i-2]+zaehl_j*step2, l3=menge_T2[zaehl_i-1]+zaehl_j*step3}, [int_x, int_z, int_y]); # print(vector_X); points_P:=[op(points_P), plots[surfdata](evalf(subs(menge_S, [[vector_A1,vector_A2,vector_A3,vector_A1],[vector_X,vector_X,vector_X,vector_X]])), style=wireframe)]; od; od; fi; elif (kin="3D2") then menge_C:={a,b}; if (nops(indets(subs(menge_S, menge_C)))>0) then error "All constants %1 must be substituted. %2 is missing.",menge_C, indets(subs(menge_S, menge_C)); fi; if (nops(menge_T) mod 3 <> 0) then error "List of parameters must contain 3*n entries."; fi; if (nops(menge_T)=3) then menge_T2:=[op(menge_T), op(menge_T)]; else menge_T2:=menge_T; fi; vector_A1:=[0,a,0]; vector_A2:=[-1*sqrt(3)*a/2, -a/2,0]; vector_A3:=[sqrt(3)*a/2, -a/2,0]; if (value=inverse) then for zaehl_i from 2 to nops(menge_T2) by 3 do if (menge_T2[zaehl_i]>=0) then error "List of parameters in y could only contain negative entries."; fi; int_temp:=menge_T2[zaehl_i]; menge_T2[zaehl_i]:=menge_T2[zaehl_i+1]; menge_T2[zaehl_i+1]:=int_temp; od; for zaehl_i from 4 to nops(menge_T2)-2 by 3 do step1:=(menge_T2[zaehl_i]-menge_T2[zaehl_i-3])/nsteps; step2:=(menge_T2[zaehl_i+1]-menge_T2[zaehl_i-2])/nsteps; step3:=(menge_T2[zaehl_i+2]-menge_T2[zaehl_i-1])/nsteps; for zaehl_j from 0 to nsteps do vector_X:=[menge_T2[zaehl_i-3]+zaehl_j*step1, menge_T2[zaehl_i-2]+zaehl_j*step2, menge_T2[zaehl_i-1]+zaehl_j*step3]; vector_B1:=[0,-l,b/sqrt(3)]; vector_B2:=[-b/2, -l, -b/sqrt(3)/2]; vector_B3:=[b/2, -l, -b/sqrt(3)/2]; matrix_M:=<<-y/d, 0, x/d>|<-x/l,-z/l, -y/l>|<-z*x/l/d,d/l,-z*y/l/d>>; vector_B1:=convert(matrix_M.convert(vector_B1, Vector), list); vector_B2:=convert(matrix_M.convert(vector_B2, Vector), list); vector_B3:=convert(matrix_M.convert(vector_B3, Vector), list); # print(vector_B1, vector_B2, vector_B3); points_P:=[op(points_P), plots[surfdata](evalf(subs(menge_S union {x=vector_X[1], y=vector_X[3], z=vector_X[2], l=sqrt(vector_X[1]^2+vector_X[2]^2+vector_X[3]^2), d=sqrt(vector_X[2]^2+vector_X[3]^2)}, [[vector_O3, vector_A1, vector_A2, vector_O3, vector_A3, vector_A1, vector_A2, vector_A3],[vector_X, vector_B1,vector_B2,vector_X,vector_B3,vector_B1,vector_B2,vector_B3]])), style=wireframe)]; od; od; else for zaehl_i from 1 to nops(menge_T2) do if (menge_T2[zaehl_i]<0) then error "List of parameters could only contain positive entries."; fi; od; userinfo(3, 'information', "Solution for kinematic transformation of 3D-Kinematik 2 is not computed yet."); return; fi; else error "Drawing of kinematik %1 not supported.",kin; fi; # print(points_P); plots[display](points_P, insequence=true, scaling=constrained, color=black, axes=normal); end proc: setBasis(false); infolevel[information]:=3: