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