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