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