![[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.6 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/lib/gr,v 1.5 2000/08/21 08:31:41 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);
257: print(["DN",B[1]]);
258: }
259: return SL;
260: }
261:
262: def hen_ttob_gsl(LHS,RHS,TERMS,M)
263: {
264: LDN = LHS[1]; RDN = RHS[1]; LCM = ilcm(LDN,RDN);
265: L1 = idiv(LCM,LDN); R1 = idiv(LCM,RDN);
266: T0 = time()[0];
267: S = henleq_gsl(RHS[0],LHS[0]*L1,M);
268: print(["henleq_gsl",time()[0]-T0]);
269: N = length(TERMS);
270: return [S[0],S[1]*R1];
271: }
272:
273: def gather_terms(GM,W,M,NPOSV)
274: {
275: N = length(W); V0 = W[NPOSV];
276: for ( T = GM; T != []; T = cdr(T) ) {
277: if ( vars(car(T)) == [V0] )
278: break;
279: }
280: U = car(T); DU = diff(U,V0);
281: R = tpoly(cdr(p_terms(U,W,2)));
282: for ( I = 0; I < N; I++ ) {
283: if ( I == NPOSV )
284: continue;
285: V1 = W[I];
286: for ( T = GM; T != []; T = cdr(T) )
287: if ( member(V1,vars(car(T))) )
288: break;
289: P = car(T);
290: R += tpoly(p_terms(srem(DU*coef(P,0,V1),U,M),W,2));
291: }
292: return p_terms(R,W,2);
293: }
294:
295: def tpoly(L)
296: {
297: for ( R = 0; L != []; L = cdr(L) )
298: R += car(L);
299: return R;
300: }
301:
302: def dptov(P,W,MB)
303: {
304: N = size(W)[0];
305: for ( I = 0; I < N; I++ )
306: W[I] = 0;
307: for ( I = 0, S = MB; P; P = dp_rest(P) ) {
308: HM = dp_hm(P); C = dp_hc(HM); T = dp_ht(HM);
309: for ( ; T != car(S); S = cdr(S), I++ );
310: W[I] = C;
311: I++; S = cdr(S);
312: }
313: }
314:
315: def tolex_main(V,O,NF,GM,M,MB)
316: {
317: DIM = length(MB);
318: DV = newvect(DIM);
319: for ( T = GM, SL = [], LCM = 1; T != []; T = cdr(T) ) {
320: S = p_terms(car(T),V,2);
321: dp_ord(O); RHS = termstomat(NF,map(dp_ptod,cdr(S),V),MB,M);
322: dp_ord(0); NHT = nf_tab_gsl(dp_ptod(LCM*car(S),V),NF);
323: dptov(NHT[0],DV,MB);
324: dp_ord(O); B = hen_ttob_gsl([DV,NHT[1]],RHS,cdr(S),M);
325: if ( !B )
326: return 0;
327: Len = length(S);
328: LCM *= B[1];
329: for ( U = LCM*car(S), I = 1; I < Len; I++ )
330: U += B[0][I-1]*S[I];
331: R = ptozp(U);
332: SL = cons(R,SL);
333: print(["DN",B[1]]);
334: }
335: return SL;
336: }
337:
338: def reduce_dn(L)
339: {
340: NM = L[0]; DN = L[1]; V = vars(NM);
341: T = remove_cont([dp_ptod(NM,V),DN]);
342: return [dp_dtop(T[0],V),T[1]];
343: }
344:
345: /* a function for computation of minimal polynomial */
346:
347: def minipoly(G0,V,O,P,V0)
348: {
349: if ( !zero_dim(hmlist(G0,V,O),V,O) )
350: error("tolex : ideal is not zero-dimensional!");
351:
352: G1 = cons(V0-P,G0);
353: O1 = [[0,1],[O,length(V)]];
354: V1 = cons(V0,V);
355: W = append(V,[V0]);
356:
357: N = length(V1);
358: dp_ord(O1);
359: HM = hmlist(G1,V1,O1);
360: MB = dp_mbase(map(dp_ptod,HM,V1));
361: dp_ord(O);
362:
363: for ( J = 0; ; J++ ) {
364: M = lprime(J);
365: if ( !valid_modulus(HM,M) )
366: continue;
367: MP = minipolym(G0,V,O,P,V0,M);
368: for ( D = deg(MP,V0), TL = [], J = 0; J <= D; J++ )
369: TL = cons(V0^J,TL);
370: NF = gennf(G1,TL,V1,O1,V0,1)[0];
371: R = tolex_main(V1,O1,NF,[MP],M,MB);
372: return R[0];
373: }
374: }
375:
376: /* subroutines */
377:
378: def gennf(G,TL,V,O,V0,FLAG)
379: {
380: N = length(V); Len = length(G); dp_ord(O); PS = newvect(Len);
381: for ( I = 0, T = G, HL = []; T != []; T = cdr(T), I++ ) {
382: PS[I] = dp_ptod(car(T),V); HL = cons(dp_ht(PS[I]),HL);
383: }
384: for ( I = 0, DTL = []; TL != []; TL = cdr(TL) )
385: DTL = cons(dp_ptod(car(TL),V),DTL);
386: for ( I = Len - 1, GI = []; I >= 0; I-- )
387: GI = cons(I,GI);
388: T = car(DTL); DTL = cdr(DTL);
389: H = [nf(GI,T,T,PS)];
390:
391: USE_TAB = (FLAG != 0);
392: if ( USE_TAB ) {
393: T0 = time()[0];
394: MB = dp_mbase(HL); DIM = length(MB);
395: U = dp_ptod(V0,V);
396: UTAB = newvect(DIM);
397: for ( I = 0; I < DIM; I++ ) {
398: UTAB[I] = [MB[I],remove_cont(dp_true_nf(GI,U*MB[I],PS,1))];
399: if ( dp_gr_print() )
400: print(".",2);
401: }
402: print("");
403: TTAB = time()[0]-T0;
404: }
405:
406: T0 = time()[0];
407: for ( LCM = 1; DTL != []; ) {
408: if ( dp_gr_print() )
409: print(".",2);
410: T = car(DTL); DTL = cdr(DTL);
411: if ( L = search_redble(T,H) ) {
412: DD = dp_subd(T,L[1]);
413: if ( USE_TAB && (DD == U) ) {
414: NF = nf_tab(L[0],UTAB);
415: NF = [NF[0],dp_hc(L[1])*NF[1]*T];
416: } else
417: NF = nf(GI,L[0]*dp_subd(T,L[1]),dp_hc(L[1])*T,PS);
418: } else
419: NF = nf(GI,T,T,PS);
420: NF = remove_cont(NF);
421: H = cons(NF,H);
422: LCM = ilcm(LCM,dp_hc(NF[1]));
423: }
424: TNF = time()[0]-T0;
425: if ( dp_gr_print() )
426: print("gennf(TAB="+rtostr(TTAB)+" NF="+rtostr(TNF)+")");
427: return [[map(adj_dn,H,LCM),LCM],PS,GI];
428: }
429:
430: def adj_dn(P,D)
431: {
432: return [(idiv(D,dp_hc(P[1])))*P[0],dp_ht(P[1])];
433: }
434:
435: def hen_ttob(T,NF,LHS,V,MOD)
436: {
437: if ( length(T) == 1 )
438: return car(T);
439: T0 = time()[0]; M = etom(leq_nf(T,NF,LHS,V)); TE = time()[0] - T0;
440: T0 = time()[0]; U = henleq(M,MOD); TH = time()[0] - T0;
441: if ( dp_gr_print() ) {
442: print("(etom="+rtostr(TE)+" hen="+rtostr(TH)+")");
443: }
444: return U ? vtop(T,U,LHS) : 0;
445: }
446:
447: def vtop(S,L,GSL)
448: {
449: U = L[0]; H = L[1];
450: if ( GSL ) {
451: for ( A = 0, I = 0; S != []; S = cdr(S), I++ )
452: A += U[I]*car(S);
453: return [A,H];
454: } else {
455: for ( A = H*car(S), S = cdr(S), I = 0; S != []; S = cdr(S), I++ )
456: A += U[I]*car(S);
457: return ptozp(A);
458: }
459: }
460:
461: def leq_nf(TL,NF,LHS,V)
462: {
463: TLen = length(NF);
464: T = newvect(TLen); M = newvect(TLen);
465: for ( I = 0; I < TLen; I++ ) {
466: T[I] = dp_ht(NF[I][1]);
467: M[I] = dp_hc(NF[I][1]);
468: }
469: Len = length(TL); INDEX = newvect(Len); COEF = newvect(Len);
470: for ( L = TL, J = 0; L != []; L = cdr(L), J++ ) {
471: D = dp_ptod(car(L),V);
472: for ( I = 0; I < TLen; I++ )
473: if ( D == T[I] )
474: break;
475: INDEX[J] = I; COEF[J] = strtov("u"+rtostr(J));
476: }
477: if ( !LHS ) {
478: COEF[0] = 1; NM = 0; DN = 1;
479: } else {
480: NM = LHS[0]; DN = LHS[1];
481: }
482: for ( J = 0, S = -NM; J < Len; J++ ) {
483: DNJ = M[INDEX[J]];
484: GCD = igcd(DN,DNJ); CS = DNJ/GCD; CJ = DN/GCD;
485: S = CS*S + CJ*NF[INDEX[J]][0]*COEF[J];
486: DN *= CS;
487: }
488: for ( D = S, E = []; D; D = dp_rest(D) )
489: E = cons(dp_hc(D),E);
490: BOUND = LHS ? 0 : 1;
491: for ( I = Len - 1, W = []; I >= BOUND; I-- )
492: W = cons(COEF[I],W);
493: return [E,W];
494: }
495:
496: def nf_tab(F,TAB)
497: {
498: for ( NM = 0, DN = 1, I = 0; F; F = dp_rest(F) ) {
499: T = dp_ht(F);
500: for ( ; TAB[I][0] != T; I++);
501: NF = TAB[I][1]; N = NF[0]; D = NF[1];
502: G = igcd(DN,D); DN1 = idiv(DN,G); D1 = idiv(D,G);
503: NM = D1*NM + DN1*dp_hc(F)*N; DN *= D1;
504: }
505: return [NM,DN];
506: }
507:
508: def nf_tab_gsl(A,NF)
509: {
510: DN = NF[1];
511: NF = NF[0];
512: TLen = length(NF);
513: for ( R = 0; A; A = dp_rest(A) ) {
514: HM = dp_hm(A); C = dp_hc(HM); T = dp_ht(HM);
515: for ( I = 0; I < TLen; I++ )
516: if ( NF[I][1] == T )
517: break;
518: R += C*NF[I][0];
519: }
520: return remove_cont([R,DN]);
521: }
522:
523: def redble(D1,D2,N)
524: {
525: for ( I = 0; I < N; I++ )
526: if ( D1[I] > D2[I] )
527: break;
528: return I == N ? 1 : 0;
529: }
530:
531: def tolexm(G,V,O,W,M)
532: {
533: N = length(V); Len = length(G);
534: dp_ord(O); setmod(M); PS = newvect(Len);
535: for ( I = 0, T = G; T != []; T = cdr(T), I++ )
536: PS[I] = dp_mod(dp_ptod(car(T),V),M,[]);
537: for ( I = Len-1, HL = []; I >= 0; I-- )
538: HL = cons(dp_ht(PS[I]),HL);
539: G2 = tolexm_main(PS,HL,V,W,M,ZeroDim);
540: L = map(dp_dtop,G2,V);
541: return L;
542: }
543:
544: def tolexm_main(PS,HL,V,W,M,FLAG)
545: {
546: N = length(W); D = newvect(N); Len = size(PS)[0];
547: for ( I = Len - 1, GI = []; I >= 0; I-- )
548: GI = cons(I,GI);
549: MB = dp_mbase(HL); DIM = length(MB);
550: U = dp_mod(dp_ptod(W[N-1],V),M,[]);
551: UTAB = newvect(DIM);
552: for ( I = 0; I < DIM; I++ ) {
553: if ( dp_gr_print() )
554: print(".",2);
555: UTAB[I] = [MB[I],dp_nf_mod(GI,U*dp_mod(MB[I],M,[]),PS,1,M)];
556: }
557: print("");
558: T = dp_mod(dp_ptod(dp_dtop(dp_vtoe(D),W),V),M,[]);
559: H = G = [[T,T]];
560: DL = []; G2 = [];
561: TNF = 0;
562: while ( 1 ) {
563: if ( dp_gr_print() )
564: print(".",2);
565: S = nextm(D,DL,N);
566: if ( !S )
567: break;
568: T = dp_mod(dp_ptod(dp_dtop(dp_vtoe(D),W),V),M,[]);
569: T0 = time()[0];
570: if ( L = search_redble(T,H) ) {
571: DD = dp_mod(dp_subd(T,L[1]),M,[]);
572: if ( DD == U )
573: NT = dp_nf_tab_mod(L[0],UTAB,M);
574: else
575: NT = dp_nf_mod(GI,L[0]*DD,PS,1,M);
576: } else
577: NT = dp_nf_mod(GI,T,PS,1,M);
578: TNF += time()[0] - T0;
579: H = cons([NT,T],H);
580: T0 = time()[0];
581: L = dp_lnf_mod([NT,T],G,M); N1 = L[0]; N2 = L[1];
582: TLNF += time()[0] - T0;
583: if ( !N1 ) {
584: G2 = cons(N2,G2);
585: if ( FLAG == MiniPoly )
586: break;
587: D1 = newvect(N);
588: for ( I = 0; I < N; I++ )
589: D1[I] = D[I];
590: DL = cons(D1,DL);
591: } else
592: G = insert(G,L);
593: }
594: if ( dp_gr_print() )
595: print("tolexm(nfm="+rtostr(TNF)+" lnfm="+rtostr(TLNF)+")");
596: return G2;
597: }
598:
599: def minipolym(G,V,O,P,V0,M)
600: {
601: N = length(V); Len = length(G);
602: dp_ord(O); setmod(M); PS = newvect(Len);
603: for ( I = 0, T = G; T != []; T = cdr(T), I++ )
604: PS[I] = dp_mod(dp_ptod(car(T),V),M,[]);
605: for ( I = Len-1, HL = []; I >= 0; I-- )
606: HL = cons(dp_ht(PS[I]),HL);
607: for ( I = Len - 1, GI = []; I >= 0; I-- )
608: GI = cons(I,GI);
609: MB = dp_mbase(HL); DIM = length(MB); UT = newvect(DIM);
610: U = dp_mod(dp_ptod(P,V),M,[]);
611: for ( I = 0; I < DIM; I++ )
612: UT[I] = [MB[I],dp_nf_mod(GI,U*dp_mod(MB[I],M,[]),PS,1,M)];
613: T = dp_mod(<<0>>,M,[]); TT = dp_mod(dp_ptod(1,V),M,[]);
614: G = H = [[TT,T]]; TNF = TLNF = 0;
615: for ( I = 1; ; I++ ) {
616: T = dp_mod(<<I>>,M,[]);
617: T0 = time()[0]; NT = dp_nf_tab_mod(H[0][0],UT,M); TNF += time()[0] - T0;
618: H = cons([NT,T],H);
619: T0 = time()[0]; L = dp_lnf_mod([NT,T],G,M); TLNF += time()[0] - T0;
620: if ( !L[0] ) {
621: if ( dp_gr_print() ) print(["nfm",TNF,"lnfm",TLNF]);
622: return dp_dtop(L[1],[V0]);
623: } else
624: G = insert(G,L);
625: }
626: }
627:
628: def nextm(D,DL,N)
629: {
630: for ( I = N-1; I >= 0; ) {
631: D[I]++;
632: for ( T = DL; T != []; T = cdr(T) )
633: if ( car(T) == D )
634: return 1;
635: else if ( redble(car(T),D,N) )
636: break;
637: if ( T != [] ) {
638: for ( J = N-1; J >= I; J-- )
639: D[J] = 0;
640: I--;
641: } else
642: break;
643: }
644: if ( I < 0 )
645: return 0;
646: else
647: return 1;
648: }
649:
650: def search_redble(T,G)
651: {
652: for ( ; G != []; G = cdr(G) )
653: if ( dp_redble(T,car(G)[1]) )
654: return car(G);
655: return 0;
656: }
657:
658: def insert(G,A)
659: {
660: if ( G == [] )
661: return [A];
662: else if ( dp_ht(car(A)) > dp_ht(car(car(G))) )
663: return cons(A,G);
664: else
665: return cons(car(G),insert(cdr(G),A));
666: }
667:
668: #if 0
669: def etom(L) {
670: E = L[0]; W = L[1];
671: LE = length(E); LW = length(W);
672: M = newmat(LE,LW+1);
673: for(J=0;J<LE;J++) {
674: for ( T = E[J]; T && (type(T) == 2); )
675: for ( V = var(T), I = 0; I < LW; I++ )
676: if ( V == W[I] ) {
677: M[J][I] = coef(T,1,V);
678: T = coef(T,0,V);
679: }
680: M[J][LW] = T;
681: }
682: return M;
683: }
684: #endif
685:
686: def etom(L) {
687: E = L[0]; W = L[1];
688: LE = length(E); LW = length(W);
689: M = newmat(LE,LW+1);
690: for(J=0;J<LE;J++) {
691: for ( I = 0, T = E[J]; I < LW; I++ ) {
692: M[J][I] = coef(T,1,W[I]); T = coef(T,0,W[I]);
693: }
694: M[J][LW] = T;
695: }
696: return M;
697: }
698:
699: def calcb_old(M) {
700: N = 2*M;
701: T = gr_sqrt(N);
702: if ( T^2 <= N && N < (T+1)^2 )
703: return idiv(T,2);
704: else
705: error("afo");
706: }
707:
708: def calcb_special(PK,P,K) { /* PK = P^K */
709: N = 2*PK;
710: T = sqrt_special(N,2,P,K);
711: if ( T^2 <= N && N < (T+1)^2 )
712: return idiv(T,2);
713: else
714: error("afo");
715: }
716:
717: def sqrt_special(A,C,M,K) { /* A = C*M^K */
718: L = idiv(K,2); B = M^L;
719: if ( K % 2 )
720: C *= M;
721: D = 2^K; X = idiv((gr_sqrt(C*D^2)+1)*B,D)+1;
722: while ( 1 )
723: if ( (Y = X^2) <= A )
724: return X;
725: else
726: X = idiv(A + Y,2*X);
727: }
728:
729: def gr_sqrt(A) {
730: for ( J = 0, T = A; T >= 2^27; J++ ) {
731: T = idiv(T,2^27)+1;
732: }
733: for ( I = 0; T >= 2; I++ ) {
734: S = idiv(T,2);
735: if ( T = S+S )
736: T = S;
737: else
738: T = S+1;
739: }
740: X = (2^27)^idiv(J,2)*2^idiv(I,2);
741: while ( 1 ) {
742: if ( (Y=X^2) < A )
743: X += X;
744: else if ( Y == A )
745: return X;
746: else
747: break;
748: }
749: while ( 1 )
750: if ( (Y = X^2) <= A )
751: return X;
752: else
753: X = idiv(A + Y,2*X);
754: }
755:
756: #define ABS(a) ((a)>=0?(a):(-a))
757:
758: def inttorat_asir(C,M,B)
759: {
760: if ( M < 0 )
761: M = -M;
762: C %= M;
763: if ( C < 0 )
764: C += M;
765: U1 = 0; U2 = M; V1 = 1; V2 = C;
766: while ( V2 >= B ) {
767: L = iqr(U2,V2); Q = L[0]; R2 = L[1];
768: R1 = U1 - Q*V1;
769: U1 = V1; U2 = V2;
770: V1 = R1; V2 = R2;
771: }
772: if ( ABS(V1) >= B )
773: return 0;
774: else
775: if ( V1 < 0 )
776: return [-V2,-V1];
777: else
778: return [V2,V1];
779: }
780:
781: def intvtoratv(V,M,B) {
782: if ( !B )
783: B = 1;
784: N = size(V)[0];
785: W = newvect(N);
786: if ( ITOR_FAIL >= 0 ) {
787: if ( V[ITOR_FAIL] ) {
788: T = inttorat(V[ITOR_FAIL],M,B);
789: if ( !T ) {
790: if ( dp_gr_print() ) {
791: print("F",2);
792: }
793: return 0;
794: }
795: }
796: }
797: for ( I = 0, DN = 1; I < N; I++ )
798: if ( V[I] ) {
799: T = inttorat((V[I]*DN) % M,M,B);
800: if ( !T ) {
801: ITOR_FAIL = I;
802: if ( dp_gr_print() ) {
803: #if 0
804: print("intvtoratv : failed at I = ",0); print(ITOR_FAIL);
805: #endif
806: print("F",2);
807: }
808: return 0;
809: } else {
810: for( J = 0; J < I; J++ )
811: W[J] *= T[1];
812: W[I] = T[0]; DN *= T[1];
813: }
814: }
815: return [W,DN];
816: }
817:
818: def nf(B,G,M,PS)
819: {
820: for ( D = 0; G; ) {
821: for ( U = 0, L = B; L != []; L = cdr(L) ) {
822: if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
823: GCD = igcd(dp_hc(G),dp_hc(R));
824: CG = idiv(dp_hc(R),GCD); CR = idiv(dp_hc(G),GCD);
825: U = CG*G-dp_subd(G,R)*CR*R;
826: if ( !U )
827: return [D,M];
828: D *= CG; M *= CG;
829: break;
830: }
831: }
832: if ( U )
833: G = U;
834: else {
835: D += dp_hm(G); G = dp_rest(G);
836: }
837: }
838: return [D,M];
839: }
840:
841: def remove_cont(L)
842: {
843: if ( type(L[1]) == 1 ) {
844: T = remove_cont([L[0],L[1]*<<0>>]);
845: return [T[0],dp_hc(T[1])];
846: } else if ( !L[0] )
847: return [0,dp_ptozp(L[1])];
848: else if ( !L[1] )
849: return [dp_ptozp(L[0]),0];
850: else {
851: A0 = dp_ptozp(L[0]); A1 = dp_ptozp(L[1]);
852: C0 = idiv(dp_hc(L[0]),dp_hc(A0)); C1 = idiv(dp_hc(L[1]),dp_hc(A1));
853: GCD = igcd(C0,C1); M0 = idiv(C0,GCD); M1 = idiv(C1,GCD);
854: return [M0*A0,M1*A1];
855: }
856: }
857:
858: def union(A,B)
859: {
860: for ( T = B; T != []; T = cdr(T) )
861: A = union1(A,car(T));
862: return A;
863: }
864:
865: def union1(A,E)
866: {
867: if ( A == [] )
868: return [E];
869: else if ( car(A) == E )
870: return A;
871: else
872: return cons(car(A),union1(cdr(A),E));
873: }
874:
875: def setminus(A,B) {
876: for ( T = reverse(A), R = []; T != []; T = cdr(T) ) {
877: for ( S = B, M = car(T); S != []; S = cdr(S) )
878: if ( car(S) == M )
879: break;
880: if ( S == [] )
881: R = cons(M,R);
882: }
883: return R;
884: }
885:
886: def member(A,L) {
887: for ( ; L != []; L = cdr(L) )
888: if ( A == car(L) )
889: return 1;
890: return 0;
891: }
892:
893: /* several functions for computation of normal forms etc. */
894:
895: def p_nf(P,B,V,O) {
896: dp_ord(O); DP = dp_ptod(P,V);
897: N = length(B); DB = newvect(N);
898: for ( I = N-1, IL = []; I >= 0; I-- ) {
899: DB[I] = dp_ptod(B[I],V);
900: IL = cons(I,IL);
901: }
902: return dp_dtop(dp_nf(IL,DP,DB,1),V);
903: }
904:
905: def p_true_nf(P,B,V,O) {
906: dp_ord(O); DP = dp_ptod(P,V);
907: N = length(B); DB = newvect(N);
908: for ( I = N-1, IL = []; I >= 0; I-- ) {
909: DB[I] = dp_ptod(B[I],V);
910: IL = cons(I,IL);
911: }
912: L = dp_true_nf(IL,DP,DB,1);
913: return [dp_dtop(L[0],V),L[1]];
914: }
915:
916: def p_terms(D,V,O)
917: {
918: dp_ord(O);
919: for ( L = [], T = dp_ptod(D,V); T; T = dp_rest(T) )
920: L = cons(dp_dtop(dp_ht(T),V),L);
921: return reverse(L);
922: }
923:
924: def dp_terms(D,V)
925: {
926: for ( L = [], T = D; T; T = dp_rest(T) )
927: L = cons(dp_dtop(dp_ht(T),V),L);
928: return reverse(L);
929: }
930:
931: def gb_comp(A,B)
932: {
933: for ( T = A; T != []; T = cdr(T) ) {
934: for ( S = B, M = car(T), N = -M; S != []; S = cdr(S) )
935: if ( car(S) == M || car(S) == N )
936: break;
937: if ( S == [] )
938: break;
939: }
940: return T == [] ? 1 : 0;
941: }
942:
943: def zero_dim(G,V,O) {
944: dp_ord(O);
945: HL = map(dp_dtop,map(dp_ht,map(dp_ptod,G,V)),V);
946: for ( L = []; HL != []; HL = cdr(HL) )
947: if ( length(vars(car(HL))) == 1 )
948: L = cons(car(HL),L);
949: return length(vars(L)) == length(V) ? 1 : 0;
950: }
951:
952: def hmlist(G,V,O) {
953: dp_ord(O);
954: return map(dp_dtop,map(dp_hm,map(dp_ptod,G,V)),V);
955: }
956:
957: def valid_modulus(HL,M) {
958: V = vars(HL);
959: for ( T = HL; T != []; T = cdr(T) )
960: if ( !dp_mod(dp_ptod(car(T),V),M,[]) )
961: break;
962: return T == [] ? 1 : 0;
963: }
964:
965: def npos_check(DL) {
966: N = size(car(DL))[0];
967: if ( length(DL) != N )
968: return [-1,0];
969: D = newvect(N);
970: for ( I = 0; I < N; I++ ) {
971: for ( J = 0; J < N; J++ )
972: D[J] = 0;
973: D[I] = 1;
974: for ( T = DL; T != []; T = cdr(T) )
975: if ( D == car(T) )
976: break;
977: if ( T != [] )
978: DL = setminus(DL,[car(T)]);
979: }
980: if ( length(DL) != 1 )
981: return [-1,0];
982: U = car(DL);
983: for ( I = 0, J = 0, I0 = -1; I < N; I++ )
984: if ( U[I] ) {
985: I0 = I; J++;
986: }
987: if ( J != 1 )
988: return [-1,0];
989: else
990: return [I0,U[I0]];
991: }
992:
993: def mult_mat(L,TAB,MB)
994: {
995: A = L[0]; DN0 = L[1];
996: for ( NM = 0, DN = 1, I = 0; A; A = dp_rest(A) ) {
997: H = dp_ht(A);
998: for ( ; MB[I] != H; I++ );
999: NM1 = TAB[I][0]; DN1 = TAB[I][1]; I++;
1000: GCD = igcd(DN,DN1); C = DN1/GCD; C1 = DN/GCD;
1001: NM = C*NM + C1*dp_hc(A)*NM1;
1002: DN *= C;
1003: }
1004: Z=remove_cont([NM,DN*DN0]);
1005: return Z;
1006: }
1007:
1008: def sepm(MAT)
1009: {
1010: S = size(MAT); N = S[0]; M = S[1]-1;
1011: A = newmat(N,M); B = newvect(N);
1012: for ( I = 0; I < N; I++ )
1013: for ( J = 0, T1 = MAT[I], T2 = A[I]; J < M; J++ )
1014: T2[J] = T1[J];
1015: for ( I = 0; I < N; I++ )
1016: B[I] = MAT[I][M];
1017: return [A,B];
1018: }
1019:
1020: def henleq(M,MOD)
1021: {
1022: SIZE = size(M); ROW = SIZE[0]; COL = SIZE[1];
1023: W = newvect(COL);
1024: L = sepm(M); A = L[0]; B = L[1];
1025: COUNT = INIT_COUNT?INIT_COUNT:idiv(max_mag(M),54);
1026: if ( !COUNT )
1027: COUNT = 1;
1028:
1029: TINV = TC = TR = TS = TM = TDIV = 0;
1030:
1031: T0 = time()[0];
1032: L = geninvm_swap(A,MOD); INV = L[0]; INDEX = L[1];
1033: TS += time()[0] - T0;
1034:
1035: COL1 = COL - 1;
1036: AA = newmat(COL1,COL1); BB = newvect(COL1);
1037: for ( I = 0; I < COL1; I++ ) {
1038: for ( J = 0, T = AA[I], S = A[INDEX[I]]; J < COL1; J++ )
1039: T[J] = S[J];
1040: BB[I] = B[INDEX[I]];
1041: }
1042: if ( COL1 != ROW ) {
1043: RESTA = newmat(ROW-COL1,COL1); RESTB = newvect(ROW-COL1);
1044: for ( ; I < ROW; I++ ) {
1045: for ( J = 0, T = RESTA[I-COL1], S = A[INDEX[I]]; J < COL1; J++ )
1046: T[J] = S[J];
1047: RESTB[I-COL1] = B[INDEX[I]];
1048: }
1049: } else
1050: RESTA = RESTB = 0;
1051:
1052: MOD2 = idiv(MOD,2);
1053: for ( I = 0, C = BB, X = 0, PK = 1, CCC = 0, ITOR_FAIL = -1; ;
1054: I++, PK *= MOD ) {
1055: if ( COUNT == CCC ) {
1056: CCC = 0;
1057: T0 = time()[0];
1058: ND = intvtoratv(X,PK,ishift(calcb_special(PK,MOD,I),32));
1059: TR += time()[0]-T0;
1060: if ( ND ) {
1061: T0 = time()[0];
1062: F = ND[0]; LCM = ND[1]; T = AA*F+LCM*BB;
1063: TM += time()[0]-T0;
1064: if ( zerovector(T) ) {
1065: T0 = time()[0]; T = RESTA*F+LCM*RESTB; TM += time()[0]-T0;
1066: if ( zerovector(T) ) {
1067: #if 0
1068: if ( dp_gr_print() ) print(["init",TS,"pinv",TINV,"c",TC,"div",TDIV,"rat",TR,"mul",TM]);
1069: #endif
1070: if ( dp_gr_print() ) print("end",2);
1071: return [F,LCM];
1072: } else
1073: return 0;
1074: }
1075: } else {
1076: #if 0
1077: if ( dp_gr_print() ) print(I);
1078: #endif
1079: }
1080: } else {
1081: #if 0
1082: if ( dp_gr_print() ) print([I,TINV,TC,TDIV]);
1083: #endif
1084: if ( dp_gr_print() ) print(".",2);
1085: CCC++;
1086: }
1087: T0 = time()[0];
1088: XT = sremainder(INV*sremainder(-C,MOD),MOD);
1089: XT = map(adj_sgn,XT,MOD,MOD2);
1090: TINV += time()[0] - T0;
1091: X += XT*PK;
1092: T0 = time()[0];
1093: C += mul_mat_vect_int(AA,XT);
1094: TC += time()[0] - T0;
1095: T0 = time()[0]; C = map(idiv,C,MOD); TDIV += time()[0] - T0;
1096: }
1097: }
1098:
1099: def henleq_prep(A,MOD)
1100: {
1101: SIZE = size(A); ROW = SIZE[0]; COL = SIZE[1];
1102: L = geninvm_swap(A,MOD); INV = L[0]; INDEX = L[1];
1103: AA = newmat(COL,COL);
1104: for ( I = 0; I < COL; I++ )
1105: for ( J = 0, T = AA[I], S = A[INDEX[I]]; J < COL; J++ )
1106: T[J] = S[J];
1107: if ( COL != ROW ) {
1108: RESTA = newmat(ROW-COL,COL);
1109: for ( ; I < ROW; I++ )
1110: for ( J = 0, T = RESTA[I-COL], S = A[INDEX[I]]; J < COL; J++ )
1111: T[J] = S[J];
1112: } else
1113: RESTA = 0;
1114: return [[A,AA,RESTA],L];
1115: }
1116:
1117: def henleq_gsl(L,B,MOD)
1118: {
1119: AL = L[0]; INVL = L[1];
1120: A = AL[0]; AA = AL[1]; RESTA = AL[2];
1121: INV = INVL[0]; INDEX = INVL[1];
1122: SIZE = size(A); ROW = SIZE[0]; COL = SIZE[1];
1123: BB = newvect(COL);
1124: for ( I = 0; I < COL; I++ )
1125: BB[I] = B[INDEX[I]];
1126: if ( COL != ROW ) {
1127: RESTB = newvect(ROW-COL);
1128: for ( ; I < ROW; I++ )
1129: RESTB[I-COL] = B[INDEX[I]];
1130: } else
1131: RESTB = 0;
1132:
1133: COUNT = INIT_COUNT?INIT_COUNT:idiv(MAX(max_mag(A),max_mag_vect(B)),54);
1134: if ( !COUNT )
1135: COUNT = 1;
1136: MOD2 = idiv(MOD,2);
1.3 noro 1137: X = newvect(size(AA)[0]);
1138: for ( I = 0, C = BB, PK = 1, CCC = 0, ITOR_FAIL = -1; ;
1.1 noro 1139: I++, PK *= MOD ) {
1140: if ( zerovector(C) )
1141: if ( zerovector(RESTA*X+RESTB) ) {
1142: if ( dp_gr_print() ) print("end",0);
1143: return [X,1];
1144: } else
1145: return 0;
1146: else if ( COUNT == CCC ) {
1147: CCC = 0;
1148: ND = intvtoratv(X,PK,ishift(calcb_special(PK,MOD,I),32));
1149: if ( ND ) {
1150: F = ND[0]; LCM = ND[1]; T = AA*F+LCM*BB;
1151: if ( zerovector(T) ) {
1152: T = RESTA*F+LCM*RESTB;
1153: if ( zerovector(T) ) {
1154: if ( dp_gr_print() ) print("end",0);
1155: return [F,LCM];
1156: } else
1157: return 0;
1158: }
1159: } else {
1160: }
1161: } else {
1162: if ( dp_gr_print() ) print(".",2);
1163: CCC++;
1164: }
1165: XT = sremainder(INV*sremainder(-C,MOD),MOD);
1166: XT = map(adj_sgn,XT,MOD,MOD2);
1167: X += XT*PK;
1168: C += mul_mat_vect_int(AA,XT);
1169: C = map(idiv,C,MOD);
1170: }
1171: }
1172:
1173: def adj_sgn(A,M,M2)
1174: {
1175: return A > M2 ? A-M : A;
1176: }
1177:
1178: def zerovector(C)
1179: {
1180: if ( !C )
1181: return 1;
1182: for ( I = size(C)[0]-1; I >= 0 && !C[I]; I-- );
1183: if ( I < 0 )
1184: return 1;
1185: else
1186: return 0;
1187: }
1188:
1189: def solvem(INV,COMP,V,MOD)
1190: {
1191: T = COMP*V;
1192: N = size(T)[0];
1193: for ( I = 0; I < N; I++ )
1194: if ( T[I] % MOD )
1195: return 0;
1196: return modvect(INV*V,MOD);
1197: }
1198:
1199: def modmat(A,MOD)
1200: {
1201: if ( !A )
1202: return 0;
1203: S = size(A); N = S[0]; M = S[1];
1204: MAT = newmat(N,M);
1205: for ( I = 0, NZ = 0; I < N; I++ )
1206: for ( J = 0, T1 = A[I], T2 = MAT[I]; J < M; J++ ) {
1207: T2[J] = T1[J] % MOD;
1208: NZ = NZ || T2[J];
1209: }
1210: return NZ?MAT:0;
1211: }
1212:
1213: def modvect(A,MOD)
1214: {
1215: if ( !A )
1216: return 0;
1217: N = size(A)[0];
1218: VECT = newvect(N);
1219: for ( I = 0, NZ = 0; I < N; I++ ) {
1220: VECT[I] = A[I] % MOD;
1221: NZ = NZ || VECT[I];
1222: }
1223: return NZ?VECT:0;
1224: }
1225:
1226: def qrmat(A,MOD)
1227: {
1228: if ( !A )
1229: return [0,0];
1230: S = size(A); N = S[0]; M = S[1];
1231: Q = newmat(N,M); R = newmat(N,M);
1232: for ( I = 0, NZQ = 0, NZR = 0; I < N; I++ )
1233: for ( J = 0, TA = A[I], TQ = Q[I], TR = R[I]; J < M; J++ ) {
1234: L = iqr(TA[J],MOD); TQ[J] = L[0]; TR[J] = L[1];
1235: NZQ = NZQ || TQ[J]; NZR = NZR || TR[J];
1236: }
1237: return [NZQ?Q:0,NZR?R:0];
1238: }
1239:
1240: def qrvect(A,MOD)
1241: {
1242: if ( !A )
1243: return [0,0];
1244: N = size(A)[0];
1245: Q = newvect(N); R = newvect(N);
1246: for ( I = 0, NZQ = 0, NZR = 0; I < N; I++ ) {
1247: L = iqr(A[I],MOD); Q[I] = L[0]; R[I] = L[1];
1248: NZQ = NZQ || Q[I]; NZR = NZR || R[I];
1249: }
1250: return [NZQ?Q:0,NZR?R:0];
1251: }
1252:
1253: def max_mag(M)
1254: {
1255: R = size(M)[0];
1256: U = 1;
1257: for ( I = 0; I < R; I++ ) {
1258: A = max_mag_vect(M[I]);
1259: U = MAX(A,U);
1260: }
1261: return U;
1262: }
1263:
1264: def max_mag_vect(V)
1265: {
1266: R = size(V)[0];
1267: U = 1;
1268: for ( I = 0; I < R; I++ ) {
1269: A = dp_mag(V[I]*<<0>>);
1270: U = MAX(A,U);
1271: }
1272: return U;
1273: }
1274:
1275: def gsl_check(B,V,S)
1276: {
1277: N = length(V);
1278: U = S[N-1]; M = U[1]; D = U[2];
1279: W = setminus(V,[var(M)]);
1280: H = uc(); VH = append(W,[H]);
1281: for ( T = B; T != []; T = cdr(T) ) {
1282: A = car(T);
1283: AH = dp_dtop(dp_homo(dp_ptod(A,W)),VH);
1284: for ( I = 0, Z = S; I < N-1; I++, Z = cdr(Z) ) {
1285: L = car(Z); AH = ptozp(subst(AH,L[0],L[1]/L[2]));
1286: }
1287: AH = ptozp(subst(AH,H,D));
1288: R = srem(AH,M);
1289: if ( dp_gr_print() )
1290: if ( !R )
1291: print([A,"ok"]);
1292: else
1293: print([A,"bad"]);
1294: if ( R )
1295: break;
1296: }
1297: return T == [] ? 1 : 0;
1298: }
1299:
1300: def vs_dim(G,V,O)
1301: {
1302: HM = hmlist(G,V,O); ZD = zero_dim(HM,V,O);
1303: if ( ZD ) {
1304: MB = dp_mbase(map(dp_ptod,HM,V));
1305: return length(MB);
1306: } else
1307: error("vs_dim : ideal is not zero-dimensional!");
1308: }
1309:
1.2 noro 1310: def dgr(G,V,O)
1.1 noro 1311: {
1.2 noro 1312: P = getopt(proc);
1313: if ( type(P) == -1 )
1314: return gr(G,V,O);
1.1 noro 1315: P0 = P[0]; P1 = P[1]; P = [P0,P1];
1.2 noro 1316: map(ox_reset,P);
1317: ox_cmo_rpc(P0,"dp_gr_main",G,V,0,1,O);
1318: ox_cmo_rpc(P1,"dp_gr_main",G,V,1,1,O);
1319: map(ox_push_cmd,P,262); /* 262 = OX_popCMO */
1320: F = ox_select(P);
1321: R = ox_get(F[0]);
1322: if ( F[0] == P0 ) {
1323: Win = "nonhomo";
1324: Lose = P1;
1325: } else {
1326: Win = "nhomo";
1327: Lose = P0;
1328: }
1329: ox_reset(Lose);
1330: return [Win,R];
1.1 noro 1331: }
1332:
1333: /* functions for rpc */
1334:
1335: def register_matrix(M)
1336: {
1337: REMOTE_MATRIX = M; return 0;
1338: }
1339:
1340: def register_nfv(L)
1341: {
1342: REMOTE_NF = L[0]; REMOTE_VARS = L[1]; return 0;
1343: }
1344:
1345: def r_ttob(T,M)
1346: {
1347: return hen_ttob(T,REMOTE_NF,0,REMOTE_VARS,M);
1348: }
1349:
1350: def r_ttob_gsl(L,M)
1351: {
1352: return cons(L[2],hen_ttob(L[0],REMOTE_NF,L[1],REMOTE_VARS,M));
1353: }
1354:
1355: def get_matrix()
1356: {
1357: REMOTE_MATRIX;
1.4 noro 1358: }
1359:
1360: extern NFArray$
1361:
1362: /*
1363: * HL = [[c,i,m,d],...]
1364: * if c != 0
1365: * g = 0
1366: * g = (c*g + m*gi)/d
1367: * ...
1368: * finally compare g with NF
1369: * if g == NF then NFArray[NFIndex] = g
1370: *
1371: * if c = 0 then HL consists of single history [0,i,0,0],
1372: * which means that dehomogenization of NFArray[i] should be
1373: * eqall to NF.
1374: */
1375:
1376: def check_trace(NF,NFIndex,HL)
1377: {
1378: if ( !car(HL)[0] ) {
1379: /* dehomogenization */
1380: DH = dp_dehomo(NFArray[car(HL)[1]]);
1381: if ( NF == DH ) {
1382: realloc_NFArray(NFIndex);
1383: NFArray[NFIndex] = NF;
1384: return 0;
1385: } else
1386: error("check_trace(dehomo)");
1387: }
1388:
1389: for ( G = 0, T = HL; T != []; T = cdr(T) ) {
1390: H = car(T);
1391:
1392: Coeff = H[0];
1393: Index = H[1];
1394: Monomial = H[2];
1395: Denominator = H[3];
1396:
1397: Reducer = NFArray[Index];
1398: G = (Coeff*G+Monomial*Reducer)/Denominator;
1399: }
1400: if ( NF == G ) {
1401: realloc_NFArray(NFIndex);
1402: NFArray[NFIndex] = NF;
1403: return 0;
1404: } else
1405: error("check_trace");
1406: }
1407:
1408: /*
1409: * realloc NFArray so that it can hold * an element as NFArray[Ind].
1410: */
1411:
1412: def realloc_NFArray(Ind)
1413: {
1414: if ( Ind == size(NFArray)[0] ) {
1415: New = newvect(Ind + 100);
1416: for ( I = 0; I < Ind; I++ )
1417: New[I] = NFArray[I];
1418: NFArray = New;
1419: }
1420: }
1421:
1422: /*
1423: * create NFArray and initialize it by List.
1424: */
1425:
1426: def register_input(List)
1427: {
1428: Len = length(List);
1429: NFArray = newvect(Len+100,List);
1.1 noro 1430: }
1431: end$