![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to dist.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine|
Return to dist.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.8 noro 1: /*
2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
3: * All rights reserved.
4: *
5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
6: * non-exclusive and royalty-free license to use, copy, modify and
7: * redistribute, solely for non-commercial and non-profit purposes, the
8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
9: * conditions of this Agreement. For the avoidance of doubt, you acquire
10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
11: * third party developer retains all rights, including but not limited to
12: * copyrights, in and to the SOFTWARE.
13: *
14: * (1) FLL does not grant you a license in any way for commercial
15: * purposes. You may use the SOFTWARE only for non-commercial and
16: * non-profit purposes only, such as academic, research and internal
17: * business use.
18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
19: * international copyright treaties. If you make copies of the SOFTWARE,
20: * with or without modification, as permitted hereunder, you shall affix
21: * to all such copies of the SOFTWARE the above copyright notice.
22: * (3) An explicit reference to this SOFTWARE and its copyright owner
23: * shall be made on your publication or presentation in any form of the
24: * results obtained by use of the SOFTWARE.
25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.9 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.8 noro 27: * for such modification or the source code of the modified part of the
28: * SOFTWARE.
29: *
30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
47: *
1.28 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.27 2004/02/03 23:31:57 noro Exp $
1.8 noro 49: */
1.1 noro 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
1.12 noro 62: #define ORD_WEYL_ELIM 10
1.13 noro 63: #define ORD_HOMO_WW_DRL 11
1.21 noro 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();
1.1 noro 68:
69: int (*cmpdl)()=cmpdl_revgradlex;
70: int (*primitive_cmpdl[3])() = {cmpdl_revgradlex,cmpdl_gradlex,cmpdl_lex};
71:
1.2 noro 72: int do_weyl;
73:
1.1 noro 74: int dp_nelim,dp_fcoeffs;
1.27 noro 75: struct order_spec *dp_current_spec;
1.1 noro 76: int *dp_dl_work;
77:
1.24 noro 78: void comm_muld_trunc(VL vl,DP p1,DP p2,DL dl,DP *pr);
79: void comm_quod(VL vl,DP p1,DP p2,DP *pr);
80: void muldm_trunc(VL vl,DP p,MP m0,DL dl,DP *pr);
81: void muldc_trunc(VL vl,DP p,P c,DL dl,DP *pr);
82:
1.22 noro 83: int has_sfcoef(DP f)
1.1 noro 84: {
85: MP t;
86:
87: if ( !f )
88: return 0;
89: for ( t = BDY(f); t; t = NEXT(t) )
1.22 noro 90: if ( has_sfcoef_p(t->c) )
1.1 noro 91: break;
92: return t ? 1 : 0;
93: }
94:
1.22 noro 95: int has_sfcoef_p(P f)
1.1 noro 96: {
97: DCP dc;
98:
99: if ( !f )
100: return 0;
101: else if ( NUM(f) )
1.22 noro 102: return (NID((Num)f) == N_GFS) ? 1 : 0;
1.1 noro 103: else {
104: for ( dc = DC(f); dc; dc = NEXT(dc) )
1.22 noro 105: if ( has_sfcoef_p(COEF(dc)) )
1.1 noro 106: return 1;
107: return 0;
108: }
109: }
110:
1.19 noro 111: void initd(struct order_spec *spec)
1.1 noro 112: {
113: switch ( spec->id ) {
1.28 ! noro 114: case 3:
! 115: cmpdl = cmpdl_composite;
! 116: dp_dl_work = (int *)MALLOC_ATOMIC(spec->nv*sizeof(int));
! 117: break;
1.1 noro 118: case 2:
119: cmpdl = cmpdl_matrix;
120: dp_dl_work = (int *)MALLOC_ATOMIC(spec->nv*sizeof(int));
121: break;
122: case 1:
123: cmpdl = cmpdl_order_pair;
124: break;
125: default:
126: switch ( spec->ord.simple ) {
127: case ORD_REVGRADLEX:
128: cmpdl = cmpdl_revgradlex; break;
129: case ORD_GRADLEX:
130: cmpdl = cmpdl_gradlex; break;
131: case ORD_BREVGRADLEX:
132: cmpdl = cmpdl_brevgradlex; break;
133: case ORD_BGRADLEX:
134: cmpdl = cmpdl_bgradlex; break;
135: case ORD_BLEX:
136: cmpdl = cmpdl_blex; break;
137: case ORD_BREVREV:
138: cmpdl = cmpdl_brevrev; break;
139: case ORD_BGRADREV:
140: cmpdl = cmpdl_bgradrev; break;
141: case ORD_BLEXREV:
142: cmpdl = cmpdl_blexrev; break;
143: case ORD_ELIM:
144: cmpdl = cmpdl_elim; break;
1.12 noro 145: case ORD_WEYL_ELIM:
146: cmpdl = cmpdl_weyl_elim; break;
1.13 noro 147: case ORD_HOMO_WW_DRL:
148: cmpdl = cmpdl_homo_ww_drl; break;
1.21 noro 149: case ORD_DRL_ZIGZAG:
150: cmpdl = cmpdl_drl_zigzag; break;
151: case ORD_HOMO_WW_DRL_ZIGZAG:
152: cmpdl = cmpdl_homo_ww_drl_zigzag; break;
1.1 noro 153: case ORD_LEX: default:
154: cmpdl = cmpdl_lex; break;
155: }
156: break;
157: }
1.27 noro 158: dp_current_spec = spec;
1.1 noro 159: }
160:
1.19 noro 161: void ptod(VL vl,VL dvl,P p,DP *pr)
1.1 noro 162: {
163: int isconst = 0;
1.16 noro 164: int n,i,j,k;
1.1 noro 165: VL tvl;
166: V v;
167: DL d;
168: MP m;
169: DCP dc;
1.16 noro 170: DCP *w;
1.1 noro 171: DP r,s,t,u;
172: P x,c;
173:
174: if ( !p )
175: *pr = 0;
176: else {
177: for ( n = 0, tvl = dvl; tvl; tvl = NEXT(tvl), n++ );
178: if ( NUM(p) ) {
179: NEWDL(d,n);
180: NEWMP(m); m->dl = d; C(m) = p; NEXT(m) = 0; MKDP(n,m,*pr); (*pr)->sugar = 0;
181: } else {
182: for ( i = 0, tvl = dvl, v = VR(p);
183: tvl && tvl->v != v; tvl = NEXT(tvl), i++ );
184: if ( !tvl ) {
1.16 noro 185: for ( dc = DC(p), k = 0; dc; dc = NEXT(dc), k++ );
186: w = (DCP *)ALLOCA(k*sizeof(DCP));
187: for ( dc = DC(p), j = 0; j < k; dc = NEXT(dc), j++ )
188: w[j] = dc;
189:
190: for ( j = k-1, s = 0, MKV(v,x); j >= 0; j-- ) {
191: ptod(vl,dvl,COEF(w[j]),&t); pwrp(vl,x,DEG(w[j]),&c);
1.1 noro 192: muldc(vl,t,c,&r); addd(vl,r,s,&t); s = t;
193: }
194: *pr = s;
195: } else {
1.16 noro 196: for ( dc = DC(p), k = 0; dc; dc = NEXT(dc), k++ );
197: w = (DCP *)ALLOCA(k*sizeof(DCP));
198: for ( dc = DC(p), j = 0; j < k; dc = NEXT(dc), j++ )
199: w[j] = dc;
200:
201: for ( j = k-1, s = 0; j >= 0; j-- ) {
202: ptod(vl,dvl,COEF(w[j]),&t);
1.20 noro 203: NEWDL(d,n); d->d[i] = QTOS(DEG(w[j]));
204: d->td = MUL_WEIGHT(d->d[i],i);
1.1 noro 205: NEWMP(m); m->dl = d; C(m) = (P)ONE; NEXT(m) = 0; MKDP(n,m,u); u->sugar = d->td;
1.2 noro 206: comm_muld(vl,t,u,&r); addd(vl,r,s,&t); s = t;
1.1 noro 207: }
208: *pr = s;
209: }
210: }
211: }
1.17 noro 212: #if 0
1.22 noro 213: if ( !dp_fcoeffs && has_sfcoef(*pr) )
214: dp_fcoeffs = N_GFS;
1.17 noro 215: #endif
1.1 noro 216: }
217:
1.19 noro 218: void dtop(VL vl,VL dvl,DP p,P *pr)
1.1 noro 219: {
1.16 noro 220: int n,i,j,k;
1.1 noro 221: DL d;
222: MP m;
1.16 noro 223: MP *a;
1.1 noro 224: P r,s,t,u,w;
225: Q q;
226: VL tvl;
227:
228: if ( !p )
229: *pr = 0;
230: else {
1.16 noro 231: for ( k = 0, m = BDY(p); m; m = NEXT(m), k++ );
232: a = (MP *)ALLOCA(k*sizeof(MP));
233: for ( j = 0, m = BDY(p); j < k; m = NEXT(m), j++ )
234: a[j] = m;
235:
236: for ( n = p->nv, j = k-1, s = 0; j >= 0; j-- ) {
237: m = a[j];
1.1 noro 238: t = C(m);
239: if ( NUM(t) && NID((Num)t) == N_M ) {
240: mptop(t,&u); t = u;
241: }
242: for ( i = 0, d = m->dl, tvl = dvl;
243: i < n; tvl = NEXT(tvl), i++ ) {
244: MKV(tvl->v,r); STOQ(d->d[i],q); pwrp(vl,r,q,&u);
245: mulp(vl,t,u,&w); t = w;
246: }
247: addp(vl,s,t,&u); s = u;
248: }
249: *pr = s;
250: }
251: }
252:
1.19 noro 253: void nodetod(NODE node,DP *dp)
1.1 noro 254: {
255: NODE t;
256: int len,i,td;
257: Q e;
258: DL d;
259: MP m;
260: DP u;
261:
262: for ( t = node, len = 0; t; t = NEXT(t), len++ );
263: NEWDL(d,len);
264: for ( t = node, i = 0, td = 0; i < len; t = NEXT(t), i++ ) {
265: e = (Q)BDY(t);
266: if ( !e )
267: d->d[i] = 0;
268: else if ( !NUM(e) || !RATN(e) || !INT(e) )
269: error("nodetod : invalid input");
270: else {
1.20 noro 271: d->d[i] = QTOS((Q)e); td += MUL_WEIGHT(d->d[i],i);
1.1 noro 272: }
273: }
274: d->td = td;
275: NEWMP(m); m->dl = d; C(m) = (P)ONE; NEXT(m) = 0;
276: MKDP(len,m,u); u->sugar = td; *dp = u;
277: }
278:
1.19 noro 279: int sugard(MP m)
1.1 noro 280: {
281: int s;
282:
283: for ( s = 0; m; m = NEXT(m) )
284: s = MAX(s,m->dl->td);
285: return s;
286: }
287:
1.19 noro 288: void addd(VL vl,DP p1,DP p2,DP *pr)
1.1 noro 289: {
290: int n;
291: MP m1,m2,mr,mr0;
292: P t;
293:
294: if ( !p1 )
295: *pr = p2;
296: else if ( !p2 )
297: *pr = p1;
298: else {
299: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
300: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
301: case 0:
302: addp(vl,C(m1),C(m2),&t);
303: if ( t ) {
304: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = t;
305: }
306: m1 = NEXT(m1); m2 = NEXT(m2); break;
307: case 1:
308: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
309: m1 = NEXT(m1); break;
310: case -1:
311: NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
312: m2 = NEXT(m2); break;
313: }
314: if ( !mr0 )
315: if ( m1 )
316: mr0 = m1;
317: else if ( m2 )
318: mr0 = m2;
319: else {
320: *pr = 0;
321: return;
322: }
323: else if ( m1 )
324: NEXT(mr) = m1;
325: else if ( m2 )
326: NEXT(mr) = m2;
327: else
328: NEXT(mr) = 0;
329: MKDP(NV(p1),mr0,*pr);
330: if ( *pr )
331: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
332: }
333: }
334:
335: /* for F4 symbolic reduction */
336:
1.19 noro 337: void symb_addd(DP p1,DP p2,DP *pr)
1.1 noro 338: {
339: int n;
340: MP m1,m2,mr,mr0;
341:
342: if ( !p1 )
343: *pr = p2;
344: else if ( !p2 )
345: *pr = p1;
346: else {
347: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
348: NEXTMP(mr0,mr); C(mr) = (P)ONE;
349: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
350: case 0:
351: mr->dl = m1->dl;
352: m1 = NEXT(m1); m2 = NEXT(m2); break;
353: case 1:
354: mr->dl = m1->dl;
355: m1 = NEXT(m1); break;
356: case -1:
357: mr->dl = m2->dl;
358: m2 = NEXT(m2); break;
359: }
360: }
361: if ( !mr0 )
362: if ( m1 )
363: mr0 = m1;
364: else if ( m2 )
365: mr0 = m2;
366: else {
367: *pr = 0;
368: return;
369: }
370: else if ( m1 )
371: NEXT(mr) = m1;
372: else if ( m2 )
373: NEXT(mr) = m2;
374: else
375: NEXT(mr) = 0;
376: MKDP(NV(p1),mr0,*pr);
377: if ( *pr )
378: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
1.3 noro 379: }
380: }
381:
382: /*
383: * destructive merge of two list
384: *
385: * p1, p2 : list of DL
386: * return : a merged list
387: */
388:
1.19 noro 389: NODE symb_merge(NODE m1,NODE m2,int n)
1.3 noro 390: {
391: NODE top,prev,cur,m,t;
1.25 noro 392: int c,i;
393: DL d1,d2;
1.3 noro 394:
395: if ( !m1 )
396: return m2;
397: else if ( !m2 )
398: return m1;
399: else {
400: switch ( (*cmpdl)(n,(DL)BDY(m1),(DL)BDY(m2)) ) {
401: case 0:
402: top = m1; m = NEXT(m2);
403: break;
404: case 1:
405: top = m1; m = m2;
406: break;
407: case -1:
408: top = m2; m = m1;
409: break;
410: }
411: prev = top; cur = NEXT(top);
412: /* BDY(prev) > BDY(m) always holds */
413: while ( cur && m ) {
1.25 noro 414: d1 = (DL)BDY(cur);
415: d2 = (DL)BDY(m);
1.26 noro 416: #if 1
417: switch ( (*cmpdl)(n,(DL)BDY(cur),(DL)BDY(m)) ) {
418: #else
419: /* XXX only valid for DRL */
1.25 noro 420: if ( d1->td > d2->td )
421: c = 1;
422: else if ( d1->td < d2->td )
423: c = -1;
424: else {
425: for ( i = n-1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
426: if ( i < 0 )
427: c = 0;
428: else if ( d1->d[i] < d2->d[i] )
429: c = 1;
430: else
431: c = -1;
432: }
433: switch ( c ) {
434: #endif
1.3 noro 435: case 0:
436: m = NEXT(m);
437: prev = cur; cur = NEXT(cur);
438: break;
439: case 1:
440: t = NEXT(cur); NEXT(cur) = m; m = t;
441: prev = cur; cur = NEXT(cur);
442: break;
443: case -1:
444: NEXT(prev) = m; m = cur;
445: prev = NEXT(prev); cur = NEXT(prev);
446: break;
1.18 noro 447: }
448: }
449: if ( !cur )
450: NEXT(prev) = m;
1.23 noro 451: return top;
452: }
453: }
454:
455: void _adddl(int n,DL d1,DL d2,DL d3)
456: {
457: int i;
458:
459: d3->td = d1->td+d2->td;
460: for ( i = 0; i < n; i++ )
461: d3->d[i] = d1->d[i]+d2->d[i];
462: }
463:
464: /* m1 <- m1 U dl*f, destructive */
465:
466: NODE mul_dllist(DL dl,DP f);
467:
468: NODE symb_mul_merge(NODE m1,DL dl,DP f,int n)
469: {
470: NODE top,prev,cur,n1;
471: DP g;
472: DL t,s;
473: MP m;
474:
475: if ( !m1 )
476: return mul_dllist(dl,f);
477: else if ( !f )
478: return m1;
479: else {
480: m = BDY(f);
481: NEWDL_NOINIT(t,n);
482: _adddl(n,m->dl,dl,t);
483: top = m1; prev = 0; cur = m1;
484: while ( m ) {
485: switch ( (*cmpdl)(n,(DL)BDY(cur),t) ) {
486: case 0:
487: prev = cur; cur = NEXT(cur);
488: if ( !cur ) {
489: MKDP(n,m,g);
490: NEXT(prev) = mul_dllist(dl,g);
491: return;
492: }
493: m = NEXT(m);
494: if ( m ) _adddl(n,m->dl,dl,t);
495: break;
496: case 1:
497: prev = cur; cur = NEXT(cur);
498: if ( !cur ) {
499: MKDP(n,m,g);
500: NEXT(prev) = mul_dllist(dl,g);
501: return;
502: }
503: break;
504: case -1:
505: NEWDL_NOINIT(s,n);
506: s->td = t->td;
507: bcopy(t->d,s->d,n*sizeof(int));
508: NEWNODE(n1);
509: n1->body = (pointer)s;
510: NEXT(n1) = cur;
511: if ( !prev ) {
512: top = n1; cur = n1;
513: } else {
514: NEXT(prev) = n1; prev = n1;
515: }
516: m = NEXT(m);
517: if ( m ) _adddl(n,m->dl,dl,t);
518: break;
519: }
520: }
1.18 noro 521: return top;
522: }
523: }
524:
1.19 noro 525: DLBUCKET symb_merge_bucket(DLBUCKET m1,DLBUCKET m2,int n)
1.18 noro 526: {
527: DLBUCKET top,prev,cur,m,t;
528:
529: if ( !m1 )
530: return m2;
531: else if ( !m2 )
532: return m1;
533: else {
534: if ( m1->td == m2->td ) {
535: top = m1;
536: BDY(top) = symb_merge(BDY(top),BDY(m2),n);
537: m = NEXT(m2);
538: } else if ( m1->td > m2->td ) {
539: top = m1; m = m2;
540: } else {
541: top = m2; m = m1;
542: }
543: prev = top; cur = NEXT(top);
544: /* prev->td > m->td always holds */
545: while ( cur && m ) {
546: if ( cur->td == m->td ) {
547: BDY(cur) = symb_merge(BDY(cur),BDY(m),n);
548: m = NEXT(m);
549: prev = cur; cur = NEXT(cur);
550: } else if ( cur->td > m->td ) {
551: t = NEXT(cur); NEXT(cur) = m; m = t;
552: prev = cur; cur = NEXT(cur);
553: } else {
554: NEXT(prev) = m; m = cur;
555: prev = NEXT(prev); cur = NEXT(prev);
1.3 noro 556: }
557: }
558: if ( !cur )
559: NEXT(prev) = m;
560: return top;
1.1 noro 561: }
562: }
563:
1.19 noro 564: void subd(VL vl,DP p1,DP p2,DP *pr)
1.1 noro 565: {
566: DP t;
567:
568: if ( !p2 )
569: *pr = p1;
570: else {
571: chsgnd(p2,&t); addd(vl,p1,t,pr);
572: }
573: }
574:
1.19 noro 575: void chsgnd(DP p,DP *pr)
1.1 noro 576: {
577: MP m,mr,mr0;
578:
579: if ( !p )
580: *pr = 0;
581: else {
582: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
583: NEXTMP(mr0,mr); chsgnp(C(m),&C(mr)); mr->dl = m->dl;
584: }
585: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
586: if ( *pr )
587: (*pr)->sugar = p->sugar;
588: }
589: }
590:
1.19 noro 591: void muld(VL vl,DP p1,DP p2,DP *pr)
1.1 noro 592: {
1.2 noro 593: if ( ! do_weyl )
594: comm_muld(vl,p1,p2,pr);
595: else
596: weyl_muld(vl,p1,p2,pr);
597: }
598:
1.19 noro 599: void comm_muld(VL vl,DP p1,DP p2,DP *pr)
1.2 noro 600: {
1.1 noro 601: MP m;
602: DP s,t,u;
1.5 noro 603: int i,l,l1;
604: static MP *w;
605: static int wlen;
1.1 noro 606:
607: if ( !p1 || !p2 )
608: *pr = 0;
609: else if ( OID(p1) <= O_P )
610: muldc(vl,p2,(P)p1,pr);
611: else if ( OID(p2) <= O_P )
612: muldc(vl,p1,(P)p2,pr);
613: else {
1.5 noro 614: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
1.4 noro 615: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
1.5 noro 616: if ( l1 < l ) {
617: t = p1; p1 = p2; p2 = t;
618: l = l1;
619: }
620: if ( l > wlen ) {
621: if ( w ) GC_free(w);
622: w = (MP *)MALLOC(l*sizeof(MP));
623: wlen = l;
624: }
1.4 noro 625: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
626: w[i] = m;
627: for ( s = 0, i = l-1; i >= 0; i-- ) {
628: muldm(vl,p1,w[i],&t); addd(vl,s,t,&u); s = u;
1.1 noro 629: }
1.5 noro 630: bzero(w,l*sizeof(MP));
1.1 noro 631: *pr = s;
632: }
633: }
634:
1.24 noro 635: /* discard terms which is not a multiple of dl */
636:
637: void comm_muld_trunc(VL vl,DP p1,DP p2,DL dl,DP *pr)
638: {
639: MP m;
640: DP s,t,u;
641: int i,l,l1;
642: static MP *w;
643: static int wlen;
644:
645: if ( !p1 || !p2 )
646: *pr = 0;
647: else if ( OID(p1) <= O_P )
648: muldc_trunc(vl,p2,(P)p1,dl,pr);
649: else if ( OID(p2) <= O_P )
650: muldc_trunc(vl,p1,(P)p2,dl,pr);
651: else {
652: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
653: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
654: if ( l1 < l ) {
655: t = p1; p1 = p2; p2 = t;
656: l = l1;
657: }
658: if ( l > wlen ) {
659: if ( w ) GC_free(w);
660: w = (MP *)MALLOC(l*sizeof(MP));
661: wlen = l;
662: }
663: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
664: w[i] = m;
665: for ( s = 0, i = l-1; i >= 0; i-- ) {
666: muldm_trunc(vl,p1,w[i],dl,&t); addd(vl,s,t,&u); s = u;
667: }
668: bzero(w,l*sizeof(MP));
669: *pr = s;
670: }
671: }
672:
673: void comm_quod(VL vl,DP p1,DP p2,DP *pr)
674: {
675: MP m,m0;
676: DP s,t;
677: int i,n,sugar;
678: DL d1,d2,d;
679: Q a,b;
680:
681: if ( !p2 )
682: error("comm_quod : invalid input");
683: if ( !p1 )
684: *pr = 0;
685: else {
686: n = NV(p1);
687: d2 = BDY(p2)->dl;
688: m0 = 0;
689: sugar = p1->sugar;
690: while ( p1 ) {
691: d1 = BDY(p1)->dl;
692: NEWDL(d,n);
693: d->td = d1->td - d2->td;
694: for ( i = 0; i < n; i++ )
695: d->d[i] = d1->d[i]-d2->d[i];
696: NEXTMP(m0,m);
697: m->dl = d;
698: divq((Q)BDY(p1)->c,(Q)BDY(p2)->c,&a); chsgnq(a,&b);
699: C(m) = (P)b;
700: muldm_trunc(vl,p2,m,d2,&t);
701: addd(vl,p1,t,&s); p1 = s;
702: C(m) = (P)a;
703: }
704: if ( m0 ) {
705: NEXT(m) = 0; MKDP(n,m0,*pr);
706: } else
707: *pr = 0;
708: /* XXX */
709: if ( *pr )
710: (*pr)->sugar = sugar - d2->td;
711: }
712: }
713:
1.19 noro 714: void muldm(VL vl,DP p,MP m0,DP *pr)
1.1 noro 715: {
716: MP m,mr,mr0;
717: P c;
718: DL d;
719: int n;
720:
721: if ( !p )
722: *pr = 0;
723: else {
724: for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
725: m; m = NEXT(m) ) {
726: NEXTMP(mr0,mr);
727: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
728: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
729: else
730: mulp(vl,C(m),c,&C(mr));
731: adddl(n,m->dl,d,&mr->dl);
732: }
733: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
734: if ( *pr )
735: (*pr)->sugar = p->sugar + m0->dl->td;
1.2 noro 736: }
737: }
738:
1.24 noro 739: void muldm_trunc(VL vl,DP p,MP m0,DL dl,DP *pr)
740: {
741: MP m,mr,mr0;
742: P c;
743: DL d,tdl;
744: int n,i;
745:
746: if ( !p )
747: *pr = 0;
748: else {
749: n = NV(p);
750: NEWDL(tdl,n);
751: for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl;
752: m; m = NEXT(m) ) {
753: _adddl(n,m->dl,d,tdl);
754: for ( i = 0; i < n; i++ )
755: if ( tdl->d[i] < dl->d[i] )
756: break;
757: if ( i < n )
758: continue;
759: NEXTMP(mr0,mr);
760: mr->dl = tdl;
761: NEWDL(tdl,n);
762: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
763: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
764: else
765: mulp(vl,C(m),c,&C(mr));
766: }
767: if ( mr0 ) {
768: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
769: } else
770: *pr = 0;
771: if ( *pr )
772: (*pr)->sugar = p->sugar + m0->dl->td;
773: }
774: }
775:
1.19 noro 776: void weyl_muld(VL vl,DP p1,DP p2,DP *pr)
1.2 noro 777: {
778: MP m;
779: DP s,t,u;
1.4 noro 780: int i,l;
1.5 noro 781: static MP *w;
782: static int wlen;
1.2 noro 783:
784: if ( !p1 || !p2 )
785: *pr = 0;
786: else if ( OID(p1) <= O_P )
787: muldc(vl,p2,(P)p1,pr);
788: else if ( OID(p2) <= O_P )
789: muldc(vl,p1,(P)p2,pr);
790: else {
1.10 noro 791: for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
1.5 noro 792: if ( l > wlen ) {
793: if ( w ) GC_free(w);
794: w = (MP *)MALLOC(l*sizeof(MP));
795: wlen = l;
796: }
1.10 noro 797: for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
1.4 noro 798: w[i] = m;
799: for ( s = 0, i = l-1; i >= 0; i-- ) {
1.10 noro 800: weyl_muldm(vl,w[i],p2,&t); addd(vl,s,t,&u); s = u;
1.2 noro 801: }
1.5 noro 802: bzero(w,l*sizeof(MP));
1.2 noro 803: *pr = s;
804: }
805: }
806:
1.10 noro 807: /* monomial * polynomial */
808:
1.19 noro 809: void weyl_muldm(VL vl,MP m0,DP p,DP *pr)
1.2 noro 810: {
811: DP r,t,t1;
812: MP m;
1.10 noro 813: DL d0;
814: int n,n2,l,i,j,tlen;
815: static MP *w,*psum;
816: static struct cdl *tab;
1.5 noro 817: static int wlen;
1.10 noro 818: static int rtlen;
1.2 noro 819:
820: if ( !p )
821: *pr = 0;
822: else {
1.4 noro 823: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
1.5 noro 824: if ( l > wlen ) {
825: if ( w ) GC_free(w);
826: w = (MP *)MALLOC(l*sizeof(MP));
827: wlen = l;
828: }
1.4 noro 829: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
830: w[i] = m;
1.10 noro 831:
832: n = NV(p); n2 = n>>1;
833: d0 = m0->dl;
834: for ( i = 0, tlen = 1; i < n2; i++ )
835: tlen *= d0->d[n2+i]+1;
836: if ( tlen > rtlen ) {
837: if ( tab ) GC_free(tab);
838: if ( psum ) GC_free(psum);
839: rtlen = tlen;
840: tab = (struct cdl *)MALLOC(rtlen*sizeof(struct cdl));
841: psum = (MP *)MALLOC(rtlen*sizeof(MP));
842: }
843: bzero(psum,tlen*sizeof(MP));
844: for ( i = l-1; i >= 0; i-- ) {
845: bzero(tab,tlen*sizeof(struct cdl));
846: weyl_mulmm(vl,m0,w[i],n,tab,tlen);
847: for ( j = 0; j < tlen; j++ ) {
848: if ( tab[j].c ) {
849: NEWMP(m); m->dl = tab[j].d; C(m) = tab[j].c; NEXT(m) = psum[j];
850: psum[j] = m;
851: }
852: }
1.2 noro 853: }
1.10 noro 854: for ( j = tlen-1, r = 0; j >= 0; j-- )
855: if ( psum[j] ) {
856: MKDP(n,psum[j],t); addd(vl,r,t,&t1); r = t1;
857: }
1.2 noro 858: if ( r )
859: r->sugar = p->sugar + m0->dl->td;
860: *pr = r;
861: }
862: }
863:
1.10 noro 864: /* m0 = x0^d0*x1^d1*... * dx0^e0*dx1^e1*... */
865: /* rtab : array of length (e0+1)*(e1+1)*... */
1.2 noro 866:
1.19 noro 867: void weyl_mulmm(VL vl,MP m0,MP m1,int n,struct cdl *rtab,int rtablen)
1.2 noro 868: {
1.19 noro 869: P c,c0,c1;
1.10 noro 870: DL d,d0,d1,dt;
871: int i,j,a,b,k,l,n2,s,min,curlen;
872: struct cdl *p;
873: static Q *ctab;
874: static struct cdl *tab;
1.5 noro 875: static int tablen;
1.10 noro 876: static struct cdl *tmptab;
877: static int tmptablen;
1.2 noro 878:
1.10 noro 879:
880: if ( !m0 || !m1 ) {
881: rtab[0].c = 0;
882: rtab[0].d = 0;
883: return;
884: }
885: c0 = C(m0); c1 = C(m1);
886: mulp(vl,c0,c1,&c);
887: d0 = m0->dl; d1 = m1->dl;
888: n2 = n>>1;
889: curlen = 1;
890: NEWDL(d,n);
891: if ( n & 1 )
892: /* offset of h-degree */
893: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
894: else
895: d->td = 0;
896: rtab[0].c = c;
897: rtab[0].d = d;
898:
899: if ( rtablen > tmptablen ) {
900: if ( tmptab ) GC_free(tmptab);
901: tmptab = (struct cdl *)MALLOC(rtablen*sizeof(struct cdl));
902: tmptablen = rtablen;
903: }
904: for ( i = 0; i < n2; i++ ) {
905: a = d0->d[i]; b = d1->d[n2+i];
906: k = d0->d[n2+i]; l = d1->d[i];
1.20 noro 907:
908: /* degree of xi^a*(Di^k*xi^l)*Di^b */
909: a += l;
910: b += k;
911: s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);
912:
1.10 noro 913: if ( !k || !l ) {
914: for ( j = 0, p = rtab; j < curlen; j++, p++ ) {
915: if ( p->c ) {
916: dt = p->d;
917: dt->d[i] = a;
918: dt->d[n2+i] = b;
919: dt->td += s;
1.5 noro 920: }
1.10 noro 921: }
922: curlen *= k+1;
923: continue;
924: }
925: if ( k+1 > tablen ) {
926: if ( tab ) GC_free(tab);
927: if ( ctab ) GC_free(ctab);
928: tablen = k+1;
929: tab = (struct cdl *)MALLOC(tablen*sizeof(struct cdl));
930: ctab = (Q *)MALLOC(tablen*sizeof(Q));
931: }
932: /* compute xi^a*(Di^k*xi^l)*Di^b */
933: min = MIN(k,l);
934: mkwc(k,l,ctab);
935: bzero(tab,(k+1)*sizeof(struct cdl));
936: if ( n & 1 )
937: for ( j = 0; j <= min; j++ ) {
938: NEWDL(d,n);
1.20 noro 939: d->d[i] = a-j; d->d[n2+i] = b-j;
1.10 noro 940: d->td = s;
1.20 noro 941: d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
1.10 noro 942: tab[j].d = d;
943: tab[j].c = (P)ctab[j];
944: }
945: else
946: for ( j = 0; j <= min; j++ ) {
947: NEWDL(d,n);
1.20 noro 948: d->d[i] = a-j; d->d[n2+i] = b-j;
949: d->td = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i); /* XXX */
1.10 noro 950: tab[j].d = d;
951: tab[j].c = (P)ctab[j];
952: }
953: bzero(ctab,(min+1)*sizeof(Q));
954: comm_muld_tab(vl,n,rtab,curlen,tab,k+1,tmptab);
955: curlen *= k+1;
956: bcopy(tmptab,rtab,curlen*sizeof(struct cdl));
957: }
958: }
959:
960: /* direct product of two cdl tables
961: rt[] = [
962: t[0]*t1[0],...,t[n-1]*t1[0],
963: t[0]*t1[1],...,t[n-1]*t1[1],
964: ...
965: t[0]*t1[n1-1],...,t[n-1]*t1[n1-1]
966: ]
967: */
968:
1.19 noro 969: void comm_muld_tab(VL vl,int nv,struct cdl *t,int n,struct cdl *t1,int n1,struct cdl *rt)
1.10 noro 970: {
971: int i,j;
972: struct cdl *p;
973: P c;
974: DL d;
975:
976: bzero(rt,n*n1*sizeof(struct cdl));
977: for ( j = 0, p = rt; j < n1; j++ ) {
978: c = t1[j].c;
979: d = t1[j].d;
980: if ( !c )
981: break;
982: for ( i = 0; i < n; i++, p++ ) {
983: if ( t[i].c ) {
984: mulp(vl,t[i].c,c,&p->c);
985: adddl(nv,t[i].d,d,&p->d);
986: }
1.6 noro 987: }
1.1 noro 988: }
989: }
990:
1.19 noro 991: void muldc(VL vl,DP p,P c,DP *pr)
1.1 noro 992: {
993: MP m,mr,mr0;
994:
995: if ( !p || !c )
996: *pr = 0;
997: else if ( NUM(c) && UNIQ((Q)c) )
998: *pr = p;
999: else if ( NUM(c) && MUNIQ((Q)c) )
1000: chsgnd(p,pr);
1001: else {
1002: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1003: NEXTMP(mr0,mr);
1004: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
1005: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
1006: else
1007: mulp(vl,C(m),c,&C(mr));
1008: mr->dl = m->dl;
1009: }
1010: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
1011: if ( *pr )
1012: (*pr)->sugar = p->sugar;
1013: }
1.24 noro 1014: }
1015:
1016: void muldc_trunc(VL vl,DP p,P c,DL dl,DP *pr)
1017: {
1018: MP m,mr,mr0;
1019: DL mdl;
1020: int i,n;
1021:
1022: if ( !p || !c ) {
1023: *pr = 0; return;
1024: }
1025: n = NV(p);
1026: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1027: mdl = m->dl;
1028: for ( i = 0; i < n; i++ )
1029: if ( mdl->d[i] < dl->d[i] )
1030: break;
1031: if ( i < n )
1032: break;
1033: NEXTMP(mr0,mr);
1034: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
1035: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
1036: else
1037: mulp(vl,C(m),c,&C(mr));
1038: mr->dl = m->dl;
1039: }
1040: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
1041: if ( *pr )
1042: (*pr)->sugar = p->sugar;
1.1 noro 1043: }
1044:
1.19 noro 1045: void divsdc(VL vl,DP p,P c,DP *pr)
1.1 noro 1046: {
1047: MP m,mr,mr0;
1048:
1049: if ( !c )
1050: error("disvsdc : division by 0");
1051: else if ( !p )
1052: *pr = 0;
1053: else {
1054: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1055: NEXTMP(mr0,mr); divsp(vl,C(m),c,&C(mr)); mr->dl = m->dl;
1056: }
1057: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
1058: if ( *pr )
1059: (*pr)->sugar = p->sugar;
1060: }
1061: }
1062:
1.19 noro 1063: void adddl(int n,DL d1,DL d2,DL *dr)
1.1 noro 1064: {
1065: DL dt;
1066: int i;
1067:
1068: if ( !d1->td )
1069: *dr = d2;
1070: else if ( !d2->td )
1071: *dr = d1;
1072: else {
1073: *dr = dt = (DL)MALLOC_ATOMIC((n+1)*sizeof(int));
1074: dt->td = d1->td + d2->td;
1075: for ( i = 0; i < n; i++ )
1076: dt->d[i] = d1->d[i]+d2->d[i];
1077: }
1.11 noro 1078: }
1079:
1080: /* d1 += d2 */
1081:
1.19 noro 1082: void adddl_destructive(int n,DL d1,DL d2)
1.11 noro 1083: {
1084: int i;
1085:
1086: d1->td += d2->td;
1087: for ( i = 0; i < n; i++ )
1088: d1->d[i] += d2->d[i];
1.1 noro 1089: }
1090:
1.19 noro 1091: int compd(VL vl,DP p1,DP p2)
1.1 noro 1092: {
1093: int n,t;
1094: MP m1,m2;
1095:
1096: if ( !p1 )
1097: return p2 ? -1 : 0;
1098: else if ( !p2 )
1099: return 1;
1100: else {
1101: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
1102: m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
1103: if ( (t = (*cmpdl)(n,m1->dl,m2->dl)) ||
1104: (t = compp(vl,C(m1),C(m2)) ) )
1105: return t;
1106: if ( m1 )
1107: return 1;
1108: else if ( m2 )
1109: return -1;
1110: else
1111: return 0;
1112: }
1113: }
1114:
1.19 noro 1115: int cmpdl_lex(int n,DL d1,DL d2)
1.1 noro 1116: {
1117: int i;
1118:
1119: for ( i = 0; i < n && d1->d[i] == d2->d[i]; i++ );
1120: return i == n ? 0 : (d1->d[i] > d2->d[i] ? 1 : -1);
1121: }
1122:
1.19 noro 1123: int cmpdl_revlex(int n,DL d1,DL d2)
1.1 noro 1124: {
1125: int i;
1126:
1127: for ( i = n - 1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
1128: return i < 0 ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1129: }
1130:
1.19 noro 1131: int cmpdl_gradlex(int n,DL d1,DL d2)
1.1 noro 1132: {
1133: if ( d1->td > d2->td )
1134: return 1;
1135: else if ( d1->td < d2->td )
1136: return -1;
1137: else
1138: return cmpdl_lex(n,d1,d2);
1139: }
1140:
1.19 noro 1141: int cmpdl_revgradlex(int n,DL d1,DL d2)
1.1 noro 1142: {
1.25 noro 1143: register int i,c;
1.7 noro 1144: register int *p1,*p2;
1145:
1.1 noro 1146: if ( d1->td > d2->td )
1147: return 1;
1148: else if ( d1->td < d2->td )
1149: return -1;
1.7 noro 1150: else {
1.25 noro 1151: i = n-1;
1152: p1 = d1->d+n-1;
1153: p2 = d2->d+n-1;
1154: while ( i >= 7 ) {
1155: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1156: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1157: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1158: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1159: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1160: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1161: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1162: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1163: i -= 8;
1164: }
1165: switch ( i ) {
1166: case 6:
1167: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1168: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1169: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1170: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1171: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1172: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1173: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1174: return 0;
1175: case 5:
1176: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1177: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1178: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1179: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1180: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1181: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1182: return 0;
1183: case 4:
1184: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1185: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1186: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1187: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1188: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1189: return 0;
1190: case 3:
1191: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1192: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1193: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1194: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1195: return 0;
1196: case 2:
1197: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1198: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1199: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1200: return 0;
1201: case 1:
1202: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1203: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1204: return 0;
1205: case 0:
1206: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1207: return 0;
1208: default:
1209: return 0;
1210: }
1211: LAST:
1212: if ( c > 0 ) return -1;
1213: else return 1;
1.7 noro 1214: }
1.1 noro 1215: }
1216:
1.19 noro 1217: int cmpdl_blex(int n,DL d1,DL d2)
1.1 noro 1218: {
1219: int c;
1220:
1221: if ( c = cmpdl_lex(n-1,d1,d2) )
1222: return c;
1223: else {
1224: c = d1->d[n-1] - d2->d[n-1];
1225: return c > 0 ? 1 : c < 0 ? -1 : 0;
1226: }
1227: }
1228:
1.19 noro 1229: int cmpdl_bgradlex(int n,DL d1,DL d2)
1.1 noro 1230: {
1231: int e1,e2,c;
1232:
1233: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
1234: if ( e1 > e2 )
1235: return 1;
1236: else if ( e1 < e2 )
1237: return -1;
1238: else {
1239: c = cmpdl_lex(n-1,d1,d2);
1240: if ( c )
1241: return c;
1242: else
1243: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
1244: }
1245: }
1246:
1.19 noro 1247: int cmpdl_brevgradlex(int n,DL d1,DL d2)
1.1 noro 1248: {
1249: int e1,e2,c;
1250:
1251: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
1252: if ( e1 > e2 )
1253: return 1;
1254: else if ( e1 < e2 )
1255: return -1;
1256: else {
1257: c = cmpdl_revlex(n-1,d1,d2);
1258: if ( c )
1259: return c;
1260: else
1261: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
1262: }
1263: }
1264:
1.19 noro 1265: int cmpdl_brevrev(int n,DL d1,DL d2)
1.1 noro 1266: {
1267: int e1,e2,f1,f2,c,i;
1268:
1269: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1270: e1 += d1->d[i]; e2 += d2->d[i];
1271: }
1272: f1 = d1->td - e1; f2 = d2->td - e2;
1273: if ( e1 > e2 )
1274: return 1;
1275: else if ( e1 < e2 )
1276: return -1;
1277: else {
1278: c = cmpdl_revlex(dp_nelim,d1,d2);
1279: if ( c )
1280: return c;
1281: else if ( f1 > f2 )
1282: return 1;
1283: else if ( f1 < f2 )
1284: return -1;
1285: else {
1286: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1287: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1288: }
1289: }
1290: }
1291:
1.19 noro 1292: int cmpdl_bgradrev(int n,DL d1,DL d2)
1.1 noro 1293: {
1294: int e1,e2,f1,f2,c,i;
1295:
1296: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1297: e1 += d1->d[i]; e2 += d2->d[i];
1298: }
1299: f1 = d1->td - e1; f2 = d2->td - e2;
1300: if ( e1 > e2 )
1301: return 1;
1302: else if ( e1 < e2 )
1303: return -1;
1304: else {
1305: c = cmpdl_lex(dp_nelim,d1,d2);
1306: if ( c )
1307: return c;
1308: else if ( f1 > f2 )
1309: return 1;
1310: else if ( f1 < f2 )
1311: return -1;
1312: else {
1313: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1314: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1315: }
1316: }
1317: }
1318:
1.19 noro 1319: int cmpdl_blexrev(int n,DL d1,DL d2)
1.1 noro 1320: {
1321: int e1,e2,f1,f2,c,i;
1322:
1323: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1324: e1 += d1->d[i]; e2 += d2->d[i];
1325: }
1326: f1 = d1->td - e1; f2 = d2->td - e2;
1327: c = cmpdl_lex(dp_nelim,d1,d2);
1328: if ( c )
1329: return c;
1330: else if ( f1 > f2 )
1331: return 1;
1332: else if ( f1 < f2 )
1333: return -1;
1334: else {
1335: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1336: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1337: }
1338: }
1339:
1.19 noro 1340: int cmpdl_elim(int n,DL d1,DL d2)
1.1 noro 1341: {
1342: int e1,e2,i;
1343:
1344: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1345: e1 += d1->d[i]; e2 += d2->d[i];
1346: }
1347: if ( e1 > e2 )
1348: return 1;
1349: else if ( e1 < e2 )
1350: return -1;
1351: else
1352: return cmpdl_revgradlex(n,d1,d2);
1.12 noro 1353: }
1354:
1.19 noro 1355: int cmpdl_weyl_elim(int n,DL d1,DL d2)
1.12 noro 1356: {
1357: int e1,e2,i;
1358:
1359: for ( i = 1, e1 = 0, e2 = 0; i <= dp_nelim; i++ ) {
1360: e1 += d1->d[n-i]; e2 += d2->d[n-i];
1361: }
1362: if ( e1 > e2 )
1363: return 1;
1364: else if ( e1 < e2 )
1365: return -1;
1366: else if ( d1->td > d2->td )
1367: return 1;
1368: else if ( d1->td < d2->td )
1369: return -1;
1370: else return -cmpdl_revlex(n,d1,d2);
1.13 noro 1371: }
1372:
1373: /*
1374: a special ordering
1375: 1. total order
1376: 2. (-w,w) for the first 2*m variables
1377: 3. DRL for the first 2*m variables
1378: */
1379:
1.20 noro 1380: extern int *current_weyl_weight_vector;
1.13 noro 1381:
1.19 noro 1382: int cmpdl_homo_ww_drl(int n,DL d1,DL d2)
1.13 noro 1383: {
1384: int e1,e2,m,i;
1385: int *p1,*p2;
1386:
1387: if ( d1->td > d2->td )
1388: return 1;
1389: else if ( d1->td < d2->td )
1390: return -1;
1391:
1392: m = n>>1;
1.21 noro 1393: for ( i = 0, e1 = e2 = 0, p1 = d1->d, p2 = d2->d; i < m; i++ ) {
1394: e1 += current_weyl_weight_vector[i]*(p1[m+i] - p1[i]);
1395: e2 += current_weyl_weight_vector[i]*(p2[m+i] - p2[i]);
1.13 noro 1396: }
1397: if ( e1 > e2 )
1398: return 1;
1399: else if ( e1 < e2 )
1400: return -1;
1401:
1402: e1 = d1->td - d1->d[n-1];
1403: e2 = d2->td - d2->d[n-1];
1404: if ( e1 > e2 )
1405: return 1;
1406: else if ( e1 < e2 )
1407: return -1;
1408:
1409: for ( i= n - 1, p1 = d1->d+n-1, p2 = d2->d+n-1;
1410: i >= 0 && *p1 == *p2; i--, p1--, p2-- );
1411: return i < 0 ? 0 : (*p1 < *p2 ? 1 : -1);
1.21 noro 1412: }
1413:
1414: int cmpdl_drl_zigzag(int n,DL d1,DL d2)
1415: {
1416: int i,t,m;
1417: int *p1,*p2;
1418:
1419: if ( d1->td > d2->td )
1420: return 1;
1421: else if ( d1->td < d2->td )
1422: return -1;
1423: else {
1424: m = n>>1;
1425: for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) {
1426: if ( t = p1[m+i] - p2[m+i] ) return t > 0 ? -1 : 1;
1427: if ( t = p1[i] - p2[i] ) return t > 0 ? -1 : 1;
1428: }
1429: return 0;
1430: }
1431: }
1432:
1433: int cmpdl_homo_ww_drl_zigzag(int n,DL d1,DL d2)
1434: {
1435: int e1,e2,m,i,t;
1436: int *p1,*p2;
1437:
1438: if ( d1->td > d2->td )
1439: return 1;
1440: else if ( d1->td < d2->td )
1441: return -1;
1442:
1443: m = n>>1;
1444: for ( i = 0, e1 = e2 = 0, p1 = d1->d, p2 = d2->d; i < m; i++ ) {
1445: e1 += current_weyl_weight_vector[i]*(p1[m+i] - p1[i]);
1446: e2 += current_weyl_weight_vector[i]*(p2[m+i] - p2[i]);
1447: }
1448: if ( e1 > e2 )
1449: return 1;
1450: else if ( e1 < e2 )
1451: return -1;
1452:
1453: e1 = d1->td - d1->d[n-1];
1454: e2 = d2->td - d2->d[n-1];
1455: if ( e1 > e2 )
1456: return 1;
1457: else if ( e1 < e2 )
1458: return -1;
1459:
1460: for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) {
1461: if ( t = p1[m+i] - p2[m+i] ) return t > 0 ? -1 : 1;
1462: if ( t = p1[i] - p2[i] ) return t > 0 ? -1 : 1;
1463: }
1464: return 0;
1.1 noro 1465: }
1466:
1.19 noro 1467: int cmpdl_order_pair(int n,DL d1,DL d2)
1.1 noro 1468: {
1469: int e1,e2,i,j,l;
1470: int *t1,*t2;
1.20 noro 1471: int len,head;
1.1 noro 1472: struct order_pair *pair;
1473:
1.27 noro 1474: len = dp_current_spec->ord.block.length;
1475: pair = dp_current_spec->ord.block.order_pair;
1.1 noro 1476:
1.20 noro 1477: head = 0;
1.1 noro 1478: for ( i = 0, t1 = d1->d, t2 = d2->d; i < len; i++ ) {
1479: l = pair[i].length;
1480: switch ( pair[i].order ) {
1481: case 0:
1482: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
1.20 noro 1483: e1 += MUL_WEIGHT(t1[j],head+j);
1484: e2 += MUL_WEIGHT(t2[j],head+j);
1.1 noro 1485: }
1486: if ( e1 > e2 )
1487: return 1;
1488: else if ( e1 < e2 )
1489: return -1;
1490: else {
1491: for ( j = l - 1; j >= 0 && t1[j] == t2[j]; j-- );
1492: if ( j >= 0 )
1493: return t1[j] < t2[j] ? 1 : -1;
1494: }
1495: break;
1496: case 1:
1497: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
1.20 noro 1498: e1 += MUL_WEIGHT(t1[j],head+j);
1499: e2 += MUL_WEIGHT(t2[j],head+j);
1.1 noro 1500: }
1501: if ( e1 > e2 )
1502: return 1;
1503: else if ( e1 < e2 )
1504: return -1;
1505: else {
1506: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
1507: if ( j < l )
1508: return t1[j] > t2[j] ? 1 : -1;
1509: }
1510: break;
1511: case 2:
1512: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
1513: if ( j < l )
1514: return t1[j] > t2[j] ? 1 : -1;
1515: break;
1516: default:
1517: error("cmpdl_order_pair : invalid order"); break;
1518: }
1.20 noro 1519: t1 += l; t2 += l; head += l;
1.28 ! noro 1520: }
! 1521: return 0;
! 1522: }
! 1523:
! 1524: int cmpdl_composite(int nv,DL d1,DL d2)
! 1525: {
! 1526: int n,i,j,k,start,s,len;
! 1527: int *dw;
! 1528: struct sparse_weight *sw;
! 1529: struct weight_or_block *worb;
! 1530: int *w,*t1,*t2;
! 1531:
! 1532: n = dp_current_spec->ord.composite.length;
! 1533: worb = dp_current_spec->ord.composite.w_or_b;
! 1534: w = dp_dl_work;
! 1535: for ( i = 0, t1 = d1->d, t2 = d2->d; i < nv; i++ )
! 1536: w[i] = t1[i]-t2[i];
! 1537: for ( i = 0; i < n; i++, worb++ ) {
! 1538: len = worb->length;
! 1539: switch ( worb->type ) {
! 1540: case IS_DENSE_WEIGHT:
! 1541: dw = worb->body.dense_weight;
! 1542: for ( j = 0, s = 0; j < len; j++ )
! 1543: s += dw[j]*w[j];
! 1544: if ( s > 0 ) return 1;
! 1545: else if ( s < 0 ) return -1;
! 1546: break;
! 1547: case IS_SPARSE_WEIGHT:
! 1548: sw = worb->body.sparse_weight;
! 1549: for ( j = 0, s = 0; j < len; j++ )
! 1550: s += sw[j].value*w[sw[j].pos];
! 1551: if ( s > 0 ) return 1;
! 1552: else if ( s < 0 ) return -1;
! 1553: break;
! 1554: case IS_BLOCK:
! 1555: start = worb->body.block.start;
! 1556: switch ( worb->body.block.order ) {
! 1557: case 0:
! 1558: for ( j = 0, k = start, s = 0; j < len; j++, k++ ) {
! 1559: s += MUL_WEIGHT(w[k],k);
! 1560: }
! 1561: if ( s > 0 ) return 1;
! 1562: else if ( s < 0 ) return -1;
! 1563: else {
! 1564: for ( j = k-1; j >= start && w[j] == 0; j-- );
! 1565: if ( j >= start )
! 1566: return w[j] < 0 ? 1 : -1;
! 1567: }
! 1568: break;
! 1569: case 1:
! 1570: for ( j = 0, k = start, s = 0; j < len; j++, k++ ) {
! 1571: s += MUL_WEIGHT(w[k],k);
! 1572: }
! 1573: if ( s > 0 ) return 1;
! 1574: else if ( s < 0 ) return -1;
! 1575: else {
! 1576: for ( j = 0, k = start; j < len && w[j] == 0; j++, k++ );
! 1577: if ( j < len )
! 1578: return w[j] > 0 ? 1 : -1;
! 1579: }
! 1580: break;
! 1581: case 2:
! 1582: for ( j = 0, k = start; j < len && w[j] == 0; j++, k++ );
! 1583: if ( j < len )
! 1584: return w[j] > 0 ? 1 : -1;
! 1585: break;
! 1586: }
! 1587: break;
! 1588: }
1.1 noro 1589: }
1590: return 0;
1591: }
1592:
1.19 noro 1593: int cmpdl_matrix(int n,DL d1,DL d2)
1.1 noro 1594: {
1595: int *v,*w,*t1,*t2;
1596: int s,i,j,len;
1597: int **matrix;
1598:
1599: for ( i = 0, t1 = d1->d, t2 = d2->d, w = dp_dl_work; i < n; i++ )
1600: w[i] = t1[i]-t2[i];
1.27 noro 1601: len = dp_current_spec->ord.matrix.row;
1602: matrix = dp_current_spec->ord.matrix.matrix;
1.1 noro 1603: for ( j = 0; j < len; j++ ) {
1604: v = matrix[j];
1605: for ( i = 0, s = 0; i < n; i++ )
1606: s += v[i]*w[i];
1607: if ( s > 0 )
1608: return 1;
1609: else if ( s < 0 )
1610: return -1;
1611: }
1612: return 0;
1.25 noro 1613: }
1614:
1615: GeoBucket create_bucket()
1616: {
1617: GeoBucket g;
1618:
1619: g = CALLOC(1,sizeof(struct oGeoBucket));
1620: g->m = 32;
1621: return g;
1622: }
1623:
1624: void add_bucket(GeoBucket g,NODE d,int nv)
1625: {
1626: int l,k,m;
1627:
1628: l = length(d);
1629: for ( k = 0, m = 1; l > m; k++, m <<= 1 );
1630: /* 2^(k-1) < l <= 2^k */
1631: d = symb_merge(g->body[k],d,nv);
1632: for ( ; length(d) > (1<<(k)); k++ ) {
1633: g->body[k] = 0;
1634: d = symb_merge(g->body[k+1],d,nv);
1635: }
1636: g->body[k] = d;
1637: g->m = MAX(g->m,k);
1638: }
1639:
1640: DL remove_head_bucket(GeoBucket g,int nv)
1641: {
1642: int j,i,c,m;
1643: DL d;
1644:
1645: j = -1;
1646: m = g->m;
1647: for ( i = 0; i <= m; i++ ) {
1648: if ( !g->body[i] )
1649: continue;
1650: if ( j < 0 ) j = i;
1651: else {
1652: c = (*cmpdl)(nv,g->body[i]->body,g->body[j]->body);
1653: if ( c > 0 )
1654: j = i;
1655: else if ( c == 0 )
1656: g->body[i] = NEXT(g->body[i]);
1657: }
1658: }
1659: if ( j < 0 )
1660: return 0;
1661: else {
1662: d = g->body[j]->body;
1663: g->body[j] = NEXT(g->body[j]);
1664: return d;
1665: }
1.1 noro 1666: }