Annotation of OpenXM/src/asir-contrib/testing/noro/module_syz.rr, Revision 1.8
1.1 noro 1: module newsyz;
2:
1.4 noro 3: localf module_syz, module_syz_old;
1.3 noro 4: localf simplify_syz, icont, mod, remove_cont,ordcheck;
5: localf complsb, complsb_sd, sortlsb, find_pos, find_pos, reduce, lres_setup, dpm_sort1, comp_pos;
6: localf fres,minres,sres,minsres,lres, create_base_ord, simplify_k, simplify_by_k, remove_k, remove_k1, extract_nonzero;
7: localf nonzero, phi, syz_check, renumber_pos, compress, compress_h;
1.8 ! noro 8: localf syz_check0,phi0,todpmlist,dpmlisttollist,comp_lex;
1.1 noro 9:
10: /* F : a list of (lists or polynomials),
11: V : a variable list, H >1=> over GF(H), H=0,1=> over Q
12: O : term order
13: return: [GS,G]
14: GS : a GB of syz(F) wrt [1,O] (POT), G: a GB of F wrt [1,O]
15: */
16:
1.3 noro 17: // return [dpmlist F,rank N]
18: def todpmlist(F,V)
1.1 noro 19: {
1.3 noro 20: K = length(F);
21: for ( I = 0; I < K; I++ ) if ( F[I] ) break;
22: if ( I == K ) return [];
23: if ( type(F[I]) <= 2 ) {
24: // F is a polynimial list
25: F = map(dp_ptod,F,V);
26: F = map(dpm_dptodpm,F,1);
27: N = 1;
28: } else if ( type(F[I]) == 9 ) {
29: // F is a dpoly list
30: F = map(dpm_dptodpm,F,1);
31: N = 1;
32: } else if ( type(F[I]) == 4 ) {
33: // F is a list of poly lists
34: N = length(F[0]);
35: F = map(dpm_ltod,F,V);
36: } else if ( type(F[I]) == 26 ) {
37: // F is a DPM list
38: for ( N = 0, T = F; T != []; T = cdr(T) ) {
39: for ( A = car(T); A; A = dpm_rest(A) ) {
40: N1 = dpm_hp(A);
41: if ( N1 > N ) N = N1;
42: }
43: }
44: } else {
45: error("todpmlist: arugument type is invalid.");
46: }
47: return [F,N];
48: }
49:
50: def module_syz(F,V,H,Ord)
51: {
52: if ( type(Weyl=getopt(weyl)) == -1 ) Weyl = 0;
53: if ( type(DP=getopt(dp)) == -1 ) DP = 0;
54: if ( type(F4=getopt(f4)) == -1 ) F4 = 0;
55: dp_ord(Ord);
56: K = length(F);
57: for ( I = 0; I < K; I++ ) if ( F[I] ) break;
58: if ( I == K ) return [[],[],[]];
59: L = todpmlist(F,V);
60: F = L[0]; N = L[1];
61: dp_ord([1,Ord]);
62: B = [];
63: for ( I = 0; I < K; I++ ) {
64: B = cons(F[I]+dpm_dptodpm(dp_ptod(1,V),N+I+1),B);
65: }
66: B = reverse(B);
1.4 noro 67: if ( H >= 2 ) {
68: // finite field
69: if ( Weyl )
70: G = nd_weyl_gr(B,V,H,[1,Ord]|dp=1);
71: else if ( F4 )
72: G = nd_f4(B,V,H,[1,Ord]|dp=1);
73: else
74: G = nd_gr(B,V,H,[1,Ord]|dp=1);
75: } else {
76: if ( Weyl )
77: G = nd_weyl_gr(B,V,0,[1,Ord]|dp=1,homo=H);
1.7 noro 78: else {
1.4 noro 79: Ind = 0;
80: while ( 1 ) {
1.7 noro 81: if ( F4 )
82: G = nd_f4_trace(B,V,H,-lprime(Ind),[1,Ord]|dp=1);
83: else
84: G = nd_gr_trace(B,V,H,-lprime(Ind),[1,Ord]|dp=1);
1.4 noro 85: if ( G ) break;
86: else Ind++;
87: }
1.7 noro 88: }
1.4 noro 89: }
1.3 noro 90: G0 = []; S0 = []; Gen0 = [];
91: for ( T = G; T != []; T = cdr(T) ) {
92: H = car(T);
93: if ( dpm_hp(H) > N ) {
94: S0 = cons(dpm_shift(H,N),S0);
95: } else {
96: L = dpm_split(H,N);
97: G0 = cons(L[0],G0);
98: Gen0 = cons(dpm_shift(L[1],N),Gen0);
1.1 noro 99: }
1.3 noro 100: }
101: #if 0
102: S0 = reverse(S0); G0 = reverse(G0); Gen0 = reverse(Gen0);
103: #endif
104: if ( !DP ) {
105: S0 = map(dpm_dtol,S0,V); G0 = map(dpm_dtol,G0,V); Gen0 = map(dpm_dtol,Gen0,V);
106: }
107: return [S0,G0,Gen0];
108: }
109:
1.4 noro 110: def module_syz_old(F,V,H,O)
111: {
112: Weyl = type(getopt(weyl)) != -1 ? 1 : 0;
113: K = length(F);
114: if ( type(F[0]) <= 2 ) {
115: for ( T = [], S = F; S != []; S = cdr(S) )
116: T = cons([car(S)],T);
117: F = reverse(T);
118: }
119: N = length(F[0]);
120: B = [];
121: for ( I = 0; I < K; I++ ) {
122: E = vector(N+K);
123: for ( J = 0; J < N; J++ ) E[J] = F[I][J];
124: E[N+I] = 1;
125: B = cons(vtol(E),B);
126: }
127: B = reverse(B);
128: if ( H >= 2 ) {
129: if ( Weyl )
130: G = nd_weyl_gr(B,V,H,[1,O]);
131: else
132: G = nd_gr(B,V,H,[1,O]);
133: } else {
134: if ( Weyl )
135: G = nd_weyl_gr_trace(B,V,H,-1,[1,O]);
136: else
137: G = nd_gr_trace(B,V,H,-1,[1,O]);
138: }
139: G0 = []; S0 = []; Gen0 = [];
140: for ( T = G; T != []; T = cdr(T) ) {
141: H = car(T);
142: for ( J = 0; J < N; J++ ) if ( H[J] ) break;
143: if ( J == N ) {
144: H1 = vector(K);
145: for ( J = 0; J < K; J++ ) H1[J] = H[N+J];
146: S0 = cons(vtol(H1),S0);
147: } else {
148: H1 = vector(N);
149: for ( J = 0; J < N; J++ ) H1[J] = H[J];
150: G0 = cons(vtol(H1),G0);
151: H1 = vector(K);
152: for ( J = 0; J < K; J++ ) H1[J] = H[N+J];
153: Gen0 = cons(vtol(H1),Gen0);
154: }
155: }
156: return [S0,G0,Gen0];
157: }
158:
1.3 noro 159: def fres(F,V,H,O)
160: {
161: if ( type(Weyl=getopt(weyl)) == -1 ) Weyl = 0;
162: if ( type(F4=getopt(f4)) == -1 ) F4 = 0;
163: if ( type(DP=getopt(dp)) == -1 ) DP = 0;
164: L = todpmlist(F,V); F = L[0];
165: R = [F];
166: while ( 1 ) {
167: if ( Weyl )
168: L = module_syz(car(R),V,H,O|weyl=1,dp=1,f4=F4);
169: else
170: L = module_syz(car(R),V,H,O|dp=1,f4=F4);
171: L = L[0];
172: L = ordcheck(L,V);
173: if ( L == [] ) {
1.5 noro 174: R = reverse(R);
1.3 noro 175: if ( DP ) return R;
176: else return map(dpmlisttollist,R,V);
177: }
178: R = cons(L,R);
179: }
180: }
181:
182: def minres(F,V,H,O)
183: {
184: if ( type(Weyl=getopt(weyl)) == -1 ) Weyl = 0;
185: if ( type(F4=getopt(f4)) == -1 ) F4 = 0;
186: if ( type(DP=getopt(dp)) == -1 ) DP = 0;
187: L = todpmlist(F,V); F = L[0];
188: R = [F];
189: while ( 1 ) {
190: if ( Weyl )
191: L = module_syz(car(R),V,H,O|weyl=1,dp=1,f4=F4);
192: else
193: L = module_syz(car(R),V,H,O|dp=1,f4=F4);
194: L = L[0];
195: L = ordcheck(L,V);
196: if ( L == [] ) break;
197: S = dpm_simplify_syz(L,R[0]);
198: if ( S[0] == [] && S[1] == [] ) {
199: R = cdr(R); break;
200: }
201: R = append([S[0],S[1]],cdr(R));
202: if ( R[0] == [] ) {
203: R = cdr(R); break;
204: }
205: }
1.5 noro 206: R = reverse(R);
1.3 noro 207: if ( DP ) return R;
208: else return map(dpmlisttollist,R,V);
209: }
210:
211: def dpmlisttollist(F,V)
212: {
213: return map(dpm_dtol,F,V);
214: }
215:
216: def sres(F,V,H,Ord)
217: {
218: if ( type(DP=getopt(dp)) == -1 ) DP = 0;
219: dpm_set_schreyer(0);
220: dp_ord(Ord);
221: K = length(F);
222: K = length(F);
223: for ( I = 0; I < K; I++ ) if ( F[I] ) break;
224: if ( I == K ) return [[],[],[]];
225: L = todpmlist(F,V);
226: F = L[0]; N = L[1];
1.6 noro 227: #if 0
1.3 noro 228: G = nd_gr(F,V,H,[0,Ord]|dp=1);
1.6 noro 229: #else
230: G = nd_gr_trace(F,V,H,1,[0,Ord]|dp=1);
231: #endif
1.3 noro 232: G = reverse(G);
233: R = [G];
234: dp_ord([0,Ord]);
235: while ( 1 ) {
236: S = dpm_schreyer_base(R[0]);
1.7 noro 237: if ( dp_gr_print() ) print(["length",length(S)]);
1.3 noro 238: if ( S == [] ) break;
239: else R = cons(S,R);
240: }
241: dp_ord([0,0]);
1.5 noro 242: R = reverse(R);
1.3 noro 243: if ( DP ) return R;
244: else return map(dpmlisttollist,R,V);
245: }
246:
247: def minsres(F,V,H,Ord)
248: {
249: if ( type(DP=getopt(dp)) == -1 ) DP = 0;
250: R = sres(F,V,H,Ord|dp=1);
1.5 noro 251: R = ltov(reverse(R));
1.3 noro 252: M = length(R);
253: for ( I = 0; I < M; I++ ) R[I] = map(dpm_sort,R[I]);
254: R = vtol(R);
255: for ( T = R, U = []; length(T) >= 2; ) {
256: Z = dpm_simplify_syz(T[0],T[1]);
257: U = cons(Z[0],U);
258: T = cons(Z[1],cdr(cdr(T)));
259: }
260: U = cons(T[0],U);
261: if ( DP ) return U;
262: else return map(dpmlisttollist,U,V);
263: }
264:
265: def ordcheck(D,V)
266: {
267: B=map(dpm_dtol,D,V);
268: dp_ord([1,0]);
269: D = map(dpm_sort,D);
270: D1 = map(dpm_ltod,B,V);
271: if ( D != D1 )
272: print("afo");
273: return D1;
274: }
275:
276: def complsb(A,B)
277: {
278: AL = A[3]; AD = A[4];
279: BL = B[3]; BD = B[4];
280: if ( AD > BD ) return 1;
281: else if ( AD < BD ) return -1;
282: else if ( AL < BL ) return 1;
283: else if ( AL > BL ) return -1;
284: else return 0;
285: }
286:
287: def complsb_sd(A,B)
288: {
289: AL = A[3]; AD = A[4];
290: BL = B[3]; BD = B[4];
291: if ( AD-AL > BD-BL ) return 1;
292: else if ( AD-AL < BD-BL ) return -1;
293: else if ( AL > BL ) return 1;
294: else if ( AL < BL ) return -1;
295: else return 0;
296: }
297:
298: def sortlsb(A)
299: {
300: S = [];
301: N = length(A);
302: for ( I = 0; I < N; I++ ) S = append(cdr(vtol(A[I])),S);
303: // S = qsort(S,newsyz.complsb_sd);
304: S = qsort(S,newsyz.complsb);
305: return S;
306: }
307: /* D=[in,sum,deg,lev,ind] */
308: /* B=[0 B1 B2 ...] */
309: /* C=[0 C1 C2 ...] */
310: /* H=[0 H1 H2 ...] */
311: /* Z=[0 Z1 Z2 ...] */
312: /* Zi=[0 L1 L2 ...] */
313: /* Lj=[I1,I2,...] */
314: /* One=<<0,...,0>> */
315:
316: extern Rtime,Stime,Ptime$
317:
318: // B is sorted wrt positon
319: def find_pos(HF,B)
320: {
321: Len = length(B);
322: Pos = dpm_hp(HF);
323: First = 0; Last = Len-1;
324: if ( B[First][0] == HF ) return B[First][2];
325: if ( B[Last][0] == HF ) return B[Last][2];
326: while ( 1 ) {
327: Mid = idiv(First+Last,2);
328: PosM = dpm_hp(B[Mid][0]);
329: if ( PosM == Pos ) break;
330: else if ( PosM < Pos )
331: First = Mid;
332: else
333: Last = Mid;
334: }
335: for ( I = Mid; I < Len && dpm_hp(B[I][0]) == Pos; I++ )
336: if ( HF == B[I][0] ) return B[I][2];
337: for ( I = Mid-1; I >= 1 && dpm_hp(B[I][0]) == Pos; I-- )
338: if ( HF == B[I][0] ) return B[I][2];
339: error("find_pos : cannot happen");
340: }
341:
342: def reduce(D,B,Bpos,C,H,Z,K,Kind,G,One,Top)
343: {
344: M = D[0]; Ind = D[2]; I = D[3];
345: if ( I == 1 ) {
346: C[I][Ind] = G[Ind];
347: dpm_insert_to_zlist(Z[1],dpm_hp(G[Ind]),Ind);
348: H[I][Ind] = G[Ind];
349: } else {
350: /* M is in F(I-1) => phi(M) is in F(I-2) */
351: /* reduction in F(I-2) */
352: T0 = time()[0];
353: dpm_set_schreyer_level(I-2);
354: CI = C[I-1]; ZI = Z[I-1]; BI = B[I-1]; BposI = Bpos[I-1];
355: Len = size(CI);
356: T=dpm_hc(M); EM=dpm_hp(M);
357: XiM = T*dpm_ht(CI[EM]);
358: EXiM = dpm_hp(XiM);
359: ZIE = ZI[EXiM];
360: for ( ZIE = ZI[EXiM]; ZIE != []; ZIE = cdr(ZIE) ) {
361: J = car(ZIE);
362: if ( J > EM && dpm_redble(XiM,dpm_ht(CI[J])) ) break;
363: }
364: Ptime += time()[0]-T0;
365: T0 = time()[0];
366: QR = dpm_sp_nf(CI,ZI,EM,J|top=Top);
367: Rtime += time()[0]-T0;
368: G = QR[0]; F = QR[1];
369: if ( F ) {
370: HF = dpm_ht(F); EF = dpm_hp(HF);
371: /* find HF in B[I-1] */
372: T0 = time()[0];
373: J = find_pos(HF,BposI);
374: Stime += time()[0]-T0;
375: /* F=Ret[0]*Ret[1] */
376: Ret = dpm_remove_cont(F);
377: CI[J] = Ret[1];
378: Kind[I-1][J] = 2;
379: dpm_insert_to_zlist(ZI,EF,J);
380: dpm_set_schreyer_level(I-1);
381: Tail = -Ret[0]*dpm_dptodpm(One,J);
382: G += Tail;
383: dpm_set_schreyer_level(0);
384: Ret = dpm_remove_cont(G); G = Ret[1];
385: C[I][Ind] = G;
386: Kind[I][Ind] = 1;
387: K[I] = cons([G,Tail],K[I]);
388: dpm_insert_to_zlist(Z[I],EM,Ind);
389: /* level <- 0 */
390: BI[J][3] = 0;
391: } else {
392: Kind[I][Ind] = 0;
393: Ret = dpm_remove_cont(G); G = Ret[1];
394: H[I][Ind] = G;
395: C[I][Ind] = G;
396: dpm_insert_to_zlist(Z[I],EM,Ind);
397: }
398: }
399: }
400:
1.8 ! noro 401: def comp_lex(A,B)
! 402: {
! 403: HA = dpm_hc(A); HB = dpm_hc(B);
! 404: if ( HA > HB ) return 1;
! 405: else if ( HA < HB ) return -1;
! 406: else return 0;
! 407: }
! 408:
1.3 noro 409: def lres_setup(F,V,H,Ord)
410: {
1.8 ! noro 411: if ( type(Lex=getopt(lex)) == -1 ) Lex = 0;
1.3 noro 412: dpm_set_schreyer(0);
413: dp_ord(Ord);
414: K = length(F);
415: if ( type(F[0]) <= 2 ) {
416: // F is a polynimial list
417: F = map(dp_ptod,F,V);
418: F = map(dpm_dptodpm,F,1);
419: N = 1;
420: } else if ( type(F[0]) == 9 ) {
421: // F is a dpoly list
422: F = map(dpm_dptodpm,F,1);
423: N = 1;
424: } else if ( type(F[0]) == 4 ) {
425: // F is a list of poly lists
426: N = length(F[0]);
427: F = map(dpm_ltod,F,V);
428: } else if ( type(F[0]) == 26 ) {
429: // F is a DPM list
430: for ( N = 0, T = F; T != []; T = cdr(T) ) {
431: for ( A = car(T); A; A = dpm_rest(A) ) {
432: N1 = dpm_hp(A);
433: if ( N1 > N ) N = N1;
434: }
435: }
436: } else {
437: error("lres_setup: arugument type is invalid.");
438: }
1.6 noro 439: dp_ord([0,Ord]);
440: F = map(dpm_sort,F);
441: #if 0
1.3 noro 442: G = nd_gr(F,V,H,[0,Ord]|dp=1);
1.6 noro 443: #else
444: G = nd_gr_trace(F,V,H,1,[0,Ord]|dp=1);
445: #endif
1.8 ! noro 446: if ( Lex ) {
! 447: dp_ord(2);
! 448: G = qsort(G,newsyz.comp_lex);
! 449: }
1.3 noro 450: G = reverse(G);
451: dp_ord([0,Ord]);
452: One = dp_ptod(1,V);
453: return [G,One];
454: }
455:
456: def dpm_sort1(L)
457: {
458: return [dpm_sort(L[0]),L[1]];
459: }
460:
461: def comp_pos(A,B)
462: {
463: PA = dpm_hp(A[0]); PB = dpm_hp(B[0]);
464: if ( PA > PB ) return 1;
465: else if ( PA < PB ) return -1;
466: else return 0;
467: }
468:
469: def lres(F,V,H,Ord)
470: {
471: T0 = time();
472: if ( type(Top=getopt(top)) == -1 ) Top = 0;
1.5 noro 473: if ( type(DP=getopt(dp)) == -1 ) DP = 0;
1.3 noro 474: if ( type(NoSimpK=getopt(nosimpk)) == -1 ) NoSimpK = 0;
475: if ( type(NoPreProj=getopt(nopreproj)) == -1 ) NoPreProj = 0;
1.8 ! noro 476: if ( type(Lex=getopt(lex)) == -1 ) Lex = 0;
1.3 noro 477: Rtime = Stime = Ptime = 0;
478: L = lres_setup(F,V,H,Ord);
479: G = L[0];
480: One = L[1];
1.8 ! noro 481: F = dpm_schreyer_frame(G|lex=Lex);
1.3 noro 482: G = ltov(cons(0,L[0]));
483: F = reverse(F);
484: F = ltov(F);
485: N = length(F);
486: for ( I = 0; I < N; I++ ) {
487: FI = F[I]; FI[0] = [];
488: FI = map(ltov,FI);
489: F[I] = FI;
490: }
491: R = sortlsb(F);
492: B = vector(N+1);
493: Bpos = vector(N+1);
494: C = vector(N+1);
495: H = vector(N+1);
496: Z = vector(N+1);
497: K = vector(N+1);
498: L = vector(N+1);
499: K = vector(N+1);
500: D = vector(N+1);
501: Kind = vector(N+1);
502:
503: for ( I = 1; I <= N; I++ ) {
504: FI = F[I-1]; Len = length(FI);
505: B[I] = FI;
506: T = vector(Len-1);
507: for ( J = 1; J < Len; J++ ) T[J-1] = FI[J];
508: Bpos[I] = qsort(T,newsyz.comp_pos);
509: C[I] = vector(Len);
510: H[I] = vector(Len);
511: Kind[I] = vector(Len);
512: K[I] = [];
513: Max = 0;
514: for ( J = 1; J < Len; J++ )
515: if ( (Pos = dpm_hp(FI[J][0])) > Max ) Max = Pos;
516: Z[I] = ZI = vector(Max+1);
517: for ( J = 1; J <= Max; J++ ) ZI[J] = [];
518: }
519: T1 = time(); Ftime = T1[0]-T0[0];
520: R = ltov(R); Len = length(R);
1.7 noro 521: if ( dp_gr_print() ) print(["Len",Len]);
1.3 noro 522: for ( I = 0, NF = 0; I < Len; I++ ) {
1.7 noro 523: if ( dp_gr_print() ) {
524: if ( !((I+1)%100) ) print(".",2);
525: if ( !((I+1)%10000) ) print(I+1);
526: }
1.3 noro 527: if ( !R[I][3] ) continue;
528: NF++;
529: reduce(R[I],B,Bpos,C,H,Z,K,Kind,G,One,Top);
530: }
1.7 noro 531: if ( dp_gr_print() ) {
532: print("");
533: print(["NF",NF]);
534: }
1.3 noro 535: T0 = time();
536: dpm_set_schreyer_level(0);
537: D[1] = map(dpm_sort,H[1]);
538: for ( I = 2; I <= N; I++ ) {
539: // HTT = [Head,Tab,TailTop]
540: HTT = create_base_ord(K[I],length(Kind[I-1]));
541: Head = HTT[0]; Tab = HTT[1]; TailTopPos = HTT[2];
542: if ( !NoPreProj )
543: Tab = map(remove_k,map(dpm_sort,Tab),Kind[I-1]);
544: else
545: Tab = map(dpm_sort,Tab);
546: TailTop = dpm_dptodpm(One,TailTopPos);
547: if ( !NoSimpK ) {
1.7 noro 548: if ( dp_gr_print() ) print("simplify_k "+rtostr(I)+"...",2);
1.3 noro 549: simplify_k(Head,Tab,TailTop,One);
1.7 noro 550: if ( dp_gr_print() ) print("done");
1.3 noro 551: }
552: HI = map(remove_k,map(dpm_sort,H[I]),Kind[I-1]);
553: Len = length(HI);
1.7 noro 554: if ( dp_gr_print() ) print("simplify_by_k "+rtostr(I)+"...",2);
1.3 noro 555: D[I] = vector(Len);
556: for ( J = 0; J < Len; J++ ) {
557: D[I][J] = simplify_by_k(HI[J],Tab,TailTop,One);
558: if ( NoPreProj )
559: D[I][J] = remove_k(D[I][J],Kind[I-1]);
560: }
1.7 noro 561: if ( dp_gr_print() ) print("done");
1.3 noro 562: }
563: dp_ord([1,0]);
564: T1 = time();
1.7 noro 565: if ( dp_gr_print() ) print(["Frame",Ftime,"Prep",Ptime,"Reduce",Rtime,"Search",Stime,"Minimalize",T1[0]-T0[0]]);
1.3 noro 566: // return [C,H,K,Kind,D];
567: D = compress_h(D);
1.5 noro 568: if ( DP ) return D;
569: else return vtol(map(dpmlisttollist,D,V));
1.3 noro 570: }
571:
572: def create_base_ord(K,N)
573: {
574: Tab = vector(N);
575: Ks = [];
576: for ( T = K; T != []; T = cdr(T) ) {
577: Ks = cons(I=dpm_hp(car(T)[1]),Ks);
578: Tab[I] = car(T)[0];
579: }
580: Others = [];
581: for ( I = N-1; I >= 1; I-- )
582: if ( !Tab[I] ) Others = cons(I,Others);
583: Ks = reverse(Ks);
584: dp_ord([4,append(Ks,Others),0]);
585: return [ltov(Ks),Tab,car(Others)];
586: }
587:
588: /* Head = [i1,i2,...], Tab[i1] = K[0], Tab[i2] = K[1],... */
589:
590: def simplify_k(Head,Tab,TailTop,One)
591: {
592: N = length(Tab);
593: Len = length(Head);
594: for ( I = Len-1; I >= 0; I-- ) {
595: M = 1;
596: R = Tab[Head[I]];
597: H = dpm_hm(R); R = dpm_rest(R);
598: while ( R && dpm_dptodpm(One,dpm_hp(R)) > TailTop ) {
599: Pos = dpm_hp(R); Red = Tab[Pos];
600: CRed = dp_hc(dpm_hc(Red));
601: CR = dpm_extract(R,Pos);
602: L = dpm_remove_cont(CRed*R-CR*Red);
603: R = L[1]; M *= CRed/L[0];
604: }
605: Tab[Head[I]] = M*H+R;
606: }
607: }
608:
609: def simplify_by_k(F,Tab,TailTop,One)
610: {
611: R = F;
612: M = 1;
613: while ( R && dpm_dptodpm(One,dpm_hp(R)) > TailTop ) {
614: Pos = dpm_hp(R); Red = Tab[Pos];
615: CRed = dp_hc(dpm_hc(Red));
616: CR = dpm_extract(R,Pos);
617: L = dpm_remove_cont(CRed*R-CR*Red);
618: R = L[1]; M *= CRed/L[0];
619: }
620: return (1/M)*R;
621: }
622:
623: /* Kind[I]=0 => phi(ei)\in H, Kind[I]=1 => phi(ei)\in K */
624: def remove_k(F,Kind)
625: {
626: R = [];
627: for ( T = F; T; T = dpm_rest(T) )
628: if ( Kind[dpm_hp(T)] != 1 ) R = cons(dpm_hm(T),R);
629: for ( S = 0; R != []; R = cdr(R) )
630: S += car(R);
631: return S;
632: }
633:
634: def remove_k1(F,Kind)
635: {
636: return [remove_k(F[0],Kind),F[1]];
637: }
638:
639: def extract_nonzero(A)
640: {
641: N = length(A);
642: R = [];
643: for ( C = I = 0; I < N; I++ )
644: if ( A[I] ) R=cons(A[I],R);
645: return reverse(R);;
646: }
647:
648: def nonzero(A)
649: {
650: N = length(A);
651: for ( C = I = 0; I < N; I++ )
652: if ( A[I] ) C++;
653: return C;
654: }
655:
656: def phi(C,F)
657: {
658: R = 0;
659: for ( T = F; T; T = dpm_rest(T) ) {
660: Coef = dpm_hc(T); Pos = dpm_hp(T);
1.6 noro 661: R += Coef*C[Pos-1];
1.3 noro 662: }
663: return R;
1.1 noro 664: }
665:
1.3 noro 666: def syz_check(H,I)
1.1 noro 667: {
1.3 noro 668: HI = H[I];
669: // dpm_set_schreyer_level(I-1);
670: HI1 = H[I+1];
671: Len = size(HI1)[0];
672: for ( J = 1; J < Len; J++ ) {
673: F = HI1[J];
674: R = phi(HI,F);
675: if ( R ) print([J,"NG"]);
676: }
1.1 noro 677: }
678:
1.3 noro 679: // for compressed C
680: def phi0(C,F)
1.1 noro 681: {
1.3 noro 682: R = 0;
683: for ( T = F; T; T = dpm_rest(T) ) {
684: Coef = dpm_hc(T); Pos = dpm_hp(T);
685: R += Coef*C[Pos-1];
686: }
687: return R;
1.1 noro 688: }
689:
1.3 noro 690: // for compressed H
691: def syz_check0(H,I)
1.1 noro 692: {
1.3 noro 693: HI = H[I];
694: // dpm_set_schreyer_level(I-1);
695: HI1 = H[I+1];
696: Len = length(HI1);
697: for ( J = 0; J < Len; J++ ) {
698: F = HI1[J];
699: R = phi0(HI,F);
700: if ( R ) print([J,"NG"]);
701: }
1.1 noro 702: }
703:
1.3 noro 704: // renumber position
705: def renumber_pos(F,Tab)
1.1 noro 706: {
1.3 noro 707: L = [];
708: for ( T = F; T; T = dpm_rest(T) )
709: L = cons(dpm_hm(T),L);
710: R = 0;
711: for ( T = L; T != []; T = cdr(T) )
712: R += dpm_dptodpm(dpm_hc(car(T)),Tab[dpm_hp(car(T))]);
713: return R;
1.1 noro 714: }
715:
1.3 noro 716: // compress H1 and renumber H
717: def compress(H,H1)
1.1 noro 718: {
1.3 noro 719: // create index table for H1
720: L1 = length(H1);
721: Tmp = vector(L1);
722: Tab = vector(L1);
723: for ( I = 1, J = 1; I < L1; I++ )
724: if ( H1[I] ) {
725: Tab[I] = J;
726: Tmp[J++] = H1[I];
727: }
728: NH1 = vector(J);
729: for ( I = 1; I < J; I++ )
730: NH1[I] = Tmp[I];
731: if ( H )
732: H = map(renumber_pos,H,Tab);
733: return [H,NH1];
1.1 noro 734: }
735:
1.3 noro 736: def compress_h(H)
1.1 noro 737: {
1.3 noro 738: // H = [0,H[1],...,H[L-1]]
739: L = length(H);
740: NH = vector(L-1);
741: H1 = H[1];
742: for ( I = 2; I < L; I++ ) {
743: R = compress(H[I],H1);
744: H1 = R[0];
745: NH[I-2] = cdr(vtol(R[1]));
746: }
747: R = compress(0,H1);
748: NH[L-2] = cdr(vtol(R[1]));
749: return NH;
1.1 noro 750: }
1.3 noro 751: endmodule$
1.1 noro 752: end$
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>