![[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.16 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.15 2001/03/29 09:49: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.1 noro 64:
65: int (*cmpdl)()=cmpdl_revgradlex;
66: int (*primitive_cmpdl[3])() = {cmpdl_revgradlex,cmpdl_gradlex,cmpdl_lex};
67:
1.2 noro 68: void comm_muld(VL,DP,DP,DP *);
69: void weyl_muld(VL,DP,DP,DP *);
1.10 noro 70: void weyl_muldm(VL,MP,DP,DP *);
71: void weyl_mulmm(VL,MP,MP,int,struct cdl *,int);
72: void comm_muld_tab(VL,int,struct cdl *,int,struct cdl *,int,struct cdl *);
73:
1.2 noro 74: void mkwc(int,int,Q *);
75:
1.12 noro 76: int cmpdl_weyl_elim();
1.13 noro 77: int cmpdl_homo_ww_drl();
1.12 noro 78:
1.2 noro 79: int do_weyl;
80:
1.1 noro 81: int dp_nelim,dp_fcoeffs;
82: struct order_spec dp_current_spec;
83: int *dp_dl_work;
84:
85: int has_fcoef(DP);
86: int has_fcoef_p(P);
87:
88: int has_fcoef(f)
89: DP f;
90: {
91: MP t;
92:
93: if ( !f )
94: return 0;
95: for ( t = BDY(f); t; t = NEXT(t) )
96: if ( has_fcoef_p(t->c) )
97: break;
98: return t ? 1 : 0;
99: }
100:
101: int has_fcoef_p(f)
102: P f;
103: {
104: DCP dc;
105:
106: if ( !f )
107: return 0;
108: else if ( NUM(f) )
1.14 noro 109: return (NID((Num)f) == N_LM
110: || NID((Num)f) == N_GF2N
1.15 noro 111: || NID((Num)f) == N_GFPN
112: || NID((Num)f) == N_GFS) ? 1 : 0;
1.1 noro 113: else {
114: for ( dc = DC(f); dc; dc = NEXT(dc) )
115: if ( has_fcoef_p(COEF(dc)) )
116: return 1;
117: return 0;
118: }
119: }
120:
121: void initd(spec)
122: struct order_spec *spec;
123: {
124: switch ( spec->id ) {
125: case 2:
126: cmpdl = cmpdl_matrix;
127: dp_dl_work = (int *)MALLOC_ATOMIC(spec->nv*sizeof(int));
128: break;
129: case 1:
130: cmpdl = cmpdl_order_pair;
131: break;
132: default:
133: switch ( spec->ord.simple ) {
134: case ORD_REVGRADLEX:
135: cmpdl = cmpdl_revgradlex; break;
136: case ORD_GRADLEX:
137: cmpdl = cmpdl_gradlex; break;
138: case ORD_BREVGRADLEX:
139: cmpdl = cmpdl_brevgradlex; break;
140: case ORD_BGRADLEX:
141: cmpdl = cmpdl_bgradlex; break;
142: case ORD_BLEX:
143: cmpdl = cmpdl_blex; break;
144: case ORD_BREVREV:
145: cmpdl = cmpdl_brevrev; break;
146: case ORD_BGRADREV:
147: cmpdl = cmpdl_bgradrev; break;
148: case ORD_BLEXREV:
149: cmpdl = cmpdl_blexrev; break;
150: case ORD_ELIM:
151: cmpdl = cmpdl_elim; break;
1.12 noro 152: case ORD_WEYL_ELIM:
153: cmpdl = cmpdl_weyl_elim; break;
1.13 noro 154: case ORD_HOMO_WW_DRL:
155: cmpdl = cmpdl_homo_ww_drl; break;
1.1 noro 156: case ORD_LEX: default:
157: cmpdl = cmpdl_lex; break;
158: }
159: break;
160: }
161: dp_current_spec = *spec;
162: }
163:
164: void ptod(vl,dvl,p,pr)
165: VL vl,dvl;
166: P p;
167: DP *pr;
168: {
169: int isconst = 0;
1.16 ! noro 170: int n,i,j,k;
1.1 noro 171: VL tvl;
172: V v;
173: DL d;
174: MP m;
175: DCP dc;
1.16 ! noro 176: DCP *w;
1.1 noro 177: DP r,s,t,u;
178: P x,c;
179:
180: if ( !p )
181: *pr = 0;
182: else {
183: for ( n = 0, tvl = dvl; tvl; tvl = NEXT(tvl), n++ );
184: if ( NUM(p) ) {
185: NEWDL(d,n);
186: NEWMP(m); m->dl = d; C(m) = p; NEXT(m) = 0; MKDP(n,m,*pr); (*pr)->sugar = 0;
187: } else {
188: for ( i = 0, tvl = dvl, v = VR(p);
189: tvl && tvl->v != v; tvl = NEXT(tvl), i++ );
190: if ( !tvl ) {
1.16 ! noro 191: for ( dc = DC(p), k = 0; dc; dc = NEXT(dc), k++ );
! 192: w = (DCP *)ALLOCA(k*sizeof(DCP));
! 193: for ( dc = DC(p), j = 0; j < k; dc = NEXT(dc), j++ )
! 194: w[j] = dc;
! 195:
! 196: for ( j = k-1, s = 0, MKV(v,x); j >= 0; j-- ) {
! 197: ptod(vl,dvl,COEF(w[j]),&t); pwrp(vl,x,DEG(w[j]),&c);
1.1 noro 198: muldc(vl,t,c,&r); addd(vl,r,s,&t); s = t;
199: }
200: *pr = s;
201: } else {
1.16 ! noro 202: for ( dc = DC(p), k = 0; dc; dc = NEXT(dc), k++ );
! 203: w = (DCP *)ALLOCA(k*sizeof(DCP));
! 204: for ( dc = DC(p), j = 0; j < k; dc = NEXT(dc), j++ )
! 205: w[j] = dc;
! 206:
! 207: for ( j = k-1, s = 0; j >= 0; j-- ) {
! 208: ptod(vl,dvl,COEF(w[j]),&t);
! 209: NEWDL(d,n); d->td = QTOS(DEG(w[j])); d->d[i] = d->td;
1.1 noro 210: NEWMP(m); m->dl = d; C(m) = (P)ONE; NEXT(m) = 0; MKDP(n,m,u); u->sugar = d->td;
1.2 noro 211: comm_muld(vl,t,u,&r); addd(vl,r,s,&t); s = t;
1.1 noro 212: }
213: *pr = s;
214: }
215: }
216: }
217: if ( !dp_fcoeffs && has_fcoef(*pr) )
218: dp_fcoeffs = 1;
219: }
220:
221: void dtop(vl,dvl,p,pr)
222: VL vl,dvl;
223: DP p;
224: P *pr;
225: {
1.16 ! noro 226: int n,i,j,k;
1.1 noro 227: DL d;
228: MP m;
1.16 ! noro 229: MP *a;
1.1 noro 230: P r,s,t,u,w;
231: Q q;
232: VL tvl;
233:
234: if ( !p )
235: *pr = 0;
236: else {
1.16 ! noro 237: for ( k = 0, m = BDY(p); m; m = NEXT(m), k++ );
! 238: a = (MP *)ALLOCA(k*sizeof(MP));
! 239: for ( j = 0, m = BDY(p); j < k; m = NEXT(m), j++ )
! 240: a[j] = m;
! 241:
! 242: for ( n = p->nv, j = k-1, s = 0; j >= 0; j-- ) {
! 243: m = a[j];
1.1 noro 244: t = C(m);
245: if ( NUM(t) && NID((Num)t) == N_M ) {
246: mptop(t,&u); t = u;
247: }
248: for ( i = 0, d = m->dl, tvl = dvl;
249: i < n; tvl = NEXT(tvl), i++ ) {
250: MKV(tvl->v,r); STOQ(d->d[i],q); pwrp(vl,r,q,&u);
251: mulp(vl,t,u,&w); t = w;
252: }
253: addp(vl,s,t,&u); s = u;
254: }
255: *pr = s;
256: }
257: }
258:
259: void nodetod(node,dp)
260: NODE node;
261: DP *dp;
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 {
279: d->d[i] = QTOS((Q)e); td += d->d[i];
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:
287: int sugard(m)
288: MP m;
289: {
290: int s;
291:
292: for ( s = 0; m; m = NEXT(m) )
293: s = MAX(s,m->dl->td);
294: return s;
295: }
296:
297: void addd(vl,p1,p2,pr)
298: VL vl;
299: DP p1,p2,*pr;
300: {
301: int n;
302: MP m1,m2,mr,mr0;
303: P t;
304:
305: if ( !p1 )
306: *pr = p2;
307: else if ( !p2 )
308: *pr = p1;
309: else {
310: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
311: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
312: case 0:
313: addp(vl,C(m1),C(m2),&t);
314: if ( t ) {
315: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = t;
316: }
317: m1 = NEXT(m1); m2 = NEXT(m2); break;
318: case 1:
319: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
320: m1 = NEXT(m1); break;
321: case -1:
322: NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
323: m2 = NEXT(m2); break;
324: }
325: if ( !mr0 )
326: if ( m1 )
327: mr0 = m1;
328: else if ( m2 )
329: mr0 = m2;
330: else {
331: *pr = 0;
332: return;
333: }
334: else if ( m1 )
335: NEXT(mr) = m1;
336: else if ( m2 )
337: NEXT(mr) = m2;
338: else
339: NEXT(mr) = 0;
340: MKDP(NV(p1),mr0,*pr);
341: if ( *pr )
342: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
343: }
344: }
345:
346: /* for F4 symbolic reduction */
347:
348: void symb_addd(p1,p2,pr)
349: DP p1,p2,*pr;
350: {
351: int n;
352: MP m1,m2,mr,mr0;
353: P t;
354:
355: if ( !p1 )
356: *pr = p2;
357: else if ( !p2 )
358: *pr = p1;
359: else {
360: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
361: NEXTMP(mr0,mr); C(mr) = (P)ONE;
362: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
363: case 0:
364: mr->dl = m1->dl;
365: m1 = NEXT(m1); m2 = NEXT(m2); break;
366: case 1:
367: mr->dl = m1->dl;
368: m1 = NEXT(m1); break;
369: case -1:
370: mr->dl = m2->dl;
371: m2 = NEXT(m2); break;
372: }
373: }
374: if ( !mr0 )
375: if ( m1 )
376: mr0 = m1;
377: else if ( m2 )
378: mr0 = m2;
379: else {
380: *pr = 0;
381: return;
382: }
383: else if ( m1 )
384: NEXT(mr) = m1;
385: else if ( m2 )
386: NEXT(mr) = m2;
387: else
388: NEXT(mr) = 0;
389: MKDP(NV(p1),mr0,*pr);
390: if ( *pr )
391: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
1.3 noro 392: }
393: }
394:
395: /*
396: * destructive merge of two list
397: *
398: * p1, p2 : list of DL
399: * return : a merged list
400: */
401:
402: NODE symb_merge(m1,m2,n)
403: NODE m1,m2;
404: int n;
405: {
406: NODE top,prev,cur,m,t;
407:
408: if ( !m1 )
409: return m2;
410: else if ( !m2 )
411: return m1;
412: else {
413: switch ( (*cmpdl)(n,(DL)BDY(m1),(DL)BDY(m2)) ) {
414: case 0:
415: top = m1; m = NEXT(m2);
416: break;
417: case 1:
418: top = m1; m = m2;
419: break;
420: case -1:
421: top = m2; m = m1;
422: break;
423: }
424: prev = top; cur = NEXT(top);
425: /* BDY(prev) > BDY(m) always holds */
426: while ( cur && m ) {
427: switch ( (*cmpdl)(n,(DL)BDY(cur),(DL)BDY(m)) ) {
428: case 0:
429: m = NEXT(m);
430: prev = cur; cur = NEXT(cur);
431: break;
432: case 1:
433: t = NEXT(cur); NEXT(cur) = m; m = t;
434: prev = cur; cur = NEXT(cur);
435: break;
436: case -1:
437: NEXT(prev) = m; m = cur;
438: prev = NEXT(prev); cur = NEXT(prev);
439: break;
440: }
441: }
442: if ( !cur )
443: NEXT(prev) = m;
444: return top;
1.1 noro 445: }
446: }
447:
448: void subd(vl,p1,p2,pr)
449: VL vl;
450: DP p1,p2,*pr;
451: {
452: DP t;
453:
454: if ( !p2 )
455: *pr = p1;
456: else {
457: chsgnd(p2,&t); addd(vl,p1,t,pr);
458: }
459: }
460:
461: void chsgnd(p,pr)
462: DP p,*pr;
463: {
464: MP m,mr,mr0;
465:
466: if ( !p )
467: *pr = 0;
468: else {
469: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
470: NEXTMP(mr0,mr); chsgnp(C(m),&C(mr)); mr->dl = m->dl;
471: }
472: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
473: if ( *pr )
474: (*pr)->sugar = p->sugar;
475: }
476: }
477:
478: void muld(vl,p1,p2,pr)
479: VL vl;
480: DP p1,p2,*pr;
481: {
1.2 noro 482: if ( ! do_weyl )
483: comm_muld(vl,p1,p2,pr);
484: else
485: weyl_muld(vl,p1,p2,pr);
486: }
487:
488: void comm_muld(vl,p1,p2,pr)
489: VL vl;
490: DP p1,p2,*pr;
491: {
1.1 noro 492: MP m;
493: DP s,t,u;
1.5 noro 494: int i,l,l1;
495: static MP *w;
496: static int wlen;
1.1 noro 497:
498: if ( !p1 || !p2 )
499: *pr = 0;
500: else if ( OID(p1) <= O_P )
501: muldc(vl,p2,(P)p1,pr);
502: else if ( OID(p2) <= O_P )
503: muldc(vl,p1,(P)p2,pr);
504: else {
1.5 noro 505: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
1.4 noro 506: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
1.5 noro 507: if ( l1 < l ) {
508: t = p1; p1 = p2; p2 = t;
509: l = l1;
510: }
511: if ( l > wlen ) {
512: if ( w ) GC_free(w);
513: w = (MP *)MALLOC(l*sizeof(MP));
514: wlen = l;
515: }
1.4 noro 516: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
517: w[i] = m;
518: for ( s = 0, i = l-1; i >= 0; i-- ) {
519: muldm(vl,p1,w[i],&t); addd(vl,s,t,&u); s = u;
1.1 noro 520: }
1.5 noro 521: bzero(w,l*sizeof(MP));
1.1 noro 522: *pr = s;
523: }
524: }
525:
526: void muldm(vl,p,m0,pr)
527: VL vl;
528: DP p;
529: MP m0;
530: DP *pr;
531: {
532: MP m,mr,mr0;
533: P c;
534: DL d;
535: int n;
536:
537: if ( !p )
538: *pr = 0;
539: else {
540: for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
541: m; m = NEXT(m) ) {
542: NEXTMP(mr0,mr);
543: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
544: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
545: else
546: mulp(vl,C(m),c,&C(mr));
547: adddl(n,m->dl,d,&mr->dl);
548: }
549: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
550: if ( *pr )
551: (*pr)->sugar = p->sugar + m0->dl->td;
1.2 noro 552: }
553: }
554:
555: void weyl_muld(vl,p1,p2,pr)
556: VL vl;
557: DP p1,p2,*pr;
558: {
559: MP m;
560: DP s,t,u;
1.4 noro 561: int i,l;
1.5 noro 562: static MP *w;
563: static int wlen;
1.2 noro 564:
565: if ( !p1 || !p2 )
566: *pr = 0;
567: else if ( OID(p1) <= O_P )
568: muldc(vl,p2,(P)p1,pr);
569: else if ( OID(p2) <= O_P )
570: muldc(vl,p1,(P)p2,pr);
571: else {
1.10 noro 572: for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
1.5 noro 573: if ( l > wlen ) {
574: if ( w ) GC_free(w);
575: w = (MP *)MALLOC(l*sizeof(MP));
576: wlen = l;
577: }
1.10 noro 578: for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
1.4 noro 579: w[i] = m;
580: for ( s = 0, i = l-1; i >= 0; i-- ) {
1.10 noro 581: weyl_muldm(vl,w[i],p2,&t); addd(vl,s,t,&u); s = u;
1.2 noro 582: }
1.5 noro 583: bzero(w,l*sizeof(MP));
1.2 noro 584: *pr = s;
585: }
586: }
587:
1.10 noro 588: /* monomial * polynomial */
589:
590: void weyl_muldm(vl,m0,p,pr)
1.2 noro 591: VL vl;
1.10 noro 592: MP m0;
1.2 noro 593: DP p;
594: DP *pr;
595: {
596: DP r,t,t1;
597: MP m;
1.10 noro 598: DL d0;
599: int n,n2,l,i,j,tlen;
600: static MP *w,*psum;
601: static struct cdl *tab;
1.5 noro 602: static int wlen;
1.10 noro 603: static int rtlen;
1.2 noro 604:
605: if ( !p )
606: *pr = 0;
607: else {
1.4 noro 608: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
1.5 noro 609: if ( l > wlen ) {
610: if ( w ) GC_free(w);
611: w = (MP *)MALLOC(l*sizeof(MP));
612: wlen = l;
613: }
1.4 noro 614: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
615: w[i] = m;
1.10 noro 616:
617: n = NV(p); n2 = n>>1;
618: d0 = m0->dl;
619: for ( i = 0, tlen = 1; i < n2; i++ )
620: tlen *= d0->d[n2+i]+1;
621: if ( tlen > rtlen ) {
622: if ( tab ) GC_free(tab);
623: if ( psum ) GC_free(psum);
624: rtlen = tlen;
625: tab = (struct cdl *)MALLOC(rtlen*sizeof(struct cdl));
626: psum = (MP *)MALLOC(rtlen*sizeof(MP));
627: }
628: bzero(psum,tlen*sizeof(MP));
629: for ( i = l-1; i >= 0; i-- ) {
630: bzero(tab,tlen*sizeof(struct cdl));
631: weyl_mulmm(vl,m0,w[i],n,tab,tlen);
632: for ( j = 0; j < tlen; j++ ) {
633: if ( tab[j].c ) {
634: NEWMP(m); m->dl = tab[j].d; C(m) = tab[j].c; NEXT(m) = psum[j];
635: psum[j] = m;
636: }
637: }
1.2 noro 638: }
1.10 noro 639: for ( j = tlen-1, r = 0; j >= 0; j-- )
640: if ( psum[j] ) {
641: MKDP(n,psum[j],t); addd(vl,r,t,&t1); r = t1;
642: }
1.2 noro 643: if ( r )
644: r->sugar = p->sugar + m0->dl->td;
645: *pr = r;
646: }
647: }
648:
1.10 noro 649: /* m0 = x0^d0*x1^d1*... * dx0^e0*dx1^e1*... */
650: /* rtab : array of length (e0+1)*(e1+1)*... */
1.2 noro 651:
1.10 noro 652: void weyl_mulmm(vl,m0,m1,n,rtab,rtablen)
1.2 noro 653: VL vl;
654: MP m0,m1;
655: int n;
1.10 noro 656: struct cdl *rtab;
657: int rtablen;
1.2 noro 658: {
659: MP m,mr,mr0;
660: DP r,t,t1;
661: P c,c0,c1,cc;
1.10 noro 662: DL d,d0,d1,dt;
663: int i,j,a,b,k,l,n2,s,min,curlen;
664: struct cdl *p;
665: static Q *ctab;
666: static struct cdl *tab;
1.5 noro 667: static int tablen;
1.10 noro 668: static struct cdl *tmptab;
669: static int tmptablen;
1.2 noro 670:
1.10 noro 671:
672: if ( !m0 || !m1 ) {
673: rtab[0].c = 0;
674: rtab[0].d = 0;
675: return;
676: }
677: c0 = C(m0); c1 = C(m1);
678: mulp(vl,c0,c1,&c);
679: d0 = m0->dl; d1 = m1->dl;
680: n2 = n>>1;
681: curlen = 1;
682: NEWDL(d,n);
683: if ( n & 1 )
684: /* offset of h-degree */
685: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
686: else
687: d->td = 0;
688: rtab[0].c = c;
689: rtab[0].d = d;
690:
691: if ( rtablen > tmptablen ) {
692: if ( tmptab ) GC_free(tmptab);
693: tmptab = (struct cdl *)MALLOC(rtablen*sizeof(struct cdl));
694: tmptablen = rtablen;
695: }
696: for ( i = 0; i < n2; i++ ) {
697: a = d0->d[i]; b = d1->d[n2+i];
698: k = d0->d[n2+i]; l = d1->d[i];
699: if ( !k || !l ) {
700: a += l;
701: b += k;
702: s = a+b;
703: for ( j = 0, p = rtab; j < curlen; j++, p++ ) {
704: if ( p->c ) {
705: dt = p->d;
706: dt->d[i] = a;
707: dt->d[n2+i] = b;
708: dt->td += s;
1.5 noro 709: }
1.10 noro 710: }
711: curlen *= k+1;
712: continue;
713: }
714: if ( k+1 > tablen ) {
715: if ( tab ) GC_free(tab);
716: if ( ctab ) GC_free(ctab);
717: tablen = k+1;
718: tab = (struct cdl *)MALLOC(tablen*sizeof(struct cdl));
719: ctab = (Q *)MALLOC(tablen*sizeof(Q));
720: }
721: /* degree of xi^a*(Di^k*xi^l)*Di^b */
722: s = a+k+l+b;
723: /* compute xi^a*(Di^k*xi^l)*Di^b */
724: min = MIN(k,l);
725: mkwc(k,l,ctab);
726: bzero(tab,(k+1)*sizeof(struct cdl));
727: if ( n & 1 )
728: for ( j = 0; j <= min; j++ ) {
729: NEWDL(d,n);
730: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
731: d->td = s;
732: d->d[n-1] = s-(d->d[i]+d->d[n2+i]);
733: tab[j].d = d;
734: tab[j].c = (P)ctab[j];
735: }
736: else
737: for ( j = 0; j <= min; j++ ) {
738: NEWDL(d,n);
739: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
740: d->td = d->d[i]+d->d[n2+i]; /* XXX */
741: tab[j].d = d;
742: tab[j].c = (P)ctab[j];
743: }
744: bzero(ctab,(min+1)*sizeof(Q));
745: comm_muld_tab(vl,n,rtab,curlen,tab,k+1,tmptab);
746: curlen *= k+1;
747: bcopy(tmptab,rtab,curlen*sizeof(struct cdl));
748: }
749: }
750:
751: /* direct product of two cdl tables
752: rt[] = [
753: t[0]*t1[0],...,t[n-1]*t1[0],
754: t[0]*t1[1],...,t[n-1]*t1[1],
755: ...
756: t[0]*t1[n1-1],...,t[n-1]*t1[n1-1]
757: ]
758: */
759:
760: void comm_muld_tab(vl,nv,t,n,t1,n1,rt)
761: VL vl;
762: int nv;
763: struct cdl *t;
764: int n;
765: struct cdl *t1;
766: int n1;
767: struct cdl *rt;
768: {
769: int i,j;
770: struct cdl *p;
771: P c;
772: DL d;
773:
774: bzero(rt,n*n1*sizeof(struct cdl));
775: for ( j = 0, p = rt; j < n1; j++ ) {
776: c = t1[j].c;
777: d = t1[j].d;
778: if ( !c )
779: break;
780: for ( i = 0; i < n; i++, p++ ) {
781: if ( t[i].c ) {
782: mulp(vl,t[i].c,c,&p->c);
783: adddl(nv,t[i].d,d,&p->d);
784: }
1.6 noro 785: }
1.1 noro 786: }
787: }
788:
789: void muldc(vl,p,c,pr)
790: VL vl;
791: DP p;
792: P c;
793: DP *pr;
794: {
795: MP m,mr,mr0;
796:
797: if ( !p || !c )
798: *pr = 0;
799: else if ( NUM(c) && UNIQ((Q)c) )
800: *pr = p;
801: else if ( NUM(c) && MUNIQ((Q)c) )
802: chsgnd(p,pr);
803: else {
804: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
805: NEXTMP(mr0,mr);
806: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
807: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
808: else
809: mulp(vl,C(m),c,&C(mr));
810: mr->dl = m->dl;
811: }
812: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
813: if ( *pr )
814: (*pr)->sugar = p->sugar;
815: }
816: }
817:
818: void divsdc(vl,p,c,pr)
819: VL vl;
820: DP p;
821: P c;
822: DP *pr;
823: {
824: MP m,mr,mr0;
825:
826: if ( !c )
827: error("disvsdc : division by 0");
828: else if ( !p )
829: *pr = 0;
830: else {
831: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
832: NEXTMP(mr0,mr); divsp(vl,C(m),c,&C(mr)); mr->dl = m->dl;
833: }
834: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
835: if ( *pr )
836: (*pr)->sugar = p->sugar;
837: }
838: }
839:
840: void adddl(n,d1,d2,dr)
841: int n;
842: DL d1,d2;
843: DL *dr;
844: {
845: DL dt;
846: int i;
847:
848: if ( !d1->td )
849: *dr = d2;
850: else if ( !d2->td )
851: *dr = d1;
852: else {
853: *dr = dt = (DL)MALLOC_ATOMIC((n+1)*sizeof(int));
854: dt->td = d1->td + d2->td;
855: for ( i = 0; i < n; i++ )
856: dt->d[i] = d1->d[i]+d2->d[i];
857: }
1.11 noro 858: }
859:
860: /* d1 += d2 */
861:
862: void adddl_destructive(n,d1,d2)
863: int n;
864: DL d1,d2;
865: {
866: DL dt;
867: int i;
868:
869: d1->td += d2->td;
870: for ( i = 0; i < n; i++ )
871: d1->d[i] += d2->d[i];
1.1 noro 872: }
873:
874: int compd(vl,p1,p2)
875: VL vl;
876: DP p1,p2;
877: {
878: int n,t;
879: MP m1,m2;
880:
881: if ( !p1 )
882: return p2 ? -1 : 0;
883: else if ( !p2 )
884: return 1;
885: else {
886: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
887: m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
888: if ( (t = (*cmpdl)(n,m1->dl,m2->dl)) ||
889: (t = compp(vl,C(m1),C(m2)) ) )
890: return t;
891: if ( m1 )
892: return 1;
893: else if ( m2 )
894: return -1;
895: else
896: return 0;
897: }
898: }
899:
900: int cmpdl_lex(n,d1,d2)
901: int n;
902: DL d1,d2;
903: {
904: int i;
905:
906: for ( i = 0; i < n && d1->d[i] == d2->d[i]; i++ );
907: return i == n ? 0 : (d1->d[i] > d2->d[i] ? 1 : -1);
908: }
909:
910: int cmpdl_revlex(n,d1,d2)
911: int n;
912: DL d1,d2;
913: {
914: int i;
915:
916: for ( i = n - 1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
917: return i < 0 ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
918: }
919:
920: int cmpdl_gradlex(n,d1,d2)
921: int n;
922: DL d1,d2;
923: {
924: if ( d1->td > d2->td )
925: return 1;
926: else if ( d1->td < d2->td )
927: return -1;
928: else
929: return cmpdl_lex(n,d1,d2);
930: }
931:
932: int cmpdl_revgradlex(n,d1,d2)
933: int n;
934: DL d1,d2;
935: {
1.7 noro 936: register int i;
937: register int *p1,*p2;
938:
1.1 noro 939: if ( d1->td > d2->td )
940: return 1;
941: else if ( d1->td < d2->td )
942: return -1;
1.7 noro 943: else {
944: for ( i= n - 1, p1 = d1->d+n-1, p2 = d2->d+n-1;
945: i >= 0 && *p1 == *p2; i--, p1--, p2-- );
946: return i < 0 ? 0 : (*p1 < *p2 ? 1 : -1);
947: }
1.1 noro 948: }
949:
950: int cmpdl_blex(n,d1,d2)
951: int n;
952: DL d1,d2;
953: {
954: int c;
955:
956: if ( c = cmpdl_lex(n-1,d1,d2) )
957: return c;
958: else {
959: c = d1->d[n-1] - d2->d[n-1];
960: return c > 0 ? 1 : c < 0 ? -1 : 0;
961: }
962: }
963:
964: int cmpdl_bgradlex(n,d1,d2)
965: int n;
966: DL d1,d2;
967: {
968: int e1,e2,c;
969:
970: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
971: if ( e1 > e2 )
972: return 1;
973: else if ( e1 < e2 )
974: return -1;
975: else {
976: c = cmpdl_lex(n-1,d1,d2);
977: if ( c )
978: return c;
979: else
980: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
981: }
982: }
983:
984: int cmpdl_brevgradlex(n,d1,d2)
985: int n;
986: DL d1,d2;
987: {
988: int e1,e2,c;
989:
990: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
991: if ( e1 > e2 )
992: return 1;
993: else if ( e1 < e2 )
994: return -1;
995: else {
996: c = cmpdl_revlex(n-1,d1,d2);
997: if ( c )
998: return c;
999: else
1000: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
1001: }
1002: }
1003:
1004: int cmpdl_brevrev(n,d1,d2)
1005: int n;
1006: DL d1,d2;
1007: {
1008: int e1,e2,f1,f2,c,i;
1009:
1010: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1011: e1 += d1->d[i]; e2 += d2->d[i];
1012: }
1013: f1 = d1->td - e1; f2 = d2->td - e2;
1014: if ( e1 > e2 )
1015: return 1;
1016: else if ( e1 < e2 )
1017: return -1;
1018: else {
1019: c = cmpdl_revlex(dp_nelim,d1,d2);
1020: if ( c )
1021: return c;
1022: else if ( f1 > f2 )
1023: return 1;
1024: else if ( f1 < f2 )
1025: return -1;
1026: else {
1027: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1028: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1029: }
1030: }
1031: }
1032:
1033: int cmpdl_bgradrev(n,d1,d2)
1034: int n;
1035: DL d1,d2;
1036: {
1037: int e1,e2,f1,f2,c,i;
1038:
1039: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1040: e1 += d1->d[i]; e2 += d2->d[i];
1041: }
1042: f1 = d1->td - e1; f2 = d2->td - e2;
1043: if ( e1 > e2 )
1044: return 1;
1045: else if ( e1 < e2 )
1046: return -1;
1047: else {
1048: c = cmpdl_lex(dp_nelim,d1,d2);
1049: if ( c )
1050: return c;
1051: else if ( f1 > f2 )
1052: return 1;
1053: else if ( f1 < f2 )
1054: return -1;
1055: else {
1056: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1057: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1058: }
1059: }
1060: }
1061:
1062: int cmpdl_blexrev(n,d1,d2)
1063: int n;
1064: DL d1,d2;
1065: {
1066: int e1,e2,f1,f2,c,i;
1067:
1068: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1069: e1 += d1->d[i]; e2 += d2->d[i];
1070: }
1071: f1 = d1->td - e1; f2 = d2->td - e2;
1072: c = cmpdl_lex(dp_nelim,d1,d2);
1073: if ( c )
1074: return c;
1075: else if ( f1 > f2 )
1076: return 1;
1077: else if ( f1 < f2 )
1078: return -1;
1079: else {
1080: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1081: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1082: }
1083: }
1084:
1085: int cmpdl_elim(n,d1,d2)
1086: int n;
1087: DL d1,d2;
1088: {
1089: int e1,e2,i;
1090:
1091: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1092: e1 += d1->d[i]; e2 += d2->d[i];
1093: }
1094: if ( e1 > e2 )
1095: return 1;
1096: else if ( e1 < e2 )
1097: return -1;
1098: else
1099: return cmpdl_revgradlex(n,d1,d2);
1.12 noro 1100: }
1101:
1102: int cmpdl_weyl_elim(n,d1,d2)
1103: int n;
1104: DL d1,d2;
1105: {
1106: int e1,e2,i;
1107:
1108: for ( i = 1, e1 = 0, e2 = 0; i <= dp_nelim; i++ ) {
1109: e1 += d1->d[n-i]; e2 += d2->d[n-i];
1110: }
1111: if ( e1 > e2 )
1112: return 1;
1113: else if ( e1 < e2 )
1114: return -1;
1115: else if ( d1->td > d2->td )
1116: return 1;
1117: else if ( d1->td < d2->td )
1118: return -1;
1119: else return -cmpdl_revlex(n,d1,d2);
1.13 noro 1120: }
1121:
1122: /*
1123: a special ordering
1124: 1. total order
1125: 2. (-w,w) for the first 2*m variables
1126: 3. DRL for the first 2*m variables
1127: */
1128:
1129: extern int *current_weight_vector;
1130:
1131: int cmpdl_homo_ww_drl(n,d1,d2)
1132: int n;
1133: DL d1,d2;
1134: {
1135: int e1,e2,m,i;
1136: int *p1,*p2;
1137:
1138: if ( d1->td > d2->td )
1139: return 1;
1140: else if ( d1->td < d2->td )
1141: return -1;
1142:
1143: m = n>>1;
1144: for ( i = 0, e1 = e2 = 0; i < m; i++ ) {
1145: e1 += current_weight_vector[i]*(d1->d[m+i] - d1->d[i]);
1146: e2 += current_weight_vector[i]*(d2->d[m+i] - d2->d[i]);
1147: }
1148: if ( e1 > e2 )
1149: return 1;
1150: else if ( e1 < e2 )
1151: return -1;
1152:
1153: e1 = d1->td - d1->d[n-1];
1154: e2 = d2->td - d2->d[n-1];
1155: if ( e1 > e2 )
1156: return 1;
1157: else if ( e1 < e2 )
1158: return -1;
1159:
1160: for ( i= n - 1, p1 = d1->d+n-1, p2 = d2->d+n-1;
1161: i >= 0 && *p1 == *p2; i--, p1--, p2-- );
1162: return i < 0 ? 0 : (*p1 < *p2 ? 1 : -1);
1.1 noro 1163: }
1164:
1165: int cmpdl_order_pair(n,d1,d2)
1166: int n;
1167: DL d1,d2;
1168: {
1169: int e1,e2,i,j,l;
1170: int *t1,*t2;
1171: int len;
1172: struct order_pair *pair;
1173:
1174: len = dp_current_spec.ord.block.length;
1175: pair = dp_current_spec.ord.block.order_pair;
1176:
1177: for ( i = 0, t1 = d1->d, t2 = d2->d; i < len; i++ ) {
1178: l = pair[i].length;
1179: switch ( pair[i].order ) {
1180: case 0:
1181: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
1182: e1 += t1[j]; e2 += t2[j];
1183: }
1184: if ( e1 > e2 )
1185: return 1;
1186: else if ( e1 < e2 )
1187: return -1;
1188: else {
1189: for ( j = l - 1; j >= 0 && t1[j] == t2[j]; j-- );
1190: if ( j >= 0 )
1191: return t1[j] < t2[j] ? 1 : -1;
1192: }
1193: break;
1194: case 1:
1195: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
1196: e1 += t1[j]; e2 += t2[j];
1197: }
1198: if ( e1 > e2 )
1199: return 1;
1200: else if ( e1 < e2 )
1201: return -1;
1202: else {
1203: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
1204: if ( j < l )
1205: return t1[j] > t2[j] ? 1 : -1;
1206: }
1207: break;
1208: case 2:
1209: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
1210: if ( j < l )
1211: return t1[j] > t2[j] ? 1 : -1;
1212: break;
1213: default:
1214: error("cmpdl_order_pair : invalid order"); break;
1215: }
1216: t1 += l; t2 += l;
1217: }
1218: return 0;
1219: }
1220:
1221: int cmpdl_matrix(n,d1,d2)
1222: int n;
1223: DL d1,d2;
1224: {
1225: int *v,*w,*t1,*t2;
1226: int s,i,j,len;
1227: int **matrix;
1228:
1229: for ( i = 0, t1 = d1->d, t2 = d2->d, w = dp_dl_work; i < n; i++ )
1230: w[i] = t1[i]-t2[i];
1231: len = dp_current_spec.ord.matrix.row;
1232: matrix = dp_current_spec.ord.matrix.matrix;
1233: for ( j = 0; j < len; j++ ) {
1234: v = matrix[j];
1235: for ( i = 0, s = 0; i < n; i++ )
1236: s += v[i]*w[i];
1237: if ( s > 0 )
1238: return 1;
1239: else if ( s < 0 )
1240: return -1;
1241: }
1242: return 0;
1243: }