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