Annotation of OpenXM_contrib2/asir2018/lib/primdec_lex, Revision 1.1
1.1 ! noro 1: /*
! 2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
! 3: * All rights reserved.
! 4: *
! 5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
! 6: * non-exclusive and royalty-free license to use, copy, modify and
! 7: * redistribute, solely for non-commercial and non-profit purposes, the
! 8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
! 9: * conditions of this Agreement. For the avoidance of doubt, you acquire
! 10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
! 11: * third party developer retains all rights, including but not limited to
! 12: * copyrights, in and to the SOFTWARE.
! 13: *
! 14: * (1) FLL does not grant you a license in any way for commercial
! 15: * purposes. You may use the SOFTWARE only for non-commercial and
! 16: * non-profit purposes only, such as academic, research and internal
! 17: * business use.
! 18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
! 19: * international copyright treaties. If you make copies of the SOFTWARE,
! 20: * with or without modification, as permitted hereunder, you shall affix
! 21: * to all such copies of the SOFTWARE the above copyright notice.
! 22: * (3) An explicit reference to this SOFTWARE and its copyright owner
! 23: * shall be made on your publication or presentation in any form of the
! 24: * results obtained by use of the SOFTWARE.
! 25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
! 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
! 27: * for such modification or the source code of the modified part of the
! 28: * SOFTWARE.
! 29: *
! 30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
! 31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
! 32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
! 33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
! 34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
! 35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
! 36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
! 37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
! 38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
! 39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
! 40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
! 41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
! 42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
! 43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
! 44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
! 45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
! 46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
! 47: *
! 48: * $OpenXM$
! 49: */
! 50: /* Primary decomposition & Radical decomposition program */
! 51: /* written by T.Shimoyama, Fujitsu Lab. Date: 1995.10.12 */
! 52:
! 53: /* comments about flags *
! 54: PRIMEORD term order of Groebner basis in primedec
! 55: PRIMAORD term order of Groebner basis in primadec
! 56: PRINTLOG print out log (1 (simple) -- 5 (pricise))
! 57: SHOWTIME if 1 : print out timings and data
! 58: ISOLATED if 1 : compute only isolated primary components
! 59: NOEMBEDDED if 1 : compute only pseudo-primary components
! 60: NOINITGB if 1 : no initial G-base computation
! 61: REDUNDANT if 1 : no redundancy check
! 62: COMPLETE if 1 : use complete criterion of redundancy
! 63: COMMONCHECK if 1 : redundancy check by intersection (in gr_fctr_sub)
! 64: SELECTFLAG selection strategy of separators (0 -- 3)
! 65: */
! 66:
! 67: if (vtype(minipoly) != 3) load("gr")$$
! 68:
! 69: #define GR(R,F,V,O) T2=newvect(4,time());R=dp_gr_main(F,V,0,0,O);GRTIME+=newvect(4,time())-T2;
! 70: #define HGRM(R,F,V,O) T2=newvect(4,time());R=dp_gr_main(F,V,1,1,O);GRTIME+=newvect(4,time())-T2;
! 71: #define NF(R,IN,F,G,O) T2=newvect(4,time());R=dp_nf(IN,F,G,O);NFTIME+=newvect(4,time())-T2;
! 72:
! 73: /* optional flags */
! 74: extern PRIMAORD,PRIMEORD,PRINTLOG,SHOWTIME,NOINITGB$
! 75: extern COMMONCHECK,ISOLATED,NOEMBEDDED,REDUNDANT,SELECTFLAG,COMPLETE$
! 76: extern EXTENDED,CONTINUE,BSAVE,MAXALGIND,NOSPECIALDEC$
! 77: def primflags() {
! 78: print("PRIMAORD,PRIMEORD,PRINTLOG,SHOWTIME,NOINITGB,ISOLATED,NOEMBEDDED,COMMONCHECK");
! 79: print("REDUNDANT,SELECTFLAG,COMPLETE,EXTENDED,CONTINUE,BSAVE,MAXALGIND,NOSPECIALDEC");}
! 80: PRIMAORD=0$ PRIMEORD=2$
! 81:
! 82: /* global variables */
! 83: extern DIRECTRY,COMMONIDEAL,NOMRADDEC,NOMISOLATE,RADDECTIME$
! 84: extern MISTIME,ISORADDEC,GRTIME,NFTIME$
! 85:
! 86:
! 87: /* primary ideal decomposition of ideal(F) */
! 88: /* return the list of [P,Q] where Q is P-primary component of ideal(F) */
! 89: def primadec(F,VL)
! 90: {
! 91: if ( !VL ) {
! 92: print("invalid variables"); return 0; }
! 93: if ( !F ) {
! 94: print("invalid argument"); return 0; }
! 95: if ( F == [] )
! 96: return [];
! 97: if ( length(F) == 1 && type(F[0]) == 1 )
! 98: return [[1],[1]];
! 99: NOMRADDEC = NOMISOLATE = 0; RADDECTIME = newvect(4); T0 = newvect(4,time());
! 100: GRTIME = newvect(4); NFTIME = newvect(4); MISTIME = newvect(4);
! 101:
! 102: R = primadec_main(F,[["begin",F,F,[1],[]]],[],VL);
! 103:
! 104: if ( PRINTLOG ) {
! 105: print(""); print("primary decomposition done.");
! 106: }
! 107: if ( PRINTLOG || SHOWTIME ) {
! 108: print("number of radical decompositions = ",0); print(NOMRADDEC);
! 109: print("number of primary components = ",0); print(length(R),0);
! 110: print(" ( isolated = ",0); print(NOMISOLATE,0); print(" )");
! 111: print("total time of m.i.s. computation = ",0);
! 112: print(MISTIME[0],0); GT = MISTIME[1];
! 113: if ( GT ) { print(" + gc ",0);print(GT); } else print("");
! 114: print("total time of n-form computation = ",0);
! 115: print(NFTIME[0],0); GT = NFTIME[1];
! 116: if ( GT ) { print(" + gc ",0);print(GT); } else print("");
! 117: print("total time of G-base computation = ",0);
! 118: print(GRTIME[0],0); GT = GRTIME[1];
! 119: if ( GT ) { print(" + gc ",0);print(GT); } else print("");
! 120: print("total time of radical decomposition = ",0);
! 121: print(RADDECTIME[0],0); GT = RADDECTIME[1];
! 122: if ( GT ) { print(" + gc ",0);print(GT,0); }
! 123: print(" (iso ",0); print(ISORADDEC[0],0); GT = ISORADDEC[1];
! 124: if ( GT ) { print(" +gc ",0);print(GT,0); } print(")");
! 125: print("total time of primary decomposition = ",0);
! 126: TT = newvect(4,time())-T0; print(TT[0],0);
! 127: if ( TT[1] ) { print(" + gc ",0);print(TT[1]); } else print("");
! 128: }
! 129: return R;
! 130: }
! 131:
! 132: /* prime ideal decomposition of radical(F) */
! 133: /* return the list of prime components of radical(F) */
! 134:
! 135: def primedec(F,VL)
! 136: {
! 137: if ( !VL ) {
! 138: print("invalid variables"); return 0; }
! 139: if ( !F ) {
! 140: print("invalid argument"); return 0; }
! 141: if ( F == [] )
! 142: return [];
! 143: GRTIME = newvect(4);
! 144: T0 = newvect(4,time());
! 145: if ( !NOINITGB ) {
! 146: if ( PRINTLOG ) {
! 147: print("G-base computation w.r.t ordering = ",0);
! 148: print(PRIMEORD);
! 149: }
! 150: T1 = newvect(4,time());
! 151: HGRM(F,F,VL,PRIMEORD);
! 152: Tg = newvect(4,time())-T1;
! 153: if ( PRINTLOG > 0 ) {
! 154: print(" G-base time = ",0); print(Tg[0],0);
! 155: if ( Tg[1] ) { print(" + gc : ",0); print(Tg[1]); }
! 156: else print("");
! 157: }
! 158: }
! 159:
! 160: R = primedec_main([F],VL);
! 161:
! 162: Ta = newvect(4,time())-T0;
! 163: if ( PRINTLOG || SHOWTIME ) {
! 164: print("number of prime components = ",0); print(length(R));
! 165: print("G-base time = ",0); print(GRTIME[0],0);
! 166: if ( GRTIME[1] ) { print(" + gc : ",0); print(GRTIME[1]); }
! 167: else print("");
! 168: print("total time = ",0); print(Ta[0],0);
! 169: if ( Ta[1] ) { print(" + gc : ",0); print(Ta[1]); } else print("");
! 170: }
! 171: return R;
! 172: }
! 173:
! 174: /* main procedure for primary decomposition. */
! 175: /* F : ideal, VL : variable list, REMS : remaining components */
! 176: def primadec_main(F,REMS,H,VL)
! 177: {
! 178: DEC = RES = [];
! 179: for (D = [], I = 0; I < length(REMS); I++) {
! 180: MARK = REMS[I][0]; WF = REMS[I][1]; WP = REMS[I][2]; SC = REMS[I][3];
! 181: if ( !NOINITGB || MARK != "begin" ) {
! 182: ORD = PRIMAORD;
! 183: if ( PRINTLOG > 1 ) {
! 184: if ( MARK != "begin" ) print("");
! 185: print("G-base computation w.r.t ordering = ",0);
! 186: print(ORD);
! 187: }
! 188: T1 = newvect(4,time());
! 189: /* G-base of ideal */
! 190: GR(GF,WF,VL,ORD);
! 191: if ( MARK != "begin" && ( COMPLETE || EXTENDED ) ) {
! 192: if ( SC!=[1] && SC!=[-1] ) {
! 193: LG = localize(WF,SC,VL,ORD); /* VR_s\cap R */
! 194: if ( EXTENDED ) { GF = LG; }
! 195: } else
! 196: LG = GF;
! 197: if ( idealinc(H,LG,VL,ORD) ) {
! 198: if ( PRINTLOG ) {
! 199: print("eliminate the remaining component ",0);
! 200: print("by complete criterion");
! 201: }
! 202: continue; /* complete criterion */
! 203: }
! 204: }
! 205: /* G-base of radical */
! 206: if ( MARK == "begin" ) {
! 207: RA = ["begin",GF];
! 208: } else if ( MARK == "ext" ) {
! 209: if ( WF==WP || idealinc(WP,GF,VL,ORD) )
! 210: RA = ["ext",GF];
! 211: else {
! 212: if ( EXTENDED ) {
! 213: GA = localize(WP,SC,VL,PRIMEORD);
! 214: } else {
! 215: GR(GA,WP,VL,PRIMEORD);
! 216: }
! 217: RA = ["ext",GA];
! 218: }
! 219: } else if ( MARK == "sep" ) {
! 220: for (U=[], T=WP,J=length(T)-1;J>=0;J--) {
! 221: if ( EXTENDED ) {
! 222: GA = localize(T[J],SC,VL,PRIMEORD);
! 223: } else {
! 224: GR(GA,T[J],VL,PRIMEORD);
! 225: }
! 226: if (GA != [1] && GA != [-1])
! 227: U = cons(GA,U);
! 228: }
! 229: RA = ["sep",U];
! 230: } else debug;
! 231: Tg = newvect(4,time())-T1;
! 232: if ( PRINTLOG > 1 ) {
! 233: print(" G-base time = ",0); print(Tg[0],0);
! 234: if ( Tg[1] ) { print(" + gc : ",0); print(Tg[1]); }
! 235: else print("");
! 236: }
! 237: } else {
! 238: GF = F; /* NOINITGB = 1 */
! 239: RA = ["begin",GF];
! 240: }
! 241: if ( PRINTLOG ) {
! 242: if ( MARK == "begin" ) {
! 243: print("primary decomposition of the ideal");
! 244: } else { /* at the begining */
! 245: print("");
! 246: print("primary decomposition of the remaining component");
! 247: }
! 248: if ( MARK == "begin" && PRINTLOG > 1 ) { /* at the begining */
! 249: print(" ideal = ",0);
! 250: print(WF);
! 251: } else if ( PRINTLOG >= 4 ) {
! 252: print(" remaining component = ",0);
! 253: print(GF);
! 254: }
! 255: }
! 256:
! 257: /* input: init, generator, G-base, radical, intersection, sep.cond.*/
! 258: /* output: primary comp. remaining comp. */
! 259: Y = isolated(F,WF,GF,RA,REMS[I][4],SC,VL);
! 260:
! 261: D = append(D,Y[0]);
! 262: if ( MARK == "begin" )
! 263: NOMISOLATE = length(Y[0]);
! 264: RES = append(RES,Y[1]);
! 265: }
! 266: if ( MARK == "begin" ) {
! 267: F = GF; /* input polynomial -> G-base of it */
! 268: }
! 269: DEC = append(DEC,D);
! 270: if ( PRINTLOG ) {
! 271: print("");
! 272: print("# number of remainig components = ",0); print(length(RES));
! 273: }
! 274: if ( !length(RES) )
! 275: return DEC;
! 276: if ( !REDUNDANT ) { /* check whether Id(F,RES) is redundant or not */
! 277: for(L = [H],I = length(D)-1; I>=0; I--)
! 278: L=append(L,[D[I][0]]);
! 279: H = idealsetcap(L,VL,ORD); /* new intersection H */
! 280: if ( idealinc(H,F,VL,ORD) ) {
! 281: if ( PRINTLOG ) {
! 282: print("");
! 283: print("all primary components are found!");
! 284: }
! 285: return DEC;
! 286: }
! 287: REMS = mksepext(RES,H,VL); /* new remainig comps. */
! 288: } else {
! 289: REMS = mksepext(RES,H,VL); /* new rmaining comps. */
! 290: }
! 291: D = primadec_main(F,REMS,H,VL);
! 292: return append(DEC,D);
! 293: }
! 294:
! 295: /* isolated components and embedded components */
! 296: /* GF is the G-base of F, RA is the radical of F */
! 297: def isolated(IP,F,GF,RA,H,SC,VL)
! 298: {
! 299: T0 = newvect(4,time());
! 300: if ( RA[0] == "begin" ) {
! 301: PD = primedec_main([RA[1]],VL);
! 302: PD = map(dp_gr_main,PD,VL,0,1,PRIMEORD); /* XXX */
! 303: } else if ( RA[0] == "ext" || RA[0] == "sep" ) {
! 304: if ( RA[0] == "sep" )
! 305: T = prime_irred(idealsav(RA[1]),VL);
! 306: else
! 307: T = [RA[1]];
! 308: if ( !NOSPECIALDEC )
! 309: PD = primedec_special(T,VL);
! 310: else
! 311: PD = primedec_main(T,VL);
! 312: } else debug;
! 313: T1 = newvect(4,time())-T0;
! 314: if ( RA[0] == "begin" ) ISORADDEC = T1;
! 315: NOMRADDEC++; RADDECTIME += T1;
! 316: if ( PRINTLOG ) {
! 317: print("number of prime components = ",0); print(length(PD),0);
! 318: print(", time = ",0); print(T1[0],0);
! 319: if ( T1[1] ) { print(" + gc : ",0); print(T1[1]); } else print("");
! 320: if ( PRINTLOG > 0 ) {
! 321: print("dimensions of primes = ",0);
! 322: for (I = 0; I < length(PD); I++) {
! 323: print(length(VL)-length(minalgdep(PD[I],VL,PRIMEORD)),0);
! 324: print(", ",0);
! 325: }
! 326: print("");
! 327: }
! 328: }
! 329: if ( RA[0] == "begin" ) { /* isolated part of initial ideal */
! 330: if ( PRINTLOG ) {
! 331: print("check 'prime decomposition = primary decomposition?'");
! 332: }
! 333: CP = idealsetcap(PD,VL,PRIMAORD);
! 334: if ( idealinc(CP,GF,VL,PRIMAORD) ) {
! 335: if ( PRINTLOG ) {
! 336: print("lucky!");
! 337: }
! 338: for (L = [],I = length(PD)-1; I >= 0; I--)
! 339: L = cons([PD[I],PD[I]],L);
! 340: return [L,[]];
! 341: }
! 342: if ( PRINTLOG ) {
! 343: print("done.");
! 344: }
! 345: }
! 346:
! 347: R = pseudo_extract(IP,F,GF,PD,H,SC,VL);
! 348:
! 349: return R;
! 350: }
! 351:
! 352: def pseudo_extract(IP,F,GF,PD,H,SC,VL)
! 353: {
! 354: REMS = DEC = PDEC = SEL = RES = [];
! 355: ZERODIM = 1; CAP = H;
! 356: for (I = 0; I < length(PD); I++) {
! 357: P = PD[I];
! 358: if ( length(minalgdep(P,VL,PRIMEORD)) != length(VL) )
! 359: ZERODIM=0;
! 360: if ( length(PD) == 1 ) { /* pseudo primary ideal case */
! 361: if ( PRINTLOG >= 1 ) { print(""); print("pseudo primary ideal"); }
! 362: DD = GF; SEL = [SC];
! 363: } else {
! 364: T0 = time();
! 365: Y = pseudodec_main(F,P,PD,VL);
! 366: T1 = time();
! 367: DD = Y[0]; SS = Y[1]; SEL = append(SEL,[SS]);
! 368: PDEC = append(PDEC,[[DD,P]]);
! 369: if ( PRINTLOG >= 1 ) {
! 370: print(""); print("pseudo primary component of ",0);
! 371: print(I,0); print(", time = ",0); print(T1[0]-T0[0]);
! 372: if ( PRINTLOG >= 4 ) { print(" = ",0); print(DD); }
! 373: }
! 374: if ( NOEMBEDDED )
! 375: continue;
! 376: }
! 377: if ( !REDUNDANT && H != [] ) { /* elim. pseu-comp. */
! 378: if ( !sepcond_ps(P,SC,VL) )
! 379: continue;
! 380: LG = localize(DD,SC,VL,PRIMAORD);
! 381: if ( idealinc(H,LG,VL,PRIMAORD)) {
! 382: if ( PRINTLOG ) {
! 383: print("eliminate the pseudo primary component ",0);
! 384: print(I);
! 385: }
! 386: continue;
! 387: }
! 388: }
! 389: if ( PRINTLOG ) {
! 390: print("extraction of the pseudo primary component ",0);
! 391: print(I);
! 392: if ( PRINTLOG > 2 ) {
! 393: print(" associated prime of pseudo primary ideal = ",0);
! 394: print(P);
! 395: }
! 396: }
! 397: U = extraction(DD,P,VL);
! 398: if ( !REDUNDANT && H != [] && idealinc(H,U[0][0],VL,PRIMAORD) ) {
! 399: if ( PRINTLOG )
! 400: print("redundant primary component!");
! 401: } else {
! 402: DEC = append(DEC,[U[0]]);
! 403: H = idealcap(H,U[0][0],VL,PRIMAORD);
! 404: if ( idealeq(IP,H) ) {
! 405: if ( PRINTLOG ) {
! 406: print("");
! 407: print("all primary components are found!");
! 408: }
! 409: return [DEC,[]];
! 410: }
! 411: }
! 412: if ( !ISOLATED && U[1] != [] )
! 413: if ( sepcond_re(U[1],SC,VL) ) {
! 414: NSC = setunion([SS,SC]); /* new separating condition */
! 415: REM = cons(DD,append(U[1],[NSC,H]));
! 416: REMS = append(REMS,[REM]);
! 417: }
! 418: }
! 419: if ( NOEMBEDDED )
! 420: DEC = PDEC;
! 421: if ( length(PD) != 1 && !NOEMBEDDED && !ISOLATED && !ZERODIM ) {
! 422: for (QD=[],I=length(PDEC)-1;I>=0;I--)
! 423: QD = cons(PDEC[I][0],QD);
! 424: RES = ["sep",PD,QD,GF];
! 425: if ( sepcond_re(append(RES,[SEL]),SC,VL) ) {
! 426: REM = cons(F,append(RES,[SC,H]));
! 427: REMS = append(REMS,[REM]);
! 428: }
! 429: }
! 430: return [DEC,REMS];
! 431: }
! 432:
! 433: /* F:input set, PD:primes of radical, E:higher dimensional ideal */
! 434: /* length(PD) > 1 */
! 435: /* output : pseudo primary comp., remaining comp., selectors */
! 436: def pseudodec_main(F,P,PD,VL)
! 437: {
! 438: ZERODIM = 1;
! 439: S = selects(P,PD,VL,SELECTFLAG);
! 440: R = localize(F,S,VL,PRIMAORD);
! 441: if ( R[0] == 1 || R[0] == -1 ) {
! 442: print("error, R = ",0); print(R);
! 443: debug;
! 444: }
! 445: R = idealnormal(R);
! 446: return [R,S];
! 447: }
! 448:
! 449: /* Id(GF) is a pseudo primary ideal. (GF must be the G-base.) */
! 450: def extraction(GF,Pr,VL)
! 451: {
! 452: TMPORD1=TMPORD2=0;
! 453: V = minalgdep(Pr,VL,PRIMEORD);
! 454: U = listminus(VL,V);
! 455: V0 = append(V,U);
! 456: #if 0
! 457: /* This may be a bug. GF is not a GB w.r.t the elimination order */
! 458: if ( V0 != VL ) {
! 459: ORD = [[TMPORD1,length(V)],[TMPORD2,length(U)]];
! 460: GR(G,GF,V0,ORD);
! 461: } else
! 462: G = GF;
! 463: #else
! 464: ORD = [[TMPORD1,length(V)],[TMPORD2,length(U)]];
! 465: GR(G,GF,V0,ORD);
! 466: #endif
! 467: dp_ord(TMPORD1);
! 468: for (LL = [],HC = 1,I = 0; I < length(G); I++)
! 469: LL = append(LL,cdr(fctr(dp_hc(dp_ptod(G[I],V)))));
! 470: for (L=[],I=0;I<length(LL);I++)
! 471: L = cons(LL[I][0],L);
! 472: L = setunion([L]);
! 473: for (S=1,SL=[],I=0;I<length(L);I++) {
! 474: S *= L[I];
! 475: if ( SELECTFLAG )
! 476: SL = cons(L[I],SL);
! 477: }
! 478: if ( !SELECTFLAG )
! 479: SL= [S];
! 480: if ( PRINTLOG > 1 ) {
! 481: print("extractor = ",0);
! 482: print(S);
! 483: }
! 484: T0 = time()[0];
! 485: Q = localize(GF,SL,VL,PRIMAORD);
! 486: Q = idealnormal(Q);
! 487: DEC = [Q,Pr];
! 488: if ( PRINTLOG ) {
! 489: print("find a primary component! time = ",0);
! 490: print(time()[0]-T0);
! 491: if (PRINTLOG >= 3){
! 492: print(" associated prime of primary component = ",0);
! 493: print(DEC[1]);
! 494: }
! 495: if (PRINTLOG >= 4){print(" primary component = ",0); print(Q);}
! 496: }
! 497: if ( !ISOLATED && !NOEMBEDDED && SL != [1]
! 498: && length(V) != length(VL) /* nonzerodim */
! 499: && (REDUNDANT || !idealinc(Q,GF,VL,PRIMAORD)) ) {
! 500: REM = ["ext",[Q,Pr],GF,S];
! 501: if ( PRINTLOG ) {
! 502: print("find the remaining component of the extraction");
! 503: }
! 504: } else {
! 505: REM = [];
! 506: if ( PRINTLOG ) {
! 507: print("eliminate the remaining component of the extraction");
! 508: }
! 509: }
! 510: return [DEC,REM];
! 511: }
! 512:
! 513: /* sep. cond. for pseudo-primary components */
! 514: def sepcond_ps(P,SC,VL)
! 515: {
! 516: for (J = 0; J < length(SC); J++) {
! 517: if ( idealinc([SC[J]],P,VL,PRIMEORD) )
! 518: break; /* separating condition */
! 519: }
! 520: if ( J != length(SC) ) {
! 521: if ( PRINTLOG ) {
! 522: print("");
! 523: print("elim. the pseudo primary comp. by separating cond.");
! 524: }
! 525: return 0;
! 526: }
! 527: return 1;
! 528: }
! 529:
! 530: /* sep. cond. for rem. components. */
! 531: /* REM = ["ext",[Q,Pr],GF,S] or ["sep",PD,QD,GF,SEL], SC : sep.cond. */
! 532: def sepcond_re(REM,SC,VL)
! 533: {
! 534: for (S=1,I=0;I<length(SC);I++)
! 535: S *= SC[I];
! 536: if (REM[0] == "ext") {
! 537: F = cons(REM[3],REM[1][1]);
! 538: L = localize(F,[S],VL,PRIMAORD);
! 539: if ( L != [1] )
! 540: return 1;
! 541: else
! 542: return 0;
! 543: } else if (REM[0] == "sep") {
! 544: PL = REM[1]; SEL = REM[4];
! 545: for (I=0;I<length(PL);I++) {
! 546: F = append(SEL[I],PL[I]);
! 547: L = localize(F,[S],VL,PRIMAORD);
! 548: if ( L != [1] )
! 549: return 1;
! 550: }
! 551: return 0;
! 552: }
! 553: }
! 554:
! 555: def minalgdep(Pr,VL,ORD)
! 556: {
! 557: T0=newvect(4,time());
! 558: if (MAXALGIND)
! 559: R = minalgdep1(Pr,VL,ORD); /* M.I.S. */
! 560: else
! 561: R = minalgdep0(Pr,VL,ORD); /* M.S.I.S. */
! 562: T1=newvect(4,time())-T0; MISTIME += T1;
! 563: if ( PRINTLOG >= 6 ) {
! 564: if (MAXALGIND) print("M.I.S. time = ",0);
! 565: else print("M.S.I.S. time = ",0);
! 566: print(T1[0]);
! 567: }
! 568: return R;
! 569: }
! 570:
! 571: def minalgdep0(Pr,VL,ORD)
! 572: {
! 573: dp_ord(ORD); N = length(Pr);
! 574: for (HD = [],I = N-1; I >= 0; I--)
! 575: HD = cons(vars(dp_dtop(dp_ht(dp_ptod(Pr[I],VL)),VL)),HD);
! 576: for (R = X = [], I = length(VL)-1; I >= 0; I--) {
! 577: V = VL[I];
! 578: TX = cons(V,X);
! 579: for (J = 0; J < N; J++) {
! 580: if ( varinc(HD[J],TX) )
! 581: break;
! 582: }
! 583: if ( J == N )
! 584: X = TX;
! 585: else
! 586: R = cons(V,R);
! 587: }
! 588: return R;
! 589: }
! 590:
! 591: def minalgdep1(Pr,VL,ORD)
! 592: {
! 593: R = all_msis(Pr,VL,ORD);
! 594: for (E=I=0,D=length(R[0]);I<length(R);I++) {
! 595: if ( length(R[I])>=D ) {
! 596: D=length(R[I]); E=I;
! 597: }
! 598: }
! 599: S = listminus(VL,R[E]);
! 600: return S;
! 601: }
! 602:
! 603: def all_msis(Pr,VL,ORD)
! 604: {
! 605: dp_ord(ORD); N = length(Pr);
! 606: for (HD = [],I = N-1; I >= 0; I--)
! 607: HD = cons(vars(dp_dtop(dp_ht(dp_ptod(Pr[I],VL)),VL)),HD);
! 608: R = dimrec(HD,0,[],[],VL);
! 609: return R;
! 610: }
! 611:
! 612: def dimrec(S,K,U,M,X)
! 613: {
! 614: MM = M;
! 615: for (I=K;I<length(X);I++) {
! 616: TX = append(U,[X[I]]);
! 617: for (J = 0; J < length(S); J++) {
! 618: if ( varinc(S[J],TX) )
! 619: break;
! 620: }
! 621: if ( J == length(S) )
! 622: MM = dimrec(S,I+1,TX,MM,X);
! 623: }
! 624: for (J = 0; J < length(MM); J++) {
! 625: if ( varinc(U,MM[J]) )
! 626: break;
! 627: }
! 628: if ( J == length(MM) )
! 629: MM = append(MM,[U]);
! 630: return MM;
! 631: }
! 632:
! 633: /* if V is min.alg.dep. return 1 */
! 634: def ismad(Pr,V,VL,ORD)
! 635: {
! 636: dp_ord(ORD); N = length(Pr);
! 637: for (HD = [],I = N-1; I >= 0; I--)
! 638: HD = cons(vars(dp_dtop(dp_ht(dp_ptod(Pr[I],VL)),VL)),HD);
! 639: W = listminus(VL,V); L = append(W,V);
! 640: for (R = TX = [],I = 0; I < length(L); I++) {
! 641: TX = cons(L[I],TX);
! 642: for (J = 0; J < N; J++) {
! 643: if ( varinc(HD[J],TX) )
! 644: break;
! 645: }
! 646: if ((I<length(W)&&J!=N) || (I>=length(W)&&J==N))
! 647: return 0;
! 648:
! 649: }
! 650: return 1;
! 651: }
! 652:
! 653: /* output:list of [flag, generator, primes, sep.cond., intersection] */
! 654: def mksepext(RES,H,VL)
! 655: {
! 656: for (R = [],I=length(RES)-1;I>=0;I--) {
! 657: W = RES[I];
! 658: if ( COMPLETE || EXTENDED )
! 659: HI = idealcap(H,W[6],VL,PRIMAORD);
! 660: else
! 661: HI = H;
! 662: if (W[1] == "sep") {
! 663: E = mksep(W[2],W[3],W[4],VL);
! 664: F = append(E[0],W[0]);
! 665: if ( ( COMPLETE || EXTENDED ) && !REDUNDANT )
! 666: HI = idealsetcap(cons(HI,W[3]),VL,PRIMAORD);
! 667: R = cons(["sep",F,E[1],W[5],HI],R);
! 668: } else if (W[1] == "ext") {
! 669: E = mkext(W[2][0],W[3],W[4],VL);
! 670: F = cons(E[0],W[0]);
! 671: P = cons(E[1],W[2][1]); /* prime component of rem. comp. */
! 672: R = cons(["ext",F,P,W[5],HI],R);
! 673: } else debug;
! 674: }
! 675: return R;
! 676: }
! 677:
! 678: /* make extractor with exponents ; F must be G-Base */
! 679: /* output : extractor, max.sqfr.of extractor */
! 680: def mkext(Q,F,S,VL)
! 681: {
! 682: if ( PRINTLOG > 1 ) {
! 683: print("make extractor = ",0);
! 684: }
! 685: if ( S == 1 || S == -1 ) {
! 686: print(1); return [1,1];
! 687: }
! 688: DD=dp_mkdist([Q],F,VL,PRIMAORD);
! 689: DQ=DD[0][0]; DF=DD[1];
! 690: for (L=fctr(S),SL=[],J=length(L)-1;J>0;J--)
! 691: SL=cons(dp_ptod(L[J][0],VL),SL);
! 692: K = dp_fexponent(dp_ptod(S,VL),DF,DQ);
! 693: E = dp_vectexp(SL,EX=dp_minexp(K,SL,DF,DQ));
! 694: if ( E != 1 ) E = dp_dtop(E,VL);
! 695: for (U=1, I=0;I<size(EX)[0];I++)
! 696: if ( EX[I] )
! 697: U *= L[I+1][0];
! 698: /* for (L=sqfr(E),U=1,J=length(L)-1;J>0;J--)
! 699: U *= L[J][0]; */
! 700: if ( PRINTLOG > 1 ) {
! 701: print(E);
! 702: }
! 703: return [E,U];
! 704: }
! 705:
! 706: /* make separators with exponents ; F must be G-Base */
! 707: /* output [separator, list of components of radical] */
! 708: def mksep(PP,QQ,F,VL)
! 709: {
! 710: if ( PRINTLOG > 1 ) {
! 711: print("make separators, ",0); print(length(PP));
! 712: }
! 713: DD=dp_mkdist(QQ,F,VL,PRIMAORD);
! 714: DQ=DD[0]; DF=DD[1];
! 715: E = dp_mksep(PP,DQ,DF,VL,PRIMAORD);
! 716: for (RA=[],I=length(E)-1;I>=0;I--) {
! 717: for (L=fctr(E[I]),J=length(L)-1;J>0;J--) {
! 718: ES=cons(L[J][0],PP[I]);
! 719: RA = cons(ES,RA); /* radical of remaining comp. */
! 720: }
! 721: }
! 722: return [E,RA];
! 723: }
! 724:
! 725: def dp_mksep(PP,DQ,DF,VL,ORD)
! 726: {
! 727: dp_ord(ORD);
! 728: for (E=[],I=0;I<length(PP);I++) {
! 729: T0=time()[0];
! 730: if ( CONTINUE ) {
! 731: E = bload("sepa");
! 732: I = length(E);
! 733: CONTINUE = 0;
! 734: }
! 735: S = selects(PP[I],PP,VL,0)[0];
! 736: for (L=fctr(S),SL=[],J=length(L)-1;J>0;J--)
! 737: SL=cons(dp_ptod(L[J][0],VL),SL);
! 738: K = dp_fexponent(dp_ptod(S,VL),DF,DQ[I]);
! 739: M = dp_vectexp(SL,dp_minexp(K,SL,DF,DQ[I]));
! 740: if (M != 1) M = dp_dtop(M,VL);
! 741: E = append(E,[M]);
! 742: if ( PRINTLOG > 1 ) {
! 743: print("separator of ",0); print(I,0);
! 744: print(" in ",0); print(length(PP),0);
! 745: print(", time = ",0); print(time()[0]-T0,0); print(", ",0);
! 746: printsep(cdr(fctr(M))); print("");
! 747: }
! 748: if (BSAVE) {
! 749: bsave(E,"sepa");
! 750: print("bsave to sepa");
! 751: }
! 752: }
! 753: return E;
! 754: }
! 755:
! 756: extern AFO$
! 757:
! 758: /* F,Id,Q is given by distributed polynomials. Id is a G-Base.*/
! 759: /* get the exponent K s.t. (Id : F^K) = Q, Id is a G-Basis */
! 760: def dp_fexponent(F,Id,Q)
! 761: {
! 762: if (!AFO)
! 763: for (IN=[],I=size(Id)[0]-1;I>=0;I--) IN = cons(I,IN);
! 764: else
! 765: for (IN=[],I=0;I<size(Id)[0];I++) IN = cons(I,IN);
! 766: NF(G,IN,F,Id,0);
! 767: if ( F != G ) { F = G; }
! 768: N=size(Q)[0]; U = newvect(N);
! 769: for (I=0;I<N;I++)
! 770: U[I]=Q[I];
! 771: for (K = 1;; K++) { /* U[I] = Q[I]*F^K mod Id */
! 772: for (FLAG = I = 0; I < N; I++) {
! 773: NF(U[I],IN,U[I]*F,Id,1);
! 774: if ( U[I] ) FLAG = 1;
! 775: }
! 776: if ( !FLAG )
! 777: return K;
! 778: }
! 779: }
! 780:
! 781: /* F must be a G-Base. SL is a set of factors of S.*/
! 782: def dp_minexp(K,SL,F,Q)
! 783: {
! 784: if (!AFO)
! 785: for (IN=[],I=size(F)[0]-1;I>=0;I--) IN = cons(I,IN);
! 786: else
! 787: for (IN=[],I=0;I<size(F)[0];I++) IN = cons(I,IN);
! 788: N = length(SL); M = size(Q)[0];
! 789: EX = newvect(N); U = newvect(M); T = newvect(M);
! 790: for (I=0;I<N;I++)
! 791: EX[I]=K;
! 792: for (I=0;I<M;I++) {
! 793: NF(Q[I],IN,Q[I],F,1);
! 794: }
! 795: for (I=0;I<N;I++) {
! 796: EX[I] = 0;
! 797: for (FF=1,II=0;II<N;II++)
! 798: for (J = 0; J < EX[II]; J++) {
! 799: NF(FF,IN,FF*SL[II],F,1);
! 800: }
! 801: for (J = 0; J < M; J++)
! 802: U[J] = Q[J]*FF;
! 803: for (; EX[I] < K; EX[I]++) {
! 804: FLAG = 0;
! 805: for (J = 0; J < M; J++) {
! 806: NF(U[J],IN,U[J],F,1);
! 807: if ( U[J] ) FLAG = 1;
! 808: }
! 809: if ( !FLAG )
! 810: break;
! 811: if ( EX[I] < K-1 )
! 812: for (J = 0; J < M; J++) U[J] = U[J]*SL[I];
! 813: }
! 814: }
! 815: return EX;
! 816: }
! 817:
! 818: def dp_vectexp(SS,EE)
! 819: {
! 820: N = length(SS);
! 821: for (A=1,I=0;I<N;I++)
! 822: for (J = 0; J < EE[I]; J++)
! 823: A *= SS[I];
! 824: return A;
! 825: }
! 826:
! 827: def dp_mkdist(QQ,F,VL,ORD)
! 828: {
! 829: dp_ord(ORD);
! 830: for (DQ=[],I=length(QQ)-1;I>=0;I--) {
! 831: for (T=[],J=length(QQ[I])-1;J>=0;J--)
! 832: T = cons(dp_ptod(QQ[I][J],VL),T);
! 833: DQ=cons(newvect(length(T),T),DQ);
! 834: }
! 835: for (T=[],J=length(F)-1;J>=0;J--)
! 836: T = cons(dp_ptod(F[J],VL),T);
! 837: DF = newvect(length(T),T);
! 838: return [DQ,DF]$
! 839: }
! 840:
! 841: /* if FLAG == 0 return singltonset */
! 842: def selects(P,PD,VL,FLAG)
! 843: {
! 844: PP=dp_mkdist([],P,VL,PRIMEORD)[1];
! 845: for (M=[],I=length(P)-1;I>=0;I--)
! 846: M = cons(I,M);
! 847: for (SL = [],I = length(PD)-1; I >= 0; I--) {
! 848: B = PD[I];
! 849: if (B == P) continue;
! 850: for (L = [],J = 0; J < length(B); J++) {
! 851: A = B[J];
! 852: NF(E,M,dp_ptod(A,VL),PP,0);
! 853: if ( E ) {
! 854: L = append(L,[A]); /* A \notin Id(P) */
! 855: }
! 856: }
! 857: SL = cons(L,SL); /* candidate of separators */
! 858: }
! 859: for (S=[],I=length(SL)-1;I>=0;I--) {/* choose an element from L */
! 860: if ( FLAG >= 2 ) {
! 861: S = append(SL[I],S);
! 862: continue;
! 863: }
! 864: C = minselect(SL[I],VL,PRIMAORD);
! 865: S = cons(C,S);
! 866: }
! 867: S = setunion([S]);
! 868: if ( FLAG == 3 || length(S) == 1 )
! 869: return S;
! 870: if ( FLAG == 2 ) {
! 871: for (R=1,I=0;I<length(S);I++)
! 872: R *= S[I];
! 873: return [R];
! 874: }
! 875: /* FLAG == 0 || FLAG == 1 */
! 876: for (T=[]; S!=[];S = cdr(S)) { /* FLAG <= 1 and length(S) != 1 */
! 877: A = car(S);
! 878: U = append(T,cdr(S));
! 879: for (R=1,I=0;I<length(U);I++)
! 880: R *= U[I];
! 881: for (I=0;I<length(PD);I++) {
! 882: if ( PD[I] != P && !idealinc([R],PD[I],VL,PRIMAORD) )
! 883: break;
! 884: }
! 885: if ( I != length(PD) )
! 886: T = append(T,[A]);
! 887: }
! 888: if ( FLAG )
! 889: return T;
! 890: for (R=1,I=0;I<length(T);I++)
! 891: R *= T[I];
! 892: return [R]; /* if FLAG == 0 */
! 893: }
! 894:
! 895: def minselect(L,VL,ORD)
! 896: {
! 897: dp_ord(ORD);
! 898: F = dp_ptod(L[0],VL); MAG=dp_mag(F); DEG=dp_td(F);
! 899: for (J = 0, I = 1; I < length(L); I++) {
! 900: A = dp_ptod(L[I],VL); M=dp_mag(A); D=dp_td(A);
! 901: if ( dp_ht(A) > dp_ht(F) )
! 902: continue;
! 903: else if ( dp_ht(A) == dp_ht(F) && (D > DEG || (D == DEG && M > MAG)))
! 904: continue;
! 905: F = A; J = I; MAG = M; DEG=D;
! 906: }
! 907: return L[J];
! 908: }
! 909:
! 910: /* localization Id(F)_MI \cap Q[VL] */
! 911: /* output is the G-base w.r.t ordering O */
! 912: def localize(F,MI,VL,O)
! 913: {
! 914: if ( MI == [1] || MI == [-1] )
! 915: return F;
! 916: for (SVL = [],R = [],I = 0; I < length(MI); I++) {
! 917: V = strtov("zz"+rtostr(I));
! 918: R = append(R,[V*MI[I]-1]);
! 919: SVL = append(SVL,[V]);
! 920: }
! 921: if ( O == 0 ) {
! 922: dp_nelim(length(SVL)); ORD = 6;
! 923: } else if ( O == 1 || O == 2 ) {
! 924: ORD = [[0,length(SVL)],[O,length(VL)]];
! 925: } else if ( O == 3 || O == 4 || O == 5 )
! 926: ORD = [[0,length(SVL)],[O-3,length(VL)-1],[2,1]];
! 927: else
! 928: error("invarid ordering");
! 929: GR(G,append(F,R),append(SVL,VL),ORD);
! 930: S = varminus(G,SVL);
! 931: return S;
! 932: }
! 933:
! 934: def varminus(G,VL)
! 935: {
! 936: for (S = [],I = 0; I < length(G); I++) {
! 937: V = vars(G[I]);
! 938: if (listminus(V,VL) == V)
! 939: S = append(S,[G[I]]);
! 940: }
! 941: return S;
! 942: }
! 943:
! 944: def idealnormal(L)
! 945: {
! 946: R=[];
! 947: for (I=length(L)-1;I>=0;I--) {
! 948: A = ptozp(L[I]);
! 949: V = vars(A);
! 950: for (B = A,J = 0;J < length(V);J++)
! 951: B = coef(B,deg(B,V[J]),V[J]);
! 952: R = cons((B>0?A:-A),R);
! 953: }
! 954: return R;
! 955: }
! 956:
! 957: def idealsav(L) /* sub procedure of welldec and normposdec */
! 958: {
! 959: if ( PRINTLOG >= 4 )
! 960: print("multiple check ",2);
! 961: for (R = [],I = 0,N=length(L); I < N; I++) {
! 962: if ( PRINTLOG >= 4 && !irem(I,10) )
! 963: print(".",2);
! 964: if ( R == [] )
! 965: R = append(R,[L[I]]);
! 966: else {
! 967: for (J = 0; J < length(R); J++)
! 968: if ( idealeq(L[I],R[J]) )
! 969: break;
! 970: if ( J == length(R) )
! 971: R = append(R,[L[I]]);
! 972: }
! 973: }
! 974: if ( PRINTLOG >= 4 )
! 975: print("done.");
! 976: return R;
! 977: }
! 978:
! 979: /* intersection of ideals in PL, PL : list of ideals */
! 980: /* VL : variable list, O : output the G-base w.r.t order O */
! 981: def idealsetcap(PL,VL,O)
! 982: {
! 983: for (U = [], I = 0; I < length(PL); I++)
! 984: if ( PL[I] != [] )
! 985: U = append(U,[PL[I]]);
! 986: if ( U == [] )
! 987: return [];
! 988: if ( PRINTLOG >= 4 ) {print("intersection of ideal ",0); print(length(U),0);}
! 989: for (L = U[0],I=1;I<length(U);I++) {
! 990: if ( PRINTLOG >=4 ) print(".",2);
! 991: L = idealcap(L,U[I],VL,O);
! 992: }
! 993: if ( PRINTLOG >=4 ) print("");
! 994: return L;
! 995: }
! 996:
! 997: /* return intersection set of ideals P and Q */
! 998: def idealcap(P,Q,VL,O)
! 999: {
! 1000: if ( P == [] )
! 1001: if ( VL )
! 1002: return Q;
! 1003: else
! 1004: return [];
! 1005: if ( Q == [] )
! 1006: return P;
! 1007: T = tmp;
! 1008: for (PP = [],I = 0; I < length(P); I++)
! 1009: PP = append(PP,[P[I]*T]);
! 1010: for (QQ = [],I = 0; I < length(Q); I++)
! 1011: QQ = append(QQ,[Q[I]*(T-1)]);
! 1012: if ( O == 0 ) {
! 1013: dp_nelim(1); ORD = 6;
! 1014: } else if ( O == 1 || O == 2 )
! 1015: ORD = [[2,1],[O,length(VL)]];
! 1016: else if ( O == 3 || O == 4 || O == 5 )
! 1017: ORD = [[2,1],[O-3,length(VL)-1],[2,1]];
! 1018: else
! 1019: error("invarid ordering");
! 1020: GR(G,append(PP,QQ),append([T],VL),ORD);
! 1021: S = varminus(G,[T]);
! 1022: return S;
! 1023: }
! 1024:
! 1025: /* return ideal P : Q */
! 1026: def idealquo(P,Q,VL,O)
! 1027: {
! 1028: for (R = [], I = 0; I < length(Q); I++) {
! 1029: F = car(Q);
! 1030: G = idealcap(P,[F],VL,O);
! 1031: for (H = [],J = 0; J < length(G); J++) {
! 1032: H = append(H,[sdiv(G[J],F)]);
! 1033: }
! 1034: R = idealcap(R,H,VL,O);
! 1035: }
! 1036: }
! 1037:
! 1038: /* if ideal Q includes P then return 1 else return 0 */
! 1039: /* Q must be a Groebner Basis w.r.t ordering ORD */
! 1040: def idealinc(P,Q,VL,ORD)
! 1041: {
! 1042: dp_ord(ORD);
! 1043: for (T=IN=[],I=length(Q)-1;I>=0;I--) {
! 1044: T=cons(dp_ptod(Q[I],VL),T);
! 1045: IN=cons(I,IN);
! 1046: }
! 1047: DQ = newvect(length(Q),T);
! 1048: for (I = 0; I < length(P); I++) {
! 1049: NF(E,IN,dp_ptod(P[I],VL),DQ,0);
! 1050: if ( E )
! 1051: return 0;
! 1052: }
! 1053: return 1;
! 1054: }
! 1055:
! 1056: /* FL : list of polynomial set */
! 1057: def primedec_main(FL,VL)
! 1058: {
! 1059: if ( PRINTLOG ) {
! 1060: print("prime decomposition of the radical");
! 1061: }
! 1062: if ( PRINTLOG >= 2 )
! 1063: print("G-base factorization");
! 1064: PP = gr_fctr(FL,VL);
! 1065: if ( PRINTLOG >= 2 )
! 1066: print("irreducible variety decomposition");
! 1067: RP = welldec(PP,VL);
! 1068: SP = normposdec(RP,VL);
! 1069: for (L=[],I=length(SP)-1;I>=0;I--) {
! 1070: L=cons(idealnormal(SP[I]),L); /* head coef. > 0 */
! 1071: }
! 1072: SP = L;
! 1073: return SP;
! 1074: }
! 1075:
! 1076: def gr_fctr(FL,VL)
! 1077: {
! 1078: for (TP = [],I = 0; I<length(FL); I++ ) {
! 1079: F = FL[I];
! 1080: SF = idealsqfr(F);
! 1081: if ( !idealeq(F,SF) ) { GR(F,SF,VL,PRIMEORD); }
! 1082: DIRECTRY = []; COMMONIDEAL=[1];
! 1083: SP = gr_fctr_sub(F,VL);
! 1084: TP = append(TP,SP);
! 1085: }
! 1086: TP = idealsav(TP);
! 1087: TP = prime_irred(TP,VL);
! 1088: return TP;
! 1089: }
! 1090:
! 1091: def gr_fctr_sub(G,VL)
! 1092: {
! 1093: CONTINUE;
! 1094: if ( length(G) == 1 && type(G[0]) == 1 )
! 1095: return [G];
! 1096: RL = [];
! 1097: for (I = 0; I < length(G); I++) {
! 1098: FL = fctr(G[I]); L = length(FL); N = FL[1][1];
! 1099: if (L > 2 || N > 1) {
! 1100: TLL = [];
! 1101: for (J = 1; J < L; J++) {
! 1102: if ( CONTINUE ) {
! 1103: JCOPP=bload("jcopp");
! 1104: J = JCOPP[0]+1;
! 1105: COMMONIDEAL = JCOPP[1];
! 1106: RL = JCOPP[2];
! 1107: CONTINUE = 0;
! 1108: }
! 1109: if ( PRINTLOG >= 5 ) {
! 1110: for (D = length(DIRECTRY)-1; D >= 0; D--)
! 1111: print(DIRECTRY[D],0);
! 1112: print([L-1,J,length(RL)]);
! 1113: /*
! 1114: L-1:a number of factors of polynomial
! 1115: J:a factor executed now [0,...L-1]
! 1116: length(RL):a number of prime components
! 1117: */
! 1118: }
! 1119: W = cons(FL[J][0],G);
! 1120: GR(NG,W,VL,PRIMEORD);
! 1121: TNG = idealsqfr(NG);
! 1122: if ( NG != TNG ) { GR(NG,TNG,VL,PRIMEORD); }
! 1123: if ( idealinc(COMMONIDEAL,NG,VL,PRIMEORD) )
! 1124: continue;
! 1125: else {
! 1126: DIRECTRY=cons(L-J-1,DIRECTRY);
! 1127: DG = gr_fctr_sub(NG,VL);
! 1128: DIRECTRY=cdr(DIRECTRY);
! 1129: RL = append(RL,DG);
! 1130: if ( J <= L-2 && !idealinc(COMMONIDEAL,NG,VL,PRIMEORD)
! 1131: && COMMONCHECK ) {
! 1132: if ( PRINTLOG >= 5 ) {
! 1133: for (D = 0; D < length(DIRECTRY); D++) print(" ",2);
! 1134: print("intersection ideal");
! 1135: }
! 1136: COMMONIDEAL=idealcap(COMMONIDEAL,NG,VL,PRIMEORD);
! 1137: }
! 1138: }
! 1139: if ( BSAVE && !length(DIRECTRY) ) {
! 1140: bsave([J,COMMONIDEAL,RL],"jcopp");
! 1141: if ( PRINTLOG >= 2 )
! 1142: print("bsave j, intersection ideal and primes to icopp ");
! 1143: }
! 1144: }
! 1145: break;
! 1146: }
! 1147: }
! 1148: if (I == length(G)) {
! 1149: RL = append([G],RL);
! 1150: if ( PRINTLOG >= 5 ) {
! 1151: for (D = 0; D < length(DIRECTRY)-1; D++) print(" ",0);
! 1152: print("prime G-base ",0);
! 1153: if ( PRINTLOG >= 6 )
! 1154: print(G);
! 1155: else
! 1156: print("");
! 1157: }
! 1158: }
! 1159: return RL;
! 1160: }
! 1161:
! 1162: def prime_irred(TP,VL)
! 1163: {
! 1164: if ( PRINTLOG >= 4 )
! 1165: print("irredundancy check for prime ideal");
! 1166: N = length(TP);
! 1167: for (P = [], I = 0; I < N; I++) {
! 1168: if ( PRINTLOG >= 4 )
! 1169: print(".",2);
! 1170: for (J = 0; J < N; J++)
! 1171: if ( I != J && idealinc(TP[J],TP[I],VL,PRIMEORD) )
! 1172: break;
! 1173: if (J == N)
! 1174: P = append(P,[TP[I]]);
! 1175: }
! 1176: if ( PRINTLOG >= 4 )
! 1177: print("done.");
! 1178: return P;
! 1179: }
! 1180:
! 1181: /* if V1 \subset V2 then return 1 else return 0 */
! 1182: def varinc(V1,V2)
! 1183: {
! 1184: N1=length(V1); N2=length(V2);
! 1185: for (I=N1-1;I>=0;I--) {
! 1186: for (J=N2-1;J>=0;J--)
! 1187: if (V1[I]==V2[J])
! 1188: break;
! 1189: if (J==-1)
! 1190: return 0;
! 1191: }
! 1192: return 1;
! 1193: }
! 1194:
! 1195: def setunion(PS)
! 1196: {
! 1197: for (R=C=[]; PS != [] && car(PS); PS=cdr(PS)) {
! 1198: for (L = car(PS); L != []; L = cdr(L)) {
! 1199: A = car(L);
! 1200: for (G = C; G != []; G = cdr(G)) {
! 1201: if ( A == car(G) || -A == car(G) )
! 1202: break;
! 1203: }
! 1204: if ( G == [] ) {
! 1205: R = append(R,[A]);
! 1206: C = cons(A,C);
! 1207: }
! 1208: }
! 1209: }
! 1210: return R;
! 1211: }
! 1212:
! 1213: def idealeq(L,M)
! 1214: {
! 1215: if ((N1=length(L)) != (N2=length(M)))
! 1216: return 0;
! 1217: for (I = 0; I < length(L); I++) {
! 1218: for (J = 0; J < length(M); J++)
! 1219: if (L[I] == M[J] || L[I] == - M[J])
! 1220: break;
! 1221: if (J == length(M))
! 1222: return 0;
! 1223: }
! 1224: return 1;
! 1225: }
! 1226:
! 1227: def listminus(L,M)
! 1228: {
! 1229: for (R = []; L != []; L = cdr(L) ) {
! 1230: A = car(L);
! 1231: for (G = M; G != []; G = cdr(G)) {
! 1232: if ( A == car(G) )
! 1233: break;
! 1234: }
! 1235: if ( G == [] ) {
! 1236: R = append(R,[A]);
! 1237: M = cons(A,M);
! 1238: }
! 1239: }
! 1240: return R;
! 1241: }
! 1242:
! 1243: def idealsqfr(G)
! 1244: {
! 1245: for(Flag=0,LL=[],I=length(G)-1;I>=0;I--) {
! 1246: for(A=1,L=sqfr(G[I]),J=1;J<length(L);J++)
! 1247: A*=L[J][0];
! 1248: LL=cons(A,LL);
! 1249: }
! 1250: return LL;
! 1251: }
! 1252:
! 1253: def welldec(PRIME,VL)
! 1254: {
! 1255: PP = PRIME; NP = [];
! 1256: while ( !Flag ) {
! 1257: LP = [];
! 1258: Flag = 1;
! 1259: for (I=0;I<length(PP);I++) {
! 1260: U = welldec_sub(PP[I],VL);
! 1261: if (length(U) >= 2) {
! 1262: Flag = 0;
! 1263: if ( PRINTLOG >= 3 ) {
! 1264: print("now decomposing to irreducible varieties ",0);
! 1265: print(I,0); print(" ",0); print(length(PP));
! 1266: }
! 1267: DIRECTRY = []; COMMONIDEAL=[1];
! 1268: PL1 = gr_fctr([U[0]],VL);
! 1269: DIRECTRY = []; COMMONIDEAL=[1];
! 1270: PL2 = gr_fctr([U[1]],VL);
! 1271: LP = append(LP,append(PL1,PL2));
! 1272: }
! 1273: else
! 1274: NP = append(NP,U);
! 1275: }
! 1276: PP = LP;
! 1277: if ( PRINTLOG > 3 ) {
! 1278: print("prime length and non-prime length = ",0);
! 1279: print(length(NP),0); print(" ",0); print(length(PP));
! 1280: }
! 1281: }
! 1282: if ( length(PRIME) != length(NP) ) {
! 1283: NP = idealsav(NP);
! 1284: NP = prime_irred(NP,VL);
! 1285: }
! 1286: return NP;
! 1287: for (I = 0; I<length(PP);I++) {
! 1288: }
! 1289: }
! 1290:
! 1291: def welldec_sub(PP,VL)
! 1292: {
! 1293: VV = minalgdep(PP,VL,PRIMEORD);
! 1294: S = wellsep(PP,VV,VL);
! 1295: if ( S == 1 )
! 1296: return [PP];
! 1297: P1 = localize(PP,[S],VL,PRIMEORD);
! 1298: if ( idealeq(P1,PP) )
! 1299: return([PP]);
! 1300: GR(P2,cons(S,PP),VL,PRIMEORD);
! 1301: return [P1,P2];
! 1302: }
! 1303:
! 1304: def wellsep(PP,VV,VL)
! 1305: {
! 1306: TMPORD = 0;
! 1307: V0 = listminus(VL,VV);
! 1308: V1 = append(VV,V0);
! 1309: /* ORD = [[TMPORD,length(VV)],[0,length(V0)]]; */
! 1310: dp_nelim(length(VV)); ORD = 6;
! 1311: GR(PP,PP,V1,ORD);
! 1312: dp_ord(TMPORD);
! 1313: for (E=1,I=0;I<length(PP);I++)
! 1314: E = lcm(E,dp_hc(dp_ptod(PP[I],VV)));
! 1315: for (F=1,L=sqfr(E),K=1;K<length(L);K++)
! 1316: F *= L[K][0];
! 1317: return F;
! 1318: }
! 1319:
! 1320: /* prime decomposition by using primitive element method */
! 1321: def normposdec(NP,VL)
! 1322: {
! 1323: if ( PRINTLOG >= 3 )
! 1324: print("radical decomposition by using normal position.");
! 1325: for (R=MP=[],I=0;I<length(NP);I++) {
! 1326: V=minalgdep(NP[I],VL,PRIMEORD);
! 1327: L=raddec(NP[I],V,VL,1);
! 1328: if ( PRINTLOG >= 3 )
! 1329: if ( length(L) == 1 ) {
! 1330: print(".",2);
! 1331: } else {
! 1332: print(" ",2); print(length(L),2); print(" ",2);
! 1333: }
! 1334: /* if (length(L)==1) {
! 1335: MP=append(MP,[NP[I]]);
! 1336: continue;
! 1337: } */
! 1338: R=append(R,L);
! 1339: }
! 1340: if ( PRINTLOG >= 3 )
! 1341: print("done");
! 1342: if ( length(R) )
! 1343: MP = idealsav(append(R,MP));
! 1344: LP = prime_irred(MP,VL);
! 1345: return LP;
! 1346: }
! 1347:
! 1348: /* radical decomposition program using by Zianni, Trager, Zacharias */
! 1349: /* R : factorized G-base in K[X], if FLAG then A is well comp. */
! 1350: def raddec(R,V,X,FLAG)
! 1351: {
! 1352: GR(A,R,V,irem(PRIMEORD,3));
! 1353: ZQ=zraddec(A,V); /* zero dimensional */
! 1354: /* if ( FLAG && length(ZQ) == 1 )
! 1355: return [R]; */
! 1356: if ( length(V) != length(X) )
! 1357: CQ=radcont(ZQ,V,X); /* contraction */
! 1358: else {
! 1359: for (CQ=[],I=length(ZQ)-1;I>=0;I--) {
! 1360: GR(R,ZQ[I],X,PRIMEORD);
! 1361: CQ=cons(R,CQ);
! 1362: }
! 1363: }
! 1364: return CQ;
! 1365: }
! 1366:
! 1367: /* radical decomposition for zero dimensional ideal */
! 1368: /* F is G-base w.r.t irem(PRIMEORD,3) */
! 1369: def zraddec(F,X)
! 1370: {
! 1371: if (length(F) == 1)
! 1372: return [F];
! 1373: Z=vvv;
! 1374: C=normposit(F,X,Z);
! 1375: /* C=[minimal polynomial, adding polynomial] */
! 1376: T=cdr(fctr(C[0]));
! 1377: if ( length(T) == 1 && T[0][1] == 1 )
! 1378: return [F];
! 1379: for (Q=[],I=0;I<length(T);I++) {
! 1380: if ( !C[1] ) {
! 1381: GR(P,cons(T[I][0],F),X,irem(PRIMEORD,3));
! 1382: } else {
! 1383: P=idealgrsub(append([T[I][0],C[1]],F),X,Z);
! 1384: }
! 1385: Q=cons(P,Q);
! 1386: }
! 1387: return Q;
! 1388: }
! 1389:
! 1390: /* contraction from V to X */
! 1391: def radcont(Q,V,X)
! 1392: {
! 1393: for (R=[],I=length(Q)-1;I>=0;I--) {
! 1394: dp_ord(irem(PRIMEORD,3));
! 1395: G=Q[I];
! 1396: for (E=1,J=0;J<length(G);J++)
! 1397: E = lcm(E,dp_hc(dp_ptod(G[J],V)));
! 1398: for (F=1,L=sqfr(E),K=1;K<length(L);K++)
! 1399: F *= L[K][0];
! 1400: H=localize(G,[F],X,PRIMEORD);
! 1401: R=cons(H,R);
! 1402: }
! 1403: return R;
! 1404: }
! 1405:
! 1406: /* A : polynomials (G-Base) */
! 1407: /* return [T,Y], T : minimal polynomial of Z, Y : append polynomial */
! 1408: def normposit(A,X,Z)
! 1409: {
! 1410: D = idim(A,X,irem(PRIMEORD,3)); /* dimension of the ideal Id(A) */
! 1411: for (I = length(A)-1;I>=0;I--) {
! 1412: for (J = length(X)-1; J>= 0; J--) {
! 1413: T=deg(A[I],X[J]);
! 1414: if ( T == D ) /* A[I] is a primitive polynomial of A */
! 1415: return [A[I],0];
! 1416: }
! 1417: }
! 1418: N=length(X);
! 1419: for (C = [],I = 0; I < N-1; I++)
! 1420: C=newvect(N-1);
! 1421: V=append(X,[Z]);
! 1422: ZD = (vars(A)==vars(X)); /* A is zero dim. over Q[X] or not */
! 1423: while( 1 ) {
! 1424: C = nextweight(C,cdr(X));
! 1425: for (Y=Z-X[0],I=0;I<N-1;I++) {
! 1426: Y-=C[I]*X[I+1]; /* new polynomial */
! 1427: }
! 1428: if ( !ZD ) {
! 1429: if ( version() == 940420 ) NOREDUCE = 1;
! 1430: else dp_gr_flags(["NoRA",1]);
! 1431: GR(G,cons(Y,A),V,2);
! 1432: if ( version() == 940420 ) NOREDUCE = 0;
! 1433: else dp_gr_flags(["NoRA",0]);
! 1434: for (T=G[length(G)-1],I = length(G)-2;I >= 0; I--)
! 1435: if ( deg(G[I],Z) > deg(T,Z) )
! 1436: T = G[I]; /* minimal polynomial */
! 1437: } else {
! 1438: T = minipoly(A,X,irem(PRIMEORD,3),Z-Y,Z);
! 1439: }
! 1440: if (deg(T,Z)==D)
! 1441: return [T,Y];
! 1442: }
! 1443: }
! 1444:
! 1445: def nextweight(C,V) /* degrevlex */
! 1446: {
! 1447: dp_ord(2);
! 1448: N = length(V);
! 1449: for (D=I=0;I<size(C)[0];I++)
! 1450: D += C[I];
! 1451: if ( C[N-1] == D ) {
! 1452: for (L=[],I=0;I<N-1;I++)
! 1453: L=cons(0,L);
! 1454: L = cons(D+1,L);
! 1455: return newvect(N,L);
! 1456: }
! 1457: CD=dp_vtoe(C);
! 1458: for (F=I=0;I<N;I++)
! 1459: F+=V[I];
! 1460: for (DF=dp_ptod(F^D,V);dp_ht(DF)!=CD;DF=dp_rest(DF));
! 1461: MD = dp_ht(dp_rest(DF));
! 1462: CC = dp_etov(MD);
! 1463: return CC;
! 1464: }
! 1465:
! 1466: def printsep(L)
! 1467: {
! 1468: for (I=0;I<length(L);I++) {
! 1469: if ( nmono(L[I][0]) == 1 )
! 1470: print(L[I][0],0);
! 1471: else {
! 1472: print("(",0); print(L[I][0],0); print(")",0);
! 1473: }
! 1474: if (L[I][1] > 1) {
! 1475: print("^",0); print(L[I][1],0);
! 1476: }
! 1477: if (I<length(L)-1)
! 1478: print("*",0);
! 1479: }
! 1480: }
! 1481:
! 1482: def idealgrsub(A,X,Z)
! 1483: {
! 1484: if ( PRIMEORD == 0 ) {
! 1485: dp_nelim(1); ORD = 8;
! 1486: } else
! 1487: ORD = [[2,1],[PRIMERORD,length(X)]];
! 1488: GR(G,A,cons(Z,X),ORD);
! 1489: for (R=[],I=length(G)-1;I>=0;I--) {
! 1490: V=vars(G[I]);
! 1491: for (J=0;J<length(V);J++)
! 1492: if (V[J]==Z)
! 1493: break;
! 1494: if (J==length(V))
! 1495: R=cons(G[I],R);
! 1496: }
! 1497: return R;
! 1498: }
! 1499:
! 1500: /* dimension of ideal */
! 1501: def idim(F,V,ORD)
! 1502: {
! 1503: return length(modbase(F,V,ORD));
! 1504: }
! 1505:
! 1506: def modbase(F,V,ORD)
! 1507: {
! 1508: dp_ord(ORD);
! 1509: for(G=[],I=length(F)-1;I>=0;I--)
! 1510: G = cons(dp_ptod(F[I],V),G);
! 1511: R = dp_modbase(G,length(V));
! 1512: for(S=[],I=length(R)-1;I>=0;I--)
! 1513: S = cons(dp_dtop(R[I],V),S);
! 1514: return S;
! 1515: }
! 1516:
! 1517: def dp_modbase(G,N)
! 1518: {
! 1519: E = newvect(N);
! 1520: R = []; I = 1;
! 1521: for (L = [];G != []; G = cdr(G))
! 1522: L = cons(dp_ht(car(G)),L);
! 1523: while ( I <= N ) {
! 1524: R = cons(dp_vtoe(E),R);
! 1525: E[0]++; I = 1;
! 1526: while( s_redble(dp_vtoe(E),L) ) {
! 1527: E[I-1] = 0;
! 1528: if (I < N)
! 1529: E[I]++;
! 1530: I++;
! 1531: }
! 1532: }
! 1533: return R;
! 1534: }
! 1535:
! 1536: def s_redble(M,L)
! 1537: {
! 1538: for (; L != []; L = cdr(L))
! 1539: if ( dp_redble(M,car(L)) )
! 1540: return 1;
! 1541: return 0;
! 1542: }
! 1543:
! 1544: /* FL : list of ideal Id(P,s) */
! 1545: def primedec_special(FL,VL)
! 1546: {
! 1547: if ( PRINTLOG ) {
! 1548: print("prime decomposition of the radical");
! 1549: }
! 1550: if ( PRINTLOG >= 2 )
! 1551: print("G-base factorization");
! 1552: PP = gr_fctr(FL,VL);
! 1553: if ( PRINTLOG >= 2 )
! 1554: print("special decomposition");
! 1555: SP = yokodec(PP,VL);
! 1556: for (L=[],I=length(SP)-1;I>=0;I--) {
! 1557: L=cons(idealnormal(SP[I]),L); /* head coef. > 0 */
! 1558: }
! 1559: SP = L;
! 1560: return SP;
! 1561: }
! 1562:
! 1563: /* PL : list of ideal Id(P,s) */
! 1564: def yokodec(PL,VL)
! 1565: {
! 1566: MSISTIME = 0; T = time()[0];
! 1567: for (R = [],I = 0; I<length(PL);I++) {
! 1568: PP = PL[I];
! 1569: if ( PRINTLOG >= 3 ) print(".",2);
! 1570: V = minalgdep(PP,VL,PRIMEORD);
! 1571: E = raddec(PP,V,VL,0);
! 1572: if ( length(E) >= 2 || !idealeq(E[0],PP) ) { /* prime check */
! 1573: T0 = time()[0];
! 1574: L=all_msis(PP,VL,PRIMEORD);
! 1575: MSISTIME += time()[0]-T0;
! 1576: E = yokodec_main(PP,L,VL);
! 1577: }
! 1578: R = append(R,E);
! 1579: }
! 1580: R = idealsav(R);
! 1581: R = prime_irred(R,VL);
! 1582: if ( PRINTLOG >= 3 ) {
! 1583: print("special dec time = ",0); print(time()[0]-T0);
! 1584: /* print(", ",0); print(MSISTIME); */
! 1585: }
! 1586: return R;
! 1587: }
! 1588:
! 1589: def yokodec_main(PP,L,VL)
! 1590: {
! 1591: AL = misset(L,VL); H = [1]; R = [];
! 1592: for (P=PP,I=0; I<length(AL); I++) {
! 1593: V = AL[I];
! 1594: if ( I && !ismad(P,AL[I],VL,PRIMEORD) )
! 1595: continue;
! 1596: U = raddec(PP,V,VL,0);
! 1597: R = append(R,U);
! 1598: for (I=0;I<length(U);I++)
! 1599: H = idealcap(H,U[I],VL,PRIMEORD);
! 1600: if ( idealeq(H,P) )
! 1601: break;
! 1602: F = wellsep(P,V,VL);
! 1603: GR(P,cons(F,P),VL,PRIMEORD);
! 1604: }
! 1605: return R;
! 1606: }
! 1607:
! 1608: /* output M.A.D. sets. AL : M.S.I.S sets */
! 1609: def misset(AL,VL)
! 1610: {
! 1611: for (L=[],D=0,I=length(AL)-1;I>=0;I--) {
! 1612: E = length(AL[I]);
! 1613: C = listminus(VL,AL[I]);
! 1614: if ( E == D )
! 1615: L = cons(C,L);
! 1616: else if ( E > D ) {
! 1617: L = [C]; D = E;
! 1618: }
! 1619: }
! 1620: return L;
! 1621: }
! 1622:
! 1623: end$
! 1624:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to primdec_lex CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2018 / lib