![[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 / asir2000 / lib|
Return to gr CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / lib |
1.1 ! noro 1: /* $OpenXM: OpenXM/src/asir99/lib/gr,v 1.1.1.1 1999/11/10 08:12:31 noro Exp $ */
! 2: extern INIT_COUNT,ITOR_FAIL$
! 3: extern REMOTE_MATRIX,REMOTE_NF,REMOTE_VARS$
! 4:
! 5: #define MAX(a,b) ((a)>(b)?(a):(b))
! 6: #define HigherDim 0
! 7: #define ZeroDim 1
! 8: #define MiniPoly 2
! 9:
! 10: /* toplevel functions for Groebner basis computation */
! 11:
! 12: def gr(B,V,O)
! 13: {
! 14: G = dp_gr_main(B,V,0,1,O);
! 15: return G;
! 16: }
! 17:
! 18: def hgr(B,V,O)
! 19: {
! 20: G = dp_gr_main(B,V,1,1,O);
! 21: return G;
! 22: }
! 23:
! 24: def gr_mod(B,V,O,M)
! 25: {
! 26: G = dp_gr_mod_main(B,V,0,M,O);
! 27: return G;
! 28: }
! 29:
! 30: def hgr_mod(B,V,O,M)
! 31: {
! 32: G = dp_gr_mod_main(B,V,1,M,O);
! 33: return G;
! 34: }
! 35:
! 36: /* toplevel functions for change-of-ordering */
! 37:
! 38: def lex_hensel(B,V,O,W,H)
! 39: {
! 40: G = dp_gr_main(B,V,H,1,O);
! 41: return tolex(G,V,O,W);
! 42: }
! 43:
! 44: def lex_hensel_gsl(B,V,O,W,H)
! 45: {
! 46: G = dp_gr_main(B,V,H,1,O);
! 47: return tolex_gsl(G,V,O,W);
! 48: }
! 49:
! 50: def gr_minipoly(B,V,O,P,V0,H)
! 51: {
! 52: G = dp_gr_main(B,V,H,1,O);
! 53: return minipoly(G,V,O,P,V0);
! 54: }
! 55:
! 56: def lex_tl(B,V,O,W,H)
! 57: {
! 58: G = dp_gr_main(B,V,H,1,O);
! 59: return tolex_tl(G,V,O,W,H);
! 60: }
! 61:
! 62: def tolex_tl(G0,V,O,W,H)
! 63: {
! 64: N = length(V); HM = hmlist(G0,V,O); ZD = zero_dim(HM,V,O);
! 65: for ( I = 0; ; I++ ) {
! 66: M = lprime(I);
! 67: if ( !valid_modulus(HM,M) )
! 68: continue;
! 69: if ( ZD ) {
! 70: if ( G3 = dp_gr_main(G0,W,H,-M,3) )
! 71: for ( J = 0; ; J++ )
! 72: if ( G2 = dp_gr_main(G3,W,0,-lprime(J),2) )
! 73: return G2;
! 74: } else if ( G2 = dp_gr_main(G0,W,H,-M,2) )
! 75: return G2;
! 76: }
! 77: }
! 78:
! 79: def tolex(G0,V,O,W)
! 80: {
! 81: TM = TE = TNF = 0;
! 82: N = length(V); HM = hmlist(G0,V,O); ZD = zero_dim(HM,V,O);
! 83: if ( !ZD )
! 84: error("tolex : ideal is not zero-dimensional!");
! 85: MB = dp_mbase(map(dp_ptod,HM,V));
! 86: for ( J = 0; ; J++ ) {
! 87: M = lprime(J);
! 88: if ( !valid_modulus(HM,M) )
! 89: continue;
! 90: T0 = time()[0]; GM = tolexm(G0,V,O,W,M); TM += time()[0] - T0;
! 91: dp_ord(2);
! 92: DL = map(dp_etov,map(dp_ht,map(dp_ptod,GM,W)));
! 93: D = newvect(N); TL = [];
! 94: do
! 95: TL = cons(dp_dtop(dp_vtoe(D),W),TL);
! 96: while ( nextm(D,DL,N) );
! 97: L = npos_check(DL); NPOSV = L[0]; DIM = L[1];
! 98: T0 = time()[0]; NF = gennf(G0,TL,V,O,W[N-1],1)[0];
! 99: TNF += time()[0] - T0;
! 100: T0 = time()[0];
! 101: R = tolex_main(V,O,NF,GM,M,MB);
! 102: TE += time()[0] - T0;
! 103: if ( R ) {
! 104: if ( dp_gr_print() )
! 105: print("mod="+rtostr(TM)+",nf="+rtostr(TNF)+",eq="+rtostr(TE));
! 106: return R;
! 107: }
! 108: }
! 109: }
! 110:
! 111: def tolex_gsl(G0,V,O,W)
! 112: {
! 113: TM = TE = TNF = 0;
! 114: N = length(V); HM = hmlist(G0,V,O); ZD = zero_dim(HM,V,O);
! 115: MB = dp_mbase(map(dp_ptod,HM,V));
! 116: if ( !ZD )
! 117: error("tolex_gsl : ideal is not zero-dimensional!");
! 118: for ( J = 0; ; J++ ) {
! 119: M = lprime(J);
! 120: if ( !valid_modulus(HM,M) )
! 121: continue;
! 122: T0 = time()[0]; GM = tolexm(G0,V,O,W,M); TM += time()[0] - T0;
! 123: dp_ord(2);
! 124: DL = map(dp_etov,map(dp_ht,map(dp_ptod,GM,W)));
! 125: D = newvect(N); TL = [];
! 126: do
! 127: TL = cons(dp_dtop(dp_vtoe(D),W),TL);
! 128: while ( nextm(D,DL,N) );
! 129: L = npos_check(DL); NPOSV = L[0]; DIM = L[1];
! 130: if ( NPOSV >= 0 ) {
! 131: V0 = W[NPOSV];
! 132: T0 = time()[0]; NFL = gennf(G0,TL,V,O,V0,1);
! 133: TNF += time()[0] - T0;
! 134: T0 = time()[0];
! 135: R = tolex_gsl_main(G0,V,O,W,NFL,NPOSV,GM,M,MB);
! 136: TE += time()[0] - T0;
! 137: } else {
! 138: T0 = time()[0]; NF = gennf(G0,TL,V,O,W[N-1],1)[0];
! 139: TNF += time()[0] - T0;
! 140: T0 = time()[0];
! 141: R = tolex_main(V,O,NF,GM,M,MB);
! 142: TE += time()[0] - T0;
! 143: }
! 144: if ( R ) {
! 145: if ( dp_gr_print() )
! 146: print("mod="+rtostr(TM)+",nf="+rtostr(TNF)+",eq="+rtostr(TE));
! 147: return R;
! 148: }
! 149: }
! 150: }
! 151:
! 152: def termstomat(NF,TERMS,MB,MOD)
! 153: {
! 154: DN = NF[1];
! 155: NF = NF[0];
! 156: N = length(MB);
! 157: M = length(TERMS);
! 158: MAT = newmat(N,M);
! 159: W = newvect(N);
! 160: Len = length(NF);
! 161: for ( I = 0; I < M; I++ ) {
! 162: T = TERMS[I];
! 163: for ( K = 0; K < Len; K++ )
! 164: if ( T == NF[K][1] )
! 165: break;
! 166: dptov(NF[K][0],W,MB);
! 167: for ( J = 0; J < N; J++ )
! 168: MAT[J][I] = W[J];
! 169: }
! 170: return [henleq_prep(MAT,MOD),DN];
! 171: }
! 172:
! 173: def tolex_gsl_main(G0,V,O,W,NFL,NPOSV,GM,M,MB)
! 174: {
! 175: NF = NFL[0]; PS = NFL[1]; GI = NFL[2];
! 176: V0 = W[NPOSV]; N = length(W);
! 177: DIM = length(MB);
! 178: DV = newvect(DIM);
! 179: TERMS = gather_terms(GM,W,M,NPOSV);
! 180: Len = length(TERMS);
! 181: dp_ord(O); RHS = termstomat(NF,map(dp_ptod,TERMS,V),MB,M);
! 182: for ( T = GM; T != []; T = cdr(T) )
! 183: if ( vars(car(T)) == [V0] )
! 184: break;
! 185: dp_ord(0); NHT = nf_tab_gsl(dp_ptod(V0^deg(car(T),V0),V),NF);
! 186: dptov(NHT[0],DV,MB);
! 187: B = hen_ttob_gsl([DV,NHT[1]],RHS,TERMS,M);
! 188: if ( !B )
! 189: return 0;
! 190: for ( I = 0, U = B[1]*V0^deg(car(T),V0); I < Len; I++ )
! 191: U += B[0][I]*TERMS[I];
! 192: DN0 = diff(U,V0);
! 193: dp_ord(O); DN0NF = nf_tab_gsl(dp_ptod(DN0,V),NF);
! 194: SL = [[V0,U,DN0]];
! 195: for ( I = N-1, LCM = 1; I >= 0; I-- ) {
! 196: if ( I == NPOSV )
! 197: continue;
! 198: V1 = W[I];
! 199: dp_ord(O); L = nf(GI,DN0NF[0]*dp_ptod(-LCM*V1,V),DN0NF[1],PS);
! 200: L = remove_cont(L);
! 201: dptov(L[0],DV,MB);
! 202: dp_ord(O); B = hen_ttob_gsl([DV,L[1]],RHS,TERMS,M);
! 203: if ( !B )
! 204: return 0;
! 205: for ( K = 0, R = 0; K < Len; K++ )
! 206: R += B[0][K]*TERMS[K];
! 207: LCM *= B[1];
! 208: SL = cons(cons(V1,[R,LCM]),SL);
! 209: print(["DN",B[1]]);
! 210: }
! 211: return SL;
! 212: }
! 213:
! 214: def hen_ttob_gsl(LHS,RHS,TERMS,M)
! 215: {
! 216: LDN = LHS[1]; RDN = RHS[1]; LCM = ilcm(LDN,RDN);
! 217: L1 = idiv(LCM,LDN); R1 = idiv(LCM,RDN);
! 218: T0 = time()[0];
! 219: S = henleq_gsl(RHS[0],LHS[0]*L1,M);
! 220: print(["henleq_gsl",time()[0]-T0]);
! 221: N = length(TERMS);
! 222: return [S[0],S[1]*R1];
! 223: }
! 224:
! 225: def gather_terms(GM,W,M,NPOSV)
! 226: {
! 227: N = length(W); V0 = W[NPOSV];
! 228: for ( T = GM; T != []; T = cdr(T) ) {
! 229: if ( vars(car(T)) == [V0] )
! 230: break;
! 231: }
! 232: U = car(T); DU = diff(U,V0);
! 233: R = tpoly(cdr(p_terms(U,W,2)));
! 234: for ( I = 0; I < N; I++ ) {
! 235: if ( I == NPOSV )
! 236: continue;
! 237: V1 = W[I];
! 238: for ( T = GM; T != []; T = cdr(T) )
! 239: if ( member(V1,vars(car(T))) )
! 240: break;
! 241: P = car(T);
! 242: R += tpoly(p_terms(srem(DU*coef(P,0,V1),U,M),W,2));
! 243: }
! 244: return p_terms(R,W,2);
! 245: }
! 246:
! 247: def tpoly(L)
! 248: {
! 249: for ( R = 0; L != []; L = cdr(L) )
! 250: R += car(L);
! 251: return R;
! 252: }
! 253:
! 254: def dptov(P,W,MB)
! 255: {
! 256: N = size(W)[0];
! 257: for ( I = 0; I < N; I++ )
! 258: W[I] = 0;
! 259: for ( I = 0, S = MB; P; P = dp_rest(P) ) {
! 260: HM = dp_hm(P); C = dp_hc(HM); T = dp_ht(HM);
! 261: for ( ; T != car(S); S = cdr(S), I++ );
! 262: W[I] = C;
! 263: I++; S = cdr(S);
! 264: }
! 265: }
! 266:
! 267: def tolex_main(V,O,NF,GM,M,MB)
! 268: {
! 269: DIM = length(MB);
! 270: DV = newvect(DIM);
! 271: for ( T = GM, SL = [], LCM = 1; T != []; T = cdr(T) ) {
! 272: S = p_terms(car(T),V,2);
! 273: dp_ord(O); RHS = termstomat(NF,map(dp_ptod,cdr(S),V),MB,M);
! 274: dp_ord(0); NHT = nf_tab_gsl(dp_ptod(LCM*car(S),V),NF);
! 275: dptov(NHT[0],DV,MB);
! 276: dp_ord(O); B = hen_ttob_gsl([DV,NHT[1]],RHS,cdr(S),M);
! 277: if ( !B )
! 278: return 0;
! 279: Len = length(S);
! 280: LCM *= B[1];
! 281: for ( U = LCM*car(S), I = 1; I < Len; I++ )
! 282: U += B[0][I-1]*S[I];
! 283: R = ptozp(U);
! 284: SL = cons(R,SL);
! 285: print(["DN",B[1]]);
! 286: }
! 287: return SL;
! 288: }
! 289:
! 290: def reduce_dn(L)
! 291: {
! 292: NM = L[0]; DN = L[1]; V = vars(NM);
! 293: T = remove_cont([dp_ptod(NM,V),DN]);
! 294: return [dp_dtop(T[0],V),T[1]];
! 295: }
! 296:
! 297: /* a function for computation of minimal polynomial */
! 298:
! 299: def minipoly(G0,V,O,P,V0)
! 300: {
! 301: if ( !zero_dim(hmlist(G0,V,O),V,O) )
! 302: error("tolex : ideal is not zero-dimensional!");
! 303:
! 304: G1 = cons(V0-P,G0);
! 305: O1 = [[0,1],[O,length(V)]];
! 306: V1 = cons(V0,V);
! 307: W = append(V,[V0]);
! 308:
! 309: N = length(V1);
! 310: dp_ord(O1);
! 311: HM = hmlist(G1,V1,O1);
! 312: MB = dp_mbase(map(dp_ptod,HM,V1));
! 313: dp_ord(O);
! 314:
! 315: for ( J = 0; ; J++ ) {
! 316: M = lprime(J);
! 317: if ( !valid_modulus(HM,M) )
! 318: continue;
! 319: MP = minipolym(G0,V,O,P,V0,M);
! 320: for ( D = deg(MP,V0), TL = [], J = 0; J <= D; J++ )
! 321: TL = cons(V0^J,TL);
! 322: NF = gennf(G1,TL,V1,O1,V0,1)[0];
! 323: R = tolex_main(V1,O1,NF,[MP],M,MB);
! 324: return R[0];
! 325: }
! 326: }
! 327:
! 328: /* subroutines */
! 329:
! 330: def gennf(G,TL,V,O,V0,FLAG)
! 331: {
! 332: N = length(V); Len = length(G); dp_ord(O); PS = newvect(Len);
! 333: for ( I = 0, T = G, HL = []; T != []; T = cdr(T), I++ ) {
! 334: PS[I] = dp_ptod(car(T),V); HL = cons(dp_ht(PS[I]),HL);
! 335: }
! 336: for ( I = 0, DTL = []; TL != []; TL = cdr(TL) )
! 337: DTL = cons(dp_ptod(car(TL),V),DTL);
! 338: for ( I = Len - 1, GI = []; I >= 0; I-- )
! 339: GI = cons(I,GI);
! 340: T = car(DTL); DTL = cdr(DTL);
! 341: H = [nf(GI,T,T,PS)];
! 342:
! 343: USE_TAB = (FLAG != 0);
! 344: if ( USE_TAB ) {
! 345: T0 = time()[0];
! 346: MB = dp_mbase(HL); DIM = length(MB);
! 347: U = dp_ptod(V0,V);
! 348: UTAB = newvect(DIM);
! 349: for ( I = 0; I < DIM; I++ ) {
! 350: UTAB[I] = [MB[I],remove_cont(dp_true_nf(GI,U*MB[I],PS,1))];
! 351: if ( dp_gr_print() )
! 352: print(".",2);
! 353: }
! 354: print("");
! 355: TTAB = time()[0]-T0;
! 356: }
! 357:
! 358: T0 = time()[0];
! 359: for ( LCM = 1; DTL != []; ) {
! 360: if ( dp_gr_print() )
! 361: print(".",2);
! 362: T = car(DTL); DTL = cdr(DTL);
! 363: if ( L = search_redble(T,H) ) {
! 364: DD = dp_subd(T,L[1]);
! 365: if ( USE_TAB && (DD == U) ) {
! 366: NF = nf_tab(L[0],UTAB);
! 367: NF = [NF[0],dp_hc(L[1])*NF[1]*T];
! 368: } else
! 369: NF = nf(GI,L[0]*dp_subd(T,L[1]),dp_hc(L[1])*T,PS);
! 370: } else
! 371: NF = nf(GI,T,T,PS);
! 372: NF = remove_cont(NF);
! 373: H = cons(NF,H);
! 374: LCM = ilcm(LCM,dp_hc(NF[1]));
! 375: }
! 376: TNF = time()[0]-T0;
! 377: if ( dp_gr_print() )
! 378: print("gennf(TAB="+rtostr(TTAB)+" NF="+rtostr(TNF)+")");
! 379: return [[map(adj_dn,H,LCM),LCM],PS,GI];
! 380: }
! 381:
! 382: def adj_dn(P,D)
! 383: {
! 384: return [(idiv(D,dp_hc(P[1])))*P[0],dp_ht(P[1])];
! 385: }
! 386:
! 387: def hen_ttob(T,NF,LHS,V,MOD)
! 388: {
! 389: if ( length(T) == 1 )
! 390: return car(T);
! 391: T0 = time()[0]; M = etom(leq_nf(T,NF,LHS,V)); TE = time()[0] - T0;
! 392: T0 = time()[0]; U = henleq(M,MOD); TH = time()[0] - T0;
! 393: if ( dp_gr_print() ) {
! 394: print("(etom="+rtostr(TE)+" hen="+rtostr(TH)+")");
! 395: }
! 396: return U ? vtop(T,U,LHS) : 0;
! 397: }
! 398:
! 399: def vtop(S,L,GSL)
! 400: {
! 401: U = L[0]; H = L[1];
! 402: if ( GSL ) {
! 403: for ( A = 0, I = 0; S != []; S = cdr(S), I++ )
! 404: A += U[I]*car(S);
! 405: return [A,H];
! 406: } else {
! 407: for ( A = H*car(S), S = cdr(S), I = 0; S != []; S = cdr(S), I++ )
! 408: A += U[I]*car(S);
! 409: return ptozp(A);
! 410: }
! 411: }
! 412:
! 413: def leq_nf(TL,NF,LHS,V)
! 414: {
! 415: TLen = length(NF);
! 416: T = newvect(TLen); M = newvect(TLen);
! 417: for ( I = 0; I < TLen; I++ ) {
! 418: T[I] = dp_ht(NF[I][1]);
! 419: M[I] = dp_hc(NF[I][1]);
! 420: }
! 421: Len = length(TL); INDEX = newvect(Len); COEF = newvect(Len);
! 422: for ( L = TL, J = 0; L != []; L = cdr(L), J++ ) {
! 423: D = dp_ptod(car(L),V);
! 424: for ( I = 0; I < TLen; I++ )
! 425: if ( D == T[I] )
! 426: break;
! 427: INDEX[J] = I; COEF[J] = strtov("u"+rtostr(J));
! 428: }
! 429: if ( !LHS ) {
! 430: COEF[0] = 1; NM = 0; DN = 1;
! 431: } else {
! 432: NM = LHS[0]; DN = LHS[1];
! 433: }
! 434: for ( J = 0, S = -NM; J < Len; J++ ) {
! 435: DNJ = M[INDEX[J]];
! 436: GCD = igcd(DN,DNJ); CS = DNJ/GCD; CJ = DN/GCD;
! 437: S = CS*S + CJ*NF[INDEX[J]][0]*COEF[J];
! 438: DN *= CS;
! 439: }
! 440: for ( D = S, E = []; D; D = dp_rest(D) )
! 441: E = cons(dp_hc(D),E);
! 442: BOUND = LHS ? 0 : 1;
! 443: for ( I = Len - 1, W = []; I >= BOUND; I-- )
! 444: W = cons(COEF[I],W);
! 445: return [E,W];
! 446: }
! 447:
! 448: def nf_tab(F,TAB)
! 449: {
! 450: for ( NM = 0, DN = 1, I = 0; F; F = dp_rest(F) ) {
! 451: T = dp_ht(F);
! 452: for ( ; TAB[I][0] != T; I++);
! 453: NF = TAB[I][1]; N = NF[0]; D = NF[1];
! 454: G = igcd(DN,D); DN1 = idiv(DN,G); D1 = idiv(D,G);
! 455: NM = D1*NM + DN1*dp_hc(F)*N; DN *= D1;
! 456: }
! 457: return [NM,DN];
! 458: }
! 459:
! 460: def nf_tab_gsl(A,NF)
! 461: {
! 462: DN = NF[1];
! 463: NF = NF[0];
! 464: TLen = length(NF);
! 465: for ( R = 0; A; A = dp_rest(A) ) {
! 466: HM = dp_hm(A); C = dp_hc(HM); T = dp_ht(HM);
! 467: for ( I = 0; I < TLen; I++ )
! 468: if ( NF[I][1] == T )
! 469: break;
! 470: R += C*NF[I][0];
! 471: }
! 472: return remove_cont([R,DN]);
! 473: }
! 474:
! 475: def redble(D1,D2,N)
! 476: {
! 477: for ( I = 0; I < N; I++ )
! 478: if ( D1[I] > D2[I] )
! 479: break;
! 480: return I == N ? 1 : 0;
! 481: }
! 482:
! 483: def tolexm(G,V,O,W,M)
! 484: {
! 485: N = length(V); Len = length(G);
! 486: dp_ord(O); setmod(M); PS = newvect(Len);
! 487: for ( I = 0, T = G; T != []; T = cdr(T), I++ )
! 488: PS[I] = dp_mod(dp_ptod(car(T),V),M,[]);
! 489: for ( I = Len-1, HL = []; I >= 0; I-- )
! 490: HL = cons(dp_ht(PS[I]),HL);
! 491: G2 = tolexm_main(PS,HL,V,W,M,ZeroDim);
! 492: L = map(dp_dtop,G2,V);
! 493: return L;
! 494: }
! 495:
! 496: def tolexm_main(PS,HL,V,W,M,FLAG)
! 497: {
! 498: N = length(W); D = newvect(N); Len = size(PS)[0];
! 499: for ( I = Len - 1, GI = []; I >= 0; I-- )
! 500: GI = cons(I,GI);
! 501: MB = dp_mbase(HL); DIM = length(MB);
! 502: U = dp_mod(dp_ptod(W[N-1],V),M,[]);
! 503: UTAB = newvect(DIM);
! 504: for ( I = 0; I < DIM; I++ ) {
! 505: if ( dp_gr_print() )
! 506: print(".",2);
! 507: UTAB[I] = [MB[I],dp_nf_mod(GI,U*dp_mod(MB[I],M,[]),PS,1,M)];
! 508: }
! 509: print("");
! 510: T = dp_mod(dp_ptod(dp_dtop(dp_vtoe(D),W),V),M,[]);
! 511: H = G = [[T,T]];
! 512: DL = []; G2 = [];
! 513: TNF = 0;
! 514: while ( 1 ) {
! 515: if ( dp_gr_print() )
! 516: print(".",2);
! 517: S = nextm(D,DL,N);
! 518: if ( !S )
! 519: break;
! 520: T = dp_mod(dp_ptod(dp_dtop(dp_vtoe(D),W),V),M,[]);
! 521: T0 = time()[0];
! 522: if ( L = search_redble(T,H) ) {
! 523: DD = dp_mod(dp_subd(T,L[1]),M,[]);
! 524: if ( DD == U )
! 525: NT = dp_nf_tab_mod(L[0],UTAB,M);
! 526: else
! 527: NT = dp_nf_mod(GI,L[0]*DD,PS,1,M);
! 528: } else
! 529: NT = dp_nf_mod(GI,T,PS,1,M);
! 530: TNF += time()[0] - T0;
! 531: H = cons([NT,T],H);
! 532: T0 = time()[0];
! 533: L = dp_lnf_mod([NT,T],G,M); N1 = L[0]; N2 = L[1];
! 534: TLNF += time()[0] - T0;
! 535: if ( !N1 ) {
! 536: G2 = cons(N2,G2);
! 537: if ( FLAG == MiniPoly )
! 538: break;
! 539: D1 = newvect(N);
! 540: for ( I = 0; I < N; I++ )
! 541: D1[I] = D[I];
! 542: DL = cons(D1,DL);
! 543: } else
! 544: G = insert(G,L);
! 545: }
! 546: if ( dp_gr_print() )
! 547: print("tolexm(nfm="+rtostr(TNF)+" lnfm="+rtostr(TLNF)+")");
! 548: return G2;
! 549: }
! 550:
! 551: def minipolym(G,V,O,P,V0,M)
! 552: {
! 553: N = length(V); Len = length(G);
! 554: dp_ord(O); setmod(M); PS = newvect(Len);
! 555: for ( I = 0, T = G; T != []; T = cdr(T), I++ )
! 556: PS[I] = dp_mod(dp_ptod(car(T),V),M,[]);
! 557: for ( I = Len-1, HL = []; I >= 0; I-- )
! 558: HL = cons(dp_ht(PS[I]),HL);
! 559: for ( I = Len - 1, GI = []; I >= 0; I-- )
! 560: GI = cons(I,GI);
! 561: MB = dp_mbase(HL); DIM = length(MB); UT = newvect(DIM);
! 562: U = dp_mod(dp_ptod(P,V),M,[]);
! 563: for ( I = 0; I < DIM; I++ )
! 564: UT[I] = [MB[I],dp_nf_mod(GI,U*dp_mod(MB[I],M,[]),PS,1,M)];
! 565: T = dp_mod(<<0>>,M,[]); TT = dp_mod(dp_ptod(1,V),M,[]);
! 566: G = H = [[TT,T]]; TNF = TLNF = 0;
! 567: for ( I = 1; ; I++ ) {
! 568: T = dp_mod(<<I>>,M,[]);
! 569: T0 = time()[0]; NT = dp_nf_tab_mod(H[0][0],UT,M); TNF += time()[0] - T0;
! 570: H = cons([NT,T],H);
! 571: T0 = time()[0]; L = dp_lnf_mod([NT,T],G,M); TLNF += time()[0] - T0;
! 572: if ( !L[0] ) {
! 573: if ( dp_gr_print() ) print(["nfm",TNF,"lnfm",TLNF]);
! 574: return dp_dtop(L[1],[V0]);
! 575: } else
! 576: G = insert(G,L);
! 577: }
! 578: }
! 579:
! 580: def nextm(D,DL,N)
! 581: {
! 582: for ( I = N-1; I >= 0; ) {
! 583: D[I]++;
! 584: for ( T = DL; T != []; T = cdr(T) )
! 585: if ( car(T) == D )
! 586: return 1;
! 587: else if ( redble(car(T),D,N) )
! 588: break;
! 589: if ( T != [] ) {
! 590: for ( J = N-1; J >= I; J-- )
! 591: D[J] = 0;
! 592: I--;
! 593: } else
! 594: break;
! 595: }
! 596: if ( I < 0 )
! 597: return 0;
! 598: else
! 599: return 1;
! 600: }
! 601:
! 602: def search_redble(T,G)
! 603: {
! 604: for ( ; G != []; G = cdr(G) )
! 605: if ( dp_redble(T,car(G)[1]) )
! 606: return car(G);
! 607: return 0;
! 608: }
! 609:
! 610: def insert(G,A)
! 611: {
! 612: if ( G == [] )
! 613: return [A];
! 614: else if ( dp_ht(car(A)) > dp_ht(car(car(G))) )
! 615: return cons(A,G);
! 616: else
! 617: return cons(car(G),insert(cdr(G),A));
! 618: }
! 619:
! 620: #if 0
! 621: def etom(L) {
! 622: E = L[0]; W = L[1];
! 623: LE = length(E); LW = length(W);
! 624: M = newmat(LE,LW+1);
! 625: for(J=0;J<LE;J++) {
! 626: for ( T = E[J]; T && (type(T) == 2); )
! 627: for ( V = var(T), I = 0; I < LW; I++ )
! 628: if ( V == W[I] ) {
! 629: M[J][I] = coef(T,1,V);
! 630: T = coef(T,0,V);
! 631: }
! 632: M[J][LW] = T;
! 633: }
! 634: return M;
! 635: }
! 636: #endif
! 637:
! 638: def etom(L) {
! 639: E = L[0]; W = L[1];
! 640: LE = length(E); LW = length(W);
! 641: M = newmat(LE,LW+1);
! 642: for(J=0;J<LE;J++) {
! 643: for ( I = 0, T = E[J]; I < LW; I++ ) {
! 644: M[J][I] = coef(T,1,W[I]); T = coef(T,0,W[I]);
! 645: }
! 646: M[J][LW] = T;
! 647: }
! 648: return M;
! 649: }
! 650:
! 651: def calcb_old(M) {
! 652: N = 2*M;
! 653: T = gr_sqrt(N);
! 654: if ( T^2 <= N && N < (T+1)^2 )
! 655: return idiv(T,2);
! 656: else
! 657: error("afo");
! 658: }
! 659:
! 660: def calcb_special(PK,P,K) { /* PK = P^K */
! 661: N = 2*PK;
! 662: T = sqrt_special(N,2,P,K);
! 663: if ( T^2 <= N && N < (T+1)^2 )
! 664: return idiv(T,2);
! 665: else
! 666: error("afo");
! 667: }
! 668:
! 669: def sqrt_special(A,C,M,K) { /* A = C*M^K */
! 670: L = idiv(K,2); B = M^L;
! 671: if ( K % 2 )
! 672: C *= M;
! 673: D = 2^K; X = idiv((gr_sqrt(C*D^2)+1)*B,D)+1;
! 674: while ( 1 )
! 675: if ( (Y = X^2) <= A )
! 676: return X;
! 677: else
! 678: X = idiv(A + Y,2*X);
! 679: }
! 680:
! 681: def gr_sqrt(A) {
! 682: for ( J = 0, T = A; T >= 2^27; J++ ) {
! 683: T = idiv(T,2^27)+1;
! 684: }
! 685: for ( I = 0; T >= 2; I++ ) {
! 686: S = idiv(T,2);
! 687: if ( T = S+S )
! 688: T = S;
! 689: else
! 690: T = S+1;
! 691: }
! 692: X = (2^27)^idiv(J,2)*2^idiv(I,2);
! 693: while ( 1 ) {
! 694: if ( (Y=X^2) < A )
! 695: X += X;
! 696: else if ( Y == A )
! 697: return X;
! 698: else
! 699: break;
! 700: }
! 701: while ( 1 )
! 702: if ( (Y = X^2) <= A )
! 703: return X;
! 704: else
! 705: X = idiv(A + Y,2*X);
! 706: }
! 707:
! 708: #define ABS(a) ((a)>=0?(a):(-a))
! 709:
! 710: def inttorat_asir(C,M,B)
! 711: {
! 712: if ( M < 0 )
! 713: M = -M;
! 714: C %= M;
! 715: if ( C < 0 )
! 716: C += M;
! 717: U1 = 0; U2 = M; V1 = 1; V2 = C;
! 718: while ( V2 >= B ) {
! 719: L = iqr(U2,V2); Q = L[0]; R2 = L[1];
! 720: R1 = U1 - Q*V1;
! 721: U1 = V1; U2 = V2;
! 722: V1 = R1; V2 = R2;
! 723: }
! 724: if ( ABS(V1) >= B )
! 725: return 0;
! 726: else
! 727: if ( V1 < 0 )
! 728: return [-V2,-V1];
! 729: else
! 730: return [V2,V1];
! 731: }
! 732:
! 733: def intvtoratv(V,M,B) {
! 734: if ( !B )
! 735: B = 1;
! 736: N = size(V)[0];
! 737: W = newvect(N);
! 738: if ( ITOR_FAIL >= 0 ) {
! 739: if ( V[ITOR_FAIL] ) {
! 740: T = inttorat(V[ITOR_FAIL],M,B);
! 741: if ( !T ) {
! 742: if ( dp_gr_print() ) {
! 743: print("F",2);
! 744: }
! 745: return 0;
! 746: }
! 747: }
! 748: }
! 749: for ( I = 0, DN = 1; I < N; I++ )
! 750: if ( V[I] ) {
! 751: T = inttorat((V[I]*DN) % M,M,B);
! 752: if ( !T ) {
! 753: ITOR_FAIL = I;
! 754: if ( dp_gr_print() ) {
! 755: #if 0
! 756: print("intvtoratv : failed at I = ",0); print(ITOR_FAIL);
! 757: #endif
! 758: print("F",2);
! 759: }
! 760: return 0;
! 761: } else {
! 762: for( J = 0; J < I; J++ )
! 763: W[J] *= T[1];
! 764: W[I] = T[0]; DN *= T[1];
! 765: }
! 766: }
! 767: return [W,DN];
! 768: }
! 769:
! 770: def nf(B,G,M,PS)
! 771: {
! 772: for ( D = 0; G; ) {
! 773: for ( U = 0, L = B; L != []; L = cdr(L) ) {
! 774: if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
! 775: GCD = igcd(dp_hc(G),dp_hc(R));
! 776: CG = idiv(dp_hc(R),GCD); CR = idiv(dp_hc(G),GCD);
! 777: U = CG*G-dp_subd(G,R)*CR*R;
! 778: if ( !U )
! 779: return [D,M];
! 780: D *= CG; M *= CG;
! 781: break;
! 782: }
! 783: }
! 784: if ( U )
! 785: G = U;
! 786: else {
! 787: D += dp_hm(G); G = dp_rest(G);
! 788: }
! 789: }
! 790: return [D,M];
! 791: }
! 792:
! 793: def remove_cont(L)
! 794: {
! 795: if ( type(L[1]) == 1 ) {
! 796: T = remove_cont([L[0],L[1]*<<0>>]);
! 797: return [T[0],dp_hc(T[1])];
! 798: } else if ( !L[0] )
! 799: return [0,dp_ptozp(L[1])];
! 800: else if ( !L[1] )
! 801: return [dp_ptozp(L[0]),0];
! 802: else {
! 803: A0 = dp_ptozp(L[0]); A1 = dp_ptozp(L[1]);
! 804: C0 = idiv(dp_hc(L[0]),dp_hc(A0)); C1 = idiv(dp_hc(L[1]),dp_hc(A1));
! 805: GCD = igcd(C0,C1); M0 = idiv(C0,GCD); M1 = idiv(C1,GCD);
! 806: return [M0*A0,M1*A1];
! 807: }
! 808: }
! 809:
! 810: def union(A,B)
! 811: {
! 812: for ( T = B; T != []; T = cdr(T) )
! 813: A = union1(A,car(T));
! 814: return A;
! 815: }
! 816:
! 817: def union1(A,E)
! 818: {
! 819: if ( A == [] )
! 820: return [E];
! 821: else if ( car(A) == E )
! 822: return A;
! 823: else
! 824: return cons(car(A),union1(cdr(A),E));
! 825: }
! 826:
! 827: def setminus(A,B) {
! 828: for ( T = reverse(A), R = []; T != []; T = cdr(T) ) {
! 829: for ( S = B, M = car(T); S != []; S = cdr(S) )
! 830: if ( car(S) == M )
! 831: break;
! 832: if ( S == [] )
! 833: R = cons(M,R);
! 834: }
! 835: return R;
! 836: }
! 837:
! 838: def member(A,L) {
! 839: for ( ; L != []; L = cdr(L) )
! 840: if ( A == car(L) )
! 841: return 1;
! 842: return 0;
! 843: }
! 844:
! 845: /* several functions for computation of normal forms etc. */
! 846:
! 847: def p_nf(P,B,V,O) {
! 848: dp_ord(O); DP = dp_ptod(P,V);
! 849: N = length(B); DB = newvect(N);
! 850: for ( I = N-1, IL = []; I >= 0; I-- ) {
! 851: DB[I] = dp_ptod(B[I],V);
! 852: IL = cons(I,IL);
! 853: }
! 854: return dp_dtop(dp_nf(IL,DP,DB,1),V);
! 855: }
! 856:
! 857: def p_true_nf(P,B,V,O) {
! 858: dp_ord(O); DP = dp_ptod(P,V);
! 859: N = length(B); DB = newvect(N);
! 860: for ( I = N-1, IL = []; I >= 0; I-- ) {
! 861: DB[I] = dp_ptod(B[I],V);
! 862: IL = cons(I,IL);
! 863: }
! 864: L = dp_true_nf(IL,DP,DB,1);
! 865: return [dp_dtop(L[0],V),L[1]];
! 866: }
! 867:
! 868: def p_terms(D,V,O)
! 869: {
! 870: dp_ord(O);
! 871: for ( L = [], T = dp_ptod(D,V); T; T = dp_rest(T) )
! 872: L = cons(dp_dtop(dp_ht(T),V),L);
! 873: return reverse(L);
! 874: }
! 875:
! 876: def dp_terms(D,V)
! 877: {
! 878: for ( L = [], T = D; T; T = dp_rest(T) )
! 879: L = cons(dp_dtop(dp_ht(T),V),L);
! 880: return reverse(L);
! 881: }
! 882:
! 883: def gb_comp(A,B)
! 884: {
! 885: for ( T = A; T != []; T = cdr(T) ) {
! 886: for ( S = B, M = car(T), N = -M; S != []; S = cdr(S) )
! 887: if ( car(S) == M || car(S) == N )
! 888: break;
! 889: if ( S == [] )
! 890: break;
! 891: }
! 892: return T == [] ? 1 : 0;
! 893: }
! 894:
! 895: def zero_dim(G,V,O) {
! 896: dp_ord(O);
! 897: HL = map(dp_dtop,map(dp_ht,map(dp_ptod,G,V)),V);
! 898: for ( L = []; HL != []; HL = cdr(HL) )
! 899: if ( length(vars(car(HL))) == 1 )
! 900: L = cons(car(HL),L);
! 901: return length(vars(L)) == length(V) ? 1 : 0;
! 902: }
! 903:
! 904: def hmlist(G,V,O) {
! 905: dp_ord(O);
! 906: return map(dp_dtop,map(dp_hm,map(dp_ptod,G,V)),V);
! 907: }
! 908:
! 909: def valid_modulus(HL,M) {
! 910: V = vars(HL);
! 911: for ( T = HL; T != []; T = cdr(T) )
! 912: if ( !dp_mod(dp_ptod(car(T),V),M,[]) )
! 913: break;
! 914: return T == [] ? 1 : 0;
! 915: }
! 916:
! 917: def npos_check(DL) {
! 918: N = size(car(DL))[0];
! 919: if ( length(DL) != N )
! 920: return [-1,0];
! 921: D = newvect(N);
! 922: for ( I = 0; I < N; I++ ) {
! 923: for ( J = 0; J < N; J++ )
! 924: D[J] = 0;
! 925: D[I] = 1;
! 926: for ( T = DL; T != []; T = cdr(T) )
! 927: if ( D == car(T) )
! 928: break;
! 929: if ( T != [] )
! 930: DL = setminus(DL,[car(T)]);
! 931: }
! 932: if ( length(DL) != 1 )
! 933: return [-1,0];
! 934: U = car(DL);
! 935: for ( I = 0, J = 0, I0 = -1; I < N; I++ )
! 936: if ( U[I] ) {
! 937: I0 = I; J++;
! 938: }
! 939: if ( J != 1 )
! 940: return [-1,0];
! 941: else
! 942: return [I0,U[I0]];
! 943: }
! 944:
! 945: def mult_mat(L,TAB,MB)
! 946: {
! 947: A = L[0]; DN0 = L[1];
! 948: for ( NM = 0, DN = 1, I = 0; A; A = dp_rest(A) ) {
! 949: H = dp_ht(A);
! 950: for ( ; MB[I] != H; I++ );
! 951: NM1 = TAB[I][0]; DN1 = TAB[I][1]; I++;
! 952: GCD = igcd(DN,DN1); C = DN1/GCD; C1 = DN/GCD;
! 953: NM = C*NM + C1*dp_hc(A)*NM1;
! 954: DN *= C;
! 955: }
! 956: Z=remove_cont([NM,DN*DN0]);
! 957: return Z;
! 958: }
! 959:
! 960: def sepm(MAT)
! 961: {
! 962: S = size(MAT); N = S[0]; M = S[1]-1;
! 963: A = newmat(N,M); B = newvect(N);
! 964: for ( I = 0; I < N; I++ )
! 965: for ( J = 0, T1 = MAT[I], T2 = A[I]; J < M; J++ )
! 966: T2[J] = T1[J];
! 967: for ( I = 0; I < N; I++ )
! 968: B[I] = MAT[I][M];
! 969: return [A,B];
! 970: }
! 971:
! 972: def henleq(M,MOD)
! 973: {
! 974: SIZE = size(M); ROW = SIZE[0]; COL = SIZE[1];
! 975: W = newvect(COL);
! 976: L = sepm(M); A = L[0]; B = L[1];
! 977: COUNT = INIT_COUNT?INIT_COUNT:idiv(max_mag(M),54);
! 978: if ( !COUNT )
! 979: COUNT = 1;
! 980:
! 981: TINV = TC = TR = TS = TM = TDIV = 0;
! 982:
! 983: T0 = time()[0];
! 984: L = geninvm_swap(A,MOD); INV = L[0]; INDEX = L[1];
! 985: TS += time()[0] - T0;
! 986:
! 987: COL1 = COL - 1;
! 988: AA = newmat(COL1,COL1); BB = newvect(COL1);
! 989: for ( I = 0; I < COL1; I++ ) {
! 990: for ( J = 0, T = AA[I], S = A[INDEX[I]]; J < COL1; J++ )
! 991: T[J] = S[J];
! 992: BB[I] = B[INDEX[I]];
! 993: }
! 994: if ( COL1 != ROW ) {
! 995: RESTA = newmat(ROW-COL1,COL1); RESTB = newvect(ROW-COL1);
! 996: for ( ; I < ROW; I++ ) {
! 997: for ( J = 0, T = RESTA[I-COL1], S = A[INDEX[I]]; J < COL1; J++ )
! 998: T[J] = S[J];
! 999: RESTB[I-COL1] = B[INDEX[I]];
! 1000: }
! 1001: } else
! 1002: RESTA = RESTB = 0;
! 1003:
! 1004: MOD2 = idiv(MOD,2);
! 1005: for ( I = 0, C = BB, X = 0, PK = 1, CCC = 0, ITOR_FAIL = -1; ;
! 1006: I++, PK *= MOD ) {
! 1007: if ( COUNT == CCC ) {
! 1008: CCC = 0;
! 1009: T0 = time()[0];
! 1010: ND = intvtoratv(X,PK,ishift(calcb_special(PK,MOD,I),32));
! 1011: TR += time()[0]-T0;
! 1012: if ( ND ) {
! 1013: T0 = time()[0];
! 1014: F = ND[0]; LCM = ND[1]; T = AA*F+LCM*BB;
! 1015: TM += time()[0]-T0;
! 1016: if ( zerovector(T) ) {
! 1017: T0 = time()[0]; T = RESTA*F+LCM*RESTB; TM += time()[0]-T0;
! 1018: if ( zerovector(T) ) {
! 1019: #if 0
! 1020: if ( dp_gr_print() ) print(["init",TS,"pinv",TINV,"c",TC,"div",TDIV,"rat",TR,"mul",TM]);
! 1021: #endif
! 1022: if ( dp_gr_print() ) print("end",2);
! 1023: return [F,LCM];
! 1024: } else
! 1025: return 0;
! 1026: }
! 1027: } else {
! 1028: #if 0
! 1029: if ( dp_gr_print() ) print(I);
! 1030: #endif
! 1031: }
! 1032: } else {
! 1033: #if 0
! 1034: if ( dp_gr_print() ) print([I,TINV,TC,TDIV]);
! 1035: #endif
! 1036: if ( dp_gr_print() ) print(".",2);
! 1037: CCC++;
! 1038: }
! 1039: T0 = time()[0];
! 1040: XT = sremainder(INV*sremainder(-C,MOD),MOD);
! 1041: XT = map(adj_sgn,XT,MOD,MOD2);
! 1042: TINV += time()[0] - T0;
! 1043: X += XT*PK;
! 1044: T0 = time()[0];
! 1045: C += mul_mat_vect_int(AA,XT);
! 1046: TC += time()[0] - T0;
! 1047: T0 = time()[0]; C = map(idiv,C,MOD); TDIV += time()[0] - T0;
! 1048: }
! 1049: }
! 1050:
! 1051: def henleq_prep(A,MOD)
! 1052: {
! 1053: SIZE = size(A); ROW = SIZE[0]; COL = SIZE[1];
! 1054: L = geninvm_swap(A,MOD); INV = L[0]; INDEX = L[1];
! 1055: AA = newmat(COL,COL);
! 1056: for ( I = 0; I < COL; I++ )
! 1057: for ( J = 0, T = AA[I], S = A[INDEX[I]]; J < COL; J++ )
! 1058: T[J] = S[J];
! 1059: if ( COL != ROW ) {
! 1060: RESTA = newmat(ROW-COL,COL);
! 1061: for ( ; I < ROW; I++ )
! 1062: for ( J = 0, T = RESTA[I-COL], S = A[INDEX[I]]; J < COL; J++ )
! 1063: T[J] = S[J];
! 1064: } else
! 1065: RESTA = 0;
! 1066: return [[A,AA,RESTA],L];
! 1067: }
! 1068:
! 1069: def henleq_gsl(L,B,MOD)
! 1070: {
! 1071: AL = L[0]; INVL = L[1];
! 1072: A = AL[0]; AA = AL[1]; RESTA = AL[2];
! 1073: INV = INVL[0]; INDEX = INVL[1];
! 1074: SIZE = size(A); ROW = SIZE[0]; COL = SIZE[1];
! 1075: BB = newvect(COL);
! 1076: for ( I = 0; I < COL; I++ )
! 1077: BB[I] = B[INDEX[I]];
! 1078: if ( COL != ROW ) {
! 1079: RESTB = newvect(ROW-COL);
! 1080: for ( ; I < ROW; I++ )
! 1081: RESTB[I-COL] = B[INDEX[I]];
! 1082: } else
! 1083: RESTB = 0;
! 1084:
! 1085: COUNT = INIT_COUNT?INIT_COUNT:idiv(MAX(max_mag(A),max_mag_vect(B)),54);
! 1086: if ( !COUNT )
! 1087: COUNT = 1;
! 1088: MOD2 = idiv(MOD,2);
! 1089: for ( I = 0, C = BB, X = 0, PK = 1, CCC = 0, ITOR_FAIL = -1; ;
! 1090: I++, PK *= MOD ) {
! 1091: if ( zerovector(C) )
! 1092: if ( zerovector(RESTA*X+RESTB) ) {
! 1093: if ( dp_gr_print() ) print("end",0);
! 1094: return [X,1];
! 1095: } else
! 1096: return 0;
! 1097: else if ( COUNT == CCC ) {
! 1098: CCC = 0;
! 1099: ND = intvtoratv(X,PK,ishift(calcb_special(PK,MOD,I),32));
! 1100: if ( ND ) {
! 1101: F = ND[0]; LCM = ND[1]; T = AA*F+LCM*BB;
! 1102: if ( zerovector(T) ) {
! 1103: T = RESTA*F+LCM*RESTB;
! 1104: if ( zerovector(T) ) {
! 1105: if ( dp_gr_print() ) print("end",0);
! 1106: return [F,LCM];
! 1107: } else
! 1108: return 0;
! 1109: }
! 1110: } else {
! 1111: }
! 1112: } else {
! 1113: if ( dp_gr_print() ) print(".",2);
! 1114: CCC++;
! 1115: }
! 1116: XT = sremainder(INV*sremainder(-C,MOD),MOD);
! 1117: XT = map(adj_sgn,XT,MOD,MOD2);
! 1118: X += XT*PK;
! 1119: C += mul_mat_vect_int(AA,XT);
! 1120: C = map(idiv,C,MOD);
! 1121: }
! 1122: }
! 1123:
! 1124: def adj_sgn(A,M,M2)
! 1125: {
! 1126: return A > M2 ? A-M : A;
! 1127: }
! 1128:
! 1129: def zerovector(C)
! 1130: {
! 1131: if ( !C )
! 1132: return 1;
! 1133: for ( I = size(C)[0]-1; I >= 0 && !C[I]; I-- );
! 1134: if ( I < 0 )
! 1135: return 1;
! 1136: else
! 1137: return 0;
! 1138: }
! 1139:
! 1140: def solvem(INV,COMP,V,MOD)
! 1141: {
! 1142: T = COMP*V;
! 1143: N = size(T)[0];
! 1144: for ( I = 0; I < N; I++ )
! 1145: if ( T[I] % MOD )
! 1146: return 0;
! 1147: return modvect(INV*V,MOD);
! 1148: }
! 1149:
! 1150: def modmat(A,MOD)
! 1151: {
! 1152: if ( !A )
! 1153: return 0;
! 1154: S = size(A); N = S[0]; M = S[1];
! 1155: MAT = newmat(N,M);
! 1156: for ( I = 0, NZ = 0; I < N; I++ )
! 1157: for ( J = 0, T1 = A[I], T2 = MAT[I]; J < M; J++ ) {
! 1158: T2[J] = T1[J] % MOD;
! 1159: NZ = NZ || T2[J];
! 1160: }
! 1161: return NZ?MAT:0;
! 1162: }
! 1163:
! 1164: def modvect(A,MOD)
! 1165: {
! 1166: if ( !A )
! 1167: return 0;
! 1168: N = size(A)[0];
! 1169: VECT = newvect(N);
! 1170: for ( I = 0, NZ = 0; I < N; I++ ) {
! 1171: VECT[I] = A[I] % MOD;
! 1172: NZ = NZ || VECT[I];
! 1173: }
! 1174: return NZ?VECT:0;
! 1175: }
! 1176:
! 1177: def qrmat(A,MOD)
! 1178: {
! 1179: if ( !A )
! 1180: return [0,0];
! 1181: S = size(A); N = S[0]; M = S[1];
! 1182: Q = newmat(N,M); R = newmat(N,M);
! 1183: for ( I = 0, NZQ = 0, NZR = 0; I < N; I++ )
! 1184: for ( J = 0, TA = A[I], TQ = Q[I], TR = R[I]; J < M; J++ ) {
! 1185: L = iqr(TA[J],MOD); TQ[J] = L[0]; TR[J] = L[1];
! 1186: NZQ = NZQ || TQ[J]; NZR = NZR || TR[J];
! 1187: }
! 1188: return [NZQ?Q:0,NZR?R:0];
! 1189: }
! 1190:
! 1191: def qrvect(A,MOD)
! 1192: {
! 1193: if ( !A )
! 1194: return [0,0];
! 1195: N = size(A)[0];
! 1196: Q = newvect(N); R = newvect(N);
! 1197: for ( I = 0, NZQ = 0, NZR = 0; I < N; I++ ) {
! 1198: L = iqr(A[I],MOD); Q[I] = L[0]; R[I] = L[1];
! 1199: NZQ = NZQ || Q[I]; NZR = NZR || R[I];
! 1200: }
! 1201: return [NZQ?Q:0,NZR?R:0];
! 1202: }
! 1203:
! 1204: def max_mag(M)
! 1205: {
! 1206: R = size(M)[0];
! 1207: U = 1;
! 1208: for ( I = 0; I < R; I++ ) {
! 1209: A = max_mag_vect(M[I]);
! 1210: U = MAX(A,U);
! 1211: }
! 1212: return U;
! 1213: }
! 1214:
! 1215: def max_mag_vect(V)
! 1216: {
! 1217: R = size(V)[0];
! 1218: U = 1;
! 1219: for ( I = 0; I < R; I++ ) {
! 1220: A = dp_mag(V[I]*<<0>>);
! 1221: U = MAX(A,U);
! 1222: }
! 1223: return U;
! 1224: }
! 1225:
! 1226: def gsl_check(B,V,S)
! 1227: {
! 1228: N = length(V);
! 1229: U = S[N-1]; M = U[1]; D = U[2];
! 1230: W = setminus(V,[var(M)]);
! 1231: H = uc(); VH = append(W,[H]);
! 1232: for ( T = B; T != []; T = cdr(T) ) {
! 1233: A = car(T);
! 1234: AH = dp_dtop(dp_homo(dp_ptod(A,W)),VH);
! 1235: for ( I = 0, Z = S; I < N-1; I++, Z = cdr(Z) ) {
! 1236: L = car(Z); AH = ptozp(subst(AH,L[0],L[1]/L[2]));
! 1237: }
! 1238: AH = ptozp(subst(AH,H,D));
! 1239: R = srem(AH,M);
! 1240: if ( dp_gr_print() )
! 1241: if ( !R )
! 1242: print([A,"ok"]);
! 1243: else
! 1244: print([A,"bad"]);
! 1245: if ( R )
! 1246: break;
! 1247: }
! 1248: return T == [] ? 1 : 0;
! 1249: }
! 1250:
! 1251: def vs_dim(G,V,O)
! 1252: {
! 1253: HM = hmlist(G,V,O); ZD = zero_dim(HM,V,O);
! 1254: if ( ZD ) {
! 1255: MB = dp_mbase(map(dp_ptod,HM,V));
! 1256: return length(MB);
! 1257: } else
! 1258: error("vs_dim : ideal is not zero-dimensional!");
! 1259: }
! 1260:
! 1261: def dgr(G,V,O,P)
! 1262: {
! 1263: P0 = P[0]; P1 = P[1]; P = [P0,P1];
! 1264: flush(P0); flush(P1);
! 1265: rpc(P0,"dp_gr_main",G,V,0,1,O);
! 1266: rpc(P1,"dp_gr_main",G,V,1,1,O);
! 1267: F = select(P);
! 1268: R = rpcrecv(F[0]); flush(P0); flush(P1);
! 1269: return R;
! 1270: }
! 1271:
! 1272: /* functions for rpc */
! 1273:
! 1274: def register_matrix(M)
! 1275: {
! 1276: REMOTE_MATRIX = M; return 0;
! 1277: }
! 1278:
! 1279: def register_nfv(L)
! 1280: {
! 1281: REMOTE_NF = L[0]; REMOTE_VARS = L[1]; return 0;
! 1282: }
! 1283:
! 1284: def r_ttob(T,M)
! 1285: {
! 1286: return hen_ttob(T,REMOTE_NF,0,REMOTE_VARS,M);
! 1287: }
! 1288:
! 1289: def r_ttob_gsl(L,M)
! 1290: {
! 1291: return cons(L[2],hen_ttob(L[0],REMOTE_NF,L[1],REMOTE_VARS,M));
! 1292: }
! 1293:
! 1294: def get_matrix()
! 1295: {
! 1296: REMOTE_MATRIX;
! 1297: }
! 1298: end$