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