Package $contrib/coala -- COmputational ALgebraic Analysis Alias Coala := $contrib/coala; Define About() PrintLn " KeyWords : Algebraic Analysis, Clifford Operator, Cauchy-Fueter Operator, Noetherian Operators, Megaforms Author : A.Damiano E-mail : alberto@tlc185.com Version : CoCoA4.4 Date : July 2005 Thanks To: The author is thankful to Daniele Struppa, Irene Sabadini and Fabrizio Colombo for working with him on this subject and allowing him to play with CoCoA " End; -- About ------[ Manual ]-------- Define Man() PrintLn " Suggested alias for this package: Alias Coala := $contrib/coala; (please note that starting with CoCoA version 4.3 this alias is automatically set) SYNTAX please visit the website www.tlc185.com/coala for an updated version of both the manual and the package " End; -- Man ------[ Matrices and Resolutions ]-------- Define MyRes(M); -- I: a MODULE or IDEAL, not necessarily homogeneous -- O: a LIST of the ranks of A free resolution, starting from the right K:=0; S:=M; BET:=[Len(S.Gens)]; While K<10 And Not(IsZero(S)) Do S:=SyzOfGens(S); -- IT USES SYZOFGENS!!! If IsZero(S) Then B:=0; Else B:=Len(S.Gens); EndIf; PrintLn(B); Append(BET,B); K:=K+1; EndWhile; Return BET; EndDefine; Define MyResMap(M); -- I: a MODULE or IDEAL, not necessarily homogeneous -- O: a LIST of the maps of a free resolution, starting from the right K:=0; S:=M; MAPS:=[S.Gens]; While K<7 And Not(IsZero(S)) Do S:=SyzOfGens(S); S:=Minimalized(S); Append(MAPS,S.Gens); K:=K+1; EndWhile; Return MAPS; EndDefine; ------------------------------ CAUCHY FUETER Define FueterMat(X); -- I: a LIST of variables of length 4N -- O: the MATrix symbol of the left CF operator in N quaternionic variables N:=Len(X)/4; Q:=[]; For J:=1 To N Do K:=4*(J-1); Q:=Coala.MatConcat(Q,Mat[[X[1+K],-X[2+K],-X[3+K],-X[4+K]],[X[2+K],X[1+K],-X[4+K],X[3+K]],[X[3+K],X[4+K],X[1+K],-X[2+K]],[X[4+K],-X[3+K],X[2+K],X[1+K]]]); EndFor; Return Q; EndDefine; -- FueterMat LeftFueterMat(X):=Coala.FueterMat(X); -- an alias of the left CF operator Define RightFueterMat(X); -- I: a LIST of variables of length 4N -- O: the MATrix symbol of the right CF operator in N quaternionic variables N:=Len(X)/4; Q:=[]; For J:=1 To N Do K:=4*(J-1); Q:=Coala.MatConcat(Q,Mat[[X[1+K],-X[2+K],-X[3+K],-X[4+K]],[X[2+K],X[1+K],X[4+K],-X[3+K]],[X[3+K],-X[4+K],X[1+K],X[2+K]],[X[4+K],X[3+K],-X[2+K],X[1+K]]]); EndFor; Return Q; EndDefine; -- RightFueterMat ----[ Complex matrix ]---- Define ComplexMat(X); -- left CF matrix written with complex notation N:=Len(X)/4; Q:=[]; For J:=1 To N Do K:=4*(J-1); Q:=MatConcat(Q,Mat[[X[1+K],X[2+K]],[X[3+K],X[4+K]]]); EndFor; Return Q; EndDefine; -------------- DIRAC OPERATOR -------------------------- ----[ auxiliary functions for the Clifford product ]---- -- an element E_{i1...ik} is represented by a list E of a sign E[1] -- and a multiindex E[2]=[i1,...,ik] Define Rearrange(E); -- take a basis element and rewrite it such that -- i1<...A[K+1] Then T:=A[K]; A[K]:=A[K+1]; A[K+1]:=T; S:=-S; K:=K+1; Else K:=K+1; EndIf; If OrdA=A Then Fine:=TRUE; Else If K=L Then K:=1; EndIf; EndIf; Until Fine; NewE:=Coala.KillPairs(A,S); Return NewE; EndDefine; Define KillPairs(A,S); DistrA:=Distrib(A); NewA:=[]; ForEach P In DistrA Do If Mod(P[2],4)=0 Then Skip; Elsif Mod(P[2],4)=1 Then Append(NewA,[P[1],1]); Elsif Mod(P[2],4)=2 Then S:=-S; Else Append(NewA,[P[1],1]); S:=-S; EndIf; EndForEach; A:=[P[1]|P In NewA]; Return [S,A]; EndDefine; Define CliffMult(E,F); --Multiplies E and F, basis elements E:=Coala.Rearrange(E); F:=Coala.Rearrange(F); S1:=E[1]; S2:=F[1]; S:=S1*S2; E:=E[2]; F:=F[2]; I:=Intersection(E,F); A:=Concat(Diff(E,I),Diff(F,I)); A:=Coala.Ordina(A); U:=Concat(E,F); L:=Len(U); For J:=1 To L-1 Do For K:=J+1 To L Do If U[J]>=U[K] Then S:=-S; EndIf; EndFor; EndFor; Return([S,A]); EndDefine; Define ListMultiVector(N); -- list of all the multivectors E_{A}, with A in S_n A:=[[]]; For I:=1 To N Do ForEach T In Subsets(1..N,I) Do Append(A,T); EndForEach; EndFor; Return A; EndDefine; Define IndexMat(N); --- this is the function that decides what basis to choose in C_n A:=Coala.ListMultiVector(N); B:=[]; I:=1; ForEach S In A Do ToApp:=[1,S,I]; Append(B,ToApp); I:=I+1; EndForEach; Return B; EndDefine; ----[ Clifford matrices ]---- Define CliffordSquareMat(X); -- I: a LIST of variables of length N -- O: the square MATrix symbol of the Dirac operator corresponding to one Clifford variable N:=Len(X); Rows:=Coala.IndexMat(N); Cols:=Rows; NN:=2^N; M:=NewMat(NN,NN,0); ForEach R In Rows Do ForEach C In Cols Do I:=R[3]; J:=C[3]; If I>J Then A:=[R[1],R[2]]; B:=[C[1],C[2]]; LA:=Len(A[2]); LB:=Len(B[2]); If Abs(LA-LB)=1 Then B:=[B[1]*((-1)^(LB+1)),Reversed(B[2])]; S:=Coala.CliffMult(A,B); Z:=1; If Len(S[2])>1 Then Z:=0; EndIf; IndVar:=S[2,1]; SignVar:=S[1]; M[I,J]:=Z*SignVar*X[IndVar]; EndIf; EndIf; EndForEach; EndForEach; Return Transposed(M)-M; EndDefine; Define CliffordMat(X,K); -- I: a LIST of variables of length K*N and an INTeger K (the number of ops) -- O: the MATrix symbol of the Dirac operator corresponding to K Clifford variables in C_N If Not(IsZero(Mod(Len(X),K))) Then Error('Number of variables and number of Operators don''t match'); EndIf; N:=Div(Len(X),K); A:=[]; For I:=1 To K Do --each step constructs the square matrix for the Ith operator MinIndex:=1+N*(I-1); MaxIndex:=N*I; V:=[X[J]| J In MinIndex..MaxIndex]; M:=Coala.CliffordSquareMat(V); Append(A,M); EndFor; --D:=A[1]; --For I:=2 To K Do -- D:=Coala.MatConcat(D,A[I]); --EndFor; Return Coala.MatConcat(A); EndDefine; --------------[ Rarita Schwinger Operator ]-------------------- -- The RS operator is represented by a n2^n by n2^n matrix whose diagonal block are -- Dirac like operators and on the off diagonal blocks there are some variations -- It is again a skew-symmetric operator so we just need the upper diagonal blocks -- the strategy to follow will be to build those blocks and then "concat" them using BlockMatrix Define RaritaSquareMat(X); --I: A LIST of INDETerminateS whose length is the dimension N of the space --O: the symbol MATrix of the Rarita-Schwinger operator with NO spinor reduction N:=Len(X); RS:=NewMat(N,N); -- the BLOCK matrix -- -- The multiindices A and B are used in each block for columns and rows respectively -- -- the index K=1..N indicate rows in the block matrix, while S indicate columns -- For K:=1 To N Do For S:=K To N Do If K=S Then RS[K,S]:=Coala.CliffordSquareMat(X); RS[K,S]:=((N-1)/N)*RS[K,S]; -- a coefficient coming from the formula Else -- S>K Rows:=Coala.IndexMat(N); Cols:=Rows; NN:=2^N; M:=NewMat(NN,NN,0); ForEach R In Rows Do ForEach C In Cols Do I:=R[3];--I; J:=C[3];--J; If I>J Then EA:=[R[1],R[2]];--EA; EB:=[C[1],C[2]];--EB; EK:=[1,[K]];--EK; ES:=[1,[S]];--ES; LA:=Len(EA[2]); LB:=Len(EB[2]); -- the actual elements to multiply are V:=Coala.CliffMult([(-1)^LA,Reversed(EA[2])],ES);--V; W:=Coala.CliffMult(EK,EB);--W; LV:=Len(V[2]); LW:=Len(W[2]); If Abs(LV-LW)=1 Then VW:=Coala.CliffMult(V,W);--VW; Z:=1; If Len(VW[2])>1 Then Z:=0; EndIf; IndVar:=VW[2,1]; SignVar:=VW[1]; M[I,J]:=Z*SignVar*X[IndVar];--M[I,J]; EndIf; EndIf; EndForEach; EndForEach; M:=Transposed(M)-M; RS[K,S]:=(1/N)*M; RS[S,K]:=-Transposed(RS[K,S]); EndIf; EndFor; EndFor; Return BlockMatrix(RS); EndDefine; Define RaritaSchwingerMat(X,K); -- I: a LIST of variables of length K*N and an INTeger K (the number of ops) -- O: the MATrix symbol of the RS operator(s) If Not(IsZero(Mod(Len(X),K))) Then Error('Number of variables and number of Operators don''t match'); EndIf; If K=1 Then Return Coala.RaritaSquareMat(X); EndIf; N:=Div(Len(X),K); A:=[]; For I:=1 To K Do --each step constructs the square matrix for the Ith operator MinIndex:=1+N*(I-1); MaxIndex:=N*I; V:=[X[J]| J In MinIndex..MaxIndex]; M:=Coala.RaritaSquareMat(V); Append(A,M); EndFor; Return Coala.MatConcat(A); EndDefine; Rarita(X,K):=Coala.RaritaSchwingerMat(X,K); RaritaSchwinger(X,K):=Coala.RaritaSchwingerMat(X,K); RaritaMat(X,K):=Coala.RaritaSchwingerMat(X,K); RaritaSchw(X,K):=Coala.RaritaSchwingerMat(X,K); RaritaSchwMat(X,K):=Coala.RaritaSchwingerMat(X,K); Define MatToSingular(M); -- prints a string that can be copied to define the mat M on Singular, supposing that the -- ring is Q[x[1..N]] on CoCoA, i.e. r=0,(x(1..N)),dp on Singular X:=Indets(); N:=NumIndets(); PrintLn;PrintLn; I:=Len(M); J:=Len(M[1]); Print('matrix m'); Print('[',I,']','[',J,']='); For A:=1 To I Do For B:=1 To J Do El:=M[A,B]; If IsZero(El) Then V:='0'; Elsif Sign(El)=1 Then P:=Coala.PositionIndet(El); V:='x('+P+')'; Else P:=Coala.PositionIndet(-El); V:='-x('+P+')'; EndIf; Print(V); If A<>I Or B<>J Then Print(','); EndIf; EndFor; EndFor; Print(';'); PrintLn(); PrintLn(); Return; EndDefine; Define PositionIndet(T); N:=NumIndets(); X:=Indets(); For I:=1 To N Do; If T=X[I] Then If I<10 Then P:=Ascii(48+I); Else --- supposing I<100 for now Pdec:=Ascii(48+Div(I,10)); Pun:=Ascii(48+Mod(I,10)); P:=Pdec+Pun; EndIf; Break; EndIf; EndFor; Return P; EndDefine; -----------------[ Noetherian Operators ]-------------------- ------- MAIN FUNCTIONS ------- Define NoetherianOps(Q); -- redirect the input to the various routines X:=Indets(); N:=NumIndets(); If Type(Q)=MODULE Then G:=Q.Gens; S:=Len(G[1]); Dimension:=Dim(CurrentRing()^S/Q); If Not(IsZero(Dimension)) Then Error('Zerodimensional Module Required'); EndIf; PrintLn(' '); PrintLn('Assuming that the input is primary and centered at the origin'); PrintLn(' '); Return [Coala.NoetherianED_Module(Q,0),X]; Elsif Type(Q)=IDEAL Then Dimension:=Dim(CurrentRing()/Q); PrintLn(' '); PrintLn('Assuming that the input is primary'); PrintLn(' '); Normalization:=Coala.IdealNormalPos(Q); IdN:=Identity(N); If Normalization[3]=IdN Then PrintLn(' '); PrintLn('The ideal is already in Noether Normal Position with respect to ',Normalization[2]); PrintLn(' '); Else PrintLn(' '); PrintLn('A change of coordinates has been performed on the ideal. The change is represented by the matrix'); PrintLn(Normalization[3]); PrintLn(' '); PrintLn('The ideal is now in Noether Normal Position with respect to ',Normalization[2]); PrintLn(' '); PrintLn('The new expression of the ideal is ',Normalization[1]); Q:=Normalization[1]; Endif; If IsZero(Dimension) Then If Coala.IsMonomialIdeal(Q) Then Return [QuotientBasis(Q),X]; EndIf; RI:=Radical(Q); M:=Ideal(Indets()); If Not(RI=M) Then Error('The zerodimensional primary ideal is not centered at the origin'); EndIf; Return [Coala.InvertList(Coala.SortedPoly(Coala.NoetherianED(Q,0))),X]; Else T:=Diff(X,Normalization[2]); X:=Normalization[2]; PrintLn(' '); PrintLn('Assuming that the variety associated to the extended ideal in K(',T,')',X,' is the origin'); PrintLn(' '); Return [Coala.InvertList(Coala.SortedPoly(Coala.NoetherianHighIdeal(Q,T,0))),X]; EndIf; Else Error('Ideal or Zerodimensional Module required'); EndIf; Return ; EndDefine; Define NoetherianOpsMMM(Q); If Not(IsZero(Dim(CurrentRing()/I))) Then Error('Zerodimensional Ideal Required'); EndIf; If Type(Q)=IDEAL Then PrintLn(' '); PrintLn('Assuming that the input is primary'); PrintLn(' '); RI:=Radical(Q); M:=Ideal(Indets()); If Not(RI=M) Then Error('The zerodimensional primary ideal is not centered at the origin'); EndIf; Return [Coala.NoetherianMMM(Q,0),Indets()]; Else Error('Zerodimensional Ideal Required'); EndIf; Return ; EndDefine; Define NoetherNormal(I); -- an alias Return Coala.IdealNormalPos(I); EndDefine; ------- PRINT FUNCTIONS -------- Define NoetherianPrint(O); -- prints the noetherian ops as differentials. The input is the output of any main function L:=O[1]; -- list of ops X:=O[2]; -- differentials Count:=0; --L:=Coala.ToDerivative(L); ForEach D In L Do --- print the operator If Coala.Grado(D)=0 Then Print(D); Else If Type(D)=VECTOR Then Print('Vector('); CountV:=0; ForEach V In D Do CountV:=CountV+1; OperatorPrint(V,X); If CountV1 Then Print(TPart); EndIf; Coala.DiffPrint(XPart,X); EndDefine; Define MonomialPrint(T,X); If IsZero(Coala.Grado(T)) Then If LC(T)>=0 Then Print('+'); EndIf; Print(T); Else Coeff:=LC(T); Ter:=T/Coeff; If Coeff<0 Then If IsZero(Coeff+1) Then Print('-'); Else Print(Coeff); EndIf; Else Print('+'); If Coeff<>1 Then Print(Coeff);EndIf; EndIf; Par:=Diff(Indets(),X); -- parameters ToSubst:=[[P,1]|P In Par]; XPart:=Subst(Ter,ToSubst); TPart:=Ter/XPart; If TPart<>1 Then Print(TPart);EndIf; Coala.DiffPrint(XPart,X); EndIf; EndDefine; Define DiffPrint(XPart,X); ForEach V In X Do If Der(XPart,V)<>0 Then Print('d',V); N:=Coala.LogVar(XPart,V); -- this is at least one If N>1 Then Print('^',N);EndIf; EndIf; EndForEach; EndDefine; Define LogVar(XPart,V); N:=-1; D:=XPart; While D<>0 Do; D:=Der(D,V); N:=N+1; EndWhile; Return N; EndDefine; --------- NORMAL POSITION ---------- --- Term ordering assumed to be degree-compatible Define IdealNormalPos(I); -- I: an IDEAL I -- O: a LIST of the following elements: -- the IDEAL I put in normal position with noether normalization -- the list of variables wrt which I is now in normal position -- the matrix representing the linear change of coordinates -- the ideal generated by the monic generators x_i^D + tail from the normalization lemma D:=Dim(CurrentRing()/I); X:=Indets(); Parameters:=[]; L:=Len(I.Gens); N:=NumIndets(); A:=Identity(N); NewI:=I; Collect:=[]; MonicFound:=FALSE; J:=1; VarToChange:=1..NumIndets(); While Not(IsZero(NewI)) And J<100 Do MonicFound:=FALSE; G:=NewI.Gens; ResultMonic:=Coala.IsThereMonic(N,G,MonicFound,X,NewI,VarToChange); If ResultMonic[1] Then MonicFound:=TRUE; XK:=ResultMonic[2]; V:=ResultMonic[3]; K:=ResultMonic[4]; Append(Parameters,XK); Append(Collect,V); NewI:=Elim(XK,NewI); VarToChange:=Diff(VarToChange,[K]); -- PrintLn(Parameters); EndIf; J:=J+1; If Not(MonicFound) Then A:=Coala.ChangeVar(A,VarToChange); -- PrintLn(A); NewI:=Coala.ApplyChange(NewI,A,VarToChange); I:=Coala.ApplyChange(I,A,VarToChange); -- PrintLn(NewI); EndIf; EndWhile; --EndFor; Return [I,Reversed(Sorted(Parameters)),A,Ideal(Collect)]; EndDefine; Define Monico(T,K); -- true if the monomial T is a power of X_K X:=Indets(); F:=Factor(T); If Len(F)=1 And X[K]=LT(F[1,1]) Then Return TRUE; EndIf; Return FALSE; EndDefine; Define ChangeVar(A,V); -- I: a MATrix representing a change of vars -- the LIST of the indices of the variables that will not be fixed -- O: a new MATrix representing a different change of coordinates, perturbing the vars in V N:=Len(A); NewA:=A; ForEach I In V Do ForEach J In V Do NewA[I,J]:=Coala.PerturbElement(A,I,J); EndForEach; EndForEach; Return NewA; EndDefine; Define PerturbElement(A,I,J); If J>=I Then Return A[I,J]; EndIf; --- still lower triangular B:=Coala.Sopra(A,I,J); C:=Coala.Davanti(A,I,J); Return A[I,J]+Rand(-2,2); EndDefine; Define Sopra(A,I,J); N:=Len(A); NewI:=Mod(I-1,N); If IsZero(NewI) Then NewI:=N;EndIf; Return A[NewI,J]; EndDefine; Define Davanti(A,I,J); N:=Len(A); NewJ:=Mod(J+1,N); If IsZero(NewJ) Then NewJ:=N;EndIf; Return A[I,NewJ]; EndDefine; Define ApplyChange(I,A,VarToChange); N:=Len(A); Mappa:=NewList(N); X:=Indets(); For J:=1 To N Do Mappa[J]:=Sum([(A[J,K]*X[K])|K In 1..N]); EndFor; --Mappa; ToSub:=[[X[J],Mappa[J]]|J In 1..N]; --ToSub; Return Subst(I,ToSub); EndDefine; Define IsThereMonic(N,G,MonicFound,X,NewI,VarToChange); ForEach K In VarToChange Do -- K; ForEach V In G Do If Coala.Monico(LT(V),K) Then -- PrintLn(LT(V)); Return ([TRUE,X[K],V,K]); EndIf; EndForEach; EndForEach; Return([FALSE,X[N],First(G),N]); EndDefine; -------- Sigma, Ro, Lambda ------------ -- by DERIVATIVE I mean a sum of elements of the type d/dx_i etc -- by DIFFERENTIAL I mean e sum of elements of the type D(i_1,...,i_n) like in [MMM] Define FactorialList(L) -- for L=[a,b] returns a! * b! Return Product([Fact(A)|A In L]); EndDefine; Define ToDerivative(F); -- transforms a differential into a derivative ---- If it's a list or vector then do it for every element -- If Type(F)=LIST Then Return [Coala.ToDerivative(V)|V In F]; Elsif Type(F)=VECTOR Then Return Vector([Coala.ToDerivative(V)|V In F]); Else -- M:=Monomials(Cast(F,POLY)); NewF:=0; Foreach T In M Do L:=Log(T); I:=Coala.FactorialList(L); NewF:=NewF+(1/I)*T; EndForEach; EndIf; Return NewF; EndDefine; Define ToDifferential(F); -- inverse of ToDerivative ---- If it's a list then do it for every element -- If Type(F)=LIST Then Return [Coala.ToDifferential(V)|V In F];EndIf; -- --- M:=Monomials(Cast(F,POLY)); NewF:=0; Foreach T In M Do L:=Log(T); I:=Coala.FactorialList(L); NewF:=NewF+(I)*T; EndForEach; Return NewF; EndDefine; ------------- sigma ------------- Define DerivationByTerm(T,V); -- derives a term V with respect to a term T: d/dT (V) If IsZero(Coala.Grado(T)) Then Return T*V; EndIf; X:=Indets(); N:=NumIndets(); L:=Log(T); C:=LC(T); For I:=1 To N Do; If Not(IsZero(L[I])) Then V:=Coala.DerivationByIndet(X[I],L[I],V); EndIf; EndFor; Return C*V; EndDefine; Define DerivationByIndet(X,M,V); -- derives M>0 times a term V with respect to an indeterminate X: d^M/dX^M (V) For I:=1 To M Do If IsZero(V) Then Break EndIf; V:=Der(V,X); EndFor; Return V; EndDefine; Define Sigma(...) -- Sigma(T,F,[optional]) is sigma_t(F) where F is a DERIVATIVE. -- If the optional variable is set, F is intended to be a DIFFERENTIAL and -- it returns a DIFFERENTIAL T:=ARGV[1]; F:=ARGV[2]; If IsZero(F) Then Return(F); EndIf; If Len(ARGV)=3 Then F:=Coala.ToDerivative(F); EndIf; NewF:=Sum([Coala.DerivationByTerm(T,V)|V In Monomials(F)]); If Len(ARGV)=3 Then NewF:=Coala.ToDifferential(NewF); EndIf; Return NewF; EndDefine; --------------- Ro --------------------- Define IntegrationByTerm(T,V); -- integrates a term V with respect to a term T: d/dT (V) X:=Indets(); N:=NumIndets(); L:=Log(T); C:=LC(T); For I:=1 To N Do; If Not(IsZero(L[I])) Then V:=Coala.IntegrationByIndet(I,L[I],V); EndIf; EndFor; Return C*V; EndDefine; Define IntegrationByIndet(I,M,V); -- integrates M>0 times a term V with respect to the indeterminate X[I] X:=Indets(); For J:=1 To M Do L:=Log(V); C:=L[I]; V:=1/(C+1)*X[I]*V; EndFor; Return V; EndDefine; Define Ro(...); -- Ro(T,F,[optional]) is ro_t(F) where F is a DERIVATIVE. -- If the optional variable is set, F is intended to be a DIFFERENTIAL and -- it returns a DIFFERENTIAL T:=ARGV[1]; F:=ARGV[2]; If IsZero(F) Then Return(F); EndIf; If Deg(F)=0 Then Return (F*T); EndIf; If Len(ARGV)=3 Then F:=Coala.ToDerivative(F); EndIf; NewF:=Sum([Coala.IntegrationByTerm(T,V)|V In Monomials(F)]); If Len(ARGV)=3 Then NewF:=Coala.ToDifferential(NewF); EndIf; Return NewF; EndDefine; --------------- Lambda = Ro(Sigma()) --------- Define Lambda(...); -- lambda_T(F) where F is a DERIVATIVE. -- If the optional variable is set, F is intended to be a DIFFERENTIAL and -- it returns a DIFFERENTIAL T:=ARGV[1]; F:=ARGV[2]; If Len(ARGV)=3 Then Return Coala.Ro(T,Coala.Sigma(T,F,'diff'),'erential'); EndIf; Return Coala.Ro(T,Coala.Sigma(T,F)); EndDefine; ----- MAIN functions for noetherian operators ----- -- [Zerodimensional ideal, Standard approach] -- -- To be implemented in a next version : Primary Decomposition Of an Ideal I into Q1..Qs -- For each Primary Ideal Q check if it is in normal position vith Vasconcelos' test -- This can be done inside the primary decomposition. -- Noether Normalization and change of coordinates: -- Find a subset of variables x1..xk in x1..xn such that the image of Q is in normal position -- and zerodimensional in the new ring K(x_1..x_k)[x_k+1 ... x_n] -- Perform a change of coordinates (and a field extension) such that V(Q)={(0,...,0)} Define NoetherianED(Q,Flag); If Flag=1 Then Talk:=TRUE; -- The Flag makes the routine verbose Else Talk:=FALSE; EndIf; -- IN: We assume that the ideal Q is primary, zerodimensional, and centered in zero If Dim(CurrentRing()/Q)>0 Then Error('ZeroDimensional Ideal Required'); EndIf; -- STEP: let M be the multiplicity of the ideal Q, i.e. the multiplicity of (0,...,0) as a point M:=Multiplicity(CurrentRing()/Q); If Talk Then PrintLn('Multiplicity of the point: ',M); EndIf; -- STEP: construct the germ of a holomorphic function H(x1..xn) such that -- it represents regular functions over Q. Coefficients are not specified -- Define a new ring with the coefficients of H as new variables NX:=NumIndets(); -- number of old unknowns X:=Indets(); -- old unkonwns MAX:=Ideal(X); -- Find the possible terms for the Taylor Series of H T:=Sum([Ideal(MAX^J)|J In 0..(M-1)]); -- the degree of a germ is at most the multiplicity! NC:=Len(T); S::=CoeffRing[x[1..NX],c[1..NC]]; --Do I need to use LEX? not here maybe Using S Do SIndets:=Indets(); SIndetsX:=[SIndets[I]|I In 1..NX]; If Talk Then PrintLn('Old Indeterminates: ',SIndetsX); EndIf; SIndetsC:=[SIndets[I]|I In (NX+1)..(NX+NC)]; ST:=Image(T,RMap(x[1]..x[NX])); SQ:=Image(Q,RMap(x[1]..x[NX])); SMAX:=Image(MAX,RMap(x[1]..x[NX])); -- Correspondence Table between the coefficients c and the terms in T CT:=[]; For I:=1 To NC Do ToApp:=[SIndetsC[I],ST.Gens[I]]; Append(CT,ToApp); EndFor; If Talk Then PrintLn(NewLine); PrintLn('Coefficients and Series Terms',CT); EndIf; -- Finally H H:=Sum([SIndetsC[I]*ST.Gens[I]|I In 1..NC]); If Talk Then PrintLn(NewLine); PrintLn('Germ of the Holomorphic Function',NewLine,H); EndIf; -- STEP: Compute a Groebner Bases of Q wrt any ordering G:=GBasis(SQ); --- (May be Skipped) -- STEP: Reduce H to his class modulo Q using Normal Form, leave coefficients as they are NFH:=NF(H,SQ); If Talk Then PrintLn(NewLine); PrintLn('Reduced Germ of the Holomorphic Function: ',NewLine,NFH); EndIf; -- FACT: Annihilation of the coefficients of NF(H) wrt the indeterminates x1..xn give noetherian operators: -- list of the terms in the x1..xn (generators of T) -- XTermList:=ST.Gens; -- Coefficients SC:=[]; -------- ADDING HERE SOME SPEED-UP TRICK ----------- MAXN:=Deg(LT(NFH))-1; -- This is the degree of the highest residual monomial in NF(h) NewSTGens:=[V|V In ST.Gens And Deg(V)<=MAXN]; -------- END OF TRICK ------------------------------ ToSubst:=[[x[I],0]|I In 1..NX]; -- ForEach V In ST.Gens ForEach V In NewSTGens Do C:=Coala.NoetherPolyPoly(V,NFH,1); C:=Subst(C,ToSubst); -- only the coefficients without x Append(SC,C); EndForEach; If Talk Then PrintLn(NewLine); PrintLn('c Coefficients of the Holomorphic Function',NewLine,SC); EndIf; -- STEP: transform each coefficient into an operator using the Taylor series formula. SC:=Subst(SC,CT); EndUsing; --Indets(); ListImage:=Concat(Indets(),NewList(NC,1)); C:=Image(SC,RMap(ListImage)); --C; -- Get rid of some zeroes !! C:=Coala.KillZero(C); If Talk Then PrintLn(NewLine); PrintLn('Dual of the operators using the original variables: '); EndIf; -- OUT: Print the operators nicely with dx1,...,dxn and also return some algebraic output (polys) Return C; EndDefine; Define KillZero(L); -- Eliminates the zero elements form a list L of any object NewL:=[]; ForEach V In L Do If Not(IsZero(V)) Then Append(NewL,V); EndIf; EndForEach; Return NewL; EndDefine; --- [ Zerodimensional Module, Standard Approach ] --- -- To be implemented in a next version : Primary Decomposition Of a Module into Q1..Qs -- For each Primary Component Q check if it is in normal position with Vasconcelos' test -- This can be done inside the primary decomposition. -- Noether Normalization and change of coordinates: -- Find a subset of variables x1..xk in x1..xn such that the image of Q is in normal position -- and zerodimensional in the new free R-module with ring R=K(x_1..x_k)[x_k+1 ... x_n] -- Perform a change of coordinates (and a field extension) such that V(Q)={(0,...,0)} Define NoetherianED_Module(Q,Flag); If Flag=1 Then Talk:=TRUE; -- The Flag makes the routine verbose Else Talk:=FALSE; EndIf; -- IN: We assume that the module Q is primary, zerodimensional, and centered in zero -- it is a submodule of R^t so we can define his dimension T:=NumComps(Q); If Dim(CurrentRing()^T/Q)>0 Then Error('ZeroDimensional Module Required'); EndIf; -- STEP: let M be the multiplicity of the ideal Q, i.e. the multiplicity of (0,...,0) as a point of -- the characteristic variety M:=Multiplicity(CurrentRing()^T/Q); If Talk Then PrintLn('Multiplicity of the zero: ',M); EndIf; -- STEP: construct the germ of a holomorphic function H(x1..xn) such that -- it represents regular functions over Q. Coefficients are not specified -- Define a new ring with the coefficients of H as new variables NX:=NumIndets(); -- number of old unknowns X:=Indets(); -- old unkonwns MAX:=Ideal(X); -- Find the possible terms for the Taylor Series of H TERMS:=Sum([Ideal(MAX^J)|J In 0..(M-1)]); -- the degree of a germ is at most the multiplicity! NC:=Len(TERMS); -- the module terms are of the form te_i with t in TERMS and i=1..T S::=CoeffRing[x[1..NX],c[1..T,1..NC]]; --Do I need to use LEX? Using S Do SIndets:=Indets(); SIndetsX:=[SIndets[I]|I In 1..NX]; If Talk Then PrintLn('Old Indeterminates: ',SIndetsX); EndIf; SIndetsC:=[SIndets[I]|I In (NX+1)..(NX+(NC*T))];--SIndetsC; ST:=Image(TERMS,RMap(x[1]..x[NX])); SQ:=Image(Q,RMap(x[1]..x[NX])); SMAX:=Image(MAX,RMap(x[1]..x[NX])); -- Correspondence Table between the coefficients c and the terms in T CT:=[]; For J:=1 To T Do For I:=1 To NC Do ToApp:=[SIndetsC[I+NC*(J-1)],ST.Gens[I]*E_(J,T)]; Append(CT,ToApp); EndFor; EndFor; If Talk Then PrintLn(NewLine); PrintLn('Coefficients and Series Terms',CT); EndIf; -- Finally H H:=Sum([CT[I][1]*CT[I][2]|I In 1..Len(CT)]); If Talk Then PrintLn(NewLine); PrintLn('Germ of the Holomorphic Function',NewLine,H); EndIf; -- STEP: Compute a Groebner Bases of Q wrt any ordering G:=GBasis(SQ); --- (May be Skipped) -- STEP: Reduce H to his class modulo Q using Normal Form, leave coefficients as they are NFH:=NF(H,SQ); If Talk Then PrintLn(NewLine); PrintLn('Reduced Germ of the Holomorphic Function: ',NewLine,NFH); EndIf; -- FACT: Annihilation of the coefficients of NF(H) wrt the Module Terms tej give noetherian operators: -- list of the terms in the x1..xn and E_i (generators of T) XEiTermList:=[CT[I][2]|I In 1..Len(CT)]; -- Coefficients SC:=[];--XEiTermList;PrintLN; ToSubst:=[[x[I],0]|I In 1..NX];--ToSubst; ForEach V In XEiTermList Do C:=Coala.NoetherPolyPoly_Module(V,NFH,1);--C; --this performs the derivation if the component is the same C:=Subst(C,ToSubst);--C; -- only the coefficients without x, but the component ej should remain Append(SC,C); EndForEach; If Talk Then PrintLn(NewLine); PrintLn('c Coefficients of the Holomorphic Function',NewLine,SC); EndIf; -- STEP: transform each coefficient into an operator using the Taylor series formula. -- warning the Correspondence Table has to be used wisely, and Subst does not work in this case -- SC:=Subst(SC,CT);PrintLn;SC; NewSC:=[]; ForEach C In SC Do NewC:=NewVector(NumComps(SQ)); ForEach X In Monomials(C) Do NewC:=NewC+LC(X)*Coala.CorrespondingOper(X/LC(X),CT); -- based on the CT EndForEach; Append(NewSC,NewC); EndForEach; SC:=NewSC; EndUsing; --Indets(); ListImage:=Concat(Indets(),NewList(NC*T,1)); SC:=Image(SC,RMap(ListImage)); --SC; -- Get rid of some zeroes !! SC:=Coala.KillZero(SC); If Talk Then PrintLn(NewLine); PrintLn('Dual operators using the original variables: ',SC); EndIf; -- OUT: Print the operators nicely with dx1,...,dxn and also return some algebraic output (polys) Return SC; EndDefine; --- [ Positive Dimensional Ideal: Standard Approach ] --- -- To be implemented in a next version : Primary Decomposition Of an Ideal I into Q1..Qs -- For each Primary Ideal Q check if it is in normal position vith Vasconcelos' test -- This can be done inside the primary decomposition. ----------------------------------------------------------------- -- Noether Normalization and change of coordinates: -- We ASSUME that the ideal is already pplace in Normal position with respect to the variables not in T=[t1..tk] -- We also assume that in Q(t1..tk)[x1..xn-k] the ideal is zeromensional, primary and that the variety is -- The origin, meaning that for every choice of the parameters t the only solution of the poly system given by I -- is just (0..0) in Q^(n-k) Define NoetherianHighIdeal(Q,T,Flag); If Flag=1 Then Talk:=TRUE; -- The Flag makes the routine verbose Else Talk:=FALSE; EndIf; -- STEP: let M be the multiplicity of the ideal Q, i.e. the multiplicity of (0,...,0) as a point M:=Multiplicity(CurrentRing()/Q); If Talk Then PrintLn('Multiplicity of the Variety: ',M); EndIf; -- STEP: construct the germ of a holomorphic function H_t(x) such that -- it represents regular functions over Q in K(t)[x]. Coefficients are not specified -- Define a new ring with the coefficients of H as new variables ALL:=Indets(); -- old unkonwns X:=Diff(ALL,T); -- the coordinates x1..xn-k MAX:=Ideal(X); NX:=Len(X); -- it corresponds to n-k NT:=Len(T); -- this is k -- Find the possible terms for the Taylor Series of H TAY:=Sum([Ideal(MAX^J)|J In 0..(M-1)]); -- the degree of a germ is at most the multiplicity! NC:=Len(TAY); S::=CoeffRing[x[1..NX],t[1..NT],c[1..NC]],Lex; --Do I need to use LEX? Or a product order? Using S Do SIndets:=Indets(); SIndetsX:=[SIndets[I]|I In 1..NX]; SIndetsT:=[SIndets[I]|I In (NX+1)..(NX+NT)]; If Talk Then PrintLn('Coordinates: ',SIndetsX); PrintLn('Parameters: ',SIndetsT); EndIf; SIndetsC:=[SIndets[I]|I In (NX+NT+1)..(NX+NT+NC)]; ToMap:=Concat(SIndetsX,SIndetsT); STAY:=Image(TAY,RMap(ToMap)); SQ:=Image(Q,RMap(ToMap)); SMAX:=Image(MAX,RMap(ToMap)); -- Correspondence Table between the coefficients c and the terms in T CT:=[]; For I:=1 To NC Do ToApp:=[SIndetsC[I],STAY.Gens[I]]; Append(CT,ToApp); EndFor; If Talk Then PrintLn(NewLine); PrintLn('Coefficients and Series Terms',CT); EndIf; -- Finally H H:=Sum([SIndetsC[I]*STAY.Gens[I]|I In 1..NC]); If Talk Then PrintLn(NewLine); PrintLn('Germ of the Holomorphic Function',NewLine,H); EndIf; -- STEP: Compute a Groebner Basis of Q wrt the ordering defined for this new ring G:=GBasis(SQ); --- may be skipped if the input is already a GB, but i think not in this case since the ordering changed ------------ HERE IS THE DIFFERENCE BETWEEN THE ZERODIMENSIONAL CASE AND THE HIGER DIMENSIONAL CASE ------------------------ -- Find the opportune term in t1..tn-k to multiply to H XToOne:=[[IND,1]|IND In SIndetsX]; TTerms:=[Subst(LT(V),XToOne)|V In G]; ToMult:=Product(TTerms); Control:=0; Repeat H:=H*ToMult; Control:=Control+1; --Control; -- STEP: Reduce H to his class modulo Q using Normal Form, leave coefficients as they are NFH:=NF(H,SQ); If Talk Then PrintLn(NewLine); PrintLn('Reduced Germ of the Holomorphic Function: ',NewLine,NFH); EndIf; -- FACT: Annihilation of the coefficients of NF(H) wrt the indeterminates x1..xn give noetherian operators: -- list of the terms in the x1..xn (generators of T) -- XTermList:=ST.Gens; -- Coefficients SC:=[]; -------- ADDING HERE SOME SPEED-UP TRICK ----------- MAXN:=Deg(LT(NFH))-1; -- This is the degree of the highest residual monomial in NF(h) NewSTGens:=[V|V In STAY.Gens And Deg(V)<=MAXN]; -------- END OF TRICK ------------------------------ ToSubst:=[[x[I],0]|I In 1..NX]; -- ForEach V In STAY.Gens ForEach V In NewSTGens Do C:=Coala.NoetherPolyPoly(V,NFH,1); C:=Subst(C,ToSubst); -- only the coefficients without x Append(SC,C); EndForEach; If Talk Then PrintLn(NewLine); PrintLn('c Coefficients of the Holomorphic Function',NewLine,SC); EndIf; -- STEP: transform each coefficient into an operator using the Taylor series formula. SC:=Subst(SC,CT); -- Get rid of some zeroes !! SC:=Coala.KillZero(SC); --PrintLn(' Found: ',SC); Until (Len(SC)=M Or Control=100); SC:=CleanT(SC,SIndetsT); EndUsing; --Indets(); ListImage:=Concat(Indets(),NewList(NC,1)); C:=Image(SC,RMap(ListImage)); --C; If Control=100 Then PrintLn('The algorithm was not able to reduce the taylor expansion in ',Control,' steps'); EndIf; If Talk Then PrintLn(NewLine); PrintLn(Diff(Indets(),T),'-dual of the operators using the original variables: '); EndIf; -- OUT: Print the operators nicely with dx1,...,dxn and also return some algebraic output (polys) Return C; EndDefine; --- [Zerodimensional Ideal: MMM approach ] --- -- To be implemented in a next version : Primary Decomposition Of an Ideal I into Q1..Qs -- For each Primary Ideal Q check if it is in normal position with Vasconcelos' test -- This can be done inside the primary decomposition. -- Noether Normalization and change of coordinates: -- Find a subset of variables x1..xk in x1..xn such that the image of Q is in normal position -- and zerodimensional in the new ring K(x_1..x_k)[x_k+1 ... x_n] -- Perform a change of coordinates (and a field extension) such that V(Q)={(0,...,0)} -------------------------------------------------------------------- -- From now on we suppose that I is Zerodimensional, Primary and that CharVar(I)={(0,...,0)} Define NoetherianMMM(I,Flag); If Flag=1 Then Talk:=TRUE; Else Talk:=FALSE; EndIf; If Not(IsZero(Dim(CurrentRing()/I))) Then Error('ZeroDimensional Ideal Required'); EndIf; --- Start with the identity NOP:=[1]; NX:=NumIndets(); X:=Indets();X; --- compute the multiplicity of R/I, call it MU MU:=Multiplicity(CurrentRing()/I); --- Stop whenever you have found MU noetherian operators If MU=1 Then Return NOP; EndIf; --- construct all the differentials dx1,...,dxn DIFF:=X; -- Remember to use the additional flag "1" for Rho and Sigma --- find the first j s.t. dxj(F)(0,...,0)=0 foreach F in I and add it to the NOP list LI:=Len(I.Gens);PrintLn(LI); JJ:=[]; For J:=1 To NX Do Found:=FALSE; K:=1;PrintLn('J=',J); While IsZero(Eval(Der(I.Gens[K],X[J]),NewList(NX,0))) And Not(Found) And K1 Then LCoeff:=Solut; Found:=TRUE; Break; EndIf; LCoeff:=Last(Sol); Found:=TRUE; EndForEach EndIf; If Not(Found) Then Error('Operator not Found'); EndIf; If Not(ZeroOp) Then -- else it is in the LCoeff; LToSub:=NewList(T,0); K:=1; ForEach J In RightDeg Do LToSub[J]:=LToSub[J]+LCoeff[K]; K:=K+1; EndForEach; PrintLn('Valori da sostituire nelle l[j] ',LToSub); ToSubst:=[[L[J],LToSub[J]]|J In 1..T]; D:=Subst(NewD,ToSubst); PrintLn(' '); PrintLn(' '); PrintLn('New Candidate: ',D); If Coala.Grado(D) foreach j If Found And Not(ZeroOp) Then --- now check condition B Found:=TRUE; For J:=1 To NX Do ----- CLOSURE TEST SJD:=Coala.Sigma(X[J],D,1); If Not(Coala.Closure(NOP,SJD)) Then Found:=FALSE; Break; EndIf; EndFor; EndIf; -- now Found is true if we passed the closure test If Found Then Break; --- Interrupting the ForEach cycle since we found the op EndIf; EndForEach; ---------- FOREACH (Pick) EndUsing; ------- USING ListImage:=Concat(Indets(),NewList(T,1));--ListImage; D:=Image(D,RMap(ListImage)); NOP:=Image(NOP,RMap(ListImage)); --- now that A and B are satisfied stop and add D to NOP Append(NOP,D); EndWhile; ----------------------------------------- END OF WHILE WHILE WHILE LOOP ------------------------------------- Return NOP; EndDefine; ---- CHOICE OF THE SOLUTION -- Define ChooseSolution(Sol); -- Sol is a list of vectors [l1,...,lt]; -- Dimens:=Len(Sol[1]); -- Choice:=NewList(Dimens,0); OneNormList:=[Coala.OneNorm(V)|V In Sol]; M:=Max(OneNormList); ForEach V In Sol Do If OneNorm(V)=M Then Return V; EndIf; EndForEach; Error('A choice of a solution for the linear system could not be made'); EndDefine; --- [ Zerodimensional Ideal: backwards approach ] --- -- AUX Define CornerMonomials(L); -- given a list of monomials find the elements in the corner Corner:=[]; IND:=Indets(); ForEach T In L Do Jump:=[T*X|X In IND]; If Len(Intersection(Jump,L))=0 Then Append(Corner,T); EndIf; EndForEach; Return Corner; EndDefine; Define MonomialDivide(T,S); -- TRUE if T divides S If IsZero(T) Then Return FALSE; EndIf; LogT:=Log(T); LogS:=Log(S); For I:=1 To NumIndets() Do If LogT[I]>LogS[I] Then Return FALSE; EndIf; EndFor; Return TRUE; EndDefine; Define PullBackCorner(T,I);-- given a corner monomial T, find the noetherian operator in which T appears as the ending term I:=Ideal(ReducedGBasis(I));--I; SI:=[[Support(F),F]|F In I.Gens And Len(Support(F))>1];--SI; NewOp:=T; NewTerms:=[]; ForEach S In SI Do EndTerm:=Last(S[1]);--EndTerm; If Coala.MonomialDivide(EndTerm,T) Then -- then we add the head of the polynomial F=S[2] ToMult:=T/EndTerm;--ToMult; Mons:=Monomials(S[2]*ToMult);--Mons; ToDivide:=1/(-LC(Last(Mons)));--ToDivide; MonsToAdd:=Coala.Testa(Mons)*ToDivide;--MonsToAdd; ToSum:=Sum(MonsToAdd); NewTerms:=Concat(NewTerms,Monomials(ToSum)); NewOp:=NewOp+ToSum; PrintLn; PrintLn(NewLine,'The term ',T,' has generated ',ToSum,' using ',S,NewLine); EndIf; EndForEach; MU:=Multiplicity(CurrentRing()/I); While Coala.Grado(NewOp)>=0 Do NewerTerms:=[]; ForEach TT In NewTerms Do ForEach S In SI Do EndTerm:=Last(S[1]);--EndTerm; If Coala.MonomialDivide(EndTerm,TT) Then Stopping:=TRUE; ToMult:=TT/EndTerm;--ToMult; Mons:=Monomials(S[2]*ToMult);--Mons; ToDivide:=1/(-LC(Last(Mons)));--ToDivide; MonsToAdd:=Coala.Testa(Mons)*ToDivide;--MonsToAdd; --PrintLn(MonsToAdd); SuppOp:=Monomials(NewOp); MonsToAdd:=Diff(MonsToAdd,SuppOp);PrintLn('now adding these monomials: ',MonsToAdd); ToSum:=Sum(MonsToAdd);--ToSum; --If Not(IsZero(ToSum)) Then NewerTerms:=Concat(NewerTerms,MonsToAdd); --EndIf; NewOp:=NewOp+ToSum; PrintLn(NewLine,'The term ',TT,' has generated ',ToSum,' using ',S,NewLine); EndIf; EndForEach; EndForEach; NewTerms:=NewerTerms; If NewTerms=[] Then Break;EndIf; EndWhile; PrintLn; ------------ now eliminating the extra terms --------------- NewMons:=Monomials(NewOp); FinalOp:=0; ForEach M In NewMons Do If T=T Then FinalOp:=FinalOp+M; EndIf; EndForEach; Return FinalOp; EndDefine; Define CornerDown(L); -- given the operator L associated to a corner monomial, derives all the ones below it IND:=Indets(); Count:=1; D:=Coala.Grado(L); Max:=Ideal(IND); Derivatives:=Gens(Sum([Max^K|K In 1..D]));--Derivatives; NewOps:=[Coala.Sigma(X,L,1)|X In Derivatives];--NewOps; Append(NewOps,L); Ops:=(Coala.KillZero(Set(NewOps))); Return Coala.InvertList(Coala.SortedPoly2(Ops)); EndDefine; Define NoetherianCorner(I); Q:=QuotientBasis(I);--Q; C:=Coala.CornerMonomials(Q);--C; L:=[]; ForEach M In C Do OP:=Coala.PullBackCorner(M,I);--OP; L:=Concat(L,Coala.CornerDown(OP));--L; EndForEach; Return [Coala.InvertList(Coala.SortedPoly(Set(L))),Indets()]; EndDefine; ---------- MAIN AUX Define CorrespondingOper(X,CT); --the module term in the taylor whose coefficient is X For I:=1 To Len(CT) Do If X=CT[I][1] Then Return CT[I][2]; EndIf; EndFor; Error('Corresponding operator for the coefficient ',X,' not found'); EndDefine; Define CleanT(C,T); --- eliminates any factor of the form t^m with t in T NewC:=[]; ForEach V In C Do Fac:=Factor(V); NewF:=1; ForEach F In Fac Do If Not(IsIn(F[1],T)) Then NewF:=NewF*(F[1]^F[2]); EndIf; EndForEach; Append(NewC,NewF);--NewF; EndForEach; Return NewC; EndDefine; --------------- action of Noether Operator(s) on poly/ideal --------------- Define NoetherAction(...); --- if type of objects is not specified; A:=ARGV[1]; B:=ARGV[2]; If Len(ARGV)=3 Then If Type(A)=POLY Then A:=Coala.ToDerivative(A); Elsif Type(A)=LIST Then A:=[Coala.ToDerivative(D)|D In A]; Else Error('First Argument: POLY or LIST of Noetherian operators expected'); EndIf; EndIf; If Type(A)=LIST Then Return Coala.NoetherSpaceIdeal(A,Ideal(B)); -- in case B is just a poly. If B is an ideal it doesn't affect B Elsif Type(A)=POLY Then If Type(B)=POLY Then Return Coala.NoetherPolyPoly(A,B); Elsif Type(B)=Ideal Then Return Coala.NoetherPolyIdeal(A,B); Else Error('Second Argument: POLY or IDEAL expected'); EndIf; Else Error('First Argument: POLY or LIST of Noetherian operators expected'); EndIf; Return; EndDefine; Define LocalNoetherAction(...); -- same as above but see if the results belongs to The Ideal by calculating NFs If Len(ARGV)=3 Then M:=Coala.NoetherAction(ARGV[1],ARGV[2],1); Else M:=Coala.NoetherAction(ARGV[1],ARGV[2]); EndIf; G:=GBasis(Ideal(ARGV[2])); G:=Ideal(G); If Type(M)=POLY Then Return NF(M,G); Elsif Type(M)=LIST Then Return [NF(P,G)|P In M]; Elsif Type(M)=MAT Then For I:=1 To Len(M) Do For J:=1 To Len(M[1]) Do M[I,J]:=NF(M[I,J],G); EndFor; EndFor; Return M; Else Error('Bad Parameters for Local Noether Action'); EndIf; EndDefine; Define RadicalNoetherAction(...); -- same as above but see if the results vanish on the variety by calculating NFs with respecto to Radical(I) If Len(ARGV)=3 Then M:=Coala.NoetherAction(ARGV[1],ARGV[2],1); Else M:=Coala.NoetherAction(ARGV[1],ARGV[2]); EndIf; RI:=Radical(Ideal(ARGV[2])); G:=GBasis(RI); G:=Ideal(G); If Type(M)=POLY Then Return NF(M,G); Elsif Type(M)=LIST Then Return [NF(P,G)|P In M]; Elsif Type(M)=MAT Then For I:=1 To Len(M) Do For J:=1 To Len(M[1]) Do M[I,J]:=NF(M[I,J],G); EndFor; EndFor; Return M; Else Error('Bad Parameters for Radical Noether Action'); EndIf; EndDefine; Define NoetherPolyPoly(...); -- Action of D on the polynomial F, where D is a Derivative, or a Differential if ARGV[3]<>null. Output is a poly D:=ARGV[1]; F:=ARGV[2]; If Len(ARGV)=3 Then D:=Coala.ToDerivative(D); EndIf; Return Sum([Coala.DerivationByTerm(T,F)|T In Monomials(Cast(D,POLY))]); EndDefine; Define NoetherPolyPoly_Module(...); -- Action of D on the vector F, where D is a Derivative, or a Differential if ARGV[3]<>null. Output is a POLY D:=ARGV[1]; F:=ARGV[2]; If Len(ARGV)=3 Then D:=Coala.ToDerivative(D); EndIf; NComp:=NumComps(F); NewV:=0; For K:=1 To NComp Do S:=Sum([Coala.DerivationByTerm(T[K],F[K])|T In Monomials(D)]); -- S:=S*E_(K,NComp); NewV:=NewV+S; EndFor; Return NewV; EndDefine; Define NoetherPolyIdeal(...); -- Action of D on every generator of the ideal I, same convention as before for Differentials D:=ARGV[1]; I:=ARGV[2]; If Len(ARGV)=3 Then D:=Coala.ToDerivative(D); EndIf; Return [Coala.NoetherPolyPoly(D,F)|F In Gens(I)]; EndDefine; Define NoetherSpaceIdeal(...); --- action of every operator of V (a LIST!!!) on every poly of the ideal I, same thing for differentials. Output is a Mat V:=ARGV[1];-- V is the list of operators I:=ARGV[2]; If Len(ARGV)=3 Then V:=[Coala.ToDerivative(D)|D In V]; EndIf; Return Mat([Coala.NoetherPolyIdeal(D,I)|D In V]); EndDefine; -------------------- CLOSURE TESTING ------------------------- Define ClosedSpace(U); -- TRUE if U is sigma-closed L:=Len(U); X:=Indets(); UNow:=[U[1]]; For I:=2 To L Do PrintLn('I=',I); ToTest:=U[I]; For J:=1 To NumIndets() Do ----- CLOSURE TEST SJD:=Coala.Sigma(X[J],ToTest,1); -- PrintLn('SJD=',SJD); -- PrintLn('J=',J); If Not(Coala.Closure(UNow,SJD)) Then Return FALSE; EndIf; EndFor; Append(UNow,ToTest); -- PrintLn('UNow=',UNow); EndFor; Return TRUE; EndDefine; Define Closure(NOP,SJD); --- TRUE if SJD is in the linear space generated by the elements of NOP, or if it is a constant SortNOP:=Ideal(Reversed(Coala.SortedPoly(NOP)));--SortNOP; L:=Len(SortNOP); X:=Indets(); For K:=1 To L Do S:=GenRepr(SJD,SortNOP);--S; If S=[] Then S:=[X[1]]; EndIf; Gradi:=[Coala.Grado(V)|V In S];--Gradi; If Coala.EachIsIn(Gradi,[0]) Then Return TRUE; EndIf; SortNOP:=Ideal(Tail(SortNOP.Gens)); EndFor; Return FALSE; EndDefine; ------[ Megaforms ]-------- /*-----------------------------------------------------------------*\ | | | Megaforms, 11-22-04 version 3.0 | | | | | | Alberto Damiano | | | \*-----------------------------------------------------------------*/ -- A MegaPoly is a list of Megamonomials -- A MegaMonomial (MegaMon) is a list [D,delta] where ----- D is a polynomial in the megaforms D_i and D_ij ----- delta is a noncommutative expression for the derivatives: [2,1,3] means delta_2 * delta_1 * delta_3 ----- at the end of the list delta there can be a one-element list [i] representing a function (g[i], h[i]....) ------------- LIKE TERMS FUNCTIONS ------------ Define LikeTermsDelta(P); NewP:=[]; D:=[P[J,1]|J In 1..Len(P)]; Delta:=[P[J,2]|J In 1..Len(P)]; Delta:=Set(Delta); ForEach T In Delta Do S:=[P[J,1]|J In 1..Len(P) And P[J,2]=T]; S:=Sum(S); If Not(IsZero(S)) Then Append(NewP,[S,T]); EndIf; EndForEach; Return NewP; EndDefine; Define LikeTermsD(P); --assuming the the functions g ,h etc are at the end of the delta list -- first separate monomials and put the sign "on the function" g NewP:=[]; ForEach V In P Do MonD:=Monomials(V[1]); ForEach M In MonD Do If Sign(M)=1 Then Append(NewP,[M,V[2]]);PrintLn([M,V[2]]); Else LD:=Len(V[2]); -- this is the position of the function g V[2,LD]:=[-V[2,LD,1]]; Append(NewP,[-M,V[2]]);PrintLn([-M,V[2]]); V[2,LD]:=[-V[2,LD,1]]; -- beacuse it does not have to be changed! EndIf; EndForEach; EndForEach; P:=NewP; NewP:=[]; D:=[P[J,1]|J In 1..Len(P)]; Delta:=[P[J,2]|J In 1..Len(P)]; D:=Set(D); ForEach T In D Do S:=[P[J,2]|J In 1..Len(P) And P[J,1]=T]; Append(NewP,[T,S]); EndForEach; Return NewP; EndDefine; ------------------------------------------------------------------------------------- -------------------- DIFFERENTIALS -------------------------------------------------- Define MegaDiff(D,N); -- construct the differential d for the D-th step with N dirac variables; -- -- Define the ring as follows: the d's are the megaforms with 2 indices, the e's with one index -- -- Use R::=Q[d[1..D,1..N,1..N],e[1..D,1..N]]; -- If D=0 Then Return [[e[1,I],[I]]|I In 1..N]; --- e[1,I] has to be interpreted as e^0_I, CoCoA cannot have 0 as index Elsif D=1 Then Ind2:=(1..N)><(1..N); Return [[d[1,I[1],I[2]],[I[1],I[2]]]|I In Ind2]; Else Ind2:=(1..N)><(1..N); Dif1:=[[e[D,I],[I]]|I In 1..N]; Dif2:=[[d[D,I[1],I[2]],[I[1],I[2]]]|I In Ind2]; Return Concat(Dif2,Dif1); EndIf; EndDefine; ------------------------ NORMAL FORMS -------------------------------------------- Define MegaNormalForm(P); -- -- INPUT: P is a MegaPoly -- -- OUTPUT: NewP is the MegaPoly obtained after the reduction of P -- NewP:=[]; ForEach M In P Do NewM:=Coala.RewriteMegaMonomial(M);--- this is a MegaPoly NewP:=Concat(NewP,NewM); EndForEach; Return Coala.LikeTermsDelta(NewP); EndDefine; Define RewriteMegaMonomial(M); -- -- INPUT: M is a MegaMonomial -- -- OUTPUT: NewP is the MegaPoly obtained after the reduction of M -- D:=M[1]; Delta:=M[2]; -- If Type(Last(Delta))=LIST Then --- in this case the delta contains some function g,h,... RealDelta:=[Delta[I]|I In 1..(Len(Delta)-1)]; RealM:=[D,RealDelta]; Optional:=[Last(Delta)]; Else RealM:=M; Optional:=[]; EndIf; ----------------------------- at this point I have in RealM the MegaMonomial I want to reduce -- P:=Coala.RadialReduce(RealM); -- this is a MegaPoly; If Optional=[] Then Return P; Else NewP:=[]; ForEach N In P Do DN:=N[1]; DeltaN:=N[2]; NewN:=[DN,Concat(DeltaN,Optional)]; Append(NewP,NewN); EndForEach; Return NewP; EndIf; Error('Something Wrong with RewriteMegaMonomial'); EndDefine; Define RadialReduce(M); -- -- INPUT: M is a MegaMonomial that does not contain any function in the Delta list -- -- OUTPUT: NewM is the MegaPoly obtained after the reduction of M. If the reduction leads to just one monomial NewM, it's the list [NewM]; -- NewM:=[M]; ------------------------------- USING A REPEAT LOOP TO AVOID RECURSION Repeat Original:=NewM; -- here we store the previous result from a one step reduction NewM:=[]; -- and this is going to be the new one. At the end of the loop, we compare the two as a stopping criterium ForEach N In Original Do P:=Coala.NewStepRadial(N); -- N is a MegaMon, P is a MegaPoly --P:=OneStepRadial(N); NewM:=Concat(NewM,P); EndForEach; Until NewM=Original; -- if this is true, the last one-step reduction did not produce a different poly, so stop. Return NewM; --Return Coala.SortPairs(NewM); -- depending on when the pair sorting is performed EndDefine; Define SortPairs(P); -- -- INPUT: A MegaPoly -- -- OUTPUT: The same MegaPoly where the pairs at the end are sorted decreasingly: [1,1,2,2] becomes [2,2,1,1] -- NewP:=[]; ForEach M In P Do D:=M[1]; Delta:=M[2]; -- NoPairs:=Delta; Pairs:=[]; L:=Len(Delta); For I:=1 To L Step 2 Do If (L-I>0) And Delta[L-I+1]=Delta[L-I] Then Pairs:=Concat(Pairs,[Delta[L-I+1],Delta[L-I]]); Else Break; EndIf; EndFor;--PrintLn(Pairs); K:=L-Len(Pairs); NoPairs:=[Delta[J]|J In 1..K]; Pairs:=Coala.InvertList(Sorted(Pairs)); Append(NewP,[D,Concat(Pairs,NoPairs)]); EndForEach; Return NewP; EndDefine; ------------------------ MODIFIED OneStepRadial ----------------------------- Define NewStepRadial(N); -- -- INPUT: N is a MegaMonomial that does not contain any function in the Delta list -- -- OUTPUT: NewP is the MegaPoly obtained after the collection of pairs in N on the right and just a ONE-STEP reduction of N. -- D:=N[1]; Delta:=N[2]; -- If Len(Delta)<3 Then Return [N]; EndIf; --- no reduction can be performed, but this should never happen! -- -- now we can assume that Delta has at least 3 elements -- NSplit:=Coala.CatchPairs(Delta); NPairs:=NSplit[1]; Other:=NSplit[2]; ------ If Len(Other)<3 Then Return [[D,Concat(NPairs,Other)]]; EndIf; ------- For I:=2 To (Len(Other)-1) Do A:=Delta[Len(Other)-I]; B:=Delta[Len(Other)-I+1]; C:=Delta[Len(Other)-I+2]; If Coala.NewNothingToDo(A,B,C) Then Skip; -- in this case no reduction can be performed so we simply go to the next I Else Heads:=[Delta[J]|J In 1..(Len(Other)-I-1)]; -- note that if I=Len(Other)-1 this is empty Tails:=[Delta[J]|J In (Len(Other)-I+3)..Len(Other)]; -- note that if I<3 (i.e. I=2) this is empty -------- -------- now assuming that the triple ABC has to be reduced and that A<>B<>C<>A -------- If Not(A>B And A>C) Then Error('OneStepRadial: a case not foreseen!');EndIf; -- and the only case left here should be the "real" reduction with the radial condition: AB and B<>C Return (A