1Package $DesignOfExperiments 2 3-- D := [0..3, 0..2]; 4-- D := [-1..1, -1..1]; 5-- FacCl := [1, x, y, x^2, y^2]; 6export FullDesignPoints; -- FullDesignPoints(D); 7export FullDesignLevels; -- FullDesignLevels(D); 8export FullDesignQB; -- FullDesignQB(R, D); 9export FullDesignPolys; -- FullDesignPolys(R, D); 10export corners; -- corners(FD_QB, SmallerQB); 11export border; -- border(FD_QB, SmallerQB); 12export GBRelations; -- GBRelations(D, SmallerQB) 13export BBRelations; -- BBRelations(D, SmallerQB) 14export AllGBases; -- AllGBases(GBRelations(D,SmallerQB)) 15export AllBBases; -- AllBBases(BBRelations(D,SmallerQB)) 16 17export FractionQB; -- FractionQB(FD_QB, [pp1, pp2,...]); 18 19 20 21------------------------ [ Procedures ] ---------------------------------- 22 23-- The number of parameters necessary for Corner, QB. 24 25Define MakeNumParams(Corners, QB) 26 Return sum([len([P in QB | P<T]) | T in Corners]); 27EndDefine; -- MakeNumParams 28 29 30-- For each term T_i in QB, gives the (parametric) polynomial G_i 31Define FillGB(Corners, QB) 32 P := RingOf(Corners[1]); 33 phi := CoeffEmbeddingHom(P); 34 params := apply(phi, indets(CoeffRing(P))); 35 Index := 0; 36 L := []; 37 For i := 1 To len(Corners) Do 38 K := [T In QB | T<Corners[i]]; 39 append(ref L, Corners[i] + sum([params[Index+j]*K[j] | j in 1..len(K)])); 40 Index := Index + len(K); 41 EndFor; 42 Return L; 43EndDefine; -- FillGB 44 45 46Define FillBB(Border, QB) 47 P := RingOf(Border[1]); 48 phi := CoeffEmbeddingHom(P); 49 params := apply(phi, indets(CoeffRing(P))); 50 Index := 0; 51 L := []; 52 For i := 1 To len(Border) Do 53 append(ref L, record[head:=Border[i], 54 tail:=sum([params[Index+j]*QB[j] | j in 1..len(QB)])]); 55 Index := Index + len(QB); 56 EndFor; 57 Return L; 58EndDefine; -- FillBB 59 60 61Define corners(BigQB, SmallQB) 62 return interreduced(diff(BigQB,SmallQB)); 63EndDefine; -- corners 64 65 66Define border(BigQB, SmallQB) 67 IsInBorder := 68 func(T) 69 foreach X in indets(RingOf(BigQB[1])) do 70 if IsDivisible(T,X) and (T/X isin SmallQB) then return true; endif; 71 endforeach; 72 return false; 73 endfunc; 74 return [ T in diff(BigQB, SmallQB) | IsInBorder(T) ]; 75EndDefine; -- corners 76 77 78define FullDesignPoints(FD) 79 return CartesianProductList(FD); 80enddefine; -- FullDesignPoints 81 82 83define FullDesignLevels(FD) 84 return [ len(D) | D in FD]; 85enddefine; -- FullDesignLevels 86 87 88Define FullDesignQB(P, FD) 89 if NumIndets(P) < len(FD) then error("Insufficient NumIndets in ring"); endif; 90 Return [MakeTerm(P, L) | L in CartesianProductList([0..len(D)-1 | D In FD])]; 91EndDefine; -- FullDesign 92 93 94Define FullDesignPolys(P, FD) 95 if NumIndets(P) < len(FD) then error("Insufficient NumIndets in ring"); endif; 96 Return [product([ indet(P,i)-a | a in FD[i]]) | i in 1..len(FD)]; 97EndDefine; -- FullDesignEq 98 99 100Define FractionQB(QB, Corners) 101 I := ideal(Corners); 102 Return [T in QB | not(T IsIn I)]; 103EndDefine; -- CornerFullDesign 104 105 106-------------[ Point Solving ] 107 108-- Gives the list of fractions described by the specialization of Sys for 109-- every element of Sol. 110Define Points(equations, GenericEQ, ParameterSols) 111 Return [$.Point(equations, GenericEQ, Sol) | Sol in ParameterSols]; 112EndDefine; 113 114 115Define Point(equations, GenericEQ, ParameterSol) 116 S := RingOf(GenericEQ[1]); 117 R := RingOf(equations[1]); 118 phiSR := PolyRingHom(S, R, PolyAlgebraHom(CoeffRing(S),RingQQ(),ParameterSol), indets(R)); 119 ActualSys := concat(equations, apply(phiSR, GenericEQ)); 120 Return RationalSolve(ActualSys); 121EndDefine; 122 123 124 125-----------------------[Buchberger procedures]----------------------------- 126 127 128------------[ Terms ]------------ 129 130Define TermsAreCoprime(S,T) 131 L1 := exponents(S); 132 L2 := exponents(T); 133 For I := 1 To NumIndets(RingOf(T)) Do 134 if L1[I]*L2[I] <> 0 then return false; endif; 135 EndFor; 136 Return true; 137EndDefine; 138 139-----------------[ Pairs ]--------------- 140 141Define BCriterion(ref Pairs, T) --> T is the LT of the new polynomial 142 Pairs := [P in Pairs | not(IsDivisible(P.lcm, T) 143 and lcm(T,LT(P.f1)) <> P.lcm 144 and lcm(T,LT(P.f2)) <> P.lcm) ]; 145EndDefine; 146 147 148Define NewPair(F,G) 149 return Record[f1:=F, f2:=G, lcm:=lcm(LT(F),LT(G)), 150 IsCoprime := $.TermsAreCoprime(LT(F),LT(G))]; 151EndDefine; 152 153 154Define InsertPairGM(ref Pairs, P) 155 LCM_P := P.lcm; 156 ToBeInserted := true; 157 for i := 1 to len(Pairs) do 158 LCM_i := Pairs[i].lcm; 159 if LCM_i = LCM_P Then 160 if P.IsCoprime then Pairs[i] := P; endif; 161 ToBeInserted := false; 162 elif IsDivisible(LCM_i, LCM_P) Then Pairs[i] := "remove"; 163 elif IsDivisible(LCM_P, LCM_i) Then ToBeInserted := false; 164 endif; 165 if not(ToBeInserted) then break; endif; 166 endfor; 167 If ToBeInserted Then append(ref Pairs, P); EndIf; 168 Pairs := [ PPP In Pairs | type(PPP) <> STRING ]; 169EndDefine; -- InsertPairGM 170 171 172-------------------[ Buchberger auxiliary ]----------------- 173 174Define SPoly(P) 175 F := P.f1; G := P.f2; lcmP := P.lcm; 176 phi := CoeffEmbeddingHom(RingOf(F)); 177 Return phi(LC(G))*(lcmP/LT(F))*F - phi(LC(F))*(lcmP/LT(G))*G; 178EndDefine; 179 180 181-------------------[ Buchberger Main ]----------------- 182 183Define BuildNewPairs(ref GB, F) 184 Pairs := []; 185 for i := 1 to len(GB) do 186 P := NewPair(GB[i], F); 187 if LT(GB[i]) = P.lcm then GB[i] := "remove"; endif; 188 $.InsertPairGM(ref Pairs, P); 189 endfor; 190 GB := [FF In GB | type(FF) <> STRING ]; 191 Return [ PP In Pairs | not(PP.IsCoprime) ]; 192EndDefine; 193 194 195Define UpdateBasisAndPairs(ref GB, ref Pairs, f) 196 $.BCriterion(ref Pairs, f); 197 Pairs_f := $.BuildNewPairs(ref GB, f); 198 append(ref GB, f); 199 Pairs := concat(Pairs, Pairs_f); 200EndDefine; 201 202 203Define relations_Buchberger(GB, GenericPolys) 204 Pairs := []; 205 Foreach f In GenericPolys Do 206 $.UpdateBasisAndPairs(ref GB, ref Pairs, f); 207 EndForeach; 208 PrintLn "len(Pairs) = ", len(Pairs); 209 PrintLn "GB = ", GB; 210 Out := [NR(SPoly(PP), GB) | PP in Pairs]; 211 return record[GB := GB, 212 rel := ConcatLists([coefficients(f) | f In Out and not(IsZero(f))])]; 213EndDefine; -- relations_Buchberger 214 215 216Define GBRelations(FullDesign, FractionQB) 217 rings := diff(RingsOf(FractionQB), CurrentTypes()); 218 if len(rings)>1 then error("more than 1 ring in FractionQB");endif; 219 P := rings[1]; 220 FD_QB := FullDesignQB(P, FullDesign); 221 PrintLn "FD_QB = ", FD_QB; 222 foreach T in FractionQB do 223 if not(T IsIn FD_QB) then error(sprint(T)+" not contained in QB"); endif; 224 endforeach; 225 FractionCorners := corners(FD_QB, FractionQB); 226 PrintLn "FractionCorners = ", FractionCorners; 227 S := BuildParamRing(P, "a", MakeNumParams(FractionCorners, FractionQB)); 228 phiPS := PolyRingHom(P,S, CanonicalHom(CoeffRing(P),CoeffRing(S)), indets(S)); 229 GenericEqns := FillGB(apply(phiPS, FractionCorners), apply(phiPS, FractionQB)); 230 PrintLn "GenericEqns = "; indent(GenericEqns); 231 FD_polys := FullDesignPolys(S, FullDesign); 232 return relations_Buchberger(FD_polys, GenericEqns); 233EndDefine; -- GBRelations 234 235 236define FindReducerBB(f, L) 237 for j:=1 to len(L) do 238 if IsDivisible(LT(f),L[j].head) then return j; endif; 239 endfor; 240 return 0; 241enddefine; -- FindReducerBB 242 243define NRBB(f, L) 244// PrintLn "NRBB(f, L) ", f; 245 TmpRes := zero(RingOf(f)); 246 if IsZero(f) then return TmpRes; endif; 247 i := FindReducerBB(f, L); 248 repeat 249 while i = 0 do 250 TmpRes := TmpRes + LM(f); f := f - LM(f); 251 if IsZero(f) then return TmpRes; endif; 252 i := FindReducerBB(f, L); 253 endwhile; 254 f := f - (LM(f)/L[i].head)*(L[i].head + L[i].tail); -- NR 255// PrintLn "f --> ", f; 256 if IsZero(f) then return TmpRes; endif; 257 i := FindReducerBB(f, L); 258 endrepeat; 259enddefine; -- NRBB 260 261 262Define MultiplicationMat(H, x, QB) 263// PrintLn "QB = ", QB; 264 SortBy(ref H, func(a,b) return a.head>b.head; endfunc); 265 return transposed(mat([coefficients(NRBB(x*T, H), QB) | T in QB])); 266EndDefine; -- MultiplicationMat 267 268 269define eval(f,M) 270 v := ZeroMat(RingOf(M[1]), NumRows(M[1]), NumRows(M[1])); 271 foreach m in monomials(f) do 272 L := log(LT(m)); 273 v := v+ RingElem(RingOf(v),LC(m)) * product([ M[i]^L[i] | i in 1..len(M)]); 274// PrintLn "m = ", m; 275// PrintLn "v = ", v; 276 endforeach; 277// PrintLn "f = ", f; 278 return v; 279enddefine; -- eval 280 281 282define relations_MultMat(FD_polys, GenericEqnsBB, FractionQB, FractionBorder) 283 H := concat([ record[head:=LT(f), tail:=f-LT(f)] | f in FD_polys], GenericEqnsBB); 284 X := [indet(RingOf(f), UnivariateIndetIndex(f)) | f in FD_polys]; 285 M := [MultiplicationMat(H, x, FractionQB) | x in X]; 286// PrintLn "M = ", M; 287 commutators := [M[i[1]]*M[i[2]]-M[i[2]]*M[i[1]] | i in subsets(1..len(M),2)]; 288// PrintLn "commutators = ", commutators; 289 ColMatrices := [ GetCol(eval(f,M),1) | f in FD_polys and IsZero(NR(LT(f), FractionBorder)) ]; 290// PrintLn "ColMatrices = ", ColMatrices; 291 AllMatEntries := flatten(concat([GetRows(N) | N in commutators],ColMatrices)); 292 return record[BB := H, 293 rel := [ g in AllMatEntries | not(IsZero(g))]]; 294enddefine; -- relations_MultMat 295 296 297Define BBRelations(FullDesign, FractionQB) 298 rings := diff(RingsOf(FractionQB), CurrentTypes()); 299 if len(rings)>1 then error("more than 1 ring in FractionQB");endif; 300 P := rings[1]; 301 FD_QB := FullDesignQB(P, FullDesign); 302// PrintLn "FD_QB = ", FD_QB; 303 foreach T in FractionQB do 304 if not(T IsIn FD_QB) then error(sprint(T)+" not contained in QB"); endif; 305 endforeach; 306 FractionBorder := border(FD_QB, FractionQB); -- C2 307// PrintLn "FractionBorder = ", FractionBorder; 308 S := BuildParamRing(P, "a", len(FractionBorder)*len(FractionQB)); 309 phiPS := PolyRingHom(P,S, CanonicalHom(CoeffRing(P),CoeffRing(S)), indets(S)); 310 GenericEqnsBB := FillBB(apply(phiPS, FractionBorder), apply(phiPS, FractionQB)); 311// PrintLn "GenericEqnsBB = "; indent(GenericEqnsBB); 312 FD_polys := FullDesignPolys(S, FullDesign); 313// PrintLn "FD_polys = "; indent(FD_polys); 314 return relations_MultMat(FD_polys, GenericEqnsBB, 315 apply(phiPS,FractionQB), apply(phiPS,FractionBorder)); 316EndDefine; 317 318 319define AllGBases(Rel) 320 TopLevel CurrentRing; 321 P := CurrentRing; 322 ParameterSol := RationalSolve(Rel.rel); 323 S := RingOf(Rel.GB[1]); 324 CoeffHoms := [PolyAlgebraHom(CoeffRing(S),CoeffRing(P), 325 [AsRAT(x)|x in sol]) | sol in ParameterSol]; 326 return [ apply(PolyRingHom(S,P,phi,indets(P)), Rel.GB) | phi in CoeffHoms]; 327enddefine; -- AllGBases 328 329 330define AllBBases(Rel) 331 TopLevel CurrentRing; 332 P := CurrentRing; 333 ParameterSol := RationalSolve(Rel.rel); 334 BB := [ f.head + f.tail | f in Rel.BB]; 335 S := RingOf(BB[1]); 336 CoeffHoms := [PolyAlgebraHom(CoeffRing(S),CoeffRing(P), 337 [AsRAT(x)|x in sol]) | sol in ParameterSol]; 338 return [ apply(PolyRingHom(S,P,phi,indets(P)), BB) | phi in CoeffHoms]; 339enddefine; -- AllBBases 340 341 342----------------------- [ Misc Procedures ] -------------- 343 344Define Delta_Var(D_Var, L) 345 L:=concat([0],L); 346 S:=1; 347 For I:=2 To len(L) Do 348 S:=S*binomial(AsINT(D_Var-L[I-1]), AsINT(L[I]-L[I-1])); 349 EndFor; 350 Return S; 351EndDefine; 352 353 354-- L is a list of terms, D are the vars max degs 355Define Delta(L,D) 356 L:=GetCols(mat([exponents(T)|T In L]) ); 357 L:=[sorted(R)|R In L]; 358 Ris:=1; 359 For I:=1 To len(L) Do 360 Ris := Ris*$.Delta_Var(D[I], L[I]); 361 EndFor; 362 Return Ris; 363EndDefine; 364 365--------------------------- [ Imaging Procedures ] ------------------------- 366 367Define BuildParamRing(P, ParamsName, NumParams) 368// PrintLn NumParams; 369 R := NewPolyRing(RingQQ(), SymbolRange(ParamsName, 1, NumParams)); 370// K := NewFractionField(R); 371 return NewPolyRing(R, IndetSymbols(P)); 372EndDefine; 373 374EndPackage; -- $DesignOfExperiments 375 376 377-- PrintLn "/**/ FD := [0..1, 0..1, 0..1]; QB := [1, x, y, z];"; 378-- PrintLn "/**/ FD := [-1..1, -1..1]; QB := [1, x, y, x^2, y^2]; // ex (h)"; 379-- PrintLn "/**/ FD := [-1..1, [-1,1]]; QB := [1, x, y]; // ex (g)"; 380-- PrintLn; 381-- PrintLn "/**/ Rel := GBRelations(FD, QB);"; 382-- PrintLn "/**/ multiplicity(RingOf(Rel.rel[1])/ideal(Rel.rel));"; 383-- PrintLn "/**/ AllGBs := AllGBases(Rel); indent(AllGBs);"; 384-- PrintLn "/**/ indent([RationalSolve(B) | B in AllGBs]);"; 385-- PrintLn; 386-- PrintLn "/**/ Rel := BBRelations(FD, QB);"; 387-- PrintLn "/**/ multiplicity(RingOf(Rel.rel[1])/ideal(Rel.rel));"; 388-- PrintLn "/**/ AllBBs := AllBBases(Rel); indent(AllBBs);"; 389-- PrintLn "/**/ indent([RationalSolve(B) | B in AllBBs]);"; 390-- PrintLn "/**/ indent(diff(subsets(FullDesignPoints(FD),3), [RationalSolve(B) | B in AllBBs]));"; 391