![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to gr CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2018 / lib|
Return to gr CVS log | Up to [local] / OpenXM_contrib2 / asir2018 / lib |
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:
! 51: module gr $
! 52: /* Empty for now. It will be used in a future. */
! 53: endmodule $
! 54:
! 55: extern INIT_COUNT,ITOR_FAIL$
! 56: extern REMOTE_MATRIX,REMOTE_NF,REMOTE_VARS$
! 57:
! 58: #define MAX(a,b) ((a)>(b)?(a):(b))
! 59: #define HigherDim 0
! 60: #define ZeroDim 1
! 61: #define MiniPoly 2
! 62:
! 63: /* toplevel functions for Groebner basis computation */
! 64:
! 65: def gr(B,V,O)
! 66: {
! 67: Ord=dp_ord();
! 68: G = dp_gr_main(B,V,0,1,O);
! 69: dp_ord(Ord);
! 70: return G;
! 71: }
! 72:
! 73: def hgr(B,V,O)
! 74: {
! 75: Ord=dp_ord();
! 76: G = dp_gr_main(B,V,1,1,O);
! 77: dp_ord(Ord);
! 78: return G;
! 79: }
! 80:
! 81: def gr_mod(B,V,O,M)
! 82: {
! 83: Ord=dp_ord();
! 84: G = dp_gr_mod_main(B,V,0,M,O);
! 85: dp_ord(Ord);
! 86: return G;
! 87: }
! 88:
! 89: def hgr_mod(B,V,O,M)
! 90: {
! 91: Ord=dp_ord();
! 92: G = dp_gr_mod_main(B,V,1,M,O);
! 93: dp_ord(Ord);
! 94: return G;
! 95: }
! 96:
! 97: /* toplevel functions for change-of-ordering */
! 98:
! 99: def lex_hensel(B,V,O,W,H)
! 100: {
! 101: G = dp_gr_main(B,V,H,1,O);
! 102: return tolex(G,V,O,W);
! 103: }
! 104:
! 105: def lex_hensel_gsl(B,V,O,W,H)
! 106: {
! 107: G = dp_gr_main(B,V,H,1,O);
! 108: return tolex_gsl(G,V,O,W);
! 109: }
! 110:
! 111: def gr_minipoly(B,V,O,P,V0,H)
! 112: {
! 113: G = dp_gr_main(B,V,H,1,O);
! 114: return minipoly(G,V,O,P,V0);
! 115: }
! 116:
! 117: def lex_tl(B,V,O,W,H)
! 118: {
! 119: G = dp_gr_main(B,V,H,1,O);
! 120: return tolex_tl(G,V,O,W,H);
! 121: }
! 122:
! 123: def tolex_tl(G0,V,O,W,H)
! 124: {
! 125: N = length(V); HM = hmlist(G0,V,O); ZD = zero_dim(HM,V,O);
! 126: for ( I = 0; ; I++ ) {
! 127: M = lprime(I);
! 128: if ( !valid_modulus(HM,M) )
! 129: continue;
! 130: if ( ZD ) {
! 131: if ( G3 = dp_gr_main(G0,W,H,-M,3) )
! 132: for ( J = 0; ; J++ )
! 133: if ( G2 = dp_gr_main(G3,W,0,-lprime(J),2) )
! 134: return G2;
! 135: } else if ( G2 = dp_gr_main(G0,W,H,-M,2) )
! 136: return G2;
! 137: }
! 138: }
! 139:
! 140: def tolex(G0,V,O,W)
! 141: {
! 142: Procs = getopt(procs);
! 143:
! 144: TM = TE = TNF = 0;
! 145: N = length(V); HM = hmlist(G0,V,O); ZD = zero_dim(HM,V,O);
! 146: if ( ZD )
! 147: MB = dp_mbase(map(dp_ptod,HM,V));
! 148: else
! 149: MB = 0;
! 150: for ( J = 0; ; J++ ) {
! 151: M = lprime(J);
! 152: if ( !valid_modulus(HM,M) )
! 153: continue;
! 154: T0 = time()[0];
! 155: if ( ZD ) {
! 156: GM = tolexm(G0,V,O,W,M);
! 157: dp_ord(2);
! 158: DL = map(dp_etov,map(dp_ht,map(dp_ptod,GM,W)));
! 159: D = newvect(N); TL = [];
! 160: do
! 161: TL = cons(dp_dtop(dp_vtoe(D),W),TL);
! 162: while ( nextm(D,DL,N) );
! 163: } else {
! 164: HVN = "h";
! 165: WN = map(rtostr,W);
! 166: while ( member(HVN,WN) ) HVN += "h";
! 167: HV = strtov(HVN);
! 168: G0H = map(homogenize,G0,W,HV);
! 169: GMH = nd_gr(G0H,append(W,[HV]),M,1);
! 170: GMH=map(subst,GMH,HV,1);
! 171: GM = nd_gr_postproc(GMH,W,M,2,0);
! 172: dp_ord(2);
! 173: for ( T = GM, S = 0; T != []; T = cdr(T) )
! 174: for ( D = dp_ptod(car(T),V); D; D = dp_rest(D) )
! 175: S += dp_ht(D);
! 176: TL = dp_terms(S,V);
! 177: }
! 178: TM += time()[0] - T0;
! 179: T0 = time()[0]; NF = gennf(G0,TL,V,O,W[N-1],ZD)[0];
! 180: TNF += time()[0] - T0;
! 181: T0 = time()[0];
! 182: if ( type(Procs) != -1 )
! 183: R = tolex_d_main(V,O,NF,GM,M,MB,Procs);
! 184: else
! 185: R = tolex_main(V,O,NF,GM,M,MB);
! 186: TE += time()[0] - T0;
! 187: if ( R ) {
! 188: if ( dp_gr_print() )
! 189: print("mod="+rtostr(TM)+",nf="+rtostr(TNF)+",eq="+rtostr(TE));
! 190: return R;
! 191: }
! 192: }
! 193: }
! 194:
! 195: def tolex_gsl(G0,V,O,W)
! 196: {
! 197: TM = TE = TNF = 0;
! 198: N = length(V); HM = hmlist(G0,V,O); ZD = zero_dim(HM,V,O);
! 199: MB = dp_mbase(map(dp_ptod,HM,V));
! 200: if ( !ZD )
! 201: error("tolex_gsl : ideal is not zero-dimensional!");
! 202: for ( J = 0; ; J++ ) {
! 203: M = lprime(J);
! 204: if ( !valid_modulus(HM,M) )
! 205: continue;
! 206: T0 = time()[0]; GM = tolexm(G0,V,O,W,M); TM += time()[0] - T0;
! 207: dp_ord(2);
! 208: DL = map(dp_etov,map(dp_ht,map(dp_ptod,GM,W)));
! 209: D = newvect(N); TL = [];
! 210: do
! 211: TL = cons(dp_dtop(dp_vtoe(D),W),TL);
! 212: while ( nextm(D,DL,N) );
! 213: L = npos_check(DL); NPOSV = L[0]; DIM = L[1];
! 214: if ( NPOSV >= 0 ) {
! 215: if ( dp_gr_print() ) print("shape base");
! 216: V0 = W[NPOSV];
! 217: T0 = time()[0]; NFL = gennf(G0,TL,V,O,V0,1);
! 218: TNF += time()[0] - T0;
! 219: T0 = time()[0];
! 220: R = tolex_gsl_main(G0,V,O,W,NFL,NPOSV,GM,M,MB);
! 221: TE += time()[0] - T0;
! 222: } else {
! 223: T0 = time()[0]; NF = gennf(G0,TL,V,O,W[N-1],1)[0];
! 224: TNF += time()[0] - T0;
! 225: T0 = time()[0];
! 226: R = tolex_main(V,O,NF,GM,M,MB);
! 227: TE += time()[0] - T0;
! 228: }
! 229: if ( R ) {
! 230: if ( dp_gr_print() )
! 231: print("mod="+rtostr(TM)+",nf="+rtostr(TNF)+",eq="+rtostr(TE));
! 232: return R;
! 233: }
! 234: }
! 235: }
! 236:
! 237: def termstomat(NF,TERMS,MB,MOD)
! 238: {
! 239: DN = NF[1];
! 240: NF = NF[0];
! 241: N = length(MB);
! 242: M = length(TERMS);
! 243: MAT = newmat(N,M);
! 244: W = newvect(N);
! 245: Len = length(NF);
! 246: for ( I = 0; I < M; I++ ) {
! 247: T = TERMS[I];
! 248: for ( K = 0; K < Len; K++ )
! 249: if ( T == NF[K][1] )
! 250: break;
! 251: dptov(NF[K][0],W,MB);
! 252: for ( J = 0; J < N; J++ )
! 253: MAT[J][I] = W[J];
! 254: }
! 255: return [henleq_prep(MAT,MOD),DN];
! 256: }
! 257:
! 258: def tolex_gsl_main(G0,V,O,W,NFL,NPOSV,GM,M,MB)
! 259: {
! 260: NF = NFL[0]; PS = NFL[1]; GI = NFL[2];
! 261: V0 = W[NPOSV]; N = length(W);
! 262: DIM = length(MB);
! 263: DV = newvect(DIM);
! 264: TERMS = gather_terms(GM,W,M,NPOSV);
! 265: Len = length(TERMS);
! 266: dp_ord(O); RHS = termstomat(NF,map(dp_ptod,TERMS,V),MB,M);
! 267: for ( T = GM; T != []; T = cdr(T) )
! 268: if ( vars(car(T)) == [V0] )
! 269: break;
! 270: dp_ord(0); NHT = nf_tab_gsl(dp_ptod(V0^deg(car(T),V0),V),NF);
! 271: dptov(NHT[0],DV,MB);
! 272: B = hen_ttob_gsl([DV,NHT[1]],RHS,TERMS,M);
! 273: if ( !B )
! 274: return 0;
! 275: for ( I = 0, U = B[1]*V0^deg(car(T),V0); I < Len; I++ )
! 276: U += B[0][I]*TERMS[I];
! 277: DN0 = diff(U,V0);
! 278: dp_ord(O); DN0NF = nf_tab_gsl(dp_ptod(DN0,V),NF);
! 279: SL = [[V0,U,DN0]];
! 280: for ( I = N-1, LCM = 1; I >= 0; I-- ) {
! 281: if ( I == NPOSV )
! 282: continue;
! 283: V1 = W[I];
! 284: dp_ord(O); L = nf(GI,DN0NF[0]*dp_ptod(-LCM*V1,V),DN0NF[1],PS);
! 285: L = remove_cont(L);
! 286: dptov(L[0],DV,MB);
! 287: dp_ord(O); B = hen_ttob_gsl([DV,L[1]],RHS,TERMS,M);
! 288: if ( !B )
! 289: return 0;
! 290: for ( K = 0, R = 0; K < Len; K++ )
! 291: R += B[0][K]*TERMS[K];
! 292: LCM *= B[1];
! 293: SL = cons(cons(V1,[R,LCM]),SL);
! 294: if ( dp_gr_print() )
! 295: print(["DN",B[1]]);
! 296: }
! 297: return SL;
! 298: }
! 299:
! 300: def hen_ttob_gsl(LHS,RHS,TERMS,M)
! 301: {
! 302: LDN = LHS[1]; RDN = RHS[1]; LCM = ilcm(LDN,RDN);
! 303: L1 = idiv(LCM,LDN); R1 = idiv(LCM,RDN);
! 304: T0 = time()[0];
! 305: S = henleq_gsl(RHS[0],LHS[0]*L1,M);
! 306: if ( dp_gr_print() )
! 307: print(["henleq_gsl",time()[0]-T0]);
! 308: N = length(TERMS);
! 309: return [S[0],S[1]*R1];
! 310: }
! 311:
! 312: def gather_terms(GM,W,M,NPOSV)
! 313: {
! 314: N = length(W); V0 = W[NPOSV];
! 315: for ( T = GM; T != []; T = cdr(T) ) {
! 316: if ( vars(car(T)) == [V0] )
! 317: break;
! 318: }
! 319: U = car(T); DU = diff(U,V0);
! 320: R = tpoly(cdr(p_terms(U,W,2)));
! 321: for ( I = 0; I < N; I++ ) {
! 322: if ( I == NPOSV )
! 323: continue;
! 324: V1 = W[I];
! 325: for ( T = GM; T != []; T = cdr(T) )
! 326: if ( member(V1,vars(car(T))) )
! 327: break;
! 328: P = car(T);
! 329: R += tpoly(p_terms(srem(DU*coef(P,0,V1),U,M),W,2));
! 330: }
! 331: return p_terms(R,W,2);
! 332: }
! 333:
! 334: def tpoly(L)
! 335: {
! 336: for ( R = 0; L != []; L = cdr(L) )
! 337: R += car(L);
! 338: return R;
! 339: }
! 340:
! 341: def dptov(P,W,MB)
! 342: {
! 343: N = size(W)[0];
! 344: for ( I = 0; I < N; I++ )
! 345: W[I] = 0;
! 346: for ( I = 0, S = MB; P; P = dp_rest(P) ) {
! 347: HM = dp_hm(P); C = dp_hc(HM); T = dp_ht(HM);
! 348: for ( ; T != car(S); S = cdr(S), I++ );
! 349: W[I] = C;
! 350: I++; S = cdr(S);
! 351: }
! 352: }
! 353:
! 354: def tolex_main(V,O,NF,GM,M,MB)
! 355: {
! 356: if ( MB ) {
! 357: PosDim = 0;
! 358: DIM = length(MB);
! 359: DV = newvect(DIM);
! 360: } else
! 361: PosDim = 1;
! 362: for ( T = GM, SL = [], LCM = 1; T != []; T = cdr(T) ) {
! 363: S = p_terms(car(T),V,2);
! 364: if ( PosDim ) {
! 365: MB = gather_nf_terms(S,NF,V,O);
! 366: DV = newvect(length(MB));
! 367: }
! 368: dp_ord(O); RHS = termstomat(NF,map(dp_ptod,cdr(S),V),MB,M);
! 369: dp_ord(O); NHT = nf_tab_gsl(dp_ptod(LCM*car(S),V),NF);
! 370: dptov(NHT[0],DV,MB);
! 371: dp_ord(O); B = hen_ttob_gsl([DV,NHT[1]],RHS,cdr(S),M);
! 372: if ( !B )
! 373: return 0;
! 374: Len = length(S);
! 375: LCM *= B[1];
! 376: for ( U = LCM*car(S), I = 1; I < Len; I++ )
! 377: U += B[0][I-1]*S[I];
! 378: R = ptozp(U);
! 379: SL = cons(R,SL);
! 380: if ( dp_gr_print() )
! 381: print(["DN",B[1]]);
! 382: }
! 383: return SL;
! 384: }
! 385:
! 386: def tolex_d_main(V,O,NF,GM,M,MB,Procs)
! 387: {
! 388: map(ox_reset,Procs);
! 389: /* register data in servers */
! 390: map(ox_cmo_rpc,Procs,"register_data_for_find_base",NF,V,O,MB,M);
! 391: /* discard return value in stack */
! 392: map(ox_pop_cmo,Procs);
! 393: Free = Procs;
! 394: Busy = [];
! 395: T = GM;
! 396: SL = [];
! 397: while ( T != [] || Busy != [] ){
! 398: if ( Free == [] || T == [] ) {
! 399: /* someone is working; wait for data */
! 400: Ready = ox_select(Busy);
! 401: Busy = setminus(Busy,Ready);
! 402: Free = append(Ready,Free);
! 403: for ( ; Ready != []; Ready = cdr(Ready) )
! 404: SL = cons(ox_get(car(Ready)),SL);
! 405: } else {
! 406: P = car(Free);
! 407: Free = cdr(Free);
! 408: Busy = cons(P,Busy);
! 409: Template = car(T);
! 410: T = cdr(T);
! 411: ox_cmo_rpc(P,"find_base",Template);
! 412: ox_push_cmd(P,262); /* 262 = OX_popCMO */
! 413: }
! 414: }
! 415: return SL;
! 416: }
! 417:
! 418: struct find_base_data { NF,V,O,MB,M,PosDim,DV }$
! 419: extern Find_base$
! 420:
! 421: def register_data_for_find_base(NF,V,O,MB,M)
! 422: {
! 423: Find_base = newstruct(find_base_data);
! 424: Find_base->NF = NF;
! 425: Find_base->V = V;
! 426: Find_base->O = O;
! 427: Find_base->M = M;
! 428: Find_base->MB = MB;
! 429:
! 430: if ( MB ) {
! 431: Find_base->PosDim = 0;
! 432: DIM = length(MB);
! 433: Find_base->DV = newvect(DIM);
! 434: } else
! 435: Find_base->PosDim = 1;
! 436: }
! 437:
! 438: def find_base(S) {
! 439: NF = Find_base->NF;
! 440: V = Find_base->V;
! 441: O = Find_base->O;
! 442: MB = Find_base->MB;
! 443: M = Find_base->M;
! 444: PosDim = Find_base->PosDim;
! 445: DV = Find_base->DV;
! 446:
! 447: S = p_terms(S,V,2);
! 448: if ( PosDim ) {
! 449: MB = gather_nf_terms(S,NF,V,O);
! 450: DV = newvect(length(MB));
! 451: }
! 452: dp_ord(O); RHS = termstomat(NF,map(dp_ptod,cdr(S),V),MB,M);
! 453: dp_ord(O); NHT = nf_tab_gsl(dp_ptod(car(S),V),NF);
! 454: dptov(NHT[0],DV,MB);
! 455: dp_ord(O); B = hen_ttob_gsl([DV,NHT[1]],RHS,cdr(S),M);
! 456: if ( !B )
! 457: return 0;
! 458: Len = length(S);
! 459: for ( U = B[1]*car(S), I = 1; I < Len; I++ )
! 460: U += B[0][I-1]*S[I];
! 461: R = ptozp(U);
! 462: return R;
! 463: }
! 464:
! 465: /*
! 466: * NF = [Pairs,DN]
! 467: * Pairs = [[NF1,T1],[NF2,T2],...]
! 468: */
! 469:
! 470: def gather_nf_terms(S,NF,V,O)
! 471: {
! 472: R = 0;
! 473: for ( T = S; T != []; T = cdr(T) ) {
! 474: DT = dp_ptod(car(T),V);
! 475: for ( U = NF[0]; U != []; U = cdr(U) )
! 476: if ( car(U)[1] == DT ) {
! 477: R += tpoly(dp_terms(car(U)[0],V));
! 478: break;
! 479: }
! 480: }
! 481: return map(dp_ptod,p_terms(R,V,O),V);
! 482: }
! 483:
! 484: def reduce_dn(L)
! 485: {
! 486: NM = L[0]; DN = L[1]; V = vars(NM);
! 487: T = remove_cont([dp_ptod(NM,V),DN]);
! 488: return [dp_dtop(T[0],V),T[1]];
! 489: }
! 490:
! 491: /* a function for computation of minimal polynomial */
! 492:
! 493: def minipoly(G0,V,O,P,V0)
! 494: {
! 495: if ( !zero_dim(hmlist(G0,V,O),V,O) )
! 496: error("tolex : ideal is not zero-dimensional!");
! 497:
! 498: Pin = P;
! 499: P = ptozp(P);
! 500: CP = sdiv(P,Pin);
! 501: G1 = cons(V0-P,G0);
! 502: O1 = [[0,1],[O,length(V)]];
! 503: V1 = cons(V0,V);
! 504: W = append(V,[V0]);
! 505:
! 506: N = length(V1);
! 507: dp_ord(O1);
! 508: HM = hmlist(G1,V1,O1);
! 509: MB = dp_mbase(map(dp_ptod,HM,V1));
! 510: dp_ord(O);
! 511:
! 512: for ( J = 0; ; J++ ) {
! 513: M = lprime(J);
! 514: if ( !valid_modulus(HM,M) )
! 515: continue;
! 516: MP = minipolym(G0,V,O,P,V0,M);
! 517: for ( D = deg(MP,V0), TL = [], J = 0; J <= D; J++ )
! 518: TL = cons(V0^J,TL);
! 519: NF = gennf(G1,TL,V1,O1,V0,1)[0];
! 520: R = tolex_main(V1,O1,NF,[MP],M,MB);
! 521: return ptozp(subst(R[0],V0,CP*V0));
! 522: }
! 523: }
! 524:
! 525: /* subroutines */
! 526:
! 527: def gennf(G,TL,V,O,V0,FLAG)
! 528: {
! 529: F = dp_gr_flags();
! 530: for ( T = F; T != []; T = cdr(T) ) {
! 531: Key = car(T); T = cdr(T);
! 532: if ( Key == "Demand" ) {
! 533: Dir = car(T); break;
! 534: }
! 535: }
! 536: if ( Dir )
! 537: return gennf_demand(G,TL,V,O,V0,FLAG,Dir);
! 538: N = length(V); Len = length(G); dp_ord(O); PS = newvect(Len);
! 539: for ( I = 0, T = G, HL = []; T != []; T = cdr(T), I++ ) {
! 540: PS[I] = dp_ptod(car(T),V); HL = cons(dp_ht(PS[I]),HL);
! 541: }
! 542: for ( I = 0, DTL = []; TL != []; TL = cdr(TL) )
! 543: DTL = cons(dp_ptod(car(TL),V),DTL);
! 544: for ( I = Len - 1, GI = []; I >= 0; I-- )
! 545: GI = cons(I,GI);
! 546: T = car(DTL); DTL = cdr(DTL);
! 547: H = [nf(GI,T,T,PS)];
! 548:
! 549: USE_TAB = (FLAG != 0);
! 550: if ( USE_TAB ) {
! 551: T0 = time()[0];
! 552: MB = dp_mbase(HL); DIM = length(MB);
! 553: U = dp_ptod(V0,V);
! 554: UTAB = newvect(DIM);
! 555: for ( I = 0; I < DIM; I++ ) {
! 556: UTAB[I] = [MB[I],remove_cont(dp_true_nf(GI,U*MB[I],PS,1))];
! 557: if ( dp_gr_print() )
! 558: print(".",2);
! 559: }
! 560: if ( dp_gr_print() )
! 561: print("");
! 562: TTAB = time()[0]-T0;
! 563: }
! 564:
! 565: T0 = time()[0];
! 566: for ( LCM = 1; DTL != []; ) {
! 567: if ( dp_gr_print() )
! 568: print(".",2);
! 569: T = car(DTL); DTL = cdr(DTL);
! 570: if ( L = search_redble(T,H) ) {
! 571: DD = dp_subd(T,L[1]);
! 572: if ( USE_TAB && (DD == U) ) {
! 573: NF = nf_tab(L[0],UTAB);
! 574: NF = [NF[0],dp_hc(L[1])*NF[1]*T];
! 575: } else
! 576: NF = nf(GI,L[0]*dp_subd(T,L[1]),dp_hc(L[1])*T,PS);
! 577: } else
! 578: NF = nf(GI,T,T,PS);
! 579: NF = remove_cont(NF);
! 580: H = cons(NF,H);
! 581: LCM = ilcm(LCM,dp_hc(NF[1]));
! 582: }
! 583: TNF = time()[0]-T0;
! 584: if ( dp_gr_print() )
! 585: print("gennf(TAB="+rtostr(TTAB)+" NF="+rtostr(TNF)+")");
! 586: return [[map(adj_dn,H,LCM),LCM],PS,GI];
! 587: }
! 588:
! 589: def gennf_demand(G,TL,V,O,V0,FLAG,Dir)
! 590: {
! 591: N = length(V); Len = length(G); dp_ord(O); PS = newvect(Len);
! 592: NTL = length(TL);
! 593: for ( I = 0, T = G, HL = []; T != []; T = cdr(T), I++ ) {
! 594: PS[I] = dp_ptod(car(T),V); HL = cons(dp_ht(PS[I]),HL);
! 595: }
! 596: for ( I = 0, DTL = []; TL != []; TL = cdr(TL) )
! 597: DTL = cons(dp_ptod(car(TL),V),DTL);
! 598: for ( I = Len - 1, GI = []; I >= 0; I-- )
! 599: GI = cons(I,GI);
! 600:
! 601: USE_TAB = (FLAG != 0);
! 602: if ( USE_TAB ) {
! 603: T0 = time()[0];
! 604: MB = dp_mbase(HL); DIM = length(MB);
! 605: U = dp_ptod(V0,V);
! 606: UTAB = newvect(DIM);
! 607: for ( I = 0; I < DIM; I++ ) {
! 608: UTAB[I] = [MB[I],remove_cont(dp_true_nf(GI,U*MB[I],PS,1))];
! 609: if ( dp_gr_print() )
! 610: print(".",2);
! 611: }
! 612: if ( dp_gr_print() )
! 613: print("");
! 614: TTAB = time()[0]-T0;
! 615: }
! 616:
! 617: T0 = time()[0];
! 618: for ( LCM = 1, Index = 0, H = []; DTL != []; Index++ ) {
! 619: if ( dp_gr_print() )
! 620: print(".",2);
! 621: T = car(DTL); DTL = cdr(DTL);
! 622: if ( L = search_redble(T,H) ) {
! 623: L = nf_load(Dir,L[0]);
! 624: DD = dp_subd(T,L[1]);
! 625: if ( USE_TAB && (DD == U) ) {
! 626: NF = nf_tab(L[0],UTAB);
! 627: NF = [NF[0],dp_hc(L[1])*NF[1]*T];
! 628: } else
! 629: NF = nf(GI,L[0]*dp_subd(T,L[1]),dp_hc(L[1])*T,PS);
! 630: } else
! 631: NF = nf(GI,T,T,PS);
! 632: NF = remove_cont(NF);
! 633: nf_save(NF,Dir,Index);
! 634: H = cons([Index,NF[1]],H);
! 635: LCM = ilcm(LCM,dp_hc(NF[1]));
! 636: }
! 637: TNF = time()[0]-T0;
! 638: if ( dp_gr_print() )
! 639: print("gennf(TAB="+rtostr(TTAB)+" NF="+rtostr(TNF)+")");
! 640:
! 641: for ( I = 0; I < NTL; I++ ) {
! 642: NF = nf_load(Dir,I);
! 643: NF = adj_dn(NF,LCM);
! 644: nf_save(NF,Dir,I);
! 645: }
! 646: for ( H = [], I = NTL-1; I >= 0; I-- )
! 647: H = cons(nf_load(Dir,I),H);
! 648: return [[H,LCM],PS,GI];
! 649: }
! 650:
! 651: def nf_load(Dir,I)
! 652: {
! 653: return bload(Dir+"/nf"+rtostr(I));
! 654: }
! 655:
! 656: def nf_save(NF,Dir,I)
! 657: {
! 658: bsave(NF,Dir+"/nf"+rtostr(I));
! 659: }
! 660:
! 661: def adj_dn(P,D)
! 662: {
! 663: return [(idiv(D,dp_hc(P[1])))*P[0],dp_ht(P[1])];
! 664: }
! 665:
! 666: def hen_ttob(T,NF,LHS,V,MOD)
! 667: {
! 668: if ( length(T) == 1 )
! 669: return car(T);
! 670: T0 = time()[0]; M = etom(leq_nf(T,NF,LHS,V)); TE = time()[0] - T0;
! 671: T0 = time()[0]; U = henleq(M,MOD); TH = time()[0] - T0;
! 672: if ( dp_gr_print() ) {
! 673: print("(etom="+rtostr(TE)+" hen="+rtostr(TH)+")");
! 674: }
! 675: return U ? vtop(T,U,LHS) : 0;
! 676: }
! 677:
! 678: def vtop(S,L,GSL)
! 679: {
! 680: U = L[0]; H = L[1];
! 681: if ( GSL ) {
! 682: for ( A = 0, I = 0; S != []; S = cdr(S), I++ )
! 683: A += U[I]*car(S);
! 684: return [A,H];
! 685: } else {
! 686: for ( A = H*car(S), S = cdr(S), I = 0; S != []; S = cdr(S), I++ )
! 687: A += U[I]*car(S);
! 688: return ptozp(A);
! 689: }
! 690: }
! 691:
! 692: /* broken */
! 693:
! 694: def leq_nf(TL,NF,LHS,V)
! 695: {
! 696: TLen = length(NF);
! 697: T = newvect(TLen); M = newvect(TLen);
! 698: for ( I = 0; I < TLen; I++ ) {
! 699: T[I] = dp_ht(NF[I][1]);
! 700: M[I] = dp_hc(NF[I][1]);
! 701: }
! 702: Len = length(TL); INDEX = newvect(Len); COEF = newvect(Len);
! 703: for ( L = TL, J = 0; L != []; L = cdr(L), J++ ) {
! 704: D = dp_ptod(car(L),V);
! 705: for ( I = 0; I < TLen; I++ )
! 706: if ( D == T[I] )
! 707: break;
! 708: INDEX[J] = I; COEF[J] = strtov("u"+rtostr(J));
! 709: }
! 710: if ( !LHS ) {
! 711: COEF[0] = 1; NM = 0; DN = 1;
! 712: } else {
! 713: NM = LHS[0]; DN = LHS[1];
! 714: }
! 715: for ( J = 0, S = -NM; J < Len; J++ ) {
! 716: DNJ = M[INDEX[J]];
! 717: GCD = igcd(DN,DNJ); CS = DNJ/GCD; CJ = DN/GCD;
! 718: S = CS*S + CJ*NF[INDEX[J]][0]*COEF[J];
! 719: DN *= CS;
! 720: }
! 721: for ( D = S, E = []; D; D = dp_rest(D) )
! 722: E = cons(dp_hc(D),E);
! 723: BOUND = LHS ? 0 : 1;
! 724: for ( I = Len - 1, W = []; I >= BOUND; I-- )
! 725: W = cons(COEF[I],W);
! 726: return [E,W];
! 727: }
! 728:
! 729: def nf_tab(F,TAB)
! 730: {
! 731: for ( NM = 0, DN = 1, I = 0; F; F = dp_rest(F) ) {
! 732: T = dp_ht(F);
! 733: for ( ; TAB[I][0] != T; I++);
! 734: NF = TAB[I][1]; N = NF[0]; D = NF[1];
! 735: G = igcd(DN,D); DN1 = idiv(DN,G); D1 = idiv(D,G);
! 736: NM = D1*NM + DN1*dp_hc(F)*N; DN *= D1;
! 737: }
! 738: return [NM,DN];
! 739: }
! 740:
! 741: def nf_tab_gsl(A,NF)
! 742: {
! 743: DN = NF[1];
! 744: NF = NF[0];
! 745: TLen = length(NF);
! 746: for ( R = 0; A; A = dp_rest(A) ) {
! 747: HM = dp_hm(A); C = dp_hc(HM); T = dp_ht(HM);
! 748: for ( I = 0; I < TLen; I++ )
! 749: if ( NF[I][1] == T )
! 750: break;
! 751: R += C*NF[I][0];
! 752: }
! 753: return remove_cont([R,DN]);
! 754: }
! 755:
! 756: def redble(D1,D2,N)
! 757: {
! 758: for ( I = 0; I < N; I++ )
! 759: if ( D1[I] > D2[I] )
! 760: break;
! 761: return I == N ? 1 : 0;
! 762: }
! 763:
! 764: def tolexm(G,V,O,W,M)
! 765: {
! 766: N = length(V); Len = length(G);
! 767: dp_ord(O); setmod(M); PS = newvect(Len);
! 768: for ( I = 0, T = G; T != []; T = cdr(T), I++ )
! 769: PS[I] = dp_mod(dp_ptod(car(T),V),M,[]);
! 770: for ( I = Len-1, HL = []; I >= 0; I-- )
! 771: HL = cons(dp_ht(PS[I]),HL);
! 772: G2 = tolexm_main(PS,HL,V,W,M,ZeroDim);
! 773: L = map(dp_dtop,G2,V);
! 774: return L;
! 775: }
! 776:
! 777: def tolexm_main(PS,HL,V,W,M,FLAG)
! 778: {
! 779: N = length(W); D = newvect(N); Len = size(PS)[0];
! 780: for ( I = Len - 1, GI = []; I >= 0; I-- )
! 781: GI = cons(I,GI);
! 782: MB = dp_mbase(HL); DIM = length(MB);
! 783: U = dp_mod(dp_ptod(W[N-1],V),M,[]);
! 784: UTAB = newvect(DIM);
! 785: for ( I = 0; I < DIM; I++ ) {
! 786: if ( dp_gr_print() )
! 787: print(".",2);
! 788: UTAB[I] = [MB[I],dp_nf_mod(GI,U*dp_mod(MB[I],M,[]),PS,1,M)];
! 789: }
! 790: if ( dp_gr_print() )
! 791: print("");
! 792: T = dp_mod(dp_ptod(dp_dtop(dp_vtoe(D),W),V),M,[]);
! 793: H = G = [[T,T]];
! 794: DL = []; G2 = [];
! 795: TNF = 0;
! 796: while ( 1 ) {
! 797: if ( dp_gr_print() )
! 798: print(".",2);
! 799: S = nextm(D,DL,N);
! 800: if ( !S )
! 801: break;
! 802: T = dp_mod(dp_ptod(dp_dtop(dp_vtoe(D),W),V),M,[]);
! 803: T0 = time()[0];
! 804: if ( L = search_redble(T,H) ) {
! 805: DD = dp_mod(dp_subd(T,L[1]),M,[]);
! 806: if ( DD == U )
! 807: NT = dp_nf_tab_mod(L[0],UTAB,M);
! 808: else
! 809: NT = dp_nf_mod(GI,L[0]*DD,PS,1,M);
! 810: } else
! 811: NT = dp_nf_mod(GI,T,PS,1,M);
! 812: TNF += time()[0] - T0;
! 813: H = cons([NT,T],H);
! 814: T0 = time()[0];
! 815: L = dp_lnf_mod([NT,T],G,M); N1 = L[0]; N2 = L[1];
! 816: TLNF += time()[0] - T0;
! 817: if ( !N1 ) {
! 818: G2 = cons(N2,G2);
! 819: if ( FLAG == MiniPoly )
! 820: break;
! 821: D1 = newvect(N);
! 822: for ( I = 0; I < N; I++ )
! 823: D1[I] = D[I];
! 824: DL = cons(D1,DL);
! 825: } else
! 826: G = insert(G,L);
! 827: }
! 828: if ( dp_gr_print() )
! 829: print("tolexm(nfm="+rtostr(TNF)+" lnfm="+rtostr(TLNF)+")");
! 830: return G2;
! 831: }
! 832:
! 833: def minipolym(G,V,O,P,V0,M)
! 834: {
! 835: N = length(V); Len = length(G);
! 836: dp_ord(O); setmod(M); PS = newvect(Len);
! 837: for ( I = 0, T = G; T != []; T = cdr(T), I++ )
! 838: PS[I] = dp_mod(dp_ptod(car(T),V),M,[]);
! 839: for ( I = Len-1, HL = []; I >= 0; I-- )
! 840: HL = cons(dp_ht(PS[I]),HL);
! 841: for ( I = Len - 1, GI = []; I >= 0; I-- )
! 842: GI = cons(I,GI);
! 843: MB = dp_mbase(HL); DIM = length(MB); UT = newvect(DIM);
! 844: U = dp_mod(dp_ptod(P,V),M,[]);
! 845: for ( I = 0; I < DIM; I++ )
! 846: UT[I] = [MB[I],dp_nf_mod(GI,U*dp_mod(MB[I],M,[]),PS,1,M)];
! 847: T = dp_mod(<<0>>,M,[]); TT = dp_mod(dp_ptod(1,V),M,[]);
! 848: G = H = [[TT,T]]; TNF = TLNF = 0;
! 849: for ( I = 1; ; I++ ) {
! 850: T = dp_mod(<<I>>,M,[]);
! 851: T0 = time()[0]; NT = dp_nf_tab_mod(H[0][0],UT,M); TNF += time()[0] - T0;
! 852: H = cons([NT,T],H);
! 853: T0 = time()[0]; L = dp_lnf_mod([NT,T],G,M); TLNF += time()[0] - T0;
! 854: if ( !L[0] ) {
! 855: if ( dp_gr_print() ) print(["nfm",TNF,"lnfm",TLNF]);
! 856: return dp_dtop(L[1],[V0]);
! 857: } else
! 858: G = insert(G,L);
! 859: }
! 860: }
! 861:
! 862: def nextm(D,DL,N)
! 863: {
! 864: for ( I = N-1; I >= 0; ) {
! 865: D[I]++;
! 866: for ( T = DL; T != []; T = cdr(T) )
! 867: if ( car(T) == D )
! 868: return 1;
! 869: else if ( redble(car(T),D,N) )
! 870: break;
! 871: if ( T != [] ) {
! 872: for ( J = N-1; J >= I; J-- )
! 873: D[J] = 0;
! 874: I--;
! 875: } else
! 876: break;
! 877: }
! 878: if ( I < 0 )
! 879: return 0;
! 880: else
! 881: return 1;
! 882: }
! 883:
! 884: def search_redble(T,G)
! 885: {
! 886: for ( ; G != []; G = cdr(G) )
! 887: if ( dp_redble(T,car(G)[1]) )
! 888: return car(G);
! 889: return 0;
! 890: }
! 891:
! 892: def insert(G,A)
! 893: {
! 894: if ( G == [] )
! 895: return [A];
! 896: else if ( dp_ht(car(A)) > dp_ht(car(car(G))) )
! 897: return cons(A,G);
! 898: else
! 899: return cons(car(G),insert(cdr(G),A));
! 900: }
! 901:
! 902: #if 0
! 903: def etom(L) {
! 904: E = L[0]; W = L[1];
! 905: LE = length(E); LW = length(W);
! 906: M = newmat(LE,LW+1);
! 907: for(J=0;J<LE;J++) {
! 908: for ( T = E[J]; T && (type(T) == 2); )
! 909: for ( V = var(T), I = 0; I < LW; I++ )
! 910: if ( V == W[I] ) {
! 911: M[J][I] = coef(T,1,V);
! 912: T = coef(T,0,V);
! 913: }
! 914: M[J][LW] = T;
! 915: }
! 916: return M;
! 917: }
! 918: #endif
! 919:
! 920: def etom(L) {
! 921: E = L[0]; W = L[1];
! 922: LE = length(E); LW = length(W);
! 923: M = newmat(LE,LW+1);
! 924: for(J=0;J<LE;J++) {
! 925: for ( I = 0, T = E[J]; I < LW; I++ ) {
! 926: M[J][I] = coef(T,1,W[I]); T = coef(T,0,W[I]);
! 927: }
! 928: M[J][LW] = T;
! 929: }
! 930: return M;
! 931: }
! 932:
! 933: def calcb_old(M) {
! 934: N = 2*M;
! 935: T = gr_sqrt(N);
! 936: if ( T^2 <= N && N < (T+1)^2 )
! 937: return idiv(T,2);
! 938: else
! 939: error("afo");
! 940: }
! 941:
! 942: def calcb_special(PK,P,K) { /* PK = P^K */
! 943: N = 2*PK;
! 944: T = sqrt_special(N,2,P,K);
! 945: if ( T^2 <= N && N < (T+1)^2 )
! 946: return idiv(T,2);
! 947: else
! 948: error("afo");
! 949: }
! 950:
! 951: def sqrt_special(A,C,M,K) { /* A = C*M^K */
! 952: L = idiv(K,2); B = M^L;
! 953: if ( K % 2 )
! 954: C *= M;
! 955: D = 2^K; X = idiv((gr_sqrt(C*D^2)+1)*B,D)+1;
! 956: while ( 1 )
! 957: if ( (Y = X^2) <= A )
! 958: return X;
! 959: else
! 960: X = idiv(A + Y,2*X);
! 961: }
! 962:
! 963: def gr_sqrt(A) {
! 964: for ( J = 0, T = A; T >= 2^27; J++ ) {
! 965: T = idiv(T,2^27)+1;
! 966: }
! 967: for ( I = 0; T >= 2; I++ ) {
! 968: S = idiv(T,2);
! 969: if ( T = S+S )
! 970: T = S;
! 971: else
! 972: T = S+1;
! 973: }
! 974: X = (2^27)^idiv(J,2)*2^idiv(I,2);
! 975: while ( 1 ) {
! 976: if ( (Y=X^2) < A )
! 977: X += X;
! 978: else if ( Y == A )
! 979: return X;
! 980: else
! 981: break;
! 982: }
! 983: while ( 1 )
! 984: if ( (Y = X^2) <= A )
! 985: return X;
! 986: else
! 987: X = idiv(A + Y,2*X);
! 988: }
! 989:
! 990: #define ABS(a) ((a)>=0?(a):(-a))
! 991:
! 992: def inttorat_asir(C,M,B)
! 993: {
! 994: if ( M < 0 )
! 995: M = -M;
! 996: C %= M;
! 997: if ( C < 0 )
! 998: C += M;
! 999: U1 = 0; U2 = M; V1 = 1; V2 = C;
! 1000: while ( V2 >= B ) {
! 1001: L = iqr(U2,V2); Q = L[0]; R2 = L[1];
! 1002: R1 = U1 - Q*V1;
! 1003: U1 = V1; U2 = V2;
! 1004: V1 = R1; V2 = R2;
! 1005: }
! 1006: if ( ABS(V1) >= B )
! 1007: return 0;
! 1008: else
! 1009: if ( V1 < 0 )
! 1010: return [-V2,-V1];
! 1011: else
! 1012: return [V2,V1];
! 1013: }
! 1014:
! 1015: def intvtoratv(V,M,B) {
! 1016: if ( !B )
! 1017: B = 1;
! 1018: N = size(V)[0];
! 1019: W = newvect(N);
! 1020: if ( ITOR_FAIL >= 0 ) {
! 1021: if ( V[ITOR_FAIL] ) {
! 1022: T = inttorat(V[ITOR_FAIL],M,B);
! 1023: if ( !T ) {
! 1024: if ( dp_gr_print() ) {
! 1025: print("F",2);
! 1026: }
! 1027: return 0;
! 1028: }
! 1029: }
! 1030: }
! 1031: for ( I = 0, DN = 1; I < N; I++ )
! 1032: if ( V[I] ) {
! 1033: T = inttorat((V[I]*DN) % M,M,B);
! 1034: if ( !T ) {
! 1035: ITOR_FAIL = I;
! 1036: if ( dp_gr_print() ) {
! 1037: #if 0
! 1038: print("intvtoratv : failed at I = ",0); print(ITOR_FAIL);
! 1039: #endif
! 1040: print("F",2);
! 1041: }
! 1042: return 0;
! 1043: } else {
! 1044: for( J = 0; J < I; J++ )
! 1045: W[J] *= T[1];
! 1046: W[I] = T[0]; DN *= T[1];
! 1047: }
! 1048: }
! 1049: return [W,DN];
! 1050: }
! 1051:
! 1052: def nf(B,G,M,PS)
! 1053: {
! 1054: for ( D = 0; G; ) {
! 1055: for ( U = 0, L = B; L != []; L = cdr(L) ) {
! 1056: if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
! 1057: GCD = igcd(dp_hc(G),dp_hc(R));
! 1058: CG = idiv(dp_hc(R),GCD); CR = idiv(dp_hc(G),GCD);
! 1059: U = CG*G-dp_subd(G,R)*CR*R;
! 1060: if ( !U )
! 1061: return [D,M];
! 1062: D *= CG; M *= CG;
! 1063: break;
! 1064: }
! 1065: }
! 1066: if ( U )
! 1067: G = U;
! 1068: else {
! 1069: D += dp_hm(G); G = dp_rest(G);
! 1070: }
! 1071: }
! 1072: return [D,M];
! 1073: }
! 1074:
! 1075: def remove_cont(L)
! 1076: {
! 1077: if ( type(L[1]) == 1 ) {
! 1078: T = remove_cont([L[0],L[1]*<<0>>]);
! 1079: return [T[0],dp_hc(T[1])];
! 1080: } else if ( !L[0] )
! 1081: return [0,dp_ptozp(L[1])];
! 1082: else if ( !L[1] )
! 1083: return [dp_ptozp(L[0]),0];
! 1084: else {
! 1085: A0 = dp_ptozp(L[0]); A1 = dp_ptozp(L[1]);
! 1086: C0 = idiv(dp_hc(L[0]),dp_hc(A0)); C1 = idiv(dp_hc(L[1]),dp_hc(A1));
! 1087: GCD = igcd(C0,C1); M0 = idiv(C0,GCD); M1 = idiv(C1,GCD);
! 1088: return [M0*A0,M1*A1];
! 1089: }
! 1090: }
! 1091:
! 1092: def union(A,B)
! 1093: {
! 1094: for ( T = B; T != []; T = cdr(T) )
! 1095: A = union1(A,car(T));
! 1096: return A;
! 1097: }
! 1098:
! 1099: def union1(A,E)
! 1100: {
! 1101: if ( A == [] )
! 1102: return [E];
! 1103: else if ( car(A) == E )
! 1104: return A;
! 1105: else
! 1106: return cons(car(A),union1(cdr(A),E));
! 1107: }
! 1108:
! 1109: def setminus(A,B) {
! 1110: for ( T = reverse(A), R = []; T != []; T = cdr(T) ) {
! 1111: for ( S = B, M = car(T); S != []; S = cdr(S) )
! 1112: if ( car(S) == M )
! 1113: break;
! 1114: if ( S == [] )
! 1115: R = cons(M,R);
! 1116: }
! 1117: return R;
! 1118: }
! 1119:
! 1120: def member(A,L) {
! 1121: for ( ; L != []; L = cdr(L) )
! 1122: if ( A == car(L) )
! 1123: return 1;
! 1124: return 0;
! 1125: }
! 1126:
! 1127: /* several functions for computation of normal forms etc. */
! 1128:
! 1129: def p_nf(P,B,V,O) {
! 1130: dp_ord(O); DP = dp_ptod(P,V);
! 1131: N = length(B); DB = newvect(N);
! 1132: for ( I = N-1, IL = []; I >= 0; I-- ) {
! 1133: DB[I] = dp_ptod(B[I],V);
! 1134: if ( DB[I] ) IL = cons(I,IL);
! 1135: }
! 1136: return dp_dtop(dp_nf(IL,DP,DB,1),V);
! 1137: }
! 1138:
! 1139: def p_true_nf(P,B,V,O) {
! 1140: dp_ord(O); DP = dp_ptod(P,V);
! 1141: N = length(B); DB = newvect(N);
! 1142: for ( I = N-1, IL = []; I >= 0; I-- ) {
! 1143: DB[I] = dp_ptod(B[I],V);
! 1144: if ( DB[I] ) IL = cons(I,IL);
! 1145: }
! 1146: L = dp_true_nf(IL,DP,DB,1);
! 1147: return [dp_dtop(L[0],V),L[1]];
! 1148: }
! 1149:
! 1150: def p_nf_mod(P,B,V,O,Mod) {
! 1151: setmod(Mod);
! 1152: dp_ord(O); DP = dp_mod(dp_ptod(P,V),Mod,[]);
! 1153: N = length(B); DB = newvect(N);
! 1154: for ( I = N-1, IL = []; I >= 0; I-- ) {
! 1155: DB[I] = dp_mod(dp_ptod(B[I],V),Mod,[]);
! 1156: IL = cons(I,IL);
! 1157: }
! 1158: return dp_dtop(dp_nf_mod(IL,DP,DB,1,Mod),V);
! 1159: }
! 1160:
! 1161: def p_terms(D,V,O)
! 1162: {
! 1163: dp_ord(O);
! 1164: for ( L = [], T = dp_ptod(D,V); T; T = dp_rest(T) )
! 1165: L = cons(dp_dtop(dp_ht(T),V),L);
! 1166: return reverse(L);
! 1167: }
! 1168:
! 1169: def dp_terms(D,V)
! 1170: {
! 1171: for ( L = [], T = D; T; T = dp_rest(T) )
! 1172: L = cons(dp_dtop(dp_ht(T),V),L);
! 1173: return reverse(L);
! 1174: }
! 1175:
! 1176: def gb_comp(A,B)
! 1177: {
! 1178: LA = length(A);
! 1179: LB = length(B);
! 1180: if ( LA != LB )
! 1181: return 0;
! 1182: A = newvect(LA,A);
! 1183: B = newvect(LB,B);
! 1184: for ( I = 0; I < LA; I++ )
! 1185: A[I] *= headsgn(A[I]);
! 1186: for ( I = 0; I < LB; I++ )
! 1187: B[I] *= headsgn(B[I]);
! 1188: A1 = qsort(A);
! 1189: B1 = qsort(B);
! 1190: for ( I = 0; I < LA; I++ )
! 1191: if ( A1[I] != B1[I] && A1[I] != -B1[I] )
! 1192: break;
! 1193: return I == LA ? 1 : 0;
! 1194: }
! 1195:
! 1196: def zero_dim(G,V,O) {
! 1197: dp_ord(O);
! 1198: HL = map(dp_dtop,map(dp_ht,map(dp_ptod,G,V)),V);
! 1199: for ( L = []; HL != []; HL = cdr(HL) )
! 1200: if ( length(vars(car(HL))) == 1 )
! 1201: L = cons(car(HL),L);
! 1202: return length(vars(L)) == length(V) ? 1 : 0;
! 1203: }
! 1204:
! 1205: def hmlist(G,V,O) {
! 1206: dp_ord(O);
! 1207: return map(dp_dtop,map(dp_hm,map(dp_ptod,G,V)),V);
! 1208: }
! 1209:
! 1210: def valid_modulus(HL,M) {
! 1211: V = vars(HL);
! 1212: for ( T = HL; T != []; T = cdr(T) )
! 1213: if ( !dp_mod(dp_ptod(car(T),V),M,[]) )
! 1214: break;
! 1215: return T == [] ? 1 : 0;
! 1216: }
! 1217:
! 1218: def npos_check(DL) {
! 1219: N = size(car(DL))[0];
! 1220: if ( length(DL) != N )
! 1221: return [-1,0];
! 1222: D = newvect(N);
! 1223: for ( I = 0; I < N; I++ ) {
! 1224: for ( J = 0; J < N; J++ )
! 1225: D[J] = 0;
! 1226: D[I] = 1;
! 1227: for ( T = DL; T != []; T = cdr(T) )
! 1228: if ( D == car(T) )
! 1229: break;
! 1230: if ( T != [] )
! 1231: DL = setminus(DL,[car(T)]);
! 1232: }
! 1233: if ( length(DL) != 1 )
! 1234: return [-1,0];
! 1235: U = car(DL);
! 1236: for ( I = 0, J = 0, I0 = -1; I < N; I++ )
! 1237: if ( U[I] ) {
! 1238: I0 = I; J++;
! 1239: }
! 1240: if ( J != 1 )
! 1241: return [-1,0];
! 1242: else
! 1243: return [I0,U[I0]];
! 1244: }
! 1245:
! 1246: def mult_mat(L,TAB,MB)
! 1247: {
! 1248: A = L[0]; DN0 = L[1];
! 1249: for ( NM = 0, DN = 1, I = 0; A; A = dp_rest(A) ) {
! 1250: H = dp_ht(A);
! 1251: for ( ; MB[I] != H; I++ );
! 1252: NM1 = TAB[I][0]; DN1 = TAB[I][1]; I++;
! 1253: GCD = igcd(DN,DN1); C = DN1/GCD; C1 = DN/GCD;
! 1254: NM = C*NM + C1*dp_hc(A)*NM1;
! 1255: DN *= C;
! 1256: }
! 1257: Z=remove_cont([NM,DN*DN0]);
! 1258: return Z;
! 1259: }
! 1260:
! 1261: def sepm(MAT)
! 1262: {
! 1263: S = size(MAT); N = S[0]; M = S[1]-1;
! 1264: A = newmat(N,M); B = newvect(N);
! 1265: for ( I = 0; I < N; I++ )
! 1266: for ( J = 0, T1 = MAT[I], T2 = A[I]; J < M; J++ )
! 1267: T2[J] = T1[J];
! 1268: for ( I = 0; I < N; I++ )
! 1269: B[I] = MAT[I][M];
! 1270: return [A,B];
! 1271: }
! 1272:
! 1273: def henleq(M,MOD)
! 1274: {
! 1275: SIZE = size(M); ROW = SIZE[0]; COL = SIZE[1];
! 1276: W = newvect(COL);
! 1277: L = sepm(M); A = L[0]; B = L[1];
! 1278: COUNT = INIT_COUNT?INIT_COUNT:idiv(max_mag(M),54);
! 1279: if ( !COUNT )
! 1280: COUNT = 1;
! 1281:
! 1282: TINV = TC = TR = TS = TM = TDIV = 0;
! 1283:
! 1284: T0 = time()[0];
! 1285: L = geninvm_swap(A,MOD); INV = L[0]; INDEX = L[1];
! 1286: TS += time()[0] - T0;
! 1287:
! 1288: COL1 = COL - 1;
! 1289: AA = newmat(COL1,COL1); BB = newvect(COL1);
! 1290: for ( I = 0; I < COL1; I++ ) {
! 1291: for ( J = 0, T = AA[I], S = A[INDEX[I]]; J < COL1; J++ )
! 1292: T[J] = S[J];
! 1293: BB[I] = B[INDEX[I]];
! 1294: }
! 1295: if ( COL1 != ROW ) {
! 1296: RESTA = newmat(ROW-COL1,COL1); RESTB = newvect(ROW-COL1);
! 1297: for ( ; I < ROW; I++ ) {
! 1298: for ( J = 0, T = RESTA[I-COL1], S = A[INDEX[I]]; J < COL1; J++ )
! 1299: T[J] = S[J];
! 1300: RESTB[I-COL1] = B[INDEX[I]];
! 1301: }
! 1302: } else
! 1303: RESTA = RESTB = 0;
! 1304:
! 1305: MOD2 = idiv(MOD,2);
! 1306: for ( I = 0, C = BB, X = 0, PK = 1, CCC = 0, ITOR_FAIL = -1; ;
! 1307: I++, PK *= MOD ) {
! 1308: if ( COUNT == CCC ) {
! 1309: CCC = 0;
! 1310: T0 = time()[0];
! 1311: ND = intvtoratv(X,PK,ishift(calcb_special(PK,MOD,I),32));
! 1312: TR += time()[0]-T0;
! 1313: if ( ND ) {
! 1314: T0 = time()[0];
! 1315: F = ND[0]; LCM = ND[1]; T = AA*F+LCM*BB;
! 1316: TM += time()[0]-T0;
! 1317: if ( zerovector(T) ) {
! 1318: T0 = time()[0]; T = RESTA*F+LCM*RESTB; TM += time()[0]-T0;
! 1319: if ( zerovector(T) ) {
! 1320: #if 0
! 1321: if ( dp_gr_print() ) print(["init",TS,"pinv",TINV,"c",TC,"div",TDIV,"rat",TR,"mul",TM]);
! 1322: #endif
! 1323: if ( dp_gr_print() ) print("end",2);
! 1324: return [F,LCM];
! 1325: } else
! 1326: return 0;
! 1327: }
! 1328: } else {
! 1329: #if 0
! 1330: if ( dp_gr_print() ) print(I);
! 1331: #endif
! 1332: }
! 1333: } else {
! 1334: #if 0
! 1335: if ( dp_gr_print() ) print([I,TINV,TC,TDIV]);
! 1336: #endif
! 1337: if ( dp_gr_print() ) print(".",2);
! 1338: CCC++;
! 1339: }
! 1340: T0 = time()[0];
! 1341: XT = sremainder(INV*sremainder(-C,MOD),MOD);
! 1342: XT = map(adj_sgn,XT,MOD,MOD2);
! 1343: TINV += time()[0] - T0;
! 1344: X += XT*PK;
! 1345: T0 = time()[0];
! 1346: C += mul_mat_vect_int(AA,XT);
! 1347: TC += time()[0] - T0;
! 1348: T0 = time()[0]; C = map(idiv,C,MOD); TDIV += time()[0] - T0;
! 1349: }
! 1350: }
! 1351:
! 1352: def henleq_prep(A,MOD)
! 1353: {
! 1354: SIZE = size(A); ROW = SIZE[0]; COL = SIZE[1];
! 1355: L = geninvm_swap(A,MOD); INV = L[0]; INDEX = L[1];
! 1356: AA = newmat(COL,COL);
! 1357: for ( I = 0; I < COL; I++ )
! 1358: for ( J = 0, T = AA[I], S = A[INDEX[I]]; J < COL; J++ )
! 1359: T[J] = S[J];
! 1360: if ( COL != ROW ) {
! 1361: RESTA = newmat(ROW-COL,COL);
! 1362: for ( ; I < ROW; I++ )
! 1363: for ( J = 0, T = RESTA[I-COL], S = A[INDEX[I]]; J < COL; J++ )
! 1364: T[J] = S[J];
! 1365: } else
! 1366: RESTA = 0;
! 1367: return [[A,AA,RESTA],L];
! 1368: }
! 1369:
! 1370: def henleq_gsl(L,B,MOD)
! 1371: {
! 1372: AL = L[0]; INVL = L[1];
! 1373: A = AL[0]; AA = AL[1]; RESTA = AL[2];
! 1374: INV = INVL[0]; INDEX = INVL[1];
! 1375: SIZE = size(A); ROW = SIZE[0]; COL = SIZE[1];
! 1376: BB = newvect(COL);
! 1377: for ( I = 0; I < COL; I++ )
! 1378: BB[I] = B[INDEX[I]];
! 1379: if ( COL != ROW ) {
! 1380: RESTB = newvect(ROW-COL);
! 1381: for ( ; I < ROW; I++ )
! 1382: RESTB[I-COL] = B[INDEX[I]];
! 1383: } else
! 1384: RESTB = 0;
! 1385:
! 1386: COUNT = INIT_COUNT?INIT_COUNT:idiv(MAX(max_mag(A),max_mag_vect(B)),54);
! 1387: if ( !COUNT )
! 1388: COUNT = 1;
! 1389: MOD2 = idiv(MOD,2);
! 1390: X = newvect(size(AA)[0]);
! 1391: for ( I = 0, C = BB, PK = 1, CCC = 0, ITOR_FAIL = -1; ;
! 1392: I++, PK *= MOD ) {
! 1393: if ( zerovector(C) )
! 1394: if ( zerovector(RESTA*X+RESTB) ) {
! 1395: if ( dp_gr_print() ) print("end",0);
! 1396: return [X,1];
! 1397: } else
! 1398: return 0;
! 1399: else if ( COUNT == CCC ) {
! 1400: CCC = 0;
! 1401: ND = intvtoratv(X,PK,ishift(calcb_special(PK,MOD,I),32));
! 1402: if ( ND ) {
! 1403: F = ND[0]; LCM = ND[1]; T = AA*F+LCM*BB;
! 1404: if ( zerovector(T) ) {
! 1405: T = RESTA*F+LCM*RESTB;
! 1406: if ( zerovector(T) ) {
! 1407: if ( dp_gr_print() ) print("end",0);
! 1408: return [F,LCM];
! 1409: } else
! 1410: return 0;
! 1411: }
! 1412: } else {
! 1413: }
! 1414: } else {
! 1415: if ( dp_gr_print() ) print(".",2);
! 1416: CCC++;
! 1417: }
! 1418: XT = sremainder(INV*sremainder(-C,MOD),MOD);
! 1419: XT = map(adj_sgn,XT,MOD,MOD2);
! 1420: X += XT*PK;
! 1421: C += mul_mat_vect_int(AA,XT);
! 1422: C = map(idiv,C,MOD);
! 1423: }
! 1424: }
! 1425:
! 1426: def adj_sgn(A,M,M2)
! 1427: {
! 1428: return A > M2 ? A-M : A;
! 1429: }
! 1430:
! 1431: def zerovector(C)
! 1432: {
! 1433: if ( !C )
! 1434: return 1;
! 1435: for ( I = size(C)[0]-1; I >= 0 && !C[I]; I-- );
! 1436: if ( I < 0 )
! 1437: return 1;
! 1438: else
! 1439: return 0;
! 1440: }
! 1441:
! 1442: def solvem(INV,COMP,V,MOD)
! 1443: {
! 1444: T = COMP*V;
! 1445: N = size(T)[0];
! 1446: for ( I = 0; I < N; I++ )
! 1447: if ( T[I] % MOD )
! 1448: return 0;
! 1449: return modvect(INV*V,MOD);
! 1450: }
! 1451:
! 1452: def modmat(A,MOD)
! 1453: {
! 1454: if ( !A )
! 1455: return 0;
! 1456: S = size(A); N = S[0]; M = S[1];
! 1457: MAT = newmat(N,M);
! 1458: for ( I = 0, NZ = 0; I < N; I++ )
! 1459: for ( J = 0, T1 = A[I], T2 = MAT[I]; J < M; J++ ) {
! 1460: T2[J] = T1[J] % MOD;
! 1461: NZ = NZ || T2[J];
! 1462: }
! 1463: return NZ?MAT:0;
! 1464: }
! 1465:
! 1466: def modvect(A,MOD)
! 1467: {
! 1468: if ( !A )
! 1469: return 0;
! 1470: N = size(A)[0];
! 1471: VECT = newvect(N);
! 1472: for ( I = 0, NZ = 0; I < N; I++ ) {
! 1473: VECT[I] = A[I] % MOD;
! 1474: NZ = NZ || VECT[I];
! 1475: }
! 1476: return NZ?VECT:0;
! 1477: }
! 1478:
! 1479: def qrmat(A,MOD)
! 1480: {
! 1481: if ( !A )
! 1482: return [0,0];
! 1483: S = size(A); N = S[0]; M = S[1];
! 1484: Q = newmat(N,M); R = newmat(N,M);
! 1485: for ( I = 0, NZQ = 0, NZR = 0; I < N; I++ )
! 1486: for ( J = 0, TA = A[I], TQ = Q[I], TR = R[I]; J < M; J++ ) {
! 1487: L = iqr(TA[J],MOD); TQ[J] = L[0]; TR[J] = L[1];
! 1488: NZQ = NZQ || TQ[J]; NZR = NZR || TR[J];
! 1489: }
! 1490: return [NZQ?Q:0,NZR?R:0];
! 1491: }
! 1492:
! 1493: def qrvect(A,MOD)
! 1494: {
! 1495: if ( !A )
! 1496: return [0,0];
! 1497: N = size(A)[0];
! 1498: Q = newvect(N); R = newvect(N);
! 1499: for ( I = 0, NZQ = 0, NZR = 0; I < N; I++ ) {
! 1500: L = iqr(A[I],MOD); Q[I] = L[0]; R[I] = L[1];
! 1501: NZQ = NZQ || Q[I]; NZR = NZR || R[I];
! 1502: }
! 1503: return [NZQ?Q:0,NZR?R:0];
! 1504: }
! 1505:
! 1506: def max_mag(M)
! 1507: {
! 1508: R = size(M)[0];
! 1509: U = 1;
! 1510: for ( I = 0; I < R; I++ ) {
! 1511: A = max_mag_vect(M[I]);
! 1512: U = MAX(A,U);
! 1513: }
! 1514: return U;
! 1515: }
! 1516:
! 1517: def max_mag_vect(V)
! 1518: {
! 1519: R = size(V)[0];
! 1520: U = 1;
! 1521: for ( I = 0; I < R; I++ ) {
! 1522: A = dp_mag(V[I]*<<0>>);
! 1523: U = MAX(A,U);
! 1524: }
! 1525: return U;
! 1526: }
! 1527:
! 1528: def gsl_check(B,V,S)
! 1529: {
! 1530: N = length(V);
! 1531: U = S[N-1]; M = U[1]; D = U[2];
! 1532: W = setminus(V,[var(M)]);
! 1533: H = uc(); VH = append(W,[H]);
! 1534: for ( T = B; T != []; T = cdr(T) ) {
! 1535: A = car(T);
! 1536: AH = dp_dtop(dp_homo(dp_ptod(A,W)),VH);
! 1537: for ( I = 0, Z = S; I < N-1; I++, Z = cdr(Z) ) {
! 1538: L = car(Z); AH = ptozp(subst(AH,L[0],L[1]/L[2]));
! 1539: }
! 1540: AH = ptozp(subst(AH,H,D));
! 1541: R = srem(AH,M);
! 1542: if ( dp_gr_print() )
! 1543: if ( !R )
! 1544: print([A,"ok"]);
! 1545: else
! 1546: print([A,"bad"]);
! 1547: if ( R )
! 1548: break;
! 1549: }
! 1550: return T == [] ? 1 : 0;
! 1551: }
! 1552:
! 1553: def vs_dim(G,V,O)
! 1554: {
! 1555: HM = hmlist(G,V,O); ZD = zero_dim(HM,V,O);
! 1556: if ( ZD ) {
! 1557: MB = dp_mbase(map(dp_ptod,HM,V));
! 1558: return length(MB);
! 1559: } else
! 1560: error("vs_dim : ideal is not zero-dimensional!");
! 1561: }
! 1562:
! 1563: def dgr(G,V,O)
! 1564: {
! 1565: P = getopt(proc);
! 1566: if ( type(P) == -1 )
! 1567: return gr(G,V,O);
! 1568: P0 = P[0]; P1 = P[1]; P = [P0,P1];
! 1569: map(ox_reset,P);
! 1570: ox_cmo_rpc(P0,"dp_gr_main",G,V,0,1,O);
! 1571: ox_cmo_rpc(P1,"dp_gr_main",G,V,1,1,O);
! 1572: map(ox_push_cmd,P,262); /* 262 = OX_popCMO */
! 1573: F = ox_select(P);
! 1574: R = ox_get(F[0]);
! 1575: if ( F[0] == P0 ) {
! 1576: Win = "nonhomo";
! 1577: Lose = P1;
! 1578: } else {
! 1579: Win = "homo";
! 1580: Lose = P0;
! 1581: }
! 1582: ox_reset(Lose);
! 1583: return [Win,R];
! 1584: }
! 1585:
! 1586: /* competitive Gbase computation : F4 vs. Bucbberger */
! 1587: /* P : process list */
! 1588:
! 1589: def dgrf4mod(G,V,M,O)
! 1590: {
! 1591: P = getopt(proc);
! 1592: if ( type(P) == -1 )
! 1593: return dp_f4_mod_main(G,V,M,O);
! 1594: P0 = P[0]; P1 = P[1]; P = [P0,P1];
! 1595: map(ox_reset,P);
! 1596: ox_cmo_rpc(P0,"dp_f4_mod_main",G,V,M,O);
! 1597: ox_cmo_rpc(P1,"dp_gr_mod_main",G,V,0,M,O);
! 1598: map(ox_push_cmd,P,262); /* 262 = OX_popCMO */
! 1599: F = ox_select(P);
! 1600: R = ox_get(F[0]);
! 1601: if ( F[0] == P0 ) {
! 1602: Win = "F4";
! 1603: Lose = P1;
! 1604: } else {
! 1605: Win = "Buchberger";
! 1606: Lose = P0;
! 1607: }
! 1608: ox_reset(Lose);
! 1609: return [Win,R];
! 1610: }
! 1611:
! 1612: /* functions for rpc */
! 1613:
! 1614: def register_matrix(M)
! 1615: {
! 1616: REMOTE_MATRIX = M; return 0;
! 1617: }
! 1618:
! 1619: def register_nfv(L)
! 1620: {
! 1621: REMOTE_NF = L[0]; REMOTE_VARS = L[1]; return 0;
! 1622: }
! 1623:
! 1624: def r_ttob(T,M)
! 1625: {
! 1626: return hen_ttob(T,REMOTE_NF,0,REMOTE_VARS,M);
! 1627: }
! 1628:
! 1629: def r_ttob_gsl(L,M)
! 1630: {
! 1631: return cons(L[2],hen_ttob(L[0],REMOTE_NF,L[1],REMOTE_VARS,M));
! 1632: }
! 1633:
! 1634: def get_matrix()
! 1635: {
! 1636: REMOTE_MATRIX;
! 1637: }
! 1638:
! 1639: extern NFArray$
! 1640:
! 1641: /*
! 1642: * HL = [[c,i,m,d],...]
! 1643: * if c != 0
! 1644: * g = 0
! 1645: * g = (c*g + m*gi)/d
! 1646: * ...
! 1647: * finally compare g with NF
! 1648: * if g == NF then NFArray[NFIndex] = g
! 1649: *
! 1650: * if c = 0 then HL consists of single history [0,i,0,0],
! 1651: * which means that dehomogenization of NFArray[i] should be
! 1652: * eqall to NF.
! 1653: */
! 1654:
! 1655: def check_trace(NF,NFIndex,HL)
! 1656: {
! 1657: if ( !car(HL)[0] ) {
! 1658: /* dehomogenization */
! 1659: DH = dp_dehomo(NFArray[car(HL)[1]]);
! 1660: if ( NF == DH ) {
! 1661: realloc_NFArray(NFIndex);
! 1662: NFArray[NFIndex] = NF;
! 1663: return 0;
! 1664: } else
! 1665: error("check_trace(dehomo)");
! 1666: }
! 1667:
! 1668: for ( G = 0, T = HL; T != []; T = cdr(T) ) {
! 1669: H = car(T);
! 1670:
! 1671: Coeff = H[0];
! 1672: Index = H[1];
! 1673: Monomial = H[2];
! 1674: Denominator = H[3];
! 1675:
! 1676: Reducer = NFArray[Index];
! 1677: G = (Coeff*G+Monomial*Reducer)/Denominator;
! 1678: }
! 1679: if ( NF == G ) {
! 1680: realloc_NFArray(NFIndex);
! 1681: NFArray[NFIndex] = NF;
! 1682: return 0;
! 1683: } else
! 1684: error("check_trace");
! 1685: }
! 1686:
! 1687: /*
! 1688: * Trace = [Input,[[j1,[[c,i,m,d],...]],[j2,[[...],...]],...]]
! 1689: * if c != 0
! 1690: * g = 0
! 1691: * g = (c*g + m*gi)/d
! 1692: * ...
! 1693: * finally fj = g
! 1694: */
! 1695:
! 1696: def show_trace(Trace,V)
! 1697: {
! 1698: Input = Trace[0];
! 1699: for ( I = 0, T = Input; T != []; T = cdr(T), I++ ) {
! 1700: print("F"+rtostr(I)+"=",0);
! 1701: print(dp_dtop(car(T),V));
! 1702: }
! 1703: Trace = cdr(Trace);
! 1704: for ( T = Trace; T != []; T = cdr(T) ) {
! 1705: HL = car(T);
! 1706: J = car(HL); HL = HL[1];
! 1707: L = length(HL);
! 1708: print("F"+rtostr(J)+"=",0);
! 1709: for ( I = 0; I < L; I++ ) print("(",0);
! 1710: for ( First = 1, S = HL; S != []; S = cdr(S) ) {
! 1711: H = car(S);
! 1712:
! 1713: Coeff = H[0];
! 1714: Index = H[1];
! 1715: Monomial = H[2];
! 1716: Denominator = H[3];
! 1717: if ( First ) {
! 1718: if ( Monomial != 1 ) {
! 1719: print("(",0);
! 1720: print(type(Monomial)==9?dp_dtop(Monomial,V):Monomial,0);
! 1721: print(")*",0);
! 1722: }
! 1723: print("F"+rtostr(Index)+")",0);
! 1724: } else {
! 1725: if ( Coeff != 1 ) {
! 1726: print("*(",0); print(Coeff,0); print(")",0);
! 1727: }
! 1728: print("+",0);
! 1729: if ( Monomial != 1 ) {
! 1730: print("(",0);
! 1731: print(type(Monomial)==9?dp_dtop(Monomial,V):Monomial,0);
! 1732: print(")*",0);
! 1733: }
! 1734: print("F"+rtostr(Index)+")",0);
! 1735: if ( Denominator != 1 ) {
! 1736: print("/",0); print(Denominator,0);
! 1737: }
! 1738: }
! 1739: if ( First ) First = 0;
! 1740: }
! 1741: print("");
! 1742: }
! 1743: }
! 1744:
! 1745: def generating_relation(Trace,V)
! 1746: {
! 1747: Trace = cdr(Trace);
! 1748: Tab = [];
! 1749: for ( T = Trace; T != []; T = cdr(T) ) {
! 1750: HL = car(T);
! 1751: J = car(HL); HL = HL[1];
! 1752: L = length(HL);
! 1753: LHS = strtov("f"+rtostr(J));
! 1754: Dn = 1;
! 1755: for ( First = 1, S = HL; S != []; S = cdr(S) ) {
! 1756: H = car(S);
! 1757:
! 1758: Coeff = H[0];
! 1759: Index = H[1];
! 1760: Monomial = type(H[2])==9?dp_dtop(H[2],V):H[2];
! 1761: Denominator = H[3];
! 1762: F = strtov("f"+rtostr(Index));
! 1763: for ( Z = Tab; Z != []; Z = cdr(Z) )
! 1764: if ( Z[0][0] == F ) break;
! 1765: if ( Z != [] ) Value = Z[0][1];
! 1766: else Value = [F,1];
! 1767: if ( First ) {
! 1768: RHS = Monomial*Value[0];
! 1769: Dn *= Value[1];
! 1770: } else {
! 1771: RHS = RHS*Coeff*Value[1]+Dn*Value[0]*Monomial;
! 1772: Dn = Value[1]*Dn*Denominator;
! 1773: }
! 1774: VVVV = tttttttt;
! 1775: P = ptozp(Dn*VVVV+RHS);
! 1776: RHS = coef(P,0,VVVV);
! 1777: Dn = coef(P,1,VVVV);
! 1778: if ( First ) First = 0;
! 1779: }
! 1780: Tab = cons([LHS,[RHS,Dn]],Tab);
! 1781: }
! 1782: return Tab;
! 1783: }
! 1784:
! 1785: def generating_relation_mod(Trace,V,M)
! 1786: {
! 1787: Trace = cdr(Trace);
! 1788: Tab = [];
! 1789: for ( T = Trace; T != []; T = cdr(T) ) {
! 1790: HL = car(T);
! 1791: J = car(HL); HL = HL[1];
! 1792: L = length(HL);
! 1793: LHS = strtov("f"+rtostr(J));
! 1794: Dn = 1;
! 1795: for ( First = 1, S = HL; S != []; S = cdr(S) ) {
! 1796: H = car(S);
! 1797:
! 1798: Coeff = H[0];
! 1799: Index = H[1];
! 1800: Monomial = type(H[2])==9?dp_dtop(H[2],V):H[2];
! 1801: F = strtov("f"+rtostr(Index));
! 1802: for ( Z = Tab; Z != []; Z = cdr(Z) )
! 1803: if ( Z[0][0] == F ) break;
! 1804: if ( Z != [] ) Value = Z[0][1];
! 1805: else Value = F;
! 1806: if ( First ) {
! 1807: RHS = (Monomial*Value)%M;
! 1808: } else {
! 1809: RHS = ((RHS*Coeff+Value*Monomial)*inv(H[3],M))%M;
! 1810: }
! 1811: if ( First ) First = 0;
! 1812: }
! 1813: Tab = cons([LHS,RHS],Tab);
! 1814: }
! 1815: return Tab;
! 1816: }
! 1817: /*
! 1818: * realloc NFArray so that it can hold * an element as NFArray[Ind].
! 1819: */
! 1820:
! 1821: def realloc_NFArray(Ind)
! 1822: {
! 1823: if ( Ind == size(NFArray)[0] ) {
! 1824: New = newvect(Ind + 100);
! 1825: for ( I = 0; I < Ind; I++ )
! 1826: New[I] = NFArray[I];
! 1827: NFArray = New;
! 1828: }
! 1829: }
! 1830:
! 1831: /*
! 1832: * create NFArray and initialize it by List.
! 1833: */
! 1834:
! 1835: def register_input(List)
! 1836: {
! 1837: Len = length(List);
! 1838: NFArray = newvect(Len+100,List);
! 1839: }
! 1840:
! 1841: /*
! 1842: tracetogen(): preliminary version
! 1843:
! 1844: dp_gr_main() returns [GB,GBIndex,Trace].
! 1845: GB : groebner basis
! 1846: GBIndex : IndexList (corresponding to Trace)
! 1847: Trace : [InputList,Trace0,Trace1,...]
! 1848: TraceI : [Index,TraceList]
! 1849: TraceList : [[Coef,Index,Monomial,Denominator],...]
! 1850: Poly <- 0
! 1851: Poly <- (Coef*Poly+Monomial*PolyList[Index])/Denominator
! 1852: */
! 1853:
! 1854: def tracetogen(G)
! 1855: {
! 1856: GB = G[0]; GBIndex = G[1]; Trace = G[2];
! 1857:
! 1858: InputList = Trace[0];
! 1859: Trace = cdr(Trace);
! 1860:
! 1861: /* number of initial basis */
! 1862: Nini = length(InputList);
! 1863:
! 1864: /* number of generated basis */
! 1865: Ngen = length(Trace);
! 1866:
! 1867: N = Nini + Ngen;
! 1868:
! 1869: /* stores traces */
! 1870: Tr = vector(N);
! 1871:
! 1872: /* stores coeffs */
! 1873: Coef = vector(N);
! 1874:
! 1875: /* XXX create dp_ptod(1,V) */
! 1876: HT = dp_ht(InputList[0]);
! 1877: One = dp_subd(HT,HT);
! 1878:
! 1879: for ( I = 0; I < Nini; I++ ) {
! 1880: Tr[I] = [1,I,One,1];
! 1881: C = vector(Nini);
! 1882: C[I] = One;
! 1883: Coef[I] = C;
! 1884: }
! 1885: for ( ; I < N; I++ )
! 1886: Tr[I] = Trace[I-Nini][1];
! 1887:
! 1888: for ( T = GBIndex; T != []; T = cdr(T) )
! 1889: compute_coef_by_trace(car(T),Tr,Coef);
! 1890: return Coef;
! 1891: }
! 1892:
! 1893: def compute_coef_by_trace(I,Tr,Coef)
! 1894: {
! 1895: if ( Coef[I] )
! 1896: return;
! 1897:
! 1898: /* XXX */
! 1899: Nini = size(Coef[0])[0];
! 1900:
! 1901: /* initialize coef vector */
! 1902: CI = vector(Nini);
! 1903:
! 1904: for ( T = Tr[I]; T != []; T = cdr(T) ) {
! 1905: /* Trace = [Coef,Index,Monomial,Denominator] */
! 1906: Trace = car(T);
! 1907: C = Trace[0];
! 1908: Ind = Trace[1];
! 1909: Mon = Trace[2];
! 1910: Den = Trace[3];
! 1911: if ( !Coef[Ind] )
! 1912: compute_coef_by_trace(Ind,Tr,Coef);
! 1913:
! 1914: /* XXX */
! 1915: CT = newvect(Nini);
! 1916: for ( J = 0; J < Nini; J++ )
! 1917: CT[J] = (C*CI[J]+Mon*Coef[Ind][J])/Den;
! 1918: CI = CT;
! 1919: }
! 1920: Coef[I] = CI;
! 1921: }
! 1922:
! 1923: extern Gbcheck_DP,Gbcheck_IL$
! 1924:
! 1925: def register_data_for_gbcheck(DPL)
! 1926: {
! 1927: for ( IL = [], I = length(DPL)-1; I >= 0; I-- )
! 1928: IL = cons(I,IL);
! 1929: Gbcheck_DP = newvect(length(DPL),DPL);
! 1930: Gbcheck_IL = IL;
! 1931: }
! 1932:
! 1933: def sp_nf_for_gbcheck(Pair)
! 1934: {
! 1935: SP = dp_sp(Gbcheck_DP[Pair[0]],Gbcheck_DP[Pair[1]]);
! 1936: return dp_nf(Gbcheck_IL,SP,Gbcheck_DP,1);
! 1937: }
! 1938:
! 1939: def gbcheck(B,V,O)
! 1940: {
! 1941: dp_ord(O);
! 1942: D = map(dp_ptod,B,V);
! 1943: L = dp_gr_checklist(D,length(V));
! 1944: DP = L[0]; Plist = L[1];
! 1945: for ( IL = [], I = size(DP)[0]-1; I >= 0; I-- )
! 1946: IL = cons(I,IL);
! 1947: Procs = getopt(proc);
! 1948: if ( type(Procs) == 4 ) {
! 1949: map(ox_reset,Procs);
! 1950: /* register DP in servers */
! 1951: map(ox_cmo_rpc,Procs,"register_data_for_gbcheck",vtol(DP));
! 1952: /* discard return value in stack */
! 1953: map(ox_pop_cmo,Procs);
! 1954: Free = Procs;
! 1955: Busy = [];
! 1956: T = Plist;
! 1957: while ( T != [] || Busy != [] ){
! 1958: if ( Free == [] || T == [] ) {
! 1959: /* someone is working; wait for data */
! 1960: Ready = ox_select(Busy);
! 1961: Busy = setminus(Busy,Ready);
! 1962: Free = append(Ready,Free);
! 1963: for ( ; Ready != []; Ready = cdr(Ready) ) {
! 1964: if ( ox_get(car(Ready)) ) {
! 1965: map(ox_reset,Procs);
! 1966: return 0;
! 1967: }
! 1968: }
! 1969: } else {
! 1970: P = car(Free);
! 1971: Free = cdr(Free);
! 1972: Busy = cons(P,Busy);
! 1973: Pair = car(T);
! 1974: T = cdr(T);
! 1975: ox_cmo_rpc(P,"sp_nf_for_gbcheck",Pair);
! 1976: ox_push_cmd(P,262); /* 262 = OX_popCMO */
! 1977: }
! 1978: }
! 1979: map(ox_reset,Procs);
! 1980: return 1;
! 1981: } else {
! 1982: for ( T = Plist; T != []; T = cdr(T) ) {
! 1983: Pair = T[0];
! 1984: SP = dp_sp(DP[Pair[0]],DP[Pair[1]]);
! 1985: if ( dp_nf(IL,SP,DP,1) )
! 1986: return 0;
! 1987: }
! 1988: return 1;
! 1989: }
! 1990: }
! 1991: end$