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