/* ZPrimDecGTZ Procedures Author: Luis David Garcia Date: 7 Nov 2002 CoCoA Version: 4.2 This algorithm works for characteristic zero. It might work For "big enough" positive characteristic. But not with this implementation. The PrimaryCheck needs To be fixed so that big prime characteristic is supported. It computes the primary decomposition of a zero dimensional ideal. */ Define LinMap(Flag) L := Indets(); N := NumIndets()-1; If Flag = 0 Then J := Rand(1,N); P := Indet(J); Elsif Flag = 1 Then P := Sum([Rand(0,1)*L[J] | J In 1..N]); Else P := Sum([Rand(1,5)*L[J] | J In 1..N]); EndIf; P := Indet(N+1) + P; Return P; End; -- LinearCoordChange Define LinInv(P) N := NumIndets(); Q := P; Q := 2*Indet(N)-P; Return Q; End; Define LexOrder(F,G) Return LT(F) < LT(G); End; -- LexOrder Define PrimaryTest(I,P) N := NumIndets(); Prime := Ideal(P); Prime.GBasis := Prime.Gens; J := SortedBy(I.Gens,Function('LexOrder')); For M:=2 To Len(J) Do Initial := LM(J[M]); If Count(Log(Initial),0) = NumIndets()-1 Then N := N-1; E := Deg(Initial); Quot := DivAlg(J[M] - Initial, [Indet(N)^(E-1)]); Quot := Quot.Quotients; Quot := Quot[1]; T := LC(Initial)*E*Indet(N) + Quot; J[M] := E^E*LC(Initial)^(E-1)*J[M]; If NF(J[M]-T^E, Prime) = 0 Then Coeff := [Num(C_) | C_ In Coefficients(T)]; Content := GCD(Coeff); T := T/Content; If LC(T) < 0 Then T := -T; EndIf; Prime := Prime + Ideal(T); Prime.GBasis := Prime.Gens; Else Return Ideal(0); EndIf; EndIf; EndFor; Return Prime; End; -- PrimaryTest Define PrimDecGTZ(I) RingCounter := 0; Return ZeroPrimDecGTZ(I,RingCounter,0,0); End; -- PrimDecGTZ --Define PrimDecGTZ(I,Flag) -- Return ZeroPrimDecGTZ(I,0,Flag); --End; Define ZeroPrimDecGTZ(I, Var RingCounter, Flag, ChildProcess) Set FullRed; N := NumIndets(); PDRingName := "RingForPrimaryDecomposition_" + Sprint(RingCounter); Var(PDRingName) ::= CoeffRing[x[1..N]], Lex; Using Var(PDRingName) Do I_ := Image(I, RMap(x)); Result := []; Rest := []; L := Indets(); L[N] := LinMap(Flag); Phi := RMap(L); L[N] := LinInv(L[N]); PhiInv := RMap(L); J := Image(I, Phi); L[N] := Indet(N); If IsHomog(Gens(J)) = 0 Then HomPDR_ ::= CoeffRing[x[1..N],w], Lex; Using HomPDR_ Do HI := Image(J, RMap(x)); HI := Homogenized(w,HI); HI.PSeries := Poincare(HomPDR_/HI); ReducedGBasis(HI); EndUsing; J.GBasis := Image(HI.GBasis,RMap(Concat(L,[1]))); Clear HomPDR_; Destroy HomPDR_; Else J.GBasis := ReducedGBasis(J); EndIf; Foreach T In J.GBasis Do V := Log(LT(T)); Remove(V,N); If IsZero(Vector(V)) Then F := T; Break; EndIf; EndForeach; Fac := [K | K In Factor(F) And Deg(K[1]) > 0]; Foreach K In Fac Do P := K[1]^K[2]; Primary := Ideal(ReducedGBasis(J+Ideal(P))); Prime := PrimaryTest(Primary, K[1]); If IsZero(Prime) Then Append(Rest,I_ + Image(Ideal(P),PhiInv)); Else Append(Result,[Ideal(ReducedGBasis(I_ + Image(Ideal(P),PhiInv))), Ideal(ReducedGBasis(Image(Prime,PhiInv)))]); EndIf; EndForeach; Foreach K In Rest Do RingCounter := RingCounter + 1; If RingCounter = N Then Flag := 1; Elsif RingCounter = 2*N Then Flag := 2; EndIf; Result := Concat(Result, ZeroPrimDecGTZ(K, RingCounter, Flag, 1)); EndForeach; EndUsing; If ChildProcess = 0 Then PrintLn "Ring counter ", RingCounter; EndIf; PrimaryDecomposition := Image(Result, RMap(Indets())); Clear Var(PDRingName); Destroy Var(PDRingName); Return PrimaryDecomposition; End; -- ZeroPrimDecGTZ