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