[35f23d] | 1 | /////////////////////////////////////////////////////////////////////////////// |
---|
[341696] | 2 | version="$Id$"; |
---|
[fd3fb7] | 3 | category="Singularities"; |
---|
[35f23d] | 4 | info=" |
---|
[8942a5] | 5 | LIBRARY: qhmoduli.lib Moduli Spaces of Semi-Quasihomogeneous Singularities |
---|
[35f23d] | 6 | AUTHOR: Thomas Bayer, email: bayert@in.tum.de |
---|
| 7 | |
---|
| 8 | PROCEDURES: |
---|
| 9 | ArnoldAction(f, [G, w]) Induced action of G_f on T_. |
---|
| 10 | ModEqn(f) Equations of the moduli space for principal part f |
---|
| 11 | QuotientEquations(G,A,I) Equations of Variety(I)/G w.r.t. action 'A' |
---|
| 12 | StabEqn(f) Equations of the stabilizer of f. |
---|
| 13 | StabEqnId(I, w) Equations of the stabilizer of the qhom. ideal I. |
---|
| 14 | StabOrder(f) Order of the stabilizer of f. |
---|
| 15 | UpperMonomials(f, [w]) Upper basis of the Milnor algebra of f. |
---|
| 16 | |
---|
| 17 | Max(data) maximal integer contained in 'data' |
---|
| 18 | Min(data) minimal integer contained in 'data' |
---|
| 19 | "; |
---|
| 20 | |
---|
| 21 | // NOTE: This library has been written in the frame of the diploma thesis |
---|
| 22 | // 'Computing moduli spaces of semiquasihomogeneous singularities and an |
---|
| 23 | // implementation in Singular', Arbeitsgruppe Algebraische Geometrie, |
---|
| 24 | // Fachbereich Mathematik, University Kaiserslautern, |
---|
| 25 | // Advisor: Prof. Gert-Martin Greuel |
---|
| 26 | |
---|
| 27 | LIB "rinvar.lib"; |
---|
| 28 | |
---|
| 29 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 30 | |
---|
| 31 | proc ModEqn(poly f, list #) |
---|
| 32 | "USAGE: ModEqn(f [, opt]); poly f; int opt; |
---|
[1f92589] | 33 | PURPOSE: compute equations of the moduli space of semiquasihomogenos hypersurface singularity with principal part f w.r.t. right equivalence |
---|
| 34 | ASSUME: f quasihomogeneous polynomial with an isolated singularity at 0 |
---|
[f04aaf] | 35 | RETURN: polynomial ring, possibly a simple extension of the ground field of |
---|
[35f23d] | 36 | the basering, containing the ideal 'modid' |
---|
| 37 | - 'modid' is the ideal of the moduli space if opt is even (> 0). |
---|
| 38 | otherwise it contains generators of the coordinate ring R of the |
---|
| 39 | moduli space (note : Spec(R) is the moduli space) |
---|
| 40 | OPTIONS: 1 compute equations of the mod. space, |
---|
[fd5013] | 41 | 2 use a primary decomposition, |
---|
| 42 | 4 compute E_f0, i.e., the image of G_f0, |
---|
| 43 | to combine options, add their value, default: opt =7 |
---|
[35f23d] | 44 | EXAMPLE: example ModEqn; shows an example |
---|
| 45 | " |
---|
| 46 | { |
---|
| 47 | int sizeOfAction, i, dimT, nonLinearQ, milnorNr, dbPrt; |
---|
| 48 | int imageQ, opt; |
---|
| 49 | intvec wt; |
---|
| 50 | ideal B; |
---|
| 51 | list Gf, tIndex, sList; |
---|
| 52 | string ringSTR; |
---|
| 53 | |
---|
| 54 | dbPrt = printlevel-voice+2; |
---|
| 55 | if(size(#) > 0) { opt = #[1]; } |
---|
| 56 | else { opt = 7; } |
---|
[1f9a84] | 57 | if(opt div 4 > 0) { imageQ = 1; opt = opt - 4;} |
---|
[35f23d] | 58 | else { imageQ = 0; } |
---|
| 59 | |
---|
| 60 | wt = weight(f); |
---|
| 61 | milnorNr = vdim(std(jacob(f))); |
---|
| 62 | if(milnorNr == -1) { |
---|
| 63 | ERROR("the polynomial " + string(f) + " has a nonisolated singularity at 0"); |
---|
| 64 | } // singularity not isolated |
---|
| 65 | |
---|
| 66 | // 1st step : compute a basis of T_ |
---|
| 67 | |
---|
| 68 | B = UpperMonomials(f, wt); |
---|
| 69 | dimT = size(B); |
---|
| 70 | dbprint(dbPrt, "moduli equations of f = " + string(f) + ", f has Milnor number = " + string(milnorNr)); |
---|
| 71 | dbprint(dbPrt, " upper basis = " + string(B)); |
---|
| 72 | if(size(B) > 1) { |
---|
| 73 | |
---|
| 74 | // 2nd step : compute the stabilizer G_f of f |
---|
| 75 | |
---|
| 76 | dbprint(dbPrt, " compute equations of the stabilizer of f, called G_f"); |
---|
| 77 | Gf = StabEqn(f); |
---|
| 78 | dbprint(dbPrt, " order of the stabilizer = " + string(StabOrder(Gf))); |
---|
| 79 | |
---|
| 80 | // 3rd step : compute the induced action of G_f on T_ by means of a theorem of Arnold |
---|
| 81 | |
---|
| 82 | dbprint(dbPrt, " compute the induced action"); |
---|
| 83 | def RME1 = ArnoldAction(f, Gf, B); |
---|
| 84 | setring(RME1); |
---|
| 85 | export(RME1); |
---|
| 86 | dbprint(dbPrt, " G_f = " + string(stabid)); |
---|
| 87 | dbprint(dbPrt, " action of G_f : " + string(actionid)); |
---|
| 88 | |
---|
| 89 | // 4th step : linearize the action of G_f |
---|
| 90 | |
---|
| 91 | sizeOfAction = size(actionid); |
---|
| 92 | def RME2 = LinearizeAction(stabid, actionid, nvars(Gf[1])); |
---|
| 93 | setring RME2; |
---|
| 94 | export(RME2); |
---|
[3b77465] | 95 | kill RME1; |
---|
[35f23d] | 96 | |
---|
| 97 | if(size(actionid) == sizeOfAction) { nonLinearQ = 0;} |
---|
| 98 | else { |
---|
| 99 | nonLinearQ = 1; |
---|
| 100 | dbprint(dbPrt, " linearized action = " + string(actionid)); |
---|
| 101 | dbprint(dbPrt, " embedding of T_ = " + string(embedid)); |
---|
| 102 | } |
---|
| 103 | |
---|
| 104 | |
---|
| 105 | |
---|
| 106 | if(!imageQ) { // do not compute the image of Gf |
---|
| 107 | // 5th step : set E_f = G_f, |
---|
| 108 | dbprint(dbPrt, " compute equations of the quotient T_/G_f"); |
---|
| 109 | def RME3 = basering; |
---|
| 110 | } |
---|
| 111 | else { |
---|
| 112 | |
---|
| 113 | // 5th step : compute the ideal and the action of E_f |
---|
| 114 | |
---|
| 115 | dbprint(dbPrt, " compute E_f"); |
---|
| 116 | def RME3 = ImageGroup(groupid, actionid); |
---|
| 117 | setring(RME3); |
---|
| 118 | ideal embedid = imap(RME2, embedid); |
---|
| 119 | dbprint(dbPrt, " E_f = (" + string(groupid) + ")"); |
---|
| 120 | dbprint(dbPrt, " action of E'f = " + string(actionid)); |
---|
| 121 | dbprint(dbPrt, " compute equations of the quotient T_/E_f"); |
---|
| 122 | } |
---|
| 123 | export(RME3); |
---|
[3b77465] | 124 | kill RME2; |
---|
[35f23d] | 125 | |
---|
| 126 | // 6th step : compute the equations of the quotient T_/E_f |
---|
| 127 | |
---|
| 128 | ideal G = groupid; ideal variety = embedid; |
---|
[3b77465] | 129 | kill groupid,embedid; |
---|
[35f23d] | 130 | def RME4 = QuotientEquations(G, actionid, variety, opt); |
---|
| 131 | setring RME4; |
---|
[393c47] | 132 | string @mPoly = string(minpoly); |
---|
[3b77465] | 133 | kill RME3; |
---|
[35f23d] | 134 | export(RME4); |
---|
| 135 | |
---|
| 136 | // simplify the ideal and create a new ring with propably less variables |
---|
| 137 | |
---|
| 138 | if(opt == 1 || opt == 3) { // equations computed ? |
---|
| 139 | sList = SimplifyIdeal(id, 0, "Y"); |
---|
| 140 | ideal newid = sList[1]; |
---|
| 141 | dbprint(dbPrt, " number of equations = " + string(size(sList[1]))); |
---|
| 142 | dbprint(dbPrt, " number of variables = " + string(size(sList[3]))); |
---|
| 143 | ringSTR = "ring RME5 = (" + charstr(basering) + "), (Y(1.." + string(size(sList[3])) + ")),dp;"; |
---|
| 144 | execute(ringSTR); |
---|
[393c47] | 145 | execute("minpoly = number(" + @mPoly + ");"); |
---|
[35f23d] | 146 | ideal modid = imap(RME4, newid); |
---|
| 147 | } |
---|
| 148 | else { |
---|
| 149 | def RME5 = RME4; |
---|
| 150 | setring(RME5); |
---|
| 151 | ideal modid = imap(RME4, id); |
---|
| 152 | } |
---|
| 153 | export(modid); |
---|
[3b77465] | 154 | kill RME4; |
---|
[35f23d] | 155 | } |
---|
| 156 | else { |
---|
| 157 | def RME5 = basering; |
---|
| 158 | ideal modid = maxideal(1); |
---|
| 159 | if(size(B) == 1) { // 1-dimensional |
---|
| 160 | modid[size(modid)] = 0; |
---|
| 161 | modid = simplify(modid,2); |
---|
| 162 | } |
---|
| 163 | export(modid); |
---|
| 164 | } |
---|
| 165 | dbprint(dbPrt, " |
---|
| 166 | // 'ModEqn' created a new ring. |
---|
| 167 | // To see the ring, type (if the name of the ring is R): |
---|
| 168 | show(R); |
---|
[1f92589] | 169 | // To access the ideal of the moduli space of semiquasihomogeneous singularities |
---|
[35f23d] | 170 | // with principal part f, type |
---|
| 171 | def R = ModEqn(f); setring R; modid; |
---|
| 172 | // 'modid' is the ideal of the moduli space. |
---|
| 173 | // if 'opt' = 0 or even, then 'modid' contains algebra generators of S s.t. |
---|
| 174 | // spec(S) = moduli space of f. |
---|
| 175 | "); |
---|
| 176 | return(RME5); |
---|
| 177 | } |
---|
| 178 | example |
---|
| 179 | {"EXAMPLE:"; echo = 2; |
---|
| 180 | ring B = 0,(x,y), ls; |
---|
| 181 | poly f = -x4 + xy5; |
---|
| 182 | def R = ModEqn(f); |
---|
| 183 | setring R; |
---|
| 184 | modid; |
---|
| 185 | } |
---|
| 186 | |
---|
| 187 | |
---|
| 188 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 189 | |
---|
| 190 | proc QuotientEquations(ideal G, ideal Gaction, ideal embedding, list#) |
---|
| 191 | "USAGE: QuotientEquations(G,action,emb [, opt]); ideal G,action,emb;int opt |
---|
| 192 | PURPOSE: compute the quotient of the variety given by the parameterization |
---|
| 193 | 'emb' by the linear action 'action' of the algebraic group G. |
---|
| 194 | ASSUME: 'action' is linear, G must be finite if the Reynolds operator is |
---|
| 195 | needed (i.e., NullCone(G,action) returns some non-invariant polys) |
---|
[f04aaf] | 196 | RETURN: polynomial ring over a simple extension of the ground field of the |
---|
[35f23d] | 197 | basering, containing the ideals 'id' and 'embedid'. |
---|
[1f92589] | 198 | - 'id' contains the equations of the quotient, if opt = 1; |
---|
| 199 | if opt = 0, 2, 'id' contains generators of the coordinate ring R |
---|
[8942a5] | 200 | of the quotient (Spec(R) is the quotient) |
---|
[1f92589] | 201 | - 'embedid' = 0, if opt = 1; |
---|
| 202 | if opt = 0, 2, it is the ideal defining the equivariant embedding |
---|
[35f23d] | 203 | OPTIONS: 1 compute equations of the quotient, |
---|
[fd5013] | 204 | 2 use a primary decomposition when computing the Reynolds operator,@* |
---|
| 205 | to combine options, add their value, default: opt =3. |
---|
[35f23d] | 206 | EXAMPLE: example QuotientEquations; shows an example |
---|
| 207 | " |
---|
| 208 | { |
---|
| 209 | int i, opt, primaryDec, relationsQ, dbPrt; |
---|
| 210 | ideal Gf, variety; |
---|
| 211 | intvec wt; |
---|
| 212 | |
---|
| 213 | dbPrt = printlevel-voice+3; |
---|
| 214 | if(size(#) > 0) { opt = #[1]; } |
---|
| 215 | else { opt = 3; } |
---|
| 216 | |
---|
[1f9a84] | 217 | if(opt div 2 > 0) { primaryDec = 1; opt = opt - 2; } |
---|
[35f23d] | 218 | else { primaryDec = 0; } |
---|
| 219 | if(opt > 0) { relationsQ = 1;} |
---|
| 220 | else { relationsQ = 0; } |
---|
| 221 | |
---|
| 222 | Gf = std(G); |
---|
| 223 | variety = EquationsOfEmbedding(embedding, nvars(basering) - size(Gaction)); |
---|
| 224 | |
---|
| 225 | if(size(variety) == 0) { // use Hilbert function ! |
---|
[62dc18e] | 226 | //for(i = 1; i <= ncols(Gaction); i ++) { wt[i] = 1;} |
---|
| 227 | for(i = 1; i <= nvars(basering); i ++) { wt[i] = 1;} |
---|
[35f23d] | 228 | } |
---|
| 229 | def RQER = InvariantRing(Gf, Gaction, primaryDec); // compute the nullcone of the linear action |
---|
| 230 | |
---|
| 231 | def RQEB = basering; |
---|
| 232 | setring(RQER); |
---|
| 233 | export(RQER); |
---|
| 234 | |
---|
| 235 | if(relationsQ > 0) { |
---|
| 236 | dbprint(dbPrt, " compute equations of the variety (" + string(size(imap(RQER, invars))) + " invariants) "); |
---|
| 237 | if(!defined(variety)) { ideal variety = imap(RQEB, variety); } |
---|
| 238 | if(wt[1] > 0) { |
---|
| 239 | def RQES = ImageVariety(variety, imap(RQER, invars), wt); |
---|
| 240 | } |
---|
| 241 | else { |
---|
| 242 | def RQES = ImageVariety(variety, imap(RQER, invars)); // forget imap |
---|
| 243 | } |
---|
| 244 | setring(RQES); |
---|
| 245 | ideal id = imageid; |
---|
| 246 | ideal embedid = 0; |
---|
| 247 | } |
---|
| 248 | else { |
---|
| 249 | def RQES = basering; |
---|
| 250 | ideal id = imap(RQER, invars); |
---|
| 251 | ideal embedid = imap(RQEB, variety); |
---|
| 252 | } |
---|
[3b77465] | 253 | kill RQER; |
---|
[35f23d] | 254 | export(id); |
---|
| 255 | export(embedid); |
---|
| 256 | return(RQES); |
---|
| 257 | } |
---|
| 258 | |
---|
| 259 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 260 | |
---|
| 261 | proc UpperMonomials(poly f, list #) |
---|
| 262 | "USAGE: UpperMonomials(poly f, [intvec w]) |
---|
| 263 | PURPOSE: compute the upper monomials of the milnor algebra of f. |
---|
[1f92589] | 264 | ASSUME: f is quasihomogeneous (w.r.t. w) |
---|
[35f23d] | 265 | RETURN: ideal |
---|
| 266 | EXAMPLE: example UpperMonomials; shows an example |
---|
| 267 | " |
---|
| 268 | { |
---|
| 269 | int i,d; |
---|
| 270 | intvec wt; |
---|
| 271 | ideal I, J; |
---|
| 272 | |
---|
| 273 | if(size(#) == 0) { wt = weight(f);} |
---|
| 274 | else { wt = #[1];} |
---|
| 275 | J = kbase(std(jacob(f))); |
---|
| 276 | d = deg(f, wt); |
---|
| 277 | for(i = 1; i <= size(J); i++) { if(deg(J[i], wt) > d) {I = I, J[i];} } |
---|
| 278 | return(simplify(I, 2)); |
---|
| 279 | } |
---|
| 280 | example |
---|
| 281 | {"EXAMPLE:"; echo = 2; |
---|
| 282 | ring B = 0,(x,y,z), ls; |
---|
| 283 | poly f = -z5+y5+x2z+x2y; |
---|
| 284 | UpperMonomials(f); |
---|
| 285 | } |
---|
| 286 | |
---|
| 287 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 288 | |
---|
| 289 | proc ArnoldAction(poly f, list #) |
---|
| 290 | "USAGE: ArnoldAction(f, [Gf, B]); poly f; list Gf, B; |
---|
| 291 | 'Gf' is a list of two rings (coming from 'StabEqn') |
---|
[1f92589] | 292 | PURPOSE: compute the induced action of the stabilizer G of f on T_, where |
---|
[35f23d] | 293 | T_ is given by the upper monomials B of the Milnor algebra of f. |
---|
[1f92589] | 294 | ASSUME: f is quasihomogeneous |
---|
[f04aaf] | 295 | RETURN: polynomial ring over the same ground field, containing the ideals |
---|
[35f23d] | 296 | 'actionid' and 'stabid'. |
---|
[1f92589] | 297 | - 'actionid' is the ideal defining the induced action of Gf on T_ @* |
---|
[35f23d] | 298 | - 'stabid' is the ideal of the stabilizer Gf in the new ring |
---|
| 299 | EXAMPLE: example ArnoldAction; shows an example |
---|
| 300 | " |
---|
| 301 | { |
---|
| 302 | int i, offset, ub, pos, nrStabVars, dbPrt; |
---|
| 303 | intvec wt = weight(f); |
---|
| 304 | ideal B; |
---|
| 305 | list Gf, parts, baseDeg; |
---|
| 306 | string ringSTR1, ringSTR2, parName, namesSTR, varSTR; |
---|
| 307 | |
---|
| 308 | dbPrt = printlevel-voice+2; |
---|
| 309 | if(size(#) == 0) { |
---|
| 310 | Gf = StabEqn(f); |
---|
| 311 | B = UpperMonomials(f, wt); |
---|
| 312 | } |
---|
| 313 | else { |
---|
| 314 | Gf = #[1]; |
---|
| 315 | if(size(#) > 1) { B = #[2];} |
---|
| 316 | else {B = UpperMonomials(f, wt);} |
---|
| 317 | } |
---|
| 318 | if(size(B) == 0) { ERROR("the principal part " + string(f) + " has no upper monomials");} |
---|
| 319 | for(i = 1; i <= size(B); i = i + 1) { |
---|
| 320 | baseDeg[i] = deg(B[i], wt); |
---|
| 321 | } |
---|
| 322 | ub = Max(baseDeg) + 1; // max degree of an upper mono. |
---|
| 323 | def RAAB = basering; |
---|
| 324 | def STR1 = Gf[1]; |
---|
| 325 | def STR2 = Gf[2]; |
---|
| 326 | nrStabVars = nvars(STR1); |
---|
| 327 | |
---|
| 328 | dbprint(dbPrt, "ArnoldAction of f = ", f, ", upper base = " + string(B)); |
---|
| 329 | |
---|
| 330 | setring STR1; |
---|
[393c47] | 331 | string @mPoly = string(minpoly); |
---|
[35f23d] | 332 | setring RAAB; |
---|
| 333 | |
---|
| 334 | // setup new ring with s(..) and t(..) as parameters |
---|
| 335 | |
---|
| 336 | varSTR = string(maxideal(1)); |
---|
| 337 | ringSTR2 = "ring RAAS = "; |
---|
| 338 | if(npars(basering) == 1) { |
---|
| 339 | parName = parstr(basering); |
---|
| 340 | ringSTR2 = ringSTR2 + "(0, " + parstr(1) + "), "; |
---|
| 341 | } |
---|
| 342 | else { |
---|
| 343 | parName = "a"; |
---|
| 344 | ringSTR2 = ringSTR2 + "0, "; |
---|
| 345 | } |
---|
| 346 | offset = 1 + nrStabVars; |
---|
| 347 | namesSTR = "s(1.." + string(nrStabVars) + "), t(1.." + string(size(B)) + ")"; |
---|
| 348 | ringSTR2 = ringSTR2 + "(" + namesSTR + "), lp;"; |
---|
| 349 | ringSTR1 = "ring RAAR = (0, " + parName + "," + namesSTR + "), (" + varSTR + "), ls;"; // lp ? |
---|
| 350 | |
---|
| 351 | execute(ringSTR1); |
---|
| 352 | export(RAAR); |
---|
[393c47] | 353 | ideal upperBasis, stabaction, action, reduceIdeal; |
---|
[35f23d] | 354 | poly f, F, monos, h; |
---|
| 355 | |
---|
[393c47] | 356 | execute("reduceIdeal = " + @mPoly + ";"); reduceIdeal = reduceIdeal, imap(STR1, stabid); |
---|
[35f23d] | 357 | f = imap(RAAB, f); |
---|
| 358 | F = f; |
---|
| 359 | upperBasis = imap(RAAB, B); |
---|
| 360 | for(i = 1; i <= size(upperBasis); i = i + 1) { |
---|
| 361 | F = F + par(i + offset)*upperBasis[i]; |
---|
| 362 | } |
---|
| 363 | monos = F - f; |
---|
| 364 | stabaction = imap(STR2, actionid); |
---|
| 365 | |
---|
| 366 | // action of the stabilizer on the semiuniversal unfolding of f |
---|
| 367 | |
---|
| 368 | F = f + APSubstitution(monos, stabaction, reduceIdeal, wt, ub, nrStabVars, size(upperBasis)); |
---|
| 369 | |
---|
| 370 | // apply the theorem of Arnold |
---|
| 371 | |
---|
| 372 | h = ArnoldFormMain(f, upperBasis, F, reduceIdeal, nrStabVars, size(upperBasis)) - f; |
---|
| 373 | |
---|
| 374 | // extract the polynomials of the action of the stabilizer on T_ |
---|
| 375 | |
---|
| 376 | parts = MonosAndTerms(h, wt, ub); |
---|
| 377 | for(i = 1; i <= size(parts[1]); i = i + 1) { |
---|
| 378 | pos = FirstEntryQHM(upperBasis, parts[1][i]); |
---|
| 379 | action[pos] = parts[2][i]/parts[1][i]; |
---|
| 380 | } |
---|
| 381 | execute(ringSTR2); |
---|
[393c47] | 382 | execute("minpoly = number(" + @mPoly + ");"); |
---|
[35f23d] | 383 | ideal actionid = imap(RAAR, action); |
---|
| 384 | ideal stabid = imap(STR1, stabid); |
---|
| 385 | export(actionid); |
---|
| 386 | export(stabid); |
---|
[3b77465] | 387 | kill RAAR; |
---|
[35f23d] | 388 | dbprint(dbPrt, " |
---|
| 389 | // 'ArnoldAction' created a new ring. |
---|
| 390 | // To see the ring, type (if the name of the ring is R): |
---|
| 391 | show(R); |
---|
| 392 | // To access the ideal of the stabilizer G of f and its group action, |
---|
[1f92589] | 393 | // where f is the quasihomogeneous principal part, type |
---|
[35f23d] | 394 | def R = ArnoldAction(f); setring R; stabid; actionid; |
---|
| 395 | // 'stabid' is the ideal of the group G and 'actionid' is the ideal defining |
---|
| 396 | // the group action of the group G on T_. Note: this action might be nonlinear |
---|
| 397 | "); |
---|
| 398 | return(RAAS); |
---|
| 399 | } |
---|
| 400 | example |
---|
| 401 | {"EXAMPLE:"; echo = 2; |
---|
| 402 | ring B = 0,(x,y,z), ls; |
---|
| 403 | poly f = -z5+y5+x2z+x2y; |
---|
| 404 | def R = ArnoldAction(f); |
---|
| 405 | setring R; |
---|
| 406 | actionid; |
---|
| 407 | stabid; |
---|
| 408 | } |
---|
| 409 | |
---|
| 410 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 411 | |
---|
| 412 | proc StabOrder(list #) |
---|
[fd5013] | 413 | "USAGE: StabOrder(f); poly f |
---|
[35f23d] | 414 | PURPOSE: compute the order of the stabilizer group of f. |
---|
[1f92589] | 415 | ASSUME: f quasihomogeneous polynomial with an isolated singularity at 0 |
---|
[35f23d] | 416 | RETURN: int |
---|
| 417 | GLOBAL: varSubsList |
---|
| 418 | " |
---|
| 419 | { |
---|
| 420 | list stab; |
---|
| 421 | |
---|
| 422 | if(size(#) == 1) { stab = StabEqn(#[1]); } |
---|
| 423 | else { stab = #;} |
---|
| 424 | |
---|
| 425 | def RSTO = stab[1]; |
---|
| 426 | setring(RSTO); |
---|
| 427 | return(vdim(std(stabid))); |
---|
| 428 | } |
---|
| 429 | |
---|
| 430 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 431 | |
---|
| 432 | proc StabEqn(poly f) |
---|
| 433 | "USAGE: StabEqn(f); f polynomial |
---|
| 434 | PURPOSE: compute the equations of the isometry group of f. |
---|
[1f92589] | 435 | ASSUME: f semiquasihomogeneous polynomial with an isolated singularity at 0 |
---|
[fd5013] | 436 | RETURN: list of two rings 'S1', 'S2' |
---|
[1f92589] | 437 | - 'S1' contians the equations of the stabilizer (ideal 'stabid') @* |
---|
[35f23d] | 438 | - 'S2' contains the action of the stabilizer (ideal 'actionid') |
---|
| 439 | EXAMPLE: example StabEqn; shows an example |
---|
| 440 | GLOBAL: varSubsList, contains the index j s.t. x(i) -> x(i)t(j) ... |
---|
| 441 | " |
---|
| 442 | { |
---|
| 443 | dbprint(dbPrt, " |
---|
| 444 | // 'StabEqn' created a list of 2 rings. |
---|
| 445 | // To see the rings, type (if the name of your list is stab): |
---|
| 446 | show(stab); |
---|
| 447 | // To access the 1-st ring and map (and similair for the others), type: |
---|
| 448 | def S1 = stab[1]; setring S1; stabid; |
---|
| 449 | // S1/stabid is the coordinate ring of the variety of the |
---|
| 450 | // stabilizer, say G. If G x K^n --> K^n is the action of G on |
---|
| 451 | // K^n, then the ideal 'actionid' in the second ring describes |
---|
| 452 | // the dual map on the ring level. |
---|
| 453 | // To access the 2-nd ring and map (and similair for the others), type: |
---|
| 454 | def S2 = stab[2]; setring S2; actionid; |
---|
| 455 | "); |
---|
| 456 | |
---|
| 457 | return(StabEqnId(ideal(f), qhweight(f))); |
---|
| 458 | } |
---|
| 459 | example |
---|
| 460 | {"EXAMPLE:"; echo = 2; |
---|
| 461 | ring B = 0,(x,y,z), ls; |
---|
| 462 | poly f = -z5+y5+x2z+x2y; |
---|
| 463 | list stab = StabEqn(f); |
---|
| 464 | def S1 = stab[1]; setring S1; stabid; |
---|
| 465 | def S2 = stab[2]; setring S2; actionid; |
---|
| 466 | } |
---|
| 467 | |
---|
| 468 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 469 | |
---|
| 470 | proc StabEqnId(ideal data, intvec wt) |
---|
| 471 | "USAGE: StabEqn(I, w); I ideal, w intvec |
---|
[fd5013] | 472 | PURPOSE: compute the equations of the isometry group of the ideal I, |
---|
[35f23d] | 473 | each generator of I is fixed by the stabilizer. |
---|
[3754ca] | 474 | ASSUME: I semiquasihomogeneous ideal w.r.t. 'w' with an isolated singularity at 0 |
---|
[fd5013] | 475 | RETURN: list of two rings 'S1', 'S2' |
---|
[1f92589] | 476 | - 'S1' contians the equations of the stabilizer (ideal 'stabid') @* |
---|
[35f23d] | 477 | - 'S2' contains the action of the stabilizer (ideal 'actionid') |
---|
| 478 | EXAMPLE: example StabEqnId; shows an example |
---|
| 479 | GLOBAL: varSubsList, contains the index j s.t. t(i) -> t(i)t(j) ... |
---|
| 480 | " |
---|
| 481 | { |
---|
[e7778a] | 482 | int i, j, c, k, r, nrVars, offset, n, sln, dbPrt; |
---|
[afd77b2] | 483 | list Variables, rd, temp, sList, varSubsList; |
---|
[35f23d] | 484 | string ringSTR, ringSTR1, varString, parString; |
---|
| 485 | |
---|
| 486 | dbPrt = printlevel-voice+2; |
---|
| 487 | dbprint(dbPrt, "StabilizerEquations of " + string(data)); |
---|
| 488 | |
---|
| 489 | export(varSubsList); |
---|
| 490 | n = nvars(basering); |
---|
[afd77b2] | 491 | Variables = StabVar(wt); // possible quasihomogeneous substitutions |
---|
[35f23d] | 492 | nrVars = 0; |
---|
[c13ff56] | 493 | for(i = 1; i <= size(wt); i++) |
---|
| 494 | { |
---|
[afd77b2] | 495 | nrVars = nrVars + size(Variables[i]); |
---|
[35f23d] | 496 | } |
---|
| 497 | |
---|
| 498 | // set the new basering needed for the substitutions |
---|
| 499 | |
---|
| 500 | varString = "s(1.." + string(nrVars) + ")"; |
---|
[c13ff56] | 501 | if(npars(basering) == 1) |
---|
| 502 | { |
---|
[35f23d] | 503 | parString = "(0, " + parstr(basering) + ")"; |
---|
| 504 | } |
---|
| 505 | else { parString = "0"; } |
---|
| 506 | |
---|
| 507 | def RSTB = basering; |
---|
[393c47] | 508 | string @mPoly = string(minpoly); |
---|
[35f23d] | 509 | ringSTR = "ring RSTR = " + parString + ", (" + varstr(basering) + ", " + varString + "), dp;"; // dp |
---|
| 510 | ringSTR1 = "ring RSTT = " + parString + ", (" + varString + ", " + varstr(basering) + "), dp;"; |
---|
| 511 | |
---|
[3b77465] | 512 | if(defined(RSTR)) { kill RSTR;} |
---|
| 513 | if(defined(RSTT)) { kill RSTT;} |
---|
[35f23d] | 514 | execute(ringSTR1); // this ring is only used for the result, where the variables |
---|
| 515 | export(RSTT); // are s(1..m),t(1..n), as needed for Derksens algorithm (NullCone) |
---|
[393c47] | 516 | execute("minpoly = number(" + @mPoly + ");"); |
---|
[35f23d] | 517 | |
---|
| 518 | execute(ringSTR); |
---|
| 519 | export(RSTR); |
---|
[393c47] | 520 | execute("minpoly = number(" + @mPoly + ");"); |
---|
[35f23d] | 521 | poly f, f1, g, h, vars, pp; // f1 is the polynomial after subs, |
---|
| 522 | ideal allEqns, qhsubs, actionid, stabid, J; |
---|
| 523 | list ringList; // all t(i)`s which do not appear in f1 |
---|
| 524 | ideal data = simplify(imap(RSTB, data), 2); |
---|
| 525 | |
---|
[1f92589] | 526 | // generate the quasihomogeneous substitution map F |
---|
[35f23d] | 527 | |
---|
| 528 | nrVars = 0; |
---|
| 529 | offset = 0; |
---|
[c13ff56] | 530 | for(i = 1; i <= size(wt); i++) |
---|
| 531 | { // build the substitution t(i) -> ... |
---|
[afd77b2] | 532 | if(i > 1) { offset = offset + size(Variables[i - 1]); } |
---|
[35f23d] | 533 | g = 0; |
---|
[c13ff56] | 534 | for(j = 1; j <= size(Variables[i]); j++) |
---|
| 535 | { |
---|
[35f23d] | 536 | pp = 1; |
---|
[c13ff56] | 537 | for(k = 2; k <= size(Variables[i][j]); k++) |
---|
| 538 | { |
---|
[afd77b2] | 539 | pp = pp * var(Variables[i][j][k]); |
---|
| 540 | if(Variables[i][j][k] == i) { varSubsList[i] = offset + j;} |
---|
[35f23d] | 541 | } |
---|
| 542 | g = g + s(offset + j) * pp; |
---|
| 543 | } |
---|
| 544 | qhsubs[i] = g; |
---|
| 545 | } |
---|
| 546 | dbprint(dbPrt, " qhasihomogenous substituion =" + string(qhsubs)); |
---|
| 547 | map F = RSTR, qhsubs; |
---|
[3b77465] | 548 | kill varSubsList; |
---|
[35f23d] | 549 | |
---|
| 550 | // get the equations of the stabilizer by comparing coefficients |
---|
| 551 | // in the equation f = F(f). |
---|
| 552 | |
---|
| 553 | vars = RingVarProduct(Table("i", "i", 1, size(wt))); |
---|
| 554 | |
---|
| 555 | allEqns = 0; |
---|
| 556 | |
---|
[e7778a] | 557 | matrix newcoMx, coMx; |
---|
| 558 | int d; |
---|
| 559 | for(r = 1; r <= ncols(data); r++) |
---|
| 560 | { |
---|
[35f23d] | 561 | |
---|
| 562 | f = data[r]; |
---|
| 563 | f1 = F(f); |
---|
[7de8e4] | 564 | d = deg(f); |
---|
| 565 | newcoMx = coef(f1, vars); // coefficients of F(f) |
---|
| 566 | coMx = coef(f, vars); // coefficients of f |
---|
[35f23d] | 567 | |
---|
[c13ff56] | 568 | for(i = 1; i <= ncols(newcoMx); i++) |
---|
| 569 | { // build the system of eqns via coeff. comp. |
---|
[35f23d] | 570 | j = 1; |
---|
| 571 | h = 0; |
---|
[c13ff56] | 572 | while(j <= ncols(coMx)) |
---|
| 573 | { // all monomials in f |
---|
[35f23d] | 574 | if(coMx[j][1] == newcoMx[i][1]) { h = coMx[j][2]; j = ncols(coMx) + 1;} |
---|
| 575 | else {j = j + 1;} |
---|
| 576 | } |
---|
| 577 | J = J, newcoMx[i][2] - h; // add equation |
---|
| 578 | } |
---|
| 579 | allEqns = allEqns, J; |
---|
| 580 | |
---|
| 581 | } |
---|
| 582 | allEqns = std(allEqns); |
---|
| 583 | |
---|
| 584 | // simplify the equations, i.e., if s(i) in J then remove s(i) from J |
---|
| 585 | // and from the basering |
---|
| 586 | |
---|
| 587 | sList = SimplifyIdeal(allEqns, n, "s"); |
---|
| 588 | stabid = sList[1]; |
---|
| 589 | map phi = basering, sList[2]; |
---|
| 590 | sln = size(sList[3]) - n; |
---|
| 591 | |
---|
| 592 | // change the substitution |
---|
| 593 | |
---|
| 594 | actionid = phi(qhsubs); |
---|
| 595 | |
---|
| 596 | // change to new ring, auxillary construction |
---|
| 597 | |
---|
| 598 | setring(RSTT); |
---|
| 599 | ideal actionid, stabid; |
---|
| 600 | |
---|
| 601 | actionid = imap(RSTR, actionid); |
---|
| 602 | stabid = imap(RSTR, stabid); |
---|
| 603 | export(stabid); |
---|
| 604 | export(actionid); |
---|
| 605 | ringList[2] = RSTT; |
---|
| 606 | |
---|
| 607 | dbprint(dbPrt, " substitution = " + string(actionid)); |
---|
| 608 | dbprint(dbPrt, " equations of stabilizer = " + string(stabid)); |
---|
| 609 | |
---|
| 610 | varString = "s(1.." + string(sln) + ")"; |
---|
| 611 | ringSTR = "ring RSTS = " + parString + ", (" + varString + "), dp;"; |
---|
| 612 | execute(ringSTR); |
---|
[393c47] | 613 | execute("minpoly = number(" + @mPoly + ");"); |
---|
[35f23d] | 614 | ideal stabid = std(imap(RSTR, stabid)); |
---|
| 615 | export(stabid); |
---|
| 616 | ringList[1] = RSTS; |
---|
| 617 | dbprint(dbPrt, " |
---|
| 618 | // 'StabEqnId' created a list of 2 rings. |
---|
| 619 | // To see the rings, type (if the name of your list is stab): |
---|
| 620 | show(stab); |
---|
| 621 | // To access the 1-st ring and map (and similair for the others), type: |
---|
| 622 | def S1 = stab[1]; setring S1; stabid; |
---|
| 623 | // S1/stabid is the coordinate ring of the variety of the |
---|
| 624 | // stabilizer, say G. If G x K^n --> K^n is the action of G on |
---|
| 625 | // K^n, then the ideal 'actionid' in the second ring describes |
---|
| 626 | // the dual map on the ring level. |
---|
| 627 | // To access the 2-nd ring and map (and similair for the others), type: |
---|
| 628 | def S2 = stab[2]; setring S2; actionid; |
---|
| 629 | "); |
---|
| 630 | return(ringList); |
---|
| 631 | } |
---|
| 632 | example |
---|
| 633 | {"EXAMPLE:"; echo = 2; |
---|
| 634 | ring B = 0,(x,y,z), ls; |
---|
| 635 | ideal I = x2,y3,z6; |
---|
| 636 | intvec w = 3,2,1; |
---|
| 637 | list stab = StabEqnId(I, w); |
---|
| 638 | def S1 = stab[1]; setring S1; stabid; |
---|
| 639 | def S2 = stab[2]; setring S2; actionid; |
---|
| 640 | } |
---|
| 641 | |
---|
| 642 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 643 | static |
---|
| 644 | proc ArnoldFormMain(poly f, B, poly Fs, ideal reduceIdeal, int nrs, int nrt) |
---|
| 645 | "USAGE: ArnoldFormMain(f, B, Fs, rI, nrs, nrt); |
---|
| 646 | poly f,Fs; ideal B, rI; int nrs, nrt |
---|
| 647 | PURPOSE: compute the induced action of 'G_f' on T_, where f is the principal |
---|
| 648 | part and 'Fs' is the semiuniversal unfolding of 'f' with x_i |
---|
| 649 | substituted by actionid[i], 'B' is a list of upper basis monomials |
---|
| 650 | for the milnor algebra of 'f', 'nrs' = number of variables for 'G_f' |
---|
| 651 | and 'nrt' = dimension of T_ |
---|
[1f92589] | 652 | ASSUME: f is quasihomogeneous with an isolated singularity at 0, |
---|
[35f23d] | 653 | s(1..r), t(1..m) are parameters of the basering |
---|
| 654 | RETURN: poly |
---|
| 655 | EXAMPLE: example ArnoldAction; shows an example |
---|
| 656 | " |
---|
| 657 | { |
---|
| 658 | int i, j, d, ub, dbPrt; |
---|
| 659 | list upperBasis, basisDegList, gmonos, common, parts; |
---|
| 660 | ideal jacobianId, jacobIdstd, mapId; // needed for phi |
---|
| 661 | intvec wt = weight(f); |
---|
| 662 | matrix gCoeffMx; // for lift command |
---|
| 663 | poly newFs, g, gred, tt; // g = sum of all monomials of degree d, gred is needed for lift |
---|
| 664 | map phi; // the map from Arnold's Theorem |
---|
| 665 | |
---|
| 666 | dbPrt = printlevel-voice+2; |
---|
| 667 | jacobianId = jacob(f); |
---|
| 668 | jacobIdstd = std(jacobianId); |
---|
| 669 | newFs = Fs; |
---|
[c13ff56] | 670 | for(i = 1; i <= size(B); i++) |
---|
| 671 | { |
---|
[35f23d] | 672 | basisDegList[i] = deg(B[i], wt); |
---|
| 673 | } |
---|
| 674 | ub = Max(basisDegList) + 1; // max degree of an upper monomial |
---|
| 675 | |
---|
| 676 | parts = MonosAndTerms(newFs - f, wt, ub); |
---|
| 677 | gmonos = parts[1]; |
---|
| 678 | d = deg(f, wt); |
---|
| 679 | |
---|
[c13ff56] | 680 | for(i = d + 1; i < ub; i++) |
---|
| 681 | { // base[1] = monomials of degree i |
---|
[35f23d] | 682 | upperBasis[i] = SelectMonos(list(B, B), wt, i); // B must not contain 0's |
---|
| 683 | } |
---|
| 684 | |
---|
| 685 | // test if each monomial of Fs is contained in B, if not, |
---|
| 686 | // compute a substitution via Arnold's theorem and substitutite |
---|
| 687 | // it into newFs |
---|
| 688 | |
---|
[c13ff56] | 689 | for(i = d + 1; i < ub; i = i + 1) |
---|
| 690 | { // ub instead of @UB |
---|
[35f23d] | 691 | dbprint(dbPrt, "-- degree = " + string(i) + " of " + string(ub - 1) + " ---------------------------"); |
---|
| 692 | if(size(newFs) < 80) { dbprint(dbPrt, " polynomial = " + string(newFs - f));} |
---|
| 693 | else { dbprint(dbPrt, " poly has deg (not weighted) " + string(deg(newFs)) + " and contains " + string(size(newFs)) + " monos");} |
---|
| 694 | |
---|
| 695 | // select monomials of degree i and intersect them with upperBasis[i] |
---|
| 696 | |
---|
| 697 | gmonos = SelectMonos(parts, wt, i); |
---|
| 698 | common = IntersectionQHM(upperBasis[i][1], gmonos[1]); |
---|
[c13ff56] | 699 | if(size(common) == size(gmonos[1])) |
---|
| 700 | { |
---|
[35f23d] | 701 | dbprint(dbPrt, " no additional monomials "); |
---|
| 702 | } |
---|
| 703 | |
---|
| 704 | // other monomials than those in upperBasis occur, compute |
---|
| 705 | // the map constructed in the proof of Arnold's theorem |
---|
| 706 | // write g = c[i] * jacobianId[i] |
---|
| 707 | |
---|
[c13ff56] | 708 | else |
---|
| 709 | { |
---|
[35f23d] | 710 | dbprint(dbPrt, " additional Monomials found, compute the map "); |
---|
| 711 | g = PSum(gmonos[2]); // sum of all monomials in g of degree i |
---|
| 712 | dbprint(dbPrt, " sum of degree " + string(i) + " is " + string(g)); |
---|
| 713 | |
---|
| 714 | gred = reduce(g, jacobIdstd); |
---|
| 715 | gCoeffMx = lift(jacobianId, g - gred); // compute c[i] |
---|
| 716 | mapId = var(1) - gCoeffMx[1][1]; // generate the map |
---|
[c13ff56] | 717 | for(j = 2; j <= size(gCoeffMx); j++) |
---|
| 718 | { |
---|
[35f23d] | 719 | mapId[j] = var(j) - gCoeffMx[1][j]; |
---|
| 720 | } |
---|
| 721 | dbprint(dbPrt, " map = " + string(mapId)); |
---|
| 722 | // apply the map to newFs |
---|
| 723 | newFs = APSubstitution(newFs, mapId, reduceIdeal, wt, ub, nrs, nrt); |
---|
| 724 | parts = MonosAndTerms(newFs - f, wt, ub); // monos and terms of deg < ub |
---|
| 725 | newFs = PSum(parts[2]) + f; // result of APS... is already reduced |
---|
| 726 | dbprint(dbPrt, " monomials of degree " + string(i)); |
---|
| 727 | } |
---|
| 728 | } |
---|
| 729 | return(newFs); |
---|
| 730 | } |
---|
| 731 | |
---|
| 732 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 733 | |
---|
[70597c] | 734 | static proc MonosAndTerms(poly f, wt, int ub) |
---|
[35f23d] | 735 | "USAGE: MonosAndTerms(f, w, ub); poly f, intvec w, int ub |
---|
| 736 | PURPOSE: returns a list of all monomials and terms occuring in f of |
---|
| 737 | weighted degree < ub |
---|
| 738 | RETURN: list |
---|
| 739 | _[1] list of monomials |
---|
[8942a5] | 740 | _[2] list of terms |
---|
[35f23d] | 741 | EXAMPLE: example MonosAndTerms shows an example |
---|
| 742 | " |
---|
| 743 | { |
---|
| 744 | int i, k; |
---|
| 745 | list monomials, terms; |
---|
| 746 | poly mono, lcInv, data; |
---|
| 747 | |
---|
| 748 | data = jet(f, ub - 1, wt); |
---|
| 749 | k = 0; |
---|
[c13ff56] | 750 | for(i = 1; i <= size(data); i++) |
---|
| 751 | { |
---|
[35f23d] | 752 | mono = lead(data[i]); |
---|
[c13ff56] | 753 | if(deg(mono, wt) < ub) |
---|
| 754 | { |
---|
[35f23d] | 755 | k = k + 1; |
---|
| 756 | lcInv = 1/leadcoef(mono); |
---|
| 757 | monomials[k] = mono * lcInv; |
---|
| 758 | terms[k] = mono; |
---|
| 759 | } |
---|
| 760 | } |
---|
| 761 | return(list(monomials, terms)); |
---|
| 762 | } |
---|
| 763 | example |
---|
| 764 | {"EXAMPLE:"; echo = 2; |
---|
| 765 | ring B = 0,(x,y,z), lp; |
---|
| 766 | poly f = 6*x2 + 2*x3 + 9*x*y2 + z*y + x*z6; |
---|
| 767 | MonosAndTerms(f, intvec(2,1,1), 5); |
---|
| 768 | } |
---|
| 769 | |
---|
| 770 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 771 | |
---|
[70597c] | 772 | static proc SelectMonos(parts, intvec wt, int d) |
---|
[35f23d] | 773 | "USAGE: SelectMonos(parts, w, d); list/ideal parts, intvec w, int d |
---|
| 774 | PURPOSE: returns a list of all monomials and terms occuring in f of |
---|
| 775 | weighted degree = d |
---|
| 776 | RETURN: list |
---|
[8942a5] | 777 | _[1] list of monomials |
---|
| 778 | _[2] list of terms |
---|
[35f23d] | 779 | EXAMPLE: example SelectMonos; shows an example |
---|
| 780 | " |
---|
| 781 | { |
---|
| 782 | int i, k; |
---|
| 783 | list monomials, terms; |
---|
| 784 | poly mono; |
---|
| 785 | |
---|
| 786 | k = 0; |
---|
[c13ff56] | 787 | for(i = 1; i <= size(parts[1]); i++) |
---|
| 788 | { |
---|
[35f23d] | 789 | mono = parts[1][i]; |
---|
[c13ff56] | 790 | if(deg(mono, wt) == d) |
---|
| 791 | { |
---|
| 792 | k++; |
---|
[35f23d] | 793 | monomials[k] = mono; |
---|
| 794 | terms[k] = parts[2][i]; |
---|
| 795 | } |
---|
| 796 | } |
---|
| 797 | return(list(monomials, terms)); |
---|
| 798 | } |
---|
| 799 | example |
---|
| 800 | {"EXAMPLE:"; echo = 2; |
---|
| 801 | ring B = 0,(x,y,z), lp; |
---|
| 802 | poly f = 6*x2 + 2*x3 + 9*x*y2 + z*y + x*z6; |
---|
| 803 | list mt = MonosAndTerms(f, intvec(2,1,1), 5); |
---|
| 804 | SelectMonos(mt, intvec(2,1,1), 4); |
---|
| 805 | } |
---|
| 806 | |
---|
| 807 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 808 | |
---|
[70597c] | 809 | static proc Expand(substitution, degVec, ideal reduceI, intvec w1, int ub, list truncated) |
---|
[35f23d] | 810 | "USAGE: Expand(substitution, degVec, reduceI, w, ub, truncated); |
---|
| 811 | ideal/list substitution, list/intvec degVec, ideal reduceI, intvec w, |
---|
| 812 | int ub, list truncated |
---|
| 813 | PURPOSE: substitute 'substitution' in the monomial given by the list of |
---|
| 814 | exponents 'degVec', omit all terms of weighted degree > ub and reduce |
---|
| 815 | the result w.r.t. 'reduceI'. If truncated[i] = 0 then the result is |
---|
| 816 | stored for later use. |
---|
| 817 | RETURN: poly |
---|
| 818 | NOTE: used by APSubstitution |
---|
| 819 | GLOBAL: computedPowers |
---|
| 820 | " |
---|
| 821 | { |
---|
| 822 | int i, minDeg; |
---|
| 823 | list powerList; |
---|
| 824 | poly g, h; |
---|
| 825 | |
---|
| 826 | // compute substitution[1]^degVec[1],...,subs[n]^degVec[n] |
---|
| 827 | |
---|
[c13ff56] | 828 | for(i = 1; i <= ncols(substitution); i++) |
---|
| 829 | { |
---|
| 830 | if(size(substitution[i]) < 3 || degVec[i] < 4) |
---|
| 831 | { |
---|
[35f23d] | 832 | powerList[i] = reduce(substitution[i]^degVec[i], reduceI); // new |
---|
| 833 | } // directly for small exponents |
---|
[c13ff56] | 834 | else |
---|
| 835 | { |
---|
[35f23d] | 836 | powerList[i] = PolyPower1(i, substitution[i], degVec[i], reduceI, w1, truncated[i], ub); |
---|
| 837 | } |
---|
| 838 | } |
---|
| 839 | // multiply the terms obtained by using PolyProduct(); |
---|
| 840 | g = powerList[1]; |
---|
| 841 | minDeg = w1[1] * degVec[1]; |
---|
[c13ff56] | 842 | for(i = 2; i <= ncols(substitution); i++) |
---|
| 843 | { |
---|
[35f23d] | 844 | g = jet(g, ub - w1[i] * degVec[i] - 1, w1); |
---|
| 845 | h = jet(powerList[i], ub - minDeg - 1, w1); |
---|
| 846 | g = PolyProduct(g, h, reduceI, w1, ub); |
---|
| 847 | if(g == 0) { Print(" g = 0 "); break;} |
---|
| 848 | minDeg = minDeg + w1[i] * degVec[i]; |
---|
| 849 | } |
---|
| 850 | return(g); |
---|
| 851 | } |
---|
| 852 | |
---|
| 853 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 854 | |
---|
[70597c] | 855 | static proc PolyProduct(poly g1, poly h1, ideal reduceI, intvec wt, int ub) |
---|
[35f23d] | 856 | "USAGE: PolyProduct(g, h, reduceI, wt, ub); poly g, h; ideal reduceI, |
---|
| 857 | intvec wt, int ub. |
---|
| 858 | PURPOSE: compute g*h and reduce it w.r.t 'reduceI' and omit terms of weighted |
---|
| 859 | degree > ub. |
---|
| 860 | RETURN: poly |
---|
| 861 | NOTE: used by 'Expand' |
---|
| 862 | " |
---|
| 863 | { |
---|
[c13ff56] | 864 | int SUBSMAXSIZE = 3000; |
---|
[35f23d] | 865 | int i, nrParts, sizeOfPart, currentPos, partSize, maxSIZE; |
---|
| 866 | poly g, h, gxh, prodComp, @g2; // replace @g2 by subst. |
---|
| 867 | |
---|
| 868 | g = g1; |
---|
| 869 | h = h1; |
---|
| 870 | |
---|
[c13ff56] | 871 | if(size(g)*size(h) > SUBSMAXSIZE) |
---|
| 872 | { |
---|
[35f23d] | 873 | // divide the polynomials with more terms in parts s.t. |
---|
| 874 | // the product of each part with the other polynomial |
---|
| 875 | // has at most SUBMAXSIZE terms |
---|
| 876 | |
---|
| 877 | if(size(g) < size(h)) { poly @h = h; h = g; g = @h;@h = 0; } |
---|
| 878 | maxSIZE = SUBSMAXSIZE / size(h); |
---|
| 879 | //print(" SUBSMAXSIZE = "+string(SUBSMAXSIZE)+" exceeded by "+string(size(g)*size(h)) + ", maxSIZE = ", string(maxSIZE)); |
---|
[1f9a84] | 880 | nrParts = size(g) div maxSIZE + 1; |
---|
| 881 | partSize = size(g) div nrParts; |
---|
[35f23d] | 882 | gxh = 0; // 'g times h' |
---|
[c13ff56] | 883 | for(i = 1; i <= nrParts; i++) |
---|
| 884 | { |
---|
[35f23d] | 885 | //print(" loop #" + string(i) + " of " + string(nrParts)); |
---|
| 886 | currentPos = (i - 1) * partSize; |
---|
| 887 | if(i < nrParts) {sizeOfPart = partSize;} |
---|
| 888 | else { sizeOfPart = size(g) - (nrParts - 1) * partSize; print(" last #" + string(sizeOfPart) + " terms ");} |
---|
| 889 | prodComp = g[currentPos + 1..sizeOfPart + currentPos] * h; // multiply a part |
---|
| 890 | @g2 = jet(prodComp, ub - 1, wt); // eventual reduce ... |
---|
| 891 | if(size(@g2) < size(prodComp)) { print(" killed " + string(size(prodComp) - size(@g2)) + " terms ");} |
---|
| 892 | gxh = reduce(gxh + @g2, reduceI); |
---|
| 893 | } |
---|
| 894 | } |
---|
[c13ff56] | 895 | else |
---|
| 896 | { |
---|
[35f23d] | 897 | gxh = reduce(jet(g * h,ub - 1, wt), reduceI); |
---|
| 898 | } // compute directly |
---|
| 899 | return(gxh); |
---|
| 900 | } |
---|
| 901 | |
---|
| 902 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 903 | |
---|
[3c4dcc] | 904 | static proc PolyPower1(int varIndex, poly f, int e, ideal reduceI, intvec wt, |
---|
[70597c] | 905 | int truncated, int ub) |
---|
[35f23d] | 906 | "USAGE: PolyPower1(i, f, e, reduceI, wt, truncated, ub);int i, e, ub;poly f; |
---|
| 907 | ideal reduceI; intvec wt; list truncated; |
---|
| 908 | PURPOSE: compute f^e, use previous computations if possible, and reduce it |
---|
| 909 | w.r.t reudecI and omit terms of weighted degree > ub. |
---|
| 910 | RETURN: poly |
---|
| 911 | NOTE: used by 'Expand' |
---|
| 912 | GLOBAL: 'computedPowers' |
---|
| 913 | " |
---|
| 914 | { |
---|
| 915 | int i, ordOfg, lb, maxPrecomputedPower; |
---|
| 916 | poly g, fn; |
---|
| 917 | |
---|
| 918 | if(e == 0) { return(1);} |
---|
| 919 | if(e == 1) { return(f);} |
---|
| 920 | if(f == 0) { return(1); } |
---|
[c13ff56] | 921 | else |
---|
| 922 | { |
---|
[35f23d] | 923 | // test if f has been computed to some power |
---|
[c13ff56] | 924 | if(computedPowers[varIndex][1] > 0) |
---|
| 925 | { |
---|
[35f23d] | 926 | maxPrecomputedPower = computedPowers[varIndex][1]; |
---|
[c13ff56] | 927 | if(maxPrecomputedPower >= e) |
---|
| 928 | { |
---|
[35f23d] | 929 | // no computation necessary, f^e has already benn computed |
---|
| 930 | g = computedPowers[varIndex][2][e - 1]; |
---|
| 931 | //Print("No computation, from list : g = elem [", varIndex, ", 2, ", e - 1, "]"); |
---|
| 932 | lb = e + 1; |
---|
| 933 | } |
---|
| 934 | else { // f^d computed, where d < e |
---|
| 935 | g = computedPowers[varIndex][2][maxPrecomputedPower - 1]; |
---|
| 936 | ordOfg = maxPrecomputedPower * wt[varIndex]; |
---|
| 937 | lb = maxPrecomputedPower + 1; |
---|
| 938 | } |
---|
| 939 | } |
---|
[c13ff56] | 940 | else |
---|
| 941 | { // no precomputed data |
---|
[35f23d] | 942 | lb = 2; |
---|
| 943 | ordOfg = wt[varIndex]; |
---|
| 944 | g = f; |
---|
| 945 | } |
---|
[c13ff56] | 946 | for(i = lb; i <= e; i++) |
---|
| 947 | { |
---|
[35f23d] | 948 | fn = jet(f, ub - ordOfg - 1, wt); // reduce w.r.t. reduceI |
---|
| 949 | g = PolyProduct(g, fn, reduceI, wt, ub); |
---|
| 950 | ordOfg = ordOfg + wt[varIndex]; |
---|
| 951 | if(g == 0) { break; } |
---|
[c13ff56] | 952 | if((i > maxPrecomputedPower) && !truncated) |
---|
| 953 | { |
---|
| 954 | if(maxPrecomputedPower == 0) |
---|
[0610f0e] | 955 | { // init computedPowers |
---|
[35f23d] | 956 | computedPowers[varIndex] = list(i, list(g)); |
---|
| 957 | } |
---|
| 958 | computedPowers[varIndex][1] = i; // new degree |
---|
| 959 | computedPowers[varIndex][2][i - 1] = g; |
---|
| 960 | maxPrecomputedPower = i; |
---|
| 961 | } |
---|
| 962 | } |
---|
| 963 | } |
---|
| 964 | return(g); |
---|
| 965 | } |
---|
| 966 | |
---|
| 967 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 968 | |
---|
[70597c] | 969 | static proc RingVarsToList(list @index) |
---|
[35f23d] | 970 | { |
---|
| 971 | int i; |
---|
| 972 | list temp; |
---|
| 973 | |
---|
[c13ff56] | 974 | for(i = 1; i <= size(@index); i++) { temp[i] = string(var(@index[i])); } |
---|
[35f23d] | 975 | return(temp); |
---|
| 976 | } |
---|
| 977 | |
---|
| 978 | /////////////////////////////////////////////////////////////////////////////// |
---|
[70597c] | 979 | static |
---|
[35f23d] | 980 | proc APSubstitution(poly f, ideal substitution, ideal reduceIdeal, intvec wt, int ub, int nrs, int nrt) |
---|
| 981 | "USAGE: APSubstitution(f, subs, reduceI, w, ub, int nrs, int nrt); poly f |
---|
| 982 | ideal subs, reduceI, intvec w, int ub, nrs, nrt; |
---|
| 983 | nrs = number of parameters s(1..nrs), |
---|
| 984 | nrt = number of parameters t(1..nrt) |
---|
| 985 | PURPOSE: substitute 'subs' in f, omit all terms with weighted degree > ub and |
---|
| 986 | reduce the result w.r.t. 'reduceI'. |
---|
| 987 | RETURN: poly |
---|
| 988 | GLOBAL: 'computedPowers' |
---|
| 989 | " |
---|
| 990 | { |
---|
| 991 | int i, j, k, d, offset; |
---|
| 992 | int n = nvars(basering); |
---|
| 993 | list coeffList, parts, degVecList, degOfMonos; |
---|
| 994 | list computedPowers, truncatedQ, degOfSubs, @temp; |
---|
| 995 | string ringSTR, @ringVars; |
---|
| 996 | |
---|
| 997 | export(computedPowers); |
---|
| 998 | |
---|
| 999 | // store arguments in strings |
---|
| 1000 | |
---|
| 1001 | def RASB = basering; |
---|
| 1002 | |
---|
| 1003 | parts = MonosAndTerms(f, wt, ub); |
---|
[c13ff56] | 1004 | for(i = 1; i <= size(parts[1]); i = i + 1) |
---|
| 1005 | { |
---|
[35f23d] | 1006 | coeffList[i] = parts[2][i]/parts[1][i]; |
---|
| 1007 | degVecList[i] = leadexp(parts[1][i]); |
---|
| 1008 | degOfMonos[i] = deg(parts[1][i], wt); |
---|
| 1009 | } |
---|
| 1010 | |
---|
| 1011 | // built new basering with no parameters and order dp ! |
---|
| 1012 | // the parameters of the basering are appended to |
---|
| 1013 | // the variables of the basering ! |
---|
| 1014 | // set ideal mpoly = minpoly for reduction ! |
---|
| 1015 | |
---|
| 1016 | @ringVars = "(" + varstr(basering) + ", " + parstr(1) + ","; // precondition |
---|
[c13ff56] | 1017 | if(nrs > 0) |
---|
| 1018 | { |
---|
[35f23d] | 1019 | @ringVars = @ringVars + "s(1.." + string(nrs) + "), "; |
---|
| 1020 | } |
---|
| 1021 | @ringVars = @ringVars + "t(1.." + string(nrt) + "))"; |
---|
| 1022 | ringSTR = "ring RASR = 0, " + @ringVars + ", dp;"; // new basering |
---|
| 1023 | |
---|
| 1024 | // built the "reduction" ring with the reduction ideal |
---|
| 1025 | |
---|
| 1026 | execute(ringSTR); |
---|
| 1027 | export(RASR); |
---|
| 1028 | ideal reduceIdeal, substitution, newSubs; |
---|
| 1029 | intvec w1, degVec; |
---|
| 1030 | list minDeg, coeffList, degList; |
---|
| 1031 | poly f, g, h, subsPoly; |
---|
| 1032 | |
---|
| 1033 | w1 = wt; // new weights |
---|
| 1034 | offset = nrs + nrt + 1; |
---|
| 1035 | for(i = n + 1; i <= offset + n; i = i + 1) { w1[i] = 0; } |
---|
| 1036 | |
---|
| 1037 | reduceIdeal = std(imap(RASB, reduceIdeal)); // omit later ! |
---|
| 1038 | coeffList = imap(RASB, coeffList); |
---|
| 1039 | substitution = imap(RASB, substitution); |
---|
| 1040 | |
---|
| 1041 | f = imap(RASB, f); |
---|
| 1042 | |
---|
[c13ff56] | 1043 | for(i = 1; i <= n; i++) |
---|
| 1044 | { // all "base" variables |
---|
[35f23d] | 1045 | computedPowers[i] = list(0); |
---|
[c13ff56] | 1046 | for(j = 1; j <= size(substitution[i]); j++) { degList[j] = deg(substitution[i][j], w1);} |
---|
[35f23d] | 1047 | degOfSubs[i] = degList; |
---|
| 1048 | } |
---|
| 1049 | |
---|
| 1050 | // substitute in each monomial seperately |
---|
| 1051 | |
---|
| 1052 | g = 0; |
---|
[c13ff56] | 1053 | for(i = 1; i <= size(degVecList); i++) |
---|
| 1054 | { |
---|
[35f23d] | 1055 | truncatedQ = Table("0", "i", 1, n); |
---|
| 1056 | newSubs = 0; |
---|
| 1057 | degVec = degVecList[i]; |
---|
| 1058 | d = degOfMonos[i]; |
---|
| 1059 | |
---|
| 1060 | // check if some terms in the substitution can be omitted |
---|
| 1061 | // degVec = list of exponents of the monomial m |
---|
| 1062 | // minDeg[j] denotes the weighted degree of the monomial m' |
---|
| 1063 | // where m' is the monomial m without the j-th variable |
---|
| 1064 | |
---|
[c13ff56] | 1065 | for(j = 1; j <= size(degVec); j++) { minDeg[j] = d - degVec[j] * wt[j]; } |
---|
[35f23d] | 1066 | |
---|
[c13ff56] | 1067 | for(j = 1; j <= size(degVec); j++) |
---|
| 1068 | { |
---|
[35f23d] | 1069 | subsPoly = 0; // set substitution to 0 |
---|
[c13ff56] | 1070 | if(degVec[j] > 0) |
---|
| 1071 | { |
---|
[35f23d] | 1072 | // if variable occurs then check if |
---|
| 1073 | // substitution[j][k] * (linear part)^(degVec[j]-1) + minDeg[j] < ub |
---|
| 1074 | // i.e. look for the smallest possible combination in subs[j]^degVec[j] |
---|
| 1075 | // which comes from the term substitution[j][k]. This term is multiplied |
---|
| 1076 | // with the rest of the monomial, which has at least degree minDeg[j]. |
---|
| 1077 | // If the degree of this product is < ub then substitution[j][k] contributes |
---|
| 1078 | // to the result and cannot be omitted |
---|
| 1079 | |
---|
[c13ff56] | 1080 | for(k = 1; k <= size(substitution[j]); k++) |
---|
[0610f0e] | 1081 | { |
---|
[c13ff56] | 1082 | if(degOfSubs[j][k] + (degVec[j] - 1) * wt[j] + minDeg[j] < ub) |
---|
[0610f0e] | 1083 | { |
---|
[35f23d] | 1084 | subsPoly = subsPoly + substitution[j][k]; |
---|
| 1085 | } |
---|
| 1086 | } |
---|
| 1087 | } |
---|
| 1088 | newSubs[j] = subsPoly; // set substitution |
---|
| 1089 | if(substitution[j] - subsPoly != 0) { truncatedQ[j] = 1;} // mark that substitution[j] is truncated |
---|
| 1090 | } |
---|
| 1091 | h = Expand(newSubs, degVec, reduceIdeal, w1, ub, truncatedQ) * coeffList[i]; // already reduced |
---|
| 1092 | g = reduce(g + h, reduceIdeal); |
---|
| 1093 | } |
---|
[3b77465] | 1094 | kill computedPowers; |
---|
[35f23d] | 1095 | setring RASB; |
---|
| 1096 | poly fnew = imap(RASR, g); |
---|
[3b77465] | 1097 | kill RASR; |
---|
[35f23d] | 1098 | return(fnew); |
---|
| 1099 | } |
---|
| 1100 | |
---|
| 1101 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1102 | |
---|
| 1103 | static proc StabVar(intvec wt) |
---|
| 1104 | "USAGE: StabVar(w); intvec w |
---|
[1f92589] | 1105 | PURPOSE: compute the indicies for quasihomogeneous substitutions of each |
---|
[35f23d] | 1106 | variable. |
---|
[1f92589] | 1107 | ASSUME: f semiquasihomogeneous polynomial with an isolated singularity at 0 |
---|
[35f23d] | 1108 | RETURN: list |
---|
| 1109 | _[i] list of combinations for var(i) (i must be appended |
---|
| 1110 | to each comb) |
---|
| 1111 | GLOBAL: 'varSubsList', contains the index j s.t. x(i) -> x(i)t(j) ... |
---|
| 1112 | " |
---|
| 1113 | { |
---|
| 1114 | int i, j, k, uw, ic; |
---|
[afd77b2] | 1115 | list varList, Variables, subs; |
---|
[35f23d] | 1116 | string str, varString; |
---|
| 1117 | |
---|
| 1118 | varList = StabVarComb(wt); |
---|
[c13ff56] | 1119 | for(i = 1; i <= size(wt); i = i + 1) |
---|
| 1120 | { |
---|
[35f23d] | 1121 | subs = 0; |
---|
| 1122 | // built linear substituitons |
---|
[c13ff56] | 1123 | for(j = 1; j <= size(varList[1][i]); j++) |
---|
| 1124 | { |
---|
[35f23d] | 1125 | subs[j] = list(i) + list(varList[1][i][j]); |
---|
| 1126 | } |
---|
[afd77b2] | 1127 | Variables[i] = subs; |
---|
[c13ff56] | 1128 | if(size(varList[2][i]) > 0) |
---|
| 1129 | { |
---|
[35f23d] | 1130 | // built nonlinear substituitons |
---|
| 1131 | subs = 0; |
---|
[c13ff56] | 1132 | for(j = 1; j <= size(varList[2][i]); j++) |
---|
| 1133 | { |
---|
[35f23d] | 1134 | subs[j] = list(i) + varList[2][i][j]; |
---|
| 1135 | } |
---|
[afd77b2] | 1136 | Variables[i] = Variables[i] + subs; |
---|
[35f23d] | 1137 | } |
---|
| 1138 | } |
---|
[afd77b2] | 1139 | return(Variables); |
---|
[35f23d] | 1140 | } |
---|
| 1141 | |
---|
| 1142 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1143 | |
---|
| 1144 | static proc StabVarComb(intvec wt) |
---|
| 1145 | "USAGE: StabVarComb(w); intvec w |
---|
| 1146 | PURPOSE: list all possible indices of indeterminates for a quasihom. subs. |
---|
| 1147 | RETURN: list |
---|
| 1148 | _[1] linear substitutions |
---|
| 1149 | _[2] nonlinear substiutions |
---|
| 1150 | GLOBAL: 'varSubsList', contains the index j s.t. x(i) -> x(i)t(j) ... |
---|
| 1151 | " |
---|
| 1152 | { |
---|
[c13ff56] | 1153 | int mmi, mma, ii, j, k, uw, ic; |
---|
[35f23d] | 1154 | list index, indices, usedWeights, combList, combinations; |
---|
| 1155 | list linearSubs, nonlinearSubs; |
---|
| 1156 | list partitions, subs, temp; // subs[i] = substitution for var(i) |
---|
| 1157 | |
---|
| 1158 | linearSubs = Table("0", "i", 1, size(wt)); |
---|
| 1159 | nonlinearSubs = Table("0", "i", 1, size(wt)); |
---|
| 1160 | |
---|
| 1161 | uw = 0; |
---|
| 1162 | ic = 0; |
---|
[c13ff56] | 1163 | mmi = Min(wt); |
---|
| 1164 | mma = Max(wt); |
---|
[35f23d] | 1165 | |
---|
[c13ff56] | 1166 | for(ii = mmi; ii <= mma; ii++) |
---|
| 1167 | { |
---|
[9ea76a] | 1168 | if(containedQ(wt, ii)) |
---|
[c13ff56] | 1169 | { // find variables of weight ii |
---|
[35f23d] | 1170 | k = 0; |
---|
| 1171 | index = 0; |
---|
| 1172 | // collect the indices of all variables of weight i |
---|
[c13ff56] | 1173 | for(j = 1; j <= size(wt); j++) |
---|
| 1174 | { |
---|
| 1175 | if(wt[j] == ii) |
---|
| 1176 | { |
---|
| 1177 | k++; |
---|
[35f23d] | 1178 | index[k] = j; |
---|
| 1179 | } |
---|
| 1180 | } |
---|
[c13ff56] | 1181 | uw++; |
---|
| 1182 | usedWeights[uw] = ii; |
---|
| 1183 | ic++; |
---|
| 1184 | indices[ii] = index; |
---|
[35f23d] | 1185 | |
---|
| 1186 | // linear part of the substitution |
---|
| 1187 | |
---|
[c13ff56] | 1188 | for(j = 1; j <= size(index); j++) |
---|
| 1189 | { |
---|
[35f23d] | 1190 | linearSubs[index[j]] = index; |
---|
| 1191 | } |
---|
| 1192 | |
---|
| 1193 | // nonlinear part of the substitution |
---|
| 1194 | |
---|
[c13ff56] | 1195 | if(uw > 1) |
---|
| 1196 | { // variables of least weight do not allow nonlinear subs. |
---|
[35f23d] | 1197 | |
---|
[c13ff56] | 1198 | partitions = Partitions(ii, delete(usedWeights, uw)); |
---|
| 1199 | for(j = 1; j <= size(partitions); j++) |
---|
| 1200 | { |
---|
| 1201 | combinations[j] = AllCombinations(partitions[j], indices); |
---|
| 1202 | } |
---|
| 1203 | for(j = 1; j <= size(index); j++) |
---|
| 1204 | { |
---|
| 1205 | nonlinearSubs[index[j]] = FlattenQHM(combinations); // flatten one level ! |
---|
| 1206 | } |
---|
[35f23d] | 1207 | } |
---|
| 1208 | } |
---|
| 1209 | } |
---|
[9ea76a] | 1210 | "type of i:", typeof(i); |
---|
[35f23d] | 1211 | combList[1] = linearSubs; |
---|
| 1212 | combList[2] = nonlinearSubs; |
---|
| 1213 | return(combList); |
---|
| 1214 | } |
---|
| 1215 | |
---|
| 1216 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1217 | |
---|
[70597c] | 1218 | static proc AllCombinations(list partition, list indices) |
---|
[35f23d] | 1219 | "USAGE: AllCombinations(partition,indices); list partition, indices) |
---|
[8942a5] | 1220 | PURPOSE: all combinations for a given partititon |
---|
[35f23d] | 1221 | RETURN: list |
---|
| 1222 | GLOBAL: varSubsList, contains the index j s.t. x(i) -> x(i)t(j) ... |
---|
| 1223 | " |
---|
| 1224 | { |
---|
| 1225 | int i, k, m, ok, p, offset; |
---|
| 1226 | list nrList, indexList; |
---|
| 1227 | |
---|
| 1228 | k = 0; |
---|
| 1229 | offset = 0; |
---|
| 1230 | i = 1; |
---|
| 1231 | ok = 1; |
---|
| 1232 | m = partition[1]; |
---|
[c13ff56] | 1233 | while(ok) |
---|
| 1234 | { |
---|
| 1235 | if(i > size(partition)) |
---|
| 1236 | { |
---|
[35f23d] | 1237 | ok = 0; |
---|
| 1238 | p = 0; |
---|
| 1239 | } |
---|
| 1240 | else { p = partition[i];} |
---|
| 1241 | if(p == m) { i = i + 1;} |
---|
[c13ff56] | 1242 | else |
---|
| 1243 | { |
---|
[35f23d] | 1244 | k = k + 1; |
---|
| 1245 | nrList[k] = i - 1 - offset; |
---|
| 1246 | offset = offset + i - 1; |
---|
| 1247 | indexList[k] = indices[m]; |
---|
| 1248 | if(ok) { m = partition[i];} |
---|
| 1249 | } |
---|
| 1250 | } |
---|
| 1251 | return(AllCombinationsAux(nrList, indexList)); |
---|
| 1252 | } |
---|
| 1253 | |
---|
| 1254 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1255 | |
---|
[70597c] | 1256 | static proc AllSingleCombinations(int n, list index) |
---|
[35f23d] | 1257 | "USAGE: AllSingleCombinations(n index); int n, list index |
---|
| 1258 | PURPOSE: all combinations for var(n) |
---|
| 1259 | RETURN: list |
---|
| 1260 | " |
---|
| 1261 | { |
---|
| 1262 | int i, j, k; |
---|
| 1263 | list comb, newC, temp, newIndex; |
---|
| 1264 | |
---|
[c13ff56] | 1265 | if(n == 1) |
---|
| 1266 | { |
---|
| 1267 | for(i = 1; i <= size(index); i++) |
---|
| 1268 | { |
---|
[35f23d] | 1269 | temp = index[i]; |
---|
| 1270 | comb[i] = temp; |
---|
| 1271 | } |
---|
| 1272 | return(comb); |
---|
| 1273 | } |
---|
[c13ff56] | 1274 | if(size(index) == 1) |
---|
| 1275 | { |
---|
[35f23d] | 1276 | temp = Table(string(index[1]), "i", 1, n); |
---|
| 1277 | comb[1] = temp; |
---|
| 1278 | return(comb); |
---|
| 1279 | } |
---|
| 1280 | newIndex = index; |
---|
[c13ff56] | 1281 | for(i = 1; i <= size(index); i = i + 1) |
---|
| 1282 | { |
---|
[35f23d] | 1283 | if(i > 1) { newIndex = delete(newIndex, 1); } |
---|
| 1284 | newC = AllSingleCombinations(n - 1, newIndex); |
---|
| 1285 | k = size(comb); |
---|
| 1286 | temp = 0; |
---|
[c13ff56] | 1287 | for(j = 1; j <= size(newC); j++) |
---|
| 1288 | { |
---|
[35f23d] | 1289 | temp[1] = index[i]; |
---|
| 1290 | temp = temp + newC[j]; |
---|
| 1291 | comb[k + j] = temp; |
---|
| 1292 | temp = 0; |
---|
| 1293 | } |
---|
| 1294 | } |
---|
| 1295 | return(comb); |
---|
| 1296 | } |
---|
| 1297 | |
---|
| 1298 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1299 | |
---|
[70597c] | 1300 | static proc AllCombinationsAux(list parts, list index) |
---|
[35f23d] | 1301 | "USAGE: AllCombinationsAux(parts ,index); list parts, index |
---|
| 1302 | PURPOSE: all compbinations for nonlinear substituiton |
---|
| 1303 | RETURN: list |
---|
| 1304 | " |
---|
| 1305 | { |
---|
| 1306 | int i, j, k; |
---|
| 1307 | list comb, firstC, restC; |
---|
| 1308 | |
---|
| 1309 | if(size(parts) == 0 || size(index) == 0) { return(comb);} |
---|
| 1310 | |
---|
| 1311 | firstC = AllSingleCombinations(parts[1], index[1]); |
---|
| 1312 | restC = AllCombinationsAux(delete(parts, 1), delete(index, 1)); |
---|
| 1313 | |
---|
| 1314 | if(size(restC) == 0) { comb = firstC;} |
---|
[c13ff56] | 1315 | else |
---|
| 1316 | { |
---|
| 1317 | for(i = 1; i <= size(firstC); i++) |
---|
| 1318 | { |
---|
[35f23d] | 1319 | k = size(comb); |
---|
[c13ff56] | 1320 | for(j = 1; j <= size(restC); j++) |
---|
| 1321 | { |
---|
[35f23d] | 1322 | //elem = firstC[i] + restC[j]; |
---|
| 1323 | // comb[k + j] = elem; |
---|
| 1324 | comb[k + j] = firstC[i] + restC[j]; |
---|
| 1325 | } |
---|
| 1326 | } |
---|
| 1327 | } |
---|
| 1328 | return(comb); |
---|
| 1329 | } |
---|
| 1330 | |
---|
| 1331 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1332 | |
---|
[70597c] | 1333 | static proc Partitions(int n, list nr) |
---|
[35f23d] | 1334 | "USAGE: Partitions(n, nr); int n, list nr |
---|
| 1335 | PURPOSE: partitions of n consisting of elements from nr |
---|
| 1336 | RETURN: list |
---|
| 1337 | " |
---|
| 1338 | { |
---|
| 1339 | int i, j, k; |
---|
| 1340 | list parts, temp, restP, newP, decP; |
---|
| 1341 | |
---|
| 1342 | if(size(nr) == 0) { return(list());} |
---|
[c13ff56] | 1343 | if(size(nr) == 1) |
---|
| 1344 | { |
---|
| 1345 | if(NumFactor(nr[1], n) > 0) |
---|
| 1346 | { |
---|
[35f23d] | 1347 | parts[1] = Table(string(nr[1]), "i", 1, NumFactor(nr[1], n)); |
---|
| 1348 | } |
---|
| 1349 | return(parts); |
---|
| 1350 | } |
---|
[c13ff56] | 1351 | else |
---|
| 1352 | { |
---|
[35f23d] | 1353 | parts = Partitions(n, nr[1]); |
---|
| 1354 | restP = Partitions(n, delete(nr, 1)); |
---|
| 1355 | |
---|
| 1356 | parts = parts + restP; |
---|
[1f9a84] | 1357 | for(i = 1; i <= n div nr[1]; i = i + 1) |
---|
[c13ff56] | 1358 | { |
---|
[35f23d] | 1359 | temp = Table(string(nr[1]), "i", 1, i); |
---|
| 1360 | decP = Partitions(n - i*nr[1], delete(nr, 1)); |
---|
| 1361 | k = size(parts); |
---|
[c13ff56] | 1362 | for(j = 1; j <= size(decP); j++) |
---|
| 1363 | { |
---|
[35f23d] | 1364 | newP = temp + decP[j]; // new partition |
---|
[9ea76a] | 1365 | if(!containedQ(parts, newP, 1)) |
---|
[0610f0e] | 1366 | { |
---|
[35f23d] | 1367 | k = k + 1; |
---|
| 1368 | parts[k] = newP; |
---|
| 1369 | } |
---|
| 1370 | } |
---|
| 1371 | } |
---|
| 1372 | } |
---|
| 1373 | return(parts); |
---|
| 1374 | } |
---|
| 1375 | |
---|
| 1376 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1377 | |
---|
[70597c] | 1378 | static proc NumFactor(int a, int b) |
---|
[35f23d] | 1379 | " USAGE: NumFactor(a, b); int a, b |
---|
| 1380 | PURPOSE: if b divides a then return b/a, else return 0 |
---|
| 1381 | RETURN: int |
---|
| 1382 | " |
---|
| 1383 | { |
---|
[1f9a84] | 1384 | int c = b div a; |
---|
[35f23d] | 1385 | if(c*a == b) { return(c); } |
---|
| 1386 | else {return(0)} |
---|
| 1387 | } |
---|
| 1388 | |
---|
| 1389 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1390 | |
---|
[70597c] | 1391 | static proc Table(string cmd, string iterator, int lb, int ub) |
---|
[35f23d] | 1392 | " USAGE: Table(cmd,i, lb, ub); string cmd, i; int lb, ub |
---|
| 1393 | PURPOSE: generate a list of size ub - lb + 1 s.t. _[i] = cmd(i) |
---|
| 1394 | RETURN: list |
---|
| 1395 | " |
---|
| 1396 | { |
---|
| 1397 | list data; |
---|
| 1398 | execute("int " + iterator + ";"); |
---|
| 1399 | |
---|
[c13ff56] | 1400 | for(int @i = lb; @i <= ub; @i++) |
---|
| 1401 | { |
---|
[35f23d] | 1402 | execute(iterator + " = " + string(@i)); |
---|
| 1403 | execute("data[" + string(@i) + "] = " + cmd + ";"); |
---|
| 1404 | } |
---|
| 1405 | return(data); |
---|
| 1406 | } |
---|
| 1407 | |
---|
| 1408 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1409 | |
---|
| 1410 | static proc FlattenQHM(list data) |
---|
| 1411 | " USAGE: FlattenQHM(n, nr); list data |
---|
| 1412 | PURPOSE: flatten the list (one level) 'data', which is a list of lists |
---|
| 1413 | RETURN: list |
---|
| 1414 | " |
---|
| 1415 | { |
---|
| 1416 | int i, j, c; |
---|
| 1417 | list fList, temp; |
---|
| 1418 | |
---|
| 1419 | c = 1; |
---|
| 1420 | |
---|
[c13ff56] | 1421 | for(i = 1; i <= size(data); i++) |
---|
| 1422 | { |
---|
| 1423 | for(j = 1; j <= size(data[i]); j++) |
---|
| 1424 | { |
---|
[35f23d] | 1425 | fList[c] = data[i][j]; |
---|
| 1426 | c = c + 1; |
---|
| 1427 | } |
---|
| 1428 | } |
---|
| 1429 | return(fList); |
---|
| 1430 | } |
---|
| 1431 | |
---|
| 1432 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1433 | |
---|
| 1434 | static proc IntersectionQHM(list l1, list l2) |
---|
| 1435 | // Type : list |
---|
| 1436 | // Purpose : Intersection of l1 and l2 |
---|
| 1437 | { |
---|
| 1438 | list l; |
---|
| 1439 | int b, c; |
---|
| 1440 | |
---|
| 1441 | c = 1; |
---|
| 1442 | |
---|
[c13ff56] | 1443 | for(int i = 1; i <= size(l1); i++) |
---|
| 1444 | { |
---|
[9ea76a] | 1445 | b = containedQ(l2, l1[i]); |
---|
[c13ff56] | 1446 | if(b == 1) |
---|
| 1447 | { |
---|
[35f23d] | 1448 | l[c] = l1[i]; |
---|
[c13ff56] | 1449 | c++; |
---|
[35f23d] | 1450 | } |
---|
| 1451 | } |
---|
| 1452 | return(l); |
---|
| 1453 | } |
---|
| 1454 | |
---|
| 1455 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1456 | |
---|
| 1457 | static proc FirstEntryQHM(data, elem) |
---|
| 1458 | // Type : int |
---|
| 1459 | // Purpose : position of first entry equal to elem in data (from left to right) |
---|
| 1460 | { |
---|
| 1461 | int i, pos; |
---|
| 1462 | |
---|
| 1463 | i = 0; |
---|
| 1464 | pos = 0; |
---|
[c13ff56] | 1465 | while(i < size(data)) |
---|
| 1466 | { |
---|
| 1467 | i++; |
---|
| 1468 | if(data[i] == elem) { pos = i; break;} |
---|
[35f23d] | 1469 | } |
---|
| 1470 | return(pos); |
---|
| 1471 | } |
---|
| 1472 | |
---|
| 1473 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1474 | |
---|
| 1475 | static proc PSum(e) |
---|
| 1476 | { |
---|
| 1477 | poly f; |
---|
[c13ff56] | 1478 | for(int i = size(e);i>=1;i--) |
---|
| 1479 | { |
---|
[35f23d] | 1480 | f = f + e[i]; |
---|
| 1481 | } |
---|
| 1482 | return(f); |
---|
| 1483 | } |
---|
| 1484 | |
---|
| 1485 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1486 | |
---|
| 1487 | proc Max(data) |
---|
[fd5013] | 1488 | "USAGE: Max(data); intvec/list of integers |
---|
[35f23d] | 1489 | PURPOSE: find the maximal integer contained in 'data' |
---|
| 1490 | RETURN: list |
---|
| 1491 | ASSUME: 'data' contians only integers and is not empty |
---|
| 1492 | " |
---|
| 1493 | { |
---|
| 1494 | int i; |
---|
| 1495 | int max = data[1]; |
---|
| 1496 | |
---|
[c13ff56] | 1497 | for(i = size(data); i>1;i--) |
---|
| 1498 | { |
---|
[35f23d] | 1499 | if(data[i] > max) { max = data[i]; } |
---|
| 1500 | } |
---|
| 1501 | return(max); |
---|
| 1502 | } |
---|
[fd5013] | 1503 | example |
---|
| 1504 | {"EXAMPLE:"; echo = 2; |
---|
| 1505 | Max(list(1,2,3)); |
---|
| 1506 | } |
---|
[35f23d] | 1507 | |
---|
| 1508 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1509 | |
---|
| 1510 | proc Min(data) |
---|
[fd5013] | 1511 | "USAGE: Min(data); intvec/list of integers |
---|
[35f23d] | 1512 | PURPOSE: find the minimal integer contained in 'data' |
---|
| 1513 | RETURN: list |
---|
| 1514 | ASSUME: 'data' contians only integers and is not empty |
---|
| 1515 | " |
---|
| 1516 | { |
---|
| 1517 | int i; |
---|
| 1518 | int min = data[1]; |
---|
| 1519 | |
---|
[c13ff56] | 1520 | for(i = size(data);i>1; i--) |
---|
| 1521 | { |
---|
[35f23d] | 1522 | if(data[i] < min) { min = data[i]; } |
---|
| 1523 | } |
---|
| 1524 | return(min); |
---|
| 1525 | } |
---|
[fd5013] | 1526 | example |
---|
| 1527 | {"EXAMPLE:"; echo = 2; |
---|
| 1528 | Min(intvec(1,2,3)); |
---|
| 1529 | } |
---|
[35f23d] | 1530 | |
---|
| 1531 | /////////////////////////////////////////////////////////////////////////////// |
---|