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