![[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.23 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.22 2003/01/04 09:06:17 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;
1.23 ! noro 420: return top;
! 421: }
! 422: }
! 423:
! 424: void _adddl(int n,DL d1,DL d2,DL d3)
! 425: {
! 426: int i;
! 427:
! 428: d3->td = d1->td+d2->td;
! 429: for ( i = 0; i < n; i++ )
! 430: d3->d[i] = d1->d[i]+d2->d[i];
! 431: }
! 432:
! 433: /* m1 <- m1 U dl*f, destructive */
! 434:
! 435: NODE mul_dllist(DL dl,DP f);
! 436:
! 437: NODE symb_mul_merge(NODE m1,DL dl,DP f,int n)
! 438: {
! 439: NODE top,prev,cur,n1;
! 440: DP g;
! 441: DL t,s;
! 442: MP m;
! 443:
! 444: if ( !m1 )
! 445: return mul_dllist(dl,f);
! 446: else if ( !f )
! 447: return m1;
! 448: else {
! 449: m = BDY(f);
! 450: NEWDL_NOINIT(t,n);
! 451: _adddl(n,m->dl,dl,t);
! 452: top = m1; prev = 0; cur = m1;
! 453: while ( m ) {
! 454: switch ( (*cmpdl)(n,(DL)BDY(cur),t) ) {
! 455: case 0:
! 456: prev = cur; cur = NEXT(cur);
! 457: if ( !cur ) {
! 458: MKDP(n,m,g);
! 459: NEXT(prev) = mul_dllist(dl,g);
! 460: return;
! 461: }
! 462: m = NEXT(m);
! 463: if ( m ) _adddl(n,m->dl,dl,t);
! 464: break;
! 465: case 1:
! 466: prev = cur; cur = NEXT(cur);
! 467: if ( !cur ) {
! 468: MKDP(n,m,g);
! 469: NEXT(prev) = mul_dllist(dl,g);
! 470: return;
! 471: }
! 472: break;
! 473: case -1:
! 474: NEWDL_NOINIT(s,n);
! 475: s->td = t->td;
! 476: bcopy(t->d,s->d,n*sizeof(int));
! 477: NEWNODE(n1);
! 478: n1->body = (pointer)s;
! 479: NEXT(n1) = cur;
! 480: if ( !prev ) {
! 481: top = n1; cur = n1;
! 482: } else {
! 483: NEXT(prev) = n1; prev = n1;
! 484: }
! 485: m = NEXT(m);
! 486: if ( m ) _adddl(n,m->dl,dl,t);
! 487: break;
! 488: }
! 489: }
1.18 noro 490: return top;
491: }
492: }
493:
1.19 noro 494: DLBUCKET symb_merge_bucket(DLBUCKET m1,DLBUCKET m2,int n)
1.18 noro 495: {
496: DLBUCKET top,prev,cur,m,t;
497:
498: if ( !m1 )
499: return m2;
500: else if ( !m2 )
501: return m1;
502: else {
503: if ( m1->td == m2->td ) {
504: top = m1;
505: BDY(top) = symb_merge(BDY(top),BDY(m2),n);
506: m = NEXT(m2);
507: } else if ( m1->td > m2->td ) {
508: top = m1; m = m2;
509: } else {
510: top = m2; m = m1;
511: }
512: prev = top; cur = NEXT(top);
513: /* prev->td > m->td always holds */
514: while ( cur && m ) {
515: if ( cur->td == m->td ) {
516: BDY(cur) = symb_merge(BDY(cur),BDY(m),n);
517: m = NEXT(m);
518: prev = cur; cur = NEXT(cur);
519: } else if ( cur->td > m->td ) {
520: t = NEXT(cur); NEXT(cur) = m; m = t;
521: prev = cur; cur = NEXT(cur);
522: } else {
523: NEXT(prev) = m; m = cur;
524: prev = NEXT(prev); cur = NEXT(prev);
1.3 noro 525: }
526: }
527: if ( !cur )
528: NEXT(prev) = m;
529: return top;
1.1 noro 530: }
531: }
532:
1.19 noro 533: void subd(VL vl,DP p1,DP p2,DP *pr)
1.1 noro 534: {
535: DP t;
536:
537: if ( !p2 )
538: *pr = p1;
539: else {
540: chsgnd(p2,&t); addd(vl,p1,t,pr);
541: }
542: }
543:
1.19 noro 544: void chsgnd(DP p,DP *pr)
1.1 noro 545: {
546: MP m,mr,mr0;
547:
548: if ( !p )
549: *pr = 0;
550: else {
551: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
552: NEXTMP(mr0,mr); chsgnp(C(m),&C(mr)); mr->dl = m->dl;
553: }
554: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
555: if ( *pr )
556: (*pr)->sugar = p->sugar;
557: }
558: }
559:
1.19 noro 560: void muld(VL vl,DP p1,DP p2,DP *pr)
1.1 noro 561: {
1.2 noro 562: if ( ! do_weyl )
563: comm_muld(vl,p1,p2,pr);
564: else
565: weyl_muld(vl,p1,p2,pr);
566: }
567:
1.19 noro 568: void comm_muld(VL vl,DP p1,DP p2,DP *pr)
1.2 noro 569: {
1.1 noro 570: MP m;
571: DP s,t,u;
1.5 noro 572: int i,l,l1;
573: static MP *w;
574: static int wlen;
1.1 noro 575:
576: if ( !p1 || !p2 )
577: *pr = 0;
578: else if ( OID(p1) <= O_P )
579: muldc(vl,p2,(P)p1,pr);
580: else if ( OID(p2) <= O_P )
581: muldc(vl,p1,(P)p2,pr);
582: else {
1.5 noro 583: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
1.4 noro 584: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
1.5 noro 585: if ( l1 < l ) {
586: t = p1; p1 = p2; p2 = t;
587: l = l1;
588: }
589: if ( l > wlen ) {
590: if ( w ) GC_free(w);
591: w = (MP *)MALLOC(l*sizeof(MP));
592: wlen = l;
593: }
1.4 noro 594: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
595: w[i] = m;
596: for ( s = 0, i = l-1; i >= 0; i-- ) {
597: muldm(vl,p1,w[i],&t); addd(vl,s,t,&u); s = u;
1.1 noro 598: }
1.5 noro 599: bzero(w,l*sizeof(MP));
1.1 noro 600: *pr = s;
601: }
602: }
603:
1.19 noro 604: void muldm(VL vl,DP p,MP m0,DP *pr)
1.1 noro 605: {
606: MP m,mr,mr0;
607: P c;
608: DL d;
609: int n;
610:
611: if ( !p )
612: *pr = 0;
613: else {
614: for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
615: m; m = NEXT(m) ) {
616: NEXTMP(mr0,mr);
617: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
618: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
619: else
620: mulp(vl,C(m),c,&C(mr));
621: adddl(n,m->dl,d,&mr->dl);
622: }
623: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
624: if ( *pr )
625: (*pr)->sugar = p->sugar + m0->dl->td;
1.2 noro 626: }
627: }
628:
1.19 noro 629: void weyl_muld(VL vl,DP p1,DP p2,DP *pr)
1.2 noro 630: {
631: MP m;
632: DP s,t,u;
1.4 noro 633: int i,l;
1.5 noro 634: static MP *w;
635: static int wlen;
1.2 noro 636:
637: if ( !p1 || !p2 )
638: *pr = 0;
639: else if ( OID(p1) <= O_P )
640: muldc(vl,p2,(P)p1,pr);
641: else if ( OID(p2) <= O_P )
642: muldc(vl,p1,(P)p2,pr);
643: else {
1.10 noro 644: for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
1.5 noro 645: if ( l > wlen ) {
646: if ( w ) GC_free(w);
647: w = (MP *)MALLOC(l*sizeof(MP));
648: wlen = l;
649: }
1.10 noro 650: for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
1.4 noro 651: w[i] = m;
652: for ( s = 0, i = l-1; i >= 0; i-- ) {
1.10 noro 653: weyl_muldm(vl,w[i],p2,&t); addd(vl,s,t,&u); s = u;
1.2 noro 654: }
1.5 noro 655: bzero(w,l*sizeof(MP));
1.2 noro 656: *pr = s;
657: }
658: }
659:
1.10 noro 660: /* monomial * polynomial */
661:
1.19 noro 662: void weyl_muldm(VL vl,MP m0,DP p,DP *pr)
1.2 noro 663: {
664: DP r,t,t1;
665: MP m;
1.10 noro 666: DL d0;
667: int n,n2,l,i,j,tlen;
668: static MP *w,*psum;
669: static struct cdl *tab;
1.5 noro 670: static int wlen;
1.10 noro 671: static int rtlen;
1.2 noro 672:
673: if ( !p )
674: *pr = 0;
675: else {
1.4 noro 676: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
1.5 noro 677: if ( l > wlen ) {
678: if ( w ) GC_free(w);
679: w = (MP *)MALLOC(l*sizeof(MP));
680: wlen = l;
681: }
1.4 noro 682: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
683: w[i] = m;
1.10 noro 684:
685: n = NV(p); n2 = n>>1;
686: d0 = m0->dl;
687: for ( i = 0, tlen = 1; i < n2; i++ )
688: tlen *= d0->d[n2+i]+1;
689: if ( tlen > rtlen ) {
690: if ( tab ) GC_free(tab);
691: if ( psum ) GC_free(psum);
692: rtlen = tlen;
693: tab = (struct cdl *)MALLOC(rtlen*sizeof(struct cdl));
694: psum = (MP *)MALLOC(rtlen*sizeof(MP));
695: }
696: bzero(psum,tlen*sizeof(MP));
697: for ( i = l-1; i >= 0; i-- ) {
698: bzero(tab,tlen*sizeof(struct cdl));
699: weyl_mulmm(vl,m0,w[i],n,tab,tlen);
700: for ( j = 0; j < tlen; j++ ) {
701: if ( tab[j].c ) {
702: NEWMP(m); m->dl = tab[j].d; C(m) = tab[j].c; NEXT(m) = psum[j];
703: psum[j] = m;
704: }
705: }
1.2 noro 706: }
1.10 noro 707: for ( j = tlen-1, r = 0; j >= 0; j-- )
708: if ( psum[j] ) {
709: MKDP(n,psum[j],t); addd(vl,r,t,&t1); r = t1;
710: }
1.2 noro 711: if ( r )
712: r->sugar = p->sugar + m0->dl->td;
713: *pr = r;
714: }
715: }
716:
1.10 noro 717: /* m0 = x0^d0*x1^d1*... * dx0^e0*dx1^e1*... */
718: /* rtab : array of length (e0+1)*(e1+1)*... */
1.2 noro 719:
1.19 noro 720: void weyl_mulmm(VL vl,MP m0,MP m1,int n,struct cdl *rtab,int rtablen)
1.2 noro 721: {
1.19 noro 722: P c,c0,c1;
1.10 noro 723: DL d,d0,d1,dt;
724: int i,j,a,b,k,l,n2,s,min,curlen;
725: struct cdl *p;
726: static Q *ctab;
727: static struct cdl *tab;
1.5 noro 728: static int tablen;
1.10 noro 729: static struct cdl *tmptab;
730: static int tmptablen;
1.2 noro 731:
1.10 noro 732:
733: if ( !m0 || !m1 ) {
734: rtab[0].c = 0;
735: rtab[0].d = 0;
736: return;
737: }
738: c0 = C(m0); c1 = C(m1);
739: mulp(vl,c0,c1,&c);
740: d0 = m0->dl; d1 = m1->dl;
741: n2 = n>>1;
742: curlen = 1;
743: NEWDL(d,n);
744: if ( n & 1 )
745: /* offset of h-degree */
746: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
747: else
748: d->td = 0;
749: rtab[0].c = c;
750: rtab[0].d = d;
751:
752: if ( rtablen > tmptablen ) {
753: if ( tmptab ) GC_free(tmptab);
754: tmptab = (struct cdl *)MALLOC(rtablen*sizeof(struct cdl));
755: tmptablen = rtablen;
756: }
757: for ( i = 0; i < n2; i++ ) {
758: a = d0->d[i]; b = d1->d[n2+i];
759: k = d0->d[n2+i]; l = d1->d[i];
1.20 noro 760:
761: /* degree of xi^a*(Di^k*xi^l)*Di^b */
762: a += l;
763: b += k;
764: s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);
765:
1.10 noro 766: if ( !k || !l ) {
767: for ( j = 0, p = rtab; j < curlen; j++, p++ ) {
768: if ( p->c ) {
769: dt = p->d;
770: dt->d[i] = a;
771: dt->d[n2+i] = b;
772: dt->td += s;
1.5 noro 773: }
1.10 noro 774: }
775: curlen *= k+1;
776: continue;
777: }
778: if ( k+1 > tablen ) {
779: if ( tab ) GC_free(tab);
780: if ( ctab ) GC_free(ctab);
781: tablen = k+1;
782: tab = (struct cdl *)MALLOC(tablen*sizeof(struct cdl));
783: ctab = (Q *)MALLOC(tablen*sizeof(Q));
784: }
785: /* compute xi^a*(Di^k*xi^l)*Di^b */
786: min = MIN(k,l);
787: mkwc(k,l,ctab);
788: bzero(tab,(k+1)*sizeof(struct cdl));
789: if ( n & 1 )
790: for ( j = 0; j <= min; j++ ) {
791: NEWDL(d,n);
1.20 noro 792: d->d[i] = a-j; d->d[n2+i] = b-j;
1.10 noro 793: d->td = s;
1.20 noro 794: d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
1.10 noro 795: tab[j].d = d;
796: tab[j].c = (P)ctab[j];
797: }
798: else
799: for ( j = 0; j <= min; j++ ) {
800: NEWDL(d,n);
1.20 noro 801: d->d[i] = a-j; d->d[n2+i] = b-j;
802: d->td = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i); /* XXX */
1.10 noro 803: tab[j].d = d;
804: tab[j].c = (P)ctab[j];
805: }
806: bzero(ctab,(min+1)*sizeof(Q));
807: comm_muld_tab(vl,n,rtab,curlen,tab,k+1,tmptab);
808: curlen *= k+1;
809: bcopy(tmptab,rtab,curlen*sizeof(struct cdl));
810: }
811: }
812:
813: /* direct product of two cdl tables
814: rt[] = [
815: t[0]*t1[0],...,t[n-1]*t1[0],
816: t[0]*t1[1],...,t[n-1]*t1[1],
817: ...
818: t[0]*t1[n1-1],...,t[n-1]*t1[n1-1]
819: ]
820: */
821:
1.19 noro 822: void comm_muld_tab(VL vl,int nv,struct cdl *t,int n,struct cdl *t1,int n1,struct cdl *rt)
1.10 noro 823: {
824: int i,j;
825: struct cdl *p;
826: P c;
827: DL d;
828:
829: bzero(rt,n*n1*sizeof(struct cdl));
830: for ( j = 0, p = rt; j < n1; j++ ) {
831: c = t1[j].c;
832: d = t1[j].d;
833: if ( !c )
834: break;
835: for ( i = 0; i < n; i++, p++ ) {
836: if ( t[i].c ) {
837: mulp(vl,t[i].c,c,&p->c);
838: adddl(nv,t[i].d,d,&p->d);
839: }
1.6 noro 840: }
1.1 noro 841: }
842: }
843:
1.19 noro 844: void muldc(VL vl,DP p,P c,DP *pr)
1.1 noro 845: {
846: MP m,mr,mr0;
847:
848: if ( !p || !c )
849: *pr = 0;
850: else if ( NUM(c) && UNIQ((Q)c) )
851: *pr = p;
852: else if ( NUM(c) && MUNIQ((Q)c) )
853: chsgnd(p,pr);
854: else {
855: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
856: NEXTMP(mr0,mr);
857: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
858: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
859: else
860: mulp(vl,C(m),c,&C(mr));
861: mr->dl = m->dl;
862: }
863: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
864: if ( *pr )
865: (*pr)->sugar = p->sugar;
866: }
867: }
868:
1.19 noro 869: void divsdc(VL vl,DP p,P c,DP *pr)
1.1 noro 870: {
871: MP m,mr,mr0;
872:
873: if ( !c )
874: error("disvsdc : division by 0");
875: else if ( !p )
876: *pr = 0;
877: else {
878: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
879: NEXTMP(mr0,mr); divsp(vl,C(m),c,&C(mr)); mr->dl = m->dl;
880: }
881: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
882: if ( *pr )
883: (*pr)->sugar = p->sugar;
884: }
885: }
886:
1.19 noro 887: void adddl(int n,DL d1,DL d2,DL *dr)
1.1 noro 888: {
889: DL dt;
890: int i;
891:
892: if ( !d1->td )
893: *dr = d2;
894: else if ( !d2->td )
895: *dr = d1;
896: else {
897: *dr = dt = (DL)MALLOC_ATOMIC((n+1)*sizeof(int));
898: dt->td = d1->td + d2->td;
899: for ( i = 0; i < n; i++ )
900: dt->d[i] = d1->d[i]+d2->d[i];
901: }
1.11 noro 902: }
903:
904: /* d1 += d2 */
905:
1.19 noro 906: void adddl_destructive(int n,DL d1,DL d2)
1.11 noro 907: {
908: int i;
909:
910: d1->td += d2->td;
911: for ( i = 0; i < n; i++ )
912: d1->d[i] += d2->d[i];
1.1 noro 913: }
914:
1.19 noro 915: int compd(VL vl,DP p1,DP p2)
1.1 noro 916: {
917: int n,t;
918: MP m1,m2;
919:
920: if ( !p1 )
921: return p2 ? -1 : 0;
922: else if ( !p2 )
923: return 1;
924: else {
925: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
926: m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
927: if ( (t = (*cmpdl)(n,m1->dl,m2->dl)) ||
928: (t = compp(vl,C(m1),C(m2)) ) )
929: return t;
930: if ( m1 )
931: return 1;
932: else if ( m2 )
933: return -1;
934: else
935: return 0;
936: }
937: }
938:
1.19 noro 939: int cmpdl_lex(int n,DL d1,DL d2)
1.1 noro 940: {
941: int i;
942:
943: for ( i = 0; i < n && d1->d[i] == d2->d[i]; i++ );
944: return i == n ? 0 : (d1->d[i] > d2->d[i] ? 1 : -1);
945: }
946:
1.19 noro 947: int cmpdl_revlex(int n,DL d1,DL d2)
1.1 noro 948: {
949: int i;
950:
951: for ( i = n - 1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
952: return i < 0 ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
953: }
954:
1.19 noro 955: int cmpdl_gradlex(int n,DL d1,DL d2)
1.1 noro 956: {
957: if ( d1->td > d2->td )
958: return 1;
959: else if ( d1->td < d2->td )
960: return -1;
961: else
962: return cmpdl_lex(n,d1,d2);
963: }
964:
1.19 noro 965: int cmpdl_revgradlex(int n,DL d1,DL d2)
1.1 noro 966: {
1.7 noro 967: register int i;
968: register int *p1,*p2;
969:
1.1 noro 970: if ( d1->td > d2->td )
971: return 1;
972: else if ( d1->td < d2->td )
973: return -1;
1.7 noro 974: else {
975: for ( i= n - 1, p1 = d1->d+n-1, p2 = d2->d+n-1;
976: i >= 0 && *p1 == *p2; i--, p1--, p2-- );
977: return i < 0 ? 0 : (*p1 < *p2 ? 1 : -1);
978: }
1.1 noro 979: }
980:
1.19 noro 981: int cmpdl_blex(int n,DL d1,DL d2)
1.1 noro 982: {
983: int c;
984:
985: if ( c = cmpdl_lex(n-1,d1,d2) )
986: return c;
987: else {
988: c = d1->d[n-1] - d2->d[n-1];
989: return c > 0 ? 1 : c < 0 ? -1 : 0;
990: }
991: }
992:
1.19 noro 993: int cmpdl_bgradlex(int n,DL d1,DL d2)
1.1 noro 994: {
995: int e1,e2,c;
996:
997: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
998: if ( e1 > e2 )
999: return 1;
1000: else if ( e1 < e2 )
1001: return -1;
1002: else {
1003: c = cmpdl_lex(n-1,d1,d2);
1004: if ( c )
1005: return c;
1006: else
1007: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
1008: }
1009: }
1010:
1.19 noro 1011: int cmpdl_brevgradlex(int n,DL d1,DL d2)
1.1 noro 1012: {
1013: int e1,e2,c;
1014:
1015: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
1016: if ( e1 > e2 )
1017: return 1;
1018: else if ( e1 < e2 )
1019: return -1;
1020: else {
1021: c = cmpdl_revlex(n-1,d1,d2);
1022: if ( c )
1023: return c;
1024: else
1025: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
1026: }
1027: }
1028:
1.19 noro 1029: int cmpdl_brevrev(int n,DL d1,DL d2)
1.1 noro 1030: {
1031: int e1,e2,f1,f2,c,i;
1032:
1033: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1034: e1 += d1->d[i]; e2 += d2->d[i];
1035: }
1036: f1 = d1->td - e1; f2 = d2->td - e2;
1037: if ( e1 > e2 )
1038: return 1;
1039: else if ( e1 < e2 )
1040: return -1;
1041: else {
1042: c = cmpdl_revlex(dp_nelim,d1,d2);
1043: if ( c )
1044: return c;
1045: else if ( f1 > f2 )
1046: return 1;
1047: else if ( f1 < f2 )
1048: return -1;
1049: else {
1050: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1051: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1052: }
1053: }
1054: }
1055:
1.19 noro 1056: int cmpdl_bgradrev(int n,DL d1,DL d2)
1.1 noro 1057: {
1058: int e1,e2,f1,f2,c,i;
1059:
1060: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1061: e1 += d1->d[i]; e2 += d2->d[i];
1062: }
1063: f1 = d1->td - e1; f2 = d2->td - e2;
1064: if ( e1 > e2 )
1065: return 1;
1066: else if ( e1 < e2 )
1067: return -1;
1068: else {
1069: c = cmpdl_lex(dp_nelim,d1,d2);
1070: if ( c )
1071: return c;
1072: else if ( f1 > f2 )
1073: return 1;
1074: else if ( f1 < f2 )
1075: return -1;
1076: else {
1077: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1078: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1079: }
1080: }
1081: }
1082:
1.19 noro 1083: int cmpdl_blexrev(int n,DL d1,DL d2)
1.1 noro 1084: {
1085: int e1,e2,f1,f2,c,i;
1086:
1087: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1088: e1 += d1->d[i]; e2 += d2->d[i];
1089: }
1090: f1 = d1->td - e1; f2 = d2->td - e2;
1091: c = cmpdl_lex(dp_nelim,d1,d2);
1092: if ( c )
1093: return c;
1094: else if ( f1 > f2 )
1095: return 1;
1096: else if ( f1 < f2 )
1097: return -1;
1098: else {
1099: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1100: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1101: }
1102: }
1103:
1.19 noro 1104: int cmpdl_elim(int n,DL d1,DL d2)
1.1 noro 1105: {
1106: int e1,e2,i;
1107:
1108: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1109: e1 += d1->d[i]; e2 += d2->d[i];
1110: }
1111: if ( e1 > e2 )
1112: return 1;
1113: else if ( e1 < e2 )
1114: return -1;
1115: else
1116: return cmpdl_revgradlex(n,d1,d2);
1.12 noro 1117: }
1118:
1.19 noro 1119: int cmpdl_weyl_elim(int n,DL d1,DL d2)
1.12 noro 1120: {
1121: int e1,e2,i;
1122:
1123: for ( i = 1, e1 = 0, e2 = 0; i <= dp_nelim; i++ ) {
1124: e1 += d1->d[n-i]; e2 += d2->d[n-i];
1125: }
1126: if ( e1 > e2 )
1127: return 1;
1128: else if ( e1 < e2 )
1129: return -1;
1130: else if ( d1->td > d2->td )
1131: return 1;
1132: else if ( d1->td < d2->td )
1133: return -1;
1134: else return -cmpdl_revlex(n,d1,d2);
1.13 noro 1135: }
1136:
1137: /*
1138: a special ordering
1139: 1. total order
1140: 2. (-w,w) for the first 2*m variables
1141: 3. DRL for the first 2*m variables
1142: */
1143:
1.20 noro 1144: extern int *current_weyl_weight_vector;
1.13 noro 1145:
1.19 noro 1146: int cmpdl_homo_ww_drl(int n,DL d1,DL d2)
1.13 noro 1147: {
1148: int e1,e2,m,i;
1149: int *p1,*p2;
1150:
1151: if ( d1->td > d2->td )
1152: return 1;
1153: else if ( d1->td < d2->td )
1154: return -1;
1155:
1156: m = n>>1;
1.21 noro 1157: for ( i = 0, e1 = e2 = 0, p1 = d1->d, p2 = d2->d; i < m; i++ ) {
1158: e1 += current_weyl_weight_vector[i]*(p1[m+i] - p1[i]);
1159: e2 += current_weyl_weight_vector[i]*(p2[m+i] - p2[i]);
1.13 noro 1160: }
1161: if ( e1 > e2 )
1162: return 1;
1163: else if ( e1 < e2 )
1164: return -1;
1165:
1166: e1 = d1->td - d1->d[n-1];
1167: e2 = d2->td - d2->d[n-1];
1168: if ( e1 > e2 )
1169: return 1;
1170: else if ( e1 < e2 )
1171: return -1;
1172:
1173: for ( i= n - 1, p1 = d1->d+n-1, p2 = d2->d+n-1;
1174: i >= 0 && *p1 == *p2; i--, p1--, p2-- );
1175: return i < 0 ? 0 : (*p1 < *p2 ? 1 : -1);
1.21 noro 1176: }
1177:
1178: int cmpdl_drl_zigzag(int n,DL d1,DL d2)
1179: {
1180: int i,t,m;
1181: int *p1,*p2;
1182:
1183: if ( d1->td > d2->td )
1184: return 1;
1185: else if ( d1->td < d2->td )
1186: return -1;
1187: else {
1188: m = n>>1;
1189: for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) {
1190: if ( t = p1[m+i] - p2[m+i] ) return t > 0 ? -1 : 1;
1191: if ( t = p1[i] - p2[i] ) return t > 0 ? -1 : 1;
1192: }
1193: return 0;
1194: }
1195: }
1196:
1197: int cmpdl_homo_ww_drl_zigzag(int n,DL d1,DL d2)
1198: {
1199: int e1,e2,m,i,t;
1200: int *p1,*p2;
1201:
1202: if ( d1->td > d2->td )
1203: return 1;
1204: else if ( d1->td < d2->td )
1205: return -1;
1206:
1207: m = n>>1;
1208: for ( i = 0, e1 = e2 = 0, p1 = d1->d, p2 = d2->d; i < m; i++ ) {
1209: e1 += current_weyl_weight_vector[i]*(p1[m+i] - p1[i]);
1210: e2 += current_weyl_weight_vector[i]*(p2[m+i] - p2[i]);
1211: }
1212: if ( e1 > e2 )
1213: return 1;
1214: else if ( e1 < e2 )
1215: return -1;
1216:
1217: e1 = d1->td - d1->d[n-1];
1218: e2 = d2->td - d2->d[n-1];
1219: if ( e1 > e2 )
1220: return 1;
1221: else if ( e1 < e2 )
1222: return -1;
1223:
1224: for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) {
1225: if ( t = p1[m+i] - p2[m+i] ) return t > 0 ? -1 : 1;
1226: if ( t = p1[i] - p2[i] ) return t > 0 ? -1 : 1;
1227: }
1228: return 0;
1.1 noro 1229: }
1230:
1.19 noro 1231: int cmpdl_order_pair(int n,DL d1,DL d2)
1.1 noro 1232: {
1233: int e1,e2,i,j,l;
1234: int *t1,*t2;
1.20 noro 1235: int len,head;
1.1 noro 1236: struct order_pair *pair;
1237:
1238: len = dp_current_spec.ord.block.length;
1239: pair = dp_current_spec.ord.block.order_pair;
1240:
1.20 noro 1241: head = 0;
1.1 noro 1242: for ( i = 0, t1 = d1->d, t2 = d2->d; i < len; i++ ) {
1243: l = pair[i].length;
1244: switch ( pair[i].order ) {
1245: case 0:
1246: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
1.20 noro 1247: e1 += MUL_WEIGHT(t1[j],head+j);
1248: e2 += MUL_WEIGHT(t2[j],head+j);
1.1 noro 1249: }
1250: if ( e1 > e2 )
1251: return 1;
1252: else if ( e1 < e2 )
1253: return -1;
1254: else {
1255: for ( j = l - 1; j >= 0 && t1[j] == t2[j]; j-- );
1256: if ( j >= 0 )
1257: return t1[j] < t2[j] ? 1 : -1;
1258: }
1259: break;
1260: case 1:
1261: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
1.20 noro 1262: e1 += MUL_WEIGHT(t1[j],head+j);
1263: e2 += MUL_WEIGHT(t2[j],head+j);
1.1 noro 1264: }
1265: if ( e1 > e2 )
1266: return 1;
1267: else if ( e1 < e2 )
1268: return -1;
1269: else {
1270: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
1271: if ( j < l )
1272: return t1[j] > t2[j] ? 1 : -1;
1273: }
1274: break;
1275: case 2:
1276: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
1277: if ( j < l )
1278: return t1[j] > t2[j] ? 1 : -1;
1279: break;
1280: default:
1281: error("cmpdl_order_pair : invalid order"); break;
1282: }
1.20 noro 1283: t1 += l; t2 += l; head += l;
1.1 noro 1284: }
1285: return 0;
1286: }
1287:
1.19 noro 1288: int cmpdl_matrix(int n,DL d1,DL d2)
1.1 noro 1289: {
1290: int *v,*w,*t1,*t2;
1291: int s,i,j,len;
1292: int **matrix;
1293:
1294: for ( i = 0, t1 = d1->d, t2 = d2->d, w = dp_dl_work; i < n; i++ )
1295: w[i] = t1[i]-t2[i];
1296: len = dp_current_spec.ord.matrix.row;
1297: matrix = dp_current_spec.ord.matrix.matrix;
1298: for ( j = 0; j < len; j++ ) {
1299: v = matrix[j];
1300: for ( i = 0, s = 0; i < n; i++ )
1301: s += v[i]*w[i];
1302: if ( s > 0 )
1303: return 1;
1304: else if ( s < 0 )
1305: return -1;
1306: }
1307: return 0;
1308: }