Annotation of OpenXM_contrib2/asir2000/engine/dist.c, Revision 1.43
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.43 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.42 2006/08/27 22:17:27 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.43 ! noro 68: int cmpdl_top_weight();
1.1 noro 69:
70: int (*cmpdl)()=cmpdl_revgradlex;
1.43 ! noro 71: int (*cmpdl_tie_breaker)();
1.1 noro 72: int (*primitive_cmpdl[3])() = {cmpdl_revgradlex,cmpdl_gradlex,cmpdl_lex};
73:
1.2 noro 74: int do_weyl;
75:
1.1 noro 76: int dp_nelim,dp_fcoeffs;
1.27 noro 77: struct order_spec *dp_current_spec;
1.31 noro 78: struct modorder_spec *dp_current_modspec;
1.1 noro 79: int *dp_dl_work;
80:
1.24 noro 81: void comm_muld_trunc(VL vl,DP p1,DP p2,DL dl,DP *pr);
82: void comm_quod(VL vl,DP p1,DP p2,DP *pr);
83: void muldm_trunc(VL vl,DP p,MP m0,DL dl,DP *pr);
84: void muldc_trunc(VL vl,DP p,P c,DL dl,DP *pr);
1.29 noro 85:
86: void order_init()
87: {
88: struct order_spec *spec;
89:
90: create_order_spec(0,0,&spec);
91: initd(spec);
1.31 noro 92: create_modorder_spec(0,0,&dp_current_modspec);
1.29 noro 93: }
1.24 noro 94:
1.22 noro 95: int has_sfcoef(DP f)
1.1 noro 96: {
97: MP t;
98:
99: if ( !f )
100: return 0;
101: for ( t = BDY(f); t; t = NEXT(t) )
1.22 noro 102: if ( has_sfcoef_p(t->c) )
1.1 noro 103: break;
104: return t ? 1 : 0;
105: }
106:
1.22 noro 107: int has_sfcoef_p(P f)
1.1 noro 108: {
109: DCP dc;
110:
111: if ( !f )
112: return 0;
113: else if ( NUM(f) )
1.22 noro 114: return (NID((Num)f) == N_GFS) ? 1 : 0;
1.1 noro 115: else {
116: for ( dc = DC(f); dc; dc = NEXT(dc) )
1.22 noro 117: if ( has_sfcoef_p(COEF(dc)) )
1.1 noro 118: return 1;
119: return 0;
120: }
121: }
122:
1.43 ! noro 123: extern N *current_top_weight_vector;
! 124: static int current_top_weight_len;
! 125:
1.19 noro 126: void initd(struct order_spec *spec)
1.1 noro 127: {
1.43 ! noro 128: int len,i;
! 129:
1.1 noro 130: switch ( spec->id ) {
1.28 noro 131: case 3:
132: cmpdl = cmpdl_composite;
133: dp_dl_work = (int *)MALLOC_ATOMIC(spec->nv*sizeof(int));
134: break;
1.1 noro 135: case 2:
136: cmpdl = cmpdl_matrix;
137: dp_dl_work = (int *)MALLOC_ATOMIC(spec->nv*sizeof(int));
138: break;
139: case 1:
140: cmpdl = cmpdl_order_pair;
141: break;
142: default:
143: switch ( spec->ord.simple ) {
144: case ORD_REVGRADLEX:
145: cmpdl = cmpdl_revgradlex; break;
146: case ORD_GRADLEX:
147: cmpdl = cmpdl_gradlex; break;
148: case ORD_BREVGRADLEX:
149: cmpdl = cmpdl_brevgradlex; break;
150: case ORD_BGRADLEX:
151: cmpdl = cmpdl_bgradlex; break;
152: case ORD_BLEX:
153: cmpdl = cmpdl_blex; break;
154: case ORD_BREVREV:
155: cmpdl = cmpdl_brevrev; break;
156: case ORD_BGRADREV:
157: cmpdl = cmpdl_bgradrev; break;
158: case ORD_BLEXREV:
159: cmpdl = cmpdl_blexrev; break;
160: case ORD_ELIM:
161: cmpdl = cmpdl_elim; break;
1.12 noro 162: case ORD_WEYL_ELIM:
163: cmpdl = cmpdl_weyl_elim; break;
1.13 noro 164: case ORD_HOMO_WW_DRL:
165: cmpdl = cmpdl_homo_ww_drl; break;
1.21 noro 166: case ORD_DRL_ZIGZAG:
167: cmpdl = cmpdl_drl_zigzag; break;
168: case ORD_HOMO_WW_DRL_ZIGZAG:
169: cmpdl = cmpdl_homo_ww_drl_zigzag; break;
1.1 noro 170: case ORD_LEX: default:
171: cmpdl = cmpdl_lex; break;
172: }
173: break;
174: }
1.43 ! noro 175: if ( current_top_weight_vector ) {
! 176: cmpdl_tie_breaker = cmpdl;
! 177: cmpdl = cmpdl_top_weight;
! 178: for ( len = 0, i = 0; i < spec->nv; i++ )
! 179: if ( current_top_weight_vector[i] )
! 180: len = MAX(PL(current_top_weight_vector[i]),len);
! 181: current_top_weight_len = len;
! 182: }
1.27 noro 183: dp_current_spec = spec;
1.1 noro 184: }
185:
1.19 noro 186: void ptod(VL vl,VL dvl,P p,DP *pr)
1.1 noro 187: {
188: int isconst = 0;
1.16 noro 189: int n,i,j,k;
1.1 noro 190: VL tvl;
191: V v;
192: DL d;
193: MP m;
194: DCP dc;
1.16 noro 195: DCP *w;
1.1 noro 196: DP r,s,t,u;
197: P x,c;
198:
199: if ( !p )
200: *pr = 0;
201: else {
202: for ( n = 0, tvl = dvl; tvl; tvl = NEXT(tvl), n++ );
203: if ( NUM(p) ) {
204: NEWDL(d,n);
205: NEWMP(m); m->dl = d; C(m) = p; NEXT(m) = 0; MKDP(n,m,*pr); (*pr)->sugar = 0;
206: } else {
207: for ( i = 0, tvl = dvl, v = VR(p);
208: tvl && tvl->v != v; tvl = NEXT(tvl), i++ );
209: if ( !tvl ) {
1.16 noro 210: for ( dc = DC(p), k = 0; dc; dc = NEXT(dc), k++ );
211: w = (DCP *)ALLOCA(k*sizeof(DCP));
212: for ( dc = DC(p), j = 0; j < k; dc = NEXT(dc), j++ )
213: w[j] = dc;
214:
215: for ( j = k-1, s = 0, MKV(v,x); j >= 0; j-- ) {
216: ptod(vl,dvl,COEF(w[j]),&t); pwrp(vl,x,DEG(w[j]),&c);
1.1 noro 217: muldc(vl,t,c,&r); addd(vl,r,s,&t); s = t;
218: }
219: *pr = s;
220: } else {
1.16 noro 221: for ( dc = DC(p), k = 0; dc; dc = NEXT(dc), k++ );
222: w = (DCP *)ALLOCA(k*sizeof(DCP));
223: for ( dc = DC(p), j = 0; j < k; dc = NEXT(dc), j++ )
224: w[j] = dc;
225:
226: for ( j = k-1, s = 0; j >= 0; j-- ) {
227: ptod(vl,dvl,COEF(w[j]),&t);
1.20 noro 228: NEWDL(d,n); d->d[i] = QTOS(DEG(w[j]));
229: d->td = MUL_WEIGHT(d->d[i],i);
1.1 noro 230: NEWMP(m); m->dl = d; C(m) = (P)ONE; NEXT(m) = 0; MKDP(n,m,u); u->sugar = d->td;
1.2 noro 231: comm_muld(vl,t,u,&r); addd(vl,r,s,&t); s = t;
1.1 noro 232: }
233: *pr = s;
234: }
235: }
236: }
1.17 noro 237: #if 0
1.22 noro 238: if ( !dp_fcoeffs && has_sfcoef(*pr) )
239: dp_fcoeffs = N_GFS;
1.17 noro 240: #endif
1.1 noro 241: }
242:
1.19 noro 243: void dtop(VL vl,VL dvl,DP p,P *pr)
1.1 noro 244: {
1.16 noro 245: int n,i,j,k;
1.1 noro 246: DL d;
247: MP m;
1.16 noro 248: MP *a;
1.1 noro 249: P r,s,t,u,w;
250: Q q;
251: VL tvl;
252:
253: if ( !p )
254: *pr = 0;
255: else {
1.16 noro 256: for ( k = 0, m = BDY(p); m; m = NEXT(m), k++ );
257: a = (MP *)ALLOCA(k*sizeof(MP));
258: for ( j = 0, m = BDY(p); j < k; m = NEXT(m), j++ )
259: a[j] = m;
260:
261: for ( n = p->nv, j = k-1, s = 0; j >= 0; j-- ) {
262: m = a[j];
1.1 noro 263: t = C(m);
264: if ( NUM(t) && NID((Num)t) == N_M ) {
265: mptop(t,&u); t = u;
266: }
267: for ( i = 0, d = m->dl, tvl = dvl;
268: i < n; tvl = NEXT(tvl), i++ ) {
269: MKV(tvl->v,r); STOQ(d->d[i],q); pwrp(vl,r,q,&u);
270: mulp(vl,t,u,&w); t = w;
271: }
272: addp(vl,s,t,&u); s = u;
273: }
274: *pr = s;
275: }
276: }
277:
1.19 noro 278: void nodetod(NODE node,DP *dp)
1.1 noro 279: {
280: NODE t;
281: int len,i,td;
282: Q e;
283: DL d;
284: MP m;
285: DP u;
286:
287: for ( t = node, len = 0; t; t = NEXT(t), len++ );
288: NEWDL(d,len);
289: for ( t = node, i = 0, td = 0; i < len; t = NEXT(t), i++ ) {
290: e = (Q)BDY(t);
291: if ( !e )
292: d->d[i] = 0;
293: else if ( !NUM(e) || !RATN(e) || !INT(e) )
294: error("nodetod : invalid input");
295: else {
1.20 noro 296: d->d[i] = QTOS((Q)e); td += MUL_WEIGHT(d->d[i],i);
1.1 noro 297: }
298: }
299: d->td = td;
300: NEWMP(m); m->dl = d; C(m) = (P)ONE; NEXT(m) = 0;
301: MKDP(len,m,u); u->sugar = td; *dp = u;
302: }
303:
1.19 noro 304: int sugard(MP m)
1.1 noro 305: {
306: int s;
307:
308: for ( s = 0; m; m = NEXT(m) )
309: s = MAX(s,m->dl->td);
310: return s;
311: }
312:
1.19 noro 313: void addd(VL vl,DP p1,DP p2,DP *pr)
1.1 noro 314: {
315: int n;
316: MP m1,m2,mr,mr0;
317: P t;
1.30 ohara 318: DL d;
1.1 noro 319:
320: if ( !p1 )
321: *pr = p2;
322: else if ( !p2 )
323: *pr = p1;
324: else {
1.30 ohara 325: if ( OID(p1) <= O_R ) {
326: n = NV(p2); NEWDL(d,n);
1.31 noro 327: NEWMP(m1); m1->dl = d; C(m1) = (P)p1; NEXT(m1) = 0;
1.30 ohara 328: MKDP(n,m1,p1); (p1)->sugar = 0;
329: }
330: if ( OID(p2) <= O_R ) {
331: n = NV(p1); NEWDL(d,n);
1.31 noro 332: NEWMP(m2); m2->dl = d; C(m2) = (P)p2; NEXT(m2) = 0;
1.30 ohara 333: MKDP(n,m2,p2); (p2)->sugar = 0;
334: }
1.1 noro 335: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
336: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
337: case 0:
338: addp(vl,C(m1),C(m2),&t);
339: if ( t ) {
340: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = t;
341: }
342: m1 = NEXT(m1); m2 = NEXT(m2); break;
343: case 1:
344: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
345: m1 = NEXT(m1); break;
346: case -1:
347: NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
348: m2 = NEXT(m2); break;
349: }
350: if ( !mr0 )
351: if ( m1 )
352: mr0 = m1;
353: else if ( m2 )
354: mr0 = m2;
355: else {
356: *pr = 0;
357: return;
358: }
359: else if ( m1 )
360: NEXT(mr) = m1;
361: else if ( m2 )
362: NEXT(mr) = m2;
363: else
364: NEXT(mr) = 0;
365: MKDP(NV(p1),mr0,*pr);
366: if ( *pr )
367: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
368: }
369: }
370:
371: /* for F4 symbolic reduction */
372:
1.19 noro 373: void symb_addd(DP p1,DP p2,DP *pr)
1.1 noro 374: {
375: int n;
376: MP m1,m2,mr,mr0;
377:
378: if ( !p1 )
379: *pr = p2;
380: else if ( !p2 )
381: *pr = p1;
382: else {
383: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
384: NEXTMP(mr0,mr); C(mr) = (P)ONE;
385: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
386: case 0:
387: mr->dl = m1->dl;
388: m1 = NEXT(m1); m2 = NEXT(m2); break;
389: case 1:
390: mr->dl = m1->dl;
391: m1 = NEXT(m1); break;
392: case -1:
393: mr->dl = m2->dl;
394: m2 = NEXT(m2); break;
395: }
396: }
397: if ( !mr0 )
398: if ( m1 )
399: mr0 = m1;
400: else if ( m2 )
401: mr0 = m2;
402: else {
403: *pr = 0;
404: return;
405: }
406: else if ( m1 )
407: NEXT(mr) = m1;
408: else if ( m2 )
409: NEXT(mr) = m2;
410: else
411: NEXT(mr) = 0;
412: MKDP(NV(p1),mr0,*pr);
413: if ( *pr )
414: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
1.3 noro 415: }
416: }
417:
418: /*
419: * destructive merge of two list
420: *
421: * p1, p2 : list of DL
422: * return : a merged list
423: */
424:
1.19 noro 425: NODE symb_merge(NODE m1,NODE m2,int n)
1.3 noro 426: {
427: NODE top,prev,cur,m,t;
1.25 noro 428: int c,i;
429: DL d1,d2;
1.3 noro 430:
431: if ( !m1 )
432: return m2;
433: else if ( !m2 )
434: return m1;
435: else {
436: switch ( (*cmpdl)(n,(DL)BDY(m1),(DL)BDY(m2)) ) {
437: case 0:
438: top = m1; m = NEXT(m2);
439: break;
440: case 1:
441: top = m1; m = m2;
442: break;
443: case -1:
444: top = m2; m = m1;
445: break;
446: }
447: prev = top; cur = NEXT(top);
448: /* BDY(prev) > BDY(m) always holds */
449: while ( cur && m ) {
1.25 noro 450: d1 = (DL)BDY(cur);
451: d2 = (DL)BDY(m);
1.26 noro 452: #if 1
453: switch ( (*cmpdl)(n,(DL)BDY(cur),(DL)BDY(m)) ) {
454: #else
455: /* XXX only valid for DRL */
1.25 noro 456: if ( d1->td > d2->td )
457: c = 1;
458: else if ( d1->td < d2->td )
459: c = -1;
460: else {
461: for ( i = n-1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
462: if ( i < 0 )
463: c = 0;
464: else if ( d1->d[i] < d2->d[i] )
465: c = 1;
466: else
467: c = -1;
468: }
469: switch ( c ) {
470: #endif
1.3 noro 471: case 0:
472: m = NEXT(m);
473: prev = cur; cur = NEXT(cur);
474: break;
475: case 1:
476: t = NEXT(cur); NEXT(cur) = m; m = t;
477: prev = cur; cur = NEXT(cur);
478: break;
479: case -1:
480: NEXT(prev) = m; m = cur;
481: prev = NEXT(prev); cur = NEXT(prev);
482: break;
1.18 noro 483: }
484: }
485: if ( !cur )
486: NEXT(prev) = m;
1.23 noro 487: return top;
488: }
489: }
490:
491: void _adddl(int n,DL d1,DL d2,DL d3)
492: {
493: int i;
494:
495: d3->td = d1->td+d2->td;
496: for ( i = 0; i < n; i++ )
497: d3->d[i] = d1->d[i]+d2->d[i];
498: }
499:
500: /* m1 <- m1 U dl*f, destructive */
501:
502: NODE mul_dllist(DL dl,DP f);
503:
504: NODE symb_mul_merge(NODE m1,DL dl,DP f,int n)
505: {
506: NODE top,prev,cur,n1;
507: DP g;
508: DL t,s;
509: MP m;
510:
511: if ( !m1 )
512: return mul_dllist(dl,f);
513: else if ( !f )
514: return m1;
515: else {
516: m = BDY(f);
517: NEWDL_NOINIT(t,n);
518: _adddl(n,m->dl,dl,t);
519: top = m1; prev = 0; cur = m1;
520: while ( m ) {
521: switch ( (*cmpdl)(n,(DL)BDY(cur),t) ) {
522: case 0:
523: prev = cur; cur = NEXT(cur);
524: if ( !cur ) {
525: MKDP(n,m,g);
526: NEXT(prev) = mul_dllist(dl,g);
527: return;
528: }
529: m = NEXT(m);
530: if ( m ) _adddl(n,m->dl,dl,t);
531: break;
532: case 1:
533: prev = cur; cur = NEXT(cur);
534: if ( !cur ) {
535: MKDP(n,m,g);
536: NEXT(prev) = mul_dllist(dl,g);
537: return;
538: }
539: break;
540: case -1:
541: NEWDL_NOINIT(s,n);
542: s->td = t->td;
543: bcopy(t->d,s->d,n*sizeof(int));
544: NEWNODE(n1);
545: n1->body = (pointer)s;
546: NEXT(n1) = cur;
547: if ( !prev ) {
548: top = n1; cur = n1;
549: } else {
550: NEXT(prev) = n1; prev = n1;
551: }
552: m = NEXT(m);
553: if ( m ) _adddl(n,m->dl,dl,t);
554: break;
555: }
556: }
1.18 noro 557: return top;
558: }
559: }
560:
1.19 noro 561: DLBUCKET symb_merge_bucket(DLBUCKET m1,DLBUCKET m2,int n)
1.18 noro 562: {
563: DLBUCKET top,prev,cur,m,t;
564:
565: if ( !m1 )
566: return m2;
567: else if ( !m2 )
568: return m1;
569: else {
570: if ( m1->td == m2->td ) {
571: top = m1;
572: BDY(top) = symb_merge(BDY(top),BDY(m2),n);
573: m = NEXT(m2);
574: } else if ( m1->td > m2->td ) {
575: top = m1; m = m2;
576: } else {
577: top = m2; m = m1;
578: }
579: prev = top; cur = NEXT(top);
580: /* prev->td > m->td always holds */
581: while ( cur && m ) {
582: if ( cur->td == m->td ) {
583: BDY(cur) = symb_merge(BDY(cur),BDY(m),n);
584: m = NEXT(m);
585: prev = cur; cur = NEXT(cur);
586: } else if ( cur->td > m->td ) {
587: t = NEXT(cur); NEXT(cur) = m; m = t;
588: prev = cur; cur = NEXT(cur);
589: } else {
590: NEXT(prev) = m; m = cur;
591: prev = NEXT(prev); cur = NEXT(prev);
1.3 noro 592: }
593: }
594: if ( !cur )
595: NEXT(prev) = m;
596: return top;
1.1 noro 597: }
598: }
599:
1.19 noro 600: void subd(VL vl,DP p1,DP p2,DP *pr)
1.1 noro 601: {
602: DP t;
603:
604: if ( !p2 )
605: *pr = p1;
606: else {
607: chsgnd(p2,&t); addd(vl,p1,t,pr);
608: }
609: }
610:
1.19 noro 611: void chsgnd(DP p,DP *pr)
1.1 noro 612: {
613: MP m,mr,mr0;
1.33 noro 614: Obj r;
1.1 noro 615:
616: if ( !p )
617: *pr = 0;
1.33 noro 618: else if ( OID(p) <= O_R ) {
619: chsgnr((Obj)p,&r); *pr = (DP)r;
620: } else {
1.1 noro 621: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
622: NEXTMP(mr0,mr); chsgnp(C(m),&C(mr)); mr->dl = m->dl;
623: }
624: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
625: if ( *pr )
626: (*pr)->sugar = p->sugar;
627: }
628: }
629:
1.19 noro 630: void muld(VL vl,DP p1,DP p2,DP *pr)
1.1 noro 631: {
1.2 noro 632: if ( ! do_weyl )
633: comm_muld(vl,p1,p2,pr);
634: else
635: weyl_muld(vl,p1,p2,pr);
636: }
637:
1.19 noro 638: void comm_muld(VL vl,DP p1,DP p2,DP *pr)
1.2 noro 639: {
1.1 noro 640: MP m;
641: DP s,t,u;
1.5 noro 642: int i,l,l1;
643: static MP *w;
644: static int wlen;
1.1 noro 645:
646: if ( !p1 || !p2 )
647: *pr = 0;
648: else if ( OID(p1) <= O_P )
649: muldc(vl,p2,(P)p1,pr);
650: else if ( OID(p2) <= O_P )
651: muldc(vl,p1,(P)p2,pr);
652: else {
1.5 noro 653: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
1.4 noro 654: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
1.5 noro 655: if ( l1 < l ) {
656: t = p1; p1 = p2; p2 = t;
657: l = l1;
658: }
659: if ( l > wlen ) {
660: if ( w ) GC_free(w);
661: w = (MP *)MALLOC(l*sizeof(MP));
662: wlen = l;
663: }
1.4 noro 664: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
665: w[i] = m;
666: for ( s = 0, i = l-1; i >= 0; i-- ) {
667: muldm(vl,p1,w[i],&t); addd(vl,s,t,&u); s = u;
1.1 noro 668: }
1.5 noro 669: bzero(w,l*sizeof(MP));
1.1 noro 670: *pr = s;
671: }
672: }
673:
1.24 noro 674: /* discard terms which is not a multiple of dl */
675:
676: void comm_muld_trunc(VL vl,DP p1,DP p2,DL dl,DP *pr)
677: {
678: MP m;
679: DP s,t,u;
680: int i,l,l1;
681: static MP *w;
682: static int wlen;
683:
684: if ( !p1 || !p2 )
685: *pr = 0;
686: else if ( OID(p1) <= O_P )
687: muldc_trunc(vl,p2,(P)p1,dl,pr);
688: else if ( OID(p2) <= O_P )
689: muldc_trunc(vl,p1,(P)p2,dl,pr);
690: else {
691: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
692: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
693: if ( l1 < l ) {
694: t = p1; p1 = p2; p2 = t;
695: l = l1;
696: }
697: if ( l > wlen ) {
698: if ( w ) GC_free(w);
699: w = (MP *)MALLOC(l*sizeof(MP));
700: wlen = l;
701: }
702: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
703: w[i] = m;
704: for ( s = 0, i = l-1; i >= 0; i-- ) {
705: muldm_trunc(vl,p1,w[i],dl,&t); addd(vl,s,t,&u); s = u;
706: }
707: bzero(w,l*sizeof(MP));
708: *pr = s;
709: }
710: }
711:
712: void comm_quod(VL vl,DP p1,DP p2,DP *pr)
713: {
714: MP m,m0;
715: DP s,t;
716: int i,n,sugar;
717: DL d1,d2,d;
718: Q a,b;
719:
720: if ( !p2 )
721: error("comm_quod : invalid input");
722: if ( !p1 )
723: *pr = 0;
724: else {
725: n = NV(p1);
726: d2 = BDY(p2)->dl;
727: m0 = 0;
728: sugar = p1->sugar;
729: while ( p1 ) {
730: d1 = BDY(p1)->dl;
731: NEWDL(d,n);
732: d->td = d1->td - d2->td;
733: for ( i = 0; i < n; i++ )
734: d->d[i] = d1->d[i]-d2->d[i];
735: NEXTMP(m0,m);
736: m->dl = d;
737: divq((Q)BDY(p1)->c,(Q)BDY(p2)->c,&a); chsgnq(a,&b);
738: C(m) = (P)b;
739: muldm_trunc(vl,p2,m,d2,&t);
740: addd(vl,p1,t,&s); p1 = s;
741: C(m) = (P)a;
742: }
743: if ( m0 ) {
744: NEXT(m) = 0; MKDP(n,m0,*pr);
745: } else
746: *pr = 0;
747: /* XXX */
748: if ( *pr )
749: (*pr)->sugar = sugar - d2->td;
750: }
751: }
752:
1.19 noro 753: void muldm(VL vl,DP p,MP m0,DP *pr)
1.1 noro 754: {
755: MP m,mr,mr0;
756: P c;
757: DL d;
758: int n;
759:
760: if ( !p )
761: *pr = 0;
762: else {
763: for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
764: m; m = NEXT(m) ) {
765: NEXTMP(mr0,mr);
766: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
767: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
768: else
769: mulp(vl,C(m),c,&C(mr));
770: adddl(n,m->dl,d,&mr->dl);
771: }
772: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
773: if ( *pr )
774: (*pr)->sugar = p->sugar + m0->dl->td;
1.2 noro 775: }
776: }
777:
1.24 noro 778: void muldm_trunc(VL vl,DP p,MP m0,DL dl,DP *pr)
779: {
780: MP m,mr,mr0;
781: P c;
782: DL d,tdl;
783: int n,i;
784:
785: if ( !p )
786: *pr = 0;
787: else {
788: n = NV(p);
789: NEWDL(tdl,n);
790: for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl;
791: m; m = NEXT(m) ) {
792: _adddl(n,m->dl,d,tdl);
793: for ( i = 0; i < n; i++ )
794: if ( tdl->d[i] < dl->d[i] )
795: break;
796: if ( i < n )
797: continue;
798: NEXTMP(mr0,mr);
799: mr->dl = tdl;
800: NEWDL(tdl,n);
801: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
802: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
803: else
804: mulp(vl,C(m),c,&C(mr));
805: }
806: if ( mr0 ) {
807: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
808: } else
809: *pr = 0;
810: if ( *pr )
811: (*pr)->sugar = p->sugar + m0->dl->td;
812: }
813: }
814:
1.19 noro 815: void weyl_muld(VL vl,DP p1,DP p2,DP *pr)
1.2 noro 816: {
817: MP m;
818: DP s,t,u;
1.4 noro 819: int i,l;
1.5 noro 820: static MP *w;
821: static int wlen;
1.2 noro 822:
823: if ( !p1 || !p2 )
824: *pr = 0;
825: else if ( OID(p1) <= O_P )
826: muldc(vl,p2,(P)p1,pr);
827: else if ( OID(p2) <= O_P )
828: muldc(vl,p1,(P)p2,pr);
829: else {
1.10 noro 830: for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
1.5 noro 831: if ( l > wlen ) {
832: if ( w ) GC_free(w);
833: w = (MP *)MALLOC(l*sizeof(MP));
834: wlen = l;
835: }
1.10 noro 836: for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
1.4 noro 837: w[i] = m;
838: for ( s = 0, i = l-1; i >= 0; i-- ) {
1.10 noro 839: weyl_muldm(vl,w[i],p2,&t); addd(vl,s,t,&u); s = u;
1.2 noro 840: }
1.5 noro 841: bzero(w,l*sizeof(MP));
1.2 noro 842: *pr = s;
843: }
844: }
845:
1.10 noro 846: /* monomial * polynomial */
847:
1.19 noro 848: void weyl_muldm(VL vl,MP m0,DP p,DP *pr)
1.2 noro 849: {
850: DP r,t,t1;
851: MP m;
1.10 noro 852: DL d0;
853: int n,n2,l,i,j,tlen;
854: static MP *w,*psum;
855: static struct cdl *tab;
1.5 noro 856: static int wlen;
1.10 noro 857: static int rtlen;
1.2 noro 858:
859: if ( !p )
860: *pr = 0;
861: else {
1.4 noro 862: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
1.5 noro 863: if ( l > wlen ) {
864: if ( w ) GC_free(w);
865: w = (MP *)MALLOC(l*sizeof(MP));
866: wlen = l;
867: }
1.4 noro 868: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
869: w[i] = m;
1.10 noro 870:
871: n = NV(p); n2 = n>>1;
872: d0 = m0->dl;
873: for ( i = 0, tlen = 1; i < n2; i++ )
874: tlen *= d0->d[n2+i]+1;
875: if ( tlen > rtlen ) {
876: if ( tab ) GC_free(tab);
877: if ( psum ) GC_free(psum);
878: rtlen = tlen;
879: tab = (struct cdl *)MALLOC(rtlen*sizeof(struct cdl));
880: psum = (MP *)MALLOC(rtlen*sizeof(MP));
881: }
882: bzero(psum,tlen*sizeof(MP));
883: for ( i = l-1; i >= 0; i-- ) {
884: bzero(tab,tlen*sizeof(struct cdl));
885: weyl_mulmm(vl,m0,w[i],n,tab,tlen);
886: for ( j = 0; j < tlen; j++ ) {
887: if ( tab[j].c ) {
888: NEWMP(m); m->dl = tab[j].d; C(m) = tab[j].c; NEXT(m) = psum[j];
889: psum[j] = m;
890: }
891: }
1.2 noro 892: }
1.10 noro 893: for ( j = tlen-1, r = 0; j >= 0; j-- )
894: if ( psum[j] ) {
895: MKDP(n,psum[j],t); addd(vl,r,t,&t1); r = t1;
896: }
1.2 noro 897: if ( r )
898: r->sugar = p->sugar + m0->dl->td;
899: *pr = r;
900: }
901: }
902:
1.10 noro 903: /* m0 = x0^d0*x1^d1*... * dx0^e0*dx1^e1*... */
904: /* rtab : array of length (e0+1)*(e1+1)*... */
1.2 noro 905:
1.19 noro 906: void weyl_mulmm(VL vl,MP m0,MP m1,int n,struct cdl *rtab,int rtablen)
1.2 noro 907: {
1.19 noro 908: P c,c0,c1;
1.10 noro 909: DL d,d0,d1,dt;
910: int i,j,a,b,k,l,n2,s,min,curlen;
911: struct cdl *p;
912: static Q *ctab;
913: static struct cdl *tab;
1.5 noro 914: static int tablen;
1.10 noro 915: static struct cdl *tmptab;
916: static int tmptablen;
1.2 noro 917:
1.10 noro 918:
919: if ( !m0 || !m1 ) {
920: rtab[0].c = 0;
921: rtab[0].d = 0;
922: return;
923: }
924: c0 = C(m0); c1 = C(m1);
925: mulp(vl,c0,c1,&c);
926: d0 = m0->dl; d1 = m1->dl;
927: n2 = n>>1;
928: curlen = 1;
929: NEWDL(d,n);
930: if ( n & 1 )
931: /* offset of h-degree */
932: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
933: else
934: d->td = 0;
935: rtab[0].c = c;
936: rtab[0].d = d;
937:
938: if ( rtablen > tmptablen ) {
939: if ( tmptab ) GC_free(tmptab);
940: tmptab = (struct cdl *)MALLOC(rtablen*sizeof(struct cdl));
941: tmptablen = rtablen;
942: }
943: for ( i = 0; i < n2; i++ ) {
944: a = d0->d[i]; b = d1->d[n2+i];
945: k = d0->d[n2+i]; l = d1->d[i];
1.20 noro 946:
947: /* degree of xi^a*(Di^k*xi^l)*Di^b */
948: a += l;
949: b += k;
950: s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);
951:
1.10 noro 952: if ( !k || !l ) {
953: for ( j = 0, p = rtab; j < curlen; j++, p++ ) {
954: if ( p->c ) {
955: dt = p->d;
956: dt->d[i] = a;
957: dt->d[n2+i] = b;
958: dt->td += s;
1.5 noro 959: }
1.10 noro 960: }
961: curlen *= k+1;
962: continue;
963: }
964: if ( k+1 > tablen ) {
965: if ( tab ) GC_free(tab);
966: if ( ctab ) GC_free(ctab);
967: tablen = k+1;
968: tab = (struct cdl *)MALLOC(tablen*sizeof(struct cdl));
969: ctab = (Q *)MALLOC(tablen*sizeof(Q));
970: }
971: /* compute xi^a*(Di^k*xi^l)*Di^b */
972: min = MIN(k,l);
973: mkwc(k,l,ctab);
974: bzero(tab,(k+1)*sizeof(struct cdl));
975: if ( n & 1 )
976: for ( j = 0; j <= min; j++ ) {
977: NEWDL(d,n);
1.20 noro 978: d->d[i] = a-j; d->d[n2+i] = b-j;
1.10 noro 979: d->td = s;
1.20 noro 980: d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
1.10 noro 981: tab[j].d = d;
982: tab[j].c = (P)ctab[j];
983: }
984: else
985: for ( j = 0; j <= min; j++ ) {
986: NEWDL(d,n);
1.20 noro 987: d->d[i] = a-j; d->d[n2+i] = b-j;
988: d->td = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i); /* XXX */
1.10 noro 989: tab[j].d = d;
990: tab[j].c = (P)ctab[j];
991: }
992: bzero(ctab,(min+1)*sizeof(Q));
993: comm_muld_tab(vl,n,rtab,curlen,tab,k+1,tmptab);
994: curlen *= k+1;
995: bcopy(tmptab,rtab,curlen*sizeof(struct cdl));
996: }
997: }
998:
999: /* direct product of two cdl tables
1000: rt[] = [
1001: t[0]*t1[0],...,t[n-1]*t1[0],
1002: t[0]*t1[1],...,t[n-1]*t1[1],
1003: ...
1004: t[0]*t1[n1-1],...,t[n-1]*t1[n1-1]
1005: ]
1006: */
1007:
1.19 noro 1008: void comm_muld_tab(VL vl,int nv,struct cdl *t,int n,struct cdl *t1,int n1,struct cdl *rt)
1.10 noro 1009: {
1010: int i,j;
1011: struct cdl *p;
1012: P c;
1013: DL d;
1014:
1015: bzero(rt,n*n1*sizeof(struct cdl));
1016: for ( j = 0, p = rt; j < n1; j++ ) {
1017: c = t1[j].c;
1018: d = t1[j].d;
1019: if ( !c )
1020: break;
1021: for ( i = 0; i < n; i++, p++ ) {
1022: if ( t[i].c ) {
1023: mulp(vl,t[i].c,c,&p->c);
1024: adddl(nv,t[i].d,d,&p->d);
1025: }
1.6 noro 1026: }
1.1 noro 1027: }
1028: }
1029:
1.19 noro 1030: void muldc(VL vl,DP p,P c,DP *pr)
1.1 noro 1031: {
1032: MP m,mr,mr0;
1033:
1034: if ( !p || !c )
1035: *pr = 0;
1036: else if ( NUM(c) && UNIQ((Q)c) )
1037: *pr = p;
1038: else if ( NUM(c) && MUNIQ((Q)c) )
1039: chsgnd(p,pr);
1040: else {
1041: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1042: NEXTMP(mr0,mr);
1043: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
1044: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
1045: else
1046: mulp(vl,C(m),c,&C(mr));
1047: mr->dl = m->dl;
1048: }
1049: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
1050: if ( *pr )
1051: (*pr)->sugar = p->sugar;
1052: }
1.24 noro 1053: }
1054:
1055: void muldc_trunc(VL vl,DP p,P c,DL dl,DP *pr)
1056: {
1057: MP m,mr,mr0;
1058: DL mdl;
1059: int i,n;
1060:
1061: if ( !p || !c ) {
1062: *pr = 0; return;
1063: }
1064: n = NV(p);
1065: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1066: mdl = m->dl;
1067: for ( i = 0; i < n; i++ )
1068: if ( mdl->d[i] < dl->d[i] )
1069: break;
1070: if ( i < n )
1071: break;
1072: NEXTMP(mr0,mr);
1073: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
1074: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
1075: else
1076: mulp(vl,C(m),c,&C(mr));
1077: mr->dl = m->dl;
1078: }
1079: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
1080: if ( *pr )
1081: (*pr)->sugar = p->sugar;
1.1 noro 1082: }
1083:
1.19 noro 1084: void divsdc(VL vl,DP p,P c,DP *pr)
1.1 noro 1085: {
1086: MP m,mr,mr0;
1087:
1088: if ( !c )
1089: error("disvsdc : division by 0");
1090: else if ( !p )
1091: *pr = 0;
1092: else {
1093: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1094: NEXTMP(mr0,mr); divsp(vl,C(m),c,&C(mr)); mr->dl = m->dl;
1095: }
1096: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
1097: if ( *pr )
1098: (*pr)->sugar = p->sugar;
1099: }
1100: }
1101:
1.19 noro 1102: void adddl(int n,DL d1,DL d2,DL *dr)
1.1 noro 1103: {
1104: DL dt;
1105: int i;
1106:
1107: if ( !d1->td )
1108: *dr = d2;
1109: else if ( !d2->td )
1110: *dr = d1;
1111: else {
1112: *dr = dt = (DL)MALLOC_ATOMIC((n+1)*sizeof(int));
1113: dt->td = d1->td + d2->td;
1114: for ( i = 0; i < n; i++ )
1115: dt->d[i] = d1->d[i]+d2->d[i];
1116: }
1.11 noro 1117: }
1118:
1119: /* d1 += d2 */
1120:
1.19 noro 1121: void adddl_destructive(int n,DL d1,DL d2)
1.11 noro 1122: {
1123: int i;
1124:
1125: d1->td += d2->td;
1126: for ( i = 0; i < n; i++ )
1127: d1->d[i] += d2->d[i];
1.1 noro 1128: }
1129:
1.19 noro 1130: int compd(VL vl,DP p1,DP p2)
1.1 noro 1131: {
1132: int n,t;
1133: MP m1,m2;
1134:
1135: if ( !p1 )
1136: return p2 ? -1 : 0;
1137: else if ( !p2 )
1138: return 1;
1.39 noro 1139: else if ( NV(p1) != NV(p2) )
1140: error("compd : size mismatch");
1.1 noro 1141: else {
1142: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
1143: m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
1144: if ( (t = (*cmpdl)(n,m1->dl,m2->dl)) ||
1145: (t = compp(vl,C(m1),C(m2)) ) )
1146: return t;
1147: if ( m1 )
1148: return 1;
1149: else if ( m2 )
1150: return -1;
1151: else
1152: return 0;
1153: }
1154: }
1155:
1.19 noro 1156: int cmpdl_lex(int n,DL d1,DL d2)
1.1 noro 1157: {
1158: int i;
1159:
1160: for ( i = 0; i < n && d1->d[i] == d2->d[i]; i++ );
1161: return i == n ? 0 : (d1->d[i] > d2->d[i] ? 1 : -1);
1162: }
1163:
1.19 noro 1164: int cmpdl_revlex(int n,DL d1,DL d2)
1.1 noro 1165: {
1166: int i;
1167:
1168: for ( i = n - 1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
1169: return i < 0 ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1170: }
1171:
1.19 noro 1172: int cmpdl_gradlex(int n,DL d1,DL d2)
1.1 noro 1173: {
1174: if ( d1->td > d2->td )
1175: return 1;
1176: else if ( d1->td < d2->td )
1177: return -1;
1178: else
1179: return cmpdl_lex(n,d1,d2);
1180: }
1181:
1.19 noro 1182: int cmpdl_revgradlex(int n,DL d1,DL d2)
1.1 noro 1183: {
1.25 noro 1184: register int i,c;
1.7 noro 1185: register int *p1,*p2;
1186:
1.1 noro 1187: if ( d1->td > d2->td )
1188: return 1;
1189: else if ( d1->td < d2->td )
1190: return -1;
1.7 noro 1191: else {
1.25 noro 1192: i = n-1;
1193: p1 = d1->d+n-1;
1194: p2 = d2->d+n-1;
1195: while ( i >= 7 ) {
1196: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1197: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1198: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1199: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1200: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1201: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1202: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1203: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1204: i -= 8;
1205: }
1206: switch ( i ) {
1207: case 6:
1208: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1209: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1210: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1211: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1212: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1213: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1214: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1215: return 0;
1216: case 5:
1217: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1218: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1219: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1220: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1221: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1222: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1223: return 0;
1224: case 4:
1225: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1226: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1227: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1228: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1229: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1230: return 0;
1231: case 3:
1232: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1233: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1234: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1235: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1236: return 0;
1237: case 2:
1238: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1239: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1240: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1241: return 0;
1242: case 1:
1243: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1244: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1245: return 0;
1246: case 0:
1247: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1248: return 0;
1249: default:
1250: return 0;
1251: }
1252: LAST:
1253: if ( c > 0 ) return -1;
1254: else return 1;
1.7 noro 1255: }
1.1 noro 1256: }
1257:
1.19 noro 1258: int cmpdl_blex(int n,DL d1,DL d2)
1.1 noro 1259: {
1260: int c;
1261:
1262: if ( c = cmpdl_lex(n-1,d1,d2) )
1263: return c;
1264: else {
1265: c = d1->d[n-1] - d2->d[n-1];
1266: return c > 0 ? 1 : c < 0 ? -1 : 0;
1267: }
1268: }
1269:
1.19 noro 1270: int cmpdl_bgradlex(int n,DL d1,DL d2)
1.1 noro 1271: {
1272: int e1,e2,c;
1273:
1274: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
1275: if ( e1 > e2 )
1276: return 1;
1277: else if ( e1 < e2 )
1278: return -1;
1279: else {
1280: c = cmpdl_lex(n-1,d1,d2);
1281: if ( c )
1282: return c;
1283: else
1284: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
1285: }
1286: }
1287:
1.19 noro 1288: int cmpdl_brevgradlex(int n,DL d1,DL d2)
1.1 noro 1289: {
1290: int e1,e2,c;
1291:
1292: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
1293: if ( e1 > e2 )
1294: return 1;
1295: else if ( e1 < e2 )
1296: return -1;
1297: else {
1298: c = cmpdl_revlex(n-1,d1,d2);
1299: if ( c )
1300: return c;
1301: else
1302: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
1303: }
1304: }
1305:
1.19 noro 1306: int cmpdl_brevrev(int n,DL d1,DL d2)
1.1 noro 1307: {
1308: int e1,e2,f1,f2,c,i;
1309:
1310: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1311: e1 += d1->d[i]; e2 += d2->d[i];
1312: }
1313: f1 = d1->td - e1; f2 = d2->td - e2;
1314: if ( e1 > e2 )
1315: return 1;
1316: else if ( e1 < e2 )
1317: return -1;
1318: else {
1319: c = cmpdl_revlex(dp_nelim,d1,d2);
1320: if ( c )
1321: return c;
1322: else if ( f1 > f2 )
1323: return 1;
1324: else if ( f1 < f2 )
1325: return -1;
1326: else {
1327: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1328: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1329: }
1330: }
1331: }
1332:
1.19 noro 1333: int cmpdl_bgradrev(int n,DL d1,DL d2)
1.1 noro 1334: {
1335: int e1,e2,f1,f2,c,i;
1336:
1337: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1338: e1 += d1->d[i]; e2 += d2->d[i];
1339: }
1340: f1 = d1->td - e1; f2 = d2->td - e2;
1341: if ( e1 > e2 )
1342: return 1;
1343: else if ( e1 < e2 )
1344: return -1;
1345: else {
1346: c = cmpdl_lex(dp_nelim,d1,d2);
1347: if ( c )
1348: return c;
1349: else if ( f1 > f2 )
1350: return 1;
1351: else if ( f1 < f2 )
1352: return -1;
1353: else {
1354: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1355: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1356: }
1357: }
1358: }
1359:
1.19 noro 1360: int cmpdl_blexrev(int n,DL d1,DL d2)
1.1 noro 1361: {
1362: int e1,e2,f1,f2,c,i;
1363:
1364: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1365: e1 += d1->d[i]; e2 += d2->d[i];
1366: }
1367: f1 = d1->td - e1; f2 = d2->td - e2;
1368: c = cmpdl_lex(dp_nelim,d1,d2);
1369: if ( c )
1370: return c;
1371: else if ( f1 > f2 )
1372: return 1;
1373: else if ( f1 < f2 )
1374: return -1;
1375: else {
1376: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1377: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1378: }
1379: }
1380:
1.19 noro 1381: int cmpdl_elim(int n,DL d1,DL d2)
1.1 noro 1382: {
1383: int e1,e2,i;
1384:
1385: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1386: e1 += d1->d[i]; e2 += d2->d[i];
1387: }
1388: if ( e1 > e2 )
1389: return 1;
1390: else if ( e1 < e2 )
1391: return -1;
1392: else
1393: return cmpdl_revgradlex(n,d1,d2);
1.12 noro 1394: }
1395:
1.19 noro 1396: int cmpdl_weyl_elim(int n,DL d1,DL d2)
1.12 noro 1397: {
1398: int e1,e2,i;
1399:
1400: for ( i = 1, e1 = 0, e2 = 0; i <= dp_nelim; i++ ) {
1401: e1 += d1->d[n-i]; e2 += d2->d[n-i];
1402: }
1403: if ( e1 > e2 )
1404: return 1;
1405: else if ( e1 < e2 )
1406: return -1;
1407: else if ( d1->td > d2->td )
1408: return 1;
1409: else if ( d1->td < d2->td )
1410: return -1;
1411: else return -cmpdl_revlex(n,d1,d2);
1.13 noro 1412: }
1413:
1414: /*
1415: a special ordering
1416: 1. total order
1417: 2. (-w,w) for the first 2*m variables
1418: 3. DRL for the first 2*m variables
1419: */
1420:
1.20 noro 1421: extern int *current_weyl_weight_vector;
1.13 noro 1422:
1.19 noro 1423: int cmpdl_homo_ww_drl(int n,DL d1,DL d2)
1.13 noro 1424: {
1425: int e1,e2,m,i;
1426: int *p1,*p2;
1427:
1428: if ( d1->td > d2->td )
1429: return 1;
1430: else if ( d1->td < d2->td )
1431: return -1;
1432:
1433: m = n>>1;
1.21 noro 1434: for ( i = 0, e1 = e2 = 0, p1 = d1->d, p2 = d2->d; i < m; i++ ) {
1435: e1 += current_weyl_weight_vector[i]*(p1[m+i] - p1[i]);
1436: e2 += current_weyl_weight_vector[i]*(p2[m+i] - p2[i]);
1.13 noro 1437: }
1438: if ( e1 > e2 )
1439: return 1;
1440: else if ( e1 < e2 )
1441: return -1;
1442:
1443: e1 = d1->td - d1->d[n-1];
1444: e2 = d2->td - d2->d[n-1];
1445: if ( e1 > e2 )
1446: return 1;
1447: else if ( e1 < e2 )
1448: return -1;
1449:
1450: for ( i= n - 1, p1 = d1->d+n-1, p2 = d2->d+n-1;
1451: i >= 0 && *p1 == *p2; i--, p1--, p2-- );
1452: return i < 0 ? 0 : (*p1 < *p2 ? 1 : -1);
1.21 noro 1453: }
1454:
1455: int cmpdl_drl_zigzag(int n,DL d1,DL d2)
1456: {
1457: int i,t,m;
1458: int *p1,*p2;
1459:
1460: if ( d1->td > d2->td )
1461: return 1;
1462: else if ( d1->td < d2->td )
1463: return -1;
1464: else {
1465: m = n>>1;
1466: for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) {
1467: if ( t = p1[m+i] - p2[m+i] ) return t > 0 ? -1 : 1;
1468: if ( t = p1[i] - p2[i] ) return t > 0 ? -1 : 1;
1469: }
1470: return 0;
1471: }
1472: }
1473:
1474: int cmpdl_homo_ww_drl_zigzag(int n,DL d1,DL d2)
1475: {
1476: int e1,e2,m,i,t;
1477: int *p1,*p2;
1478:
1479: if ( d1->td > d2->td )
1480: return 1;
1481: else if ( d1->td < d2->td )
1482: return -1;
1483:
1484: m = n>>1;
1485: for ( i = 0, e1 = e2 = 0, p1 = d1->d, p2 = d2->d; i < m; i++ ) {
1486: e1 += current_weyl_weight_vector[i]*(p1[m+i] - p1[i]);
1487: e2 += current_weyl_weight_vector[i]*(p2[m+i] - p2[i]);
1488: }
1489: if ( e1 > e2 )
1490: return 1;
1491: else if ( e1 < e2 )
1492: return -1;
1493:
1494: e1 = d1->td - d1->d[n-1];
1495: e2 = d2->td - d2->d[n-1];
1496: if ( e1 > e2 )
1497: return 1;
1498: else if ( e1 < e2 )
1499: return -1;
1500:
1501: for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) {
1502: if ( t = p1[m+i] - p2[m+i] ) return t > 0 ? -1 : 1;
1503: if ( t = p1[i] - p2[i] ) return t > 0 ? -1 : 1;
1504: }
1505: return 0;
1.1 noro 1506: }
1507:
1.19 noro 1508: int cmpdl_order_pair(int n,DL d1,DL d2)
1.1 noro 1509: {
1510: int e1,e2,i,j,l;
1511: int *t1,*t2;
1.20 noro 1512: int len,head;
1.1 noro 1513: struct order_pair *pair;
1514:
1.27 noro 1515: len = dp_current_spec->ord.block.length;
1.39 noro 1516: if ( n != dp_current_spec->nv )
1517: error("cmpdl_order_pair : incompatible order specification");
1.27 noro 1518: pair = dp_current_spec->ord.block.order_pair;
1.1 noro 1519:
1.20 noro 1520: head = 0;
1.1 noro 1521: for ( i = 0, t1 = d1->d, t2 = d2->d; i < len; i++ ) {
1522: l = pair[i].length;
1523: switch ( pair[i].order ) {
1524: case 0:
1525: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
1.20 noro 1526: e1 += MUL_WEIGHT(t1[j],head+j);
1527: e2 += MUL_WEIGHT(t2[j],head+j);
1.1 noro 1528: }
1529: if ( e1 > e2 )
1530: return 1;
1531: else if ( e1 < e2 )
1532: return -1;
1533: else {
1534: for ( j = l - 1; j >= 0 && t1[j] == t2[j]; j-- );
1535: if ( j >= 0 )
1536: return t1[j] < t2[j] ? 1 : -1;
1537: }
1538: break;
1539: case 1:
1540: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
1.20 noro 1541: e1 += MUL_WEIGHT(t1[j],head+j);
1542: e2 += MUL_WEIGHT(t2[j],head+j);
1.1 noro 1543: }
1544: if ( e1 > e2 )
1545: return 1;
1546: else if ( e1 < e2 )
1547: return -1;
1548: else {
1549: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
1550: if ( j < l )
1551: return t1[j] > t2[j] ? 1 : -1;
1552: }
1553: break;
1554: case 2:
1555: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
1556: if ( j < l )
1557: return t1[j] > t2[j] ? 1 : -1;
1558: break;
1559: default:
1560: error("cmpdl_order_pair : invalid order"); break;
1561: }
1.20 noro 1562: t1 += l; t2 += l; head += l;
1.28 noro 1563: }
1564: return 0;
1565: }
1566:
1567: int cmpdl_composite(int nv,DL d1,DL d2)
1568: {
1569: int n,i,j,k,start,s,len;
1570: int *dw;
1571: struct sparse_weight *sw;
1572: struct weight_or_block *worb;
1573: int *w,*t1,*t2;
1574:
1575: n = dp_current_spec->ord.composite.length;
1576: worb = dp_current_spec->ord.composite.w_or_b;
1577: w = dp_dl_work;
1578: for ( i = 0, t1 = d1->d, t2 = d2->d; i < nv; i++ )
1579: w[i] = t1[i]-t2[i];
1580: for ( i = 0; i < n; i++, worb++ ) {
1581: len = worb->length;
1582: switch ( worb->type ) {
1583: case IS_DENSE_WEIGHT:
1584: dw = worb->body.dense_weight;
1585: for ( j = 0, s = 0; j < len; j++ )
1586: s += dw[j]*w[j];
1587: if ( s > 0 ) return 1;
1588: else if ( s < 0 ) return -1;
1589: break;
1590: case IS_SPARSE_WEIGHT:
1591: sw = worb->body.sparse_weight;
1592: for ( j = 0, s = 0; j < len; j++ )
1593: s += sw[j].value*w[sw[j].pos];
1594: if ( s > 0 ) return 1;
1595: else if ( s < 0 ) return -1;
1596: break;
1597: case IS_BLOCK:
1598: start = worb->body.block.start;
1599: switch ( worb->body.block.order ) {
1600: case 0:
1601: for ( j = 0, k = start, s = 0; j < len; j++, k++ ) {
1602: s += MUL_WEIGHT(w[k],k);
1603: }
1604: if ( s > 0 ) return 1;
1605: else if ( s < 0 ) return -1;
1606: else {
1607: for ( j = k-1; j >= start && w[j] == 0; j-- );
1608: if ( j >= start )
1609: return w[j] < 0 ? 1 : -1;
1610: }
1611: break;
1612: case 1:
1613: for ( j = 0, k = start, s = 0; j < len; j++, k++ ) {
1614: s += MUL_WEIGHT(w[k],k);
1615: }
1616: if ( s > 0 ) return 1;
1617: else if ( s < 0 ) return -1;
1618: else {
1619: for ( j = 0, k = start; j < len && w[j] == 0; j++, k++ );
1620: if ( j < len )
1621: return w[j] > 0 ? 1 : -1;
1622: }
1623: break;
1624: case 2:
1625: for ( j = 0, k = start; j < len && w[j] == 0; j++, k++ );
1626: if ( j < len )
1627: return w[j] > 0 ? 1 : -1;
1628: break;
1629: }
1630: break;
1631: }
1.1 noro 1632: }
1633: return 0;
1634: }
1635:
1.19 noro 1636: int cmpdl_matrix(int n,DL d1,DL d2)
1.1 noro 1637: {
1638: int *v,*w,*t1,*t2;
1639: int s,i,j,len;
1640: int **matrix;
1641:
1642: for ( i = 0, t1 = d1->d, t2 = d2->d, w = dp_dl_work; i < n; i++ )
1643: w[i] = t1[i]-t2[i];
1.27 noro 1644: len = dp_current_spec->ord.matrix.row;
1645: matrix = dp_current_spec->ord.matrix.matrix;
1.1 noro 1646: for ( j = 0; j < len; j++ ) {
1647: v = matrix[j];
1648: for ( i = 0, s = 0; i < n; i++ )
1649: s += v[i]*w[i];
1650: if ( s > 0 )
1651: return 1;
1652: else if ( s < 0 )
1653: return -1;
1654: }
1655: return 0;
1.25 noro 1656: }
1657:
1.43 ! noro 1658: int cmpdl_top_weight(int n,DL d1,DL d2)
! 1659: {
! 1660: int *w;
! 1661: N sum,wm,wma,t;
! 1662: N *a;
! 1663: struct oN tn;
! 1664: int len,i,sgn,tsgn;
! 1665: int *t1,*t2;
! 1666:
! 1667: w = (int *)ALLOCA(n*sizeof(int));
! 1668: len = current_top_weight_len+3;
! 1669: t1 = d1->d; t2 = d2->d;
! 1670: for ( i = 0; i < n; i++ ) w[i] = t1[i]-t2[i];
! 1671: sum = (N)W_ALLOC(len); sgn = 0;
! 1672: wm = (N)W_ALLOC(len);
! 1673: wma = (N)W_ALLOC(len);
! 1674: a = current_top_weight_vector;
! 1675: for ( i = 0; i < n; i++ ) {
! 1676: if ( !a[i] || !w[i] ) continue;
! 1677: tn.p = 1;
! 1678: if ( w[i] > 0 ) {
! 1679: tn.b[0] = w[i]; tsgn = 1;
! 1680: } else {
! 1681: tn.b[0] = -w[i]; tsgn = -1;
! 1682: }
! 1683: _muln(a[i],&tn,wm);
! 1684: if ( !sgn ) {
! 1685: sgn = tsgn;
! 1686: t = wm; wm = sum; sum = t;
! 1687: } else if ( sgn == tsgn ) {
! 1688: _addn(sum,wm,wma);
! 1689: if ( !PL(wma) )
! 1690: sgn = 0;
! 1691: t = wma; wma = sum; sum = t;
! 1692: } else {
! 1693: sgn *= _subn(sum,wm,wma);
! 1694: t = wma; wma = sum; sum = t;
! 1695: }
! 1696: }
! 1697: if ( sgn > 0 ) return 1;
! 1698: else if ( sgn < 0 ) return -1;
! 1699: else return (*cmpdl_tie_breaker)(n,d1,d2);
! 1700: }
! 1701:
1.25 noro 1702: GeoBucket create_bucket()
1703: {
1704: GeoBucket g;
1705:
1706: g = CALLOC(1,sizeof(struct oGeoBucket));
1707: g->m = 32;
1708: return g;
1709: }
1710:
1711: void add_bucket(GeoBucket g,NODE d,int nv)
1712: {
1713: int l,k,m;
1714:
1715: l = length(d);
1716: for ( k = 0, m = 1; l > m; k++, m <<= 1 );
1717: /* 2^(k-1) < l <= 2^k */
1718: d = symb_merge(g->body[k],d,nv);
1719: for ( ; length(d) > (1<<(k)); k++ ) {
1720: g->body[k] = 0;
1721: d = symb_merge(g->body[k+1],d,nv);
1722: }
1723: g->body[k] = d;
1724: g->m = MAX(g->m,k);
1725: }
1726:
1727: DL remove_head_bucket(GeoBucket g,int nv)
1728: {
1729: int j,i,c,m;
1730: DL d;
1731:
1732: j = -1;
1733: m = g->m;
1734: for ( i = 0; i <= m; i++ ) {
1735: if ( !g->body[i] )
1736: continue;
1737: if ( j < 0 ) j = i;
1738: else {
1739: c = (*cmpdl)(nv,g->body[i]->body,g->body[j]->body);
1740: if ( c > 0 )
1741: j = i;
1742: else if ( c == 0 )
1743: g->body[i] = NEXT(g->body[i]);
1744: }
1745: }
1746: if ( j < 0 )
1747: return 0;
1748: else {
1749: d = g->body[j]->body;
1750: g->body[j] = NEXT(g->body[j]);
1751: return d;
1.31 noro 1752: }
1753: }
1754:
1755: /* DPV functions */
1756:
1757: void adddv(VL vl,DPV p1,DPV p2,DPV *pr)
1758: {
1759: int i,len;
1760: DP *e;
1761:
1762: if ( !p1 || !p2 )
1763: error("adddv : invalid argument");
1764: else if ( p1->len != p2->len )
1765: error("adddv : size mismatch");
1766: else {
1767: len = p1->len;
1768: e = (DP *)MALLOC(p1->len*sizeof(DP));
1769: for ( i = 0; i < len; i++ )
1770: addd(vl,p1->body[i],p2->body[i],&e[i]);
1771: MKDPV(len,e,*pr);
1772: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
1773: }
1774: }
1775:
1776: void subdv(VL vl,DPV p1,DPV p2,DPV *pr)
1777: {
1778: int i,len;
1779: DP *e;
1780:
1781: if ( !p1 || !p2 )
1782: error("subdv : invalid argument");
1783: else if ( p1->len != p2->len )
1784: error("subdv : size mismatch");
1785: else {
1786: len = p1->len;
1787: e = (DP *)MALLOC(p1->len*sizeof(DP));
1788: for ( i = 0; i < len; i++ )
1789: subd(vl,p1->body[i],p2->body[i],&e[i]);
1790: MKDPV(len,e,*pr);
1791: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
1792: }
1793: }
1794:
1795: void chsgndv(DPV p1,DPV *pr)
1796: {
1797: int i,len;
1798: DP *e;
1799:
1800: if ( !p1 )
1801: error("subdv : invalid argument");
1802: else {
1803: len = p1->len;
1804: e = (DP *)MALLOC(p1->len*sizeof(DP));
1805: for ( i = 0; i < len; i++ )
1806: chsgnd(p1->body[i],&e[i]);
1807: MKDPV(len,e,*pr);
1808: (*pr)->sugar = p1->sugar;
1809: }
1810: }
1811:
1812: void muldv(VL vl,DP p1,DPV p2,DPV *pr)
1813: {
1814: int i,len;
1815: DP *e;
1816:
1817: len = p2->len;
1818: e = (DP *)MALLOC(p2->len*sizeof(DP));
1819: if ( !p1 ) {
1820: MKDPV(len,e,*pr);
1821: (*pr)->sugar = 0;
1822: } else {
1823: for ( i = 0; i < len; i++ )
1824: muld(vl,p1,p2->body[i],&e[i]);
1825: MKDPV(len,e,*pr);
1826: (*pr)->sugar = p1->sugar + p2->sugar;
1827: }
1828: }
1829:
1830: int compdv(VL vl,DPV p1,DPV p2)
1831: {
1832: int i,t,len;
1833:
1834: if ( p1->len != p2->len )
1835: error("compdv : size mismatch");
1836: else {
1837: len = p1->len;
1838: for ( i = 0; i < len; i++ )
1839: if ( t = compd(vl,p1->body[i],p2->body[i]) )
1840: return t;
1841: return 0;
1.33 noro 1842: }
1843: }
1844:
1845: int ni_next(int *a,int n)
1846: {
1847: int i,j,k,kj;
1848:
1849: /* find the first nonzero a[j] */
1.35 noro 1850: for ( j = 0; j < n && a[j] == 0; j++ );
1.33 noro 1851: /* find the first zero a[k] after a[j] */
1852: for ( k = j; k < n && a[k] == 1; k++ );
1853: if ( k == n ) return 0;
1854: /* a[0] = 0, ... , a[j-1] = 0, a[j] = 1, ..., a[k-1] = 1, a[k] = 0 */
1855: /* a[0] = 1,..., a[k-j-2] = 1, a[k-j-1] = 0, ..., a[k-1] = 0, a[k] = 1 */
1856: kj = k-j-1;
1857: for ( i = 0; i < kj; i++ ) a[i] = 1;
1858: for ( ; i < k; i++ ) a[i] = 0;
1859: a[k] = 1;
1860: return 1;
1861: }
1862:
1863: int comp_nbm(NBM a,NBM b)
1864: {
1.41 noro 1865: int d,i,w,ai,bi;
1.33 noro 1866: int *ab,*bb;
1867:
1868: if ( a->d > b->d ) return 1;
1869: else if ( a->d < b->d ) return -1;
1870: else {
1871: d = a->d; ab = a->b; bb = b->b;
1.41 noro 1872: #if 0
1.33 noro 1873: w = (d+31)/32;
1874: for ( i = 0; i < w; i++ )
1875: if ( ab[i] > bb[i] ) return 1;
1876: else if ( ab[i] < bb[i] ) return -1;
1.41 noro 1877: #else
1878: for ( i = 0; i < d; i++ ) {
1879: ai = NBM_GET(ab,i);
1880: bi = NBM_GET(bb,i);
1881: if ( ai > bi ) return 1;
1882: else if ( ai < bi ) return -1;
1883: }
1884: #endif
1.33 noro 1885: return 0;
1886: }
1887: }
1888:
1889: NBM mul_nbm(NBM a,NBM b)
1890: {
1891: int ad,bd,d,i,j;
1892: int *ab,*bb,*mb;
1893: NBM m;
1894: NODE r;
1895: NBP u;
1896:
1897: ad = a->d; bd = b->d; ab = a->b; bb = b->b;
1898: d = ad + bd;
1899: NEWNBM(m); NEWNBMBDY(m,d);
1.40 noro 1900: m->d = d; mulp(CO,a->c,b->c,&m->c); mb = m->b;
1.33 noro 1901: j = 0;
1902: for ( i = 0; i < ad; i++, j++ )
1903: if ( NBM_GET(ab,i) ) NBM_SET(mb,j);
1904: else NBM_CLR(mb,j);
1905: for ( i = 0; i < bd; i++, j++ )
1906: if ( NBM_GET(bb,i) ) NBM_SET(mb,j);
1907: else NBM_CLR(mb,j);
1908: return m;
1909: }
1910:
1.37 noro 1911: NBP nbmtonbp(NBM m)
1912: {
1913: NODE n;
1914: NBP u;
1915:
1916: MKNODE(n,m,0);
1917: MKNBP(u,n);
1918: return u;
1919: }
1920:
1921: /* a=c*x*rest -> a0= x*rest, ah=x, ar=rest */
1922:
1.40 noro 1923: P separate_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar)
1.37 noro 1924: {
1925: int i,d1;
1926: NBM t;
1927:
1928: if ( !a->d ) error("separate_nbm : invalid argument");
1929:
1.38 noro 1930: if ( a0 ) {
1.40 noro 1931: NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
1.38 noro 1932: *a0 = nbmtonbp(t);
1933: }
1934:
1935: if ( ah ) {
1.40 noro 1936: NEWNBM(t); NEWNBMBDY(t,1); t->d = 1; t->c = (P)ONE;
1.38 noro 1937: if ( NBM_GET(a->b,0) ) NBM_SET(t->b,0);
1938: else NBM_CLR(t->b,0);
1939: *ah = nbmtonbp(t);
1940: }
1941:
1942: if ( ar ) {
1943: d1 = a->d-1;
1.40 noro 1944: NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
1.38 noro 1945: for ( i = 0; i < d1; i++ ) {
1946: if ( NBM_GET(a->b,i+1) ) NBM_SET(t->b,i);
1.42 noro 1947: else NBM_CLR(t->b,i);
1948: }
1949: *ar = nbmtonbp(t);
1950: }
1951:
1952: return a->c;
1953: }
1954:
1955: /* a=c*rest*x -> a0= rest*x, ar=rest, at=x */
1956:
1957: P separate_tail_nbm(NBM a,NBP *a0,NBP *ar,NBP *at)
1958: {
1959: int i,d,d1;
1960: NBM t;
1961:
1962: if ( !(d=a->d) ) error("separate_tail_nbm : invalid argument");
1963:
1964: if ( a0 ) {
1965: NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
1966: *a0 = nbmtonbp(t);
1967: }
1968:
1969: d1 = a->d-1;
1970: if ( at ) {
1971: NEWNBM(t); NEWNBMBDY(t,1); t->d = 1; t->c = (P)ONE;
1972: if ( NBM_GET(a->b,d1) ) NBM_SET(t->b,0);
1973: else NBM_CLR(t->b,0);
1974: *at = nbmtonbp(t);
1975: }
1976:
1977: if ( ar ) {
1978: NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
1979: for ( i = 0; i < d1; i++ ) {
1980: if ( NBM_GET(a->b,i) ) NBM_SET(t->b,i);
1.38 noro 1981: else NBM_CLR(t->b,i);
1982: }
1983: *ar = nbmtonbp(t);
1984: }
1.37 noro 1985:
1986: return a->c;
1987: }
1988:
1989: NBP make_xky(int k)
1990: {
1991: int k1,i;
1992: NBM t;
1993:
1.40 noro 1994: NEWNBM(t); NEWNBMBDY(t,k); t->d = k; t->c = (P)ONE;
1.37 noro 1995: k1 = k-1;
1996: for ( i = 0; i < k1; i++ ) NBM_SET(t->b,i);
1997: NBM_CLR(t->b,i);
1998: return nbmtonbp(t);
1999: }
2000:
2001: /* a=c*x^(k-1)*y*rest -> a0= x^(k-1)*y*rest, ah=x^(k-1)*y, ar=rest */
2002:
1.40 noro 2003: P separate_xky_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar)
1.37 noro 2004: {
2005: int i,d1,k,k1;
2006: NBM t;
2007:
2008: if ( !a->d )
2009: error("separate_nbm : invalid argument");
2010: for ( i = 0; i < a->d && NBM_GET(a->b,i); i++ );
2011: if ( i == a->d )
2012: error("separate_nbm : invalid argument");
2013: k1 = i;
2014: k = i+1;
2015:
1.38 noro 2016: if ( a0 ) {
1.40 noro 2017: NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
1.38 noro 2018: *a0 = nbmtonbp(t);
2019: }
2020:
2021: if ( ah ) {
1.40 noro 2022: NEWNBM(t); NEWNBMBDY(t,k); t->d = k; t->c = (P)ONE;
1.38 noro 2023: for ( i = 0; i < k1; i++ ) NBM_SET(t->b,i);
2024: NBM_CLR(t->b,i);
2025: *ah = nbmtonbp(t);
2026: }
2027:
2028: if ( ar ) {
2029: d1 = a->d-k;
1.40 noro 2030: NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
1.38 noro 2031: for ( i = 0; i < d1; i++ ) {
2032: if ( NBM_GET(a->b,i+k) ) NBM_SET(t->b,i);
2033: else NBM_CLR(t->b,i);
2034: }
2035: *ar = nbmtonbp(t);
1.37 noro 2036: }
2037:
2038: return a->c;
2039: }
1.33 noro 2040:
1.37 noro 2041: void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp);
2042: void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp);
1.38 noro 2043: void mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp);
2044: void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp);
1.37 noro 2045:
2046: NBP shuffle_mul_nbm(NBM a,NBM b)
2047: {
2048: NBP u,a0,ah,ar,b0,bh,br,a1,b1,t;
1.40 noro 2049: P ac,bc,c;
1.37 noro 2050:
2051: if ( !a->d || !b->d )
2052: u = nbmtonbp(mul_nbm(a,b));
2053: else {
2054: ac = separate_nbm(a,&a0,&ah,&ar);
2055: bc = separate_nbm(b,&b0,&bh,&br);
1.40 noro 2056: mulp(CO,ac,bc,&c);
1.37 noro 2057: shuffle_mulnbp(CO,ar,b0,&t); mulnbp(CO,ah,t,&a1);
2058: shuffle_mulnbp(CO,a0,br,&t); mulnbp(CO,bh,t,&b1);
2059: addnbp(CO,a1,b1,&t); mulnbp(CO,(NBP)c,t,&u);
2060: }
2061: return u;
2062: }
1.33 noro 2063:
1.37 noro 2064: NBP harmonic_mul_nbm(NBM a,NBM b)
2065: {
2066: NBP u,a0,ah,ar,b0,bh,br,a1,b1,t,s,abk,ab1;
1.40 noro 2067: P ac,bc,c;
1.37 noro 2068:
2069: if ( !a->d || !b->d )
2070: u = nbmtonbp(mul_nbm(a,b));
2071: else {
1.40 noro 2072: mulp(CO,a->c,b->c,&c);
1.37 noro 2073: ac = separate_xky_nbm(a,&a0,&ah,&ar);
2074: bc = separate_xky_nbm(b,&b0,&bh,&br);
1.40 noro 2075: mulp(CO,ac,bc,&c);
1.37 noro 2076: harmonic_mulnbp(CO,ar,b0,&t); mulnbp(CO,ah,t,&a1);
2077: harmonic_mulnbp(CO,a0,br,&t); mulnbp(CO,bh,t,&b1);
2078: abk = make_xky(((NBM)BDY(BDY(ah)))->d+((NBM)BDY(BDY(bh)))->d);
2079: harmonic_mulnbp(CO,ar,br,&t); mulnbp(CO,abk,t,&ab1);
2080: addnbp(CO,a1,b1,&t); addnbp(CO,t,ab1,&s); mulnbp(CO,(NBP)c,s,&u);
2081: }
2082: return u;
2083:
2084: }
1.34 noro 2085:
1.33 noro 2086: void addnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2087: {
2088: NODE b1,b2,br,br0;
2089: NBM m1,m2,m;
1.40 noro 2090: P c;
1.33 noro 2091:
2092: if ( !p1 )
2093: *rp = p2;
2094: else if ( !p2 )
2095: *rp = p1;
2096: else {
2097: for ( b1 = BDY(p1), b2 = BDY(p2), br0 = 0; b1 && b2; ) {
2098: m1 = (NBM)BDY(b1); m2 = (NBM)BDY(b2);
2099: switch ( comp_nbm(m1,m2) ) {
2100: case 0:
1.40 noro 2101: addp(CO,m1->c,m2->c,&c);
1.33 noro 2102: if ( c ) {
2103: NEXTNODE(br0,br);
2104: NEWNBM(m); m->d = m1->d; m->c = c; m->b = m1->b;
2105: BDY(br) = (pointer)m;
2106: }
2107: b1 = NEXT(b1); b2 = NEXT(b2); break;
2108: case 1:
2109: NEXTNODE(br0,br); BDY(br) = BDY(b1);
2110: b1 = NEXT(b1); break;
2111: case -1:
2112: NEXTNODE(br0,br); BDY(br) = BDY(b2);
2113: b2 = NEXT(b2); break;
2114: }
1.34 noro 2115: }
2116: if ( !br0 )
2117: if ( b1 )
2118: br0 = b1;
1.33 noro 2119: else if ( b2 )
1.34 noro 2120: br0 = b2;
2121: else {
2122: *rp = 0;
2123: return;
2124: }
2125: else if ( b1 )
2126: NEXT(br) = b1;
2127: else if ( b2 )
1.33 noro 2128: NEXT(br) = b2;
1.34 noro 2129: else
2130: NEXT(br) = 0;
2131: MKNBP(*rp,br0);
1.33 noro 2132: }
2133: }
2134:
2135: void subnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2136: {
2137: NBP t;
2138:
2139: chsgnnbp(p2,&t);
2140: addnbp(vl,p1,t,rp);
2141: }
2142:
2143: void chsgnnbp(NBP p,NBP *rp)
2144: {
2145: NODE r0,r,b;
2146: NBM m,m1;
2147:
2148: for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
2149: NEXTNODE(r0,r);
2150: m = (NBM)BDY(b);
1.40 noro 2151: NEWNBM(m1); m1->d = m->d; m1->b = m->b; chsgnp(m->c,&m1->c);
1.34 noro 2152: BDY(r) = m1;
1.33 noro 2153: }
2154: if ( r0 ) NEXT(r) = 0;
2155: MKNBP(*rp,r0);
2156: }
2157:
2158: void mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2159: {
1.36 noro 2160: NODE b,n;
1.33 noro 2161: NBP r,t,s;
1.36 noro 2162: NBM m;
1.33 noro 2163:
1.36 noro 2164: if ( !p1 || !p2 ) {
2165: *rp = 0; return;
2166: }
2167: if ( OID(p1) != O_NBP ) {
1.40 noro 2168: if ( !POLY(p1) )
1.37 noro 2169: error("mulnbp : invalid argument");
1.40 noro 2170: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
1.36 noro 2171: MKNODE(n,m,0); MKNBP(p1,n);
2172: }
2173: if ( OID(p2) != O_NBP ) {
1.40 noro 2174: if ( !POLY(p2) )
1.37 noro 2175: error("mulnbp : invalid argument");
1.40 noro 2176: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
1.36 noro 2177: MKNODE(n,m,0); MKNBP(p2,n);
2178: }
2179: if ( length(BDY(p1)) < length(BDY(p2)) ) {
1.33 noro 2180: for ( r = 0, b = BDY(p1); b; b = NEXT(b) ) {
2181: mulnbmnbp(vl,(NBM)BDY(b),p2,&t);
2182: addnbp(vl,r,t,&s); r = s;
2183: }
2184: *rp = r;
2185: } else {
2186: for ( r = 0, b = BDY(p2); b; b = NEXT(b) ) {
2187: mulnbpnbm(vl,p1,(NBM)BDY(b),&t);
2188: addnbp(vl,r,t,&s); r = s;
2189: }
2190: *rp = r;
2191: }
2192: }
2193:
2194: void mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp)
2195: {
2196: NODE b,r0,r;
2197:
2198: if ( !p ) *rp = 0;
2199: else {
2200: for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
2201: NEXTNODE(r0,r);
2202: BDY(r) = mul_nbm(m,(NBM)BDY(b));
2203: }
2204: if ( r0 ) NEXT(r) = 0;
2205: MKNBP(*rp,r0);
2206: }
2207: }
2208:
2209: void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp)
2210: {
2211: NODE b,r0,r;
2212:
2213: if ( !p ) *rp = 0;
2214: else {
2215: for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
2216: NEXTNODE(r0,r);
2217: BDY(r) = mul_nbm((NBM)BDY(b),m);
2218: }
2219: if ( r0 ) NEXT(r) = 0;
2220: MKNBP(*rp,r0);
2221: }
2222: }
2223:
2224: void pwrnbp(VL vl,NBP a,Q q,NBP *c)
2225: {
2226: int t;
2227: NBP a1,a2;
2228: N n1;
2229: Q q1;
2230: NBM m;
2231: NODE r;
2232:
2233: if ( !q ) {
1.40 noro 2234: NEWNBM(m); m->d = 0; m->c = (P)ONE; m->b = 0;
1.33 noro 2235: MKNODE(r,m,0); MKNBP(*c,r);
2236: } else if ( !a )
2237: *c = 0;
2238: else if ( UNIQ(q) )
2239: *c = a;
2240: else {
2241: t = divin(NM(q),2,&n1); NTOQ(n1,1,q1);
2242: pwrnbp(vl,a,q1,&a1);
2243: mulnbp(vl,a1,a1,&a2);
2244: if ( t )
2245: mulnbp(vl,a2,a,c);
2246: else
2247: *c = a2;
2248: }
2249: }
2250:
1.38 noro 2251: int compnbp(VL vl,NBP p1,NBP p2)
2252: {
2253: NODE n1,n2;
2254: NBM m1,m2;
2255: int t;
2256:
2257: if ( !p1 )
2258: return p2 ? -1 : 0;
2259: else if ( !p2 )
2260: return 1;
2261: else {
2262: for ( n1 = BDY(p1), n2 = BDY(p2);
2263: n1 && n2; n1 = NEXT(n1), n2 = NEXT(n2) ) {
2264: m1 = (NBM)BDY(n1); m2 = (NBM)BDY(n2);
1.40 noro 2265: if ( (t = comp_nbm(m1,m2)) || (t = compp(CO,m1->c,m2->c) ) )
1.38 noro 2266: return t;
2267: }
2268: if ( n1 )
2269: return 1;
2270: else if ( n2 )
2271: return -1;
2272: else
2273: return 0;
2274: }
2275: }
2276:
1.33 noro 2277: void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2278: {
1.36 noro 2279: NODE b1,b2,n;
1.33 noro 2280: NBP r,t,s;
1.34 noro 2281: NBM m;
1.33 noro 2282:
1.36 noro 2283: if ( !p1 || !p2 ) {
2284: *rp = 0; return;
1.33 noro 2285: }
1.36 noro 2286: if ( OID(p1) != O_NBP ) {
1.40 noro 2287: if ( !POLY(p1) )
1.37 noro 2288: error("shuffle_mulnbp : invalid argument");
1.40 noro 2289: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
1.36 noro 2290: MKNODE(n,m,0); MKNBP(p1,n);
2291: }
2292: if ( OID(p2) != O_NBP ) {
1.40 noro 2293: if ( !POLY(p2) )
1.37 noro 2294: error("shuffle_mulnbp : invalid argument");
1.40 noro 2295: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
1.36 noro 2296: MKNODE(n,m,0); MKNBP(p2,n);
2297: }
2298: for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) )
2299: for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) {
2300: t = shuffle_mul_nbm(m,(NBM)BDY(b2));
2301: addnbp(vl,r,t,&s); r = s;
2302: }
2303: *rp = r;
1.33 noro 2304: }
2305:
1.34 noro 2306: void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
1.33 noro 2307: {
1.36 noro 2308: NODE b1,b2,n;
1.34 noro 2309: NBP r,t,s;
2310: NBM m;
1.33 noro 2311:
1.36 noro 2312: if ( !p1 || !p2 ) {
2313: *rp = 0; return;
1.25 noro 2314: }
1.36 noro 2315: if ( OID(p1) != O_NBP ) {
1.40 noro 2316: if ( !POLY(p1) )
1.37 noro 2317: error("harmonic_mulnbp : invalid argument");
1.40 noro 2318: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
1.36 noro 2319: MKNODE(n,m,0); MKNBP(p1,n);
2320: }
2321: if ( OID(p2) != O_NBP ) {
1.40 noro 2322: if ( !POLY(p2) )
1.37 noro 2323: error("harmonic_mulnbp : invalid argument");
1.40 noro 2324: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
1.36 noro 2325: MKNODE(n,m,0); MKNBP(p2,n);
2326: }
2327: for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) )
2328: for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) {
2329: t = harmonic_mul_nbm(m,(NBM)BDY(b2));
2330: addnbp(vl,r,t,&s); r = s;
2331: }
2332: *rp = r;
1.1 noro 2333: }
1.38 noro 2334:
2335: #if 0
2336: NBP shuffle_mul_nbm(NBM a,NBM b)
2337: {
2338: int ad,bd,d,i,ai,bi,bit,s;
2339: int *ab,*bb,*wmb,*w;
2340: NBM wm,tm;
1.40 noro 2341: P c,c1;
1.38 noro 2342: NODE r,t,t1,p;
2343: NBP u;
2344:
2345: ad = a->d; bd = b->d; ab = a->b; bb = b->b;
2346: d = ad + bd;
2347: w = (int *)ALLOCA(d*sizeof(int));
2348: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2349: for ( i = 0; i < ad; i++ ) w[i] = 1;
2350: for ( ; i < d; i++ ) w[i] = 0;
1.40 noro 2351: mulp(CO,a->c,b->c,&c);
1.38 noro 2352: r = 0;
2353: do {
2354: wm->d = d; wm->c = c;
2355: ai = 0; bi = 0;
2356: for ( i = 0; i < d; i++ ) {
2357: if ( w[i] ) { bit = NBM_GET(ab,ai); ai++; }
2358: else { bit = NBM_GET(bb,bi); bi++; }
2359: if ( bit ) NBM_SET(wmb,i);
2360: else NBM_CLR(wmb,i);
2361: }
2362: for ( p = 0, t = r; t; p = t, t = NEXT(t) ) {
2363: tm = (NBM)BDY(t);
2364: s = comp_nbm(tm,wm);
2365: if ( s < 0 ) {
2366: /* insert */
2367: MKNODE(t1,wm,t);
2368: if ( !p ) r = t1;
2369: else NEXT(p) = t1;
2370: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2371: break;
2372: } else if ( s == 0 ) {
2373: /* add coefs */
1.40 noro 2374: addp(CO,tm->c,c,&c1);
1.38 noro 2375: if ( c1 ) tm->c = c1;
2376: else NEXT(p) = NEXT(t);
2377: break;
2378: }
2379: }
2380: if ( !t ) {
2381: /* append */
2382: MKNODE(t1,wm,t);
2383: if ( !p ) r = t1;
2384: else NEXT(p) = t1;
2385: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2386: }
2387: } while ( ni_next(w,d) );
2388: MKNBP(u,r);
2389: return u;
2390: }
2391:
2392: int nbmtoxky(NBM a,int *b)
2393: {
2394: int d,i,j,k;
2395: int *p;
2396:
2397: d = a->d; p = a->b;
2398: for ( i = j = 0, k = 1; i < d; i++ ) {
2399: if ( !NBM_GET(p,i) ) {
2400: b[j++] = k;
2401: k = 1;
2402: } else k++;
2403: }
2404: return j;
2405: }
2406:
2407: NBP harmonic_mul_nbm(NBM a,NBM b)
2408: {
2409: int da,db,d,la,lb,lmax,lmin,l,lab,la1,lb1,lab1;
2410: int i,j,k,ia,ib,s;
2411: int *wa,*wb,*w,*wab,*wa1,*wmb;
1.40 noro 2412: P c,c1;
1.38 noro 2413: NBM wm,tm;
2414: NODE r,t1,t,p;
2415: NBP u;
2416:
2417: da = a->d; db = b->d; d = da+db;
2418: wa = (int *)ALLOCA(da*sizeof(int));
2419: wb = (int *)ALLOCA(db*sizeof(int));
2420: la = nbmtoxky(a,wa);
2421: lb = nbmtoxky(b,wb);
1.40 noro 2422: mulp(CO,a->c,b->c,&c);
1.38 noro 2423: /* wa[0],..,wa[la-1] <-> x^wa[0]y x^wa[1]y .. */
2424: /* lmax : total length */
2425: lmax = la+lb;
2426: lmin = la>lb?la:lb;
2427: w = (int *)ALLOCA(lmax*sizeof(int));
2428: /* position of a+b */
2429: wab = (int *)ALLOCA(lmax*sizeof(int));
2430: /* position of a */
2431: wa1 = (int *)ALLOCA(lmax*sizeof(int));
2432: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2433: for ( l = lmin, r = 0; l <= lmax; l++ ) {
2434: lab = lmax - l;
2435: la1 = la - lab;
2436: lb1 = lb - lab;
2437: lab1 = l-lab;
2438: /* partion l into three parts: a, b, a+b */
2439: /* initialize wab */
2440: for ( i = 0; i < lab; i++ ) wab[i] = 1;
2441: for ( ; i < l; i++ ) wab[i] = 0;
2442: do {
2443: /* initialize wa1 */
2444: for ( i = 0; i < la1; i++ ) wa1[i] = 1;
2445: for ( ; i < lab1; i++ ) wa1[i] = 0;
2446: do {
2447: ia = 0; ib = 0;
2448: for ( i = j = 0; i < l; i++ )
2449: if ( wab[i] ) w[i] = wa[ia++]+wb[ib++];
2450: else if ( wa1[j++] ) w[i] = wa[ia++];
2451: else w[i] = wb[ib++];
2452: for ( i = j = 0; i < l; i++ ) {
2453: for ( k = w[i]-1; k > 0; k--, j++ ) NBM_SET(wmb,j);
2454: NBM_CLR(wmb,j); j++;
2455: }
2456: wm->d = j; wm->c = c;
2457: for ( p = 0, t = r; t; p = t, t = NEXT(t) ) {
2458: tm = (NBM)BDY(t);
2459: s = comp_nbm(tm,wm);
2460: if ( s < 0 ) {
2461: /* insert */
2462: MKNODE(t1,wm,t);
2463: if ( !p ) r = t1;
2464: else NEXT(p) = t1;
2465: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2466: break;
2467: } else if ( s == 0 ) {
2468: /* add coefs */
1.40 noro 2469: addp(CO,tm->c,c,&c1);
1.38 noro 2470: if ( c1 ) tm->c = c1;
2471: else NEXT(p) = NEXT(t);
2472: break;
2473: }
2474: }
2475: if ( !t ) {
2476: /* append */
2477: MKNODE(t1,wm,t);
2478: if ( !p ) r = t1;
2479: else NEXT(p) = t1;
2480: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2481: }
2482: } while ( ni_next(wa1,lab1) );
2483: } while ( ni_next(wab,l) );
2484: }
2485: MKNBP(u,r);
2486: return u;
2487: }
2488: #endif
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>