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