Annotation of OpenXM_contrib2/asir2018/engine/dist.c, Revision 1.27
1.1 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
26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
27: * for such modification or the source code of the modified part of the
28: * SOFTWARE.
29: *
30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
47: *
1.27 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2018/engine/dist.c,v 1.26 2021/12/23 08:55:33 noro Exp $
1.1 noro 49: */
50: #include "ca.h"
51:
52: #define ORD_REVGRADLEX 0
53: #define ORD_GRADLEX 1
54: #define ORD_LEX 2
55: #define ORD_BREVGRADLEX 3
56: #define ORD_BGRADLEX 4
57: #define ORD_BLEX 5
58: #define ORD_BREVREV 6
59: #define ORD_BGRADREV 7
60: #define ORD_BLEXREV 8
61: #define ORD_ELIM 9
62: #define ORD_WEYL_ELIM 10
63: #define ORD_HOMO_WW_DRL 11
64: #define ORD_DRL_ZIGZAG 12
65: #define ORD_HOMO_WW_DRL_ZIGZAG 13
66:
67: int cmpdl_drl_zigzag(), cmpdl_homo_ww_drl_zigzag();
68: int cmpdl_top_weight();
69:
70: int (*cmpdl)()=cmpdl_revgradlex;
71: int (*cmpdl_tie_breaker)();
72: int (*primitive_cmpdl[3])() = {cmpdl_revgradlex,cmpdl_gradlex,cmpdl_lex};
73:
74: Obj current_top_weight;
75: int current_top_weight_len;
76:
77: int do_weyl;
78:
79: int dp_nelim,dp_fcoeffs;
80: struct order_spec *dp_current_spec;
81: struct modorder_spec *dp_current_modspec;
82: int *dp_dl_work;
83:
84: void comm_muld_trunc(VL vl,DP p1,DP p2,DL dl,DP *pr);
85: void comm_quod(VL vl,DP p1,DP p2,DP *pr);
86: void muldm_trunc(VL vl,DP p,MP m0,DL dl,DP *pr);
87: void muldc_trunc(VL vl,DP p,Obj c,DL dl,DP *pr);
88: int create_order_spec(VL vl,Obj obj,struct order_spec **specp);
89: void create_modorder_spec(int id,LIST shift,struct modorder_spec **s);
90:
91: void order_init()
92: {
93: struct order_spec *spec;
94:
95: create_order_spec(0,0,&spec);
96: initd(spec);
97: create_modorder_spec(0,0,&dp_current_modspec);
98: }
99:
100: int has_sfcoef_p(Obj f);
101:
102: int has_sfcoef(DP f)
103: {
104: MP t;
105:
106: if ( !f )
107: return 0;
108: for ( t = BDY(f); t; t = NEXT(t) )
109: if ( has_sfcoef_p(t->c) )
110: break;
111: return t ? 1 : 0;
112: }
113:
114: int has_sfcoef_p(Obj f)
115: {
116: DCP dc;
117:
118: if ( !f )
119: return 0;
120: else if ( NUM(f) )
121: return (NID((Num)f) == N_GFS) ? 1 : 0;
122: else if ( POLY(f) ) {
123: for ( dc = DC((P)f); dc; dc = NEXT(dc) )
124: if ( has_sfcoef_p((Obj)COEF(dc)) )
125: return 1;
126: return 0;
127: } else
128: return 0;
129: }
130:
131: extern Obj nd_top_weight;
132:
133: void reset_top_weight()
134: {
135: cmpdl = cmpdl_tie_breaker;
136: cmpdl_tie_breaker = 0;
137: nd_top_weight = 0;
138: current_top_weight = 0;
139: current_top_weight_len = 0;
140: }
141:
142: void initd(struct order_spec *spec)
143: {
144: int len,i,k,row;
145: Q **mat;
146:
147: switch ( spec->id ) {
148: case 3:
149: cmpdl = cmpdl_composite;
150: dp_dl_work = (int *)MALLOC_ATOMIC(spec->nv*sizeof(int));
151: break;
152: case 2:
153: cmpdl = cmpdl_matrix;
154: dp_dl_work = (int *)MALLOC_ATOMIC(spec->nv*sizeof(int));
155: break;
156: case 1:
157: cmpdl = cmpdl_order_pair;
158: break;
159: default:
160: switch ( spec->ord.simple ) {
161: case ORD_REVGRADLEX:
162: cmpdl = cmpdl_revgradlex; break;
163: case ORD_GRADLEX:
164: cmpdl = cmpdl_gradlex; break;
165: case ORD_BREVGRADLEX:
166: cmpdl = cmpdl_brevgradlex; break;
167: case ORD_BGRADLEX:
168: cmpdl = cmpdl_bgradlex; break;
169: case ORD_BLEX:
170: cmpdl = cmpdl_blex; break;
171: case ORD_BREVREV:
172: cmpdl = cmpdl_brevrev; break;
173: case ORD_BGRADREV:
174: cmpdl = cmpdl_bgradrev; break;
175: case ORD_BLEXREV:
176: cmpdl = cmpdl_blexrev; break;
177: case ORD_ELIM:
178: cmpdl = cmpdl_elim; break;
179: case ORD_WEYL_ELIM:
180: cmpdl = cmpdl_weyl_elim; break;
181: case ORD_HOMO_WW_DRL:
182: cmpdl = cmpdl_homo_ww_drl; break;
183: case ORD_DRL_ZIGZAG:
184: cmpdl = cmpdl_drl_zigzag; break;
185: case ORD_HOMO_WW_DRL_ZIGZAG:
186: cmpdl = cmpdl_homo_ww_drl_zigzag; break;
187: case ORD_LEX: default:
188: cmpdl = cmpdl_lex; break;
189: }
190: break;
191: }
192: if ( current_top_weight ) {
193: cmpdl_tie_breaker = cmpdl;
194: cmpdl = cmpdl_top_weight;
195: if ( OID(current_top_weight) == O_VECT ) {
196: mat = (Q **)&BDY((VECT)current_top_weight);
197: row = 1;
198: } else {
199: mat = (Q **)BDY((MAT)current_top_weight);
200: row = ((MAT)current_top_weight)->row;
201: }
202: for ( k = 0, len = 0; k < row; k++ )
203: for ( i = 0; i < spec->nv; i++ )
204: if ( mat[k][i] )
205: len = MAX(mpz_size(BDY((Z)mat[k][i])),len);
206: current_top_weight_len = len;
207: }
208: dp_current_spec = spec;
209: }
210:
1.3 noro 211: int dpm_ordtype;
1.1 noro 212:
1.27 ! noro 213: #if 1
1.1 noro 214: void ptod(VL vl,VL dvl,P p,DP *pr)
215: {
216: int n,i,j,k;
217: VL tvl;
218: V v;
219: DL d;
220: MP m;
221: DCP dc;
222: DCP *w;
223: DP r,s,t,u;
224: P x,c;
225:
226: if ( !p )
227: *pr = 0;
228: else if ( OID(p) > O_P )
229: error("ptod : only polynomials can be converted.");
230: else {
231: for ( n = 0, tvl = dvl; tvl; tvl = NEXT(tvl), n++ );
232: if ( NUM(p) ) {
233: NEWDL(d,n);
234: NEWMP(m); m->dl = d; C(m) = (Obj)p; NEXT(m) = 0; MKDP(n,m,*pr); (*pr)->sugar = 0;
235: } else {
236: for ( i = 0, tvl = dvl, v = VR(p);
237: tvl && tvl->v != v; tvl = NEXT(tvl), i++ );
238: if ( !tvl ) {
239: for ( dc = DC(p), k = 0; dc; dc = NEXT(dc), k++ );
240: w = (DCP *)ALLOCA(k*sizeof(DCP));
241: for ( dc = DC(p), j = 0; j < k; dc = NEXT(dc), j++ )
242: w[j] = dc;
243:
244: for ( j = k-1, s = 0, MKV(v,x); j >= 0; j-- ) {
245: ptod(vl,dvl,COEF(w[j]),&t); pwrp(vl,x,DEG(w[j]),&c);
246: muldc(vl,t,(Obj)c,&r); addd(vl,r,s,&t); s = t;
247: }
248: *pr = s;
249: } else {
250: for ( dc = DC(p), k = 0; dc; dc = NEXT(dc), k++ );
251: w = (DCP *)ALLOCA(k*sizeof(DCP));
252: for ( dc = DC(p), j = 0; j < k; dc = NEXT(dc), j++ )
253: w[j] = dc;
254:
255: for ( j = k-1, s = 0; j >= 0; j-- ) {
256: ptod(vl,dvl,COEF(w[j]),&t);
1.2 noro 257: NEWDL(d,n); d->d[i] = ZTOS(DEG(w[j]));
1.1 noro 258: d->td = MUL_WEIGHT(d->d[i],i);
259: NEWMP(m); m->dl = d; C(m) = (Obj)ONE; NEXT(m) = 0; MKDP(n,m,u); u->sugar = d->td;
260: comm_muld(vl,t,u,&r); addd(vl,r,s,&t); s = t;
261: }
262: *pr = s;
263: }
264: }
265: }
266: #if 0
267: if ( !dp_fcoeffs && has_sfcoef(*pr) )
268: dp_fcoeffs = N_GFS;
269: #endif
270: }
1.26 noro 271: #else
272: void ptod(VL vl,VL dvl,P p,DP *pr)
273: {
274: int n,i,j,k;
275: VL tvl;
276: V v;
277: DL d;
278: MP m;
279: DCP dc;
280: DP *y;
281: DP r,s,t,u;
282: P x,c;
283:
284: if ( !p )
285: *pr = 0;
286: else if ( OID(p) > O_P )
287: error("ptod : only polynomials can be converted.");
288: else {
289: for ( n = 0, tvl = dvl; tvl; tvl = NEXT(tvl), n++ );
290: if ( NUM(p) ) {
291: NEWDL(d,n);
292: NEWMP(m); m->dl = d; C(m) = (Obj)p; NEXT(m) = 0; MKDP(n,m,*pr); (*pr)->sugar = 0;
293: } else {
294: for ( i = 0, tvl = dvl, v = VR(p); tvl && tvl->v != v; tvl = NEXT(tvl), i++ );
295: for ( dc = DC(p), k = 0; dc; dc = NEXT(dc), k++ );
296: y = (DP *)ALLOCA(k*sizeof(DP));
297: if ( !tvl ) {
298: MKV(v,x);
299: for ( dc = DC(p), j = 0; dc; dc = NEXT(dc), j++ ) {
300: ptod(vl,dvl,COEF(dc),&t); pwrp(vl,x,DEG(dc),&c);
301: muldc(vl,t,(Obj)c,&y[j]);
302: }
303: } else {
304: for ( dc = DC(p), j = 0; dc; dc = NEXT(dc), j++ ) {
305: ptod(vl,dvl,COEF(dc),&t);
306: NEWDL(d,n);
307: d->d[i] = ZTOS(DEG(dc));
308: d->td = MUL_WEIGHT(d->d[i],i);
309: NEWMP(m); m->dl = d; C(m) = (Obj)ONE; NEXT(m) = 0; MKDP(n,m,u); u->sugar = d->td;
310: comm_muld(vl,t,u,&y[j]);
311: }
312: }
313: for ( j = k-1, s = 0; j >= 0; j-- ) {
314: addd(vl,y[j],s,&t); s = t;
315: }
316: *pr = s;
317: }
318: }
319: #if 0
320: if ( !dp_fcoeffs && has_sfcoef(*pr) )
321: dp_fcoeffs = N_GFS;
322: #endif
323: }
324: #endif
1.1 noro 325:
326: void dtop(VL vl,VL dvl,DP p,Obj *pr)
327: {
328: int n,i,j,k;
329: DL d;
330: MP m;
331: MP *a;
332: P r;
333: Obj t,w,s,u;
334: Z q;
335: VL tvl;
336:
337: if ( !p )
338: *pr = 0;
339: else {
340: for ( k = 0, m = BDY(p); m; m = NEXT(m), k++ );
341: a = (MP *)ALLOCA(k*sizeof(MP));
342: for ( j = 0, m = BDY(p); j < k; m = NEXT(m), j++ )
343: a[j] = m;
344:
345: for ( n = p->nv, j = k-1, s = 0; j >= 0; j-- ) {
346: m = a[j];
347: t = C(m);
348: if ( NUM(t) && NID((Num)t) == N_M ) {
349: mptop((P)t,(P *)&u); t = u;
350: }
351: for ( i = 0, d = m->dl, tvl = dvl;
352: i < n; tvl = NEXT(tvl), i++ ) {
1.2 noro 353: MKV(tvl->v,r); STOZ(d->d[i],q); pwrp(vl,r,q,(P *)&u);
1.1 noro 354: arf_mul(vl,t,(Obj)u,&w); t = w;
355: }
356: arf_add(vl,s,t,&u); s = u;
357: }
358: *pr = s;
359: }
360: }
361:
362: void nodetod(NODE node,DP *dp)
363: {
364: NODE t;
365: int len,i,td;
366: Q e;
367: DL d;
368: MP m;
369: DP u;
370:
371: for ( t = node, len = 0; t; t = NEXT(t), len++ );
372: NEWDL(d,len);
373: for ( t = node, i = 0, td = 0; i < len; t = NEXT(t), i++ ) {
374: e = (Q)BDY(t);
375: if ( !e )
376: d->d[i] = 0;
377: else if ( !NUM(e) || !RATN(e) || !INT(e) )
378: error("nodetod : invalid input");
379: else {
1.2 noro 380: d->d[i] = ZTOS((Q)e); td += MUL_WEIGHT(d->d[i],i);
1.1 noro 381: }
382: }
383: d->td = td;
384: NEWMP(m); m->dl = d; C(m) = (Obj)ONE; NEXT(m) = 0;
385: MKDP(len,m,u); u->sugar = td; *dp = u;
386: }
387:
388: void nodetodpm(NODE node,Obj pos,DPM *dp)
389: {
390: NODE t;
1.9 noro 391: int len,i,td,p;
1.1 noro 392: Q e;
393: DL d;
394: DMM m;
395: DPM u;
396:
397: for ( t = node, len = 0; t; t = NEXT(t), len++ );
398: NEWDL(d,len);
399: for ( t = node, i = 0, td = 0; i < len; t = NEXT(t), i++ ) {
400: e = (Q)BDY(t);
401: if ( !e )
402: d->d[i] = 0;
403: else if ( !NUM(e) || !RATN(e) || !INT(e) )
404: error("nodetodpm : invalid input");
405: else {
1.2 noro 406: d->d[i] = ZTOS((Q)e); td += MUL_WEIGHT(d->d[i],i);
1.1 noro 407: }
408: }
409: d->td = td;
1.9 noro 410: p = ZTOS((Q)pos);
1.19 noro 411: if ( dp_current_spec->module_top_weight ) {
1.9 noro 412: if ( p > dp_current_spec->module_rank )
413: error("nodetodpm : inconsistent order spec");
414: d->td += dp_current_spec->module_top_weight[p-1];
415: }
416: NEWDMM(m); m->dl = d; m->pos = p; C(m) = (Obj)ONE; NEXT(m) = 0;
1.1 noro 417: MKDPM(len,m,u); u->sugar = td; *dp = u;
418: }
419:
420: void dtodpm(DP d,int pos,DPM *dp)
421: {
422: DMM mr0,mr;
423: MP m;
1.9 noro 424: int shift;
1.1 noro 425:
426: if ( !d ) *dp = 0;
427: else {
1.9 noro 428: shift = 0;
1.19 noro 429: if ( dp_current_spec->module_top_weight ) {
1.9 noro 430: if ( pos > dp_current_spec->module_rank )
431: error("nodetodpm : inconsistent order spec");
432: shift = dp_current_spec->module_top_weight[pos-1];
433: }
1.1 noro 434: for ( m = BDY(d), mr0 = 0; m; m = NEXT(m) ) {
435: NEXTDMM(mr0,mr);
436: mr->dl = m->dl;
1.9 noro 437: mr->dl->td += shift;
1.1 noro 438: mr->pos = pos;
439: C(mr) = C(m);
440: }
441: MKDPM(d->nv,mr0,*dp); (*dp)->sugar = d->sugar;
442: }
443: }
444:
445: int sugard(MP m)
446: {
447: int s;
448:
449: for ( s = 0; m; m = NEXT(m) )
450: s = MAX(s,m->dl->td);
451: return s;
452: }
453:
454: void addd(VL vl,DP p1,DP p2,DP *pr)
455: {
456: int n;
457: MP m1,m2,mr=0,mr0;
458: Obj t;
459: DL d;
460:
461: if ( !p1 )
462: *pr = p2;
463: else if ( !p2 )
464: *pr = p1;
465: else {
466: if ( OID(p1) <= O_R ) {
467: n = NV(p2); NEWDL(d,n);
468: NEWMP(m1); m1->dl = d; C(m1) = (Obj)p1; NEXT(m1) = 0;
469: MKDP(n,m1,p1); (p1)->sugar = 0;
470: }
471: if ( OID(p2) <= O_R ) {
472: n = NV(p1); NEWDL(d,n);
473: NEWMP(m2); m2->dl = d; C(m2) = (Obj)p2; NEXT(m2) = 0;
474: MKDP(n,m2,p2); (p2)->sugar = 0;
475: }
476: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
477: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
478: case 0:
479: arf_add(vl,C(m1),C(m2),&t);
480: if ( t ) {
481: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = t;
482: }
483: m1 = NEXT(m1); m2 = NEXT(m2); break;
484: case 1:
485: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
486: m1 = NEXT(m1); break;
487: case -1:
488: NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
489: m2 = NEXT(m2); break;
490: }
491: if ( !mr0 )
492: if ( m1 )
493: mr0 = m1;
494: else if ( m2 )
495: mr0 = m2;
496: else {
497: *pr = 0;
498: return;
499: }
500: else if ( m1 )
501: NEXT(mr) = m1;
502: else if ( m2 )
503: NEXT(mr) = m2;
504: else
505: NEXT(mr) = 0;
506: MKDP(NV(p1),mr0,*pr);
507: if ( *pr )
508: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
509: }
510: }
511:
512: /* for F4 symbolic reduction */
513:
514: void symb_addd(DP p1,DP p2,DP *pr)
515: {
516: int n;
517: MP m1,m2,mr=0,mr0;
518:
519: if ( !p1 )
520: *pr = p2;
521: else if ( !p2 )
522: *pr = p1;
523: else {
524: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
525: NEXTMP(mr0,mr); C(mr) = (Obj)ONE;
526: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
527: case 0:
528: mr->dl = m1->dl;
529: m1 = NEXT(m1); m2 = NEXT(m2); break;
530: case 1:
531: mr->dl = m1->dl;
532: m1 = NEXT(m1); break;
533: case -1:
534: mr->dl = m2->dl;
535: m2 = NEXT(m2); break;
536: }
537: }
538: if ( !mr0 )
539: if ( m1 )
540: mr0 = m1;
541: else if ( m2 )
542: mr0 = m2;
543: else {
544: *pr = 0;
545: return;
546: }
547: else if ( m1 )
548: NEXT(mr) = m1;
549: else if ( m2 )
550: NEXT(mr) = m2;
551: else
552: NEXT(mr) = 0;
553: MKDP(NV(p1),mr0,*pr);
554: if ( *pr )
555: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
556: }
557: }
558:
559: /*
560: * destructive merge of two list
561: *
562: * p1, p2 : list of DL
563: * return : a merged list
564: */
565:
566: NODE symb_merge(NODE m1,NODE m2,int n)
567: {
568: NODE top=0,prev,cur,m=0,t;
569: DL d1,d2;
570:
571: if ( !m1 )
572: return m2;
573: else if ( !m2 )
574: return m1;
575: else {
576: switch ( (*cmpdl)(n,(DL)BDY(m1),(DL)BDY(m2)) ) {
577: case 0:
578: top = m1; m = NEXT(m2);
579: break;
580: case 1:
581: top = m1; m = m2;
582: break;
583: case -1:
584: top = m2; m = m1;
585: break;
586: }
587: prev = top; cur = NEXT(top);
588: /* BDY(prev) > BDY(m) always holds */
589: while ( cur && m ) {
590: d1 = (DL)BDY(cur);
591: d2 = (DL)BDY(m);
592: #if 1
593: switch ( (*cmpdl)(n,(DL)BDY(cur),(DL)BDY(m)) ) {
594: #else
595: /* XXX only valid for DRL */
596: if ( d1->td > d2->td )
597: c = 1;
598: else if ( d1->td < d2->td )
599: c = -1;
600: else {
601: for ( i = n-1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
602: if ( i < 0 )
603: c = 0;
604: else if ( d1->d[i] < d2->d[i] )
605: c = 1;
606: else
607: c = -1;
608: }
609: switch ( c ) {
610: #endif
611: case 0:
612: m = NEXT(m);
613: prev = cur; cur = NEXT(cur);
614: break;
615: case 1:
616: t = NEXT(cur); NEXT(cur) = m; m = t;
617: prev = cur; cur = NEXT(cur);
618: break;
619: case -1:
620: NEXT(prev) = m; m = cur;
621: prev = NEXT(prev); cur = NEXT(prev);
622: break;
623: }
624: }
625: if ( !cur )
626: NEXT(prev) = m;
627: return top;
628: }
629: }
630:
631: void _adddl(int n,DL d1,DL d2,DL d3)
632: {
633: int i;
634:
635: d3->td = d1->td+d2->td;
636: for ( i = 0; i < n; i++ )
637: d3->d[i] = d1->d[i]+d2->d[i];
638: }
639:
1.25 noro 640: void _subdl(int n,DL d1,DL d2,DL d3)
641: {
642: int i;
643:
644: d3->td = d1->td-d2->td;
645: for ( i = 0; i < n; i++ )
646: d3->d[i] = d1->d[i]-d2->d[i];
647: }
648:
649:
1.3 noro 650: void _addtodl(int n,DL d1,DL d2)
651: {
652: int i;
653:
654: d2->td += d1->td;
655: for ( i = 0; i < n; i++ )
656: d2->d[i] += d1->d[i];
657: }
658:
1.23 noro 659: void _subfromdl(int n,DL d1,DL d2)
660: {
661: int i;
662:
663: d2->td -= d1->td;
664: for ( i = 0; i < n; i++ )
665: d2->d[i] -= d1->d[i];
666: }
667:
1.3 noro 668: void _copydl(int n,DL d1,DL d2)
669: {
670: int i;
671:
672: d2->td = d1->td;
673: for ( i = 0; i < n; i++ )
674: d2->d[i] = d1->d[i];
675: }
676:
677: int _eqdl(int n,DL d1,DL d2)
678: {
679: int i;
680:
681: if ( d2->td != d1->td ) return 0;
682: for ( i = 0; i < n; i++ )
683: if ( d2->d[i] != d1->d[i] ) return 0;
684: return 1;
685: }
686:
1.1 noro 687: /* m1 <- m1 U dl*f, destructive */
688:
689: NODE mul_dllist(DL dl,DP f);
690:
691: NODE symb_mul_merge(NODE m1,DL dl,DP f,int n)
692: {
693: NODE top,prev,cur,n1;
694: DP g;
695: DL t,s;
696: MP m;
697:
698: if ( !m1 )
699: return mul_dllist(dl,f);
700: else if ( !f )
701: return m1;
702: else {
703: m = BDY(f);
704: NEWDL_NOINIT(t,n);
705: _adddl(n,m->dl,dl,t);
706: top = m1; prev = 0; cur = m1;
707: while ( m ) {
708: switch ( (*cmpdl)(n,(DL)BDY(cur),t) ) {
709: case 0:
710: prev = cur; cur = NEXT(cur);
711: if ( !cur ) {
712: MKDP(n,m,g);
713: NEXT(prev) = mul_dllist(dl,g);
714: return top;
715: }
716: m = NEXT(m);
717: if ( m ) _adddl(n,m->dl,dl,t);
718: break;
719: case 1:
720: prev = cur; cur = NEXT(cur);
721: if ( !cur ) {
722: MKDP(n,m,g);
723: NEXT(prev) = mul_dllist(dl,g);
724: return top;
725: }
726: break;
727: case -1:
728: NEWDL_NOINIT(s,n);
729: s->td = t->td;
730: bcopy(t->d,s->d,n*sizeof(int));
731: NEWNODE(n1);
732: n1->body = (pointer)s;
733: NEXT(n1) = cur;
734: if ( !prev ) {
735: top = n1; cur = n1;
736: } else {
737: NEXT(prev) = n1; prev = n1;
738: }
739: m = NEXT(m);
740: if ( m ) _adddl(n,m->dl,dl,t);
741: break;
742: }
743: }
744: return top;
745: }
746: }
747:
748: DLBUCKET symb_merge_bucket(DLBUCKET m1,DLBUCKET m2,int n)
749: {
750: DLBUCKET top,prev,cur,m,t;
751:
752: if ( !m1 )
753: return m2;
754: else if ( !m2 )
755: return m1;
756: else {
757: if ( m1->td == m2->td ) {
758: top = m1;
759: BDY(top) = symb_merge(BDY(top),BDY(m2),n);
760: m = NEXT(m2);
761: } else if ( m1->td > m2->td ) {
762: top = m1; m = m2;
763: } else {
764: top = m2; m = m1;
765: }
766: prev = top; cur = NEXT(top);
767: /* prev->td > m->td always holds */
768: while ( cur && m ) {
769: if ( cur->td == m->td ) {
770: BDY(cur) = symb_merge(BDY(cur),BDY(m),n);
771: m = NEXT(m);
772: prev = cur; cur = NEXT(cur);
773: } else if ( cur->td > m->td ) {
774: t = NEXT(cur); NEXT(cur) = m; m = t;
775: prev = cur; cur = NEXT(cur);
776: } else {
777: NEXT(prev) = m; m = cur;
778: prev = NEXT(prev); cur = NEXT(prev);
779: }
780: }
781: if ( !cur )
782: NEXT(prev) = m;
783: return top;
784: }
785: }
786:
787: void subd(VL vl,DP p1,DP p2,DP *pr)
788: {
789: DP t;
790:
791: if ( !p2 )
792: *pr = p1;
793: else {
794: chsgnd(p2,&t); addd(vl,p1,t,pr);
795: }
796: }
797:
798: void chsgnd(DP p,DP *pr)
799: {
800: MP m,mr=0,mr0;
801: Obj r;
802:
803: if ( !p )
804: *pr = 0;
805: else if ( OID(p) <= O_R ) {
806: arf_chsgn((Obj)p,&r); *pr = (DP)r;
807: } else {
808: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
809: NEXTMP(mr0,mr); arf_chsgn(C(m),&C(mr)); mr->dl = m->dl;
810: }
811: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
812: if ( *pr )
813: (*pr)->sugar = p->sugar;
814: }
815: }
816:
817: void muld(VL vl,DP p1,DP p2,DP *pr)
818: {
819: if ( ! do_weyl )
820: comm_muld(vl,p1,p2,pr);
821: else
822: weyl_muld(vl,p1,p2,pr);
823: }
824:
825: void comm_muld(VL vl,DP p1,DP p2,DP *pr)
826: {
827: MP m;
828: DP s,t,u;
829: int i,l,l1;
830: static MP *w;
831: static int wlen;
832:
833: if ( !p1 || !p2 )
834: *pr = 0;
835: else if ( OID(p1) != O_DP )
836: muldc(vl,p2,(Obj)p1,pr);
837: else if ( OID(p2) != O_DP )
838: muldc(vl,p1,(Obj)p2,pr);
839: else {
840: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
841: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
842: if ( l1 < l ) {
843: t = p1; p1 = p2; p2 = t;
844: l = l1;
845: }
846: if ( l > wlen ) {
847: if ( w ) GCFREE(w);
848: w = (MP *)MALLOC(l*sizeof(MP));
849: wlen = l;
850: }
851: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
852: w[i] = m;
853: for ( s = 0, i = l-1; i >= 0; i-- ) {
854: muldm(vl,p1,w[i],&t); addd(vl,s,t,&u); s = u;
855: }
856: bzero(w,l*sizeof(MP));
857: *pr = s;
858: }
859: }
860:
861: /* discard terms which is not a multiple of dl */
862:
863: void comm_muld_trunc(VL vl,DP p1,DP p2,DL dl,DP *pr)
864: {
865: MP m;
866: DP s,t,u;
867: int i,l,l1;
868: static MP *w;
869: static int wlen;
870:
871: if ( !p1 || !p2 )
872: *pr = 0;
873: else if ( OID(p1) != O_DP )
874: muldc_trunc(vl,p2,(Obj)p1,dl,pr);
875: else if ( OID(p2) != O_DP )
876: muldc_trunc(vl,p1,(Obj)p2,dl,pr);
877: else {
878: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
879: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
880: if ( l1 < l ) {
881: t = p1; p1 = p2; p2 = t;
882: l = l1;
883: }
884: if ( l > wlen ) {
885: if ( w ) GCFREE(w);
886: w = (MP *)MALLOC(l*sizeof(MP));
887: wlen = l;
888: }
889: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
890: w[i] = m;
891: for ( s = 0, i = l-1; i >= 0; i-- ) {
892: muldm_trunc(vl,p1,w[i],dl,&t); addd(vl,s,t,&u); s = u;
893: }
894: bzero(w,l*sizeof(MP));
895: *pr = s;
896: }
897: }
898:
899: void comm_quod(VL vl,DP p1,DP p2,DP *pr)
900: {
901: MP m=0,m0;
902: DP s,t;
903: int i,n,sugar;
904: DL d1,d2,d;
905: Q a,b;
906:
907: if ( !p2 )
908: error("comm_quod : invalid input");
909: if ( !p1 )
910: *pr = 0;
911: else {
912: n = NV(p1);
913: d2 = BDY(p2)->dl;
914: m0 = 0;
915: sugar = p1->sugar;
916: while ( p1 ) {
917: d1 = BDY(p1)->dl;
918: NEWDL(d,n);
919: d->td = d1->td - d2->td;
920: for ( i = 0; i < n; i++ )
921: d->d[i] = d1->d[i]-d2->d[i];
922: NEXTMP(m0,m);
923: m->dl = d;
924: divq((Q)BDY(p1)->c,(Q)BDY(p2)->c,&a); chsgnq(a,&b);
925: C(m) = (Obj)b;
926: muldm_trunc(vl,p2,m,d2,&t);
927: addd(vl,p1,t,&s); p1 = s;
928: C(m) = (Obj)a;
929: }
930: if ( m0 ) {
931: NEXT(m) = 0; MKDP(n,m0,*pr);
932: } else
933: *pr = 0;
934: /* XXX */
935: if ( *pr )
936: (*pr)->sugar = sugar - d2->td;
937: }
938: }
939:
940: void muldm(VL vl,DP p,MP m0,DP *pr)
941: {
942: MP m,mr=0,mr0;
943: Obj c;
944: DL d;
945: int n;
946:
947: if ( !p )
948: *pr = 0;
949: else {
950: for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
951: m; m = NEXT(m) ) {
952: NEXTMP(mr0,mr);
953: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
954: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
955: else
956: arf_mul(vl,C(m),c,&C(mr));
957: adddl(n,m->dl,d,&mr->dl);
958: }
959: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
960: if ( *pr )
961: (*pr)->sugar = p->sugar + m0->dl->td;
962: }
963: }
964:
965: void muldm_trunc(VL vl,DP p,MP m0,DL dl,DP *pr)
966: {
967: MP m,mr=0,mr0;
968: Obj c;
969: DL d,tdl;
970: int n,i;
971:
972: if ( !p )
973: *pr = 0;
974: else {
975: n = NV(p);
976: NEWDL(tdl,n);
977: for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl;
978: m; m = NEXT(m) ) {
979: _adddl(n,m->dl,d,tdl);
980: for ( i = 0; i < n; i++ )
981: if ( tdl->d[i] < dl->d[i] )
982: break;
983: if ( i < n )
984: continue;
985: NEXTMP(mr0,mr);
986: mr->dl = tdl;
987: NEWDL(tdl,n);
988: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
989: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
990: else
991: arf_mul(vl,C(m),(Obj)c,&C(mr));
992: }
993: if ( mr0 ) {
994: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
995: } else
996: *pr = 0;
997: if ( *pr )
998: (*pr)->sugar = p->sugar + m0->dl->td;
999: }
1000: }
1001:
1002: void weyl_muld(VL vl,DP p1,DP p2,DP *pr)
1003: {
1004: MP m;
1005: DP s,t,u;
1006: int i,l;
1007: static MP *w;
1008: static int wlen;
1009:
1010: if ( !p1 || !p2 )
1011: *pr = 0;
1012: else if ( OID(p1) != O_DP )
1013: muldc(vl,p2,(Obj)p1,pr);
1014: else if ( OID(p2) != O_DP )
1015: muldc(vl,p1,(Obj)p2,pr);
1016: else {
1017: for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
1018: if ( l > wlen ) {
1019: if ( w ) GCFREE(w);
1020: w = (MP *)MALLOC(l*sizeof(MP));
1021: wlen = l;
1022: }
1023: for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
1024: w[i] = m;
1025: for ( s = 0, i = l-1; i >= 0; i-- ) {
1026: weyl_muldm(vl,w[i],p2,&t); addd(vl,s,t,&u); s = u;
1027: }
1028: bzero(w,l*sizeof(MP));
1029: *pr = s;
1030: }
1031: }
1032:
1033: void actm(VL vl,int nv,MP m1,MP m2,DP *pr)
1034: {
1035: DL d1,d2,d;
1036: int n2,i,j,k;
1037: Z jq,c,c1;
1038: MP m;
1039: Obj t;
1040:
1041: d1 = m1->dl;
1042: d2 = m2->dl;
1043: for ( i = 0; i < nv; i++ )
1044: if ( d1->d[i] > d2->d[i] ) {
1045: *pr = 0; return;
1046: }
1047: NEWDL(d,nv);
1048: c = ONE;
1049: for ( i = 0; i < nv; i++ ) {
1050: for ( j = d2->d[i], k = d1->d[i]; k > 0; k--, j-- ) {
1.2 noro 1051: STOZ(j,jq); mulz(c,jq,&c1); c = c1;
1.1 noro 1052: }
1053: d->d[i] = d2->d[i]-d1->d[i];
1054: }
1055: arf_mul(vl,C(m1),C(m2),&t);
1056: NEWMP(m);
1057: arf_mul(vl,(Obj)c,t,&C(m));
1058: m->dl = d;
1059: MKDP(nv,m,*pr);
1060: }
1061:
1062: void weyl_actd(VL vl,DP p1,DP p2,DP *pr)
1063: {
1064: int n;
1065: MP m1,m2;
1066: DP d,r,s;
1067:
1068: if ( !p1 || !p2 ) *pr = 0;
1069: else {
1070: n = NV(p1);
1071: r = 0;
1072: for ( m1 = BDY(p1); m1; m1 = NEXT(m1) )
1073: for ( m2 = BDY(p2); m2; m2 = NEXT(m2) ) {
1074: actm(vl,n,m1,m2,&d);
1075: addd(vl,r,d,&s); r = s;
1076: }
1077: *pr = r;
1078: }
1079: }
1080:
1081: /* monomial * polynomial */
1082:
1083: void weyl_muldm(VL vl,MP m0,DP p,DP *pr)
1084: {
1085: DP r,t,t1;
1086: MP m;
1087: DL d0;
1088: int n,n2,l,i,j,tlen;
1089: static MP *w,*psum;
1090: static struct cdl *tab;
1091: static int wlen;
1092: static int rtlen;
1093:
1094: if ( !p )
1095: *pr = 0;
1096: else {
1097: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
1098: if ( l > wlen ) {
1099: if ( w ) GCFREE(w);
1100: w = (MP *)MALLOC(l*sizeof(MP));
1101: wlen = l;
1102: }
1103: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
1104: w[i] = m;
1105:
1106: n = NV(p); n2 = n>>1;
1107: d0 = m0->dl;
1108: for ( i = 0, tlen = 1; i < n2; i++ )
1109: tlen *= d0->d[n2+i]+1;
1110: if ( tlen > rtlen ) {
1111: if ( tab ) GCFREE(tab);
1112: if ( psum ) GCFREE(psum);
1113: rtlen = tlen;
1114: tab = (struct cdl *)MALLOC(rtlen*sizeof(struct cdl));
1115: psum = (MP *)MALLOC(rtlen*sizeof(MP));
1116: }
1117: bzero(psum,tlen*sizeof(MP));
1118: for ( i = l-1; i >= 0; i-- ) {
1119: bzero(tab,tlen*sizeof(struct cdl));
1120: weyl_mulmm(vl,m0,w[i],n,tab,tlen);
1121: for ( j = 0; j < tlen; j++ ) {
1122: if ( tab[j].c ) {
1123: NEWMP(m); m->dl = tab[j].d; C(m) = (Obj)tab[j].c; NEXT(m) = psum[j];
1124: psum[j] = m;
1125: }
1126: }
1127: }
1128: for ( j = tlen-1, r = 0; j >= 0; j-- )
1129: if ( psum[j] ) {
1130: MKDP(n,psum[j],t); addd(vl,r,t,&t1); r = t1;
1131: }
1132: if ( r )
1133: r->sugar = p->sugar + m0->dl->td;
1134: *pr = r;
1135: }
1136: }
1137:
1138: /* m0 = x0^d0*x1^d1*... * dx0^e0*dx1^e1*... */
1139: /* rtab : array of length (e0+1)*(e1+1)*... */
1140:
1141: void weyl_mulmm(VL vl,MP m0,MP m1,int n,struct cdl *rtab,int rtablen)
1142: {
1143: Obj c,c0,c1;
1144: DL d,d0,d1,dt;
1145: int i,j,a,b,k,l,n2,s,min,curlen;
1146: struct cdl *p;
1147: static Z *ctab;
1148: static struct cdl *tab;
1149: static int tablen;
1150: static struct cdl *tmptab;
1151: static int tmptablen;
1152:
1153:
1154: if ( !m0 || !m1 ) {
1155: rtab[0].c = 0;
1156: rtab[0].d = 0;
1157: return;
1158: }
1159: c0 = C(m0); c1 = C(m1);
1160: arf_mul(vl,c0,c1,&c);
1161: d0 = m0->dl; d1 = m1->dl;
1162: n2 = n>>1;
1163: curlen = 1;
1164: NEWDL(d,n);
1165: if ( n & 1 )
1166: /* offset of h-degree */
1167: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
1168: else
1169: d->td = 0;
1170: rtab[0].c = c;
1171: rtab[0].d = d;
1172:
1173: if ( rtablen > tmptablen ) {
1174: if ( tmptab ) GCFREE(tmptab);
1175: tmptab = (struct cdl *)MALLOC(rtablen*sizeof(struct cdl));
1176: tmptablen = rtablen;
1177: }
1178: for ( i = 0; i < n2; i++ ) {
1179: a = d0->d[i]; b = d1->d[n2+i];
1180: k = d0->d[n2+i]; l = d1->d[i];
1181:
1182: /* degree of xi^a*(Di^k*xi^l)*Di^b */
1183: a += l;
1184: b += k;
1185: s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);
1186:
1187: if ( !k || !l ) {
1188: for ( j = 0, p = rtab; j < curlen; j++, p++ ) {
1189: if ( p->c ) {
1190: dt = p->d;
1191: dt->d[i] = a;
1192: dt->d[n2+i] = b;
1193: dt->td += s;
1194: }
1195: }
1196: curlen *= k+1;
1197: continue;
1198: }
1199: if ( k+1 > tablen ) {
1200: if ( tab ) GCFREE(tab);
1201: if ( ctab ) GCFREE(ctab);
1202: tablen = k+1;
1203: tab = (struct cdl *)MALLOC(tablen*sizeof(struct cdl));
1204: ctab = (Z *)MALLOC(tablen*sizeof(Q));
1205: }
1206: /* compute xi^a*(Di^k*xi^l)*Di^b */
1207: min = MIN(k,l);
1208: mkwc(k,l,ctab);
1209: bzero(tab,(k+1)*sizeof(struct cdl));
1210: if ( n & 1 )
1211: for ( j = 0; j <= min; j++ ) {
1212: NEWDL(d,n);
1213: d->d[i] = a-j; d->d[n2+i] = b-j;
1214: d->td = s;
1215: d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
1216: tab[j].d = d;
1217: tab[j].c = (Obj)ctab[j];
1218: }
1219: else
1220: for ( j = 0; j <= min; j++ ) {
1221: NEWDL(d,n);
1222: d->d[i] = a-j; d->d[n2+i] = b-j;
1223: d->td = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i); /* XXX */
1224: tab[j].d = d;
1225: tab[j].c = (Obj)ctab[j];
1226: }
1227: bzero(ctab,(min+1)*sizeof(Q));
1228: comm_muld_tab(vl,n,rtab,curlen,tab,k+1,tmptab);
1229: curlen *= k+1;
1230: bcopy(tmptab,rtab,curlen*sizeof(struct cdl));
1231: }
1232: }
1233:
1234: /* direct product of two cdl tables
1235: rt[] = [
1236: t[0]*t1[0],...,t[n-1]*t1[0],
1237: t[0]*t1[1],...,t[n-1]*t1[1],
1238: ...
1239: t[0]*t1[n1-1],...,t[n-1]*t1[n1-1]
1240: ]
1241: */
1242:
1243: void comm_muld_tab(VL vl,int nv,struct cdl *t,int n,struct cdl *t1,int n1,struct cdl *rt)
1244: {
1245: int i,j;
1246: struct cdl *p;
1247: Obj c;
1248: DL d;
1249:
1250: bzero(rt,n*n1*sizeof(struct cdl));
1251: for ( j = 0, p = rt; j < n1; j++ ) {
1252: c = (Obj)t1[j].c;
1253: d = t1[j].d;
1254: if ( !c )
1255: break;
1256: for ( i = 0; i < n; i++, p++ ) {
1257: if ( t[i].c ) {
1258: arf_mul(vl,(Obj)t[i].c,c,(Obj *)&p->c);
1259: adddl(nv,t[i].d,d,&p->d);
1260: }
1261: }
1262: }
1263: }
1264:
1265: void muldc(VL vl,DP p,Obj c,DP *pr)
1266: {
1267: MP m,mr=0,mr0;
1268:
1269: if ( !p || !c )
1270: *pr = 0;
1271: else if ( NUM(c) && UNIQ((Q)c) )
1272: *pr = p;
1273: else if ( NUM(c) && MUNIQ((Q)c) )
1274: chsgnd(p,pr);
1275: else {
1276: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1277: NEXTMP(mr0,mr);
1278: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
1279: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
1280: else
1281: arf_mul(vl,C(m),c,&C(mr));
1282: mr->dl = m->dl;
1283: }
1284: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
1285: if ( *pr )
1286: (*pr)->sugar = p->sugar;
1287: }
1288: }
1289:
1290: void divdc(VL vl,DP p,Obj c,DP *pr)
1291: {
1292: Obj inv;
1293:
1294: arf_div(vl,(Obj)ONE,c,&inv);
1295: muld(vl,p,(DP)inv,pr);
1296: }
1297:
1298: void muldc_trunc(VL vl,DP p,Obj c,DL dl,DP *pr)
1299: {
1300: MP m,mr=0,mr0;
1301: DL mdl;
1302: int i,n;
1303:
1304: if ( !p || !c ) {
1305: *pr = 0; return;
1306: }
1307: n = NV(p);
1308: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1309: mdl = m->dl;
1310: for ( i = 0; i < n; i++ )
1311: if ( mdl->d[i] < dl->d[i] )
1312: break;
1313: if ( i < n )
1314: break;
1315: NEXTMP(mr0,mr);
1316: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
1317: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
1318: else
1319: arf_mul(vl,C(m),c,&C(mr));
1320: mr->dl = m->dl;
1321: }
1322: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
1323: if ( *pr )
1324: (*pr)->sugar = p->sugar;
1325: }
1326:
1327: void divsdc(VL vl,DP p,P c,DP *pr)
1328: {
1329: MP m,mr=0,mr0;
1330:
1331: if ( !c )
1332: error("disvsdc : division by 0");
1333: else if ( !p )
1334: *pr = 0;
1335: else if ( OID(c) > O_P )
1336: error("divsdc : invalid argument");
1337: else {
1338: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1339: NEXTMP(mr0,mr); divsp(vl,(P)C(m),c,(P *)&C(mr)); mr->dl = m->dl;
1340: }
1341: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
1342: if ( *pr )
1343: (*pr)->sugar = p->sugar;
1344: }
1345: }
1346:
1347: void adddl(int n,DL d1,DL d2,DL *dr)
1348: {
1349: DL dt;
1350: int i;
1351:
1352: *dr = dt = (DL)MALLOC_ATOMIC((n+1)*sizeof(int));
1353: dt->td = d1->td + d2->td;
1354: for ( i = 0; i < n; i++ )
1355: dt->d[i] = d1->d[i]+d2->d[i];
1356: }
1357:
1358: /* d1 += d2 */
1359:
1360: void adddl_destructive(int n,DL d1,DL d2)
1361: {
1362: int i;
1363:
1364: d1->td += d2->td;
1365: for ( i = 0; i < n; i++ )
1366: d1->d[i] += d2->d[i];
1367: }
1368:
1369: int compd(VL vl,DP p1,DP p2)
1370: {
1371: int n,t;
1372: MP m1,m2;
1373:
1374: if ( !p1 )
1375: return p2 ? -1 : 0;
1376: else if ( !p2 )
1377: return 1;
1378: else if ( NV(p1) != NV(p2) ) {
1379: error("compd : size mismatch");
1380: return 0; /* XXX */
1381: } else {
1382: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
1383: m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
1384: if ( (t = (*cmpdl)(n,m1->dl,m2->dl)) ||
1385: (t = arf_comp(vl,C(m1),C(m2)) ) )
1386: return t;
1387: if ( m1 )
1388: return 1;
1389: else if ( m2 )
1390: return -1;
1391: else
1392: return 0;
1393: }
1394: }
1395:
1396: int cmpdl_lex(int n,DL d1,DL d2)
1397: {
1398: int i;
1399:
1400: for ( i = 0; i < n && d1->d[i] == d2->d[i]; i++ );
1401: return i == n ? 0 : (d1->d[i] > d2->d[i] ? 1 : -1);
1402: }
1403:
1404: int cmpdl_revlex(int n,DL d1,DL d2)
1405: {
1406: int i;
1407:
1408: for ( i = n - 1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
1409: return i < 0 ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1410: }
1411:
1412: int cmpdl_gradlex(int n,DL d1,DL d2)
1413: {
1414: if ( d1->td > d2->td )
1415: return 1;
1416: else if ( d1->td < d2->td )
1417: return -1;
1418: else
1419: return cmpdl_lex(n,d1,d2);
1420: }
1421:
1422: int cmpdl_revgradlex(int n,DL d1,DL d2)
1423: {
1424: register int i,c;
1425: register int *p1,*p2;
1426:
1427: if ( d1->td > d2->td )
1428: return 1;
1429: else if ( d1->td < d2->td )
1430: return -1;
1431: else {
1432: i = n-1;
1433: p1 = d1->d+n-1;
1434: p2 = d2->d+n-1;
1435: while ( i >= 7 ) {
1436: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1437: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1438: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1439: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1440: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1441: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1442: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1443: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1444: i -= 8;
1445: }
1446: switch ( i ) {
1447: case 6:
1448: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1449: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1450: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1451: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1452: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1453: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1454: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1455: return 0;
1456: case 5:
1457: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1458: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1459: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1460: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1461: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1462: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1463: return 0;
1464: case 4:
1465: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1466: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1467: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1468: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1469: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1470: return 0;
1471: case 3:
1472: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1473: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1474: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1475: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1476: return 0;
1477: case 2:
1478: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1479: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1480: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1481: return 0;
1482: case 1:
1483: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1484: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1485: return 0;
1486: case 0:
1487: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1488: return 0;
1489: default:
1490: return 0;
1491: }
1492: LAST:
1493: if ( c > 0 ) return -1;
1494: else return 1;
1495: }
1496: }
1497:
1498: int cmpdl_blex(int n,DL d1,DL d2)
1499: {
1500: int c;
1501:
1502: if ( (c = cmpdl_lex(n-1,d1,d2)) )
1503: return c;
1504: else {
1505: c = d1->d[n-1] - d2->d[n-1];
1506: return c > 0 ? 1 : c < 0 ? -1 : 0;
1507: }
1508: }
1509:
1510: int cmpdl_bgradlex(int n,DL d1,DL d2)
1511: {
1512: int e1,e2,c;
1513:
1514: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
1515: if ( e1 > e2 )
1516: return 1;
1517: else if ( e1 < e2 )
1518: return -1;
1519: else {
1520: c = cmpdl_lex(n-1,d1,d2);
1521: if ( c )
1522: return c;
1523: else
1524: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
1525: }
1526: }
1527:
1528: int cmpdl_brevgradlex(int n,DL d1,DL d2)
1529: {
1530: int e1,e2,c;
1531:
1532: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
1533: if ( e1 > e2 )
1534: return 1;
1535: else if ( e1 < e2 )
1536: return -1;
1537: else {
1538: c = cmpdl_revlex(n-1,d1,d2);
1539: if ( c )
1540: return c;
1541: else
1542: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
1543: }
1544: }
1545:
1546: int cmpdl_brevrev(int n,DL d1,DL d2)
1547: {
1548: int e1,e2,f1,f2,c,i;
1549:
1550: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1551: e1 += d1->d[i]; e2 += d2->d[i];
1552: }
1553: f1 = d1->td - e1; f2 = d2->td - e2;
1554: if ( e1 > e2 )
1555: return 1;
1556: else if ( e1 < e2 )
1557: return -1;
1558: else {
1559: c = cmpdl_revlex(dp_nelim,d1,d2);
1560: if ( c )
1561: return c;
1562: else if ( f1 > f2 )
1563: return 1;
1564: else if ( f1 < f2 )
1565: return -1;
1566: else {
1567: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1568: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1569: }
1570: }
1571: }
1572:
1573: int cmpdl_bgradrev(int n,DL d1,DL d2)
1574: {
1575: int e1,e2,f1,f2,c,i;
1576:
1577: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1578: e1 += d1->d[i]; e2 += d2->d[i];
1579: }
1580: f1 = d1->td - e1; f2 = d2->td - e2;
1581: if ( e1 > e2 )
1582: return 1;
1583: else if ( e1 < e2 )
1584: return -1;
1585: else {
1586: c = cmpdl_lex(dp_nelim,d1,d2);
1587: if ( c )
1588: return c;
1589: else if ( f1 > f2 )
1590: return 1;
1591: else if ( f1 < f2 )
1592: return -1;
1593: else {
1594: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1595: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1596: }
1597: }
1598: }
1599:
1600: int cmpdl_blexrev(int n,DL d1,DL d2)
1601: {
1602: int e1,e2,f1,f2,c,i;
1603:
1604: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1605: e1 += d1->d[i]; e2 += d2->d[i];
1606: }
1607: f1 = d1->td - e1; f2 = d2->td - e2;
1608: c = cmpdl_lex(dp_nelim,d1,d2);
1609: if ( c )
1610: return c;
1611: else if ( f1 > f2 )
1612: return 1;
1613: else if ( f1 < f2 )
1614: return -1;
1615: else {
1616: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1617: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1618: }
1619: }
1620:
1621: int cmpdl_elim(int n,DL d1,DL d2)
1622: {
1623: int e1,e2,i;
1624:
1625: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1626: e1 += d1->d[i]; e2 += d2->d[i];
1627: }
1628: if ( e1 > e2 )
1629: return 1;
1630: else if ( e1 < e2 )
1631: return -1;
1632: else
1633: return cmpdl_revgradlex(n,d1,d2);
1634: }
1635:
1636: int cmpdl_weyl_elim(int n,DL d1,DL d2)
1637: {
1638: int e1,e2,i;
1639:
1640: for ( i = 1, e1 = 0, e2 = 0; i <= dp_nelim; i++ ) {
1641: e1 += d1->d[n-i]; e2 += d2->d[n-i];
1642: }
1643: if ( e1 > e2 )
1644: return 1;
1645: else if ( e1 < e2 )
1646: return -1;
1647: else if ( d1->td > d2->td )
1648: return 1;
1649: else if ( d1->td < d2->td )
1650: return -1;
1651: else return -cmpdl_revlex(n,d1,d2);
1652: }
1653:
1654: /*
1655: a special ordering
1656: 1. total order
1657: 2. (-w,w) for the first 2*m variables
1658: 3. DRL for the first 2*m variables
1659: */
1660:
1661: extern int *current_weyl_weight_vector;
1662:
1663: int cmpdl_homo_ww_drl(int n,DL d1,DL d2)
1664: {
1665: int e1,e2,m,i;
1666: int *p1,*p2;
1667:
1668: if ( d1->td > d2->td )
1669: return 1;
1670: else if ( d1->td < d2->td )
1671: return -1;
1672:
1673: m = n>>1;
1674: for ( i = 0, e1 = e2 = 0, p1 = d1->d, p2 = d2->d; i < m; i++ ) {
1675: e1 += current_weyl_weight_vector[i]*(p1[m+i] - p1[i]);
1676: e2 += current_weyl_weight_vector[i]*(p2[m+i] - p2[i]);
1677: }
1678: if ( e1 > e2 )
1679: return 1;
1680: else if ( e1 < e2 )
1681: return -1;
1682:
1683: e1 = d1->td - d1->d[n-1];
1684: e2 = d2->td - d2->d[n-1];
1685: if ( e1 > e2 )
1686: return 1;
1687: else if ( e1 < e2 )
1688: return -1;
1689:
1690: for ( i= n - 1, p1 = d1->d+n-1, p2 = d2->d+n-1;
1691: i >= 0 && *p1 == *p2; i--, p1--, p2-- );
1692: return i < 0 ? 0 : (*p1 < *p2 ? 1 : -1);
1693: }
1694:
1695: int cmpdl_drl_zigzag(int n,DL d1,DL d2)
1696: {
1697: int i,t,m;
1698: int *p1,*p2;
1699:
1700: if ( d1->td > d2->td )
1701: return 1;
1702: else if ( d1->td < d2->td )
1703: return -1;
1704: else {
1705: m = n>>1;
1706: for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) {
1707: if ( (t = p1[m+i] - p2[m+i]) ) return t > 0 ? -1 : 1;
1708: if ( (t = p1[i] - p2[i]) ) return t > 0 ? -1 : 1;
1709: }
1710: return 0;
1711: }
1712: }
1713:
1714: int cmpdl_homo_ww_drl_zigzag(int n,DL d1,DL d2)
1715: {
1716: int e1,e2,m,i,t;
1717: int *p1,*p2;
1718:
1719: if ( d1->td > d2->td )
1720: return 1;
1721: else if ( d1->td < d2->td )
1722: return -1;
1723:
1724: m = n>>1;
1725: for ( i = 0, e1 = e2 = 0, p1 = d1->d, p2 = d2->d; i < m; i++ ) {
1726: e1 += current_weyl_weight_vector[i]*(p1[m+i] - p1[i]);
1727: e2 += current_weyl_weight_vector[i]*(p2[m+i] - p2[i]);
1728: }
1729: if ( e1 > e2 )
1730: return 1;
1731: else if ( e1 < e2 )
1732: return -1;
1733:
1734: e1 = d1->td - d1->d[n-1];
1735: e2 = d2->td - d2->d[n-1];
1736: if ( e1 > e2 )
1737: return 1;
1738: else if ( e1 < e2 )
1739: return -1;
1740:
1741: for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) {
1742: if ( (t = p1[m+i] - p2[m+i]) ) return t > 0 ? -1 : 1;
1743: if ( (t = p1[i] - p2[i]) ) return t > 0 ? -1 : 1;
1744: }
1745: return 0;
1746: }
1747:
1748: int cmpdl_order_pair(int n,DL d1,DL d2)
1749: {
1750: int e1,e2,i,j,l;
1751: int *t1,*t2;
1752: int len,head;
1753: struct order_pair *pair;
1754:
1755: len = dp_current_spec->ord.block.length;
1756: if ( n != dp_current_spec->nv )
1757: error("cmpdl_order_pair : incompatible order specification");
1758: pair = dp_current_spec->ord.block.order_pair;
1759:
1760: head = 0;
1761: for ( i = 0, t1 = d1->d, t2 = d2->d; i < len; i++ ) {
1762: l = pair[i].length;
1763: switch ( pair[i].order ) {
1764: case 0:
1765: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
1766: e1 += MUL_WEIGHT(t1[j],head+j);
1767: e2 += MUL_WEIGHT(t2[j],head+j);
1768: }
1769: if ( e1 > e2 )
1770: return 1;
1771: else if ( e1 < e2 )
1772: return -1;
1773: else {
1774: for ( j = l - 1; j >= 0 && t1[j] == t2[j]; j-- );
1775: if ( j >= 0 )
1776: return t1[j] < t2[j] ? 1 : -1;
1777: }
1778: break;
1779: case 1:
1780: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
1781: e1 += MUL_WEIGHT(t1[j],head+j);
1782: e2 += MUL_WEIGHT(t2[j],head+j);
1783: }
1784: if ( e1 > e2 )
1785: return 1;
1786: else if ( e1 < e2 )
1787: return -1;
1788: else {
1789: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
1790: if ( j < l )
1791: return t1[j] > t2[j] ? 1 : -1;
1792: }
1793: break;
1794: case 2:
1795: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
1796: if ( j < l )
1797: return t1[j] > t2[j] ? 1 : -1;
1798: break;
1799: default:
1800: error("cmpdl_order_pair : invalid order"); break;
1801: }
1802: t1 += l; t2 += l; head += l;
1803: }
1804: return 0;
1805: }
1806:
1807: int cmpdl_composite(int nv,DL d1,DL d2)
1808: {
1809: int n,i,j,k,start,s,len;
1810: int *dw;
1811: struct sparse_weight *sw;
1812: struct weight_or_block *worb;
1813: int *w,*t1,*t2;
1814:
1815: n = dp_current_spec->ord.composite.length;
1816: worb = dp_current_spec->ord.composite.w_or_b;
1817: w = dp_dl_work;
1818: for ( i = 0, t1 = d1->d, t2 = d2->d; i < nv; i++ )
1819: w[i] = t1[i]-t2[i];
1820: for ( i = 0; i < n; i++, worb++ ) {
1821: len = worb->length;
1822: switch ( worb->type ) {
1823: case IS_DENSE_WEIGHT:
1824: dw = worb->body.dense_weight;
1825: for ( j = 0, s = 0; j < len; j++ )
1826: s += dw[j]*w[j];
1827: if ( s > 0 ) return 1;
1828: else if ( s < 0 ) return -1;
1829: break;
1830: case IS_SPARSE_WEIGHT:
1831: sw = worb->body.sparse_weight;
1832: for ( j = 0, s = 0; j < len; j++ )
1833: s += sw[j].value*w[sw[j].pos];
1834: if ( s > 0 ) return 1;
1835: else if ( s < 0 ) return -1;
1836: break;
1837: case IS_BLOCK:
1838: start = worb->body.block.start;
1839: switch ( worb->body.block.order ) {
1840: case 0:
1841: for ( j = 0, k = start, s = 0; j < len; j++, k++ ) {
1842: s += MUL_WEIGHT(w[k],k);
1843: }
1844: if ( s > 0 ) return 1;
1845: else if ( s < 0 ) return -1;
1846: else {
1847: for ( j = k-1; j >= start && w[j] == 0; j-- );
1848: if ( j >= start )
1849: return w[j] < 0 ? 1 : -1;
1850: }
1851: break;
1852: case 1:
1853: for ( j = 0, k = start, s = 0; j < len; j++, k++ ) {
1854: s += MUL_WEIGHT(w[k],k);
1855: }
1856: if ( s > 0 ) return 1;
1857: else if ( s < 0 ) return -1;
1858: else {
1859: for ( j = 0, k = start; j < len && w[j] == 0; j++, k++ );
1860: if ( j < len )
1861: return w[j] > 0 ? 1 : -1;
1862: }
1863: break;
1864: case 2:
1865: for ( j = 0, k = start; j < len && w[j] == 0; j++, k++ );
1866: if ( j < len )
1867: return w[j] > 0 ? 1 : -1;
1868: break;
1869: }
1870: break;
1871: }
1872: }
1873: return 0;
1874: }
1875:
1876: int cmpdl_matrix(int n,DL d1,DL d2)
1877: {
1878: int *v,*w,*t1,*t2;
1879: int s,i,j,len;
1880: int **matrix;
1881:
1882: for ( i = 0, t1 = d1->d, t2 = d2->d, w = dp_dl_work; i < n; i++ )
1883: w[i] = t1[i]-t2[i];
1884: len = dp_current_spec->ord.matrix.row;
1885: matrix = dp_current_spec->ord.matrix.matrix;
1886: for ( j = 0; j < len; j++ ) {
1887: v = matrix[j];
1888: for ( i = 0, s = 0; i < n; i++ )
1889: s += v[i]*w[i];
1890: if ( s > 0 )
1891: return 1;
1892: else if ( s < 0 )
1893: return -1;
1894: }
1895: return 0;
1896: }
1897:
1898: int cmpdl_top_weight(int n,DL d1,DL d2)
1899: {
1900: int *w;
1901: Z **mat;
1902: Z *a;
1903: mpz_t sum;
1904: int len,i,sgn,tsgn,row,k;
1905: int *t1,*t2;
1906:
1907: w = (int *)ALLOCA(n*sizeof(int));
1908: len = current_top_weight_len+3;
1909: t1 = d1->d; t2 = d2->d;
1910: for ( i = 0; i < n; i++ ) w[i] = t1[i]-t2[i];
1911: mpz_init_set_ui(sum,0);
1912: if ( OID(current_top_weight) == O_VECT ) {
1913: mat = (Z **)&BDY((VECT)current_top_weight);
1914: row = 1;
1915: } else {
1916: mat = (Z **)BDY((MAT)current_top_weight);
1917: row = ((MAT)current_top_weight)->row;
1918: }
1919: for ( k = 0; k < row; k++ ) {
1920: a = mat[k];
1921: for ( i = 0; i < n; i++ ) {
1922: if ( !a[i] || !w[i] ) continue;
1923: if ( w[i] > 0 )
1924: mpz_addmul_ui(sum,BDY(a[i]),(unsigned int)w[i]);
1925: else
1926: mpz_submul_ui(sum,BDY(a[i]),(unsigned int)(-w[i]));
1927: }
1928: sgn = mpz_sgn(sum);
1929: if ( sgn > 0 ) return 1;
1930: else if ( sgn < 0 ) return -1;
1931: }
1932: return (*cmpdl_tie_breaker)(n,d1,d2);
1933: }
1934:
1935: GeoBucket create_bucket()
1936: {
1937: GeoBucket g;
1938:
1939: g = CALLOC(1,sizeof(struct oGeoBucket));
1940: g->m = 32;
1941: return g;
1942: }
1943:
1944: int length(NODE d);
1945:
1946: void add_bucket(GeoBucket g,NODE d,int nv)
1947: {
1948: int l,k,m;
1949:
1950: l = length(d);
1951: for ( k = 0, m = 1; l > m; k++, m <<= 1 );
1952: /* 2^(k-1) < l <= 2^k */
1953: d = symb_merge(g->body[k],d,nv);
1954: for ( ; length(d) > (1<<(k)); k++ ) {
1955: g->body[k] = 0;
1956: d = symb_merge(g->body[k+1],d,nv);
1957: }
1958: g->body[k] = d;
1959: g->m = MAX(g->m,k);
1960: }
1961:
1962: DL remove_head_bucket(GeoBucket g,int nv)
1963: {
1964: int j,i,c,m;
1965: DL d;
1966:
1967: j = -1;
1968: m = g->m;
1969: for ( i = 0; i <= m; i++ ) {
1970: if ( !g->body[i] )
1971: continue;
1972: if ( j < 0 ) j = i;
1973: else {
1974: c = (*cmpdl)(nv,g->body[i]->body,g->body[j]->body);
1975: if ( c > 0 )
1976: j = i;
1977: else if ( c == 0 )
1978: g->body[i] = NEXT(g->body[i]);
1979: }
1980: }
1981: if ( j < 0 )
1982: return 0;
1983: else {
1984: d = g->body[j]->body;
1985: g->body[j] = NEXT(g->body[j]);
1986: return d;
1987: }
1988: }
1989:
1990: /* DPV functions */
1991:
1992: void adddv(VL vl,DPV p1,DPV p2,DPV *pr)
1993: {
1994: int i,len;
1995: DP *e;
1996:
1997: if ( !p1 || !p2 )
1998: error("adddv : invalid argument");
1999: else if ( p1->len != p2->len )
2000: error("adddv : size mismatch");
2001: else {
2002: len = p1->len;
2003: e = (DP *)MALLOC(p1->len*sizeof(DP));
2004: for ( i = 0; i < len; i++ )
2005: addd(vl,p1->body[i],p2->body[i],&e[i]);
2006: MKDPV(len,e,*pr);
2007: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
2008: }
2009: }
2010:
2011: void subdv(VL vl,DPV p1,DPV p2,DPV *pr)
2012: {
2013: int i,len;
2014: DP *e;
2015:
2016: if ( !p1 || !p2 )
2017: error("subdv : invalid argument");
2018: else if ( p1->len != p2->len )
2019: error("subdv : size mismatch");
2020: else {
2021: len = p1->len;
2022: e = (DP *)MALLOC(p1->len*sizeof(DP));
2023: for ( i = 0; i < len; i++ )
2024: subd(vl,p1->body[i],p2->body[i],&e[i]);
2025: MKDPV(len,e,*pr);
2026: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
2027: }
2028: }
2029:
2030: void chsgndv(DPV p1,DPV *pr)
2031: {
2032: int i,len;
2033: DP *e;
2034:
2035: if ( !p1 )
2036: error("subdv : invalid argument");
2037: else {
2038: len = p1->len;
2039: e = (DP *)MALLOC(p1->len*sizeof(DP));
2040: for ( i = 0; i < len; i++ )
2041: chsgnd(p1->body[i],&e[i]);
2042: MKDPV(len,e,*pr);
2043: (*pr)->sugar = p1->sugar;
2044: }
2045: }
2046:
2047: void muldv(VL vl,DP p1,DPV p2,DPV *pr)
2048: {
2049: int i,len;
2050: DP *e;
2051:
2052: len = p2->len;
2053: e = (DP *)MALLOC(p2->len*sizeof(DP));
2054: if ( !p1 ) {
2055: MKDPV(len,e,*pr);
2056: (*pr)->sugar = 0;
2057: } else {
2058: for ( i = 0; i < len; i++ )
2059: muld(vl,p1,p2->body[i],&e[i]);
2060: MKDPV(len,e,*pr);
2061: (*pr)->sugar = p1->sugar + p2->sugar;
2062: }
2063: }
2064:
2065: int compdv(VL vl,DPV p1,DPV p2)
2066: {
2067: int i,t,len;
2068:
2069: if ( p1->len != p2->len ) {
2070: error("compdv : size mismatch");
2071: return 0; /* XXX */
2072: } else {
2073: len = p1->len;
2074: for ( i = 0; i < len; i++ )
2075: if ( (t = compd(vl,p1->body[i],p2->body[i])) )
2076: return t;
2077: return 0;
2078: }
2079: }
2080:
2081: int ni_next(int *a,int n)
2082: {
2083: int i,j,k,kj;
2084:
2085: /* find the first nonzero a[j] */
2086: for ( j = 0; j < n && a[j] == 0; j++ );
2087: /* find the first zero a[k] after a[j] */
2088: for ( k = j; k < n && a[k] == 1; k++ );
2089: if ( k == n ) return 0;
2090: /* a[0] = 0, ... , a[j-1] = 0, a[j] = 1, ..., a[k-1] = 1, a[k] = 0 */
2091: /* a[0] = 1,..., a[k-j-2] = 1, a[k-j-1] = 0, ..., a[k-1] = 0, a[k] = 1 */
2092: kj = k-j-1;
2093: for ( i = 0; i < kj; i++ ) a[i] = 1;
2094: for ( ; i < k; i++ ) a[i] = 0;
2095: a[k] = 1;
2096: return 1;
2097: }
2098:
2099: int comp_nbm(NBM a,NBM b)
2100: {
2101: int d,i,ai,bi;
2102: unsigned int *ab,*bb;
2103:
2104: if ( a->d > b->d ) return 1;
2105: else if ( a->d < b->d ) return -1;
2106: else {
2107: d = a->d; ab = a->b; bb = b->b;
2108: #if 0
2109: w = (d+31)/32;
2110: for ( i = 0; i < w; i++ )
2111: if ( ab[i] > bb[i] ) return 1;
2112: else if ( ab[i] < bb[i] ) return -1;
2113: #else
2114: for ( i = 0; i < d; i++ ) {
2115: ai = NBM_GET(ab,i);
2116: bi = NBM_GET(bb,i);
2117: if ( ai > bi ) return 1;
2118: else if ( ai < bi ) return -1;
2119: }
2120: #endif
2121: return 0;
2122: }
2123: }
2124:
2125: NBM mul_nbm(NBM a,NBM b)
2126: {
2127: int ad,bd,d,i,j;
2128: unsigned int *ab,*bb,*mb;
2129: NBM m;
2130:
2131: ad = a->d; bd = b->d; ab = a->b; bb = b->b;
2132: d = ad + bd;
2133: NEWNBM(m); NEWNBMBDY(m,d);
2134: m->d = d; mulp(CO,a->c,b->c,&m->c); mb = m->b;
2135: j = 0;
2136: for ( i = 0; i < ad; i++, j++ )
2137: if ( NBM_GET(ab,i) ) NBM_SET(mb,j);
2138: else NBM_CLR(mb,j);
2139: for ( i = 0; i < bd; i++, j++ )
2140: if ( NBM_GET(bb,i) ) NBM_SET(mb,j);
2141: else NBM_CLR(mb,j);
2142: return m;
2143: }
2144:
2145: NBP nbmtonbp(NBM m)
2146: {
2147: NODE n;
2148: NBP u;
2149:
2150: MKNODE(n,m,0);
2151: MKNBP(u,n);
2152: return u;
2153: }
2154:
2155: /* a=c*x*rest -> a0= x*rest, ah=x, ar=rest */
2156:
2157: P separate_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar)
2158: {
2159: int i,d1;
2160: NBM t;
2161:
2162: if ( !a->d ) error("separate_nbm : invalid argument");
2163:
2164: if ( a0 ) {
2165: NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
2166: *a0 = nbmtonbp(t);
2167: }
2168:
2169: if ( ah ) {
2170: NEWNBM(t); NEWNBMBDY(t,1); t->d = 1; t->c = (P)ONE;
2171: if ( NBM_GET(a->b,0) ) NBM_SET(t->b,0);
2172: else NBM_CLR(t->b,0);
2173: *ah = nbmtonbp(t);
2174: }
2175:
2176: if ( ar ) {
2177: d1 = a->d-1;
2178: NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
2179: for ( i = 0; i < d1; i++ ) {
2180: if ( NBM_GET(a->b,i+1) ) NBM_SET(t->b,i);
2181: else NBM_CLR(t->b,i);
2182: }
2183: *ar = nbmtonbp(t);
2184: }
2185:
2186: return a->c;
2187: }
2188:
2189: /* a=c*rest*x -> a0= rest*x, ar=rest, at=x */
2190:
2191: P separate_tail_nbm(NBM a,NBP *a0,NBP *ar,NBP *at)
2192: {
2193: int i,d,d1;
2194: NBM t;
2195:
2196: if ( !(d=a->d) ) error("separate_tail_nbm : invalid argument");
2197:
2198: if ( a0 ) {
2199: NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
2200: *a0 = nbmtonbp(t);
2201: }
2202:
2203: d1 = a->d-1;
2204: if ( at ) {
2205: NEWNBM(t); NEWNBMBDY(t,1); t->d = 1; t->c = (P)ONE;
2206: if ( NBM_GET(a->b,d1) ) NBM_SET(t->b,0);
2207: else NBM_CLR(t->b,0);
2208: *at = nbmtonbp(t);
2209: }
2210:
2211: if ( ar ) {
2212: NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
2213: for ( i = 0; i < d1; i++ ) {
2214: if ( NBM_GET(a->b,i) ) NBM_SET(t->b,i);
2215: else NBM_CLR(t->b,i);
2216: }
2217: *ar = nbmtonbp(t);
2218: }
2219:
2220: return a->c;
2221: }
2222:
2223: NBP make_xky(int k)
2224: {
2225: int k1,i;
2226: NBM t;
2227:
2228: NEWNBM(t); NEWNBMBDY(t,k); t->d = k; t->c = (P)ONE;
2229: k1 = k-1;
2230: for ( i = 0; i < k1; i++ ) NBM_SET(t->b,i);
2231: NBM_CLR(t->b,i);
2232: return nbmtonbp(t);
2233: }
2234:
2235: /* a=c*x^(k-1)*y*rest -> a0= x^(k-1)*y*rest, ah=x^(k-1)*y, ar=rest */
2236:
2237: P separate_xky_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar)
2238: {
2239: int i,d1,k,k1;
2240: NBM t;
2241:
2242: if ( !a->d )
2243: error("separate_nbm : invalid argument");
2244: for ( i = 0; i < a->d && NBM_GET(a->b,i); i++ );
2245: if ( i == a->d )
2246: error("separate_nbm : invalid argument");
2247: k1 = i;
2248: k = i+1;
2249:
2250: if ( a0 ) {
2251: NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
2252: *a0 = nbmtonbp(t);
2253: }
2254:
2255: if ( ah ) {
2256: NEWNBM(t); NEWNBMBDY(t,k); t->d = k; t->c = (P)ONE;
2257: for ( i = 0; i < k1; i++ ) NBM_SET(t->b,i);
2258: NBM_CLR(t->b,i);
2259: *ah = nbmtonbp(t);
2260: }
2261:
2262: if ( ar ) {
2263: d1 = a->d-k;
2264: NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
2265: for ( i = 0; i < d1; i++ ) {
2266: if ( NBM_GET(a->b,i+k) ) NBM_SET(t->b,i);
2267: else NBM_CLR(t->b,i);
2268: }
2269: *ar = nbmtonbp(t);
2270: }
2271:
2272: return a->c;
2273: }
2274:
2275: void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp);
2276: void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp);
2277: void mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp);
2278: void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp);
2279:
2280: NBP shuffle_mul_nbm(NBM a,NBM b)
2281: {
2282: NBP u,a0,ah,ar,b0,bh,br,a1,b1,t;
2283: P ac,bc,c;
2284:
2285: if ( !a->d || !b->d )
2286: u = nbmtonbp(mul_nbm(a,b));
2287: else {
2288: ac = separate_nbm(a,&a0,&ah,&ar);
2289: bc = separate_nbm(b,&b0,&bh,&br);
2290: mulp(CO,ac,bc,&c);
2291: shuffle_mulnbp(CO,ar,b0,&t); mulnbp(CO,ah,t,&a1);
2292: shuffle_mulnbp(CO,a0,br,&t); mulnbp(CO,bh,t,&b1);
2293: addnbp(CO,a1,b1,&t); mulnbp(CO,(NBP)c,t,&u);
2294: }
2295: return u;
2296: }
2297:
2298: NBP harmonic_mul_nbm(NBM a,NBM b)
2299: {
2300: NBP u,a0,ah,ar,b0,bh,br,a1,b1,t,s,abk,ab1;
2301: P ac,bc,c;
2302:
2303: if ( !a->d || !b->d )
2304: u = nbmtonbp(mul_nbm(a,b));
2305: else {
2306: mulp(CO,a->c,b->c,&c);
2307: ac = separate_xky_nbm(a,&a0,&ah,&ar);
2308: bc = separate_xky_nbm(b,&b0,&bh,&br);
2309: mulp(CO,ac,bc,&c);
2310: harmonic_mulnbp(CO,ar,b0,&t); mulnbp(CO,ah,t,&a1);
2311: harmonic_mulnbp(CO,a0,br,&t); mulnbp(CO,bh,t,&b1);
2312: abk = make_xky(((NBM)BDY(BDY(ah)))->d+((NBM)BDY(BDY(bh)))->d);
2313: harmonic_mulnbp(CO,ar,br,&t); mulnbp(CO,abk,t,&ab1);
2314: addnbp(CO,a1,b1,&t); addnbp(CO,t,ab1,&s); mulnbp(CO,(NBP)c,s,&u);
2315: }
2316: return u;
2317:
2318: }
2319:
2320: void addnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2321: {
2322: NODE b1,b2,br=0,br0;
2323: NBM m1,m2,m;
2324: P c;
2325:
2326: if ( !p1 )
2327: *rp = p2;
2328: else if ( !p2 )
2329: *rp = p1;
2330: else {
2331: for ( b1 = BDY(p1), b2 = BDY(p2), br0 = 0; b1 && b2; ) {
2332: m1 = (NBM)BDY(b1); m2 = (NBM)BDY(b2);
2333: switch ( comp_nbm(m1,m2) ) {
2334: case 0:
2335: addp(CO,m1->c,m2->c,&c);
2336: if ( c ) {
2337: NEXTNODE(br0,br);
2338: NEWNBM(m); m->d = m1->d; m->c = c; m->b = m1->b;
2339: BDY(br) = (pointer)m;
2340: }
2341: b1 = NEXT(b1); b2 = NEXT(b2); break;
2342: case 1:
2343: NEXTNODE(br0,br); BDY(br) = BDY(b1);
2344: b1 = NEXT(b1); break;
2345: case -1:
2346: NEXTNODE(br0,br); BDY(br) = BDY(b2);
2347: b2 = NEXT(b2); break;
2348: }
2349: }
2350: if ( !br0 )
2351: if ( b1 )
2352: br0 = b1;
2353: else if ( b2 )
2354: br0 = b2;
2355: else {
2356: *rp = 0;
2357: return;
2358: }
2359: else if ( b1 )
2360: NEXT(br) = b1;
2361: else if ( b2 )
2362: NEXT(br) = b2;
2363: else
2364: NEXT(br) = 0;
2365: MKNBP(*rp,br0);
2366: }
2367: }
2368:
2369: void subnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2370: {
2371: NBP t;
2372:
2373: chsgnnbp(p2,&t);
2374: addnbp(vl,p1,t,rp);
2375: }
2376:
2377: void chsgnnbp(NBP p,NBP *rp)
2378: {
2379: NODE r0,r=0,b;
2380: NBM m,m1;
2381:
2382: for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
2383: NEXTNODE(r0,r);
2384: m = (NBM)BDY(b);
2385: NEWNBM(m1); m1->d = m->d; m1->b = m->b; chsgnp(m->c,&m1->c);
2386: BDY(r) = m1;
2387: }
2388: if ( r0 ) NEXT(r) = 0;
2389: MKNBP(*rp,r0);
2390: }
2391:
2392: void mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2393: {
2394: NODE b,n;
2395: NBP r,t,s;
2396: NBM m;
2397:
2398: if ( !p1 || !p2 ) {
2399: *rp = 0; return;
2400: }
2401: if ( OID(p1) != O_NBP ) {
2402: if ( !POLY(p1) )
2403: error("mulnbp : invalid argument");
2404: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
2405: MKNODE(n,m,0); MKNBP(p1,n);
2406: }
2407: if ( OID(p2) != O_NBP ) {
2408: if ( !POLY(p2) )
2409: error("mulnbp : invalid argument");
2410: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
2411: MKNODE(n,m,0); MKNBP(p2,n);
2412: }
2413: if ( length(BDY(p1)) < length(BDY(p2)) ) {
2414: for ( r = 0, b = BDY(p1); b; b = NEXT(b) ) {
2415: mulnbmnbp(vl,(NBM)BDY(b),p2,&t);
2416: addnbp(vl,r,t,&s); r = s;
2417: }
2418: *rp = r;
2419: } else {
2420: for ( r = 0, b = BDY(p2); b; b = NEXT(b) ) {
2421: mulnbpnbm(vl,p1,(NBM)BDY(b),&t);
2422: addnbp(vl,r,t,&s); r = s;
2423: }
2424: *rp = r;
2425: }
2426: }
2427:
2428: void mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp)
2429: {
2430: NODE b,r0,r=0;
2431:
2432: if ( !p ) *rp = 0;
2433: else {
2434: for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
2435: NEXTNODE(r0,r);
2436: BDY(r) = mul_nbm(m,(NBM)BDY(b));
2437: }
2438: if ( r0 ) NEXT(r) = 0;
2439: MKNBP(*rp,r0);
2440: }
2441: }
2442:
2443: void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp)
2444: {
2445: NODE b,r0,r=0;
2446:
2447: if ( !p ) *rp = 0;
2448: else {
2449: for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
2450: NEXTNODE(r0,r);
2451: BDY(r) = mul_nbm((NBM)BDY(b),m);
2452: }
2453: if ( r0 ) NEXT(r) = 0;
2454: MKNBP(*rp,r0);
2455: }
2456: }
2457:
2458: void pwrnbp(VL vl,NBP a,Z q,NBP *c)
2459: {
2460: NBP a1,a2;
2461: Z q1,r1,two;
2462: NBM m;
2463: NODE r;
2464:
2465: if ( !q ) {
2466: NEWNBM(m); m->d = 0; m->c = (P)ONE; m->b = 0;
2467: MKNODE(r,m,0); MKNBP(*c,r);
2468: } else if ( !a )
2469: *c = 0;
2470: else if ( UNIQ(q) )
2471: *c = a;
2472: else {
1.2 noro 2473: STOZ(2,two);
1.1 noro 2474: divqrz(q,two,&q1,&r1);
2475: pwrnbp(vl,a,q1,&a1);
2476: mulnbp(vl,a1,a1,&a2);
2477: if ( r1 )
2478: mulnbp(vl,a2,a,c);
2479: else
2480: *c = a2;
2481: }
2482: }
2483:
2484: int compnbp(VL vl,NBP p1,NBP p2)
2485: {
2486: NODE n1,n2;
2487: NBM m1,m2;
2488: int t;
2489:
2490: if ( !p1 )
2491: return p2 ? -1 : 0;
2492: else if ( !p2 )
2493: return 1;
2494: else {
2495: for ( n1 = BDY(p1), n2 = BDY(p2);
2496: n1 && n2; n1 = NEXT(n1), n2 = NEXT(n2) ) {
2497: m1 = (NBM)BDY(n1); m2 = (NBM)BDY(n2);
2498: if ( (t = comp_nbm(m1,m2)) || (t = compp(CO,m1->c,m2->c) ) )
2499: return t;
2500: }
2501: if ( n1 )
2502: return 1;
2503: else if ( n2 )
2504: return -1;
2505: else
2506: return 0;
2507: }
2508: }
2509:
2510: void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2511: {
2512: NODE b1,b2,n;
2513: NBP r,t,s;
2514: NBM m;
2515:
2516: if ( !p1 || !p2 ) {
2517: *rp = 0; return;
2518: }
2519: if ( OID(p1) != O_NBP ) {
2520: if ( !POLY(p1) )
2521: error("shuffle_mulnbp : invalid argument");
2522: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
2523: MKNODE(n,m,0); MKNBP(p1,n);
2524: }
2525: if ( OID(p2) != O_NBP ) {
2526: if ( !POLY(p2) )
2527: error("shuffle_mulnbp : invalid argument");
2528: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
2529: MKNODE(n,m,0); MKNBP(p2,n);
2530: }
2531: for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) )
2532: for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) {
2533: t = shuffle_mul_nbm(m,(NBM)BDY(b2));
2534: addnbp(vl,r,t,&s); r = s;
2535: }
2536: *rp = r;
2537: }
2538:
2539: void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2540: {
2541: NODE b1,b2,n;
2542: NBP r,t,s;
2543: NBM m;
2544:
2545: if ( !p1 || !p2 ) {
2546: *rp = 0; return;
2547: }
2548: if ( OID(p1) != O_NBP ) {
2549: if ( !POLY(p1) )
2550: error("harmonic_mulnbp : invalid argument");
2551: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
2552: MKNODE(n,m,0); MKNBP(p1,n);
2553: }
2554: if ( OID(p2) != O_NBP ) {
2555: if ( !POLY(p2) )
2556: error("harmonic_mulnbp : invalid argument");
2557: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
2558: MKNODE(n,m,0); MKNBP(p2,n);
2559: }
2560: for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) )
2561: for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) {
2562: t = harmonic_mul_nbm(m,(NBM)BDY(b2));
2563: addnbp(vl,r,t,&s); r = s;
2564: }
2565: *rp = r;
2566: }
2567:
2568: #if 0
2569: NBP shuffle_mul_nbm(NBM a,NBM b)
2570: {
2571: int ad,bd,d,i,ai,bi,bit,s;
2572: int *ab,*bb,*wmb,*w;
2573: NBM wm,tm;
2574: P c,c1;
2575: NODE r,t,t1,p;
2576: NBP u;
2577:
2578: ad = a->d; bd = b->d; ab = a->b; bb = b->b;
2579: d = ad + bd;
2580: w = (int *)ALLOCA(d*sizeof(int));
2581: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2582: for ( i = 0; i < ad; i++ ) w[i] = 1;
2583: for ( ; i < d; i++ ) w[i] = 0;
2584: mulp(CO,a->c,b->c,&c);
2585: r = 0;
2586: do {
2587: wm->d = d; wm->c = c;
2588: ai = 0; bi = 0;
2589: for ( i = 0; i < d; i++ ) {
2590: if ( w[i] ) { bit = NBM_GET(ab,ai); ai++; }
2591: else { bit = NBM_GET(bb,bi); bi++; }
2592: if ( bit ) NBM_SET(wmb,i);
2593: else NBM_CLR(wmb,i);
2594: }
2595: for ( p = 0, t = r; t; p = t, t = NEXT(t) ) {
2596: tm = (NBM)BDY(t);
2597: s = comp_nbm(tm,wm);
2598: if ( s < 0 ) {
2599: /* insert */
2600: MKNODE(t1,wm,t);
2601: if ( !p ) r = t1;
2602: else NEXT(p) = t1;
2603: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2604: break;
2605: } else if ( s == 0 ) {
2606: /* add coefs */
2607: addp(CO,tm->c,c,&c1);
2608: if ( c1 ) tm->c = c1;
2609: else NEXT(p) = NEXT(t);
2610: break;
2611: }
2612: }
2613: if ( !t ) {
2614: /* append */
2615: MKNODE(t1,wm,t);
2616: if ( !p ) r = t1;
2617: else NEXT(p) = t1;
2618: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2619: }
2620: } while ( ni_next(w,d) );
2621: MKNBP(u,r);
2622: return u;
2623: }
2624:
2625: int nbmtoxky(NBM a,int *b)
2626: {
2627: int d,i,j,k;
2628: int *p;
2629:
2630: d = a->d; p = a->b;
2631: for ( i = j = 0, k = 1; i < d; i++ ) {
2632: if ( !NBM_GET(p,i) ) {
2633: b[j++] = k;
2634: k = 1;
2635: } else k++;
2636: }
2637: return j;
2638: }
2639:
2640: NBP harmonic_mul_nbm(NBM a,NBM b)
2641: {
2642: int da,db,d,la,lb,lmax,lmin,l,lab,la1,lb1,lab1;
2643: int i,j,k,ia,ib,s;
2644: int *wa,*wb,*w,*wab,*wa1,*wmb;
2645: P c,c1;
2646: NBM wm,tm;
2647: NODE r,t1,t,p;
2648: NBP u;
2649:
2650: da = a->d; db = b->d; d = da+db;
2651: wa = (int *)ALLOCA(da*sizeof(int));
2652: wb = (int *)ALLOCA(db*sizeof(int));
2653: la = nbmtoxky(a,wa);
2654: lb = nbmtoxky(b,wb);
2655: mulp(CO,a->c,b->c,&c);
2656: /* wa[0],..,wa[la-1] <-> x^wa[0]y x^wa[1]y .. */
2657: /* lmax : total length */
2658: lmax = la+lb;
2659: lmin = la>lb?la:lb;
2660: w = (int *)ALLOCA(lmax*sizeof(int));
2661: /* position of a+b */
2662: wab = (int *)ALLOCA(lmax*sizeof(int));
2663: /* position of a */
2664: wa1 = (int *)ALLOCA(lmax*sizeof(int));
2665: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2666: for ( l = lmin, r = 0; l <= lmax; l++ ) {
2667: lab = lmax - l;
2668: la1 = la - lab;
2669: lb1 = lb - lab;
2670: lab1 = l-lab;
2671: /* partion l into three parts: a, b, a+b */
2672: /* initialize wab */
2673: for ( i = 0; i < lab; i++ ) wab[i] = 1;
2674: for ( ; i < l; i++ ) wab[i] = 0;
2675: do {
2676: /* initialize wa1 */
2677: for ( i = 0; i < la1; i++ ) wa1[i] = 1;
2678: for ( ; i < lab1; i++ ) wa1[i] = 0;
2679: do {
2680: ia = 0; ib = 0;
2681: for ( i = j = 0; i < l; i++ )
2682: if ( wab[i] ) w[i] = wa[ia++]+wb[ib++];
2683: else if ( wa1[j++] ) w[i] = wa[ia++];
2684: else w[i] = wb[ib++];
2685: for ( i = j = 0; i < l; i++ ) {
2686: for ( k = w[i]-1; k > 0; k--, j++ ) NBM_SET(wmb,j);
2687: NBM_CLR(wmb,j); j++;
2688: }
2689: wm->d = j; wm->c = c;
2690: for ( p = 0, t = r; t; p = t, t = NEXT(t) ) {
2691: tm = (NBM)BDY(t);
2692: s = comp_nbm(tm,wm);
2693: if ( s < 0 ) {
2694: /* insert */
2695: MKNODE(t1,wm,t);
2696: if ( !p ) r = t1;
2697: else NEXT(p) = t1;
2698: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2699: break;
2700: } else if ( s == 0 ) {
2701: /* add coefs */
2702: addp(CO,tm->c,c,&c1);
2703: if ( c1 ) tm->c = c1;
2704: else NEXT(p) = NEXT(t);
2705: break;
2706: }
2707: }
2708: if ( !t ) {
2709: /* append */
2710: MKNODE(t1,wm,t);
2711: if ( !p ) r = t1;
2712: else NEXT(p) = t1;
2713: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2714: }
2715: } while ( ni_next(wa1,lab1) );
2716: } while ( ni_next(wab,l) );
2717: }
2718: MKNBP(u,r);
2719: return u;
2720: }
2721: #endif
2722:
2723: /* DPM functions */
2724:
1.3 noro 2725: DMMstack dmm_stack;
1.9 noro 2726: int dpm_base_ordtype;;
1.3 noro 2727:
1.9 noro 2728: DMMstack push_schreyer_order(LIST data,DMMstack stack)
1.3 noro 2729: {
2730: DMMstack t;
1.14 noro 2731: DP dp;
2732: MP mp;
2733: DMM m0,m1;
2734: DPM dpm0,dpm1;
2735: int len,i,nv;
1.7 noro 2736: NODE in,t1;
1.9 noro 2737: LIST l;
1.4 noro 2738:
2739: /* data = [DPM,...,DPM] */
1.9 noro 2740: if ( !dmm_stack && ( !dp_current_spec || dp_current_spec->id < 256 ) )
2741: error("push_schreyer_order : base module order is not set");
1.4 noro 2742: in = BDY(data);
2743: len = length(in);
2744: NEWDMMstack(t);
2745: t->rank = len;
2746: t->in = (DMM *)MALLOC((len+1)*sizeof(DMM));
1.14 noro 2747: t->sum = (DMM *)MALLOC((len+1)*sizeof(DMM));
1.9 noro 2748: if ( stack ) {
2749: MKNODE(t1,data,BDY(stack->obj)); MKLIST(l,t1); t->obj = l;
1.14 noro 2750: for ( i = 1; i <= len; i++, in = NEXT(in) ) {
2751: m1 = t->in[i] = BDY((DPM)BDY(in));
2752: NEWMP(mp); mp->dl = m1->dl; mp->c = m1->c; NEXT(mp) = 0;
2753: nv = ((DPM)BDY(in))->nv;
2754: MKDP(nv,mp,dp); dp->sugar = mp->dl->td;
2755: m0 = stack->sum[m1->pos]; MKDPM(nv,m0,dpm0);
2756: mulobjdpm(CO,(Obj)dp,dpm0,&dpm1);
2757: t->sum[i] = BDY(dpm1);
2758: }
1.9 noro 2759: } else {
2760: MKNODE(t1,data,0); MKLIST(l,t1); t->obj = l;
1.14 noro 2761: for ( i = 1; i <= len; i++, in = NEXT(in) ) {
2762: t->sum[i] = t->in[i] = BDY((DPM)BDY(in));
2763: }
1.4 noro 2764: }
1.9 noro 2765: t->next = stack;;
2766: dpm_ordtype = 3;
2767: return t;
1.4 noro 2768: }
2769:
2770: // data=[Ink,...,In0]
2771: // Ini = a list of module monomials
2772:
2773: void set_schreyer_order(LIST data)
2774: {
2775: NODE in;
2776: LIST *w;
2777: int i,len;
1.3 noro 2778:
2779: if ( !data ) {
2780: dmm_stack = 0;
2781: if ( dp_current_spec && dp_current_spec->id >= 256 )
1.9 noro 2782: dpm_ordtype = dp_current_spec->module_ordtype;
1.3 noro 2783: else
2784: dpm_ordtype = 0;
2785: return;
2786: } else {
1.9 noro 2787: if ( !dp_current_spec || dp_current_spec->id < 256 )
2788: error("set_schreyer_order : base module order is not set");
1.4 noro 2789: dmm_stack = 0;
1.9 noro 2790: dpm_base_ordtype = dp_current_spec->module_ordtype;
1.4 noro 2791: in = BDY(data);
1.3 noro 2792: len = length(in);
1.4 noro 2793: w = (LIST *)MALLOC(len*sizeof(LIST));
2794: for ( i = 0; i < len; i++, in = NEXT(in) ) w[i] = (LIST)BDY(in);
1.9 noro 2795: for ( i = len-1; i >= 0; i-- ) dmm_stack = push_schreyer_order(w[i],dmm_stack);
1.3 noro 2796: }
2797: }
2798:
1.15 noro 2799: void set_schreyer_level(DMMstack_array array,int level)
2800: {
2801: if ( !level ) {
2802: dmm_stack = 0;
2803: if ( dp_current_spec && dp_current_spec->id >= 256 )
2804: dpm_ordtype = dp_current_spec->module_ordtype;
2805: else
2806: dpm_ordtype = 0;
2807: return;
2808: } else {
2809: if ( !dp_current_spec || dp_current_spec->id < 256 )
2810: error("set_schreyer_level : base module order is not set");
2811: if ( level > array->len )
2812: error("set_schreyer_level : invalid level");
2813: dmm_stack = array->body[level-1];
2814: dpm_base_ordtype = dp_current_spec->module_ordtype;
2815: dpm_ordtype = 3;
2816: }
2817: }
2818:
1.4 noro 2819: // construct a base of syz(g)
2820: // assuming the schrerer order is properly set
2821:
2822: DP dpm_sp_hm(DPM p1,DPM p2);
2823: void dpm_sp(DPM p1,DPM p2,DPM *sp,DP *t1,DP *t2);
1.12 noro 2824: DPM dpm_nf_and_quotient3(DPM sp,VECT psv,DPM *nf,P *dn);
2825: DPM dpm_nf_and_quotient4(DPM sp,DPM *ps,VECT psiv,DPM head,DPM *nf,P *dn);
2826: DPM dpm_sp_nf(VECT psv,VECT psiv,int i,int j,DPM *nf);
1.19 noro 2827: DPM dpm_sp_nf_zlist(VECT psv,VECT psiv,int i,int j,int top,DPM *nf);
1.13 noro 2828: DPM dpm_sp_nf_asir(VECT psv,int i,int j,DPM *nf);
1.4 noro 2829: void dpm_sort(DPM p,DPM *r);
2830:
2831: extern int DP_Multiple;
1.22 noro 2832: extern int DP_Print;
1.4 noro 2833:
2834: void dpm_nf_z(NODE b,DPM g,VECT psv,int full,int multiple,DPM *rp);
2835: NODE dpm_sort_list(NODE l);
2836: void dpm_ptozp(DPM p,Z *cont,DPM *r);
2837:
2838: NODE dpm_reduceall(NODE in)
2839: {
2840: int n,i;
2841: VECT psv;
2842: DPM *ps;
2843: NODE t,t1;
2844: DPM g,r;
2845: Z cont;
2846:
2847: n = length(in);
2848: MKVECT(psv,n);
2849: ps = (DPM *)BDY(psv);
2850: for ( i = 0, t = in; i < n; i++, t = NEXT(t) ) ps[i] = BDY(t);
2851: for ( i = 0; i < n; i++ ) {
2852: g = ps[i]; ps[i] = 0;
1.8 noro 2853: // dpm_nf_z(0,g,psv,1,DP_Multiple,&ps[i]);
2854: dpm_nf_z(0,g,psv,1,0,&ps[i]);
1.4 noro 2855: }
2856: t = 0;
2857: for ( i = n-1; i >= 0; i-- ) {
2858: dpm_ptozp(ps[i],&cont,&r);
2859: MKNODE(t1,r,t); t = t1;
2860: }
2861: return t;
2862: }
2863:
1.8 noro 2864: struct oEGT egra;
2865:
1.9 noro 2866: void dpm_ht(DPM d,DPM *r);
1.16 noro 2867: NODE reverse_node(NODE b);
1.9 noro 2868:
1.12 noro 2869: void dpm_schreyer_base(LIST g,LIST *s)
2870: {
2871: NODE nd,t0,t,b0,b;
2872: int n,i,j,k,nv,max,pos;
2873: LIST l;
2874: DP h,t1,t2;
2875: MP d;
2876: DMM r0,r,r1;
2877: DPM sp,nf,dpm;
2878: DPM *ps;
2879: VECT psv,psiv;
2880: DPM quo;
1.16 noro 2881: DP *mm;
1.12 noro 2882: NODE *psi;
1.13 noro 2883: NODE n1,n2,n3;
2884: int p1,p2,p3;
2885: struct oEGT eg0,eg1,egsp,egnf;
1.12 noro 2886: extern struct oEGT egred;
1.15 noro 2887: extern int sch_count,schrec_count,schlast_count;
1.12 noro 2888:
1.15 noro 2889: sch_count = schlast_count= 0;
1.12 noro 2890: init_eg(&egra);
1.13 noro 2891: init_eg(&egsp);
2892: init_eg(&egnf);
1.12 noro 2893: nd = BDY(g);
2894: n = length(nd);
2895: MKVECT(psv,n+1);
2896: ps = (DPM *)BDY(psv);
2897: for ( i = 1, t = nd; i <= n; i++, t = NEXT(t) ) ps[i] = (DPM)BDY(t);
2898: for ( i = 1, max = 0; i <= n; i++ )
2899: if ( (pos=BDY(ps[i])->pos) > max ) max = pos;
2900: MKVECT(psiv,max+1);
2901: psi = (NODE *)BDY(psiv);
2902: for ( i = n; i >= 1; i-- ) {
2903: pos = BDY(ps[i])->pos;
2904: MKNODE(nd,(long)i,psi[pos]); psi[pos] = nd;
2905: }
2906: nv = ps[1]->nv;
1.16 noro 2907: mm = (DP *)MALLOC((n+1)*sizeof(DP));
1.12 noro 2908: b0 = 0;
1.13 noro 2909: get_eg(&eg0);
2910: for ( i = 1; i <= max; i++ ) {
1.16 noro 2911: memset(mm,0,(n+1)*sizeof(DP));
1.13 noro 2912: for ( n1 = psi[i]; n1; n1 = NEXT(n1) ) {
2913: p1 = (long)BDY(n1);
2914: for ( n2 = NEXT(n1); n2; n2 = NEXT(n2) ) {
2915: p2 = (long)BDY(n2);
1.16 noro 2916: mm[p2] = dpm_sp_hm(ps[p1],ps[p2]);
1.13 noro 2917: }
2918: for ( n2 = NEXT(n1); n2; n2 = NEXT(n2) ) {
2919: p2 = (long)BDY(n2);
1.16 noro 2920: if ( !mm[p2] ) continue;
2921: for ( h = mm[p2], n3 = NEXT(n1); n3; n3 = NEXT(n3) ) {
1.13 noro 2922: p3 = (long)BDY(n3);
1.16 noro 2923: if ( n3 != n2 && mm[p3] && dp_redble(mm[p3],h) ) mm[p3] = 0;
2924: }
2925: }
2926: for ( j = p1+1; j <= n; j++ ) {
2927: if ( mm[j] ) {
2928: quo = dpm_sp_nf(psv,psiv,p1,j,&nf);
2929: if ( nf )
2930: error("dpm_schreyer_base : cannot happen");
2931: NEXTNODE(b0,b); BDY(b) = (pointer)quo;
1.13 noro 2932: }
2933: }
2934: }
2935: }
2936: get_eg(&eg1); add_eg(&egsp,&eg0,&eg1); print_eg("SP",&egsp);
2937: get_eg(&eg0);
2938: get_eg(&eg1); add_eg(&egnf,&eg0,&eg1); print_eg("NF",&egnf); printf("\n");
1.12 noro 2939: if ( b0 ) NEXT(b) = 0;
1.17 noro 2940: for ( t0 = 0, nd = BDY(g); nd; nd = NEXT(nd) ) {
1.12 noro 2941: dpm_ht((DPM)BDY(nd),&dpm); NEXTNODE(t0,t); BDY(t) = (pointer)dpm;
2942: }
2943: if ( t0 ) NEXT(t) = 0;
2944: MKLIST(l,t0);
2945: dmm_stack = push_schreyer_order(l,dmm_stack);
1.16 noro 2946: // b0 = dpm_sort_list(b0);
1.12 noro 2947: // get_eg(&eg0);
2948: // b0 = dpm_reduceall(b0);
2949: // get_eg(&eg1); add_eg(&egra,&eg0,&eg1); print_eg("RA",&egra);
1.4 noro 2950: MKLIST(*s,b0);
1.8 noro 2951: // print_eg("red",&egred); printf("\n");
1.15 noro 2952: printf("sch_count=%d, schlast_count=%d\n",sch_count,schlast_count);
2953: }
2954:
1.17 noro 2955: void dpm_list_to_array(LIST g,VECT *psvect,VECT *psindvect)
2956: {
2957: NODE nd,t;
2958: int n;
2959: VECT psv,psiv;
2960: DPM *ps;
2961: int i,max,pos;
2962: LIST *psi;
2963: LIST l;
2964: Z iz;
2965:
2966: nd = BDY(g);
2967: n = length(nd);
2968: MKVECT(psv,n+1);
2969: ps = (DPM *)BDY(psv);
2970: for ( i = 1, t = nd; i <= n; i++, t = NEXT(t) ) ps[i] = (DPM)BDY(t);
2971: for ( i = 1, max = 0; i <= n; i++ )
2972: if ( (pos=BDY(ps[i])->pos) > max ) max = pos;
2973: MKVECT(psiv,max+1);
2974: psi = (LIST *)BDY(psiv);
2975: for ( i = 0; i <= max; i++ ) {
2976: MKLIST(l,0); psi[i] = l;
2977: }
2978: for ( i = n; i >= 1; i-- ) {
2979: pos = BDY(ps[i])->pos;
2980: STOZ(i,iz); MKNODE(nd,iz,BDY(psi[pos])); BDY(psi[pos]) = nd;
2981: }
2982: *psvect = psv; *psindvect = psiv;
2983: }
2984:
1.19 noro 2985: #if 0
1.17 noro 2986: void dpm_insert_to_zlist(VECT psiv,int pos,int i)
2987: {
2988: LIST l;
2989: NODE prev,cur,nd;
2990: int curi;
2991: Z iz;
2992:
2993: l = (LIST)BDY(psiv)[pos];
2994: for ( prev = 0, cur = BDY(l); cur; cur = NEXT(cur) ) {
2995: curi = ZTOS((Q)BDY(cur));
2996: if ( curi == i )
2997: error("dpm_insert_to_list : invalid index");
2998: if ( i < curi ) break;
2999: prev = cur;
3000: }
3001: STOZ(i,iz); MKNODE(nd,iz,cur);
3002: if ( prev == 0 ) BDY(l) = nd;
3003: else NEXT(prev) = nd;
3004: }
1.19 noro 3005: #else
3006: void dpm_insert_to_zlist(VECT psiv,int pos,int i)
3007: {
3008: LIST l;
3009: NODE prev,cur,nd;
3010: int curi;
3011: Z iz;
3012:
3013: l = (LIST)BDY(psiv)[pos];
3014: for ( prev = 0, cur = BDY(l); cur; cur = NEXT(cur) ) {
3015: curi = ZTOS((Q)BDY(cur));
3016: if ( curi == i )
3017: error("dpm_insert_to_list : invalid index");
3018: prev = cur;
3019: }
3020: STOZ(i,iz); MKNODE(nd,iz,cur);
3021: if ( prev == 0 ) BDY(l) = nd;
3022: else NEXT(prev) = nd;
3023: }
3024: #endif
1.17 noro 3025:
3026: void dpm_schreyer_base_zlist(LIST g,LIST *s)
3027: {
3028: NODE nd,t0,t,b0,b;
3029: int n,i,j,k,nv,max,pos;
3030: LIST l;
3031: DP h,t1,t2;
3032: MP d;
3033: DMM r0,r,r1;
3034: DPM sp,nf,dpm;
3035: DPM *ps;
3036: VECT psv,psiv;
3037: DPM quo;
3038: DP *mm;
3039: LIST *psi;
3040: NODE n1,n2,n3;
3041: int p1,p2,p3;
3042: Z iz;
3043: struct oEGT eg0,eg1,egsp,egnf;
3044: extern struct oEGT egred;
3045: extern int sch_count,schrec_count,schlast_count;
3046:
3047: sch_count = schlast_count= 0;
3048: init_eg(&egra);
3049: init_eg(&egsp);
3050: init_eg(&egnf);
3051: dpm_list_to_array(g,&psv,&psiv);
3052: ps = (DPM *)BDY(psv);
3053: psi = (LIST *)BDY(psiv);
3054: nv = ps[1]->nv;
3055: n = psv->len-1;
3056: max = psiv->len-1;
3057: mm = (DP *)MALLOC((n+1)*sizeof(DP));
3058: b0 = 0;
3059: get_eg(&eg0);
3060: for ( i = 1; i <= max; i++ ) {
3061: memset(mm,0,(n+1)*sizeof(DP));
3062: for ( n1 = BDY((LIST)psi[i]); n1; n1 = NEXT(n1) ) {
3063: p1 = ZTOS((Q)BDY(n1));
3064: for ( n2 = NEXT(n1); n2; n2 = NEXT(n2) ) {
3065: p2 = ZTOS((Q)BDY(n2));
3066: mm[p2] = dpm_sp_hm(ps[p1],ps[p2]);
3067: }
3068: for ( n2 = NEXT(n1); n2; n2 = NEXT(n2) ) {
3069: p2 = ZTOS((Q)BDY(n2));
3070: if ( !mm[p2] ) continue;
3071: for ( h = mm[p2], n3 = NEXT(n1); n3; n3 = NEXT(n3) ) {
3072: p3 = ZTOS((Q)BDY(n3));
3073: if ( n3 != n2 && mm[p3] && dp_redble(mm[p3],h) ) mm[p3] = 0;
3074: }
3075: }
3076: for ( j = p1+1; j <= n; j++ ) {
3077: if ( mm[j] ) {
1.19 noro 3078: quo = dpm_sp_nf_zlist(psv,psiv,p1,j,0,&nf);
1.17 noro 3079: if ( nf )
3080: error("dpm_schreyer_base : cannot happen");
3081: NEXTNODE(b0,b); BDY(b) = (pointer)quo;
3082: }
3083: }
3084: }
3085: }
3086: get_eg(&eg1); add_eg(&egsp,&eg0,&eg1); print_eg("SP",&egsp);
3087: get_eg(&eg0);
3088: get_eg(&eg1); add_eg(&egnf,&eg0,&eg1); print_eg("NF",&egnf); printf("\n");
3089: if ( b0 ) NEXT(b) = 0;
3090: for ( t0 = 0, nd = BDY(g); nd; nd = NEXT(nd) ) {
3091: dpm_ht((DPM)BDY(nd),&dpm); NEXTNODE(t0,t); BDY(t) = (pointer)dpm;
3092: }
3093: if ( t0 ) NEXT(t) = 0;
3094: MKLIST(l,t0);
3095: dmm_stack = push_schreyer_order(l,dmm_stack);
3096: // b0 = dpm_sort_list(b0);
3097: // get_eg(&eg0);
3098: // b0 = dpm_reduceall(b0);
3099: // get_eg(&eg1); add_eg(&egra,&eg0,&eg1); print_eg("RA",&egra);
3100: MKLIST(*s,b0);
3101: // print_eg("red",&egred); printf("\n");
3102: printf("sch_count=%d, schlast_count=%d\n",sch_count,schlast_count);
3103: }
1.15 noro 3104:
1.19 noro 3105: static int compdp_nv;
3106:
3107: int compdp_revgradlex(DP *a,DP *b)
3108: {
3109: return -cmpdl_revgradlex(compdp_nv,BDY(*a)->dl,BDY(*b)->dl);
3110: }
3111:
3112: int compdp_lex(DP *a,DP *b)
3113: {
3114: return -cmpdl_lex(compdp_nv,BDY(*a)->dl,BDY(*b)->dl);
3115: }
3116:
1.22 noro 3117: DMMstack_array dpm_schreyer_frame(NODE g,int lex)
1.15 noro 3118: {
3119: LIST l;
3120: NODE nd,in,b0,b,n1,n2,n3,t;
3121: NODE *psi;
3122: long p1,p2,p3;
1.19 noro 3123: int nv,n,i,j,k,max,pos,level;
1.15 noro 3124: DMMstack s,s1;
3125: DMM m1,m0,dmm;
3126: MP mp;
3127: DP dp,h;
1.19 noro 3128: DP *m,*m2;
1.15 noro 3129: DPM dpm,dpm0,dpm1;
3130: VECT psv,psiv;
3131: DPM *ps;
3132: DMMstack_array dmmstack_array;
3133:
3134: nd = g;
1.19 noro 3135: compdp_nv = nv = ((DPM)BDY(nd))->nv;
1.15 noro 3136: s = 0;
3137: level = 0;
3138: while ( 1 ) {
3139: /* store the current nd to a DMMstack */
3140: n = length(nd);
3141: NEWDMMstack(s1);
3142: MKLIST(l,nd); s1->obj = l;
3143: s1->rank = n;
3144: s1->in = (DMM *)MALLOC((n+1)*sizeof(DMM));
3145: s1->sum = (DMM *)MALLOC((n+1)*sizeof(DMM));
3146: NEXT(s1) = s;
3147: if ( s ) {
3148: for ( i = 1, in = nd; i <= n; i++, in = NEXT(in) ) {
3149: m1 = s1->in[i] = BDY((DPM)BDY(in));
3150: NEWMP(mp); mp->dl = m1->dl; mp->c = m1->c; NEXT(mp) = 0;
3151: MKDP(nv,mp,dp); dp->sugar = mp->dl->td;
3152: m0 = s->sum[m1->pos]; MKDPM(nv,m0,dpm0);
3153: mulobjdpm(CO,(Obj)dp,dpm0,&dpm1);
3154: s1->sum[i] = BDY(dpm1);
3155: }
3156: } else {
1.18 noro 3157: for ( i = 1, in = nd; i <= n; i++, in = NEXT(in) ) {
3158: m0 = BDY((DPM)BDY(in));
3159: NEWDMM(m1); *m1 = *m0; m1->c = (Obj)ONE; NEXT(m1) = 0;
3160: s1->sum[i] = s1->in[i] = m1;
3161: }
1.15 noro 3162: }
3163: s = s1;
3164: level++;
1.22 noro 3165: if ( DP_Print ) printf("level=%d,len=%d\n",level,n);
1.15 noro 3166:
3167: /* create new list */
3168: MKVECT(psv,n+1);
3169: ps = (DPM *)BDY(psv);
3170: for ( i = 1, t = nd; i <= n; i++, t = NEXT(t) ) ps[i] = (DPM)BDY(t);
3171: for ( i = 1, max = 0; i <= n; i++ )
3172: if ( (pos=BDY(ps[i])->pos) > max ) max = pos;
3173: MKVECT(psiv,max+1);
3174: psi = (NODE *)BDY(psiv);
3175: for ( i = n; i >= 1; i-- ) {
3176: pos = BDY(ps[i])->pos;
3177: MKNODE(nd,(long)i,psi[pos]); psi[pos] = nd;
3178: }
3179: m = (DP *)MALLOC((n+1)*sizeof(DP));
1.19 noro 3180: m2 = (DP *)MALLOC((n+1)*sizeof(DP));
1.15 noro 3181: b0 = 0;
3182: for ( i = 1; i <= max; i++ ) {
3183: for ( n1 = psi[i]; n1; n1 = NEXT(n1) ) {
3184: bzero(m,(n+1)*sizeof(DP));
3185: p1 = (long)BDY(n1);
3186: for ( n2 = NEXT(n1); n2; n2 = NEXT(n2) ) {
3187: p2 = (long)BDY(n2);
3188: h = dpm_sp_hm(ps[p1],ps[p2]);
3189: for ( n3 = NEXT(n1); n3 != n2; n3 = NEXT(n3) ) {
3190: p3 = (long)BDY(n3);
3191: if ( m[p3] ) {
3192: if ( dp_redble(h,m[p3]) ) {
3193: h = 0; break;
3194: }
3195: if ( dp_redble(m[p3],h) ) m[p3] = 0;
3196: }
3197: }
3198: if ( h ) m[p2] = h;
3199: }
1.22 noro 3200: if ( lex ) {
3201: // compress m to m2
3202: for ( j = 0, n2 = NEXT(n1); n2; n2 = NEXT(n2) ) {
3203: p2 = (long)BDY(n2);
3204: if ( m[p2] ) m2[j++] = m[p2];
3205: }
3206: qsort(m2,j,sizeof(DP),(int (*)(const void *,const void *))compdp_lex);
3207: for ( k = 0; k < j; k++ ) {
3208: NEWDMM(dmm); dmm->dl = BDY(m2[k])->dl; dmm->pos = p1; dmm->c = (Obj)ONE;
1.15 noro 3209: MKDPM(nv,dmm,dpm);
3210: NEXTNODE(b0,b); BDY(b) = (pointer)dpm;
3211: }
1.22 noro 3212: } else {
3213: for ( n2 = NEXT(n1); n2; n2 = NEXT(n2) ) {
3214: p2 = (long)BDY(n2);
3215: if ( m[p2] ) {
3216: NEWDMM(dmm); dmm->dl = BDY(m[p2])->dl; dmm->pos = p1; dmm->c = (Obj)ONE;
3217: MKDPM(nv,dmm,dpm);
3218: NEXTNODE(b0,b); BDY(b) = (pointer)dpm;
3219: }
3220: }
1.15 noro 3221: }
3222: }
3223: }
3224: if ( !b0 ) {
3225: NEWDMMstack_array(dmmstack_array);
3226: dmmstack_array->len = level;
3227: dmmstack_array->body = (DMMstack *)MALLOC(level*sizeof(DMMstack));
3228: for ( i = level-1, s1 = s; i >= 0; i--, s1 = NEXT(s1) )
3229: dmmstack_array->body[i] = s1;
3230: return dmmstack_array;
3231: } else {
3232: NEXT(b) = 0;
3233: nd = b0;
3234: }
3235: }
1.4 noro 3236: }
1.15 noro 3237:
3238: int sch_count,schlast_count;
1.4 noro 3239:
1.3 noro 3240: int compdmm_schreyer(int n,DMM m1,DMM m2)
3241: {
1.15 noro 3242: int pos1,pos2,t,npos1,npos2;
3243: DMM *in,*sum;
3244: DMMstack s;
3245: static DL d1=0,d2=0;
3246: static int dlen=0;
3247:
3248: sch_count++;
3249: pos1 = m1->pos; pos2 = m2->pos;
3250: if ( pos1 == pos2 ) return (*cmpdl)(n,m1->dl,m2->dl);
3251: if ( n > dlen ) {
3252: NEWDL(d1,n); NEWDL(d2,n); dlen = n;
3253: }
3254: sum = dmm_stack->sum;
3255: _adddl(n,m1->dl,sum[pos1]->dl,d1);
3256: _adddl(n,m2->dl,sum[pos2]->dl,d2);
3257: t = (*cmpdl)(n,d1,d2);
3258: if ( sum[pos1]->pos == sum[pos2]->pos && t == 0 ) {
3259: for ( s = dmm_stack; s; s = NEXT(s) ) {
3260: in = s->in;
3261: npos1 = in[pos1]->pos;
3262: npos2 = in[pos2]->pos;
3263: if ( npos1 == npos2 ) break;
3264: else {
3265: pos1 = npos1;
3266: pos2 = npos2;
3267: }
3268: }
3269: // (pos1,pos2) = the last pair s.t. pos1 != pos2
3270: // simply compare pos1 and pos2
1.3 noro 3271: if ( pos1 < pos2 ) return 1;
3272: else if ( pos1 > pos2 ) return -1;
3273: } else {
1.15 noro 3274: schlast_count++;
3275: pos1 = sum[pos1]->pos;
3276: pos2 = sum[pos2]->pos;
3277: if ( dpm_base_ordtype == 1 ) {
3278: if ( pos1 < pos2 ) return 1;
3279: else if ( pos1 > pos2 ) return -1;
3280: else return t;
3281: } else {
3282: if ( t ) return t;
3283: else if ( pos1 < pos2 ) return 1;
3284: else if ( pos1 > pos2 ) return -1;
3285: else return 0;
3286: }
1.3 noro 3287: }
1.24 noro 3288: /* XXX */
3289: return 0;
1.3 noro 3290: }
3291:
1.20 noro 3292: int compdmm_schreyer_old(int n,DMM m1,DMM m2)
3293: {
3294: int pos1,pos2,t,npos1,npos2;
3295: DMM *in,*sum;
3296: DMMstack s;
3297: static DL d1=0,d2=0;
3298: static int dlen=0;
3299:
3300: sch_count++;
3301: pos1 = m1->pos; pos2 = m2->pos;
3302: if ( pos1 == pos2 ) return (*cmpdl)(n,m1->dl,m2->dl);
3303: if ( n > dlen ) {
3304: NEWDL(d1,n); NEWDL(d2,n); dlen = n;
3305: }
3306: sum = dmm_stack->sum;
3307: _copydl(n,m1->dl,d1);
3308: _copydl(n,m2->dl,d2);
3309: for ( s = dmm_stack; s; s = NEXT(s) ) {
3310: in = s->in;
3311: _addtodl(n,in[pos1]->dl,d1);
3312: _addtodl(n,in[pos2]->dl,d2);
3313: if ( in[pos1]->pos == in[pos2]->pos && _eqdl(n,d1,d2)) {
3314: if ( pos1 < pos2 ) return 1;
3315: else if ( pos1 > pos2 ) return -1;
3316: else return 0;
3317: }
3318: pos1 = in[pos1]->pos;
3319: pos2 = in[pos2]->pos;
3320: if ( pos1 == pos2 ) return (*cmpdl)(n,d1,d2);
3321: }
3322: if ( dpm_base_ordtype == 1 ) {
3323: if ( pos1 < pos2 ) return 1;
3324: else if ( pos1 > pos2 ) return -1;
3325: else return (*cmpdl)(n,d1,d2);
3326: } else {
3327: t = (*cmpdl)(n,d1,d2);
3328: if ( t ) return t;
3329: else if ( pos1 < pos2 ) return 1;
3330: else if ( pos1 > pos2 ) return -1;
3331: else return 0;
3332: }
3333: }
3334:
3335: extern int NaiveSchreyer;
3336:
1.4 noro 3337: int compdmm(int n,DMM m1,DMM m2)
3338: {
1.20 noro 3339: int t,t1;
1.19 noro 3340: int *base_ord;
1.4 noro 3341:
3342: switch ( dpm_ordtype ) {
1.9 noro 3343: case 0: /* TOP ord->pos */
1.4 noro 3344: t = (*cmpdl)(n,m1->dl,m2->dl);
3345: if ( t ) return t;
3346: else if ( m1->pos < m2->pos ) return 1;
3347: else if ( m1->pos > m2->pos ) return -1;
3348: else return 0;
1.9 noro 3349: case 1: /* POT : pos->ord */
1.4 noro 3350: if ( m1->pos < m2->pos ) return 1;
3351: else if ( m1->pos > m2->pos ) return -1;
3352: else return (*cmpdl)(n,m1->dl,m2->dl);
1.9 noro 3353: case 2: /* wPOT: weight->pos->ord */
3354: if ( m1->dl->td > m2->dl->td ) return 1;
3355: else if ( m1->dl->td < m2->dl->td ) return 1;
3356: else if ( m1->pos < m2->pos ) return 1;
3357: else if ( m1->pos > m2->pos ) return -1;
3358: else return (*cmpdl)(n,m1->dl,m2->dl);
3359: case 3: /* Schreyer */
1.20 noro 3360: if ( NaiveSchreyer )
3361: t = compdmm_schreyer_old(n,m1,m2);
3362: else
3363: t = compdmm_schreyer(n,m1,m2);
3364: return t;
1.19 noro 3365: case 4: /* POT with base_ord */
3366: base_ord = dp_current_spec->module_base_ord;
3367: if ( base_ord[m1->pos] < base_ord[m2->pos] ) return 1;
3368: else if ( base_ord[m1->pos] > base_ord[m2->pos] ) return -1;
3369: else return (*cmpdl)(n,m1->dl,m2->dl);
1.4 noro 3370: default:
3371: error("compdmm : invalid dpm_ordtype");
1.24 noro 3372: return 0;
1.4 noro 3373: }
3374: }
1.1 noro 3375:
3376: void adddpm(VL vl,DPM p1,DPM p2,DPM *pr)
3377: {
1.4 noro 3378: int n,s;
1.1 noro 3379: DMM m1,m2,mr=0,mr0;
3380: Obj t;
3381: DL d;
3382:
3383: if ( !p1 )
3384: *pr = p2;
3385: else if ( !p2 )
3386: *pr = p1;
3387: else {
1.4 noro 3388: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
3389: s = compdmm(n,m1,m2);
3390: switch ( s ) {
1.1 noro 3391: case 0:
3392: arf_add(vl,C(m1),C(m2),&t);
3393: if ( t ) {
3394: NEXTDMM(mr0,mr); mr->pos = m1->pos; mr->dl = m1->dl; C(mr) = t;
3395: }
3396: m1 = NEXT(m1); m2 = NEXT(m2); break;
3397: case 1:
3398: NEXTDMM(mr0,mr); mr->pos = m1->pos; mr->dl = m1->dl; C(mr) = C(m1);
3399: m1 = NEXT(m1); break;
3400: case -1:
3401: NEXTDMM(mr0,mr); mr->pos = m2->pos; mr->dl = m2->dl; C(mr) = C(m2);
3402: m2 = NEXT(m2); break;
3403: }
1.4 noro 3404: }
1.1 noro 3405: if ( !mr0 )
3406: if ( m1 )
3407: mr0 = m1;
3408: else if ( m2 )
3409: mr0 = m2;
3410: else {
3411: *pr = 0;
3412: return;
3413: }
3414: else if ( m1 )
3415: NEXT(mr) = m1;
3416: else if ( m2 )
3417: NEXT(mr) = m2;
3418: else
3419: NEXT(mr) = 0;
3420: MKDPM(NV(p1),mr0,*pr);
3421: if ( *pr )
3422: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
3423: }
3424: }
3425:
3426: void subdpm(VL vl,DPM p1,DPM p2,DPM *pr)
3427: {
3428: DPM t;
3429:
3430: if ( !p2 )
3431: *pr = p1;
3432: else {
3433: chsgndpm(p2,&t); adddpm(vl,p1,t,pr);
3434: }
3435: }
3436:
3437: void chsgndpm(DPM p,DPM *pr)
3438: {
3439: DMM m,mr=0,mr0;
3440: Obj r;
3441:
3442: if ( !p )
3443: *pr = 0;
3444: else {
3445: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
3446: NEXTDMM(mr0,mr); arf_chsgn(C(m),&C(mr)); mr->pos = m->pos; mr->dl = m->dl;
3447: }
3448: NEXT(mr) = 0; MKDPM(NV(p),mr0,*pr);
3449: if ( *pr )
3450: (*pr)->sugar = p->sugar;
3451: }
3452: }
3453:
1.8 noro 3454: void mulcmp(Obj c,MP m)
3455: {
3456: MP t;
3457: Obj c1;
3458:
3459: for ( t = m; t; t = NEXT(t) ) {
3460: arf_mul(CO,c,C(t),&c1); C(t) = c1;
3461: }
3462: }
3463:
3464: void mulcdmm(Obj c,DMM m)
3465: {
3466: DMM t;
3467: Obj c1;
3468:
3469: for ( t = m; t; t = NEXT(t) ) {
3470: arf_mul(CO,c,C(t),&c1); C(t) = c1;
3471: }
3472: }
3473:
1.1 noro 3474: void mulcdpm(VL vl,Obj c,DPM p,DPM *pr)
3475: {
3476: DMM m,mr=0,mr0;
3477:
3478: if ( !p || !c )
3479: *pr = 0;
3480: else if ( NUM(c) && UNIQ((Q)c) )
3481: *pr = p;
3482: else if ( NUM(c) && MUNIQ((Q)c) )
3483: chsgndpm(p,pr);
3484: else {
3485: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
3486: NEXTDMM(mr0,mr);
3487: arf_mul(vl,C(m),c,&C(mr));
3488: mr->pos = m->pos;
3489: mr->dl = m->dl;
3490: }
3491: NEXT(mr) = 0; MKDPM(NV(p),mr0,*pr);
3492: if ( *pr )
3493: (*pr)->sugar = p->sugar;
3494: }
3495: }
3496:
3497: void comm_mulmpdpm(VL vl,MP m0,DPM p,DPM *pr)
3498: {
3499: DMM m,mr=0,mr0;
3500: DL d;
3501: Obj c;
3502: int n;
3503:
3504: if ( !p )
3505: *pr = 0;
3506: else {
3507: n = NV(p);
3508: d = m0->dl;
3509: c = C(m0);
3510: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
3511: NEXTDMM(mr0,mr);
3512: arf_mul(vl,C(m),c,&C(mr));
3513: mr->pos = m->pos;
3514: adddl(n,m->dl,d,&mr->dl);
3515: }
3516: NEXT(mr) = 0; MKDPM(NV(p),mr0,*pr);
3517: if ( *pr )
3518: (*pr)->sugar = p->sugar;
3519: }
3520: }
3521:
3522: void weyl_mulmpdpm(VL vl,MP m0,DPM p,DPM *pr)
3523: {
3524: DPM r,t,t1;
3525: DMM m;
3526: DL d0;
3527: int n,n2,l,i,j,tlen;
3528: struct oMP mp;
3529: static DMM *w,*psum;
3530: static struct cdl *tab;
3531: static int wlen;
3532: static int rtlen;
3533:
3534: if ( !p )
3535: *pr = 0;
3536: else {
3537: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
3538: if ( l > wlen ) {
3539: if ( w ) GCFREE(w);
3540: w = (DMM *)MALLOC(l*sizeof(DMM));
3541: wlen = l;
3542: }
3543: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
3544: w[i] = m;
3545:
3546: n = NV(p); n2 = n>>1;
3547: d0 = m0->dl;
3548: for ( i = 0, tlen = 1; i < n2; i++ )
3549: tlen *= d0->d[n2+i]+1;
3550: if ( tlen > rtlen ) {
3551: if ( tab ) GCFREE(tab);
3552: if ( psum ) GCFREE(psum);
3553: rtlen = tlen;
3554: tab = (struct cdl *)MALLOC(rtlen*sizeof(struct cdl));
3555: psum = (DMM *)MALLOC(rtlen*sizeof(DMM));
3556: }
3557: bzero(psum,tlen*sizeof(DMM));
3558: for ( i = l-1; i >= 0; i-- ) {
3559: bzero(tab,tlen*sizeof(struct cdl));
3560: mp.dl = w[i]->dl; mp.c = C(w[i]); mp.next = 0;
3561: weyl_mulmm(vl,m0,&mp,n,tab,tlen);
3562: for ( j = 0; j < tlen; j++ ) {
3563: if ( tab[j].c ) {
3564: NEWDMM(m); m->dl = tab[j].d; m->pos = w[i]->pos; C(m) = (Obj)tab[j].c; NEXT(m) = psum[j];
3565: psum[j] = m;
3566: }
3567: }
3568: }
3569: for ( j = tlen-1, r = 0; j >= 0; j-- )
3570: if ( psum[j] ) {
3571: MKDPM(n,psum[j],t); adddpm(vl,r,t,&t1); r = t1;
3572: }
3573: if ( r )
3574: r->sugar = p->sugar + m0->dl->td;
3575: *pr = r;
3576: }
3577: }
3578:
3579: void mulobjdpm(VL vl,Obj p1,DPM p2,DPM *pr)
3580: {
3581: MP m;
3582: DPM s,t,u;
3583:
3584: if ( !p1 || !p2 )
3585: *pr = 0;
3586: else if ( OID(p1) != O_DP )
3587: mulcdpm(vl,p1,p2,pr);
3588: else {
3589: s = 0;
3590: for ( m = BDY((DP)p1); m; m = NEXT(m) ) {
3591: if ( do_weyl )
3592: weyl_mulmpdpm(vl,m,p2,&t);
3593: else
3594: comm_mulmpdpm(vl,m,p2,&t);
3595: adddpm(vl,s,t,&u); s = u;
3596: }
3597: *pr = s;
3598: }
3599: }
3600:
3601: int compdpm(VL vl,DPM p1,DPM p2)
3602: {
3603: int n,t;
3604: DMM m1,m2;
3605:
3606: if ( !p1 )
3607: return p2 ? -1 : 0;
3608: else if ( !p2 )
3609: return 1;
3610: else if ( NV(p1) != NV(p2) ) {
3611: error("compdpm : size mismatch");
3612: return 0; /* XXX */
3613: } else {
3614: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
3615: m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
3616: if ( (t = compdmm(n,m1,m2)) ||
3617: (t = arf_comp(vl,C(m1),C(m2)) ) )
3618: return t;
3619: if ( m1 )
3620: return 1;
3621: else if ( m2 )
3622: return -1;
3623: else
3624: return 0;
3625: }
3626: }
3627:
1.11 noro 3628: void dpm_removecont2(DPM p1,DPM p2,DPM *r1p,DPM *r2p,Z *contp);
3629:
1.5 noro 3630: // p = ...+c*<<0,...0:pos>>+...
3631: DPM dpm_eliminate_term(DPM a,DPM p,Obj c,int pos)
3632: {
3633: MP d0,d;
1.10 noro 3634: DL dl;
1.5 noro 3635: DMM m;
3636: DP f;
1.11 noro 3637: DPM a1,p1,r,r1,dmy;
3638: Z dmyz;
1.5 noro 3639:
3640: if ( !a ) return 0;
3641: d0 = 0;
3642: for ( m = BDY(a); m; m = NEXT(m) )
3643: if ( m->pos == pos ) {
1.10 noro 3644: NEXTMP(d0,d);
3645: arf_chsgn(m->c,&d->c);
1.19 noro 3646: if ( !dp_current_spec || !dp_current_spec->module_top_weight )
1.11 noro 3647: d->dl = m->dl;
1.10 noro 3648: else {
3649: NEWDL(dl,NV(a));
3650: _copydl(NV(a),m->dl,dl);
3651: dl->td -= dp_current_spec->module_top_weight[pos-1];
3652: d->dl = dl;
3653: }
1.5 noro 3654: }
3655: if ( d0 ) {
1.6 noro 3656: NEXT(d) = 0; MKDP(NV(a),d0,f);
3657: mulcdpm(CO,c,a,&a1);
3658: mulobjdpm(CO,(Obj)f,p,&p1);
3659: adddpm(CO,a1,p1,&r);
1.11 noro 3660: dpm_removecont2(0,r,&dmy,&r1,&dmyz);
3661: return r1;
1.5 noro 3662: } else
1.6 noro 3663: return a;
1.5 noro 3664: }
3665:
3666: // <<...:i>> -> <<...:tab[i]>>
3667: DPM dpm_compress(DPM p,int *tab)
3668: {
3669: DMM m,mr0,mr;
3670: DPM t;
3671:
3672: if ( !p ) return 0;
3673: else {
3674: for ( m = BDY(p), mr0 = 0; m; m = NEXT(m) ) {
3675: NEXTDMM(mr0,mr);
1.6 noro 3676: mr->dl = m->dl; mr->c = m->c; mr->pos = tab[m->pos];
1.5 noro 3677: if ( mr->pos <= 0 )
3678: error("dpm_compress : invalid position");
3679: }
3680: NEXT(mr) = 0;
3681: MKDPM(p->nv,mr0,t); t->sugar = p->sugar;
3682: return t;
3683: }
3684: }
3685:
3686: // input : s, s = syz(m) output simplified s, m
1.21 noro 3687: // assuming the term order is POT
1.10 noro 3688: void dpm_simplify_syz(LIST s,LIST m,LIST *s1,LIST *m1,LIST *w1)
1.5 noro 3689: {
1.9 noro 3690: int lm,ls,i,j,k,pos,nv;
1.5 noro 3691: DPM *am,*as;
3692: DPM p;
3693: DMM d;
3694: Obj c;
1.10 noro 3695: Z q;
3696: int *tab,*dd,*new_w;
1.5 noro 3697: NODE t,t1;
3698:
3699: lm = length(BDY(m));
3700: am = (DPM *)MALLOC((lm+1)*sizeof(DPM));
3701: ls = length(BDY(s));
3702: as = (DPM *)MALLOC(ls*sizeof(DPM));
3703: for ( i = 1, t = BDY(m); i <= lm; i++, t = NEXT(t) ) am[i] = (DPM)BDY(t);
1.6 noro 3704: for ( i = 0, t = BDY(s); i < ls; i++, t = NEXT(t) ) as[i] = (DPM)BDY(t);
1.5 noro 3705:
3706: for ( i = 0; i < ls; i++ ) {
3707: p = as[i];
1.6 noro 3708: if ( p == 0 ) continue;
1.9 noro 3709: nv = NV(p);
1.21 noro 3710: for ( d = BDY(p); d; ) {
1.9 noro 3711: dd = d->dl->d;
3712: for ( k = 0; k < nv; k++ ) if ( dd[k] ) break;
3713: if ( k == nv ) break;
1.21 noro 3714: pos = d->pos;
3715: while ( d && d->pos == pos ) d = NEXT(d);
1.9 noro 3716: }
1.5 noro 3717: if ( d ) {
3718: c = d->c; pos = d->pos;
3719: for ( j = 0; j < ls; j++ )
1.6 noro 3720: if ( j != i ) {
1.5 noro 3721: as[j] = dpm_eliminate_term(as[j],p,c,pos);
1.6 noro 3722: }
1.5 noro 3723: // remove m[i]
3724: am[pos] = 0;
3725: as[i] = 0;
3726: }
3727: }
3728: // compress s
3729: // create index table from am[]
3730: // (0 0 * 0 * ...) -> (0 0 1 0 2 ... ) which means 2->1, 4->2, ...
3731: tab = (int *)MALLOC((lm+1)*sizeof(int));
3732: for ( j = 0, i = 1; i <= lm; i++ ) {
3733: if ( am[i] ) { j++; tab[i] = j; }
1.10 noro 3734: else { tab[i] = 0; }
1.5 noro 3735: }
3736: t = 0;
3737: for ( i = ls-1; i >= 0; i-- )
3738: if ( as[i] ) {
3739: p = dpm_compress(as[i],tab);
3740: MKNODE(t1,(pointer)p,t); t = t1;
3741: }
3742: MKLIST(*s1,t);
1.10 noro 3743:
1.19 noro 3744: if ( dp_current_spec && dp_current_spec->module_top_weight ) {
1.10 noro 3745: new_w = (int *)MALLOC(j*sizeof(int));
3746: for ( j = 0, i = 1; i <= lm; i++ )
3747: if ( tab[i] ) { new_w[j++] = dp_current_spec->module_top_weight[i-1]; }
3748: t = 0;
3749: for ( i = j-1; i >= 0; i-- ) {
3750: STOZ(new_w[i],q);
3751: MKNODE(t1,q,t); t = t1;
3752: }
3753: } else
3754: t = 0;
3755: MKLIST(*w1,t);
3756:
1.5 noro 3757: t = 0;
3758: for ( i = lm; i >= 1; i-- )
3759: if ( am[i] ) {
3760: MKNODE(t1,(pointer)am[i],t); t = t1;
3761: }
3762: MKLIST(*m1,t);
3763: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>