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