[35f23d] | 1 | /////////////////////////////////////////////////////////////////////////////// |
---|
[e7778a] | 2 | version="$Id: qhmoduli.lib,v 1.14 2009-04-08 14:08:25 Singular Exp $"; |
---|
[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; } |
---|
| 57 | if(opt / 4 > 0) { imageQ = 1; opt = opt - 4;} |
---|
| 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; |
---|
| 132 | poly mPoly = 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); |
---|
| 145 | minpoly = number(imap(RME4, mPoly)); |
---|
| 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 | |
---|
| 217 | if(opt / 2 > 0) { primaryDec = 1; opt = opt - 2; } |
---|
| 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 ! |
---|
[9e9ec3] | 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; |
---|
| 331 | poly mPoly = minpoly; |
---|
| 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); |
---|
| 353 | ideal upperBasis, stabaction, action, mPoly, reduceIdeal; |
---|
| 354 | poly f, F, monos, h; |
---|
| 355 | |
---|
| 356 | reduceIdeal = imap(STR1, mPoly), imap(STR1, stabid); |
---|
| 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); |
---|
| 382 | minpoly = number(imap(STR1, mPoly)); |
---|
| 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. |
---|
[1f92589] | 474 | ASSUME: I semiquasihomogeneous ideal wrt '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 | poly mPoly; |
---|
| 485 | string ringSTR, ringSTR1, varString, parString; |
---|
| 486 | |
---|
| 487 | dbPrt = printlevel-voice+2; |
---|
| 488 | dbprint(dbPrt, "StabilizerEquations of " + string(data)); |
---|
| 489 | |
---|
| 490 | export(varSubsList); |
---|
| 491 | n = nvars(basering); |
---|
[afd77b2] | 492 | Variables = StabVar(wt); // possible quasihomogeneous substitutions |
---|
[35f23d] | 493 | nrVars = 0; |
---|
| 494 | for(i = 1; i <= size(wt); i = i + 1) { |
---|
[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) + ")"; |
---|
| 501 | if(npars(basering) == 1) { |
---|
| 502 | parString = "(0, " + parstr(basering) + ")"; |
---|
| 503 | } |
---|
| 504 | else { parString = "0"; } |
---|
| 505 | |
---|
| 506 | def RSTB = basering; |
---|
| 507 | mPoly = minpoly; |
---|
| 508 | ringSTR = "ring RSTR = " + parString + ", (" + varstr(basering) + ", " + varString + "), dp;"; // dp |
---|
| 509 | ringSTR1 = "ring RSTT = " + parString + ", (" + varString + ", " + varstr(basering) + "), dp;"; |
---|
| 510 | |
---|
[3b77465] | 511 | if(defined(RSTR)) { kill RSTR;} |
---|
| 512 | if(defined(RSTT)) { kill RSTT;} |
---|
[35f23d] | 513 | execute(ringSTR1); // this ring is only used for the result, where the variables |
---|
| 514 | export(RSTT); // are s(1..m),t(1..n), as needed for Derksens algorithm (NullCone) |
---|
| 515 | minpoly = number(imap(RSTB, mPoly)); |
---|
| 516 | |
---|
| 517 | execute(ringSTR); |
---|
| 518 | export(RSTR); |
---|
| 519 | minpoly = number(imap(RSTB, mPoly)); |
---|
| 520 | poly f, f1, g, h, vars, pp; // f1 is the polynomial after subs, |
---|
| 521 | ideal allEqns, qhsubs, actionid, stabid, J; |
---|
| 522 | list ringList; // all t(i)`s which do not appear in f1 |
---|
| 523 | ideal data = simplify(imap(RSTB, data), 2); |
---|
| 524 | |
---|
[1f92589] | 525 | // generate the quasihomogeneous substitution map F |
---|
[35f23d] | 526 | |
---|
| 527 | nrVars = 0; |
---|
| 528 | offset = 0; |
---|
| 529 | for(i = 1; i <= size(wt); i = i + 1) { // build the substitution t(i) -> ... |
---|
[afd77b2] | 530 | if(i > 1) { offset = offset + size(Variables[i - 1]); } |
---|
[35f23d] | 531 | g = 0; |
---|
[afd77b2] | 532 | for(j = 1; j <= size(Variables[i]); j = j + 1) { |
---|
[35f23d] | 533 | pp = 1; |
---|
[afd77b2] | 534 | for(k = 2; k <= size(Variables[i][j]); k = k + 1) { |
---|
| 535 | pp = pp * var(Variables[i][j][k]); |
---|
| 536 | if(Variables[i][j][k] == i) { varSubsList[i] = offset + j;} |
---|
[35f23d] | 537 | } |
---|
| 538 | g = g + s(offset + j) * pp; |
---|
| 539 | } |
---|
| 540 | qhsubs[i] = g; |
---|
| 541 | } |
---|
| 542 | dbprint(dbPrt, " qhasihomogenous substituion =" + string(qhsubs)); |
---|
| 543 | map F = RSTR, qhsubs; |
---|
[3b77465] | 544 | kill varSubsList; |
---|
[35f23d] | 545 | |
---|
| 546 | // get the equations of the stabilizer by comparing coefficients |
---|
| 547 | // in the equation f = F(f). |
---|
| 548 | |
---|
| 549 | vars = RingVarProduct(Table("i", "i", 1, size(wt))); |
---|
| 550 | |
---|
| 551 | allEqns = 0; |
---|
| 552 | |
---|
[e7778a] | 553 | matrix newcoMx, coMx; |
---|
| 554 | int d; |
---|
| 555 | for(r = 1; r <= ncols(data); r++) |
---|
| 556 | { |
---|
[35f23d] | 557 | |
---|
| 558 | f = data[r]; |
---|
| 559 | f1 = F(f); |
---|
[7de8e4] | 560 | d = deg(f); |
---|
| 561 | newcoMx = coef(f1, vars); // coefficients of F(f) |
---|
| 562 | coMx = coef(f, vars); // coefficients of f |
---|
[35f23d] | 563 | |
---|
| 564 | for(i = 1; i <= ncols(newcoMx); i = i + 1) { // build the system of eqns via coeff. comp. |
---|
| 565 | j = 1; |
---|
| 566 | h = 0; |
---|
| 567 | while(j <= ncols(coMx)) { // all monomials in f |
---|
| 568 | if(coMx[j][1] == newcoMx[i][1]) { h = coMx[j][2]; j = ncols(coMx) + 1;} |
---|
| 569 | else {j = j + 1;} |
---|
| 570 | } |
---|
| 571 | J = J, newcoMx[i][2] - h; // add equation |
---|
| 572 | } |
---|
| 573 | allEqns = allEqns, J; |
---|
| 574 | |
---|
| 575 | } |
---|
| 576 | allEqns = std(allEqns); |
---|
| 577 | |
---|
| 578 | // simplify the equations, i.e., if s(i) in J then remove s(i) from J |
---|
| 579 | // and from the basering |
---|
| 580 | |
---|
| 581 | sList = SimplifyIdeal(allEqns, n, "s"); |
---|
| 582 | stabid = sList[1]; |
---|
| 583 | map phi = basering, sList[2]; |
---|
| 584 | sln = size(sList[3]) - n; |
---|
| 585 | |
---|
| 586 | // change the substitution |
---|
| 587 | |
---|
| 588 | actionid = phi(qhsubs); |
---|
| 589 | |
---|
| 590 | // change to new ring, auxillary construction |
---|
| 591 | |
---|
| 592 | setring(RSTT); |
---|
| 593 | ideal actionid, stabid; |
---|
| 594 | |
---|
| 595 | actionid = imap(RSTR, actionid); |
---|
| 596 | stabid = imap(RSTR, stabid); |
---|
| 597 | export(stabid); |
---|
| 598 | export(actionid); |
---|
| 599 | ringList[2] = RSTT; |
---|
| 600 | |
---|
| 601 | dbprint(dbPrt, " substitution = " + string(actionid)); |
---|
| 602 | dbprint(dbPrt, " equations of stabilizer = " + string(stabid)); |
---|
| 603 | |
---|
| 604 | varString = "s(1.." + string(sln) + ")"; |
---|
| 605 | ringSTR = "ring RSTS = " + parString + ", (" + varString + "), dp;"; |
---|
| 606 | execute(ringSTR); |
---|
| 607 | minpoly = number(imap(RSTB, mPoly)); |
---|
| 608 | ideal stabid = std(imap(RSTR, stabid)); |
---|
| 609 | export(stabid); |
---|
| 610 | ringList[1] = RSTS; |
---|
| 611 | dbprint(dbPrt, " |
---|
| 612 | // 'StabEqnId' created a list of 2 rings. |
---|
| 613 | // To see the rings, type (if the name of your list is stab): |
---|
| 614 | show(stab); |
---|
| 615 | // To access the 1-st ring and map (and similair for the others), type: |
---|
| 616 | def S1 = stab[1]; setring S1; stabid; |
---|
| 617 | // S1/stabid is the coordinate ring of the variety of the |
---|
| 618 | // stabilizer, say G. If G x K^n --> K^n is the action of G on |
---|
| 619 | // K^n, then the ideal 'actionid' in the second ring describes |
---|
| 620 | // the dual map on the ring level. |
---|
| 621 | // To access the 2-nd ring and map (and similair for the others), type: |
---|
| 622 | def S2 = stab[2]; setring S2; actionid; |
---|
| 623 | "); |
---|
| 624 | return(ringList); |
---|
| 625 | } |
---|
| 626 | example |
---|
| 627 | {"EXAMPLE:"; echo = 2; |
---|
| 628 | ring B = 0,(x,y,z), ls; |
---|
| 629 | ideal I = x2,y3,z6; |
---|
| 630 | intvec w = 3,2,1; |
---|
| 631 | list stab = StabEqnId(I, w); |
---|
| 632 | def S1 = stab[1]; setring S1; stabid; |
---|
| 633 | def S2 = stab[2]; setring S2; actionid; |
---|
| 634 | } |
---|
| 635 | |
---|
| 636 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 637 | static |
---|
| 638 | proc ArnoldFormMain(poly f, B, poly Fs, ideal reduceIdeal, int nrs, int nrt) |
---|
| 639 | "USAGE: ArnoldFormMain(f, B, Fs, rI, nrs, nrt); |
---|
| 640 | poly f,Fs; ideal B, rI; int nrs, nrt |
---|
| 641 | PURPOSE: compute the induced action of 'G_f' on T_, where f is the principal |
---|
| 642 | part and 'Fs' is the semiuniversal unfolding of 'f' with x_i |
---|
| 643 | substituted by actionid[i], 'B' is a list of upper basis monomials |
---|
| 644 | for the milnor algebra of 'f', 'nrs' = number of variables for 'G_f' |
---|
| 645 | and 'nrt' = dimension of T_ |
---|
[1f92589] | 646 | ASSUME: f is quasihomogeneous with an isolated singularity at 0, |
---|
[35f23d] | 647 | s(1..r), t(1..m) are parameters of the basering |
---|
| 648 | RETURN: poly |
---|
| 649 | EXAMPLE: example ArnoldAction; shows an example |
---|
| 650 | " |
---|
| 651 | { |
---|
| 652 | int i, j, d, ub, dbPrt; |
---|
| 653 | list upperBasis, basisDegList, gmonos, common, parts; |
---|
| 654 | ideal jacobianId, jacobIdstd, mapId; // needed for phi |
---|
| 655 | intvec wt = weight(f); |
---|
| 656 | matrix gCoeffMx; // for lift command |
---|
| 657 | poly newFs, g, gred, tt; // g = sum of all monomials of degree d, gred is needed for lift |
---|
| 658 | map phi; // the map from Arnold's Theorem |
---|
| 659 | |
---|
| 660 | dbPrt = printlevel-voice+2; |
---|
| 661 | jacobianId = jacob(f); |
---|
| 662 | jacobIdstd = std(jacobianId); |
---|
| 663 | newFs = Fs; |
---|
| 664 | for(i = 1; i <= size(B); i = i + 1) { |
---|
| 665 | basisDegList[i] = deg(B[i], wt); |
---|
| 666 | } |
---|
| 667 | ub = Max(basisDegList) + 1; // max degree of an upper monomial |
---|
| 668 | |
---|
| 669 | parts = MonosAndTerms(newFs - f, wt, ub); |
---|
| 670 | gmonos = parts[1]; |
---|
| 671 | d = deg(f, wt); |
---|
| 672 | |
---|
| 673 | for(i = d + 1; i < ub; i = i + 1) { // base[1] = monomials of degree i |
---|
| 674 | upperBasis[i] = SelectMonos(list(B, B), wt, i); // B must not contain 0's |
---|
| 675 | } |
---|
| 676 | |
---|
| 677 | // test if each monomial of Fs is contained in B, if not, |
---|
| 678 | // compute a substitution via Arnold's theorem and substitutite |
---|
| 679 | // it into newFs |
---|
| 680 | |
---|
| 681 | for(i = d + 1; i < ub; i = i + 1) { // ub instead of @UB |
---|
| 682 | dbprint(dbPrt, "-- degree = " + string(i) + " of " + string(ub - 1) + " ---------------------------"); |
---|
| 683 | if(size(newFs) < 80) { dbprint(dbPrt, " polynomial = " + string(newFs - f));} |
---|
| 684 | else { dbprint(dbPrt, " poly has deg (not weighted) " + string(deg(newFs)) + " and contains " + string(size(newFs)) + " monos");} |
---|
| 685 | |
---|
| 686 | // select monomials of degree i and intersect them with upperBasis[i] |
---|
| 687 | |
---|
| 688 | gmonos = SelectMonos(parts, wt, i); |
---|
| 689 | common = IntersectionQHM(upperBasis[i][1], gmonos[1]); |
---|
| 690 | if(size(common) == size(gmonos[1])) { |
---|
| 691 | dbprint(dbPrt, " no additional monomials "); |
---|
| 692 | } |
---|
| 693 | |
---|
| 694 | // other monomials than those in upperBasis occur, compute |
---|
| 695 | // the map constructed in the proof of Arnold's theorem |
---|
| 696 | // write g = c[i] * jacobianId[i] |
---|
| 697 | |
---|
| 698 | else { |
---|
| 699 | dbprint(dbPrt, " additional Monomials found, compute the map "); |
---|
| 700 | g = PSum(gmonos[2]); // sum of all monomials in g of degree i |
---|
| 701 | dbprint(dbPrt, " sum of degree " + string(i) + " is " + string(g)); |
---|
| 702 | |
---|
| 703 | gred = reduce(g, jacobIdstd); |
---|
| 704 | gCoeffMx = lift(jacobianId, g - gred); // compute c[i] |
---|
| 705 | mapId = var(1) - gCoeffMx[1][1]; // generate the map |
---|
| 706 | for(j = 2; j <= size(gCoeffMx); j = j + 1) { |
---|
| 707 | mapId[j] = var(j) - gCoeffMx[1][j]; |
---|
| 708 | } |
---|
| 709 | dbprint(dbPrt, " map = " + string(mapId)); |
---|
| 710 | // apply the map to newFs |
---|
| 711 | newFs = APSubstitution(newFs, mapId, reduceIdeal, wt, ub, nrs, nrt); |
---|
| 712 | parts = MonosAndTerms(newFs - f, wt, ub); // monos and terms of deg < ub |
---|
| 713 | newFs = PSum(parts[2]) + f; // result of APS... is already reduced |
---|
| 714 | dbprint(dbPrt, " monomials of degree " + string(i)); |
---|
| 715 | } |
---|
| 716 | } |
---|
| 717 | return(newFs); |
---|
| 718 | } |
---|
| 719 | |
---|
| 720 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 721 | |
---|
[70597c] | 722 | static proc MonosAndTerms(poly f, wt, int ub) |
---|
[35f23d] | 723 | "USAGE: MonosAndTerms(f, w, ub); poly f, intvec w, int ub |
---|
| 724 | PURPOSE: returns a list of all monomials and terms occuring in f of |
---|
| 725 | weighted degree < ub |
---|
| 726 | RETURN: list |
---|
| 727 | _[1] list of monomials |
---|
[8942a5] | 728 | _[2] list of terms |
---|
[35f23d] | 729 | EXAMPLE: example MonosAndTerms shows an example |
---|
| 730 | " |
---|
| 731 | { |
---|
| 732 | int i, k; |
---|
| 733 | list monomials, terms; |
---|
| 734 | poly mono, lcInv, data; |
---|
| 735 | |
---|
| 736 | data = jet(f, ub - 1, wt); |
---|
| 737 | k = 0; |
---|
| 738 | for(i = 1; i <= size(data); i = i + 1) { |
---|
| 739 | mono = lead(data[i]); |
---|
| 740 | if(deg(mono, wt) < ub) { |
---|
| 741 | k = k + 1; |
---|
| 742 | lcInv = 1/leadcoef(mono); |
---|
| 743 | monomials[k] = mono * lcInv; |
---|
| 744 | terms[k] = mono; |
---|
| 745 | } |
---|
| 746 | } |
---|
| 747 | return(list(monomials, terms)); |
---|
| 748 | } |
---|
| 749 | example |
---|
| 750 | {"EXAMPLE:"; echo = 2; |
---|
| 751 | ring B = 0,(x,y,z), lp; |
---|
| 752 | poly f = 6*x2 + 2*x3 + 9*x*y2 + z*y + x*z6; |
---|
| 753 | MonosAndTerms(f, intvec(2,1,1), 5); |
---|
| 754 | } |
---|
| 755 | |
---|
| 756 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 757 | |
---|
[70597c] | 758 | static proc SelectMonos(parts, intvec wt, int d) |
---|
[35f23d] | 759 | "USAGE: SelectMonos(parts, w, d); list/ideal parts, intvec w, int d |
---|
| 760 | PURPOSE: returns a list of all monomials and terms occuring in f of |
---|
| 761 | weighted degree = d |
---|
| 762 | RETURN: list |
---|
[8942a5] | 763 | _[1] list of monomials |
---|
| 764 | _[2] list of terms |
---|
[35f23d] | 765 | EXAMPLE: example SelectMonos; shows an example |
---|
| 766 | " |
---|
| 767 | { |
---|
| 768 | int i, k; |
---|
| 769 | list monomials, terms; |
---|
| 770 | poly mono; |
---|
| 771 | |
---|
| 772 | k = 0; |
---|
| 773 | for(i = 1; i <= size(parts[1]); i = i + 1) { |
---|
| 774 | mono = parts[1][i]; |
---|
| 775 | if(deg(mono, wt) == d) { |
---|
| 776 | k = k + 1; |
---|
| 777 | monomials[k] = mono; |
---|
| 778 | terms[k] = parts[2][i]; |
---|
| 779 | } |
---|
| 780 | } |
---|
| 781 | return(list(monomials, terms)); |
---|
| 782 | } |
---|
| 783 | example |
---|
| 784 | {"EXAMPLE:"; echo = 2; |
---|
| 785 | ring B = 0,(x,y,z), lp; |
---|
| 786 | poly f = 6*x2 + 2*x3 + 9*x*y2 + z*y + x*z6; |
---|
| 787 | list mt = MonosAndTerms(f, intvec(2,1,1), 5); |
---|
| 788 | SelectMonos(mt, intvec(2,1,1), 4); |
---|
| 789 | } |
---|
| 790 | |
---|
| 791 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 792 | |
---|
[70597c] | 793 | static proc Expand(substitution, degVec, ideal reduceI, intvec w1, int ub, list truncated) |
---|
[35f23d] | 794 | "USAGE: Expand(substitution, degVec, reduceI, w, ub, truncated); |
---|
| 795 | ideal/list substitution, list/intvec degVec, ideal reduceI, intvec w, |
---|
| 796 | int ub, list truncated |
---|
| 797 | PURPOSE: substitute 'substitution' in the monomial given by the list of |
---|
| 798 | exponents 'degVec', omit all terms of weighted degree > ub and reduce |
---|
| 799 | the result w.r.t. 'reduceI'. If truncated[i] = 0 then the result is |
---|
| 800 | stored for later use. |
---|
| 801 | RETURN: poly |
---|
| 802 | NOTE: used by APSubstitution |
---|
| 803 | GLOBAL: computedPowers |
---|
| 804 | " |
---|
| 805 | { |
---|
| 806 | int i, minDeg; |
---|
| 807 | list powerList; |
---|
| 808 | poly g, h; |
---|
| 809 | |
---|
| 810 | // compute substitution[1]^degVec[1],...,subs[n]^degVec[n] |
---|
| 811 | |
---|
| 812 | for(i = 1; i <= ncols(substitution); i = i + 1) { |
---|
| 813 | if(size(substitution[i]) < 3 || degVec[i] < 4) { |
---|
| 814 | powerList[i] = reduce(substitution[i]^degVec[i], reduceI); // new |
---|
| 815 | } // directly for small exponents |
---|
| 816 | else { |
---|
| 817 | powerList[i] = PolyPower1(i, substitution[i], degVec[i], reduceI, w1, truncated[i], ub); |
---|
| 818 | } |
---|
| 819 | } |
---|
| 820 | // multiply the terms obtained by using PolyProduct(); |
---|
| 821 | g = powerList[1]; |
---|
| 822 | minDeg = w1[1] * degVec[1]; |
---|
| 823 | for(i = 2; i <= ncols(substitution); i = i + 1) { |
---|
| 824 | g = jet(g, ub - w1[i] * degVec[i] - 1, w1); |
---|
| 825 | h = jet(powerList[i], ub - minDeg - 1, w1); |
---|
| 826 | g = PolyProduct(g, h, reduceI, w1, ub); |
---|
| 827 | if(g == 0) { Print(" g = 0 "); break;} |
---|
| 828 | minDeg = minDeg + w1[i] * degVec[i]; |
---|
| 829 | } |
---|
| 830 | return(g); |
---|
| 831 | } |
---|
| 832 | |
---|
| 833 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 834 | |
---|
[70597c] | 835 | static proc PolyProduct(poly g1, poly h1, ideal reduceI, intvec wt, int ub) |
---|
[35f23d] | 836 | "USAGE: PolyProduct(g, h, reduceI, wt, ub); poly g, h; ideal reduceI, |
---|
| 837 | intvec wt, int ub. |
---|
| 838 | PURPOSE: compute g*h and reduce it w.r.t 'reduceI' and omit terms of weighted |
---|
| 839 | degree > ub. |
---|
| 840 | RETURN: poly |
---|
| 841 | NOTE: used by 'Expand' |
---|
| 842 | " |
---|
| 843 | { |
---|
| 844 | int SUBSMAXSIZE = 3000; // |
---|
| 845 | int i, nrParts, sizeOfPart, currentPos, partSize, maxSIZE; |
---|
| 846 | poly g, h, gxh, prodComp, @g2; // replace @g2 by subst. |
---|
| 847 | |
---|
| 848 | g = g1; |
---|
| 849 | h = h1; |
---|
| 850 | |
---|
| 851 | if(size(g)*size(h) > SUBSMAXSIZE) { |
---|
| 852 | |
---|
| 853 | // divide the polynomials with more terms in parts s.t. |
---|
| 854 | // the product of each part with the other polynomial |
---|
| 855 | // has at most SUBMAXSIZE terms |
---|
| 856 | |
---|
| 857 | if(size(g) < size(h)) { poly @h = h; h = g; g = @h;@h = 0; } |
---|
| 858 | maxSIZE = SUBSMAXSIZE / size(h); |
---|
| 859 | //print(" SUBSMAXSIZE = "+string(SUBSMAXSIZE)+" exceeded by "+string(size(g)*size(h)) + ", maxSIZE = ", string(maxSIZE)); |
---|
| 860 | nrParts = size(g) / maxSIZE + 1; |
---|
| 861 | partSize = size(g) / nrParts; |
---|
| 862 | gxh = 0; // 'g times h' |
---|
| 863 | for(i = 1; i <= nrParts; i = i + 1){ |
---|
| 864 | //print(" loop #" + string(i) + " of " + string(nrParts)); |
---|
| 865 | currentPos = (i - 1) * partSize; |
---|
| 866 | if(i < nrParts) {sizeOfPart = partSize;} |
---|
| 867 | else { sizeOfPart = size(g) - (nrParts - 1) * partSize; print(" last #" + string(sizeOfPart) + " terms ");} |
---|
| 868 | prodComp = g[currentPos + 1..sizeOfPart + currentPos] * h; // multiply a part |
---|
| 869 | @g2 = jet(prodComp, ub - 1, wt); // eventual reduce ... |
---|
| 870 | if(size(@g2) < size(prodComp)) { print(" killed " + string(size(prodComp) - size(@g2)) + " terms ");} |
---|
| 871 | gxh = reduce(gxh + @g2, reduceI); |
---|
| 872 | |
---|
| 873 | } |
---|
| 874 | } |
---|
| 875 | else { |
---|
| 876 | gxh = reduce(jet(g * h,ub - 1, wt), reduceI); |
---|
| 877 | } // compute directly |
---|
| 878 | return(gxh); |
---|
| 879 | } |
---|
| 880 | |
---|
| 881 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 882 | |
---|
[3c4dcc] | 883 | static proc PolyPower1(int varIndex, poly f, int e, ideal reduceI, intvec wt, |
---|
[70597c] | 884 | int truncated, int ub) |
---|
[35f23d] | 885 | "USAGE: PolyPower1(i, f, e, reduceI, wt, truncated, ub);int i, e, ub;poly f; |
---|
| 886 | ideal reduceI; intvec wt; list truncated; |
---|
| 887 | PURPOSE: compute f^e, use previous computations if possible, and reduce it |
---|
| 888 | w.r.t reudecI and omit terms of weighted degree > ub. |
---|
| 889 | RETURN: poly |
---|
| 890 | NOTE: used by 'Expand' |
---|
| 891 | GLOBAL: 'computedPowers' |
---|
| 892 | " |
---|
| 893 | { |
---|
| 894 | int i, ordOfg, lb, maxPrecomputedPower; |
---|
| 895 | poly g, fn; |
---|
| 896 | |
---|
| 897 | if(e == 0) { return(1);} |
---|
| 898 | if(e == 1) { return(f);} |
---|
| 899 | if(f == 0) { return(1); } |
---|
| 900 | else { |
---|
| 901 | |
---|
| 902 | // test if f has been computed to some power |
---|
| 903 | |
---|
| 904 | if(computedPowers[varIndex][1] > 0) { |
---|
| 905 | maxPrecomputedPower = computedPowers[varIndex][1]; |
---|
| 906 | if(maxPrecomputedPower >= e) { |
---|
| 907 | // no computation necessary, f^e has already benn computed |
---|
| 908 | g = computedPowers[varIndex][2][e - 1]; |
---|
| 909 | //Print("No computation, from list : g = elem [", varIndex, ", 2, ", e - 1, "]"); |
---|
| 910 | lb = e + 1; |
---|
| 911 | } |
---|
| 912 | else { // f^d computed, where d < e |
---|
| 913 | g = computedPowers[varIndex][2][maxPrecomputedPower - 1]; |
---|
| 914 | ordOfg = maxPrecomputedPower * wt[varIndex]; |
---|
| 915 | lb = maxPrecomputedPower + 1; |
---|
| 916 | } |
---|
| 917 | } |
---|
| 918 | else { // no precomputed data |
---|
| 919 | lb = 2; |
---|
| 920 | ordOfg = wt[varIndex]; |
---|
| 921 | g = f; |
---|
| 922 | } |
---|
| 923 | for(i = lb; i <= e; i = i + 1) { |
---|
| 924 | fn = jet(f, ub - ordOfg - 1, wt); // reduce w.r.t. reduceI |
---|
| 925 | g = PolyProduct(g, fn, reduceI, wt, ub); |
---|
| 926 | ordOfg = ordOfg + wt[varIndex]; |
---|
| 927 | if(g == 0) { break; } |
---|
| 928 | if((i > maxPrecomputedPower) && !truncated) { |
---|
| 929 | if(maxPrecomputedPower == 0) { // init computedPowers |
---|
| 930 | computedPowers[varIndex] = list(i, list(g)); |
---|
| 931 | } |
---|
| 932 | computedPowers[varIndex][1] = i; // new degree |
---|
| 933 | computedPowers[varIndex][2][i - 1] = g; |
---|
| 934 | maxPrecomputedPower = i; |
---|
| 935 | } |
---|
| 936 | } |
---|
| 937 | } |
---|
| 938 | return(g); |
---|
| 939 | } |
---|
| 940 | |
---|
| 941 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 942 | |
---|
[70597c] | 943 | static proc RingVarsToList(list @index) |
---|
[35f23d] | 944 | { |
---|
| 945 | int i; |
---|
| 946 | list temp; |
---|
| 947 | |
---|
| 948 | for(i = 1; i <= size(@index); i = i + 1) { temp[i] = string(var(@index[i])); } |
---|
| 949 | return(temp); |
---|
| 950 | } |
---|
| 951 | |
---|
| 952 | /////////////////////////////////////////////////////////////////////////////// |
---|
[70597c] | 953 | static |
---|
[35f23d] | 954 | proc APSubstitution(poly f, ideal substitution, ideal reduceIdeal, intvec wt, int ub, int nrs, int nrt) |
---|
| 955 | "USAGE: APSubstitution(f, subs, reduceI, w, ub, int nrs, int nrt); poly f |
---|
| 956 | ideal subs, reduceI, intvec w, int ub, nrs, nrt; |
---|
| 957 | nrs = number of parameters s(1..nrs), |
---|
| 958 | nrt = number of parameters t(1..nrt) |
---|
| 959 | PURPOSE: substitute 'subs' in f, omit all terms with weighted degree > ub and |
---|
| 960 | reduce the result w.r.t. 'reduceI'. |
---|
| 961 | RETURN: poly |
---|
| 962 | GLOBAL: 'computedPowers' |
---|
| 963 | " |
---|
| 964 | { |
---|
| 965 | int i, j, k, d, offset; |
---|
| 966 | int n = nvars(basering); |
---|
| 967 | list coeffList, parts, degVecList, degOfMonos; |
---|
| 968 | list computedPowers, truncatedQ, degOfSubs, @temp; |
---|
| 969 | string ringSTR, @ringVars; |
---|
| 970 | |
---|
| 971 | export(computedPowers); |
---|
| 972 | |
---|
| 973 | // store arguments in strings |
---|
| 974 | |
---|
| 975 | def RASB = basering; |
---|
| 976 | |
---|
| 977 | parts = MonosAndTerms(f, wt, ub); |
---|
| 978 | for(i = 1; i <= size(parts[1]); i = i + 1) { |
---|
| 979 | coeffList[i] = parts[2][i]/parts[1][i]; |
---|
| 980 | degVecList[i] = leadexp(parts[1][i]); |
---|
| 981 | degOfMonos[i] = deg(parts[1][i], wt); |
---|
| 982 | } |
---|
| 983 | |
---|
| 984 | // built new basering with no parameters and order dp ! |
---|
| 985 | // the parameters of the basering are appended to |
---|
| 986 | // the variables of the basering ! |
---|
| 987 | // set ideal mpoly = minpoly for reduction ! |
---|
| 988 | |
---|
| 989 | @ringVars = "(" + varstr(basering) + ", " + parstr(1) + ","; // precondition |
---|
| 990 | if(nrs > 0) { |
---|
| 991 | @ringVars = @ringVars + "s(1.." + string(nrs) + "), "; |
---|
| 992 | } |
---|
| 993 | @ringVars = @ringVars + "t(1.." + string(nrt) + "))"; |
---|
| 994 | ringSTR = "ring RASR = 0, " + @ringVars + ", dp;"; // new basering |
---|
| 995 | |
---|
| 996 | // built the "reduction" ring with the reduction ideal |
---|
| 997 | |
---|
| 998 | execute(ringSTR); |
---|
| 999 | export(RASR); |
---|
| 1000 | ideal reduceIdeal, substitution, newSubs; |
---|
| 1001 | intvec w1, degVec; |
---|
| 1002 | list minDeg, coeffList, degList; |
---|
| 1003 | poly f, g, h, subsPoly; |
---|
| 1004 | |
---|
| 1005 | w1 = wt; // new weights |
---|
| 1006 | offset = nrs + nrt + 1; |
---|
| 1007 | for(i = n + 1; i <= offset + n; i = i + 1) { w1[i] = 0; } |
---|
| 1008 | |
---|
| 1009 | reduceIdeal = std(imap(RASB, reduceIdeal)); // omit later ! |
---|
| 1010 | coeffList = imap(RASB, coeffList); |
---|
| 1011 | substitution = imap(RASB, substitution); |
---|
| 1012 | |
---|
| 1013 | f = imap(RASB, f); |
---|
| 1014 | |
---|
| 1015 | for(i = 1; i <= n; i = i + 1) { // all "base" variables |
---|
| 1016 | computedPowers[i] = list(0); |
---|
| 1017 | for(j = 1; j <= size(substitution[i]); j = j + 1) { degList[j] = deg(substitution[i][j], w1);} |
---|
| 1018 | degOfSubs[i] = degList; |
---|
| 1019 | } |
---|
| 1020 | |
---|
| 1021 | // substitute in each monomial seperately |
---|
| 1022 | |
---|
| 1023 | g = 0; |
---|
| 1024 | for(i = 1; i <= size(degVecList); i = i + 1) { |
---|
| 1025 | truncatedQ = Table("0", "i", 1, n); |
---|
| 1026 | newSubs = 0; |
---|
| 1027 | degVec = degVecList[i]; |
---|
| 1028 | d = degOfMonos[i]; |
---|
| 1029 | |
---|
| 1030 | // check if some terms in the substitution can be omitted |
---|
| 1031 | // degVec = list of exponents of the monomial m |
---|
| 1032 | // minDeg[j] denotes the weighted degree of the monomial m' |
---|
| 1033 | // where m' is the monomial m without the j-th variable |
---|
| 1034 | |
---|
| 1035 | for(j = 1; j <= size(degVec); j = j + 1) { minDeg[j] = d - degVec[j] * wt[j]; } |
---|
| 1036 | |
---|
| 1037 | for(j = 1; j <= size(degVec); j = j + 1) { |
---|
| 1038 | subsPoly = 0; // set substitution to 0 |
---|
| 1039 | if(degVec[j] > 0) { |
---|
| 1040 | |
---|
| 1041 | // if variable occurs then check if |
---|
| 1042 | // substitution[j][k] * (linear part)^(degVec[j]-1) + minDeg[j] < ub |
---|
| 1043 | // i.e. look for the smallest possible combination in subs[j]^degVec[j] |
---|
| 1044 | // which comes from the term substitution[j][k]. This term is multiplied |
---|
| 1045 | // with the rest of the monomial, which has at least degree minDeg[j]. |
---|
| 1046 | // If the degree of this product is < ub then substitution[j][k] contributes |
---|
| 1047 | // to the result and cannot be omitted |
---|
| 1048 | |
---|
| 1049 | for(k = 1; k <= size(substitution[j]); k = k + 1) { |
---|
| 1050 | if(degOfSubs[j][k] + (degVec[j] - 1) * wt[j] + minDeg[j] < ub) { |
---|
| 1051 | subsPoly = subsPoly + substitution[j][k]; |
---|
| 1052 | } |
---|
| 1053 | } |
---|
| 1054 | } |
---|
| 1055 | newSubs[j] = subsPoly; // set substitution |
---|
| 1056 | if(substitution[j] - subsPoly != 0) { truncatedQ[j] = 1;} // mark that substitution[j] is truncated |
---|
| 1057 | } |
---|
| 1058 | h = Expand(newSubs, degVec, reduceIdeal, w1, ub, truncatedQ) * coeffList[i]; // already reduced |
---|
| 1059 | g = reduce(g + h, reduceIdeal); |
---|
| 1060 | } |
---|
[3b77465] | 1061 | kill computedPowers; |
---|
[35f23d] | 1062 | setring RASB; |
---|
| 1063 | poly fnew = imap(RASR, g); |
---|
[3b77465] | 1064 | kill RASR; |
---|
[35f23d] | 1065 | return(fnew); |
---|
| 1066 | } |
---|
| 1067 | |
---|
| 1068 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1069 | |
---|
| 1070 | static proc StabVar(intvec wt) |
---|
| 1071 | "USAGE: StabVar(w); intvec w |
---|
[1f92589] | 1072 | PURPOSE: compute the indicies for quasihomogeneous substitutions of each |
---|
[35f23d] | 1073 | variable. |
---|
[1f92589] | 1074 | ASSUME: f semiquasihomogeneous polynomial with an isolated singularity at 0 |
---|
[35f23d] | 1075 | RETURN: list |
---|
| 1076 | _[i] list of combinations for var(i) (i must be appended |
---|
| 1077 | to each comb) |
---|
| 1078 | GLOBAL: 'varSubsList', contains the index j s.t. x(i) -> x(i)t(j) ... |
---|
| 1079 | " |
---|
| 1080 | { |
---|
| 1081 | int i, j, k, uw, ic; |
---|
[afd77b2] | 1082 | list varList, Variables, subs; |
---|
[35f23d] | 1083 | string str, varString; |
---|
| 1084 | |
---|
| 1085 | varList = StabVarComb(wt); |
---|
| 1086 | for(i = 1; i <= size(wt); i = i + 1) { |
---|
| 1087 | subs = 0; |
---|
| 1088 | |
---|
| 1089 | // built linear substituitons |
---|
| 1090 | for(j = 1; j <= size(varList[1][i]); j = j + 1) { |
---|
| 1091 | subs[j] = list(i) + list(varList[1][i][j]); |
---|
| 1092 | } |
---|
[afd77b2] | 1093 | Variables[i] = subs; |
---|
[35f23d] | 1094 | if(size(varList[2][i]) > 0) { |
---|
| 1095 | |
---|
| 1096 | // built nonlinear substituitons |
---|
| 1097 | subs = 0; |
---|
| 1098 | for(j = 1; j <= size(varList[2][i]); j = j + 1) { |
---|
| 1099 | subs[j] = list(i) + varList[2][i][j]; |
---|
| 1100 | } |
---|
[afd77b2] | 1101 | Variables[i] = Variables[i] + subs; |
---|
[35f23d] | 1102 | } |
---|
| 1103 | } |
---|
[afd77b2] | 1104 | return(Variables); |
---|
[35f23d] | 1105 | } |
---|
| 1106 | |
---|
| 1107 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1108 | |
---|
| 1109 | static proc StabVarComb(intvec wt) |
---|
| 1110 | "USAGE: StabVarComb(w); intvec w |
---|
| 1111 | PURPOSE: list all possible indices of indeterminates for a quasihom. subs. |
---|
| 1112 | RETURN: list |
---|
| 1113 | _[1] linear substitutions |
---|
| 1114 | _[2] nonlinear substiutions |
---|
| 1115 | GLOBAL: 'varSubsList', contains the index j s.t. x(i) -> x(i)t(j) ... |
---|
| 1116 | " |
---|
| 1117 | { |
---|
| 1118 | int mi, ma, i, j, k, uw, ic; |
---|
| 1119 | list index, indices, usedWeights, combList, combinations; |
---|
| 1120 | list linearSubs, nonlinearSubs; |
---|
| 1121 | list partitions, subs, temp; // subs[i] = substitution for var(i) |
---|
| 1122 | |
---|
| 1123 | linearSubs = Table("0", "i", 1, size(wt)); |
---|
| 1124 | nonlinearSubs = Table("0", "i", 1, size(wt)); |
---|
| 1125 | |
---|
| 1126 | uw = 0; |
---|
| 1127 | ic = 0; |
---|
| 1128 | mi = Min(wt); |
---|
| 1129 | ma = Max(wt); |
---|
| 1130 | |
---|
| 1131 | for(i = mi; i <= ma; i = i + 1) { |
---|
| 1132 | if(ContainedQ(wt, i)) { // find variables of weight i |
---|
| 1133 | k = 0; |
---|
| 1134 | index = 0; |
---|
| 1135 | // collect the indices of all variables of weight i |
---|
| 1136 | for(j = 1; j <= size(wt); j = j + 1) { |
---|
| 1137 | if(wt[j] == i) { |
---|
| 1138 | k = k + 1; |
---|
| 1139 | index[k] = j; |
---|
| 1140 | } |
---|
| 1141 | } |
---|
| 1142 | uw = uw + 1; |
---|
| 1143 | usedWeights[uw] = i; |
---|
| 1144 | ic = ic + 1; |
---|
| 1145 | indices[i] = index; |
---|
| 1146 | |
---|
| 1147 | // linear part of the substitution |
---|
| 1148 | |
---|
| 1149 | for(j = 1; j <= size(index); j = j + 1) { |
---|
| 1150 | linearSubs[index[j]] = index; |
---|
| 1151 | } |
---|
| 1152 | |
---|
| 1153 | // nonlinear part of the substitution |
---|
| 1154 | |
---|
| 1155 | if(uw > 1) { // variables of least weight do not allow nonlinear subs. |
---|
| 1156 | |
---|
| 1157 | partitions = Partitions(i, delete(usedWeights, uw)); |
---|
| 1158 | for(j = 1; j <= size(partitions); j = j + 1) { |
---|
| 1159 | combinations[j] = AllCombinations(partitions[j], indices); |
---|
| 1160 | } |
---|
| 1161 | for(j = 1; j <= size(index); j = j + 1) { |
---|
| 1162 | nonlinearSubs[index[j]] = FlattenQHM(combinations); // flatten one level ! |
---|
| 1163 | } |
---|
| 1164 | |
---|
| 1165 | } // end if |
---|
| 1166 | |
---|
| 1167 | } |
---|
| 1168 | } |
---|
| 1169 | combList[1] = linearSubs; |
---|
| 1170 | combList[2] = nonlinearSubs; |
---|
| 1171 | return(combList); |
---|
| 1172 | } |
---|
| 1173 | |
---|
| 1174 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1175 | |
---|
[70597c] | 1176 | static proc AllCombinations(list partition, list indices) |
---|
[35f23d] | 1177 | "USAGE: AllCombinations(partition,indices); list partition, indices) |
---|
[8942a5] | 1178 | PURPOSE: all combinations for a given partititon |
---|
[35f23d] | 1179 | RETURN: list |
---|
| 1180 | GLOBAL: varSubsList, contains the index j s.t. x(i) -> x(i)t(j) ... |
---|
| 1181 | " |
---|
| 1182 | { |
---|
| 1183 | int i, k, m, ok, p, offset; |
---|
| 1184 | list nrList, indexList; |
---|
| 1185 | |
---|
| 1186 | k = 0; |
---|
| 1187 | offset = 0; |
---|
| 1188 | i = 1; |
---|
| 1189 | ok = 1; |
---|
| 1190 | m = partition[1]; |
---|
| 1191 | while(ok) { |
---|
| 1192 | if(i > size(partition)) { |
---|
| 1193 | ok = 0; |
---|
| 1194 | p = 0; |
---|
| 1195 | } |
---|
| 1196 | else { p = partition[i];} |
---|
| 1197 | if(p == m) { i = i + 1;} |
---|
| 1198 | else { |
---|
| 1199 | k = k + 1; |
---|
| 1200 | nrList[k] = i - 1 - offset; |
---|
| 1201 | offset = offset + i - 1; |
---|
| 1202 | indexList[k] = indices[m]; |
---|
| 1203 | if(ok) { m = partition[i];} |
---|
| 1204 | } |
---|
| 1205 | } |
---|
| 1206 | return(AllCombinationsAux(nrList, indexList)); |
---|
| 1207 | } |
---|
| 1208 | |
---|
| 1209 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1210 | |
---|
[70597c] | 1211 | static proc AllSingleCombinations(int n, list index) |
---|
[35f23d] | 1212 | "USAGE: AllSingleCombinations(n index); int n, list index |
---|
| 1213 | PURPOSE: all combinations for var(n) |
---|
| 1214 | RETURN: list |
---|
| 1215 | " |
---|
| 1216 | { |
---|
| 1217 | int i, j, k; |
---|
| 1218 | list comb, newC, temp, newIndex; |
---|
| 1219 | |
---|
| 1220 | if(n == 1) { |
---|
| 1221 | for(i = 1; i <= size(index); i = i + 1) { |
---|
| 1222 | temp = index[i]; |
---|
| 1223 | comb[i] = temp; |
---|
| 1224 | } |
---|
| 1225 | return(comb); |
---|
| 1226 | } |
---|
| 1227 | if(size(index) == 1) { |
---|
| 1228 | temp = Table(string(index[1]), "i", 1, n); |
---|
| 1229 | comb[1] = temp; |
---|
| 1230 | return(comb); |
---|
| 1231 | } |
---|
| 1232 | newIndex = index; |
---|
| 1233 | for(i = 1; i <= size(index); i = i + 1) { |
---|
| 1234 | if(i > 1) { newIndex = delete(newIndex, 1); } |
---|
| 1235 | newC = AllSingleCombinations(n - 1, newIndex); |
---|
| 1236 | k = size(comb); |
---|
| 1237 | temp = 0; |
---|
| 1238 | for(j = 1; j <= size(newC); j = j + 1) { |
---|
| 1239 | temp[1] = index[i]; |
---|
| 1240 | temp = temp + newC[j]; |
---|
| 1241 | comb[k + j] = temp; |
---|
| 1242 | temp = 0; |
---|
| 1243 | } |
---|
| 1244 | } |
---|
| 1245 | return(comb); |
---|
| 1246 | } |
---|
| 1247 | |
---|
| 1248 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1249 | |
---|
[70597c] | 1250 | static proc AllCombinationsAux(list parts, list index) |
---|
[35f23d] | 1251 | "USAGE: AllCombinationsAux(parts ,index); list parts, index |
---|
| 1252 | PURPOSE: all compbinations for nonlinear substituiton |
---|
| 1253 | RETURN: list |
---|
| 1254 | " |
---|
| 1255 | { |
---|
| 1256 | int i, j, k; |
---|
| 1257 | list comb, firstC, restC; |
---|
| 1258 | |
---|
| 1259 | if(size(parts) == 0 || size(index) == 0) { return(comb);} |
---|
| 1260 | |
---|
| 1261 | firstC = AllSingleCombinations(parts[1], index[1]); |
---|
| 1262 | restC = AllCombinationsAux(delete(parts, 1), delete(index, 1)); |
---|
| 1263 | |
---|
| 1264 | if(size(restC) == 0) { comb = firstC;} |
---|
| 1265 | else { |
---|
| 1266 | for(i = 1; i <= size(firstC); i = i + 1) { |
---|
| 1267 | k = size(comb); |
---|
| 1268 | for(j = 1; j <= size(restC); j = j + 1) { |
---|
| 1269 | //elem = firstC[i] + restC[j]; |
---|
| 1270 | // comb[k + j] = elem; |
---|
| 1271 | comb[k + j] = firstC[i] + restC[j]; |
---|
| 1272 | } |
---|
| 1273 | } |
---|
| 1274 | } |
---|
| 1275 | return(comb); |
---|
| 1276 | } |
---|
| 1277 | |
---|
| 1278 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1279 | |
---|
[70597c] | 1280 | static proc Partitions(int n, list nr) |
---|
[35f23d] | 1281 | "USAGE: Partitions(n, nr); int n, list nr |
---|
| 1282 | PURPOSE: partitions of n consisting of elements from nr |
---|
| 1283 | RETURN: list |
---|
| 1284 | " |
---|
| 1285 | { |
---|
| 1286 | int i, j, k; |
---|
| 1287 | list parts, temp, restP, newP, decP; |
---|
| 1288 | |
---|
| 1289 | if(size(nr) == 0) { return(list());} |
---|
| 1290 | if(size(nr) == 1) { |
---|
| 1291 | if(NumFactor(nr[1], n) > 0) { |
---|
| 1292 | parts[1] = Table(string(nr[1]), "i", 1, NumFactor(nr[1], n)); |
---|
| 1293 | } |
---|
| 1294 | return(parts); |
---|
| 1295 | } |
---|
| 1296 | else { |
---|
| 1297 | parts = Partitions(n, nr[1]); |
---|
| 1298 | restP = Partitions(n, delete(nr, 1)); |
---|
| 1299 | |
---|
| 1300 | parts = parts + restP; |
---|
| 1301 | for(i = 1; i <= n / nr[1]; i = i + 1) { |
---|
| 1302 | temp = Table(string(nr[1]), "i", 1, i); |
---|
| 1303 | decP = Partitions(n - i*nr[1], delete(nr, 1)); |
---|
| 1304 | k = size(parts); |
---|
| 1305 | for(j = 1; j <= size(decP); j = j + 1) { |
---|
| 1306 | newP = temp + decP[j]; // new partition |
---|
| 1307 | if(!ContainedQ(parts, newP, 1)) { |
---|
| 1308 | k = k + 1; |
---|
| 1309 | parts[k] = newP; |
---|
| 1310 | } |
---|
| 1311 | } |
---|
| 1312 | } |
---|
| 1313 | } |
---|
| 1314 | return(parts); |
---|
| 1315 | } |
---|
| 1316 | |
---|
| 1317 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1318 | |
---|
[70597c] | 1319 | static proc NumFactor(int a, int b) |
---|
[35f23d] | 1320 | " USAGE: NumFactor(a, b); int a, b |
---|
| 1321 | PURPOSE: if b divides a then return b/a, else return 0 |
---|
| 1322 | RETURN: int |
---|
| 1323 | " |
---|
| 1324 | { |
---|
| 1325 | int c = int(b/a); |
---|
| 1326 | if(c*a == b) { return(c); } |
---|
| 1327 | else {return(0)} |
---|
| 1328 | } |
---|
| 1329 | |
---|
| 1330 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1331 | |
---|
[70597c] | 1332 | static proc Table(string cmd, string iterator, int lb, int ub) |
---|
[35f23d] | 1333 | " USAGE: Table(cmd,i, lb, ub); string cmd, i; int lb, ub |
---|
| 1334 | PURPOSE: generate a list of size ub - lb + 1 s.t. _[i] = cmd(i) |
---|
| 1335 | RETURN: list |
---|
| 1336 | " |
---|
| 1337 | { |
---|
| 1338 | list data; |
---|
| 1339 | execute("int " + iterator + ";"); |
---|
| 1340 | |
---|
| 1341 | for(int @i = lb; @i <= ub; @i++) { |
---|
| 1342 | execute(iterator + " = " + string(@i)); |
---|
| 1343 | execute("data[" + string(@i) + "] = " + cmd + ";"); |
---|
| 1344 | } |
---|
| 1345 | return(data); |
---|
| 1346 | } |
---|
| 1347 | |
---|
| 1348 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1349 | |
---|
| 1350 | static proc FlattenQHM(list data) |
---|
| 1351 | " USAGE: FlattenQHM(n, nr); list data |
---|
| 1352 | PURPOSE: flatten the list (one level) 'data', which is a list of lists |
---|
| 1353 | RETURN: list |
---|
| 1354 | " |
---|
| 1355 | { |
---|
| 1356 | int i, j, c; |
---|
| 1357 | list fList, temp; |
---|
| 1358 | |
---|
| 1359 | c = 1; |
---|
| 1360 | |
---|
| 1361 | for(i = 1; i <= size(data); i = i + 1) { |
---|
| 1362 | for(j = 1; j <= size(data[i]); j = j + 1) { |
---|
| 1363 | fList[c] = data[i][j]; |
---|
| 1364 | c = c + 1; |
---|
| 1365 | } |
---|
| 1366 | } |
---|
| 1367 | return(fList); |
---|
| 1368 | } |
---|
| 1369 | |
---|
| 1370 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1371 | |
---|
| 1372 | static proc IntersectionQHM(list l1, list l2) |
---|
| 1373 | // Type : list |
---|
| 1374 | // Purpose : Intersection of l1 and l2 |
---|
| 1375 | { |
---|
| 1376 | list l; |
---|
| 1377 | int b, c; |
---|
| 1378 | |
---|
| 1379 | c = 1; |
---|
| 1380 | |
---|
| 1381 | for(int i = 1; i <= size(l1); i = i + 1) { |
---|
| 1382 | b = ContainedQ(l2, l1[i]); |
---|
| 1383 | if(b == 1) { |
---|
| 1384 | l[c] = l1[i]; |
---|
| 1385 | c = c + 1; |
---|
| 1386 | } |
---|
| 1387 | } |
---|
| 1388 | return(l); |
---|
| 1389 | } |
---|
| 1390 | |
---|
| 1391 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1392 | |
---|
| 1393 | static proc FirstEntryQHM(data, elem) |
---|
| 1394 | // Type : int |
---|
| 1395 | // Purpose : position of first entry equal to elem in data (from left to right) |
---|
| 1396 | { |
---|
| 1397 | int i, pos; |
---|
| 1398 | |
---|
| 1399 | i = 0; |
---|
| 1400 | pos = 0; |
---|
| 1401 | while(!pos && i < size(data)) { |
---|
| 1402 | i = i + 1; |
---|
| 1403 | if(data[i] == elem) { pos = i;} |
---|
| 1404 | } |
---|
| 1405 | return(pos); |
---|
| 1406 | } |
---|
| 1407 | |
---|
| 1408 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1409 | |
---|
| 1410 | static proc PSum(e) |
---|
| 1411 | { |
---|
| 1412 | poly f; |
---|
| 1413 | for(int i = 1; i <= size(e); i = i + 1) { |
---|
| 1414 | f = f + e[i]; |
---|
| 1415 | } |
---|
| 1416 | return(f); |
---|
| 1417 | |
---|
| 1418 | } |
---|
| 1419 | |
---|
| 1420 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1421 | |
---|
| 1422 | proc Max(data) |
---|
[fd5013] | 1423 | "USAGE: Max(data); intvec/list of integers |
---|
[35f23d] | 1424 | PURPOSE: find the maximal integer contained in 'data' |
---|
| 1425 | RETURN: list |
---|
| 1426 | ASSUME: 'data' contians only integers and is not empty |
---|
| 1427 | " |
---|
| 1428 | { |
---|
| 1429 | int i; |
---|
| 1430 | int max = data[1]; |
---|
| 1431 | |
---|
| 1432 | for(i = 1; i <= size(data); i = i + 1) { |
---|
| 1433 | if(data[i] > max) { max = data[i]; } |
---|
| 1434 | } |
---|
| 1435 | return(max); |
---|
| 1436 | } |
---|
[fd5013] | 1437 | example |
---|
| 1438 | {"EXAMPLE:"; echo = 2; |
---|
| 1439 | Max(list(1,2,3)); |
---|
| 1440 | } |
---|
[35f23d] | 1441 | |
---|
| 1442 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 1443 | |
---|
| 1444 | proc Min(data) |
---|
[fd5013] | 1445 | "USAGE: Min(data); intvec/list of integers |
---|
[35f23d] | 1446 | PURPOSE: find the minimal integer contained in 'data' |
---|
| 1447 | RETURN: list |
---|
| 1448 | ASSUME: 'data' contians only integers and is not empty |
---|
| 1449 | " |
---|
| 1450 | { |
---|
| 1451 | int i; |
---|
| 1452 | int min = data[1]; |
---|
| 1453 | |
---|
| 1454 | for(i = 1; i <= size(data); i = i + 1) { |
---|
| 1455 | if(data[i] < min) { min = data[i]; } |
---|
| 1456 | } |
---|
| 1457 | return(min); |
---|
| 1458 | } |
---|
[fd5013] | 1459 | example |
---|
| 1460 | {"EXAMPLE:"; echo = 2; |
---|
| 1461 | Min(intvec(1,2,3)); |
---|
| 1462 | } |
---|
[35f23d] | 1463 | |
---|
| 1464 | /////////////////////////////////////////////////////////////////////////////// |
---|