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