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