![[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.2 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2000/lib/gr,v 1.1.1.1 1999/12/03 07:39:11 noro Exp $ */
1.1 noro 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:
1.2 ! noro 1261: def dgr(G,V,O)
1.1 noro 1262: {
1.2 ! noro 1263: P = getopt(proc);
! 1264: if ( type(P) == -1 )
! 1265: return gr(G,V,O);
1.1 noro 1266: P0 = P[0]; P1 = P[1]; P = [P0,P1];
1.2 ! noro 1267: map(ox_reset,P);
! 1268: ox_cmo_rpc(P0,"dp_gr_main",G,V,0,1,O);
! 1269: ox_cmo_rpc(P1,"dp_gr_main",G,V,1,1,O);
! 1270: map(ox_push_cmd,P,262); /* 262 = OX_popCMO */
! 1271: F = ox_select(P);
! 1272: R = ox_get(F[0]);
! 1273: if ( F[0] == P0 ) {
! 1274: Win = "nonhomo";
! 1275: Lose = P1;
! 1276: } else {
! 1277: Win = "nhomo";
! 1278: Lose = P0;
! 1279: }
! 1280: ox_reset(Lose);
! 1281: return [Win,R];
1.1 noro 1282: }
1283:
1284: /* functions for rpc */
1285:
1286: def register_matrix(M)
1287: {
1288: REMOTE_MATRIX = M; return 0;
1289: }
1290:
1291: def register_nfv(L)
1292: {
1293: REMOTE_NF = L[0]; REMOTE_VARS = L[1]; return 0;
1294: }
1295:
1296: def r_ttob(T,M)
1297: {
1298: return hen_ttob(T,REMOTE_NF,0,REMOTE_VARS,M);
1299: }
1300:
1301: def r_ttob_gsl(L,M)
1302: {
1303: return cons(L[2],hen_ttob(L[0],REMOTE_NF,L[1],REMOTE_VARS,M));
1304: }
1305:
1306: def get_matrix()
1307: {
1308: REMOTE_MATRIX;
1309: }
1310: end$