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