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