![[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.46 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.45 2012/12/17 07:20:44 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);
1.46 ! noro 527: return top;
1.23 noro 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);
1.46 ! noro 537: return top;
1.23 noro 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 ) {
1.45 noro 660: if ( w ) GCFREE(w);
1.5 noro 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 ) {
1.45 noro 698: if ( w ) GCFREE(w);
1.24 noro 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 ) {
1.45 noro 832: if ( w ) GCFREE(w);
1.5 noro 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 ) {
1.45 noro 864: if ( w ) GCFREE(w);
1.5 noro 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 ) {
1.45 noro 876: if ( tab ) GCFREE(tab);
877: if ( psum ) GCFREE(psum);
1.10 noro 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 ) {
1.45 noro 939: if ( tmptab ) GCFREE(tmptab);
1.10 noro 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 ) {
1.45 noro 965: if ( tab ) GCFREE(tab);
966: if ( ctab ) GCFREE(ctab);
1.10 noro 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:
1.44 noro 1107: *dr = dt = (DL)MALLOC_ATOMIC((n+1)*sizeof(int));
1108: dt->td = d1->td + d2->td;
1109: for ( i = 0; i < n; i++ )
1110: dt->d[i] = d1->d[i]+d2->d[i];
1.11 noro 1111: }
1112:
1113: /* d1 += d2 */
1114:
1.19 noro 1115: void adddl_destructive(int n,DL d1,DL d2)
1.11 noro 1116: {
1117: int i;
1118:
1119: d1->td += d2->td;
1120: for ( i = 0; i < n; i++ )
1121: d1->d[i] += d2->d[i];
1.1 noro 1122: }
1123:
1.19 noro 1124: int compd(VL vl,DP p1,DP p2)
1.1 noro 1125: {
1126: int n,t;
1127: MP m1,m2;
1128:
1129: if ( !p1 )
1130: return p2 ? -1 : 0;
1131: else if ( !p2 )
1132: return 1;
1.39 noro 1133: else if ( NV(p1) != NV(p2) )
1134: error("compd : size mismatch");
1.1 noro 1135: else {
1136: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
1137: m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
1138: if ( (t = (*cmpdl)(n,m1->dl,m2->dl)) ||
1139: (t = compp(vl,C(m1),C(m2)) ) )
1140: return t;
1141: if ( m1 )
1142: return 1;
1143: else if ( m2 )
1144: return -1;
1145: else
1146: return 0;
1147: }
1148: }
1149:
1.19 noro 1150: int cmpdl_lex(int n,DL d1,DL d2)
1.1 noro 1151: {
1152: int i;
1153:
1154: for ( i = 0; i < n && d1->d[i] == d2->d[i]; i++ );
1155: return i == n ? 0 : (d1->d[i] > d2->d[i] ? 1 : -1);
1156: }
1157:
1.19 noro 1158: int cmpdl_revlex(int n,DL d1,DL d2)
1.1 noro 1159: {
1160: int i;
1161:
1162: for ( i = n - 1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
1163: return i < 0 ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1164: }
1165:
1.19 noro 1166: int cmpdl_gradlex(int n,DL d1,DL d2)
1.1 noro 1167: {
1168: if ( d1->td > d2->td )
1169: return 1;
1170: else if ( d1->td < d2->td )
1171: return -1;
1172: else
1173: return cmpdl_lex(n,d1,d2);
1174: }
1175:
1.19 noro 1176: int cmpdl_revgradlex(int n,DL d1,DL d2)
1.1 noro 1177: {
1.25 noro 1178: register int i,c;
1.7 noro 1179: register int *p1,*p2;
1180:
1.1 noro 1181: if ( d1->td > d2->td )
1182: return 1;
1183: else if ( d1->td < d2->td )
1184: return -1;
1.7 noro 1185: else {
1.25 noro 1186: i = n-1;
1187: p1 = d1->d+n-1;
1188: p2 = d2->d+n-1;
1189: while ( i >= 7 ) {
1190: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1191: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1192: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1193: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1194: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1195: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1196: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1197: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1198: i -= 8;
1199: }
1200: switch ( i ) {
1201: case 6:
1202: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1203: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1204: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1205: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1206: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1207: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1208: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1209: return 0;
1210: case 5:
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: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1216: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1217: return 0;
1218: case 4:
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: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1224: return 0;
1225: case 3:
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 2:
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: return 0;
1236: case 1:
1237: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1238: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1239: return 0;
1240: case 0:
1241: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1242: return 0;
1243: default:
1244: return 0;
1245: }
1246: LAST:
1247: if ( c > 0 ) return -1;
1248: else return 1;
1.7 noro 1249: }
1.1 noro 1250: }
1251:
1.19 noro 1252: int cmpdl_blex(int n,DL d1,DL d2)
1.1 noro 1253: {
1254: int c;
1255:
1256: if ( c = cmpdl_lex(n-1,d1,d2) )
1257: return c;
1258: else {
1259: c = d1->d[n-1] - d2->d[n-1];
1260: return c > 0 ? 1 : c < 0 ? -1 : 0;
1261: }
1262: }
1263:
1.19 noro 1264: int cmpdl_bgradlex(int n,DL d1,DL d2)
1.1 noro 1265: {
1266: int e1,e2,c;
1267:
1268: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
1269: if ( e1 > e2 )
1270: return 1;
1271: else if ( e1 < e2 )
1272: return -1;
1273: else {
1274: c = cmpdl_lex(n-1,d1,d2);
1275: if ( c )
1276: return c;
1277: else
1278: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
1279: }
1280: }
1281:
1.19 noro 1282: int cmpdl_brevgradlex(int n,DL d1,DL d2)
1.1 noro 1283: {
1284: int e1,e2,c;
1285:
1286: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
1287: if ( e1 > e2 )
1288: return 1;
1289: else if ( e1 < e2 )
1290: return -1;
1291: else {
1292: c = cmpdl_revlex(n-1,d1,d2);
1293: if ( c )
1294: return c;
1295: else
1296: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
1297: }
1298: }
1299:
1.19 noro 1300: int cmpdl_brevrev(int n,DL d1,DL d2)
1.1 noro 1301: {
1302: int e1,e2,f1,f2,c,i;
1303:
1304: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1305: e1 += d1->d[i]; e2 += d2->d[i];
1306: }
1307: f1 = d1->td - e1; f2 = d2->td - e2;
1308: if ( e1 > e2 )
1309: return 1;
1310: else if ( e1 < e2 )
1311: return -1;
1312: else {
1313: c = cmpdl_revlex(dp_nelim,d1,d2);
1314: if ( c )
1315: return c;
1316: else if ( f1 > f2 )
1317: return 1;
1318: else if ( f1 < f2 )
1319: return -1;
1320: else {
1321: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1322: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1323: }
1324: }
1325: }
1326:
1.19 noro 1327: int cmpdl_bgradrev(int n,DL d1,DL d2)
1.1 noro 1328: {
1329: int e1,e2,f1,f2,c,i;
1330:
1331: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1332: e1 += d1->d[i]; e2 += d2->d[i];
1333: }
1334: f1 = d1->td - e1; f2 = d2->td - e2;
1335: if ( e1 > e2 )
1336: return 1;
1337: else if ( e1 < e2 )
1338: return -1;
1339: else {
1340: c = cmpdl_lex(dp_nelim,d1,d2);
1341: if ( c )
1342: return c;
1343: else if ( f1 > f2 )
1344: return 1;
1345: else if ( f1 < f2 )
1346: return -1;
1347: else {
1348: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1349: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1350: }
1351: }
1352: }
1353:
1.19 noro 1354: int cmpdl_blexrev(int n,DL d1,DL d2)
1.1 noro 1355: {
1356: int e1,e2,f1,f2,c,i;
1357:
1358: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1359: e1 += d1->d[i]; e2 += d2->d[i];
1360: }
1361: f1 = d1->td - e1; f2 = d2->td - e2;
1362: c = cmpdl_lex(dp_nelim,d1,d2);
1363: if ( c )
1364: return c;
1365: else if ( f1 > f2 )
1366: return 1;
1367: else if ( f1 < f2 )
1368: return -1;
1369: else {
1370: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1371: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1372: }
1373: }
1374:
1.19 noro 1375: int cmpdl_elim(int n,DL d1,DL d2)
1.1 noro 1376: {
1377: int e1,e2,i;
1378:
1379: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1380: e1 += d1->d[i]; e2 += d2->d[i];
1381: }
1382: if ( e1 > e2 )
1383: return 1;
1384: else if ( e1 < e2 )
1385: return -1;
1386: else
1387: return cmpdl_revgradlex(n,d1,d2);
1.12 noro 1388: }
1389:
1.19 noro 1390: int cmpdl_weyl_elim(int n,DL d1,DL d2)
1.12 noro 1391: {
1392: int e1,e2,i;
1393:
1394: for ( i = 1, e1 = 0, e2 = 0; i <= dp_nelim; i++ ) {
1395: e1 += d1->d[n-i]; e2 += d2->d[n-i];
1396: }
1397: if ( e1 > e2 )
1398: return 1;
1399: else if ( e1 < e2 )
1400: return -1;
1401: else if ( d1->td > d2->td )
1402: return 1;
1403: else if ( d1->td < d2->td )
1404: return -1;
1405: else return -cmpdl_revlex(n,d1,d2);
1.13 noro 1406: }
1407:
1408: /*
1409: a special ordering
1410: 1. total order
1411: 2. (-w,w) for the first 2*m variables
1412: 3. DRL for the first 2*m variables
1413: */
1414:
1.20 noro 1415: extern int *current_weyl_weight_vector;
1.13 noro 1416:
1.19 noro 1417: int cmpdl_homo_ww_drl(int n,DL d1,DL d2)
1.13 noro 1418: {
1419: int e1,e2,m,i;
1420: int *p1,*p2;
1421:
1422: if ( d1->td > d2->td )
1423: return 1;
1424: else if ( d1->td < d2->td )
1425: return -1;
1426:
1427: m = n>>1;
1.21 noro 1428: for ( i = 0, e1 = e2 = 0, p1 = d1->d, p2 = d2->d; i < m; i++ ) {
1429: e1 += current_weyl_weight_vector[i]*(p1[m+i] - p1[i]);
1430: e2 += current_weyl_weight_vector[i]*(p2[m+i] - p2[i]);
1.13 noro 1431: }
1432: if ( e1 > e2 )
1433: return 1;
1434: else if ( e1 < e2 )
1435: return -1;
1436:
1437: e1 = d1->td - d1->d[n-1];
1438: e2 = d2->td - d2->d[n-1];
1439: if ( e1 > e2 )
1440: return 1;
1441: else if ( e1 < e2 )
1442: return -1;
1443:
1444: for ( i= n - 1, p1 = d1->d+n-1, p2 = d2->d+n-1;
1445: i >= 0 && *p1 == *p2; i--, p1--, p2-- );
1446: return i < 0 ? 0 : (*p1 < *p2 ? 1 : -1);
1.21 noro 1447: }
1448:
1449: int cmpdl_drl_zigzag(int n,DL d1,DL d2)
1450: {
1451: int i,t,m;
1452: int *p1,*p2;
1453:
1454: if ( d1->td > d2->td )
1455: return 1;
1456: else if ( d1->td < d2->td )
1457: return -1;
1458: else {
1459: m = n>>1;
1460: for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) {
1461: if ( t = p1[m+i] - p2[m+i] ) return t > 0 ? -1 : 1;
1462: if ( t = p1[i] - p2[i] ) return t > 0 ? -1 : 1;
1463: }
1464: return 0;
1465: }
1466: }
1467:
1468: int cmpdl_homo_ww_drl_zigzag(int n,DL d1,DL d2)
1469: {
1470: int e1,e2,m,i,t;
1471: int *p1,*p2;
1472:
1473: if ( d1->td > d2->td )
1474: return 1;
1475: else if ( d1->td < d2->td )
1476: return -1;
1477:
1478: m = n>>1;
1479: for ( i = 0, e1 = e2 = 0, p1 = d1->d, p2 = d2->d; i < m; i++ ) {
1480: e1 += current_weyl_weight_vector[i]*(p1[m+i] - p1[i]);
1481: e2 += current_weyl_weight_vector[i]*(p2[m+i] - p2[i]);
1482: }
1483: if ( e1 > e2 )
1484: return 1;
1485: else if ( e1 < e2 )
1486: return -1;
1487:
1488: e1 = d1->td - d1->d[n-1];
1489: e2 = d2->td - d2->d[n-1];
1490: if ( e1 > e2 )
1491: return 1;
1492: else if ( e1 < e2 )
1493: return -1;
1494:
1495: for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) {
1496: if ( t = p1[m+i] - p2[m+i] ) return t > 0 ? -1 : 1;
1497: if ( t = p1[i] - p2[i] ) return t > 0 ? -1 : 1;
1498: }
1499: return 0;
1.1 noro 1500: }
1501:
1.19 noro 1502: int cmpdl_order_pair(int n,DL d1,DL d2)
1.1 noro 1503: {
1504: int e1,e2,i,j,l;
1505: int *t1,*t2;
1.20 noro 1506: int len,head;
1.1 noro 1507: struct order_pair *pair;
1508:
1.27 noro 1509: len = dp_current_spec->ord.block.length;
1.39 noro 1510: if ( n != dp_current_spec->nv )
1511: error("cmpdl_order_pair : incompatible order specification");
1.27 noro 1512: pair = dp_current_spec->ord.block.order_pair;
1.1 noro 1513:
1.20 noro 1514: head = 0;
1.1 noro 1515: for ( i = 0, t1 = d1->d, t2 = d2->d; i < len; i++ ) {
1516: l = pair[i].length;
1517: switch ( pair[i].order ) {
1518: case 0:
1519: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
1.20 noro 1520: e1 += MUL_WEIGHT(t1[j],head+j);
1521: e2 += MUL_WEIGHT(t2[j],head+j);
1.1 noro 1522: }
1523: if ( e1 > e2 )
1524: return 1;
1525: else if ( e1 < e2 )
1526: return -1;
1527: else {
1528: for ( j = l - 1; j >= 0 && t1[j] == t2[j]; j-- );
1529: if ( j >= 0 )
1530: return t1[j] < t2[j] ? 1 : -1;
1531: }
1532: break;
1533: case 1:
1534: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
1.20 noro 1535: e1 += MUL_WEIGHT(t1[j],head+j);
1536: e2 += MUL_WEIGHT(t2[j],head+j);
1.1 noro 1537: }
1538: if ( e1 > e2 )
1539: return 1;
1540: else if ( e1 < e2 )
1541: return -1;
1542: else {
1543: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
1544: if ( j < l )
1545: return t1[j] > t2[j] ? 1 : -1;
1546: }
1547: break;
1548: case 2:
1549: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
1550: if ( j < l )
1551: return t1[j] > t2[j] ? 1 : -1;
1552: break;
1553: default:
1554: error("cmpdl_order_pair : invalid order"); break;
1555: }
1.20 noro 1556: t1 += l; t2 += l; head += l;
1.28 noro 1557: }
1558: return 0;
1559: }
1560:
1561: int cmpdl_composite(int nv,DL d1,DL d2)
1562: {
1563: int n,i,j,k,start,s,len;
1564: int *dw;
1565: struct sparse_weight *sw;
1566: struct weight_or_block *worb;
1567: int *w,*t1,*t2;
1568:
1569: n = dp_current_spec->ord.composite.length;
1570: worb = dp_current_spec->ord.composite.w_or_b;
1571: w = dp_dl_work;
1572: for ( i = 0, t1 = d1->d, t2 = d2->d; i < nv; i++ )
1573: w[i] = t1[i]-t2[i];
1574: for ( i = 0; i < n; i++, worb++ ) {
1575: len = worb->length;
1576: switch ( worb->type ) {
1577: case IS_DENSE_WEIGHT:
1578: dw = worb->body.dense_weight;
1579: for ( j = 0, s = 0; j < len; j++ )
1580: s += dw[j]*w[j];
1581: if ( s > 0 ) return 1;
1582: else if ( s < 0 ) return -1;
1583: break;
1584: case IS_SPARSE_WEIGHT:
1585: sw = worb->body.sparse_weight;
1586: for ( j = 0, s = 0; j < len; j++ )
1587: s += sw[j].value*w[sw[j].pos];
1588: if ( s > 0 ) return 1;
1589: else if ( s < 0 ) return -1;
1590: break;
1591: case IS_BLOCK:
1592: start = worb->body.block.start;
1593: switch ( worb->body.block.order ) {
1594: case 0:
1595: for ( j = 0, k = start, s = 0; j < len; j++, k++ ) {
1596: s += MUL_WEIGHT(w[k],k);
1597: }
1598: if ( s > 0 ) return 1;
1599: else if ( s < 0 ) return -1;
1600: else {
1601: for ( j = k-1; j >= start && w[j] == 0; j-- );
1602: if ( j >= start )
1603: return w[j] < 0 ? 1 : -1;
1604: }
1605: break;
1606: case 1:
1607: for ( j = 0, k = start, s = 0; j < len; j++, k++ ) {
1608: s += MUL_WEIGHT(w[k],k);
1609: }
1610: if ( s > 0 ) return 1;
1611: else if ( s < 0 ) return -1;
1612: else {
1613: for ( j = 0, k = start; j < len && w[j] == 0; j++, k++ );
1614: if ( j < len )
1615: return w[j] > 0 ? 1 : -1;
1616: }
1617: break;
1618: case 2:
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: break;
1623: }
1624: break;
1625: }
1.1 noro 1626: }
1627: return 0;
1628: }
1629:
1.19 noro 1630: int cmpdl_matrix(int n,DL d1,DL d2)
1.1 noro 1631: {
1632: int *v,*w,*t1,*t2;
1633: int s,i,j,len;
1634: int **matrix;
1635:
1636: for ( i = 0, t1 = d1->d, t2 = d2->d, w = dp_dl_work; i < n; i++ )
1637: w[i] = t1[i]-t2[i];
1.27 noro 1638: len = dp_current_spec->ord.matrix.row;
1639: matrix = dp_current_spec->ord.matrix.matrix;
1.1 noro 1640: for ( j = 0; j < len; j++ ) {
1641: v = matrix[j];
1642: for ( i = 0, s = 0; i < n; i++ )
1643: s += v[i]*w[i];
1644: if ( s > 0 )
1645: return 1;
1646: else if ( s < 0 )
1647: return -1;
1648: }
1649: return 0;
1.25 noro 1650: }
1651:
1.43 noro 1652: int cmpdl_top_weight(int n,DL d1,DL d2)
1653: {
1654: int *w;
1655: N sum,wm,wma,t;
1656: N *a;
1657: struct oN tn;
1658: int len,i,sgn,tsgn;
1659: int *t1,*t2;
1660:
1661: w = (int *)ALLOCA(n*sizeof(int));
1662: len = current_top_weight_len+3;
1663: t1 = d1->d; t2 = d2->d;
1664: for ( i = 0; i < n; i++ ) w[i] = t1[i]-t2[i];
1665: sum = (N)W_ALLOC(len); sgn = 0;
1666: wm = (N)W_ALLOC(len);
1667: wma = (N)W_ALLOC(len);
1668: a = current_top_weight_vector;
1669: for ( i = 0; i < n; i++ ) {
1670: if ( !a[i] || !w[i] ) continue;
1671: tn.p = 1;
1672: if ( w[i] > 0 ) {
1673: tn.b[0] = w[i]; tsgn = 1;
1674: } else {
1675: tn.b[0] = -w[i]; tsgn = -1;
1676: }
1677: _muln(a[i],&tn,wm);
1678: if ( !sgn ) {
1679: sgn = tsgn;
1680: t = wm; wm = sum; sum = t;
1681: } else if ( sgn == tsgn ) {
1682: _addn(sum,wm,wma);
1683: if ( !PL(wma) )
1684: sgn = 0;
1685: t = wma; wma = sum; sum = t;
1686: } else {
1687: sgn *= _subn(sum,wm,wma);
1688: t = wma; wma = sum; sum = t;
1689: }
1690: }
1691: if ( sgn > 0 ) return 1;
1692: else if ( sgn < 0 ) return -1;
1693: else return (*cmpdl_tie_breaker)(n,d1,d2);
1694: }
1695:
1.25 noro 1696: GeoBucket create_bucket()
1697: {
1698: GeoBucket g;
1699:
1700: g = CALLOC(1,sizeof(struct oGeoBucket));
1701: g->m = 32;
1702: return g;
1703: }
1704:
1705: void add_bucket(GeoBucket g,NODE d,int nv)
1706: {
1707: int l,k,m;
1708:
1709: l = length(d);
1710: for ( k = 0, m = 1; l > m; k++, m <<= 1 );
1711: /* 2^(k-1) < l <= 2^k */
1712: d = symb_merge(g->body[k],d,nv);
1713: for ( ; length(d) > (1<<(k)); k++ ) {
1714: g->body[k] = 0;
1715: d = symb_merge(g->body[k+1],d,nv);
1716: }
1717: g->body[k] = d;
1718: g->m = MAX(g->m,k);
1719: }
1720:
1721: DL remove_head_bucket(GeoBucket g,int nv)
1722: {
1723: int j,i,c,m;
1724: DL d;
1725:
1726: j = -1;
1727: m = g->m;
1728: for ( i = 0; i <= m; i++ ) {
1729: if ( !g->body[i] )
1730: continue;
1731: if ( j < 0 ) j = i;
1732: else {
1733: c = (*cmpdl)(nv,g->body[i]->body,g->body[j]->body);
1734: if ( c > 0 )
1735: j = i;
1736: else if ( c == 0 )
1737: g->body[i] = NEXT(g->body[i]);
1738: }
1739: }
1740: if ( j < 0 )
1741: return 0;
1742: else {
1743: d = g->body[j]->body;
1744: g->body[j] = NEXT(g->body[j]);
1745: return d;
1.31 noro 1746: }
1747: }
1748:
1749: /* DPV functions */
1750:
1751: void adddv(VL vl,DPV p1,DPV p2,DPV *pr)
1752: {
1753: int i,len;
1754: DP *e;
1755:
1756: if ( !p1 || !p2 )
1757: error("adddv : invalid argument");
1758: else if ( p1->len != p2->len )
1759: error("adddv : size mismatch");
1760: else {
1761: len = p1->len;
1762: e = (DP *)MALLOC(p1->len*sizeof(DP));
1763: for ( i = 0; i < len; i++ )
1764: addd(vl,p1->body[i],p2->body[i],&e[i]);
1765: MKDPV(len,e,*pr);
1766: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
1767: }
1768: }
1769:
1770: void subdv(VL vl,DPV p1,DPV p2,DPV *pr)
1771: {
1772: int i,len;
1773: DP *e;
1774:
1775: if ( !p1 || !p2 )
1776: error("subdv : invalid argument");
1777: else if ( p1->len != p2->len )
1778: error("subdv : size mismatch");
1779: else {
1780: len = p1->len;
1781: e = (DP *)MALLOC(p1->len*sizeof(DP));
1782: for ( i = 0; i < len; i++ )
1783: subd(vl,p1->body[i],p2->body[i],&e[i]);
1784: MKDPV(len,e,*pr);
1785: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
1786: }
1787: }
1788:
1789: void chsgndv(DPV p1,DPV *pr)
1790: {
1791: int i,len;
1792: DP *e;
1793:
1794: if ( !p1 )
1795: error("subdv : invalid argument");
1796: else {
1797: len = p1->len;
1798: e = (DP *)MALLOC(p1->len*sizeof(DP));
1799: for ( i = 0; i < len; i++ )
1800: chsgnd(p1->body[i],&e[i]);
1801: MKDPV(len,e,*pr);
1802: (*pr)->sugar = p1->sugar;
1803: }
1804: }
1805:
1806: void muldv(VL vl,DP p1,DPV p2,DPV *pr)
1807: {
1808: int i,len;
1809: DP *e;
1810:
1811: len = p2->len;
1812: e = (DP *)MALLOC(p2->len*sizeof(DP));
1813: if ( !p1 ) {
1814: MKDPV(len,e,*pr);
1815: (*pr)->sugar = 0;
1816: } else {
1817: for ( i = 0; i < len; i++ )
1818: muld(vl,p1,p2->body[i],&e[i]);
1819: MKDPV(len,e,*pr);
1820: (*pr)->sugar = p1->sugar + p2->sugar;
1821: }
1822: }
1823:
1824: int compdv(VL vl,DPV p1,DPV p2)
1825: {
1826: int i,t,len;
1827:
1828: if ( p1->len != p2->len )
1829: error("compdv : size mismatch");
1830: else {
1831: len = p1->len;
1832: for ( i = 0; i < len; i++ )
1833: if ( t = compd(vl,p1->body[i],p2->body[i]) )
1834: return t;
1835: return 0;
1.33 noro 1836: }
1837: }
1838:
1839: int ni_next(int *a,int n)
1840: {
1841: int i,j,k,kj;
1842:
1843: /* find the first nonzero a[j] */
1.35 noro 1844: for ( j = 0; j < n && a[j] == 0; j++ );
1.33 noro 1845: /* find the first zero a[k] after a[j] */
1846: for ( k = j; k < n && a[k] == 1; k++ );
1847: if ( k == n ) return 0;
1848: /* a[0] = 0, ... , a[j-1] = 0, a[j] = 1, ..., a[k-1] = 1, a[k] = 0 */
1849: /* a[0] = 1,..., a[k-j-2] = 1, a[k-j-1] = 0, ..., a[k-1] = 0, a[k] = 1 */
1850: kj = k-j-1;
1851: for ( i = 0; i < kj; i++ ) a[i] = 1;
1852: for ( ; i < k; i++ ) a[i] = 0;
1853: a[k] = 1;
1854: return 1;
1855: }
1856:
1857: int comp_nbm(NBM a,NBM b)
1858: {
1.41 noro 1859: int d,i,w,ai,bi;
1.33 noro 1860: int *ab,*bb;
1861:
1862: if ( a->d > b->d ) return 1;
1863: else if ( a->d < b->d ) return -1;
1864: else {
1865: d = a->d; ab = a->b; bb = b->b;
1.41 noro 1866: #if 0
1.33 noro 1867: w = (d+31)/32;
1868: for ( i = 0; i < w; i++ )
1869: if ( ab[i] > bb[i] ) return 1;
1870: else if ( ab[i] < bb[i] ) return -1;
1.41 noro 1871: #else
1872: for ( i = 0; i < d; i++ ) {
1873: ai = NBM_GET(ab,i);
1874: bi = NBM_GET(bb,i);
1875: if ( ai > bi ) return 1;
1876: else if ( ai < bi ) return -1;
1877: }
1878: #endif
1.33 noro 1879: return 0;
1880: }
1881: }
1882:
1883: NBM mul_nbm(NBM a,NBM b)
1884: {
1885: int ad,bd,d,i,j;
1886: int *ab,*bb,*mb;
1887: NBM m;
1888: NODE r;
1889: NBP u;
1890:
1891: ad = a->d; bd = b->d; ab = a->b; bb = b->b;
1892: d = ad + bd;
1893: NEWNBM(m); NEWNBMBDY(m,d);
1.40 noro 1894: m->d = d; mulp(CO,a->c,b->c,&m->c); mb = m->b;
1.33 noro 1895: j = 0;
1896: for ( i = 0; i < ad; i++, j++ )
1897: if ( NBM_GET(ab,i) ) NBM_SET(mb,j);
1898: else NBM_CLR(mb,j);
1899: for ( i = 0; i < bd; i++, j++ )
1900: if ( NBM_GET(bb,i) ) NBM_SET(mb,j);
1901: else NBM_CLR(mb,j);
1902: return m;
1903: }
1904:
1.37 noro 1905: NBP nbmtonbp(NBM m)
1906: {
1907: NODE n;
1908: NBP u;
1909:
1910: MKNODE(n,m,0);
1911: MKNBP(u,n);
1912: return u;
1913: }
1914:
1915: /* a=c*x*rest -> a0= x*rest, ah=x, ar=rest */
1916:
1.40 noro 1917: P separate_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar)
1.37 noro 1918: {
1919: int i,d1;
1920: NBM t;
1921:
1922: if ( !a->d ) error("separate_nbm : invalid argument");
1923:
1.38 noro 1924: if ( a0 ) {
1.40 noro 1925: NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
1.38 noro 1926: *a0 = nbmtonbp(t);
1927: }
1928:
1929: if ( ah ) {
1.40 noro 1930: NEWNBM(t); NEWNBMBDY(t,1); t->d = 1; t->c = (P)ONE;
1.38 noro 1931: if ( NBM_GET(a->b,0) ) NBM_SET(t->b,0);
1932: else NBM_CLR(t->b,0);
1933: *ah = nbmtonbp(t);
1934: }
1935:
1936: if ( ar ) {
1937: d1 = a->d-1;
1.40 noro 1938: NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
1.38 noro 1939: for ( i = 0; i < d1; i++ ) {
1940: if ( NBM_GET(a->b,i+1) ) NBM_SET(t->b,i);
1.42 noro 1941: else NBM_CLR(t->b,i);
1942: }
1943: *ar = nbmtonbp(t);
1944: }
1945:
1946: return a->c;
1947: }
1948:
1949: /* a=c*rest*x -> a0= rest*x, ar=rest, at=x */
1950:
1951: P separate_tail_nbm(NBM a,NBP *a0,NBP *ar,NBP *at)
1952: {
1953: int i,d,d1;
1954: NBM t;
1955:
1956: if ( !(d=a->d) ) error("separate_tail_nbm : invalid argument");
1957:
1958: if ( a0 ) {
1959: NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
1960: *a0 = nbmtonbp(t);
1961: }
1962:
1963: d1 = a->d-1;
1964: if ( at ) {
1965: NEWNBM(t); NEWNBMBDY(t,1); t->d = 1; t->c = (P)ONE;
1966: if ( NBM_GET(a->b,d1) ) NBM_SET(t->b,0);
1967: else NBM_CLR(t->b,0);
1968: *at = nbmtonbp(t);
1969: }
1970:
1971: if ( ar ) {
1972: NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
1973: for ( i = 0; i < d1; i++ ) {
1974: if ( NBM_GET(a->b,i) ) NBM_SET(t->b,i);
1.38 noro 1975: else NBM_CLR(t->b,i);
1976: }
1977: *ar = nbmtonbp(t);
1978: }
1.37 noro 1979:
1980: return a->c;
1981: }
1982:
1983: NBP make_xky(int k)
1984: {
1985: int k1,i;
1986: NBM t;
1987:
1.40 noro 1988: NEWNBM(t); NEWNBMBDY(t,k); t->d = k; t->c = (P)ONE;
1.37 noro 1989: k1 = k-1;
1990: for ( i = 0; i < k1; i++ ) NBM_SET(t->b,i);
1991: NBM_CLR(t->b,i);
1992: return nbmtonbp(t);
1993: }
1994:
1995: /* a=c*x^(k-1)*y*rest -> a0= x^(k-1)*y*rest, ah=x^(k-1)*y, ar=rest */
1996:
1.40 noro 1997: P separate_xky_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar)
1.37 noro 1998: {
1999: int i,d1,k,k1;
2000: NBM t;
2001:
2002: if ( !a->d )
2003: error("separate_nbm : invalid argument");
2004: for ( i = 0; i < a->d && NBM_GET(a->b,i); i++ );
2005: if ( i == a->d )
2006: error("separate_nbm : invalid argument");
2007: k1 = i;
2008: k = i+1;
2009:
1.38 noro 2010: if ( a0 ) {
1.40 noro 2011: NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
1.38 noro 2012: *a0 = nbmtonbp(t);
2013: }
2014:
2015: if ( ah ) {
1.40 noro 2016: NEWNBM(t); NEWNBMBDY(t,k); t->d = k; t->c = (P)ONE;
1.38 noro 2017: for ( i = 0; i < k1; i++ ) NBM_SET(t->b,i);
2018: NBM_CLR(t->b,i);
2019: *ah = nbmtonbp(t);
2020: }
2021:
2022: if ( ar ) {
2023: d1 = a->d-k;
1.40 noro 2024: NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
1.38 noro 2025: for ( i = 0; i < d1; i++ ) {
2026: if ( NBM_GET(a->b,i+k) ) NBM_SET(t->b,i);
2027: else NBM_CLR(t->b,i);
2028: }
2029: *ar = nbmtonbp(t);
1.37 noro 2030: }
2031:
2032: return a->c;
2033: }
1.33 noro 2034:
1.37 noro 2035: void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp);
2036: void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp);
1.38 noro 2037: void mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp);
2038: void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp);
1.37 noro 2039:
2040: NBP shuffle_mul_nbm(NBM a,NBM b)
2041: {
2042: NBP u,a0,ah,ar,b0,bh,br,a1,b1,t;
1.40 noro 2043: P ac,bc,c;
1.37 noro 2044:
2045: if ( !a->d || !b->d )
2046: u = nbmtonbp(mul_nbm(a,b));
2047: else {
2048: ac = separate_nbm(a,&a0,&ah,&ar);
2049: bc = separate_nbm(b,&b0,&bh,&br);
1.40 noro 2050: mulp(CO,ac,bc,&c);
1.37 noro 2051: shuffle_mulnbp(CO,ar,b0,&t); mulnbp(CO,ah,t,&a1);
2052: shuffle_mulnbp(CO,a0,br,&t); mulnbp(CO,bh,t,&b1);
2053: addnbp(CO,a1,b1,&t); mulnbp(CO,(NBP)c,t,&u);
2054: }
2055: return u;
2056: }
1.33 noro 2057:
1.37 noro 2058: NBP harmonic_mul_nbm(NBM a,NBM b)
2059: {
2060: NBP u,a0,ah,ar,b0,bh,br,a1,b1,t,s,abk,ab1;
1.40 noro 2061: P ac,bc,c;
1.37 noro 2062:
2063: if ( !a->d || !b->d )
2064: u = nbmtonbp(mul_nbm(a,b));
2065: else {
1.40 noro 2066: mulp(CO,a->c,b->c,&c);
1.37 noro 2067: ac = separate_xky_nbm(a,&a0,&ah,&ar);
2068: bc = separate_xky_nbm(b,&b0,&bh,&br);
1.40 noro 2069: mulp(CO,ac,bc,&c);
1.37 noro 2070: harmonic_mulnbp(CO,ar,b0,&t); mulnbp(CO,ah,t,&a1);
2071: harmonic_mulnbp(CO,a0,br,&t); mulnbp(CO,bh,t,&b1);
2072: abk = make_xky(((NBM)BDY(BDY(ah)))->d+((NBM)BDY(BDY(bh)))->d);
2073: harmonic_mulnbp(CO,ar,br,&t); mulnbp(CO,abk,t,&ab1);
2074: addnbp(CO,a1,b1,&t); addnbp(CO,t,ab1,&s); mulnbp(CO,(NBP)c,s,&u);
2075: }
2076: return u;
2077:
2078: }
1.34 noro 2079:
1.33 noro 2080: void addnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2081: {
2082: NODE b1,b2,br,br0;
2083: NBM m1,m2,m;
1.40 noro 2084: P c;
1.33 noro 2085:
2086: if ( !p1 )
2087: *rp = p2;
2088: else if ( !p2 )
2089: *rp = p1;
2090: else {
2091: for ( b1 = BDY(p1), b2 = BDY(p2), br0 = 0; b1 && b2; ) {
2092: m1 = (NBM)BDY(b1); m2 = (NBM)BDY(b2);
2093: switch ( comp_nbm(m1,m2) ) {
2094: case 0:
1.40 noro 2095: addp(CO,m1->c,m2->c,&c);
1.33 noro 2096: if ( c ) {
2097: NEXTNODE(br0,br);
2098: NEWNBM(m); m->d = m1->d; m->c = c; m->b = m1->b;
2099: BDY(br) = (pointer)m;
2100: }
2101: b1 = NEXT(b1); b2 = NEXT(b2); break;
2102: case 1:
2103: NEXTNODE(br0,br); BDY(br) = BDY(b1);
2104: b1 = NEXT(b1); break;
2105: case -1:
2106: NEXTNODE(br0,br); BDY(br) = BDY(b2);
2107: b2 = NEXT(b2); break;
2108: }
1.34 noro 2109: }
2110: if ( !br0 )
2111: if ( b1 )
2112: br0 = b1;
1.33 noro 2113: else if ( b2 )
1.34 noro 2114: br0 = b2;
2115: else {
2116: *rp = 0;
2117: return;
2118: }
2119: else if ( b1 )
2120: NEXT(br) = b1;
2121: else if ( b2 )
1.33 noro 2122: NEXT(br) = b2;
1.34 noro 2123: else
2124: NEXT(br) = 0;
2125: MKNBP(*rp,br0);
1.33 noro 2126: }
2127: }
2128:
2129: void subnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2130: {
2131: NBP t;
2132:
2133: chsgnnbp(p2,&t);
2134: addnbp(vl,p1,t,rp);
2135: }
2136:
2137: void chsgnnbp(NBP p,NBP *rp)
2138: {
2139: NODE r0,r,b;
2140: NBM m,m1;
2141:
2142: for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
2143: NEXTNODE(r0,r);
2144: m = (NBM)BDY(b);
1.40 noro 2145: NEWNBM(m1); m1->d = m->d; m1->b = m->b; chsgnp(m->c,&m1->c);
1.34 noro 2146: BDY(r) = m1;
1.33 noro 2147: }
2148: if ( r0 ) NEXT(r) = 0;
2149: MKNBP(*rp,r0);
2150: }
2151:
2152: void mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2153: {
1.36 noro 2154: NODE b,n;
1.33 noro 2155: NBP r,t,s;
1.36 noro 2156: NBM m;
1.33 noro 2157:
1.36 noro 2158: if ( !p1 || !p2 ) {
2159: *rp = 0; return;
2160: }
2161: if ( OID(p1) != O_NBP ) {
1.40 noro 2162: if ( !POLY(p1) )
1.37 noro 2163: error("mulnbp : invalid argument");
1.40 noro 2164: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
1.36 noro 2165: MKNODE(n,m,0); MKNBP(p1,n);
2166: }
2167: if ( OID(p2) != O_NBP ) {
1.40 noro 2168: if ( !POLY(p2) )
1.37 noro 2169: error("mulnbp : invalid argument");
1.40 noro 2170: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
1.36 noro 2171: MKNODE(n,m,0); MKNBP(p2,n);
2172: }
2173: if ( length(BDY(p1)) < length(BDY(p2)) ) {
1.33 noro 2174: for ( r = 0, b = BDY(p1); b; b = NEXT(b) ) {
2175: mulnbmnbp(vl,(NBM)BDY(b),p2,&t);
2176: addnbp(vl,r,t,&s); r = s;
2177: }
2178: *rp = r;
2179: } else {
2180: for ( r = 0, b = BDY(p2); b; b = NEXT(b) ) {
2181: mulnbpnbm(vl,p1,(NBM)BDY(b),&t);
2182: addnbp(vl,r,t,&s); r = s;
2183: }
2184: *rp = r;
2185: }
2186: }
2187:
2188: void mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp)
2189: {
2190: NODE b,r0,r;
2191:
2192: if ( !p ) *rp = 0;
2193: else {
2194: for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
2195: NEXTNODE(r0,r);
2196: BDY(r) = mul_nbm(m,(NBM)BDY(b));
2197: }
2198: if ( r0 ) NEXT(r) = 0;
2199: MKNBP(*rp,r0);
2200: }
2201: }
2202:
2203: void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp)
2204: {
2205: NODE b,r0,r;
2206:
2207: if ( !p ) *rp = 0;
2208: else {
2209: for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
2210: NEXTNODE(r0,r);
2211: BDY(r) = mul_nbm((NBM)BDY(b),m);
2212: }
2213: if ( r0 ) NEXT(r) = 0;
2214: MKNBP(*rp,r0);
2215: }
2216: }
2217:
2218: void pwrnbp(VL vl,NBP a,Q q,NBP *c)
2219: {
2220: int t;
2221: NBP a1,a2;
2222: N n1;
2223: Q q1;
2224: NBM m;
2225: NODE r;
2226:
2227: if ( !q ) {
1.40 noro 2228: NEWNBM(m); m->d = 0; m->c = (P)ONE; m->b = 0;
1.33 noro 2229: MKNODE(r,m,0); MKNBP(*c,r);
2230: } else if ( !a )
2231: *c = 0;
2232: else if ( UNIQ(q) )
2233: *c = a;
2234: else {
2235: t = divin(NM(q),2,&n1); NTOQ(n1,1,q1);
2236: pwrnbp(vl,a,q1,&a1);
2237: mulnbp(vl,a1,a1,&a2);
2238: if ( t )
2239: mulnbp(vl,a2,a,c);
2240: else
2241: *c = a2;
2242: }
2243: }
2244:
1.38 noro 2245: int compnbp(VL vl,NBP p1,NBP p2)
2246: {
2247: NODE n1,n2;
2248: NBM m1,m2;
2249: int t;
2250:
2251: if ( !p1 )
2252: return p2 ? -1 : 0;
2253: else if ( !p2 )
2254: return 1;
2255: else {
2256: for ( n1 = BDY(p1), n2 = BDY(p2);
2257: n1 && n2; n1 = NEXT(n1), n2 = NEXT(n2) ) {
2258: m1 = (NBM)BDY(n1); m2 = (NBM)BDY(n2);
1.40 noro 2259: if ( (t = comp_nbm(m1,m2)) || (t = compp(CO,m1->c,m2->c) ) )
1.38 noro 2260: return t;
2261: }
2262: if ( n1 )
2263: return 1;
2264: else if ( n2 )
2265: return -1;
2266: else
2267: return 0;
2268: }
2269: }
2270:
1.33 noro 2271: void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2272: {
1.36 noro 2273: NODE b1,b2,n;
1.33 noro 2274: NBP r,t,s;
1.34 noro 2275: NBM m;
1.33 noro 2276:
1.36 noro 2277: if ( !p1 || !p2 ) {
2278: *rp = 0; return;
1.33 noro 2279: }
1.36 noro 2280: if ( OID(p1) != O_NBP ) {
1.40 noro 2281: if ( !POLY(p1) )
1.37 noro 2282: error("shuffle_mulnbp : invalid argument");
1.40 noro 2283: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
1.36 noro 2284: MKNODE(n,m,0); MKNBP(p1,n);
2285: }
2286: if ( OID(p2) != O_NBP ) {
1.40 noro 2287: if ( !POLY(p2) )
1.37 noro 2288: error("shuffle_mulnbp : invalid argument");
1.40 noro 2289: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
1.36 noro 2290: MKNODE(n,m,0); MKNBP(p2,n);
2291: }
2292: for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) )
2293: for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) {
2294: t = shuffle_mul_nbm(m,(NBM)BDY(b2));
2295: addnbp(vl,r,t,&s); r = s;
2296: }
2297: *rp = r;
1.33 noro 2298: }
2299:
1.34 noro 2300: void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
1.33 noro 2301: {
1.36 noro 2302: NODE b1,b2,n;
1.34 noro 2303: NBP r,t,s;
2304: NBM m;
1.33 noro 2305:
1.36 noro 2306: if ( !p1 || !p2 ) {
2307: *rp = 0; return;
1.25 noro 2308: }
1.36 noro 2309: if ( OID(p1) != O_NBP ) {
1.40 noro 2310: if ( !POLY(p1) )
1.37 noro 2311: error("harmonic_mulnbp : invalid argument");
1.40 noro 2312: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
1.36 noro 2313: MKNODE(n,m,0); MKNBP(p1,n);
2314: }
2315: if ( OID(p2) != O_NBP ) {
1.40 noro 2316: if ( !POLY(p2) )
1.37 noro 2317: error("harmonic_mulnbp : invalid argument");
1.40 noro 2318: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
1.36 noro 2319: MKNODE(n,m,0); MKNBP(p2,n);
2320: }
2321: for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) )
2322: for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) {
2323: t = harmonic_mul_nbm(m,(NBM)BDY(b2));
2324: addnbp(vl,r,t,&s); r = s;
2325: }
2326: *rp = r;
1.1 noro 2327: }
1.38 noro 2328:
2329: #if 0
2330: NBP shuffle_mul_nbm(NBM a,NBM b)
2331: {
2332: int ad,bd,d,i,ai,bi,bit,s;
2333: int *ab,*bb,*wmb,*w;
2334: NBM wm,tm;
1.40 noro 2335: P c,c1;
1.38 noro 2336: NODE r,t,t1,p;
2337: NBP u;
2338:
2339: ad = a->d; bd = b->d; ab = a->b; bb = b->b;
2340: d = ad + bd;
2341: w = (int *)ALLOCA(d*sizeof(int));
2342: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2343: for ( i = 0; i < ad; i++ ) w[i] = 1;
2344: for ( ; i < d; i++ ) w[i] = 0;
1.40 noro 2345: mulp(CO,a->c,b->c,&c);
1.38 noro 2346: r = 0;
2347: do {
2348: wm->d = d; wm->c = c;
2349: ai = 0; bi = 0;
2350: for ( i = 0; i < d; i++ ) {
2351: if ( w[i] ) { bit = NBM_GET(ab,ai); ai++; }
2352: else { bit = NBM_GET(bb,bi); bi++; }
2353: if ( bit ) NBM_SET(wmb,i);
2354: else NBM_CLR(wmb,i);
2355: }
2356: for ( p = 0, t = r; t; p = t, t = NEXT(t) ) {
2357: tm = (NBM)BDY(t);
2358: s = comp_nbm(tm,wm);
2359: if ( s < 0 ) {
2360: /* insert */
2361: MKNODE(t1,wm,t);
2362: if ( !p ) r = t1;
2363: else NEXT(p) = t1;
2364: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2365: break;
2366: } else if ( s == 0 ) {
2367: /* add coefs */
1.40 noro 2368: addp(CO,tm->c,c,&c1);
1.38 noro 2369: if ( c1 ) tm->c = c1;
2370: else NEXT(p) = NEXT(t);
2371: break;
2372: }
2373: }
2374: if ( !t ) {
2375: /* append */
2376: MKNODE(t1,wm,t);
2377: if ( !p ) r = t1;
2378: else NEXT(p) = t1;
2379: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2380: }
2381: } while ( ni_next(w,d) );
2382: MKNBP(u,r);
2383: return u;
2384: }
2385:
2386: int nbmtoxky(NBM a,int *b)
2387: {
2388: int d,i,j,k;
2389: int *p;
2390:
2391: d = a->d; p = a->b;
2392: for ( i = j = 0, k = 1; i < d; i++ ) {
2393: if ( !NBM_GET(p,i) ) {
2394: b[j++] = k;
2395: k = 1;
2396: } else k++;
2397: }
2398: return j;
2399: }
2400:
2401: NBP harmonic_mul_nbm(NBM a,NBM b)
2402: {
2403: int da,db,d,la,lb,lmax,lmin,l,lab,la1,lb1,lab1;
2404: int i,j,k,ia,ib,s;
2405: int *wa,*wb,*w,*wab,*wa1,*wmb;
1.40 noro 2406: P c,c1;
1.38 noro 2407: NBM wm,tm;
2408: NODE r,t1,t,p;
2409: NBP u;
2410:
2411: da = a->d; db = b->d; d = da+db;
2412: wa = (int *)ALLOCA(da*sizeof(int));
2413: wb = (int *)ALLOCA(db*sizeof(int));
2414: la = nbmtoxky(a,wa);
2415: lb = nbmtoxky(b,wb);
1.40 noro 2416: mulp(CO,a->c,b->c,&c);
1.38 noro 2417: /* wa[0],..,wa[la-1] <-> x^wa[0]y x^wa[1]y .. */
2418: /* lmax : total length */
2419: lmax = la+lb;
2420: lmin = la>lb?la:lb;
2421: w = (int *)ALLOCA(lmax*sizeof(int));
2422: /* position of a+b */
2423: wab = (int *)ALLOCA(lmax*sizeof(int));
2424: /* position of a */
2425: wa1 = (int *)ALLOCA(lmax*sizeof(int));
2426: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2427: for ( l = lmin, r = 0; l <= lmax; l++ ) {
2428: lab = lmax - l;
2429: la1 = la - lab;
2430: lb1 = lb - lab;
2431: lab1 = l-lab;
2432: /* partion l into three parts: a, b, a+b */
2433: /* initialize wab */
2434: for ( i = 0; i < lab; i++ ) wab[i] = 1;
2435: for ( ; i < l; i++ ) wab[i] = 0;
2436: do {
2437: /* initialize wa1 */
2438: for ( i = 0; i < la1; i++ ) wa1[i] = 1;
2439: for ( ; i < lab1; i++ ) wa1[i] = 0;
2440: do {
2441: ia = 0; ib = 0;
2442: for ( i = j = 0; i < l; i++ )
2443: if ( wab[i] ) w[i] = wa[ia++]+wb[ib++];
2444: else if ( wa1[j++] ) w[i] = wa[ia++];
2445: else w[i] = wb[ib++];
2446: for ( i = j = 0; i < l; i++ ) {
2447: for ( k = w[i]-1; k > 0; k--, j++ ) NBM_SET(wmb,j);
2448: NBM_CLR(wmb,j); j++;
2449: }
2450: wm->d = j; wm->c = c;
2451: for ( p = 0, t = r; t; p = t, t = NEXT(t) ) {
2452: tm = (NBM)BDY(t);
2453: s = comp_nbm(tm,wm);
2454: if ( s < 0 ) {
2455: /* insert */
2456: MKNODE(t1,wm,t);
2457: if ( !p ) r = t1;
2458: else NEXT(p) = t1;
2459: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2460: break;
2461: } else if ( s == 0 ) {
2462: /* add coefs */
1.40 noro 2463: addp(CO,tm->c,c,&c1);
1.38 noro 2464: if ( c1 ) tm->c = c1;
2465: else NEXT(p) = NEXT(t);
2466: break;
2467: }
2468: }
2469: if ( !t ) {
2470: /* append */
2471: MKNODE(t1,wm,t);
2472: if ( !p ) r = t1;
2473: else NEXT(p) = t1;
2474: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2475: }
2476: } while ( ni_next(wa1,lab1) );
2477: } while ( ni_next(wab,l) );
2478: }
2479: MKNBP(u,r);
2480: return u;
2481: }
2482: #endif