![[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.51 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/dist.c,v 1.50 2014/10/10 09:02:25 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.51 ! noro 868: void actm(VL vl,int nv,MP m1,MP m2,DP *pr)
! 869: {
! 870: DL d1,d2,d;
! 871: int n2,i,j,k;
! 872: Q jq,c,c1;
! 873: MP m;
! 874: P t;
! 875:
! 876: d1 = m1->dl;
! 877: d2 = m2->dl;
! 878: for ( i = 0; i < nv; i++ )
! 879: if ( d1->d[i] > d2->d[i] ) {
! 880: *pr = 0; return;
! 881: }
! 882: NEWDL(d,nv);
! 883: c = ONE;
! 884: for ( i = 0; i < nv; i++ ) {
! 885: for ( j = d2->d[i], k = d1->d[i]; k > 0; k--, j-- ) {
! 886: STOQ(j,jq); mulq(c,jq,&c1); c = c1;
! 887: }
! 888: d->d[i] = d2->d[i]-d1->d[i];
! 889: }
! 890: mulp(vl,C(m1),C(m2),&t);
! 891: NEWMP(m);
! 892: mulp(vl,(P)c,t,&C(m));
! 893: m->dl = d;
! 894: MKDP(nv,m,*pr);
! 895: }
! 896:
! 897: void weyl_actd(VL vl,DP p1,DP p2,DP *pr)
! 898: {
! 899: int n;
! 900: MP m1,m2;
! 901: DP d,r,s;
! 902:
! 903: if ( !p1 || !p2 ) *pr = 0;
! 904: else {
! 905: n = NV(p1);
! 906: r = 0;
! 907: for ( m1 = BDY(p1); m1; m1 = NEXT(m1) )
! 908: for ( m2 = BDY(p2); m2; m2 = NEXT(m2) ) {
! 909: actm(vl,n,m1,m2,&d);
! 910: addd(vl,r,d,&s); r = s;
! 911: }
! 912: *pr = r;
! 913: }
! 914: }
! 915:
1.10 noro 916: /* monomial * polynomial */
917:
1.19 noro 918: void weyl_muldm(VL vl,MP m0,DP p,DP *pr)
1.2 noro 919: {
920: DP r,t,t1;
921: MP m;
1.10 noro 922: DL d0;
923: int n,n2,l,i,j,tlen;
924: static MP *w,*psum;
925: static struct cdl *tab;
1.5 noro 926: static int wlen;
1.10 noro 927: static int rtlen;
1.2 noro 928:
929: if ( !p )
930: *pr = 0;
931: else {
1.4 noro 932: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
1.5 noro 933: if ( l > wlen ) {
1.45 noro 934: if ( w ) GCFREE(w);
1.5 noro 935: w = (MP *)MALLOC(l*sizeof(MP));
936: wlen = l;
937: }
1.4 noro 938: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
939: w[i] = m;
1.10 noro 940:
941: n = NV(p); n2 = n>>1;
942: d0 = m0->dl;
943: for ( i = 0, tlen = 1; i < n2; i++ )
944: tlen *= d0->d[n2+i]+1;
945: if ( tlen > rtlen ) {
1.45 noro 946: if ( tab ) GCFREE(tab);
947: if ( psum ) GCFREE(psum);
1.10 noro 948: rtlen = tlen;
949: tab = (struct cdl *)MALLOC(rtlen*sizeof(struct cdl));
950: psum = (MP *)MALLOC(rtlen*sizeof(MP));
951: }
952: bzero(psum,tlen*sizeof(MP));
953: for ( i = l-1; i >= 0; i-- ) {
954: bzero(tab,tlen*sizeof(struct cdl));
955: weyl_mulmm(vl,m0,w[i],n,tab,tlen);
956: for ( j = 0; j < tlen; j++ ) {
957: if ( tab[j].c ) {
958: NEWMP(m); m->dl = tab[j].d; C(m) = tab[j].c; NEXT(m) = psum[j];
959: psum[j] = m;
960: }
961: }
1.2 noro 962: }
1.10 noro 963: for ( j = tlen-1, r = 0; j >= 0; j-- )
964: if ( psum[j] ) {
965: MKDP(n,psum[j],t); addd(vl,r,t,&t1); r = t1;
966: }
1.2 noro 967: if ( r )
968: r->sugar = p->sugar + m0->dl->td;
969: *pr = r;
970: }
971: }
972:
1.10 noro 973: /* m0 = x0^d0*x1^d1*... * dx0^e0*dx1^e1*... */
974: /* rtab : array of length (e0+1)*(e1+1)*... */
1.2 noro 975:
1.19 noro 976: void weyl_mulmm(VL vl,MP m0,MP m1,int n,struct cdl *rtab,int rtablen)
1.2 noro 977: {
1.19 noro 978: P c,c0,c1;
1.10 noro 979: DL d,d0,d1,dt;
980: int i,j,a,b,k,l,n2,s,min,curlen;
981: struct cdl *p;
982: static Q *ctab;
983: static struct cdl *tab;
1.5 noro 984: static int tablen;
1.10 noro 985: static struct cdl *tmptab;
986: static int tmptablen;
1.2 noro 987:
1.10 noro 988:
989: if ( !m0 || !m1 ) {
990: rtab[0].c = 0;
991: rtab[0].d = 0;
992: return;
993: }
994: c0 = C(m0); c1 = C(m1);
995: mulp(vl,c0,c1,&c);
996: d0 = m0->dl; d1 = m1->dl;
997: n2 = n>>1;
998: curlen = 1;
999: NEWDL(d,n);
1000: if ( n & 1 )
1001: /* offset of h-degree */
1002: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
1003: else
1004: d->td = 0;
1005: rtab[0].c = c;
1006: rtab[0].d = d;
1007:
1008: if ( rtablen > tmptablen ) {
1.45 noro 1009: if ( tmptab ) GCFREE(tmptab);
1.10 noro 1010: tmptab = (struct cdl *)MALLOC(rtablen*sizeof(struct cdl));
1011: tmptablen = rtablen;
1012: }
1013: for ( i = 0; i < n2; i++ ) {
1014: a = d0->d[i]; b = d1->d[n2+i];
1015: k = d0->d[n2+i]; l = d1->d[i];
1.20 noro 1016:
1017: /* degree of xi^a*(Di^k*xi^l)*Di^b */
1018: a += l;
1019: b += k;
1020: s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);
1021:
1.10 noro 1022: if ( !k || !l ) {
1023: for ( j = 0, p = rtab; j < curlen; j++, p++ ) {
1024: if ( p->c ) {
1025: dt = p->d;
1026: dt->d[i] = a;
1027: dt->d[n2+i] = b;
1028: dt->td += s;
1.5 noro 1029: }
1.10 noro 1030: }
1031: curlen *= k+1;
1032: continue;
1033: }
1034: if ( k+1 > tablen ) {
1.45 noro 1035: if ( tab ) GCFREE(tab);
1036: if ( ctab ) GCFREE(ctab);
1.10 noro 1037: tablen = k+1;
1038: tab = (struct cdl *)MALLOC(tablen*sizeof(struct cdl));
1039: ctab = (Q *)MALLOC(tablen*sizeof(Q));
1040: }
1041: /* compute xi^a*(Di^k*xi^l)*Di^b */
1042: min = MIN(k,l);
1043: mkwc(k,l,ctab);
1044: bzero(tab,(k+1)*sizeof(struct cdl));
1045: if ( n & 1 )
1046: for ( j = 0; j <= min; j++ ) {
1047: NEWDL(d,n);
1.20 noro 1048: d->d[i] = a-j; d->d[n2+i] = b-j;
1.10 noro 1049: d->td = s;
1.20 noro 1050: d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
1.10 noro 1051: tab[j].d = d;
1052: tab[j].c = (P)ctab[j];
1053: }
1054: else
1055: for ( j = 0; j <= min; j++ ) {
1056: NEWDL(d,n);
1.20 noro 1057: d->d[i] = a-j; d->d[n2+i] = b-j;
1058: d->td = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i); /* XXX */
1.10 noro 1059: tab[j].d = d;
1060: tab[j].c = (P)ctab[j];
1061: }
1062: bzero(ctab,(min+1)*sizeof(Q));
1063: comm_muld_tab(vl,n,rtab,curlen,tab,k+1,tmptab);
1064: curlen *= k+1;
1065: bcopy(tmptab,rtab,curlen*sizeof(struct cdl));
1066: }
1067: }
1068:
1069: /* direct product of two cdl tables
1070: rt[] = [
1071: t[0]*t1[0],...,t[n-1]*t1[0],
1072: t[0]*t1[1],...,t[n-1]*t1[1],
1073: ...
1074: t[0]*t1[n1-1],...,t[n-1]*t1[n1-1]
1075: ]
1076: */
1077:
1.19 noro 1078: void comm_muld_tab(VL vl,int nv,struct cdl *t,int n,struct cdl *t1,int n1,struct cdl *rt)
1.10 noro 1079: {
1080: int i,j;
1081: struct cdl *p;
1082: P c;
1083: DL d;
1084:
1085: bzero(rt,n*n1*sizeof(struct cdl));
1086: for ( j = 0, p = rt; j < n1; j++ ) {
1087: c = t1[j].c;
1088: d = t1[j].d;
1089: if ( !c )
1090: break;
1091: for ( i = 0; i < n; i++, p++ ) {
1092: if ( t[i].c ) {
1093: mulp(vl,t[i].c,c,&p->c);
1094: adddl(nv,t[i].d,d,&p->d);
1095: }
1.6 noro 1096: }
1.1 noro 1097: }
1098: }
1099:
1.19 noro 1100: void muldc(VL vl,DP p,P c,DP *pr)
1.1 noro 1101: {
1.47 noro 1102: MP m,mr=0,mr0;
1.1 noro 1103:
1104: if ( !p || !c )
1105: *pr = 0;
1106: else if ( NUM(c) && UNIQ((Q)c) )
1107: *pr = p;
1108: else if ( NUM(c) && MUNIQ((Q)c) )
1109: chsgnd(p,pr);
1110: else {
1111: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1112: NEXTMP(mr0,mr);
1113: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
1114: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
1115: else
1116: mulp(vl,C(m),c,&C(mr));
1117: mr->dl = m->dl;
1118: }
1119: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
1120: if ( *pr )
1121: (*pr)->sugar = p->sugar;
1122: }
1.24 noro 1123: }
1124:
1125: void muldc_trunc(VL vl,DP p,P c,DL dl,DP *pr)
1126: {
1.47 noro 1127: MP m,mr=0,mr0;
1.24 noro 1128: DL mdl;
1129: int i,n;
1130:
1131: if ( !p || !c ) {
1132: *pr = 0; return;
1133: }
1134: n = NV(p);
1135: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1136: mdl = m->dl;
1137: for ( i = 0; i < n; i++ )
1138: if ( mdl->d[i] < dl->d[i] )
1139: break;
1140: if ( i < n )
1141: break;
1142: NEXTMP(mr0,mr);
1143: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
1144: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
1145: else
1146: mulp(vl,C(m),c,&C(mr));
1147: mr->dl = m->dl;
1148: }
1149: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
1150: if ( *pr )
1151: (*pr)->sugar = p->sugar;
1.1 noro 1152: }
1153:
1.19 noro 1154: void divsdc(VL vl,DP p,P c,DP *pr)
1.1 noro 1155: {
1.47 noro 1156: MP m,mr=0,mr0;
1.1 noro 1157:
1158: if ( !c )
1159: error("disvsdc : division by 0");
1160: else if ( !p )
1161: *pr = 0;
1162: else {
1163: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1164: NEXTMP(mr0,mr); divsp(vl,C(m),c,&C(mr)); mr->dl = m->dl;
1165: }
1166: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
1167: if ( *pr )
1168: (*pr)->sugar = p->sugar;
1169: }
1170: }
1171:
1.19 noro 1172: void adddl(int n,DL d1,DL d2,DL *dr)
1.1 noro 1173: {
1174: DL dt;
1175: int i;
1176:
1.44 noro 1177: *dr = dt = (DL)MALLOC_ATOMIC((n+1)*sizeof(int));
1178: dt->td = d1->td + d2->td;
1179: for ( i = 0; i < n; i++ )
1180: dt->d[i] = d1->d[i]+d2->d[i];
1.11 noro 1181: }
1182:
1183: /* d1 += d2 */
1184:
1.19 noro 1185: void adddl_destructive(int n,DL d1,DL d2)
1.11 noro 1186: {
1187: int i;
1188:
1189: d1->td += d2->td;
1190: for ( i = 0; i < n; i++ )
1191: d1->d[i] += d2->d[i];
1.1 noro 1192: }
1193:
1.19 noro 1194: int compd(VL vl,DP p1,DP p2)
1.1 noro 1195: {
1196: int n,t;
1197: MP m1,m2;
1198:
1199: if ( !p1 )
1200: return p2 ? -1 : 0;
1201: else if ( !p2 )
1202: return 1;
1.47 noro 1203: else if ( NV(p1) != NV(p2) ) {
1.39 noro 1204: error("compd : size mismatch");
1.47 noro 1205: return 0; /* XXX */
1206: } else {
1.1 noro 1207: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
1208: m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
1209: if ( (t = (*cmpdl)(n,m1->dl,m2->dl)) ||
1210: (t = compp(vl,C(m1),C(m2)) ) )
1211: return t;
1212: if ( m1 )
1213: return 1;
1214: else if ( m2 )
1215: return -1;
1216: else
1217: return 0;
1218: }
1219: }
1220:
1.19 noro 1221: int cmpdl_lex(int n,DL d1,DL d2)
1.1 noro 1222: {
1223: int i;
1224:
1225: for ( i = 0; i < n && d1->d[i] == d2->d[i]; i++ );
1226: return i == n ? 0 : (d1->d[i] > d2->d[i] ? 1 : -1);
1227: }
1228:
1.19 noro 1229: int cmpdl_revlex(int n,DL d1,DL d2)
1.1 noro 1230: {
1231: int i;
1232:
1233: for ( i = n - 1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
1234: return i < 0 ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1235: }
1236:
1.19 noro 1237: int cmpdl_gradlex(int n,DL d1,DL d2)
1.1 noro 1238: {
1239: if ( d1->td > d2->td )
1240: return 1;
1241: else if ( d1->td < d2->td )
1242: return -1;
1243: else
1244: return cmpdl_lex(n,d1,d2);
1245: }
1246:
1.19 noro 1247: int cmpdl_revgradlex(int n,DL d1,DL d2)
1.1 noro 1248: {
1.25 noro 1249: register int i,c;
1.7 noro 1250: register int *p1,*p2;
1251:
1.1 noro 1252: if ( d1->td > d2->td )
1253: return 1;
1254: else if ( d1->td < d2->td )
1255: return -1;
1.7 noro 1256: else {
1.25 noro 1257: i = n-1;
1258: p1 = d1->d+n-1;
1259: p2 = d2->d+n-1;
1260: while ( i >= 7 ) {
1261: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1262: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1263: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1264: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1265: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1266: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1267: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1268: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1269: i -= 8;
1270: }
1271: switch ( i ) {
1272: case 6:
1273: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1274: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1275: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1276: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1277: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1278: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1279: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1280: return 0;
1281: case 5:
1282: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1283: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1284: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1285: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1286: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1287: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1288: return 0;
1289: case 4:
1290: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1291: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1292: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1293: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1294: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1295: return 0;
1296: case 3:
1297: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1298: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1299: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1300: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1301: return 0;
1302: case 2:
1303: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1304: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1305: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1306: return 0;
1307: case 1:
1308: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1309: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1310: return 0;
1311: case 0:
1312: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1313: return 0;
1314: default:
1315: return 0;
1316: }
1317: LAST:
1318: if ( c > 0 ) return -1;
1319: else return 1;
1.7 noro 1320: }
1.1 noro 1321: }
1322:
1.19 noro 1323: int cmpdl_blex(int n,DL d1,DL d2)
1.1 noro 1324: {
1325: int c;
1326:
1.47 noro 1327: if ( (c = cmpdl_lex(n-1,d1,d2)) )
1.1 noro 1328: return c;
1329: else {
1330: c = d1->d[n-1] - d2->d[n-1];
1331: return c > 0 ? 1 : c < 0 ? -1 : 0;
1332: }
1333: }
1334:
1.19 noro 1335: int cmpdl_bgradlex(int n,DL d1,DL d2)
1.1 noro 1336: {
1337: int e1,e2,c;
1338:
1339: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
1340: if ( e1 > e2 )
1341: return 1;
1342: else if ( e1 < e2 )
1343: return -1;
1344: else {
1345: c = cmpdl_lex(n-1,d1,d2);
1346: if ( c )
1347: return c;
1348: else
1349: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
1350: }
1351: }
1352:
1.19 noro 1353: int cmpdl_brevgradlex(int n,DL d1,DL d2)
1.1 noro 1354: {
1355: int e1,e2,c;
1356:
1357: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
1358: if ( e1 > e2 )
1359: return 1;
1360: else if ( e1 < e2 )
1361: return -1;
1362: else {
1363: c = cmpdl_revlex(n-1,d1,d2);
1364: if ( c )
1365: return c;
1366: else
1367: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
1368: }
1369: }
1370:
1.19 noro 1371: int cmpdl_brevrev(int n,DL d1,DL d2)
1.1 noro 1372: {
1373: int e1,e2,f1,f2,c,i;
1374:
1375: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1376: e1 += d1->d[i]; e2 += d2->d[i];
1377: }
1378: f1 = d1->td - e1; f2 = d2->td - e2;
1379: if ( e1 > e2 )
1380: return 1;
1381: else if ( e1 < e2 )
1382: return -1;
1383: else {
1384: c = cmpdl_revlex(dp_nelim,d1,d2);
1385: if ( c )
1386: return c;
1387: else if ( f1 > f2 )
1388: return 1;
1389: else if ( f1 < f2 )
1390: return -1;
1391: else {
1392: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1393: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1394: }
1395: }
1396: }
1397:
1.19 noro 1398: int cmpdl_bgradrev(int n,DL d1,DL d2)
1.1 noro 1399: {
1400: int e1,e2,f1,f2,c,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: f1 = d1->td - e1; f2 = d2->td - e2;
1406: if ( e1 > e2 )
1407: return 1;
1408: else if ( e1 < e2 )
1409: return -1;
1410: else {
1411: c = cmpdl_lex(dp_nelim,d1,d2);
1412: if ( c )
1413: return c;
1414: else if ( f1 > f2 )
1415: return 1;
1416: else if ( f1 < f2 )
1417: return -1;
1418: else {
1419: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1420: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1421: }
1422: }
1423: }
1424:
1.19 noro 1425: int cmpdl_blexrev(int n,DL d1,DL d2)
1.1 noro 1426: {
1427: int e1,e2,f1,f2,c,i;
1428:
1429: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1430: e1 += d1->d[i]; e2 += d2->d[i];
1431: }
1432: f1 = d1->td - e1; f2 = d2->td - e2;
1433: c = cmpdl_lex(dp_nelim,d1,d2);
1434: if ( c )
1435: return c;
1436: else if ( f1 > f2 )
1437: return 1;
1438: else if ( f1 < f2 )
1439: return -1;
1440: else {
1441: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1442: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1443: }
1444: }
1445:
1.19 noro 1446: int cmpdl_elim(int n,DL d1,DL d2)
1.1 noro 1447: {
1448: int e1,e2,i;
1449:
1450: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1451: e1 += d1->d[i]; e2 += d2->d[i];
1452: }
1453: if ( e1 > e2 )
1454: return 1;
1455: else if ( e1 < e2 )
1456: return -1;
1457: else
1458: return cmpdl_revgradlex(n,d1,d2);
1.12 noro 1459: }
1460:
1.19 noro 1461: int cmpdl_weyl_elim(int n,DL d1,DL d2)
1.12 noro 1462: {
1463: int e1,e2,i;
1464:
1465: for ( i = 1, e1 = 0, e2 = 0; i <= dp_nelim; i++ ) {
1466: e1 += d1->d[n-i]; e2 += d2->d[n-i];
1467: }
1468: if ( e1 > e2 )
1469: return 1;
1470: else if ( e1 < e2 )
1471: return -1;
1472: else if ( d1->td > d2->td )
1473: return 1;
1474: else if ( d1->td < d2->td )
1475: return -1;
1476: else return -cmpdl_revlex(n,d1,d2);
1.13 noro 1477: }
1478:
1479: /*
1480: a special ordering
1481: 1. total order
1482: 2. (-w,w) for the first 2*m variables
1483: 3. DRL for the first 2*m variables
1484: */
1485:
1.20 noro 1486: extern int *current_weyl_weight_vector;
1.13 noro 1487:
1.19 noro 1488: int cmpdl_homo_ww_drl(int n,DL d1,DL d2)
1.13 noro 1489: {
1490: int e1,e2,m,i;
1491: int *p1,*p2;
1492:
1493: if ( d1->td > d2->td )
1494: return 1;
1495: else if ( d1->td < d2->td )
1496: return -1;
1497:
1498: m = n>>1;
1.21 noro 1499: for ( i = 0, e1 = e2 = 0, p1 = d1->d, p2 = d2->d; i < m; i++ ) {
1500: e1 += current_weyl_weight_vector[i]*(p1[m+i] - p1[i]);
1501: e2 += current_weyl_weight_vector[i]*(p2[m+i] - p2[i]);
1.13 noro 1502: }
1503: if ( e1 > e2 )
1504: return 1;
1505: else if ( e1 < e2 )
1506: return -1;
1507:
1508: e1 = d1->td - d1->d[n-1];
1509: e2 = d2->td - d2->d[n-1];
1510: if ( e1 > e2 )
1511: return 1;
1512: else if ( e1 < e2 )
1513: return -1;
1514:
1515: for ( i= n - 1, p1 = d1->d+n-1, p2 = d2->d+n-1;
1516: i >= 0 && *p1 == *p2; i--, p1--, p2-- );
1517: return i < 0 ? 0 : (*p1 < *p2 ? 1 : -1);
1.21 noro 1518: }
1519:
1520: int cmpdl_drl_zigzag(int n,DL d1,DL d2)
1521: {
1522: int i,t,m;
1523: int *p1,*p2;
1524:
1525: if ( d1->td > d2->td )
1526: return 1;
1527: else if ( d1->td < d2->td )
1528: return -1;
1529: else {
1530: m = n>>1;
1531: for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) {
1.47 noro 1532: if ( (t = p1[m+i] - p2[m+i]) ) return t > 0 ? -1 : 1;
1533: if ( (t = p1[i] - p2[i]) ) return t > 0 ? -1 : 1;
1.21 noro 1534: }
1535: return 0;
1536: }
1537: }
1538:
1539: int cmpdl_homo_ww_drl_zigzag(int n,DL d1,DL d2)
1540: {
1541: int e1,e2,m,i,t;
1542: int *p1,*p2;
1543:
1544: if ( d1->td > d2->td )
1545: return 1;
1546: else if ( d1->td < d2->td )
1547: return -1;
1548:
1549: m = n>>1;
1550: for ( i = 0, e1 = e2 = 0, p1 = d1->d, p2 = d2->d; i < m; i++ ) {
1551: e1 += current_weyl_weight_vector[i]*(p1[m+i] - p1[i]);
1552: e2 += current_weyl_weight_vector[i]*(p2[m+i] - p2[i]);
1553: }
1554: if ( e1 > e2 )
1555: return 1;
1556: else if ( e1 < e2 )
1557: return -1;
1558:
1559: e1 = d1->td - d1->d[n-1];
1560: e2 = d2->td - d2->d[n-1];
1561: if ( e1 > e2 )
1562: return 1;
1563: else if ( e1 < e2 )
1564: return -1;
1565:
1566: for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) {
1.47 noro 1567: if ( (t = p1[m+i] - p2[m+i]) ) return t > 0 ? -1 : 1;
1568: if ( (t = p1[i] - p2[i]) ) return t > 0 ? -1 : 1;
1.21 noro 1569: }
1570: return 0;
1.1 noro 1571: }
1572:
1.19 noro 1573: int cmpdl_order_pair(int n,DL d1,DL d2)
1.1 noro 1574: {
1575: int e1,e2,i,j,l;
1576: int *t1,*t2;
1.20 noro 1577: int len,head;
1.1 noro 1578: struct order_pair *pair;
1579:
1.27 noro 1580: len = dp_current_spec->ord.block.length;
1.39 noro 1581: if ( n != dp_current_spec->nv )
1582: error("cmpdl_order_pair : incompatible order specification");
1.27 noro 1583: pair = dp_current_spec->ord.block.order_pair;
1.1 noro 1584:
1.20 noro 1585: head = 0;
1.1 noro 1586: for ( i = 0, t1 = d1->d, t2 = d2->d; i < len; i++ ) {
1587: l = pair[i].length;
1588: switch ( pair[i].order ) {
1589: case 0:
1590: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
1.20 noro 1591: e1 += MUL_WEIGHT(t1[j],head+j);
1592: e2 += MUL_WEIGHT(t2[j],head+j);
1.1 noro 1593: }
1594: if ( e1 > e2 )
1595: return 1;
1596: else if ( e1 < e2 )
1597: return -1;
1598: else {
1599: for ( j = l - 1; j >= 0 && t1[j] == t2[j]; j-- );
1600: if ( j >= 0 )
1601: return t1[j] < t2[j] ? 1 : -1;
1602: }
1603: break;
1604: case 1:
1605: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
1.20 noro 1606: e1 += MUL_WEIGHT(t1[j],head+j);
1607: e2 += MUL_WEIGHT(t2[j],head+j);
1.1 noro 1608: }
1609: if ( e1 > e2 )
1610: return 1;
1611: else if ( e1 < e2 )
1612: return -1;
1613: else {
1614: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
1615: if ( j < l )
1616: return t1[j] > t2[j] ? 1 : -1;
1617: }
1618: break;
1619: case 2:
1620: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
1621: if ( j < l )
1622: return t1[j] > t2[j] ? 1 : -1;
1623: break;
1624: default:
1625: error("cmpdl_order_pair : invalid order"); break;
1626: }
1.20 noro 1627: t1 += l; t2 += l; head += l;
1.28 noro 1628: }
1629: return 0;
1630: }
1631:
1632: int cmpdl_composite(int nv,DL d1,DL d2)
1633: {
1634: int n,i,j,k,start,s,len;
1635: int *dw;
1636: struct sparse_weight *sw;
1637: struct weight_or_block *worb;
1638: int *w,*t1,*t2;
1639:
1640: n = dp_current_spec->ord.composite.length;
1641: worb = dp_current_spec->ord.composite.w_or_b;
1642: w = dp_dl_work;
1643: for ( i = 0, t1 = d1->d, t2 = d2->d; i < nv; i++ )
1644: w[i] = t1[i]-t2[i];
1645: for ( i = 0; i < n; i++, worb++ ) {
1646: len = worb->length;
1647: switch ( worb->type ) {
1648: case IS_DENSE_WEIGHT:
1649: dw = worb->body.dense_weight;
1650: for ( j = 0, s = 0; j < len; j++ )
1651: s += dw[j]*w[j];
1652: if ( s > 0 ) return 1;
1653: else if ( s < 0 ) return -1;
1654: break;
1655: case IS_SPARSE_WEIGHT:
1656: sw = worb->body.sparse_weight;
1657: for ( j = 0, s = 0; j < len; j++ )
1658: s += sw[j].value*w[sw[j].pos];
1659: if ( s > 0 ) return 1;
1660: else if ( s < 0 ) return -1;
1661: break;
1662: case IS_BLOCK:
1663: start = worb->body.block.start;
1664: switch ( worb->body.block.order ) {
1665: case 0:
1666: for ( j = 0, k = start, s = 0; j < len; j++, k++ ) {
1667: s += MUL_WEIGHT(w[k],k);
1668: }
1669: if ( s > 0 ) return 1;
1670: else if ( s < 0 ) return -1;
1671: else {
1672: for ( j = k-1; j >= start && w[j] == 0; j-- );
1673: if ( j >= start )
1674: return w[j] < 0 ? 1 : -1;
1675: }
1676: break;
1677: case 1:
1678: for ( j = 0, k = start, s = 0; j < len; j++, k++ ) {
1679: s += MUL_WEIGHT(w[k],k);
1680: }
1681: if ( s > 0 ) return 1;
1682: else if ( s < 0 ) return -1;
1683: else {
1684: for ( j = 0, k = start; j < len && w[j] == 0; j++, k++ );
1685: if ( j < len )
1686: return w[j] > 0 ? 1 : -1;
1687: }
1688: break;
1689: case 2:
1690: for ( j = 0, k = start; j < len && w[j] == 0; j++, k++ );
1691: if ( j < len )
1692: return w[j] > 0 ? 1 : -1;
1693: break;
1694: }
1695: break;
1696: }
1.1 noro 1697: }
1698: return 0;
1699: }
1700:
1.19 noro 1701: int cmpdl_matrix(int n,DL d1,DL d2)
1.1 noro 1702: {
1703: int *v,*w,*t1,*t2;
1704: int s,i,j,len;
1705: int **matrix;
1706:
1707: for ( i = 0, t1 = d1->d, t2 = d2->d, w = dp_dl_work; i < n; i++ )
1708: w[i] = t1[i]-t2[i];
1.27 noro 1709: len = dp_current_spec->ord.matrix.row;
1710: matrix = dp_current_spec->ord.matrix.matrix;
1.1 noro 1711: for ( j = 0; j < len; j++ ) {
1712: v = matrix[j];
1713: for ( i = 0, s = 0; i < n; i++ )
1714: s += v[i]*w[i];
1715: if ( s > 0 )
1716: return 1;
1717: else if ( s < 0 )
1718: return -1;
1719: }
1720: return 0;
1.25 noro 1721: }
1722:
1.43 noro 1723: int cmpdl_top_weight(int n,DL d1,DL d2)
1724: {
1725: int *w;
1726: N sum,wm,wma,t;
1.48 noro 1727: Q **mat;
1728: Q *a;
1.43 noro 1729: struct oN tn;
1.48 noro 1730: int len,i,sgn,tsgn,row,k;
1.43 noro 1731: int *t1,*t2;
1732:
1733: w = (int *)ALLOCA(n*sizeof(int));
1734: len = current_top_weight_len+3;
1735: t1 = d1->d; t2 = d2->d;
1736: for ( i = 0; i < n; i++ ) w[i] = t1[i]-t2[i];
1737: sum = (N)W_ALLOC(len); sgn = 0;
1738: wm = (N)W_ALLOC(len);
1739: wma = (N)W_ALLOC(len);
1.48 noro 1740: if ( OID(current_top_weight) == O_VECT ) {
1741: mat = (Q **)&BDY((VECT)current_top_weight);
1742: row = 1;
1743: } else {
1744: mat = (Q **)BDY((MAT)current_top_weight);
1745: row = ((MAT)current_top_weight)->row;
1746: }
1747: for ( k = 0; k < row; k++ ) {
1748: a = mat[k];
1749: for ( i = 0; i < n; i++ ) {
1750: if ( !a[i] || !w[i] ) continue;
1751: tn.p = 1;
1752: if ( w[i] > 0 ) {
1753: tn.b[0] = w[i]; tsgn = 1;
1754: } else {
1755: tn.b[0] = -w[i]; tsgn = -1;
1756: }
1757: _muln(NM(a[i]),&tn,wm);
1758: if ( !sgn ) {
1759: sgn = tsgn;
1760: t = wm; wm = sum; sum = t;
1761: } else if ( sgn == tsgn ) {
1762: _addn(sum,wm,wma);
1763: if ( !PL(wma) )
1764: sgn = 0;
1765: t = wma; wma = sum; sum = t;
1766: } else {
1767: sgn *= _subn(sum,wm,wma);
1768: t = wma; wma = sum; sum = t;
1769: }
1770: }
1771: if ( sgn > 0 ) return 1;
1772: else if ( sgn < 0 ) return -1;
1.43 noro 1773: }
1.48 noro 1774: return (*cmpdl_tie_breaker)(n,d1,d2);
1.43 noro 1775: }
1776:
1.25 noro 1777: GeoBucket create_bucket()
1778: {
1779: GeoBucket g;
1780:
1781: g = CALLOC(1,sizeof(struct oGeoBucket));
1782: g->m = 32;
1783: return g;
1784: }
1785:
1.47 noro 1786: int length(NODE d);
1787:
1.25 noro 1788: void add_bucket(GeoBucket g,NODE d,int nv)
1789: {
1790: int l,k,m;
1791:
1792: l = length(d);
1793: for ( k = 0, m = 1; l > m; k++, m <<= 1 );
1794: /* 2^(k-1) < l <= 2^k */
1795: d = symb_merge(g->body[k],d,nv);
1796: for ( ; length(d) > (1<<(k)); k++ ) {
1797: g->body[k] = 0;
1798: d = symb_merge(g->body[k+1],d,nv);
1799: }
1800: g->body[k] = d;
1801: g->m = MAX(g->m,k);
1802: }
1803:
1804: DL remove_head_bucket(GeoBucket g,int nv)
1805: {
1806: int j,i,c,m;
1807: DL d;
1808:
1809: j = -1;
1810: m = g->m;
1811: for ( i = 0; i <= m; i++ ) {
1812: if ( !g->body[i] )
1813: continue;
1814: if ( j < 0 ) j = i;
1815: else {
1816: c = (*cmpdl)(nv,g->body[i]->body,g->body[j]->body);
1817: if ( c > 0 )
1818: j = i;
1819: else if ( c == 0 )
1820: g->body[i] = NEXT(g->body[i]);
1821: }
1822: }
1823: if ( j < 0 )
1824: return 0;
1825: else {
1826: d = g->body[j]->body;
1827: g->body[j] = NEXT(g->body[j]);
1828: return d;
1.31 noro 1829: }
1830: }
1831:
1832: /* DPV functions */
1833:
1834: void adddv(VL vl,DPV p1,DPV p2,DPV *pr)
1835: {
1836: int i,len;
1837: DP *e;
1838:
1839: if ( !p1 || !p2 )
1840: error("adddv : invalid argument");
1841: else if ( p1->len != p2->len )
1842: error("adddv : size mismatch");
1843: else {
1844: len = p1->len;
1845: e = (DP *)MALLOC(p1->len*sizeof(DP));
1846: for ( i = 0; i < len; i++ )
1847: addd(vl,p1->body[i],p2->body[i],&e[i]);
1848: MKDPV(len,e,*pr);
1849: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
1850: }
1851: }
1852:
1853: void subdv(VL vl,DPV p1,DPV p2,DPV *pr)
1854: {
1855: int i,len;
1856: DP *e;
1857:
1858: if ( !p1 || !p2 )
1859: error("subdv : invalid argument");
1860: else if ( p1->len != p2->len )
1861: error("subdv : size mismatch");
1862: else {
1863: len = p1->len;
1864: e = (DP *)MALLOC(p1->len*sizeof(DP));
1865: for ( i = 0; i < len; i++ )
1866: subd(vl,p1->body[i],p2->body[i],&e[i]);
1867: MKDPV(len,e,*pr);
1868: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
1869: }
1870: }
1871:
1872: void chsgndv(DPV p1,DPV *pr)
1873: {
1874: int i,len;
1875: DP *e;
1876:
1877: if ( !p1 )
1878: error("subdv : invalid argument");
1879: else {
1880: len = p1->len;
1881: e = (DP *)MALLOC(p1->len*sizeof(DP));
1882: for ( i = 0; i < len; i++ )
1883: chsgnd(p1->body[i],&e[i]);
1884: MKDPV(len,e,*pr);
1885: (*pr)->sugar = p1->sugar;
1886: }
1887: }
1888:
1889: void muldv(VL vl,DP p1,DPV p2,DPV *pr)
1890: {
1891: int i,len;
1892: DP *e;
1893:
1894: len = p2->len;
1895: e = (DP *)MALLOC(p2->len*sizeof(DP));
1896: if ( !p1 ) {
1897: MKDPV(len,e,*pr);
1898: (*pr)->sugar = 0;
1899: } else {
1900: for ( i = 0; i < len; i++ )
1901: muld(vl,p1,p2->body[i],&e[i]);
1902: MKDPV(len,e,*pr);
1903: (*pr)->sugar = p1->sugar + p2->sugar;
1904: }
1905: }
1906:
1907: int compdv(VL vl,DPV p1,DPV p2)
1908: {
1909: int i,t,len;
1910:
1.47 noro 1911: if ( p1->len != p2->len ) {
1.31 noro 1912: error("compdv : size mismatch");
1.47 noro 1913: return 0; /* XXX */
1914: } else {
1.31 noro 1915: len = p1->len;
1916: for ( i = 0; i < len; i++ )
1.47 noro 1917: if ( (t = compd(vl,p1->body[i],p2->body[i])) )
1.31 noro 1918: return t;
1919: return 0;
1.33 noro 1920: }
1921: }
1922:
1923: int ni_next(int *a,int n)
1924: {
1925: int i,j,k,kj;
1926:
1927: /* find the first nonzero a[j] */
1.35 noro 1928: for ( j = 0; j < n && a[j] == 0; j++ );
1.33 noro 1929: /* find the first zero a[k] after a[j] */
1930: for ( k = j; k < n && a[k] == 1; k++ );
1931: if ( k == n ) return 0;
1932: /* a[0] = 0, ... , a[j-1] = 0, a[j] = 1, ..., a[k-1] = 1, a[k] = 0 */
1933: /* a[0] = 1,..., a[k-j-2] = 1, a[k-j-1] = 0, ..., a[k-1] = 0, a[k] = 1 */
1934: kj = k-j-1;
1935: for ( i = 0; i < kj; i++ ) a[i] = 1;
1936: for ( ; i < k; i++ ) a[i] = 0;
1937: a[k] = 1;
1938: return 1;
1939: }
1940:
1941: int comp_nbm(NBM a,NBM b)
1942: {
1.47 noro 1943: int d,i,ai,bi;
1.33 noro 1944: int *ab,*bb;
1945:
1946: if ( a->d > b->d ) return 1;
1947: else if ( a->d < b->d ) return -1;
1948: else {
1949: d = a->d; ab = a->b; bb = b->b;
1.41 noro 1950: #if 0
1.33 noro 1951: w = (d+31)/32;
1952: for ( i = 0; i < w; i++ )
1953: if ( ab[i] > bb[i] ) return 1;
1954: else if ( ab[i] < bb[i] ) return -1;
1.41 noro 1955: #else
1956: for ( i = 0; i < d; i++ ) {
1957: ai = NBM_GET(ab,i);
1958: bi = NBM_GET(bb,i);
1959: if ( ai > bi ) return 1;
1960: else if ( ai < bi ) return -1;
1961: }
1962: #endif
1.33 noro 1963: return 0;
1964: }
1965: }
1966:
1967: NBM mul_nbm(NBM a,NBM b)
1968: {
1969: int ad,bd,d,i,j;
1970: int *ab,*bb,*mb;
1971: NBM m;
1972:
1973: ad = a->d; bd = b->d; ab = a->b; bb = b->b;
1974: d = ad + bd;
1975: NEWNBM(m); NEWNBMBDY(m,d);
1.40 noro 1976: m->d = d; mulp(CO,a->c,b->c,&m->c); mb = m->b;
1.33 noro 1977: j = 0;
1978: for ( i = 0; i < ad; i++, j++ )
1979: if ( NBM_GET(ab,i) ) NBM_SET(mb,j);
1980: else NBM_CLR(mb,j);
1981: for ( i = 0; i < bd; i++, j++ )
1982: if ( NBM_GET(bb,i) ) NBM_SET(mb,j);
1983: else NBM_CLR(mb,j);
1984: return m;
1985: }
1986:
1.37 noro 1987: NBP nbmtonbp(NBM m)
1988: {
1989: NODE n;
1990: NBP u;
1991:
1992: MKNODE(n,m,0);
1993: MKNBP(u,n);
1994: return u;
1995: }
1996:
1997: /* a=c*x*rest -> a0= x*rest, ah=x, ar=rest */
1998:
1.40 noro 1999: P separate_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar)
1.37 noro 2000: {
2001: int i,d1;
2002: NBM t;
2003:
2004: if ( !a->d ) error("separate_nbm : invalid argument");
2005:
1.38 noro 2006: if ( a0 ) {
1.40 noro 2007: NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
1.38 noro 2008: *a0 = nbmtonbp(t);
2009: }
2010:
2011: if ( ah ) {
1.40 noro 2012: NEWNBM(t); NEWNBMBDY(t,1); t->d = 1; t->c = (P)ONE;
1.38 noro 2013: if ( NBM_GET(a->b,0) ) NBM_SET(t->b,0);
2014: else NBM_CLR(t->b,0);
2015: *ah = nbmtonbp(t);
2016: }
2017:
2018: if ( ar ) {
2019: d1 = a->d-1;
1.40 noro 2020: NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
1.38 noro 2021: for ( i = 0; i < d1; i++ ) {
2022: if ( NBM_GET(a->b,i+1) ) NBM_SET(t->b,i);
1.42 noro 2023: else NBM_CLR(t->b,i);
2024: }
2025: *ar = nbmtonbp(t);
2026: }
2027:
2028: return a->c;
2029: }
2030:
2031: /* a=c*rest*x -> a0= rest*x, ar=rest, at=x */
2032:
2033: P separate_tail_nbm(NBM a,NBP *a0,NBP *ar,NBP *at)
2034: {
2035: int i,d,d1;
2036: NBM t;
2037:
2038: if ( !(d=a->d) ) error("separate_tail_nbm : invalid argument");
2039:
2040: if ( a0 ) {
2041: NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
2042: *a0 = nbmtonbp(t);
2043: }
2044:
2045: d1 = a->d-1;
2046: if ( at ) {
2047: NEWNBM(t); NEWNBMBDY(t,1); t->d = 1; t->c = (P)ONE;
2048: if ( NBM_GET(a->b,d1) ) NBM_SET(t->b,0);
2049: else NBM_CLR(t->b,0);
2050: *at = nbmtonbp(t);
2051: }
2052:
2053: if ( ar ) {
2054: NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
2055: for ( i = 0; i < d1; i++ ) {
2056: if ( NBM_GET(a->b,i) ) NBM_SET(t->b,i);
1.38 noro 2057: else NBM_CLR(t->b,i);
2058: }
2059: *ar = nbmtonbp(t);
2060: }
1.37 noro 2061:
2062: return a->c;
2063: }
2064:
2065: NBP make_xky(int k)
2066: {
2067: int k1,i;
2068: NBM t;
2069:
1.40 noro 2070: NEWNBM(t); NEWNBMBDY(t,k); t->d = k; t->c = (P)ONE;
1.37 noro 2071: k1 = k-1;
2072: for ( i = 0; i < k1; i++ ) NBM_SET(t->b,i);
2073: NBM_CLR(t->b,i);
2074: return nbmtonbp(t);
2075: }
2076:
2077: /* a=c*x^(k-1)*y*rest -> a0= x^(k-1)*y*rest, ah=x^(k-1)*y, ar=rest */
2078:
1.40 noro 2079: P separate_xky_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar)
1.37 noro 2080: {
2081: int i,d1,k,k1;
2082: NBM t;
2083:
2084: if ( !a->d )
2085: error("separate_nbm : invalid argument");
2086: for ( i = 0; i < a->d && NBM_GET(a->b,i); i++ );
2087: if ( i == a->d )
2088: error("separate_nbm : invalid argument");
2089: k1 = i;
2090: k = i+1;
2091:
1.38 noro 2092: if ( a0 ) {
1.40 noro 2093: NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
1.38 noro 2094: *a0 = nbmtonbp(t);
2095: }
2096:
2097: if ( ah ) {
1.40 noro 2098: NEWNBM(t); NEWNBMBDY(t,k); t->d = k; t->c = (P)ONE;
1.38 noro 2099: for ( i = 0; i < k1; i++ ) NBM_SET(t->b,i);
2100: NBM_CLR(t->b,i);
2101: *ah = nbmtonbp(t);
2102: }
2103:
2104: if ( ar ) {
2105: d1 = a->d-k;
1.40 noro 2106: NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
1.38 noro 2107: for ( i = 0; i < d1; i++ ) {
2108: if ( NBM_GET(a->b,i+k) ) NBM_SET(t->b,i);
2109: else NBM_CLR(t->b,i);
2110: }
2111: *ar = nbmtonbp(t);
1.37 noro 2112: }
2113:
2114: return a->c;
2115: }
1.33 noro 2116:
1.37 noro 2117: void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp);
2118: void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp);
1.38 noro 2119: void mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp);
2120: void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp);
1.37 noro 2121:
2122: NBP shuffle_mul_nbm(NBM a,NBM b)
2123: {
2124: NBP u,a0,ah,ar,b0,bh,br,a1,b1,t;
1.40 noro 2125: P ac,bc,c;
1.37 noro 2126:
2127: if ( !a->d || !b->d )
2128: u = nbmtonbp(mul_nbm(a,b));
2129: else {
2130: ac = separate_nbm(a,&a0,&ah,&ar);
2131: bc = separate_nbm(b,&b0,&bh,&br);
1.40 noro 2132: mulp(CO,ac,bc,&c);
1.37 noro 2133: shuffle_mulnbp(CO,ar,b0,&t); mulnbp(CO,ah,t,&a1);
2134: shuffle_mulnbp(CO,a0,br,&t); mulnbp(CO,bh,t,&b1);
2135: addnbp(CO,a1,b1,&t); mulnbp(CO,(NBP)c,t,&u);
2136: }
2137: return u;
2138: }
1.33 noro 2139:
1.37 noro 2140: NBP harmonic_mul_nbm(NBM a,NBM b)
2141: {
2142: NBP u,a0,ah,ar,b0,bh,br,a1,b1,t,s,abk,ab1;
1.40 noro 2143: P ac,bc,c;
1.37 noro 2144:
2145: if ( !a->d || !b->d )
2146: u = nbmtonbp(mul_nbm(a,b));
2147: else {
1.40 noro 2148: mulp(CO,a->c,b->c,&c);
1.37 noro 2149: ac = separate_xky_nbm(a,&a0,&ah,&ar);
2150: bc = separate_xky_nbm(b,&b0,&bh,&br);
1.40 noro 2151: mulp(CO,ac,bc,&c);
1.37 noro 2152: harmonic_mulnbp(CO,ar,b0,&t); mulnbp(CO,ah,t,&a1);
2153: harmonic_mulnbp(CO,a0,br,&t); mulnbp(CO,bh,t,&b1);
2154: abk = make_xky(((NBM)BDY(BDY(ah)))->d+((NBM)BDY(BDY(bh)))->d);
2155: harmonic_mulnbp(CO,ar,br,&t); mulnbp(CO,abk,t,&ab1);
2156: addnbp(CO,a1,b1,&t); addnbp(CO,t,ab1,&s); mulnbp(CO,(NBP)c,s,&u);
2157: }
2158: return u;
2159:
2160: }
1.34 noro 2161:
1.33 noro 2162: void addnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2163: {
1.47 noro 2164: NODE b1,b2,br=0,br0;
1.33 noro 2165: NBM m1,m2,m;
1.40 noro 2166: P c;
1.33 noro 2167:
2168: if ( !p1 )
2169: *rp = p2;
2170: else if ( !p2 )
2171: *rp = p1;
2172: else {
2173: for ( b1 = BDY(p1), b2 = BDY(p2), br0 = 0; b1 && b2; ) {
2174: m1 = (NBM)BDY(b1); m2 = (NBM)BDY(b2);
2175: switch ( comp_nbm(m1,m2) ) {
2176: case 0:
1.40 noro 2177: addp(CO,m1->c,m2->c,&c);
1.33 noro 2178: if ( c ) {
2179: NEXTNODE(br0,br);
2180: NEWNBM(m); m->d = m1->d; m->c = c; m->b = m1->b;
2181: BDY(br) = (pointer)m;
2182: }
2183: b1 = NEXT(b1); b2 = NEXT(b2); break;
2184: case 1:
2185: NEXTNODE(br0,br); BDY(br) = BDY(b1);
2186: b1 = NEXT(b1); break;
2187: case -1:
2188: NEXTNODE(br0,br); BDY(br) = BDY(b2);
2189: b2 = NEXT(b2); break;
2190: }
1.34 noro 2191: }
2192: if ( !br0 )
2193: if ( b1 )
2194: br0 = b1;
1.33 noro 2195: else if ( b2 )
1.34 noro 2196: br0 = b2;
2197: else {
2198: *rp = 0;
2199: return;
2200: }
2201: else if ( b1 )
2202: NEXT(br) = b1;
2203: else if ( b2 )
1.33 noro 2204: NEXT(br) = b2;
1.34 noro 2205: else
2206: NEXT(br) = 0;
2207: MKNBP(*rp,br0);
1.33 noro 2208: }
2209: }
2210:
2211: void subnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2212: {
2213: NBP t;
2214:
2215: chsgnnbp(p2,&t);
2216: addnbp(vl,p1,t,rp);
2217: }
2218:
2219: void chsgnnbp(NBP p,NBP *rp)
2220: {
1.47 noro 2221: NODE r0,r=0,b;
1.33 noro 2222: NBM m,m1;
2223:
2224: for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
2225: NEXTNODE(r0,r);
2226: m = (NBM)BDY(b);
1.40 noro 2227: NEWNBM(m1); m1->d = m->d; m1->b = m->b; chsgnp(m->c,&m1->c);
1.34 noro 2228: BDY(r) = m1;
1.33 noro 2229: }
2230: if ( r0 ) NEXT(r) = 0;
2231: MKNBP(*rp,r0);
2232: }
2233:
2234: void mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2235: {
1.36 noro 2236: NODE b,n;
1.33 noro 2237: NBP r,t,s;
1.36 noro 2238: NBM m;
1.33 noro 2239:
1.36 noro 2240: if ( !p1 || !p2 ) {
2241: *rp = 0; return;
2242: }
2243: if ( OID(p1) != O_NBP ) {
1.40 noro 2244: if ( !POLY(p1) )
1.37 noro 2245: error("mulnbp : invalid argument");
1.40 noro 2246: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
1.36 noro 2247: MKNODE(n,m,0); MKNBP(p1,n);
2248: }
2249: if ( OID(p2) != O_NBP ) {
1.40 noro 2250: if ( !POLY(p2) )
1.37 noro 2251: error("mulnbp : invalid argument");
1.40 noro 2252: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
1.36 noro 2253: MKNODE(n,m,0); MKNBP(p2,n);
2254: }
2255: if ( length(BDY(p1)) < length(BDY(p2)) ) {
1.33 noro 2256: for ( r = 0, b = BDY(p1); b; b = NEXT(b) ) {
2257: mulnbmnbp(vl,(NBM)BDY(b),p2,&t);
2258: addnbp(vl,r,t,&s); r = s;
2259: }
2260: *rp = r;
2261: } else {
2262: for ( r = 0, b = BDY(p2); b; b = NEXT(b) ) {
2263: mulnbpnbm(vl,p1,(NBM)BDY(b),&t);
2264: addnbp(vl,r,t,&s); r = s;
2265: }
2266: *rp = r;
2267: }
2268: }
2269:
2270: void mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp)
2271: {
1.47 noro 2272: NODE b,r0,r=0;
1.33 noro 2273:
2274: if ( !p ) *rp = 0;
2275: else {
2276: for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
2277: NEXTNODE(r0,r);
2278: BDY(r) = mul_nbm(m,(NBM)BDY(b));
2279: }
2280: if ( r0 ) NEXT(r) = 0;
2281: MKNBP(*rp,r0);
2282: }
2283: }
2284:
2285: void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp)
2286: {
1.47 noro 2287: NODE b,r0,r=0;
1.33 noro 2288:
2289: if ( !p ) *rp = 0;
2290: else {
2291: for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
2292: NEXTNODE(r0,r);
2293: BDY(r) = mul_nbm((NBM)BDY(b),m);
2294: }
2295: if ( r0 ) NEXT(r) = 0;
2296: MKNBP(*rp,r0);
2297: }
2298: }
2299:
2300: void pwrnbp(VL vl,NBP a,Q q,NBP *c)
2301: {
2302: int t;
2303: NBP a1,a2;
2304: N n1;
2305: Q q1;
2306: NBM m;
2307: NODE r;
2308:
2309: if ( !q ) {
1.40 noro 2310: NEWNBM(m); m->d = 0; m->c = (P)ONE; m->b = 0;
1.33 noro 2311: MKNODE(r,m,0); MKNBP(*c,r);
2312: } else if ( !a )
2313: *c = 0;
2314: else if ( UNIQ(q) )
2315: *c = a;
2316: else {
2317: t = divin(NM(q),2,&n1); NTOQ(n1,1,q1);
2318: pwrnbp(vl,a,q1,&a1);
2319: mulnbp(vl,a1,a1,&a2);
2320: if ( t )
2321: mulnbp(vl,a2,a,c);
2322: else
2323: *c = a2;
2324: }
2325: }
2326:
1.38 noro 2327: int compnbp(VL vl,NBP p1,NBP p2)
2328: {
2329: NODE n1,n2;
2330: NBM m1,m2;
2331: int t;
2332:
2333: if ( !p1 )
2334: return p2 ? -1 : 0;
2335: else if ( !p2 )
2336: return 1;
2337: else {
2338: for ( n1 = BDY(p1), n2 = BDY(p2);
2339: n1 && n2; n1 = NEXT(n1), n2 = NEXT(n2) ) {
2340: m1 = (NBM)BDY(n1); m2 = (NBM)BDY(n2);
1.40 noro 2341: if ( (t = comp_nbm(m1,m2)) || (t = compp(CO,m1->c,m2->c) ) )
1.38 noro 2342: return t;
2343: }
2344: if ( n1 )
2345: return 1;
2346: else if ( n2 )
2347: return -1;
2348: else
2349: return 0;
2350: }
2351: }
2352:
1.33 noro 2353: void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2354: {
1.36 noro 2355: NODE b1,b2,n;
1.33 noro 2356: NBP r,t,s;
1.34 noro 2357: NBM m;
1.33 noro 2358:
1.36 noro 2359: if ( !p1 || !p2 ) {
2360: *rp = 0; return;
1.33 noro 2361: }
1.36 noro 2362: if ( OID(p1) != O_NBP ) {
1.40 noro 2363: if ( !POLY(p1) )
1.37 noro 2364: error("shuffle_mulnbp : invalid argument");
1.40 noro 2365: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
1.36 noro 2366: MKNODE(n,m,0); MKNBP(p1,n);
2367: }
2368: if ( OID(p2) != O_NBP ) {
1.40 noro 2369: if ( !POLY(p2) )
1.37 noro 2370: error("shuffle_mulnbp : invalid argument");
1.40 noro 2371: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
1.36 noro 2372: MKNODE(n,m,0); MKNBP(p2,n);
2373: }
2374: for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) )
2375: for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) {
2376: t = shuffle_mul_nbm(m,(NBM)BDY(b2));
2377: addnbp(vl,r,t,&s); r = s;
2378: }
2379: *rp = r;
1.33 noro 2380: }
2381:
1.34 noro 2382: void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
1.33 noro 2383: {
1.36 noro 2384: NODE b1,b2,n;
1.34 noro 2385: NBP r,t,s;
2386: NBM m;
1.33 noro 2387:
1.36 noro 2388: if ( !p1 || !p2 ) {
2389: *rp = 0; return;
1.25 noro 2390: }
1.36 noro 2391: if ( OID(p1) != O_NBP ) {
1.40 noro 2392: if ( !POLY(p1) )
1.37 noro 2393: error("harmonic_mulnbp : invalid argument");
1.40 noro 2394: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
1.36 noro 2395: MKNODE(n,m,0); MKNBP(p1,n);
2396: }
2397: if ( OID(p2) != O_NBP ) {
1.40 noro 2398: if ( !POLY(p2) )
1.37 noro 2399: error("harmonic_mulnbp : invalid argument");
1.40 noro 2400: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
1.36 noro 2401: MKNODE(n,m,0); MKNBP(p2,n);
2402: }
2403: for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) )
2404: for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) {
2405: t = harmonic_mul_nbm(m,(NBM)BDY(b2));
2406: addnbp(vl,r,t,&s); r = s;
2407: }
2408: *rp = r;
1.1 noro 2409: }
1.38 noro 2410:
2411: #if 0
2412: NBP shuffle_mul_nbm(NBM a,NBM b)
2413: {
2414: int ad,bd,d,i,ai,bi,bit,s;
2415: int *ab,*bb,*wmb,*w;
2416: NBM wm,tm;
1.40 noro 2417: P c,c1;
1.38 noro 2418: NODE r,t,t1,p;
2419: NBP u;
2420:
2421: ad = a->d; bd = b->d; ab = a->b; bb = b->b;
2422: d = ad + bd;
2423: w = (int *)ALLOCA(d*sizeof(int));
2424: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2425: for ( i = 0; i < ad; i++ ) w[i] = 1;
2426: for ( ; i < d; i++ ) w[i] = 0;
1.40 noro 2427: mulp(CO,a->c,b->c,&c);
1.38 noro 2428: r = 0;
2429: do {
2430: wm->d = d; wm->c = c;
2431: ai = 0; bi = 0;
2432: for ( i = 0; i < d; i++ ) {
2433: if ( w[i] ) { bit = NBM_GET(ab,ai); ai++; }
2434: else { bit = NBM_GET(bb,bi); bi++; }
2435: if ( bit ) NBM_SET(wmb,i);
2436: else NBM_CLR(wmb,i);
2437: }
2438: for ( p = 0, t = r; t; p = t, t = NEXT(t) ) {
2439: tm = (NBM)BDY(t);
2440: s = comp_nbm(tm,wm);
2441: if ( s < 0 ) {
2442: /* insert */
2443: MKNODE(t1,wm,t);
2444: if ( !p ) r = t1;
2445: else NEXT(p) = t1;
2446: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2447: break;
2448: } else if ( s == 0 ) {
2449: /* add coefs */
1.40 noro 2450: addp(CO,tm->c,c,&c1);
1.38 noro 2451: if ( c1 ) tm->c = c1;
2452: else NEXT(p) = NEXT(t);
2453: break;
2454: }
2455: }
2456: if ( !t ) {
2457: /* append */
2458: MKNODE(t1,wm,t);
2459: if ( !p ) r = t1;
2460: else NEXT(p) = t1;
2461: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2462: }
2463: } while ( ni_next(w,d) );
2464: MKNBP(u,r);
2465: return u;
2466: }
2467:
2468: int nbmtoxky(NBM a,int *b)
2469: {
2470: int d,i,j,k;
2471: int *p;
2472:
2473: d = a->d; p = a->b;
2474: for ( i = j = 0, k = 1; i < d; i++ ) {
2475: if ( !NBM_GET(p,i) ) {
2476: b[j++] = k;
2477: k = 1;
2478: } else k++;
2479: }
2480: return j;
2481: }
2482:
2483: NBP harmonic_mul_nbm(NBM a,NBM b)
2484: {
2485: int da,db,d,la,lb,lmax,lmin,l,lab,la1,lb1,lab1;
2486: int i,j,k,ia,ib,s;
2487: int *wa,*wb,*w,*wab,*wa1,*wmb;
1.40 noro 2488: P c,c1;
1.38 noro 2489: NBM wm,tm;
2490: NODE r,t1,t,p;
2491: NBP u;
2492:
2493: da = a->d; db = b->d; d = da+db;
2494: wa = (int *)ALLOCA(da*sizeof(int));
2495: wb = (int *)ALLOCA(db*sizeof(int));
2496: la = nbmtoxky(a,wa);
2497: lb = nbmtoxky(b,wb);
1.40 noro 2498: mulp(CO,a->c,b->c,&c);
1.38 noro 2499: /* wa[0],..,wa[la-1] <-> x^wa[0]y x^wa[1]y .. */
2500: /* lmax : total length */
2501: lmax = la+lb;
2502: lmin = la>lb?la:lb;
2503: w = (int *)ALLOCA(lmax*sizeof(int));
2504: /* position of a+b */
2505: wab = (int *)ALLOCA(lmax*sizeof(int));
2506: /* position of a */
2507: wa1 = (int *)ALLOCA(lmax*sizeof(int));
2508: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2509: for ( l = lmin, r = 0; l <= lmax; l++ ) {
2510: lab = lmax - l;
2511: la1 = la - lab;
2512: lb1 = lb - lab;
2513: lab1 = l-lab;
2514: /* partion l into three parts: a, b, a+b */
2515: /* initialize wab */
2516: for ( i = 0; i < lab; i++ ) wab[i] = 1;
2517: for ( ; i < l; i++ ) wab[i] = 0;
2518: do {
2519: /* initialize wa1 */
2520: for ( i = 0; i < la1; i++ ) wa1[i] = 1;
2521: for ( ; i < lab1; i++ ) wa1[i] = 0;
2522: do {
2523: ia = 0; ib = 0;
2524: for ( i = j = 0; i < l; i++ )
2525: if ( wab[i] ) w[i] = wa[ia++]+wb[ib++];
2526: else if ( wa1[j++] ) w[i] = wa[ia++];
2527: else w[i] = wb[ib++];
2528: for ( i = j = 0; i < l; i++ ) {
2529: for ( k = w[i]-1; k > 0; k--, j++ ) NBM_SET(wmb,j);
2530: NBM_CLR(wmb,j); j++;
2531: }
2532: wm->d = j; wm->c = c;
2533: for ( p = 0, t = r; t; p = t, t = NEXT(t) ) {
2534: tm = (NBM)BDY(t);
2535: s = comp_nbm(tm,wm);
2536: if ( s < 0 ) {
2537: /* insert */
2538: MKNODE(t1,wm,t);
2539: if ( !p ) r = t1;
2540: else NEXT(p) = t1;
2541: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2542: break;
2543: } else if ( s == 0 ) {
2544: /* add coefs */
1.40 noro 2545: addp(CO,tm->c,c,&c1);
1.38 noro 2546: if ( c1 ) tm->c = c1;
2547: else NEXT(p) = NEXT(t);
2548: break;
2549: }
2550: }
2551: if ( !t ) {
2552: /* append */
2553: MKNODE(t1,wm,t);
2554: if ( !p ) r = t1;
2555: else NEXT(p) = t1;
2556: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2557: }
2558: } while ( ni_next(wa1,lab1) );
2559: } while ( ni_next(wab,l) );
2560: }
2561: MKNBP(u,r);
2562: return u;
2563: }
2564: #endif