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