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