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