Annotation of OpenXM_contrib2/asir2018/engine/dist.c, Revision 1.23
1.1 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
26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
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.23 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2018/engine/dist.c,v 1.22 2020/02/11 01:43:57 noro Exp $
1.1 noro 49: */
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
62: #define ORD_WEYL_ELIM 10
63: #define ORD_HOMO_WW_DRL 11
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();
68: int cmpdl_top_weight();
69:
70: int (*cmpdl)()=cmpdl_revgradlex;
71: int (*cmpdl_tie_breaker)();
72: int (*primitive_cmpdl[3])() = {cmpdl_revgradlex,cmpdl_gradlex,cmpdl_lex};
73:
74: Obj current_top_weight;
75: int current_top_weight_len;
76:
77: int do_weyl;
78:
79: int dp_nelim,dp_fcoeffs;
80: struct order_spec *dp_current_spec;
81: struct modorder_spec *dp_current_modspec;
82: int *dp_dl_work;
83:
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,Obj c,DL dl,DP *pr);
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);
90:
91: void order_init()
92: {
93: struct order_spec *spec;
94:
95: create_order_spec(0,0,&spec);
96: initd(spec);
97: create_modorder_spec(0,0,&dp_current_modspec);
98: }
99:
100: int has_sfcoef_p(Obj f);
101:
102: int has_sfcoef(DP f)
103: {
104: MP t;
105:
106: if ( !f )
107: return 0;
108: for ( t = BDY(f); t; t = NEXT(t) )
109: if ( has_sfcoef_p(t->c) )
110: break;
111: return t ? 1 : 0;
112: }
113:
114: int has_sfcoef_p(Obj f)
115: {
116: DCP dc;
117:
118: if ( !f )
119: return 0;
120: else if ( NUM(f) )
121: return (NID((Num)f) == N_GFS) ? 1 : 0;
122: else if ( POLY(f) ) {
123: for ( dc = DC((P)f); dc; dc = NEXT(dc) )
124: if ( has_sfcoef_p((Obj)COEF(dc)) )
125: return 1;
126: return 0;
127: } else
128: return 0;
129: }
130:
131: extern Obj nd_top_weight;
132:
133: void reset_top_weight()
134: {
135: cmpdl = cmpdl_tie_breaker;
136: cmpdl_tie_breaker = 0;
137: nd_top_weight = 0;
138: current_top_weight = 0;
139: current_top_weight_len = 0;
140: }
141:
142: void initd(struct order_spec *spec)
143: {
144: int len,i,k,row;
145: Q **mat;
146:
147: switch ( spec->id ) {
148: case 3:
149: cmpdl = cmpdl_composite;
150: dp_dl_work = (int *)MALLOC_ATOMIC(spec->nv*sizeof(int));
151: break;
152: case 2:
153: cmpdl = cmpdl_matrix;
154: dp_dl_work = (int *)MALLOC_ATOMIC(spec->nv*sizeof(int));
155: break;
156: case 1:
157: cmpdl = cmpdl_order_pair;
158: break;
159: default:
160: switch ( spec->ord.simple ) {
161: case ORD_REVGRADLEX:
162: cmpdl = cmpdl_revgradlex; break;
163: case ORD_GRADLEX:
164: cmpdl = cmpdl_gradlex; break;
165: case ORD_BREVGRADLEX:
166: cmpdl = cmpdl_brevgradlex; break;
167: case ORD_BGRADLEX:
168: cmpdl = cmpdl_bgradlex; break;
169: case ORD_BLEX:
170: cmpdl = cmpdl_blex; break;
171: case ORD_BREVREV:
172: cmpdl = cmpdl_brevrev; break;
173: case ORD_BGRADREV:
174: cmpdl = cmpdl_bgradrev; break;
175: case ORD_BLEXREV:
176: cmpdl = cmpdl_blexrev; break;
177: case ORD_ELIM:
178: cmpdl = cmpdl_elim; break;
179: case ORD_WEYL_ELIM:
180: cmpdl = cmpdl_weyl_elim; break;
181: case ORD_HOMO_WW_DRL:
182: cmpdl = cmpdl_homo_ww_drl; break;
183: case ORD_DRL_ZIGZAG:
184: cmpdl = cmpdl_drl_zigzag; break;
185: case ORD_HOMO_WW_DRL_ZIGZAG:
186: cmpdl = cmpdl_homo_ww_drl_zigzag; break;
187: case ORD_LEX: default:
188: cmpdl = cmpdl_lex; break;
189: }
190: break;
191: }
192: if ( current_top_weight ) {
193: cmpdl_tie_breaker = cmpdl;
194: cmpdl = cmpdl_top_weight;
195: if ( OID(current_top_weight) == O_VECT ) {
196: mat = (Q **)&BDY((VECT)current_top_weight);
197: row = 1;
198: } else {
199: mat = (Q **)BDY((MAT)current_top_weight);
200: row = ((MAT)current_top_weight)->row;
201: }
202: for ( k = 0, len = 0; k < row; k++ )
203: for ( i = 0; i < spec->nv; i++ )
204: if ( mat[k][i] )
205: len = MAX(mpz_size(BDY((Z)mat[k][i])),len);
206: current_top_weight_len = len;
207: }
208: dp_current_spec = spec;
209: }
210:
1.3 noro 211: int dpm_ordtype;
1.1 noro 212:
213: void ptod(VL vl,VL dvl,P p,DP *pr)
214: {
215: int n,i,j,k;
216: VL tvl;
217: V v;
218: DL d;
219: MP m;
220: DCP dc;
221: DCP *w;
222: DP r,s,t,u;
223: P x,c;
224:
225: if ( !p )
226: *pr = 0;
227: else if ( OID(p) > O_P )
228: error("ptod : only polynomials can be converted.");
229: else {
230: for ( n = 0, tvl = dvl; tvl; tvl = NEXT(tvl), n++ );
231: if ( NUM(p) ) {
232: NEWDL(d,n);
233: NEWMP(m); m->dl = d; C(m) = (Obj)p; NEXT(m) = 0; MKDP(n,m,*pr); (*pr)->sugar = 0;
234: } else {
235: for ( i = 0, tvl = dvl, v = VR(p);
236: tvl && tvl->v != v; tvl = NEXT(tvl), i++ );
237: if ( !tvl ) {
238: for ( dc = DC(p), k = 0; dc; dc = NEXT(dc), k++ );
239: w = (DCP *)ALLOCA(k*sizeof(DCP));
240: for ( dc = DC(p), j = 0; j < k; dc = NEXT(dc), j++ )
241: w[j] = dc;
242:
243: for ( j = k-1, s = 0, MKV(v,x); j >= 0; j-- ) {
244: ptod(vl,dvl,COEF(w[j]),&t); pwrp(vl,x,DEG(w[j]),&c);
245: muldc(vl,t,(Obj)c,&r); addd(vl,r,s,&t); s = t;
246: }
247: *pr = s;
248: } else {
249: for ( dc = DC(p), k = 0; dc; dc = NEXT(dc), k++ );
250: w = (DCP *)ALLOCA(k*sizeof(DCP));
251: for ( dc = DC(p), j = 0; j < k; dc = NEXT(dc), j++ )
252: w[j] = dc;
253:
254: for ( j = k-1, s = 0; j >= 0; j-- ) {
255: ptod(vl,dvl,COEF(w[j]),&t);
1.2 noro 256: NEWDL(d,n); d->d[i] = ZTOS(DEG(w[j]));
1.1 noro 257: d->td = MUL_WEIGHT(d->d[i],i);
258: NEWMP(m); m->dl = d; C(m) = (Obj)ONE; NEXT(m) = 0; MKDP(n,m,u); u->sugar = d->td;
259: comm_muld(vl,t,u,&r); addd(vl,r,s,&t); s = t;
260: }
261: *pr = s;
262: }
263: }
264: }
265: #if 0
266: if ( !dp_fcoeffs && has_sfcoef(*pr) )
267: dp_fcoeffs = N_GFS;
268: #endif
269: }
270:
271: void dtop(VL vl,VL dvl,DP p,Obj *pr)
272: {
273: int n,i,j,k;
274: DL d;
275: MP m;
276: MP *a;
277: P r;
278: Obj t,w,s,u;
279: Z q;
280: VL tvl;
281:
282: if ( !p )
283: *pr = 0;
284: else {
285: for ( k = 0, m = BDY(p); m; m = NEXT(m), k++ );
286: a = (MP *)ALLOCA(k*sizeof(MP));
287: for ( j = 0, m = BDY(p); j < k; m = NEXT(m), j++ )
288: a[j] = m;
289:
290: for ( n = p->nv, j = k-1, s = 0; j >= 0; j-- ) {
291: m = a[j];
292: t = C(m);
293: if ( NUM(t) && NID((Num)t) == N_M ) {
294: mptop((P)t,(P *)&u); t = u;
295: }
296: for ( i = 0, d = m->dl, tvl = dvl;
297: i < n; tvl = NEXT(tvl), i++ ) {
1.2 noro 298: MKV(tvl->v,r); STOZ(d->d[i],q); pwrp(vl,r,q,(P *)&u);
1.1 noro 299: arf_mul(vl,t,(Obj)u,&w); t = w;
300: }
301: arf_add(vl,s,t,&u); s = u;
302: }
303: *pr = s;
304: }
305: }
306:
307: void nodetod(NODE node,DP *dp)
308: {
309: NODE t;
310: int len,i,td;
311: Q e;
312: DL d;
313: MP m;
314: DP u;
315:
316: for ( t = node, len = 0; t; t = NEXT(t), len++ );
317: NEWDL(d,len);
318: for ( t = node, i = 0, td = 0; i < len; t = NEXT(t), i++ ) {
319: e = (Q)BDY(t);
320: if ( !e )
321: d->d[i] = 0;
322: else if ( !NUM(e) || !RATN(e) || !INT(e) )
323: error("nodetod : invalid input");
324: else {
1.2 noro 325: d->d[i] = ZTOS((Q)e); td += MUL_WEIGHT(d->d[i],i);
1.1 noro 326: }
327: }
328: d->td = td;
329: NEWMP(m); m->dl = d; C(m) = (Obj)ONE; NEXT(m) = 0;
330: MKDP(len,m,u); u->sugar = td; *dp = u;
331: }
332:
333: void nodetodpm(NODE node,Obj pos,DPM *dp)
334: {
335: NODE t;
1.9 noro 336: int len,i,td,p;
1.1 noro 337: Q e;
338: DL d;
339: DMM m;
340: DPM u;
341:
342: for ( t = node, len = 0; t; t = NEXT(t), len++ );
343: NEWDL(d,len);
344: for ( t = node, i = 0, td = 0; i < len; t = NEXT(t), i++ ) {
345: e = (Q)BDY(t);
346: if ( !e )
347: d->d[i] = 0;
348: else if ( !NUM(e) || !RATN(e) || !INT(e) )
349: error("nodetodpm : invalid input");
350: else {
1.2 noro 351: d->d[i] = ZTOS((Q)e); td += MUL_WEIGHT(d->d[i],i);
1.1 noro 352: }
353: }
354: d->td = td;
1.9 noro 355: p = ZTOS((Q)pos);
1.19 noro 356: if ( dp_current_spec->module_top_weight ) {
1.9 noro 357: if ( p > dp_current_spec->module_rank )
358: error("nodetodpm : inconsistent order spec");
359: d->td += dp_current_spec->module_top_weight[p-1];
360: }
361: NEWDMM(m); m->dl = d; m->pos = p; C(m) = (Obj)ONE; NEXT(m) = 0;
1.1 noro 362: MKDPM(len,m,u); u->sugar = td; *dp = u;
363: }
364:
365: void dtodpm(DP d,int pos,DPM *dp)
366: {
367: DMM mr0,mr;
368: MP m;
1.9 noro 369: int shift;
1.1 noro 370:
371: if ( !d ) *dp = 0;
372: else {
1.9 noro 373: shift = 0;
1.19 noro 374: if ( dp_current_spec->module_top_weight ) {
1.9 noro 375: if ( pos > dp_current_spec->module_rank )
376: error("nodetodpm : inconsistent order spec");
377: shift = dp_current_spec->module_top_weight[pos-1];
378: }
1.1 noro 379: for ( m = BDY(d), mr0 = 0; m; m = NEXT(m) ) {
380: NEXTDMM(mr0,mr);
381: mr->dl = m->dl;
1.9 noro 382: mr->dl->td += shift;
1.1 noro 383: mr->pos = pos;
384: C(mr) = C(m);
385: }
386: MKDPM(d->nv,mr0,*dp); (*dp)->sugar = d->sugar;
387: }
388: }
389:
390: int sugard(MP m)
391: {
392: int s;
393:
394: for ( s = 0; m; m = NEXT(m) )
395: s = MAX(s,m->dl->td);
396: return s;
397: }
398:
399: void addd(VL vl,DP p1,DP p2,DP *pr)
400: {
401: int n;
402: MP m1,m2,mr=0,mr0;
403: Obj t;
404: DL d;
405:
406: if ( !p1 )
407: *pr = p2;
408: else if ( !p2 )
409: *pr = p1;
410: else {
411: if ( OID(p1) <= O_R ) {
412: n = NV(p2); NEWDL(d,n);
413: NEWMP(m1); m1->dl = d; C(m1) = (Obj)p1; NEXT(m1) = 0;
414: MKDP(n,m1,p1); (p1)->sugar = 0;
415: }
416: if ( OID(p2) <= O_R ) {
417: n = NV(p1); NEWDL(d,n);
418: NEWMP(m2); m2->dl = d; C(m2) = (Obj)p2; NEXT(m2) = 0;
419: MKDP(n,m2,p2); (p2)->sugar = 0;
420: }
421: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
422: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
423: case 0:
424: arf_add(vl,C(m1),C(m2),&t);
425: if ( t ) {
426: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = t;
427: }
428: m1 = NEXT(m1); m2 = NEXT(m2); break;
429: case 1:
430: NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1);
431: m1 = NEXT(m1); break;
432: case -1:
433: NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2);
434: m2 = NEXT(m2); break;
435: }
436: if ( !mr0 )
437: if ( m1 )
438: mr0 = m1;
439: else if ( m2 )
440: mr0 = m2;
441: else {
442: *pr = 0;
443: return;
444: }
445: else if ( m1 )
446: NEXT(mr) = m1;
447: else if ( m2 )
448: NEXT(mr) = m2;
449: else
450: NEXT(mr) = 0;
451: MKDP(NV(p1),mr0,*pr);
452: if ( *pr )
453: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
454: }
455: }
456:
457: /* for F4 symbolic reduction */
458:
459: void symb_addd(DP p1,DP p2,DP *pr)
460: {
461: int n;
462: MP m1,m2,mr=0,mr0;
463:
464: if ( !p1 )
465: *pr = p2;
466: else if ( !p2 )
467: *pr = p1;
468: else {
469: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
470: NEXTMP(mr0,mr); C(mr) = (Obj)ONE;
471: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
472: case 0:
473: mr->dl = m1->dl;
474: m1 = NEXT(m1); m2 = NEXT(m2); break;
475: case 1:
476: mr->dl = m1->dl;
477: m1 = NEXT(m1); break;
478: case -1:
479: mr->dl = m2->dl;
480: m2 = NEXT(m2); break;
481: }
482: }
483: if ( !mr0 )
484: if ( m1 )
485: mr0 = m1;
486: else if ( m2 )
487: mr0 = m2;
488: else {
489: *pr = 0;
490: return;
491: }
492: else if ( m1 )
493: NEXT(mr) = m1;
494: else if ( m2 )
495: NEXT(mr) = m2;
496: else
497: NEXT(mr) = 0;
498: MKDP(NV(p1),mr0,*pr);
499: if ( *pr )
500: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
501: }
502: }
503:
504: /*
505: * destructive merge of two list
506: *
507: * p1, p2 : list of DL
508: * return : a merged list
509: */
510:
511: NODE symb_merge(NODE m1,NODE m2,int n)
512: {
513: NODE top=0,prev,cur,m=0,t;
514: DL d1,d2;
515:
516: if ( !m1 )
517: return m2;
518: else if ( !m2 )
519: return m1;
520: else {
521: switch ( (*cmpdl)(n,(DL)BDY(m1),(DL)BDY(m2)) ) {
522: case 0:
523: top = m1; m = NEXT(m2);
524: break;
525: case 1:
526: top = m1; m = m2;
527: break;
528: case -1:
529: top = m2; m = m1;
530: break;
531: }
532: prev = top; cur = NEXT(top);
533: /* BDY(prev) > BDY(m) always holds */
534: while ( cur && m ) {
535: d1 = (DL)BDY(cur);
536: d2 = (DL)BDY(m);
537: #if 1
538: switch ( (*cmpdl)(n,(DL)BDY(cur),(DL)BDY(m)) ) {
539: #else
540: /* XXX only valid for DRL */
541: if ( d1->td > d2->td )
542: c = 1;
543: else if ( d1->td < d2->td )
544: c = -1;
545: else {
546: for ( i = n-1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
547: if ( i < 0 )
548: c = 0;
549: else if ( d1->d[i] < d2->d[i] )
550: c = 1;
551: else
552: c = -1;
553: }
554: switch ( c ) {
555: #endif
556: case 0:
557: m = NEXT(m);
558: prev = cur; cur = NEXT(cur);
559: break;
560: case 1:
561: t = NEXT(cur); NEXT(cur) = m; m = t;
562: prev = cur; cur = NEXT(cur);
563: break;
564: case -1:
565: NEXT(prev) = m; m = cur;
566: prev = NEXT(prev); cur = NEXT(prev);
567: break;
568: }
569: }
570: if ( !cur )
571: NEXT(prev) = m;
572: return top;
573: }
574: }
575:
576: void _adddl(int n,DL d1,DL d2,DL d3)
577: {
578: int i;
579:
580: d3->td = d1->td+d2->td;
581: for ( i = 0; i < n; i++ )
582: d3->d[i] = d1->d[i]+d2->d[i];
583: }
584:
1.3 noro 585: void _addtodl(int n,DL d1,DL d2)
586: {
587: int i;
588:
589: d2->td += d1->td;
590: for ( i = 0; i < n; i++ )
591: d2->d[i] += d1->d[i];
592: }
593:
1.23 ! noro 594: void _subfromdl(int n,DL d1,DL d2)
! 595: {
! 596: int i;
! 597:
! 598: d2->td -= d1->td;
! 599: for ( i = 0; i < n; i++ )
! 600: d2->d[i] -= d1->d[i];
! 601: }
! 602:
1.3 noro 603: void _copydl(int n,DL d1,DL d2)
604: {
605: int i;
606:
607: d2->td = d1->td;
608: for ( i = 0; i < n; i++ )
609: d2->d[i] = d1->d[i];
610: }
611:
612: int _eqdl(int n,DL d1,DL d2)
613: {
614: int i;
615:
616: if ( d2->td != d1->td ) return 0;
617: for ( i = 0; i < n; i++ )
618: if ( d2->d[i] != d1->d[i] ) return 0;
619: return 1;
620: }
621:
1.1 noro 622: /* m1 <- m1 U dl*f, destructive */
623:
624: NODE mul_dllist(DL dl,DP f);
625:
626: NODE symb_mul_merge(NODE m1,DL dl,DP f,int n)
627: {
628: NODE top,prev,cur,n1;
629: DP g;
630: DL t,s;
631: MP m;
632:
633: if ( !m1 )
634: return mul_dllist(dl,f);
635: else if ( !f )
636: return m1;
637: else {
638: m = BDY(f);
639: NEWDL_NOINIT(t,n);
640: _adddl(n,m->dl,dl,t);
641: top = m1; prev = 0; cur = m1;
642: while ( m ) {
643: switch ( (*cmpdl)(n,(DL)BDY(cur),t) ) {
644: case 0:
645: prev = cur; cur = NEXT(cur);
646: if ( !cur ) {
647: MKDP(n,m,g);
648: NEXT(prev) = mul_dllist(dl,g);
649: return top;
650: }
651: m = NEXT(m);
652: if ( m ) _adddl(n,m->dl,dl,t);
653: break;
654: case 1:
655: prev = cur; cur = NEXT(cur);
656: if ( !cur ) {
657: MKDP(n,m,g);
658: NEXT(prev) = mul_dllist(dl,g);
659: return top;
660: }
661: break;
662: case -1:
663: NEWDL_NOINIT(s,n);
664: s->td = t->td;
665: bcopy(t->d,s->d,n*sizeof(int));
666: NEWNODE(n1);
667: n1->body = (pointer)s;
668: NEXT(n1) = cur;
669: if ( !prev ) {
670: top = n1; cur = n1;
671: } else {
672: NEXT(prev) = n1; prev = n1;
673: }
674: m = NEXT(m);
675: if ( m ) _adddl(n,m->dl,dl,t);
676: break;
677: }
678: }
679: return top;
680: }
681: }
682:
683: DLBUCKET symb_merge_bucket(DLBUCKET m1,DLBUCKET m2,int n)
684: {
685: DLBUCKET top,prev,cur,m,t;
686:
687: if ( !m1 )
688: return m2;
689: else if ( !m2 )
690: return m1;
691: else {
692: if ( m1->td == m2->td ) {
693: top = m1;
694: BDY(top) = symb_merge(BDY(top),BDY(m2),n);
695: m = NEXT(m2);
696: } else if ( m1->td > m2->td ) {
697: top = m1; m = m2;
698: } else {
699: top = m2; m = m1;
700: }
701: prev = top; cur = NEXT(top);
702: /* prev->td > m->td always holds */
703: while ( cur && m ) {
704: if ( cur->td == m->td ) {
705: BDY(cur) = symb_merge(BDY(cur),BDY(m),n);
706: m = NEXT(m);
707: prev = cur; cur = NEXT(cur);
708: } else if ( cur->td > m->td ) {
709: t = NEXT(cur); NEXT(cur) = m; m = t;
710: prev = cur; cur = NEXT(cur);
711: } else {
712: NEXT(prev) = m; m = cur;
713: prev = NEXT(prev); cur = NEXT(prev);
714: }
715: }
716: if ( !cur )
717: NEXT(prev) = m;
718: return top;
719: }
720: }
721:
722: void subd(VL vl,DP p1,DP p2,DP *pr)
723: {
724: DP t;
725:
726: if ( !p2 )
727: *pr = p1;
728: else {
729: chsgnd(p2,&t); addd(vl,p1,t,pr);
730: }
731: }
732:
733: void chsgnd(DP p,DP *pr)
734: {
735: MP m,mr=0,mr0;
736: Obj r;
737:
738: if ( !p )
739: *pr = 0;
740: else if ( OID(p) <= O_R ) {
741: arf_chsgn((Obj)p,&r); *pr = (DP)r;
742: } else {
743: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
744: NEXTMP(mr0,mr); arf_chsgn(C(m),&C(mr)); mr->dl = m->dl;
745: }
746: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
747: if ( *pr )
748: (*pr)->sugar = p->sugar;
749: }
750: }
751:
752: void muld(VL vl,DP p1,DP p2,DP *pr)
753: {
754: if ( ! do_weyl )
755: comm_muld(vl,p1,p2,pr);
756: else
757: weyl_muld(vl,p1,p2,pr);
758: }
759:
760: void comm_muld(VL vl,DP p1,DP p2,DP *pr)
761: {
762: MP m;
763: DP s,t,u;
764: int i,l,l1;
765: static MP *w;
766: static int wlen;
767:
768: if ( !p1 || !p2 )
769: *pr = 0;
770: else if ( OID(p1) != O_DP )
771: muldc(vl,p2,(Obj)p1,pr);
772: else if ( OID(p2) != O_DP )
773: muldc(vl,p1,(Obj)p2,pr);
774: else {
775: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
776: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
777: if ( l1 < l ) {
778: t = p1; p1 = p2; p2 = t;
779: l = l1;
780: }
781: if ( l > wlen ) {
782: if ( w ) GCFREE(w);
783: w = (MP *)MALLOC(l*sizeof(MP));
784: wlen = l;
785: }
786: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
787: w[i] = m;
788: for ( s = 0, i = l-1; i >= 0; i-- ) {
789: muldm(vl,p1,w[i],&t); addd(vl,s,t,&u); s = u;
790: }
791: bzero(w,l*sizeof(MP));
792: *pr = s;
793: }
794: }
795:
796: /* discard terms which is not a multiple of dl */
797:
798: void comm_muld_trunc(VL vl,DP p1,DP p2,DL dl,DP *pr)
799: {
800: MP m;
801: DP s,t,u;
802: int i,l,l1;
803: static MP *w;
804: static int wlen;
805:
806: if ( !p1 || !p2 )
807: *pr = 0;
808: else if ( OID(p1) != O_DP )
809: muldc_trunc(vl,p2,(Obj)p1,dl,pr);
810: else if ( OID(p2) != O_DP )
811: muldc_trunc(vl,p1,(Obj)p2,dl,pr);
812: else {
813: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
814: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
815: if ( l1 < l ) {
816: t = p1; p1 = p2; p2 = t;
817: l = l1;
818: }
819: if ( l > wlen ) {
820: if ( w ) GCFREE(w);
821: w = (MP *)MALLOC(l*sizeof(MP));
822: wlen = l;
823: }
824: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
825: w[i] = m;
826: for ( s = 0, i = l-1; i >= 0; i-- ) {
827: muldm_trunc(vl,p1,w[i],dl,&t); addd(vl,s,t,&u); s = u;
828: }
829: bzero(w,l*sizeof(MP));
830: *pr = s;
831: }
832: }
833:
834: void comm_quod(VL vl,DP p1,DP p2,DP *pr)
835: {
836: MP m=0,m0;
837: DP s,t;
838: int i,n,sugar;
839: DL d1,d2,d;
840: Q a,b;
841:
842: if ( !p2 )
843: error("comm_quod : invalid input");
844: if ( !p1 )
845: *pr = 0;
846: else {
847: n = NV(p1);
848: d2 = BDY(p2)->dl;
849: m0 = 0;
850: sugar = p1->sugar;
851: while ( p1 ) {
852: d1 = BDY(p1)->dl;
853: NEWDL(d,n);
854: d->td = d1->td - d2->td;
855: for ( i = 0; i < n; i++ )
856: d->d[i] = d1->d[i]-d2->d[i];
857: NEXTMP(m0,m);
858: m->dl = d;
859: divq((Q)BDY(p1)->c,(Q)BDY(p2)->c,&a); chsgnq(a,&b);
860: C(m) = (Obj)b;
861: muldm_trunc(vl,p2,m,d2,&t);
862: addd(vl,p1,t,&s); p1 = s;
863: C(m) = (Obj)a;
864: }
865: if ( m0 ) {
866: NEXT(m) = 0; MKDP(n,m0,*pr);
867: } else
868: *pr = 0;
869: /* XXX */
870: if ( *pr )
871: (*pr)->sugar = sugar - d2->td;
872: }
873: }
874:
875: void muldm(VL vl,DP p,MP m0,DP *pr)
876: {
877: MP m,mr=0,mr0;
878: Obj c;
879: DL d;
880: int n;
881:
882: if ( !p )
883: *pr = 0;
884: else {
885: for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl, n = NV(p);
886: m; m = NEXT(m) ) {
887: NEXTMP(mr0,mr);
888: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
889: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
890: else
891: arf_mul(vl,C(m),c,&C(mr));
892: adddl(n,m->dl,d,&mr->dl);
893: }
894: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
895: if ( *pr )
896: (*pr)->sugar = p->sugar + m0->dl->td;
897: }
898: }
899:
900: void muldm_trunc(VL vl,DP p,MP m0,DL dl,DP *pr)
901: {
902: MP m,mr=0,mr0;
903: Obj c;
904: DL d,tdl;
905: int n,i;
906:
907: if ( !p )
908: *pr = 0;
909: else {
910: n = NV(p);
911: NEWDL(tdl,n);
912: for ( mr0 = 0, m = BDY(p), c = C(m0), d = m0->dl;
913: m; m = NEXT(m) ) {
914: _adddl(n,m->dl,d,tdl);
915: for ( i = 0; i < n; i++ )
916: if ( tdl->d[i] < dl->d[i] )
917: break;
918: if ( i < n )
919: continue;
920: NEXTMP(mr0,mr);
921: mr->dl = tdl;
922: NEWDL(tdl,n);
923: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
924: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
925: else
926: arf_mul(vl,C(m),(Obj)c,&C(mr));
927: }
928: if ( mr0 ) {
929: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
930: } else
931: *pr = 0;
932: if ( *pr )
933: (*pr)->sugar = p->sugar + m0->dl->td;
934: }
935: }
936:
937: void weyl_muld(VL vl,DP p1,DP p2,DP *pr)
938: {
939: MP m;
940: DP s,t,u;
941: int i,l;
942: static MP *w;
943: static int wlen;
944:
945: if ( !p1 || !p2 )
946: *pr = 0;
947: else if ( OID(p1) != O_DP )
948: muldc(vl,p2,(Obj)p1,pr);
949: else if ( OID(p2) != O_DP )
950: muldc(vl,p1,(Obj)p2,pr);
951: else {
952: for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
953: if ( l > wlen ) {
954: if ( w ) GCFREE(w);
955: w = (MP *)MALLOC(l*sizeof(MP));
956: wlen = l;
957: }
958: for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
959: w[i] = m;
960: for ( s = 0, i = l-1; i >= 0; i-- ) {
961: weyl_muldm(vl,w[i],p2,&t); addd(vl,s,t,&u); s = u;
962: }
963: bzero(w,l*sizeof(MP));
964: *pr = s;
965: }
966: }
967:
968: void actm(VL vl,int nv,MP m1,MP m2,DP *pr)
969: {
970: DL d1,d2,d;
971: int n2,i,j,k;
972: Z jq,c,c1;
973: MP m;
974: Obj t;
975:
976: d1 = m1->dl;
977: d2 = m2->dl;
978: for ( i = 0; i < nv; i++ )
979: if ( d1->d[i] > d2->d[i] ) {
980: *pr = 0; return;
981: }
982: NEWDL(d,nv);
983: c = ONE;
984: for ( i = 0; i < nv; i++ ) {
985: for ( j = d2->d[i], k = d1->d[i]; k > 0; k--, j-- ) {
1.2 noro 986: STOZ(j,jq); mulz(c,jq,&c1); c = c1;
1.1 noro 987: }
988: d->d[i] = d2->d[i]-d1->d[i];
989: }
990: arf_mul(vl,C(m1),C(m2),&t);
991: NEWMP(m);
992: arf_mul(vl,(Obj)c,t,&C(m));
993: m->dl = d;
994: MKDP(nv,m,*pr);
995: }
996:
997: void weyl_actd(VL vl,DP p1,DP p2,DP *pr)
998: {
999: int n;
1000: MP m1,m2;
1001: DP d,r,s;
1002:
1003: if ( !p1 || !p2 ) *pr = 0;
1004: else {
1005: n = NV(p1);
1006: r = 0;
1007: for ( m1 = BDY(p1); m1; m1 = NEXT(m1) )
1008: for ( m2 = BDY(p2); m2; m2 = NEXT(m2) ) {
1009: actm(vl,n,m1,m2,&d);
1010: addd(vl,r,d,&s); r = s;
1011: }
1012: *pr = r;
1013: }
1014: }
1015:
1016: /* monomial * polynomial */
1017:
1018: void weyl_muldm(VL vl,MP m0,DP p,DP *pr)
1019: {
1020: DP r,t,t1;
1021: MP m;
1022: DL d0;
1023: int n,n2,l,i,j,tlen;
1024: static MP *w,*psum;
1025: static struct cdl *tab;
1026: static int wlen;
1027: static int rtlen;
1028:
1029: if ( !p )
1030: *pr = 0;
1031: else {
1032: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
1033: if ( l > wlen ) {
1034: if ( w ) GCFREE(w);
1035: w = (MP *)MALLOC(l*sizeof(MP));
1036: wlen = l;
1037: }
1038: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
1039: w[i] = m;
1040:
1041: n = NV(p); n2 = n>>1;
1042: d0 = m0->dl;
1043: for ( i = 0, tlen = 1; i < n2; i++ )
1044: tlen *= d0->d[n2+i]+1;
1045: if ( tlen > rtlen ) {
1046: if ( tab ) GCFREE(tab);
1047: if ( psum ) GCFREE(psum);
1048: rtlen = tlen;
1049: tab = (struct cdl *)MALLOC(rtlen*sizeof(struct cdl));
1050: psum = (MP *)MALLOC(rtlen*sizeof(MP));
1051: }
1052: bzero(psum,tlen*sizeof(MP));
1053: for ( i = l-1; i >= 0; i-- ) {
1054: bzero(tab,tlen*sizeof(struct cdl));
1055: weyl_mulmm(vl,m0,w[i],n,tab,tlen);
1056: for ( j = 0; j < tlen; j++ ) {
1057: if ( tab[j].c ) {
1058: NEWMP(m); m->dl = tab[j].d; C(m) = (Obj)tab[j].c; NEXT(m) = psum[j];
1059: psum[j] = m;
1060: }
1061: }
1062: }
1063: for ( j = tlen-1, r = 0; j >= 0; j-- )
1064: if ( psum[j] ) {
1065: MKDP(n,psum[j],t); addd(vl,r,t,&t1); r = t1;
1066: }
1067: if ( r )
1068: r->sugar = p->sugar + m0->dl->td;
1069: *pr = r;
1070: }
1071: }
1072:
1073: /* m0 = x0^d0*x1^d1*... * dx0^e0*dx1^e1*... */
1074: /* rtab : array of length (e0+1)*(e1+1)*... */
1075:
1076: void weyl_mulmm(VL vl,MP m0,MP m1,int n,struct cdl *rtab,int rtablen)
1077: {
1078: Obj c,c0,c1;
1079: DL d,d0,d1,dt;
1080: int i,j,a,b,k,l,n2,s,min,curlen;
1081: struct cdl *p;
1082: static Z *ctab;
1083: static struct cdl *tab;
1084: static int tablen;
1085: static struct cdl *tmptab;
1086: static int tmptablen;
1087:
1088:
1089: if ( !m0 || !m1 ) {
1090: rtab[0].c = 0;
1091: rtab[0].d = 0;
1092: return;
1093: }
1094: c0 = C(m0); c1 = C(m1);
1095: arf_mul(vl,c0,c1,&c);
1096: d0 = m0->dl; d1 = m1->dl;
1097: n2 = n>>1;
1098: curlen = 1;
1099: NEWDL(d,n);
1100: if ( n & 1 )
1101: /* offset of h-degree */
1102: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
1103: else
1104: d->td = 0;
1105: rtab[0].c = c;
1106: rtab[0].d = d;
1107:
1108: if ( rtablen > tmptablen ) {
1109: if ( tmptab ) GCFREE(tmptab);
1110: tmptab = (struct cdl *)MALLOC(rtablen*sizeof(struct cdl));
1111: tmptablen = rtablen;
1112: }
1113: for ( i = 0; i < n2; i++ ) {
1114: a = d0->d[i]; b = d1->d[n2+i];
1115: k = d0->d[n2+i]; l = d1->d[i];
1116:
1117: /* degree of xi^a*(Di^k*xi^l)*Di^b */
1118: a += l;
1119: b += k;
1120: s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i);
1121:
1122: if ( !k || !l ) {
1123: for ( j = 0, p = rtab; j < curlen; j++, p++ ) {
1124: if ( p->c ) {
1125: dt = p->d;
1126: dt->d[i] = a;
1127: dt->d[n2+i] = b;
1128: dt->td += s;
1129: }
1130: }
1131: curlen *= k+1;
1132: continue;
1133: }
1134: if ( k+1 > tablen ) {
1135: if ( tab ) GCFREE(tab);
1136: if ( ctab ) GCFREE(ctab);
1137: tablen = k+1;
1138: tab = (struct cdl *)MALLOC(tablen*sizeof(struct cdl));
1139: ctab = (Z *)MALLOC(tablen*sizeof(Q));
1140: }
1141: /* compute xi^a*(Di^k*xi^l)*Di^b */
1142: min = MIN(k,l);
1143: mkwc(k,l,ctab);
1144: bzero(tab,(k+1)*sizeof(struct cdl));
1145: if ( n & 1 )
1146: for ( j = 0; j <= min; j++ ) {
1147: NEWDL(d,n);
1148: d->d[i] = a-j; d->d[n2+i] = b-j;
1149: d->td = s;
1150: d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i));
1151: tab[j].d = d;
1152: tab[j].c = (Obj)ctab[j];
1153: }
1154: else
1155: for ( j = 0; j <= min; j++ ) {
1156: NEWDL(d,n);
1157: d->d[i] = a-j; d->d[n2+i] = b-j;
1158: d->td = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i); /* XXX */
1159: tab[j].d = d;
1160: tab[j].c = (Obj)ctab[j];
1161: }
1162: bzero(ctab,(min+1)*sizeof(Q));
1163: comm_muld_tab(vl,n,rtab,curlen,tab,k+1,tmptab);
1164: curlen *= k+1;
1165: bcopy(tmptab,rtab,curlen*sizeof(struct cdl));
1166: }
1167: }
1168:
1169: /* direct product of two cdl tables
1170: rt[] = [
1171: t[0]*t1[0],...,t[n-1]*t1[0],
1172: t[0]*t1[1],...,t[n-1]*t1[1],
1173: ...
1174: t[0]*t1[n1-1],...,t[n-1]*t1[n1-1]
1175: ]
1176: */
1177:
1178: void comm_muld_tab(VL vl,int nv,struct cdl *t,int n,struct cdl *t1,int n1,struct cdl *rt)
1179: {
1180: int i,j;
1181: struct cdl *p;
1182: Obj c;
1183: DL d;
1184:
1185: bzero(rt,n*n1*sizeof(struct cdl));
1186: for ( j = 0, p = rt; j < n1; j++ ) {
1187: c = (Obj)t1[j].c;
1188: d = t1[j].d;
1189: if ( !c )
1190: break;
1191: for ( i = 0; i < n; i++, p++ ) {
1192: if ( t[i].c ) {
1193: arf_mul(vl,(Obj)t[i].c,c,(Obj *)&p->c);
1194: adddl(nv,t[i].d,d,&p->d);
1195: }
1196: }
1197: }
1198: }
1199:
1200: void muldc(VL vl,DP p,Obj c,DP *pr)
1201: {
1202: MP m,mr=0,mr0;
1203:
1204: if ( !p || !c )
1205: *pr = 0;
1206: else if ( NUM(c) && UNIQ((Q)c) )
1207: *pr = p;
1208: else if ( NUM(c) && MUNIQ((Q)c) )
1209: chsgnd(p,pr);
1210: else {
1211: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1212: NEXTMP(mr0,mr);
1213: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
1214: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
1215: else
1216: arf_mul(vl,C(m),c,&C(mr));
1217: mr->dl = m->dl;
1218: }
1219: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
1220: if ( *pr )
1221: (*pr)->sugar = p->sugar;
1222: }
1223: }
1224:
1225: void divdc(VL vl,DP p,Obj c,DP *pr)
1226: {
1227: Obj inv;
1228:
1229: arf_div(vl,(Obj)ONE,c,&inv);
1230: muld(vl,p,(DP)inv,pr);
1231: }
1232:
1233: void muldc_trunc(VL vl,DP p,Obj c,DL dl,DP *pr)
1234: {
1235: MP m,mr=0,mr0;
1236: DL mdl;
1237: int i,n;
1238:
1239: if ( !p || !c ) {
1240: *pr = 0; return;
1241: }
1242: n = NV(p);
1243: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1244: mdl = m->dl;
1245: for ( i = 0; i < n; i++ )
1246: if ( mdl->d[i] < dl->d[i] )
1247: break;
1248: if ( i < n )
1249: break;
1250: NEXTMP(mr0,mr);
1251: if ( NUM(C(m)) && RATN(C(m)) && NUM(c) && RATN(c) )
1252: mulq((Q)C(m),(Q)c,(Q *)&C(mr));
1253: else
1254: arf_mul(vl,C(m),c,&C(mr));
1255: mr->dl = m->dl;
1256: }
1257: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
1258: if ( *pr )
1259: (*pr)->sugar = p->sugar;
1260: }
1261:
1262: void divsdc(VL vl,DP p,P c,DP *pr)
1263: {
1264: MP m,mr=0,mr0;
1265:
1266: if ( !c )
1267: error("disvsdc : division by 0");
1268: else if ( !p )
1269: *pr = 0;
1270: else if ( OID(c) > O_P )
1271: error("divsdc : invalid argument");
1272: else {
1273: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
1274: NEXTMP(mr0,mr); divsp(vl,(P)C(m),c,(P *)&C(mr)); mr->dl = m->dl;
1275: }
1276: NEXT(mr) = 0; MKDP(NV(p),mr0,*pr);
1277: if ( *pr )
1278: (*pr)->sugar = p->sugar;
1279: }
1280: }
1281:
1282: void adddl(int n,DL d1,DL d2,DL *dr)
1283: {
1284: DL dt;
1285: int i;
1286:
1287: *dr = dt = (DL)MALLOC_ATOMIC((n+1)*sizeof(int));
1288: dt->td = d1->td + d2->td;
1289: for ( i = 0; i < n; i++ )
1290: dt->d[i] = d1->d[i]+d2->d[i];
1291: }
1292:
1293: /* d1 += d2 */
1294:
1295: void adddl_destructive(int n,DL d1,DL d2)
1296: {
1297: int i;
1298:
1299: d1->td += d2->td;
1300: for ( i = 0; i < n; i++ )
1301: d1->d[i] += d2->d[i];
1302: }
1303:
1304: int compd(VL vl,DP p1,DP p2)
1305: {
1306: int n,t;
1307: MP m1,m2;
1308:
1309: if ( !p1 )
1310: return p2 ? -1 : 0;
1311: else if ( !p2 )
1312: return 1;
1313: else if ( NV(p1) != NV(p2) ) {
1314: error("compd : size mismatch");
1315: return 0; /* XXX */
1316: } else {
1317: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
1318: m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
1319: if ( (t = (*cmpdl)(n,m1->dl,m2->dl)) ||
1320: (t = arf_comp(vl,C(m1),C(m2)) ) )
1321: return t;
1322: if ( m1 )
1323: return 1;
1324: else if ( m2 )
1325: return -1;
1326: else
1327: return 0;
1328: }
1329: }
1330:
1331: int cmpdl_lex(int n,DL d1,DL d2)
1332: {
1333: int i;
1334:
1335: for ( i = 0; i < n && d1->d[i] == d2->d[i]; i++ );
1336: return i == n ? 0 : (d1->d[i] > d2->d[i] ? 1 : -1);
1337: }
1338:
1339: int cmpdl_revlex(int n,DL d1,DL d2)
1340: {
1341: int i;
1342:
1343: for ( i = n - 1; i >= 0 && d1->d[i] == d2->d[i]; i-- );
1344: return i < 0 ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1345: }
1346:
1347: int cmpdl_gradlex(int n,DL d1,DL d2)
1348: {
1349: if ( d1->td > d2->td )
1350: return 1;
1351: else if ( d1->td < d2->td )
1352: return -1;
1353: else
1354: return cmpdl_lex(n,d1,d2);
1355: }
1356:
1357: int cmpdl_revgradlex(int n,DL d1,DL d2)
1358: {
1359: register int i,c;
1360: register int *p1,*p2;
1361:
1362: if ( d1->td > d2->td )
1363: return 1;
1364: else if ( d1->td < d2->td )
1365: return -1;
1366: else {
1367: i = n-1;
1368: p1 = d1->d+n-1;
1369: p2 = d2->d+n-1;
1370: while ( i >= 7 ) {
1371: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1372: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1373: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1374: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1375: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1376: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1377: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1378: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1379: i -= 8;
1380: }
1381: switch ( i ) {
1382: case 6:
1383: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1384: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1385: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1386: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1387: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1388: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1389: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1390: return 0;
1391: case 5:
1392: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1393: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1394: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1395: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1396: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1397: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1398: return 0;
1399: case 4:
1400: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1401: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1402: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1403: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1404: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1405: return 0;
1406: case 3:
1407: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1408: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1409: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1410: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1411: return 0;
1412: case 2:
1413: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1414: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1415: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1416: return 0;
1417: case 1:
1418: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1419: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1420: return 0;
1421: case 0:
1422: c = (*p1--) - (*p2--); if ( c ) goto LAST;
1423: return 0;
1424: default:
1425: return 0;
1426: }
1427: LAST:
1428: if ( c > 0 ) return -1;
1429: else return 1;
1430: }
1431: }
1432:
1433: int cmpdl_blex(int n,DL d1,DL d2)
1434: {
1435: int c;
1436:
1437: if ( (c = cmpdl_lex(n-1,d1,d2)) )
1438: return c;
1439: else {
1440: c = d1->d[n-1] - d2->d[n-1];
1441: return c > 0 ? 1 : c < 0 ? -1 : 0;
1442: }
1443: }
1444:
1445: int cmpdl_bgradlex(int n,DL d1,DL d2)
1446: {
1447: int e1,e2,c;
1448:
1449: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
1450: if ( e1 > e2 )
1451: return 1;
1452: else if ( e1 < e2 )
1453: return -1;
1454: else {
1455: c = cmpdl_lex(n-1,d1,d2);
1456: if ( c )
1457: return c;
1458: else
1459: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
1460: }
1461: }
1462:
1463: int cmpdl_brevgradlex(int n,DL d1,DL d2)
1464: {
1465: int e1,e2,c;
1466:
1467: e1 = d1->td - d1->d[n-1]; e2 = d2->td - d2->d[n-1];
1468: if ( e1 > e2 )
1469: return 1;
1470: else if ( e1 < e2 )
1471: return -1;
1472: else {
1473: c = cmpdl_revlex(n-1,d1,d2);
1474: if ( c )
1475: return c;
1476: else
1477: return d1->td > d2->td ? 1 : d1->td < d2->td ? -1 : 0;
1478: }
1479: }
1480:
1481: int cmpdl_brevrev(int n,DL d1,DL d2)
1482: {
1483: int e1,e2,f1,f2,c,i;
1484:
1485: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1486: e1 += d1->d[i]; e2 += d2->d[i];
1487: }
1488: f1 = d1->td - e1; f2 = d2->td - e2;
1489: if ( e1 > e2 )
1490: return 1;
1491: else if ( e1 < e2 )
1492: return -1;
1493: else {
1494: c = cmpdl_revlex(dp_nelim,d1,d2);
1495: if ( c )
1496: return c;
1497: else if ( f1 > f2 )
1498: return 1;
1499: else if ( f1 < f2 )
1500: return -1;
1501: else {
1502: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1503: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1504: }
1505: }
1506: }
1507:
1508: int cmpdl_bgradrev(int n,DL d1,DL d2)
1509: {
1510: int e1,e2,f1,f2,c,i;
1511:
1512: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1513: e1 += d1->d[i]; e2 += d2->d[i];
1514: }
1515: f1 = d1->td - e1; f2 = d2->td - e2;
1516: if ( e1 > e2 )
1517: return 1;
1518: else if ( e1 < e2 )
1519: return -1;
1520: else {
1521: c = cmpdl_lex(dp_nelim,d1,d2);
1522: if ( c )
1523: return c;
1524: else if ( f1 > f2 )
1525: return 1;
1526: else if ( f1 < f2 )
1527: return -1;
1528: else {
1529: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1530: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1531: }
1532: }
1533: }
1534:
1535: int cmpdl_blexrev(int n,DL d1,DL d2)
1536: {
1537: int e1,e2,f1,f2,c,i;
1538:
1539: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1540: e1 += d1->d[i]; e2 += d2->d[i];
1541: }
1542: f1 = d1->td - e1; f2 = d2->td - e2;
1543: c = cmpdl_lex(dp_nelim,d1,d2);
1544: if ( c )
1545: return c;
1546: else if ( f1 > f2 )
1547: return 1;
1548: else if ( f1 < f2 )
1549: return -1;
1550: else {
1551: for ( i = n - 1; i >= dp_nelim && d1->d[i] == d2->d[i]; i-- );
1552: return i < dp_nelim ? 0 : (d1->d[i] < d2->d[i] ? 1 : -1);
1553: }
1554: }
1555:
1556: int cmpdl_elim(int n,DL d1,DL d2)
1557: {
1558: int e1,e2,i;
1559:
1560: for ( i = 0, e1 = 0, e2 = 0; i < dp_nelim; i++ ) {
1561: e1 += d1->d[i]; e2 += d2->d[i];
1562: }
1563: if ( e1 > e2 )
1564: return 1;
1565: else if ( e1 < e2 )
1566: return -1;
1567: else
1568: return cmpdl_revgradlex(n,d1,d2);
1569: }
1570:
1571: int cmpdl_weyl_elim(int n,DL d1,DL d2)
1572: {
1573: int e1,e2,i;
1574:
1575: for ( i = 1, e1 = 0, e2 = 0; i <= dp_nelim; i++ ) {
1576: e1 += d1->d[n-i]; e2 += d2->d[n-i];
1577: }
1578: if ( e1 > e2 )
1579: return 1;
1580: else if ( e1 < e2 )
1581: return -1;
1582: else if ( d1->td > d2->td )
1583: return 1;
1584: else if ( d1->td < d2->td )
1585: return -1;
1586: else return -cmpdl_revlex(n,d1,d2);
1587: }
1588:
1589: /*
1590: a special ordering
1591: 1. total order
1592: 2. (-w,w) for the first 2*m variables
1593: 3. DRL for the first 2*m variables
1594: */
1595:
1596: extern int *current_weyl_weight_vector;
1597:
1598: int cmpdl_homo_ww_drl(int n,DL d1,DL d2)
1599: {
1600: int e1,e2,m,i;
1601: int *p1,*p2;
1602:
1603: if ( d1->td > d2->td )
1604: return 1;
1605: else if ( d1->td < d2->td )
1606: return -1;
1607:
1608: m = n>>1;
1609: for ( i = 0, e1 = e2 = 0, p1 = d1->d, p2 = d2->d; i < m; i++ ) {
1610: e1 += current_weyl_weight_vector[i]*(p1[m+i] - p1[i]);
1611: e2 += current_weyl_weight_vector[i]*(p2[m+i] - p2[i]);
1612: }
1613: if ( e1 > e2 )
1614: return 1;
1615: else if ( e1 < e2 )
1616: return -1;
1617:
1618: e1 = d1->td - d1->d[n-1];
1619: e2 = d2->td - d2->d[n-1];
1620: if ( e1 > e2 )
1621: return 1;
1622: else if ( e1 < e2 )
1623: return -1;
1624:
1625: for ( i= n - 1, p1 = d1->d+n-1, p2 = d2->d+n-1;
1626: i >= 0 && *p1 == *p2; i--, p1--, p2-- );
1627: return i < 0 ? 0 : (*p1 < *p2 ? 1 : -1);
1628: }
1629:
1630: int cmpdl_drl_zigzag(int n,DL d1,DL d2)
1631: {
1632: int i,t,m;
1633: int *p1,*p2;
1634:
1635: if ( d1->td > d2->td )
1636: return 1;
1637: else if ( d1->td < d2->td )
1638: return -1;
1639: else {
1640: m = n>>1;
1641: for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) {
1642: if ( (t = p1[m+i] - p2[m+i]) ) return t > 0 ? -1 : 1;
1643: if ( (t = p1[i] - p2[i]) ) return t > 0 ? -1 : 1;
1644: }
1645: return 0;
1646: }
1647: }
1648:
1649: int cmpdl_homo_ww_drl_zigzag(int n,DL d1,DL d2)
1650: {
1651: int e1,e2,m,i,t;
1652: int *p1,*p2;
1653:
1654: if ( d1->td > d2->td )
1655: return 1;
1656: else if ( d1->td < d2->td )
1657: return -1;
1658:
1659: m = n>>1;
1660: for ( i = 0, e1 = e2 = 0, p1 = d1->d, p2 = d2->d; i < m; i++ ) {
1661: e1 += current_weyl_weight_vector[i]*(p1[m+i] - p1[i]);
1662: e2 += current_weyl_weight_vector[i]*(p2[m+i] - p2[i]);
1663: }
1664: if ( e1 > e2 )
1665: return 1;
1666: else if ( e1 < e2 )
1667: return -1;
1668:
1669: e1 = d1->td - d1->d[n-1];
1670: e2 = d2->td - d2->d[n-1];
1671: if ( e1 > e2 )
1672: return 1;
1673: else if ( e1 < e2 )
1674: return -1;
1675:
1676: for ( i= m - 1, p1 = d1->d, p2 = d2->d; i >= 0; i-- ) {
1677: if ( (t = p1[m+i] - p2[m+i]) ) return t > 0 ? -1 : 1;
1678: if ( (t = p1[i] - p2[i]) ) return t > 0 ? -1 : 1;
1679: }
1680: return 0;
1681: }
1682:
1683: int cmpdl_order_pair(int n,DL d1,DL d2)
1684: {
1685: int e1,e2,i,j,l;
1686: int *t1,*t2;
1687: int len,head;
1688: struct order_pair *pair;
1689:
1690: len = dp_current_spec->ord.block.length;
1691: if ( n != dp_current_spec->nv )
1692: error("cmpdl_order_pair : incompatible order specification");
1693: pair = dp_current_spec->ord.block.order_pair;
1694:
1695: head = 0;
1696: for ( i = 0, t1 = d1->d, t2 = d2->d; i < len; i++ ) {
1697: l = pair[i].length;
1698: switch ( pair[i].order ) {
1699: case 0:
1700: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
1701: e1 += MUL_WEIGHT(t1[j],head+j);
1702: e2 += MUL_WEIGHT(t2[j],head+j);
1703: }
1704: if ( e1 > e2 )
1705: return 1;
1706: else if ( e1 < e2 )
1707: return -1;
1708: else {
1709: for ( j = l - 1; j >= 0 && t1[j] == t2[j]; j-- );
1710: if ( j >= 0 )
1711: return t1[j] < t2[j] ? 1 : -1;
1712: }
1713: break;
1714: case 1:
1715: for ( j = 0, e1 = e2 = 0; j < l; j++ ) {
1716: e1 += MUL_WEIGHT(t1[j],head+j);
1717: e2 += MUL_WEIGHT(t2[j],head+j);
1718: }
1719: if ( e1 > e2 )
1720: return 1;
1721: else if ( e1 < e2 )
1722: return -1;
1723: else {
1724: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
1725: if ( j < l )
1726: return t1[j] > t2[j] ? 1 : -1;
1727: }
1728: break;
1729: case 2:
1730: for ( j = 0; j < l && t1[j] == t2[j]; j++ );
1731: if ( j < l )
1732: return t1[j] > t2[j] ? 1 : -1;
1733: break;
1734: default:
1735: error("cmpdl_order_pair : invalid order"); break;
1736: }
1737: t1 += l; t2 += l; head += l;
1738: }
1739: return 0;
1740: }
1741:
1742: int cmpdl_composite(int nv,DL d1,DL d2)
1743: {
1744: int n,i,j,k,start,s,len;
1745: int *dw;
1746: struct sparse_weight *sw;
1747: struct weight_or_block *worb;
1748: int *w,*t1,*t2;
1749:
1750: n = dp_current_spec->ord.composite.length;
1751: worb = dp_current_spec->ord.composite.w_or_b;
1752: w = dp_dl_work;
1753: for ( i = 0, t1 = d1->d, t2 = d2->d; i < nv; i++ )
1754: w[i] = t1[i]-t2[i];
1755: for ( i = 0; i < n; i++, worb++ ) {
1756: len = worb->length;
1757: switch ( worb->type ) {
1758: case IS_DENSE_WEIGHT:
1759: dw = worb->body.dense_weight;
1760: for ( j = 0, s = 0; j < len; j++ )
1761: s += dw[j]*w[j];
1762: if ( s > 0 ) return 1;
1763: else if ( s < 0 ) return -1;
1764: break;
1765: case IS_SPARSE_WEIGHT:
1766: sw = worb->body.sparse_weight;
1767: for ( j = 0, s = 0; j < len; j++ )
1768: s += sw[j].value*w[sw[j].pos];
1769: if ( s > 0 ) return 1;
1770: else if ( s < 0 ) return -1;
1771: break;
1772: case IS_BLOCK:
1773: start = worb->body.block.start;
1774: switch ( worb->body.block.order ) {
1775: case 0:
1776: for ( j = 0, k = start, s = 0; j < len; j++, k++ ) {
1777: s += MUL_WEIGHT(w[k],k);
1778: }
1779: if ( s > 0 ) return 1;
1780: else if ( s < 0 ) return -1;
1781: else {
1782: for ( j = k-1; j >= start && w[j] == 0; j-- );
1783: if ( j >= start )
1784: return w[j] < 0 ? 1 : -1;
1785: }
1786: break;
1787: case 1:
1788: for ( j = 0, k = start, s = 0; j < len; j++, k++ ) {
1789: s += MUL_WEIGHT(w[k],k);
1790: }
1791: if ( s > 0 ) return 1;
1792: else if ( s < 0 ) return -1;
1793: else {
1794: for ( j = 0, k = start; j < len && w[j] == 0; j++, k++ );
1795: if ( j < len )
1796: return w[j] > 0 ? 1 : -1;
1797: }
1798: break;
1799: case 2:
1800: for ( j = 0, k = start; j < len && w[j] == 0; j++, k++ );
1801: if ( j < len )
1802: return w[j] > 0 ? 1 : -1;
1803: break;
1804: }
1805: break;
1806: }
1807: }
1808: return 0;
1809: }
1810:
1811: int cmpdl_matrix(int n,DL d1,DL d2)
1812: {
1813: int *v,*w,*t1,*t2;
1814: int s,i,j,len;
1815: int **matrix;
1816:
1817: for ( i = 0, t1 = d1->d, t2 = d2->d, w = dp_dl_work; i < n; i++ )
1818: w[i] = t1[i]-t2[i];
1819: len = dp_current_spec->ord.matrix.row;
1820: matrix = dp_current_spec->ord.matrix.matrix;
1821: for ( j = 0; j < len; j++ ) {
1822: v = matrix[j];
1823: for ( i = 0, s = 0; i < n; i++ )
1824: s += v[i]*w[i];
1825: if ( s > 0 )
1826: return 1;
1827: else if ( s < 0 )
1828: return -1;
1829: }
1830: return 0;
1831: }
1832:
1833: int cmpdl_top_weight(int n,DL d1,DL d2)
1834: {
1835: int *w;
1836: Z **mat;
1837: Z *a;
1838: mpz_t sum;
1839: int len,i,sgn,tsgn,row,k;
1840: int *t1,*t2;
1841:
1842: w = (int *)ALLOCA(n*sizeof(int));
1843: len = current_top_weight_len+3;
1844: t1 = d1->d; t2 = d2->d;
1845: for ( i = 0; i < n; i++ ) w[i] = t1[i]-t2[i];
1846: mpz_init_set_ui(sum,0);
1847: if ( OID(current_top_weight) == O_VECT ) {
1848: mat = (Z **)&BDY((VECT)current_top_weight);
1849: row = 1;
1850: } else {
1851: mat = (Z **)BDY((MAT)current_top_weight);
1852: row = ((MAT)current_top_weight)->row;
1853: }
1854: for ( k = 0; k < row; k++ ) {
1855: a = mat[k];
1856: for ( i = 0; i < n; i++ ) {
1857: if ( !a[i] || !w[i] ) continue;
1858: if ( w[i] > 0 )
1859: mpz_addmul_ui(sum,BDY(a[i]),(unsigned int)w[i]);
1860: else
1861: mpz_submul_ui(sum,BDY(a[i]),(unsigned int)(-w[i]));
1862: }
1863: sgn = mpz_sgn(sum);
1864: if ( sgn > 0 ) return 1;
1865: else if ( sgn < 0 ) return -1;
1866: }
1867: return (*cmpdl_tie_breaker)(n,d1,d2);
1868: }
1869:
1870: GeoBucket create_bucket()
1871: {
1872: GeoBucket g;
1873:
1874: g = CALLOC(1,sizeof(struct oGeoBucket));
1875: g->m = 32;
1876: return g;
1877: }
1878:
1879: int length(NODE d);
1880:
1881: void add_bucket(GeoBucket g,NODE d,int nv)
1882: {
1883: int l,k,m;
1884:
1885: l = length(d);
1886: for ( k = 0, m = 1; l > m; k++, m <<= 1 );
1887: /* 2^(k-1) < l <= 2^k */
1888: d = symb_merge(g->body[k],d,nv);
1889: for ( ; length(d) > (1<<(k)); k++ ) {
1890: g->body[k] = 0;
1891: d = symb_merge(g->body[k+1],d,nv);
1892: }
1893: g->body[k] = d;
1894: g->m = MAX(g->m,k);
1895: }
1896:
1897: DL remove_head_bucket(GeoBucket g,int nv)
1898: {
1899: int j,i,c,m;
1900: DL d;
1901:
1902: j = -1;
1903: m = g->m;
1904: for ( i = 0; i <= m; i++ ) {
1905: if ( !g->body[i] )
1906: continue;
1907: if ( j < 0 ) j = i;
1908: else {
1909: c = (*cmpdl)(nv,g->body[i]->body,g->body[j]->body);
1910: if ( c > 0 )
1911: j = i;
1912: else if ( c == 0 )
1913: g->body[i] = NEXT(g->body[i]);
1914: }
1915: }
1916: if ( j < 0 )
1917: return 0;
1918: else {
1919: d = g->body[j]->body;
1920: g->body[j] = NEXT(g->body[j]);
1921: return d;
1922: }
1923: }
1924:
1925: /* DPV functions */
1926:
1927: void adddv(VL vl,DPV p1,DPV p2,DPV *pr)
1928: {
1929: int i,len;
1930: DP *e;
1931:
1932: if ( !p1 || !p2 )
1933: error("adddv : invalid argument");
1934: else if ( p1->len != p2->len )
1935: error("adddv : size mismatch");
1936: else {
1937: len = p1->len;
1938: e = (DP *)MALLOC(p1->len*sizeof(DP));
1939: for ( i = 0; i < len; i++ )
1940: addd(vl,p1->body[i],p2->body[i],&e[i]);
1941: MKDPV(len,e,*pr);
1942: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
1943: }
1944: }
1945:
1946: void subdv(VL vl,DPV p1,DPV p2,DPV *pr)
1947: {
1948: int i,len;
1949: DP *e;
1950:
1951: if ( !p1 || !p2 )
1952: error("subdv : invalid argument");
1953: else if ( p1->len != p2->len )
1954: error("subdv : size mismatch");
1955: else {
1956: len = p1->len;
1957: e = (DP *)MALLOC(p1->len*sizeof(DP));
1958: for ( i = 0; i < len; i++ )
1959: subd(vl,p1->body[i],p2->body[i],&e[i]);
1960: MKDPV(len,e,*pr);
1961: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
1962: }
1963: }
1964:
1965: void chsgndv(DPV p1,DPV *pr)
1966: {
1967: int i,len;
1968: DP *e;
1969:
1970: if ( !p1 )
1971: error("subdv : invalid argument");
1972: else {
1973: len = p1->len;
1974: e = (DP *)MALLOC(p1->len*sizeof(DP));
1975: for ( i = 0; i < len; i++ )
1976: chsgnd(p1->body[i],&e[i]);
1977: MKDPV(len,e,*pr);
1978: (*pr)->sugar = p1->sugar;
1979: }
1980: }
1981:
1982: void muldv(VL vl,DP p1,DPV p2,DPV *pr)
1983: {
1984: int i,len;
1985: DP *e;
1986:
1987: len = p2->len;
1988: e = (DP *)MALLOC(p2->len*sizeof(DP));
1989: if ( !p1 ) {
1990: MKDPV(len,e,*pr);
1991: (*pr)->sugar = 0;
1992: } else {
1993: for ( i = 0; i < len; i++ )
1994: muld(vl,p1,p2->body[i],&e[i]);
1995: MKDPV(len,e,*pr);
1996: (*pr)->sugar = p1->sugar + p2->sugar;
1997: }
1998: }
1999:
2000: int compdv(VL vl,DPV p1,DPV p2)
2001: {
2002: int i,t,len;
2003:
2004: if ( p1->len != p2->len ) {
2005: error("compdv : size mismatch");
2006: return 0; /* XXX */
2007: } else {
2008: len = p1->len;
2009: for ( i = 0; i < len; i++ )
2010: if ( (t = compd(vl,p1->body[i],p2->body[i])) )
2011: return t;
2012: return 0;
2013: }
2014: }
2015:
2016: int ni_next(int *a,int n)
2017: {
2018: int i,j,k,kj;
2019:
2020: /* find the first nonzero a[j] */
2021: for ( j = 0; j < n && a[j] == 0; j++ );
2022: /* find the first zero a[k] after a[j] */
2023: for ( k = j; k < n && a[k] == 1; k++ );
2024: if ( k == n ) return 0;
2025: /* a[0] = 0, ... , a[j-1] = 0, a[j] = 1, ..., a[k-1] = 1, a[k] = 0 */
2026: /* a[0] = 1,..., a[k-j-2] = 1, a[k-j-1] = 0, ..., a[k-1] = 0, a[k] = 1 */
2027: kj = k-j-1;
2028: for ( i = 0; i < kj; i++ ) a[i] = 1;
2029: for ( ; i < k; i++ ) a[i] = 0;
2030: a[k] = 1;
2031: return 1;
2032: }
2033:
2034: int comp_nbm(NBM a,NBM b)
2035: {
2036: int d,i,ai,bi;
2037: unsigned int *ab,*bb;
2038:
2039: if ( a->d > b->d ) return 1;
2040: else if ( a->d < b->d ) return -1;
2041: else {
2042: d = a->d; ab = a->b; bb = b->b;
2043: #if 0
2044: w = (d+31)/32;
2045: for ( i = 0; i < w; i++ )
2046: if ( ab[i] > bb[i] ) return 1;
2047: else if ( ab[i] < bb[i] ) return -1;
2048: #else
2049: for ( i = 0; i < d; i++ ) {
2050: ai = NBM_GET(ab,i);
2051: bi = NBM_GET(bb,i);
2052: if ( ai > bi ) return 1;
2053: else if ( ai < bi ) return -1;
2054: }
2055: #endif
2056: return 0;
2057: }
2058: }
2059:
2060: NBM mul_nbm(NBM a,NBM b)
2061: {
2062: int ad,bd,d,i,j;
2063: unsigned int *ab,*bb,*mb;
2064: NBM m;
2065:
2066: ad = a->d; bd = b->d; ab = a->b; bb = b->b;
2067: d = ad + bd;
2068: NEWNBM(m); NEWNBMBDY(m,d);
2069: m->d = d; mulp(CO,a->c,b->c,&m->c); mb = m->b;
2070: j = 0;
2071: for ( i = 0; i < ad; i++, j++ )
2072: if ( NBM_GET(ab,i) ) NBM_SET(mb,j);
2073: else NBM_CLR(mb,j);
2074: for ( i = 0; i < bd; i++, j++ )
2075: if ( NBM_GET(bb,i) ) NBM_SET(mb,j);
2076: else NBM_CLR(mb,j);
2077: return m;
2078: }
2079:
2080: NBP nbmtonbp(NBM m)
2081: {
2082: NODE n;
2083: NBP u;
2084:
2085: MKNODE(n,m,0);
2086: MKNBP(u,n);
2087: return u;
2088: }
2089:
2090: /* a=c*x*rest -> a0= x*rest, ah=x, ar=rest */
2091:
2092: P separate_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar)
2093: {
2094: int i,d1;
2095: NBM t;
2096:
2097: if ( !a->d ) error("separate_nbm : invalid argument");
2098:
2099: if ( a0 ) {
2100: NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
2101: *a0 = nbmtonbp(t);
2102: }
2103:
2104: if ( ah ) {
2105: NEWNBM(t); NEWNBMBDY(t,1); t->d = 1; t->c = (P)ONE;
2106: if ( NBM_GET(a->b,0) ) NBM_SET(t->b,0);
2107: else NBM_CLR(t->b,0);
2108: *ah = nbmtonbp(t);
2109: }
2110:
2111: if ( ar ) {
2112: d1 = a->d-1;
2113: NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
2114: for ( i = 0; i < d1; i++ ) {
2115: if ( NBM_GET(a->b,i+1) ) NBM_SET(t->b,i);
2116: else NBM_CLR(t->b,i);
2117: }
2118: *ar = nbmtonbp(t);
2119: }
2120:
2121: return a->c;
2122: }
2123:
2124: /* a=c*rest*x -> a0= rest*x, ar=rest, at=x */
2125:
2126: P separate_tail_nbm(NBM a,NBP *a0,NBP *ar,NBP *at)
2127: {
2128: int i,d,d1;
2129: NBM t;
2130:
2131: if ( !(d=a->d) ) error("separate_tail_nbm : invalid argument");
2132:
2133: if ( a0 ) {
2134: NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
2135: *a0 = nbmtonbp(t);
2136: }
2137:
2138: d1 = a->d-1;
2139: if ( at ) {
2140: NEWNBM(t); NEWNBMBDY(t,1); t->d = 1; t->c = (P)ONE;
2141: if ( NBM_GET(a->b,d1) ) NBM_SET(t->b,0);
2142: else NBM_CLR(t->b,0);
2143: *at = nbmtonbp(t);
2144: }
2145:
2146: if ( ar ) {
2147: NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
2148: for ( i = 0; i < d1; i++ ) {
2149: if ( NBM_GET(a->b,i) ) NBM_SET(t->b,i);
2150: else NBM_CLR(t->b,i);
2151: }
2152: *ar = nbmtonbp(t);
2153: }
2154:
2155: return a->c;
2156: }
2157:
2158: NBP make_xky(int k)
2159: {
2160: int k1,i;
2161: NBM t;
2162:
2163: NEWNBM(t); NEWNBMBDY(t,k); t->d = k; t->c = (P)ONE;
2164: k1 = k-1;
2165: for ( i = 0; i < k1; i++ ) NBM_SET(t->b,i);
2166: NBM_CLR(t->b,i);
2167: return nbmtonbp(t);
2168: }
2169:
2170: /* a=c*x^(k-1)*y*rest -> a0= x^(k-1)*y*rest, ah=x^(k-1)*y, ar=rest */
2171:
2172: P separate_xky_nbm(NBM a,NBP *a0,NBP *ah,NBP *ar)
2173: {
2174: int i,d1,k,k1;
2175: NBM t;
2176:
2177: if ( !a->d )
2178: error("separate_nbm : invalid argument");
2179: for ( i = 0; i < a->d && NBM_GET(a->b,i); i++ );
2180: if ( i == a->d )
2181: error("separate_nbm : invalid argument");
2182: k1 = i;
2183: k = i+1;
2184:
2185: if ( a0 ) {
2186: NEWNBM(t); t->d = a->d; t->b = a->b; t->c = (P)ONE;
2187: *a0 = nbmtonbp(t);
2188: }
2189:
2190: if ( ah ) {
2191: NEWNBM(t); NEWNBMBDY(t,k); t->d = k; t->c = (P)ONE;
2192: for ( i = 0; i < k1; i++ ) NBM_SET(t->b,i);
2193: NBM_CLR(t->b,i);
2194: *ah = nbmtonbp(t);
2195: }
2196:
2197: if ( ar ) {
2198: d1 = a->d-k;
2199: NEWNBM(t); NEWNBMBDY(t,d1); t->d = d1; t->c = (P)ONE;
2200: for ( i = 0; i < d1; i++ ) {
2201: if ( NBM_GET(a->b,i+k) ) NBM_SET(t->b,i);
2202: else NBM_CLR(t->b,i);
2203: }
2204: *ar = nbmtonbp(t);
2205: }
2206:
2207: return a->c;
2208: }
2209:
2210: void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp);
2211: void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp);
2212: void mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp);
2213: void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp);
2214:
2215: NBP shuffle_mul_nbm(NBM a,NBM b)
2216: {
2217: NBP u,a0,ah,ar,b0,bh,br,a1,b1,t;
2218: P ac,bc,c;
2219:
2220: if ( !a->d || !b->d )
2221: u = nbmtonbp(mul_nbm(a,b));
2222: else {
2223: ac = separate_nbm(a,&a0,&ah,&ar);
2224: bc = separate_nbm(b,&b0,&bh,&br);
2225: mulp(CO,ac,bc,&c);
2226: shuffle_mulnbp(CO,ar,b0,&t); mulnbp(CO,ah,t,&a1);
2227: shuffle_mulnbp(CO,a0,br,&t); mulnbp(CO,bh,t,&b1);
2228: addnbp(CO,a1,b1,&t); mulnbp(CO,(NBP)c,t,&u);
2229: }
2230: return u;
2231: }
2232:
2233: NBP harmonic_mul_nbm(NBM a,NBM b)
2234: {
2235: NBP u,a0,ah,ar,b0,bh,br,a1,b1,t,s,abk,ab1;
2236: P ac,bc,c;
2237:
2238: if ( !a->d || !b->d )
2239: u = nbmtonbp(mul_nbm(a,b));
2240: else {
2241: mulp(CO,a->c,b->c,&c);
2242: ac = separate_xky_nbm(a,&a0,&ah,&ar);
2243: bc = separate_xky_nbm(b,&b0,&bh,&br);
2244: mulp(CO,ac,bc,&c);
2245: harmonic_mulnbp(CO,ar,b0,&t); mulnbp(CO,ah,t,&a1);
2246: harmonic_mulnbp(CO,a0,br,&t); mulnbp(CO,bh,t,&b1);
2247: abk = make_xky(((NBM)BDY(BDY(ah)))->d+((NBM)BDY(BDY(bh)))->d);
2248: harmonic_mulnbp(CO,ar,br,&t); mulnbp(CO,abk,t,&ab1);
2249: addnbp(CO,a1,b1,&t); addnbp(CO,t,ab1,&s); mulnbp(CO,(NBP)c,s,&u);
2250: }
2251: return u;
2252:
2253: }
2254:
2255: void addnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2256: {
2257: NODE b1,b2,br=0,br0;
2258: NBM m1,m2,m;
2259: P c;
2260:
2261: if ( !p1 )
2262: *rp = p2;
2263: else if ( !p2 )
2264: *rp = p1;
2265: else {
2266: for ( b1 = BDY(p1), b2 = BDY(p2), br0 = 0; b1 && b2; ) {
2267: m1 = (NBM)BDY(b1); m2 = (NBM)BDY(b2);
2268: switch ( comp_nbm(m1,m2) ) {
2269: case 0:
2270: addp(CO,m1->c,m2->c,&c);
2271: if ( c ) {
2272: NEXTNODE(br0,br);
2273: NEWNBM(m); m->d = m1->d; m->c = c; m->b = m1->b;
2274: BDY(br) = (pointer)m;
2275: }
2276: b1 = NEXT(b1); b2 = NEXT(b2); break;
2277: case 1:
2278: NEXTNODE(br0,br); BDY(br) = BDY(b1);
2279: b1 = NEXT(b1); break;
2280: case -1:
2281: NEXTNODE(br0,br); BDY(br) = BDY(b2);
2282: b2 = NEXT(b2); break;
2283: }
2284: }
2285: if ( !br0 )
2286: if ( b1 )
2287: br0 = b1;
2288: else if ( b2 )
2289: br0 = b2;
2290: else {
2291: *rp = 0;
2292: return;
2293: }
2294: else if ( b1 )
2295: NEXT(br) = b1;
2296: else if ( b2 )
2297: NEXT(br) = b2;
2298: else
2299: NEXT(br) = 0;
2300: MKNBP(*rp,br0);
2301: }
2302: }
2303:
2304: void subnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2305: {
2306: NBP t;
2307:
2308: chsgnnbp(p2,&t);
2309: addnbp(vl,p1,t,rp);
2310: }
2311:
2312: void chsgnnbp(NBP p,NBP *rp)
2313: {
2314: NODE r0,r=0,b;
2315: NBM m,m1;
2316:
2317: for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
2318: NEXTNODE(r0,r);
2319: m = (NBM)BDY(b);
2320: NEWNBM(m1); m1->d = m->d; m1->b = m->b; chsgnp(m->c,&m1->c);
2321: BDY(r) = m1;
2322: }
2323: if ( r0 ) NEXT(r) = 0;
2324: MKNBP(*rp,r0);
2325: }
2326:
2327: void mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2328: {
2329: NODE b,n;
2330: NBP r,t,s;
2331: NBM m;
2332:
2333: if ( !p1 || !p2 ) {
2334: *rp = 0; return;
2335: }
2336: if ( OID(p1) != O_NBP ) {
2337: if ( !POLY(p1) )
2338: error("mulnbp : invalid argument");
2339: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
2340: MKNODE(n,m,0); MKNBP(p1,n);
2341: }
2342: if ( OID(p2) != O_NBP ) {
2343: if ( !POLY(p2) )
2344: error("mulnbp : invalid argument");
2345: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
2346: MKNODE(n,m,0); MKNBP(p2,n);
2347: }
2348: if ( length(BDY(p1)) < length(BDY(p2)) ) {
2349: for ( r = 0, b = BDY(p1); b; b = NEXT(b) ) {
2350: mulnbmnbp(vl,(NBM)BDY(b),p2,&t);
2351: addnbp(vl,r,t,&s); r = s;
2352: }
2353: *rp = r;
2354: } else {
2355: for ( r = 0, b = BDY(p2); b; b = NEXT(b) ) {
2356: mulnbpnbm(vl,p1,(NBM)BDY(b),&t);
2357: addnbp(vl,r,t,&s); r = s;
2358: }
2359: *rp = r;
2360: }
2361: }
2362:
2363: void mulnbmnbp(VL vl,NBM m,NBP p, NBP *rp)
2364: {
2365: NODE b,r0,r=0;
2366:
2367: if ( !p ) *rp = 0;
2368: else {
2369: for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
2370: NEXTNODE(r0,r);
2371: BDY(r) = mul_nbm(m,(NBM)BDY(b));
2372: }
2373: if ( r0 ) NEXT(r) = 0;
2374: MKNBP(*rp,r0);
2375: }
2376: }
2377:
2378: void mulnbpnbm(VL vl,NBP p,NBM m, NBP *rp)
2379: {
2380: NODE b,r0,r=0;
2381:
2382: if ( !p ) *rp = 0;
2383: else {
2384: for ( r0 = 0, b = BDY(p); b; b = NEXT(b) ) {
2385: NEXTNODE(r0,r);
2386: BDY(r) = mul_nbm((NBM)BDY(b),m);
2387: }
2388: if ( r0 ) NEXT(r) = 0;
2389: MKNBP(*rp,r0);
2390: }
2391: }
2392:
2393: void pwrnbp(VL vl,NBP a,Z q,NBP *c)
2394: {
2395: NBP a1,a2;
2396: Z q1,r1,two;
2397: NBM m;
2398: NODE r;
2399:
2400: if ( !q ) {
2401: NEWNBM(m); m->d = 0; m->c = (P)ONE; m->b = 0;
2402: MKNODE(r,m,0); MKNBP(*c,r);
2403: } else if ( !a )
2404: *c = 0;
2405: else if ( UNIQ(q) )
2406: *c = a;
2407: else {
1.2 noro 2408: STOZ(2,two);
1.1 noro 2409: divqrz(q,two,&q1,&r1);
2410: pwrnbp(vl,a,q1,&a1);
2411: mulnbp(vl,a1,a1,&a2);
2412: if ( r1 )
2413: mulnbp(vl,a2,a,c);
2414: else
2415: *c = a2;
2416: }
2417: }
2418:
2419: int compnbp(VL vl,NBP p1,NBP p2)
2420: {
2421: NODE n1,n2;
2422: NBM m1,m2;
2423: int t;
2424:
2425: if ( !p1 )
2426: return p2 ? -1 : 0;
2427: else if ( !p2 )
2428: return 1;
2429: else {
2430: for ( n1 = BDY(p1), n2 = BDY(p2);
2431: n1 && n2; n1 = NEXT(n1), n2 = NEXT(n2) ) {
2432: m1 = (NBM)BDY(n1); m2 = (NBM)BDY(n2);
2433: if ( (t = comp_nbm(m1,m2)) || (t = compp(CO,m1->c,m2->c) ) )
2434: return t;
2435: }
2436: if ( n1 )
2437: return 1;
2438: else if ( n2 )
2439: return -1;
2440: else
2441: return 0;
2442: }
2443: }
2444:
2445: void shuffle_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2446: {
2447: NODE b1,b2,n;
2448: NBP r,t,s;
2449: NBM m;
2450:
2451: if ( !p1 || !p2 ) {
2452: *rp = 0; return;
2453: }
2454: if ( OID(p1) != O_NBP ) {
2455: if ( !POLY(p1) )
2456: error("shuffle_mulnbp : invalid argument");
2457: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
2458: MKNODE(n,m,0); MKNBP(p1,n);
2459: }
2460: if ( OID(p2) != O_NBP ) {
2461: if ( !POLY(p2) )
2462: error("shuffle_mulnbp : invalid argument");
2463: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
2464: MKNODE(n,m,0); MKNBP(p2,n);
2465: }
2466: for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) )
2467: for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) {
2468: t = shuffle_mul_nbm(m,(NBM)BDY(b2));
2469: addnbp(vl,r,t,&s); r = s;
2470: }
2471: *rp = r;
2472: }
2473:
2474: void harmonic_mulnbp(VL vl,NBP p1,NBP p2, NBP *rp)
2475: {
2476: NODE b1,b2,n;
2477: NBP r,t,s;
2478: NBM m;
2479:
2480: if ( !p1 || !p2 ) {
2481: *rp = 0; return;
2482: }
2483: if ( OID(p1) != O_NBP ) {
2484: if ( !POLY(p1) )
2485: error("harmonic_mulnbp : invalid argument");
2486: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p1;
2487: MKNODE(n,m,0); MKNBP(p1,n);
2488: }
2489: if ( OID(p2) != O_NBP ) {
2490: if ( !POLY(p2) )
2491: error("harmonic_mulnbp : invalid argument");
2492: NEWNBM(m); m->d = 0; m->b = 0; m->c = (P)p2;
2493: MKNODE(n,m,0); MKNBP(p2,n);
2494: }
2495: for ( r = 0, b1 = BDY(p1); b1; b1 = NEXT(b1) )
2496: for ( m = BDY(b1), b2 = BDY(p2); b2; b2 = NEXT(b2) ) {
2497: t = harmonic_mul_nbm(m,(NBM)BDY(b2));
2498: addnbp(vl,r,t,&s); r = s;
2499: }
2500: *rp = r;
2501: }
2502:
2503: #if 0
2504: NBP shuffle_mul_nbm(NBM a,NBM b)
2505: {
2506: int ad,bd,d,i,ai,bi,bit,s;
2507: int *ab,*bb,*wmb,*w;
2508: NBM wm,tm;
2509: P c,c1;
2510: NODE r,t,t1,p;
2511: NBP u;
2512:
2513: ad = a->d; bd = b->d; ab = a->b; bb = b->b;
2514: d = ad + bd;
2515: w = (int *)ALLOCA(d*sizeof(int));
2516: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2517: for ( i = 0; i < ad; i++ ) w[i] = 1;
2518: for ( ; i < d; i++ ) w[i] = 0;
2519: mulp(CO,a->c,b->c,&c);
2520: r = 0;
2521: do {
2522: wm->d = d; wm->c = c;
2523: ai = 0; bi = 0;
2524: for ( i = 0; i < d; i++ ) {
2525: if ( w[i] ) { bit = NBM_GET(ab,ai); ai++; }
2526: else { bit = NBM_GET(bb,bi); bi++; }
2527: if ( bit ) NBM_SET(wmb,i);
2528: else NBM_CLR(wmb,i);
2529: }
2530: for ( p = 0, t = r; t; p = t, t = NEXT(t) ) {
2531: tm = (NBM)BDY(t);
2532: s = comp_nbm(tm,wm);
2533: if ( s < 0 ) {
2534: /* insert */
2535: MKNODE(t1,wm,t);
2536: if ( !p ) r = t1;
2537: else NEXT(p) = t1;
2538: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2539: break;
2540: } else if ( s == 0 ) {
2541: /* add coefs */
2542: addp(CO,tm->c,c,&c1);
2543: if ( c1 ) tm->c = c1;
2544: else NEXT(p) = NEXT(t);
2545: break;
2546: }
2547: }
2548: if ( !t ) {
2549: /* append */
2550: MKNODE(t1,wm,t);
2551: if ( !p ) r = t1;
2552: else NEXT(p) = t1;
2553: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2554: }
2555: } while ( ni_next(w,d) );
2556: MKNBP(u,r);
2557: return u;
2558: }
2559:
2560: int nbmtoxky(NBM a,int *b)
2561: {
2562: int d,i,j,k;
2563: int *p;
2564:
2565: d = a->d; p = a->b;
2566: for ( i = j = 0, k = 1; i < d; i++ ) {
2567: if ( !NBM_GET(p,i) ) {
2568: b[j++] = k;
2569: k = 1;
2570: } else k++;
2571: }
2572: return j;
2573: }
2574:
2575: NBP harmonic_mul_nbm(NBM a,NBM b)
2576: {
2577: int da,db,d,la,lb,lmax,lmin,l,lab,la1,lb1,lab1;
2578: int i,j,k,ia,ib,s;
2579: int *wa,*wb,*w,*wab,*wa1,*wmb;
2580: P c,c1;
2581: NBM wm,tm;
2582: NODE r,t1,t,p;
2583: NBP u;
2584:
2585: da = a->d; db = b->d; d = da+db;
2586: wa = (int *)ALLOCA(da*sizeof(int));
2587: wb = (int *)ALLOCA(db*sizeof(int));
2588: la = nbmtoxky(a,wa);
2589: lb = nbmtoxky(b,wb);
2590: mulp(CO,a->c,b->c,&c);
2591: /* wa[0],..,wa[la-1] <-> x^wa[0]y x^wa[1]y .. */
2592: /* lmax : total length */
2593: lmax = la+lb;
2594: lmin = la>lb?la:lb;
2595: w = (int *)ALLOCA(lmax*sizeof(int));
2596: /* position of a+b */
2597: wab = (int *)ALLOCA(lmax*sizeof(int));
2598: /* position of a */
2599: wa1 = (int *)ALLOCA(lmax*sizeof(int));
2600: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2601: for ( l = lmin, r = 0; l <= lmax; l++ ) {
2602: lab = lmax - l;
2603: la1 = la - lab;
2604: lb1 = lb - lab;
2605: lab1 = l-lab;
2606: /* partion l into three parts: a, b, a+b */
2607: /* initialize wab */
2608: for ( i = 0; i < lab; i++ ) wab[i] = 1;
2609: for ( ; i < l; i++ ) wab[i] = 0;
2610: do {
2611: /* initialize wa1 */
2612: for ( i = 0; i < la1; i++ ) wa1[i] = 1;
2613: for ( ; i < lab1; i++ ) wa1[i] = 0;
2614: do {
2615: ia = 0; ib = 0;
2616: for ( i = j = 0; i < l; i++ )
2617: if ( wab[i] ) w[i] = wa[ia++]+wb[ib++];
2618: else if ( wa1[j++] ) w[i] = wa[ia++];
2619: else w[i] = wb[ib++];
2620: for ( i = j = 0; i < l; i++ ) {
2621: for ( k = w[i]-1; k > 0; k--, j++ ) NBM_SET(wmb,j);
2622: NBM_CLR(wmb,j); j++;
2623: }
2624: wm->d = j; wm->c = c;
2625: for ( p = 0, t = r; t; p = t, t = NEXT(t) ) {
2626: tm = (NBM)BDY(t);
2627: s = comp_nbm(tm,wm);
2628: if ( s < 0 ) {
2629: /* insert */
2630: MKNODE(t1,wm,t);
2631: if ( !p ) r = t1;
2632: else NEXT(p) = t1;
2633: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2634: break;
2635: } else if ( s == 0 ) {
2636: /* add coefs */
2637: addp(CO,tm->c,c,&c1);
2638: if ( c1 ) tm->c = c1;
2639: else NEXT(p) = NEXT(t);
2640: break;
2641: }
2642: }
2643: if ( !t ) {
2644: /* append */
2645: MKNODE(t1,wm,t);
2646: if ( !p ) r = t1;
2647: else NEXT(p) = t1;
2648: NEWNBM(wm); NEWNBMBDY(wm,d); wmb = wm->b;
2649: }
2650: } while ( ni_next(wa1,lab1) );
2651: } while ( ni_next(wab,l) );
2652: }
2653: MKNBP(u,r);
2654: return u;
2655: }
2656: #endif
2657:
2658: /* DPM functions */
2659:
1.3 noro 2660: DMMstack dmm_stack;
1.9 noro 2661: int dpm_base_ordtype;;
1.3 noro 2662:
1.9 noro 2663: DMMstack push_schreyer_order(LIST data,DMMstack stack)
1.3 noro 2664: {
2665: DMMstack t;
1.14 noro 2666: DP dp;
2667: MP mp;
2668: DMM m0,m1;
2669: DPM dpm0,dpm1;
2670: int len,i,nv;
1.7 noro 2671: NODE in,t1;
1.9 noro 2672: LIST l;
1.4 noro 2673:
2674: /* data = [DPM,...,DPM] */
1.9 noro 2675: if ( !dmm_stack && ( !dp_current_spec || dp_current_spec->id < 256 ) )
2676: error("push_schreyer_order : base module order is not set");
1.4 noro 2677: in = BDY(data);
2678: len = length(in);
2679: NEWDMMstack(t);
2680: t->rank = len;
2681: t->in = (DMM *)MALLOC((len+1)*sizeof(DMM));
1.14 noro 2682: t->sum = (DMM *)MALLOC((len+1)*sizeof(DMM));
1.9 noro 2683: if ( stack ) {
2684: MKNODE(t1,data,BDY(stack->obj)); MKLIST(l,t1); t->obj = l;
1.14 noro 2685: for ( i = 1; i <= len; i++, in = NEXT(in) ) {
2686: m1 = t->in[i] = BDY((DPM)BDY(in));
2687: NEWMP(mp); mp->dl = m1->dl; mp->c = m1->c; NEXT(mp) = 0;
2688: nv = ((DPM)BDY(in))->nv;
2689: MKDP(nv,mp,dp); dp->sugar = mp->dl->td;
2690: m0 = stack->sum[m1->pos]; MKDPM(nv,m0,dpm0);
2691: mulobjdpm(CO,(Obj)dp,dpm0,&dpm1);
2692: t->sum[i] = BDY(dpm1);
2693: }
1.9 noro 2694: } else {
2695: MKNODE(t1,data,0); MKLIST(l,t1); t->obj = l;
1.14 noro 2696: for ( i = 1; i <= len; i++, in = NEXT(in) ) {
2697: t->sum[i] = t->in[i] = BDY((DPM)BDY(in));
2698: }
1.4 noro 2699: }
1.9 noro 2700: t->next = stack;;
2701: dpm_ordtype = 3;
2702: return t;
1.4 noro 2703: }
2704:
2705: // data=[Ink,...,In0]
2706: // Ini = a list of module monomials
2707:
2708: void set_schreyer_order(LIST data)
2709: {
2710: NODE in;
2711: LIST *w;
2712: int i,len;
1.3 noro 2713:
2714: if ( !data ) {
2715: dmm_stack = 0;
2716: if ( dp_current_spec && dp_current_spec->id >= 256 )
1.9 noro 2717: dpm_ordtype = dp_current_spec->module_ordtype;
1.3 noro 2718: else
2719: dpm_ordtype = 0;
2720: return;
2721: } else {
1.9 noro 2722: if ( !dp_current_spec || dp_current_spec->id < 256 )
2723: error("set_schreyer_order : base module order is not set");
1.4 noro 2724: dmm_stack = 0;
1.9 noro 2725: dpm_base_ordtype = dp_current_spec->module_ordtype;
1.4 noro 2726: in = BDY(data);
1.3 noro 2727: len = length(in);
1.4 noro 2728: w = (LIST *)MALLOC(len*sizeof(LIST));
2729: for ( i = 0; i < len; i++, in = NEXT(in) ) w[i] = (LIST)BDY(in);
1.9 noro 2730: for ( i = len-1; i >= 0; i-- ) dmm_stack = push_schreyer_order(w[i],dmm_stack);
1.3 noro 2731: }
2732: }
2733:
1.15 noro 2734: void set_schreyer_level(DMMstack_array array,int level)
2735: {
2736: if ( !level ) {
2737: dmm_stack = 0;
2738: if ( dp_current_spec && dp_current_spec->id >= 256 )
2739: dpm_ordtype = dp_current_spec->module_ordtype;
2740: else
2741: dpm_ordtype = 0;
2742: return;
2743: } else {
2744: if ( !dp_current_spec || dp_current_spec->id < 256 )
2745: error("set_schreyer_level : base module order is not set");
2746: if ( level > array->len )
2747: error("set_schreyer_level : invalid level");
2748: dmm_stack = array->body[level-1];
2749: dpm_base_ordtype = dp_current_spec->module_ordtype;
2750: dpm_ordtype = 3;
2751: }
2752: }
2753:
1.4 noro 2754: // construct a base of syz(g)
2755: // assuming the schrerer order is properly set
2756:
2757: DP dpm_sp_hm(DPM p1,DPM p2);
2758: void dpm_sp(DPM p1,DPM p2,DPM *sp,DP *t1,DP *t2);
1.12 noro 2759: DPM dpm_nf_and_quotient3(DPM sp,VECT psv,DPM *nf,P *dn);
2760: DPM dpm_nf_and_quotient4(DPM sp,DPM *ps,VECT psiv,DPM head,DPM *nf,P *dn);
2761: DPM dpm_sp_nf(VECT psv,VECT psiv,int i,int j,DPM *nf);
1.19 noro 2762: DPM dpm_sp_nf_zlist(VECT psv,VECT psiv,int i,int j,int top,DPM *nf);
1.13 noro 2763: DPM dpm_sp_nf_asir(VECT psv,int i,int j,DPM *nf);
1.4 noro 2764: void dpm_sort(DPM p,DPM *r);
2765:
2766: extern int DP_Multiple;
1.22 noro 2767: extern int DP_Print;
1.4 noro 2768:
2769: void dpm_nf_z(NODE b,DPM g,VECT psv,int full,int multiple,DPM *rp);
2770: NODE dpm_sort_list(NODE l);
2771: void dpm_ptozp(DPM p,Z *cont,DPM *r);
2772:
2773: NODE dpm_reduceall(NODE in)
2774: {
2775: int n,i;
2776: VECT psv;
2777: DPM *ps;
2778: NODE t,t1;
2779: DPM g,r;
2780: Z cont;
2781:
2782: n = length(in);
2783: MKVECT(psv,n);
2784: ps = (DPM *)BDY(psv);
2785: for ( i = 0, t = in; i < n; i++, t = NEXT(t) ) ps[i] = BDY(t);
2786: for ( i = 0; i < n; i++ ) {
2787: g = ps[i]; ps[i] = 0;
1.8 noro 2788: // dpm_nf_z(0,g,psv,1,DP_Multiple,&ps[i]);
2789: dpm_nf_z(0,g,psv,1,0,&ps[i]);
1.4 noro 2790: }
2791: t = 0;
2792: for ( i = n-1; i >= 0; i-- ) {
2793: dpm_ptozp(ps[i],&cont,&r);
2794: MKNODE(t1,r,t); t = t1;
2795: }
2796: return t;
2797: }
2798:
1.8 noro 2799: struct oEGT egra;
2800:
1.9 noro 2801: void dpm_ht(DPM d,DPM *r);
1.16 noro 2802: NODE reverse_node(NODE b);
1.9 noro 2803:
1.12 noro 2804: void dpm_schreyer_base(LIST g,LIST *s)
2805: {
2806: NODE nd,t0,t,b0,b;
2807: int n,i,j,k,nv,max,pos;
2808: LIST l;
2809: DP h,t1,t2;
2810: MP d;
2811: DMM r0,r,r1;
2812: DPM sp,nf,dpm;
2813: DPM *ps;
2814: VECT psv,psiv;
2815: DPM quo;
1.16 noro 2816: DP *mm;
1.12 noro 2817: NODE *psi;
1.13 noro 2818: NODE n1,n2,n3;
2819: int p1,p2,p3;
2820: struct oEGT eg0,eg1,egsp,egnf;
1.12 noro 2821: extern struct oEGT egred;
1.15 noro 2822: extern int sch_count,schrec_count,schlast_count;
1.12 noro 2823:
1.15 noro 2824: sch_count = schlast_count= 0;
1.12 noro 2825: init_eg(&egra);
1.13 noro 2826: init_eg(&egsp);
2827: init_eg(&egnf);
1.12 noro 2828: nd = BDY(g);
2829: n = length(nd);
2830: MKVECT(psv,n+1);
2831: ps = (DPM *)BDY(psv);
2832: for ( i = 1, t = nd; i <= n; i++, t = NEXT(t) ) ps[i] = (DPM)BDY(t);
2833: for ( i = 1, max = 0; i <= n; i++ )
2834: if ( (pos=BDY(ps[i])->pos) > max ) max = pos;
2835: MKVECT(psiv,max+1);
2836: psi = (NODE *)BDY(psiv);
2837: for ( i = n; i >= 1; i-- ) {
2838: pos = BDY(ps[i])->pos;
2839: MKNODE(nd,(long)i,psi[pos]); psi[pos] = nd;
2840: }
2841: nv = ps[1]->nv;
1.16 noro 2842: mm = (DP *)MALLOC((n+1)*sizeof(DP));
1.12 noro 2843: b0 = 0;
1.13 noro 2844: get_eg(&eg0);
2845: for ( i = 1; i <= max; i++ ) {
1.16 noro 2846: memset(mm,0,(n+1)*sizeof(DP));
1.13 noro 2847: for ( n1 = psi[i]; n1; n1 = NEXT(n1) ) {
2848: p1 = (long)BDY(n1);
2849: for ( n2 = NEXT(n1); n2; n2 = NEXT(n2) ) {
2850: p2 = (long)BDY(n2);
1.16 noro 2851: mm[p2] = dpm_sp_hm(ps[p1],ps[p2]);
1.13 noro 2852: }
2853: for ( n2 = NEXT(n1); n2; n2 = NEXT(n2) ) {
2854: p2 = (long)BDY(n2);
1.16 noro 2855: if ( !mm[p2] ) continue;
2856: for ( h = mm[p2], n3 = NEXT(n1); n3; n3 = NEXT(n3) ) {
1.13 noro 2857: p3 = (long)BDY(n3);
1.16 noro 2858: if ( n3 != n2 && mm[p3] && dp_redble(mm[p3],h) ) mm[p3] = 0;
2859: }
2860: }
2861: for ( j = p1+1; j <= n; j++ ) {
2862: if ( mm[j] ) {
2863: quo = dpm_sp_nf(psv,psiv,p1,j,&nf);
2864: if ( nf )
2865: error("dpm_schreyer_base : cannot happen");
2866: NEXTNODE(b0,b); BDY(b) = (pointer)quo;
1.13 noro 2867: }
2868: }
2869: }
2870: }
2871: get_eg(&eg1); add_eg(&egsp,&eg0,&eg1); print_eg("SP",&egsp);
2872: get_eg(&eg0);
2873: get_eg(&eg1); add_eg(&egnf,&eg0,&eg1); print_eg("NF",&egnf); printf("\n");
1.12 noro 2874: if ( b0 ) NEXT(b) = 0;
1.17 noro 2875: for ( t0 = 0, nd = BDY(g); nd; nd = NEXT(nd) ) {
1.12 noro 2876: dpm_ht((DPM)BDY(nd),&dpm); NEXTNODE(t0,t); BDY(t) = (pointer)dpm;
2877: }
2878: if ( t0 ) NEXT(t) = 0;
2879: MKLIST(l,t0);
2880: dmm_stack = push_schreyer_order(l,dmm_stack);
1.16 noro 2881: // b0 = dpm_sort_list(b0);
1.12 noro 2882: // get_eg(&eg0);
2883: // b0 = dpm_reduceall(b0);
2884: // get_eg(&eg1); add_eg(&egra,&eg0,&eg1); print_eg("RA",&egra);
1.4 noro 2885: MKLIST(*s,b0);
1.8 noro 2886: // print_eg("red",&egred); printf("\n");
1.15 noro 2887: printf("sch_count=%d, schlast_count=%d\n",sch_count,schlast_count);
2888: }
2889:
1.17 noro 2890: void dpm_list_to_array(LIST g,VECT *psvect,VECT *psindvect)
2891: {
2892: NODE nd,t;
2893: int n;
2894: VECT psv,psiv;
2895: DPM *ps;
2896: int i,max,pos;
2897: LIST *psi;
2898: LIST l;
2899: Z iz;
2900:
2901: nd = BDY(g);
2902: n = length(nd);
2903: MKVECT(psv,n+1);
2904: ps = (DPM *)BDY(psv);
2905: for ( i = 1, t = nd; i <= n; i++, t = NEXT(t) ) ps[i] = (DPM)BDY(t);
2906: for ( i = 1, max = 0; i <= n; i++ )
2907: if ( (pos=BDY(ps[i])->pos) > max ) max = pos;
2908: MKVECT(psiv,max+1);
2909: psi = (LIST *)BDY(psiv);
2910: for ( i = 0; i <= max; i++ ) {
2911: MKLIST(l,0); psi[i] = l;
2912: }
2913: for ( i = n; i >= 1; i-- ) {
2914: pos = BDY(ps[i])->pos;
2915: STOZ(i,iz); MKNODE(nd,iz,BDY(psi[pos])); BDY(psi[pos]) = nd;
2916: }
2917: *psvect = psv; *psindvect = psiv;
2918: }
2919:
1.19 noro 2920: #if 0
1.17 noro 2921: void dpm_insert_to_zlist(VECT psiv,int pos,int i)
2922: {
2923: LIST l;
2924: NODE prev,cur,nd;
2925: int curi;
2926: Z iz;
2927:
2928: l = (LIST)BDY(psiv)[pos];
2929: for ( prev = 0, cur = BDY(l); cur; cur = NEXT(cur) ) {
2930: curi = ZTOS((Q)BDY(cur));
2931: if ( curi == i )
2932: error("dpm_insert_to_list : invalid index");
2933: if ( i < curi ) break;
2934: prev = cur;
2935: }
2936: STOZ(i,iz); MKNODE(nd,iz,cur);
2937: if ( prev == 0 ) BDY(l) = nd;
2938: else NEXT(prev) = nd;
2939: }
1.19 noro 2940: #else
2941: void dpm_insert_to_zlist(VECT psiv,int pos,int i)
2942: {
2943: LIST l;
2944: NODE prev,cur,nd;
2945: int curi;
2946: Z iz;
2947:
2948: l = (LIST)BDY(psiv)[pos];
2949: for ( prev = 0, cur = BDY(l); cur; cur = NEXT(cur) ) {
2950: curi = ZTOS((Q)BDY(cur));
2951: if ( curi == i )
2952: error("dpm_insert_to_list : invalid index");
2953: prev = cur;
2954: }
2955: STOZ(i,iz); MKNODE(nd,iz,cur);
2956: if ( prev == 0 ) BDY(l) = nd;
2957: else NEXT(prev) = nd;
2958: }
2959: #endif
1.17 noro 2960:
2961: void dpm_schreyer_base_zlist(LIST g,LIST *s)
2962: {
2963: NODE nd,t0,t,b0,b;
2964: int n,i,j,k,nv,max,pos;
2965: LIST l;
2966: DP h,t1,t2;
2967: MP d;
2968: DMM r0,r,r1;
2969: DPM sp,nf,dpm;
2970: DPM *ps;
2971: VECT psv,psiv;
2972: DPM quo;
2973: DP *mm;
2974: LIST *psi;
2975: NODE n1,n2,n3;
2976: int p1,p2,p3;
2977: Z iz;
2978: struct oEGT eg0,eg1,egsp,egnf;
2979: extern struct oEGT egred;
2980: extern int sch_count,schrec_count,schlast_count;
2981:
2982: sch_count = schlast_count= 0;
2983: init_eg(&egra);
2984: init_eg(&egsp);
2985: init_eg(&egnf);
2986: dpm_list_to_array(g,&psv,&psiv);
2987: ps = (DPM *)BDY(psv);
2988: psi = (LIST *)BDY(psiv);
2989: nv = ps[1]->nv;
2990: n = psv->len-1;
2991: max = psiv->len-1;
2992: mm = (DP *)MALLOC((n+1)*sizeof(DP));
2993: b0 = 0;
2994: get_eg(&eg0);
2995: for ( i = 1; i <= max; i++ ) {
2996: memset(mm,0,(n+1)*sizeof(DP));
2997: for ( n1 = BDY((LIST)psi[i]); n1; n1 = NEXT(n1) ) {
2998: p1 = ZTOS((Q)BDY(n1));
2999: for ( n2 = NEXT(n1); n2; n2 = NEXT(n2) ) {
3000: p2 = ZTOS((Q)BDY(n2));
3001: mm[p2] = dpm_sp_hm(ps[p1],ps[p2]);
3002: }
3003: for ( n2 = NEXT(n1); n2; n2 = NEXT(n2) ) {
3004: p2 = ZTOS((Q)BDY(n2));
3005: if ( !mm[p2] ) continue;
3006: for ( h = mm[p2], n3 = NEXT(n1); n3; n3 = NEXT(n3) ) {
3007: p3 = ZTOS((Q)BDY(n3));
3008: if ( n3 != n2 && mm[p3] && dp_redble(mm[p3],h) ) mm[p3] = 0;
3009: }
3010: }
3011: for ( j = p1+1; j <= n; j++ ) {
3012: if ( mm[j] ) {
1.19 noro 3013: quo = dpm_sp_nf_zlist(psv,psiv,p1,j,0,&nf);
1.17 noro 3014: if ( nf )
3015: error("dpm_schreyer_base : cannot happen");
3016: NEXTNODE(b0,b); BDY(b) = (pointer)quo;
3017: }
3018: }
3019: }
3020: }
3021: get_eg(&eg1); add_eg(&egsp,&eg0,&eg1); print_eg("SP",&egsp);
3022: get_eg(&eg0);
3023: get_eg(&eg1); add_eg(&egnf,&eg0,&eg1); print_eg("NF",&egnf); printf("\n");
3024: if ( b0 ) NEXT(b) = 0;
3025: for ( t0 = 0, nd = BDY(g); nd; nd = NEXT(nd) ) {
3026: dpm_ht((DPM)BDY(nd),&dpm); NEXTNODE(t0,t); BDY(t) = (pointer)dpm;
3027: }
3028: if ( t0 ) NEXT(t) = 0;
3029: MKLIST(l,t0);
3030: dmm_stack = push_schreyer_order(l,dmm_stack);
3031: // b0 = dpm_sort_list(b0);
3032: // get_eg(&eg0);
3033: // b0 = dpm_reduceall(b0);
3034: // get_eg(&eg1); add_eg(&egra,&eg0,&eg1); print_eg("RA",&egra);
3035: MKLIST(*s,b0);
3036: // print_eg("red",&egred); printf("\n");
3037: printf("sch_count=%d, schlast_count=%d\n",sch_count,schlast_count);
3038: }
1.15 noro 3039:
1.19 noro 3040: static int compdp_nv;
3041:
3042: int compdp_revgradlex(DP *a,DP *b)
3043: {
3044: return -cmpdl_revgradlex(compdp_nv,BDY(*a)->dl,BDY(*b)->dl);
3045: }
3046:
3047: int compdp_lex(DP *a,DP *b)
3048: {
3049: return -cmpdl_lex(compdp_nv,BDY(*a)->dl,BDY(*b)->dl);
3050: }
3051:
1.22 noro 3052: DMMstack_array dpm_schreyer_frame(NODE g,int lex)
1.15 noro 3053: {
3054: LIST l;
3055: NODE nd,in,b0,b,n1,n2,n3,t;
3056: NODE *psi;
3057: long p1,p2,p3;
1.19 noro 3058: int nv,n,i,j,k,max,pos,level;
1.15 noro 3059: DMMstack s,s1;
3060: DMM m1,m0,dmm;
3061: MP mp;
3062: DP dp,h;
1.19 noro 3063: DP *m,*m2;
1.15 noro 3064: DPM dpm,dpm0,dpm1;
3065: VECT psv,psiv;
3066: DPM *ps;
3067: DMMstack_array dmmstack_array;
3068:
3069: nd = g;
1.19 noro 3070: compdp_nv = nv = ((DPM)BDY(nd))->nv;
1.15 noro 3071: s = 0;
3072: level = 0;
3073: while ( 1 ) {
3074: /* store the current nd to a DMMstack */
3075: n = length(nd);
3076: NEWDMMstack(s1);
3077: MKLIST(l,nd); s1->obj = l;
3078: s1->rank = n;
3079: s1->in = (DMM *)MALLOC((n+1)*sizeof(DMM));
3080: s1->sum = (DMM *)MALLOC((n+1)*sizeof(DMM));
3081: NEXT(s1) = s;
3082: if ( s ) {
3083: for ( i = 1, in = nd; i <= n; i++, in = NEXT(in) ) {
3084: m1 = s1->in[i] = BDY((DPM)BDY(in));
3085: NEWMP(mp); mp->dl = m1->dl; mp->c = m1->c; NEXT(mp) = 0;
3086: MKDP(nv,mp,dp); dp->sugar = mp->dl->td;
3087: m0 = s->sum[m1->pos]; MKDPM(nv,m0,dpm0);
3088: mulobjdpm(CO,(Obj)dp,dpm0,&dpm1);
3089: s1->sum[i] = BDY(dpm1);
3090: }
3091: } else {
1.18 noro 3092: for ( i = 1, in = nd; i <= n; i++, in = NEXT(in) ) {
3093: m0 = BDY((DPM)BDY(in));
3094: NEWDMM(m1); *m1 = *m0; m1->c = (Obj)ONE; NEXT(m1) = 0;
3095: s1->sum[i] = s1->in[i] = m1;
3096: }
1.15 noro 3097: }
3098: s = s1;
3099: level++;
1.22 noro 3100: if ( DP_Print ) printf("level=%d,len=%d\n",level,n);
1.15 noro 3101:
3102: /* create new list */
3103: MKVECT(psv,n+1);
3104: ps = (DPM *)BDY(psv);
3105: for ( i = 1, t = nd; i <= n; i++, t = NEXT(t) ) ps[i] = (DPM)BDY(t);
3106: for ( i = 1, max = 0; i <= n; i++ )
3107: if ( (pos=BDY(ps[i])->pos) > max ) max = pos;
3108: MKVECT(psiv,max+1);
3109: psi = (NODE *)BDY(psiv);
3110: for ( i = n; i >= 1; i-- ) {
3111: pos = BDY(ps[i])->pos;
3112: MKNODE(nd,(long)i,psi[pos]); psi[pos] = nd;
3113: }
3114: m = (DP *)MALLOC((n+1)*sizeof(DP));
1.19 noro 3115: m2 = (DP *)MALLOC((n+1)*sizeof(DP));
1.15 noro 3116: b0 = 0;
3117: for ( i = 1; i <= max; i++ ) {
3118: for ( n1 = psi[i]; n1; n1 = NEXT(n1) ) {
3119: bzero(m,(n+1)*sizeof(DP));
3120: p1 = (long)BDY(n1);
3121: for ( n2 = NEXT(n1); n2; n2 = NEXT(n2) ) {
3122: p2 = (long)BDY(n2);
3123: h = dpm_sp_hm(ps[p1],ps[p2]);
3124: for ( n3 = NEXT(n1); n3 != n2; n3 = NEXT(n3) ) {
3125: p3 = (long)BDY(n3);
3126: if ( m[p3] ) {
3127: if ( dp_redble(h,m[p3]) ) {
3128: h = 0; break;
3129: }
3130: if ( dp_redble(m[p3],h) ) m[p3] = 0;
3131: }
3132: }
3133: if ( h ) m[p2] = h;
3134: }
1.22 noro 3135: if ( lex ) {
3136: // compress m to m2
3137: for ( j = 0, n2 = NEXT(n1); n2; n2 = NEXT(n2) ) {
3138: p2 = (long)BDY(n2);
3139: if ( m[p2] ) m2[j++] = m[p2];
3140: }
3141: qsort(m2,j,sizeof(DP),(int (*)(const void *,const void *))compdp_lex);
3142: for ( k = 0; k < j; k++ ) {
3143: NEWDMM(dmm); dmm->dl = BDY(m2[k])->dl; dmm->pos = p1; dmm->c = (Obj)ONE;
1.15 noro 3144: MKDPM(nv,dmm,dpm);
3145: NEXTNODE(b0,b); BDY(b) = (pointer)dpm;
3146: }
1.22 noro 3147: } else {
3148: for ( n2 = NEXT(n1); n2; n2 = NEXT(n2) ) {
3149: p2 = (long)BDY(n2);
3150: if ( m[p2] ) {
3151: NEWDMM(dmm); dmm->dl = BDY(m[p2])->dl; dmm->pos = p1; dmm->c = (Obj)ONE;
3152: MKDPM(nv,dmm,dpm);
3153: NEXTNODE(b0,b); BDY(b) = (pointer)dpm;
3154: }
3155: }
1.15 noro 3156: }
3157: }
3158: }
3159: if ( !b0 ) {
3160: NEWDMMstack_array(dmmstack_array);
3161: dmmstack_array->len = level;
3162: dmmstack_array->body = (DMMstack *)MALLOC(level*sizeof(DMMstack));
3163: for ( i = level-1, s1 = s; i >= 0; i--, s1 = NEXT(s1) )
3164: dmmstack_array->body[i] = s1;
3165: return dmmstack_array;
3166: } else {
3167: NEXT(b) = 0;
3168: nd = b0;
3169: }
3170: }
1.4 noro 3171: }
1.15 noro 3172:
3173: int sch_count,schlast_count;
1.4 noro 3174:
1.3 noro 3175: int compdmm_schreyer(int n,DMM m1,DMM m2)
3176: {
1.15 noro 3177: int pos1,pos2,t,npos1,npos2;
3178: DMM *in,*sum;
3179: DMMstack s;
3180: static DL d1=0,d2=0;
3181: static int dlen=0;
3182:
3183: sch_count++;
3184: pos1 = m1->pos; pos2 = m2->pos;
3185: if ( pos1 == pos2 ) return (*cmpdl)(n,m1->dl,m2->dl);
3186: if ( n > dlen ) {
3187: NEWDL(d1,n); NEWDL(d2,n); dlen = n;
3188: }
3189: sum = dmm_stack->sum;
3190: _adddl(n,m1->dl,sum[pos1]->dl,d1);
3191: _adddl(n,m2->dl,sum[pos2]->dl,d2);
3192: t = (*cmpdl)(n,d1,d2);
3193: if ( sum[pos1]->pos == sum[pos2]->pos && t == 0 ) {
3194: for ( s = dmm_stack; s; s = NEXT(s) ) {
3195: in = s->in;
3196: npos1 = in[pos1]->pos;
3197: npos2 = in[pos2]->pos;
3198: if ( npos1 == npos2 ) break;
3199: else {
3200: pos1 = npos1;
3201: pos2 = npos2;
3202: }
3203: }
3204: // (pos1,pos2) = the last pair s.t. pos1 != pos2
3205: // simply compare pos1 and pos2
1.3 noro 3206: if ( pos1 < pos2 ) return 1;
3207: else if ( pos1 > pos2 ) return -1;
3208: } else {
1.15 noro 3209: schlast_count++;
3210: pos1 = sum[pos1]->pos;
3211: pos2 = sum[pos2]->pos;
3212: if ( dpm_base_ordtype == 1 ) {
3213: if ( pos1 < pos2 ) return 1;
3214: else if ( pos1 > pos2 ) return -1;
3215: else return t;
3216: } else {
3217: if ( t ) return t;
3218: else if ( pos1 < pos2 ) return 1;
3219: else if ( pos1 > pos2 ) return -1;
3220: else return 0;
3221: }
1.3 noro 3222: }
3223: }
3224:
1.20 noro 3225: int compdmm_schreyer_old(int n,DMM m1,DMM m2)
3226: {
3227: int pos1,pos2,t,npos1,npos2;
3228: DMM *in,*sum;
3229: DMMstack s;
3230: static DL d1=0,d2=0;
3231: static int dlen=0;
3232:
3233: sch_count++;
3234: pos1 = m1->pos; pos2 = m2->pos;
3235: if ( pos1 == pos2 ) return (*cmpdl)(n,m1->dl,m2->dl);
3236: if ( n > dlen ) {
3237: NEWDL(d1,n); NEWDL(d2,n); dlen = n;
3238: }
3239: sum = dmm_stack->sum;
3240: _copydl(n,m1->dl,d1);
3241: _copydl(n,m2->dl,d2);
3242: for ( s = dmm_stack; s; s = NEXT(s) ) {
3243: in = s->in;
3244: _addtodl(n,in[pos1]->dl,d1);
3245: _addtodl(n,in[pos2]->dl,d2);
3246: if ( in[pos1]->pos == in[pos2]->pos && _eqdl(n,d1,d2)) {
3247: if ( pos1 < pos2 ) return 1;
3248: else if ( pos1 > pos2 ) return -1;
3249: else return 0;
3250: }
3251: pos1 = in[pos1]->pos;
3252: pos2 = in[pos2]->pos;
3253: if ( pos1 == pos2 ) return (*cmpdl)(n,d1,d2);
3254: }
3255: if ( dpm_base_ordtype == 1 ) {
3256: if ( pos1 < pos2 ) return 1;
3257: else if ( pos1 > pos2 ) return -1;
3258: else return (*cmpdl)(n,d1,d2);
3259: } else {
3260: t = (*cmpdl)(n,d1,d2);
3261: if ( t ) return t;
3262: else if ( pos1 < pos2 ) return 1;
3263: else if ( pos1 > pos2 ) return -1;
3264: else return 0;
3265: }
3266: }
3267:
3268: extern int NaiveSchreyer;
3269:
1.4 noro 3270: int compdmm(int n,DMM m1,DMM m2)
3271: {
1.20 noro 3272: int t,t1;
1.19 noro 3273: int *base_ord;
1.4 noro 3274:
3275: switch ( dpm_ordtype ) {
1.9 noro 3276: case 0: /* TOP ord->pos */
1.4 noro 3277: t = (*cmpdl)(n,m1->dl,m2->dl);
3278: if ( t ) return t;
3279: else if ( m1->pos < m2->pos ) return 1;
3280: else if ( m1->pos > m2->pos ) return -1;
3281: else return 0;
1.9 noro 3282: case 1: /* POT : pos->ord */
1.4 noro 3283: if ( m1->pos < m2->pos ) return 1;
3284: else if ( m1->pos > m2->pos ) return -1;
3285: else return (*cmpdl)(n,m1->dl,m2->dl);
1.9 noro 3286: case 2: /* wPOT: weight->pos->ord */
3287: if ( m1->dl->td > m2->dl->td ) return 1;
3288: else if ( m1->dl->td < m2->dl->td ) return 1;
3289: else if ( m1->pos < m2->pos ) return 1;
3290: else if ( m1->pos > m2->pos ) return -1;
3291: else return (*cmpdl)(n,m1->dl,m2->dl);
3292: case 3: /* Schreyer */
1.20 noro 3293: if ( NaiveSchreyer )
3294: t = compdmm_schreyer_old(n,m1,m2);
3295: else
3296: t = compdmm_schreyer(n,m1,m2);
3297: return t;
1.19 noro 3298: case 4: /* POT with base_ord */
3299: base_ord = dp_current_spec->module_base_ord;
3300: if ( base_ord[m1->pos] < base_ord[m2->pos] ) return 1;
3301: else if ( base_ord[m1->pos] > base_ord[m2->pos] ) return -1;
3302: else return (*cmpdl)(n,m1->dl,m2->dl);
1.4 noro 3303: default:
3304: error("compdmm : invalid dpm_ordtype");
3305: }
3306: }
1.1 noro 3307:
3308: void adddpm(VL vl,DPM p1,DPM p2,DPM *pr)
3309: {
1.4 noro 3310: int n,s;
1.1 noro 3311: DMM m1,m2,mr=0,mr0;
3312: Obj t;
3313: DL d;
3314:
3315: if ( !p1 )
3316: *pr = p2;
3317: else if ( !p2 )
3318: *pr = p1;
3319: else {
1.4 noro 3320: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) {
3321: s = compdmm(n,m1,m2);
3322: switch ( s ) {
1.1 noro 3323: case 0:
3324: arf_add(vl,C(m1),C(m2),&t);
3325: if ( t ) {
3326: NEXTDMM(mr0,mr); mr->pos = m1->pos; mr->dl = m1->dl; C(mr) = t;
3327: }
3328: m1 = NEXT(m1); m2 = NEXT(m2); break;
3329: case 1:
3330: NEXTDMM(mr0,mr); mr->pos = m1->pos; mr->dl = m1->dl; C(mr) = C(m1);
3331: m1 = NEXT(m1); break;
3332: case -1:
3333: NEXTDMM(mr0,mr); mr->pos = m2->pos; mr->dl = m2->dl; C(mr) = C(m2);
3334: m2 = NEXT(m2); break;
3335: }
1.4 noro 3336: }
1.1 noro 3337: if ( !mr0 )
3338: if ( m1 )
3339: mr0 = m1;
3340: else if ( m2 )
3341: mr0 = m2;
3342: else {
3343: *pr = 0;
3344: return;
3345: }
3346: else if ( m1 )
3347: NEXT(mr) = m1;
3348: else if ( m2 )
3349: NEXT(mr) = m2;
3350: else
3351: NEXT(mr) = 0;
3352: MKDPM(NV(p1),mr0,*pr);
3353: if ( *pr )
3354: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
3355: }
3356: }
3357:
3358: void subdpm(VL vl,DPM p1,DPM p2,DPM *pr)
3359: {
3360: DPM t;
3361:
3362: if ( !p2 )
3363: *pr = p1;
3364: else {
3365: chsgndpm(p2,&t); adddpm(vl,p1,t,pr);
3366: }
3367: }
3368:
3369: void chsgndpm(DPM p,DPM *pr)
3370: {
3371: DMM m,mr=0,mr0;
3372: Obj r;
3373:
3374: if ( !p )
3375: *pr = 0;
3376: else {
3377: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
3378: NEXTDMM(mr0,mr); arf_chsgn(C(m),&C(mr)); mr->pos = m->pos; mr->dl = m->dl;
3379: }
3380: NEXT(mr) = 0; MKDPM(NV(p),mr0,*pr);
3381: if ( *pr )
3382: (*pr)->sugar = p->sugar;
3383: }
3384: }
3385:
1.8 noro 3386: void mulcmp(Obj c,MP m)
3387: {
3388: MP t;
3389: Obj c1;
3390:
3391: for ( t = m; t; t = NEXT(t) ) {
3392: arf_mul(CO,c,C(t),&c1); C(t) = c1;
3393: }
3394: }
3395:
3396: void mulcdmm(Obj c,DMM m)
3397: {
3398: DMM t;
3399: Obj c1;
3400:
3401: for ( t = m; t; t = NEXT(t) ) {
3402: arf_mul(CO,c,C(t),&c1); C(t) = c1;
3403: }
3404: }
3405:
1.1 noro 3406: void mulcdpm(VL vl,Obj c,DPM p,DPM *pr)
3407: {
3408: DMM m,mr=0,mr0;
3409:
3410: if ( !p || !c )
3411: *pr = 0;
3412: else if ( NUM(c) && UNIQ((Q)c) )
3413: *pr = p;
3414: else if ( NUM(c) && MUNIQ((Q)c) )
3415: chsgndpm(p,pr);
3416: else {
3417: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
3418: NEXTDMM(mr0,mr);
3419: arf_mul(vl,C(m),c,&C(mr));
3420: mr->pos = m->pos;
3421: mr->dl = m->dl;
3422: }
3423: NEXT(mr) = 0; MKDPM(NV(p),mr0,*pr);
3424: if ( *pr )
3425: (*pr)->sugar = p->sugar;
3426: }
3427: }
3428:
3429: void comm_mulmpdpm(VL vl,MP m0,DPM p,DPM *pr)
3430: {
3431: DMM m,mr=0,mr0;
3432: DL d;
3433: Obj c;
3434: int n;
3435:
3436: if ( !p )
3437: *pr = 0;
3438: else {
3439: n = NV(p);
3440: d = m0->dl;
3441: c = C(m0);
3442: for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) {
3443: NEXTDMM(mr0,mr);
3444: arf_mul(vl,C(m),c,&C(mr));
3445: mr->pos = m->pos;
3446: adddl(n,m->dl,d,&mr->dl);
3447: }
3448: NEXT(mr) = 0; MKDPM(NV(p),mr0,*pr);
3449: if ( *pr )
3450: (*pr)->sugar = p->sugar;
3451: }
3452: }
3453:
3454: void weyl_mulmpdpm(VL vl,MP m0,DPM p,DPM *pr)
3455: {
3456: DPM r,t,t1;
3457: DMM m;
3458: DL d0;
3459: int n,n2,l,i,j,tlen;
3460: struct oMP mp;
3461: static DMM *w,*psum;
3462: static struct cdl *tab;
3463: static int wlen;
3464: static int rtlen;
3465:
3466: if ( !p )
3467: *pr = 0;
3468: else {
3469: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
3470: if ( l > wlen ) {
3471: if ( w ) GCFREE(w);
3472: w = (DMM *)MALLOC(l*sizeof(DMM));
3473: wlen = l;
3474: }
3475: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
3476: w[i] = m;
3477:
3478: n = NV(p); n2 = n>>1;
3479: d0 = m0->dl;
3480: for ( i = 0, tlen = 1; i < n2; i++ )
3481: tlen *= d0->d[n2+i]+1;
3482: if ( tlen > rtlen ) {
3483: if ( tab ) GCFREE(tab);
3484: if ( psum ) GCFREE(psum);
3485: rtlen = tlen;
3486: tab = (struct cdl *)MALLOC(rtlen*sizeof(struct cdl));
3487: psum = (DMM *)MALLOC(rtlen*sizeof(DMM));
3488: }
3489: bzero(psum,tlen*sizeof(DMM));
3490: for ( i = l-1; i >= 0; i-- ) {
3491: bzero(tab,tlen*sizeof(struct cdl));
3492: mp.dl = w[i]->dl; mp.c = C(w[i]); mp.next = 0;
3493: weyl_mulmm(vl,m0,&mp,n,tab,tlen);
3494: for ( j = 0; j < tlen; j++ ) {
3495: if ( tab[j].c ) {
3496: NEWDMM(m); m->dl = tab[j].d; m->pos = w[i]->pos; C(m) = (Obj)tab[j].c; NEXT(m) = psum[j];
3497: psum[j] = m;
3498: }
3499: }
3500: }
3501: for ( j = tlen-1, r = 0; j >= 0; j-- )
3502: if ( psum[j] ) {
3503: MKDPM(n,psum[j],t); adddpm(vl,r,t,&t1); r = t1;
3504: }
3505: if ( r )
3506: r->sugar = p->sugar + m0->dl->td;
3507: *pr = r;
3508: }
3509: }
3510:
3511: void mulobjdpm(VL vl,Obj p1,DPM p2,DPM *pr)
3512: {
3513: MP m;
3514: DPM s,t,u;
3515:
3516: if ( !p1 || !p2 )
3517: *pr = 0;
3518: else if ( OID(p1) != O_DP )
3519: mulcdpm(vl,p1,p2,pr);
3520: else {
3521: s = 0;
3522: for ( m = BDY((DP)p1); m; m = NEXT(m) ) {
3523: if ( do_weyl )
3524: weyl_mulmpdpm(vl,m,p2,&t);
3525: else
3526: comm_mulmpdpm(vl,m,p2,&t);
3527: adddpm(vl,s,t,&u); s = u;
3528: }
3529: *pr = s;
3530: }
3531: }
3532:
3533: int compdpm(VL vl,DPM p1,DPM p2)
3534: {
3535: int n,t;
3536: DMM m1,m2;
3537:
3538: if ( !p1 )
3539: return p2 ? -1 : 0;
3540: else if ( !p2 )
3541: return 1;
3542: else if ( NV(p1) != NV(p2) ) {
3543: error("compdpm : size mismatch");
3544: return 0; /* XXX */
3545: } else {
3546: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2);
3547: m1 && m2; m1 = NEXT(m1), m2 = NEXT(m2) )
3548: if ( (t = compdmm(n,m1,m2)) ||
3549: (t = arf_comp(vl,C(m1),C(m2)) ) )
3550: return t;
3551: if ( m1 )
3552: return 1;
3553: else if ( m2 )
3554: return -1;
3555: else
3556: return 0;
3557: }
3558: }
3559:
1.11 noro 3560: void dpm_removecont2(DPM p1,DPM p2,DPM *r1p,DPM *r2p,Z *contp);
3561:
1.5 noro 3562: // p = ...+c*<<0,...0:pos>>+...
3563: DPM dpm_eliminate_term(DPM a,DPM p,Obj c,int pos)
3564: {
3565: MP d0,d;
1.10 noro 3566: DL dl;
1.5 noro 3567: DMM m;
3568: DP f;
1.11 noro 3569: DPM a1,p1,r,r1,dmy;
3570: Z dmyz;
1.5 noro 3571:
3572: if ( !a ) return 0;
3573: d0 = 0;
3574: for ( m = BDY(a); m; m = NEXT(m) )
3575: if ( m->pos == pos ) {
1.10 noro 3576: NEXTMP(d0,d);
3577: arf_chsgn(m->c,&d->c);
1.19 noro 3578: if ( !dp_current_spec || !dp_current_spec->module_top_weight )
1.11 noro 3579: d->dl = m->dl;
1.10 noro 3580: else {
3581: NEWDL(dl,NV(a));
3582: _copydl(NV(a),m->dl,dl);
3583: dl->td -= dp_current_spec->module_top_weight[pos-1];
3584: d->dl = dl;
3585: }
1.5 noro 3586: }
3587: if ( d0 ) {
1.6 noro 3588: NEXT(d) = 0; MKDP(NV(a),d0,f);
3589: mulcdpm(CO,c,a,&a1);
3590: mulobjdpm(CO,(Obj)f,p,&p1);
3591: adddpm(CO,a1,p1,&r);
1.11 noro 3592: dpm_removecont2(0,r,&dmy,&r1,&dmyz);
3593: return r1;
1.5 noro 3594: } else
1.6 noro 3595: return a;
1.5 noro 3596: }
3597:
3598: // <<...:i>> -> <<...:tab[i]>>
3599: DPM dpm_compress(DPM p,int *tab)
3600: {
3601: DMM m,mr0,mr;
3602: DPM t;
3603:
3604: if ( !p ) return 0;
3605: else {
3606: for ( m = BDY(p), mr0 = 0; m; m = NEXT(m) ) {
3607: NEXTDMM(mr0,mr);
1.6 noro 3608: mr->dl = m->dl; mr->c = m->c; mr->pos = tab[m->pos];
1.5 noro 3609: if ( mr->pos <= 0 )
3610: error("dpm_compress : invalid position");
3611: }
3612: NEXT(mr) = 0;
3613: MKDPM(p->nv,mr0,t); t->sugar = p->sugar;
3614: return t;
3615: }
3616: }
3617:
3618: // input : s, s = syz(m) output simplified s, m
1.21 noro 3619: // assuming the term order is POT
1.10 noro 3620: void dpm_simplify_syz(LIST s,LIST m,LIST *s1,LIST *m1,LIST *w1)
1.5 noro 3621: {
1.9 noro 3622: int lm,ls,i,j,k,pos,nv;
1.5 noro 3623: DPM *am,*as;
3624: DPM p;
3625: DMM d;
3626: Obj c;
1.10 noro 3627: Z q;
3628: int *tab,*dd,*new_w;
1.5 noro 3629: NODE t,t1;
3630:
3631: lm = length(BDY(m));
3632: am = (DPM *)MALLOC((lm+1)*sizeof(DPM));
3633: ls = length(BDY(s));
3634: as = (DPM *)MALLOC(ls*sizeof(DPM));
3635: for ( i = 1, t = BDY(m); i <= lm; i++, t = NEXT(t) ) am[i] = (DPM)BDY(t);
1.6 noro 3636: for ( i = 0, t = BDY(s); i < ls; i++, t = NEXT(t) ) as[i] = (DPM)BDY(t);
1.5 noro 3637:
3638: for ( i = 0; i < ls; i++ ) {
3639: p = as[i];
1.6 noro 3640: if ( p == 0 ) continue;
1.9 noro 3641: nv = NV(p);
1.21 noro 3642: for ( d = BDY(p); d; ) {
1.9 noro 3643: dd = d->dl->d;
3644: for ( k = 0; k < nv; k++ ) if ( dd[k] ) break;
3645: if ( k == nv ) break;
1.21 noro 3646: pos = d->pos;
3647: while ( d && d->pos == pos ) d = NEXT(d);
1.9 noro 3648: }
1.5 noro 3649: if ( d ) {
3650: c = d->c; pos = d->pos;
3651: for ( j = 0; j < ls; j++ )
1.6 noro 3652: if ( j != i ) {
1.5 noro 3653: as[j] = dpm_eliminate_term(as[j],p,c,pos);
1.6 noro 3654: }
1.5 noro 3655: // remove m[i]
3656: am[pos] = 0;
3657: as[i] = 0;
3658: }
3659: }
3660: // compress s
3661: // create index table from am[]
3662: // (0 0 * 0 * ...) -> (0 0 1 0 2 ... ) which means 2->1, 4->2, ...
3663: tab = (int *)MALLOC((lm+1)*sizeof(int));
3664: for ( j = 0, i = 1; i <= lm; i++ ) {
3665: if ( am[i] ) { j++; tab[i] = j; }
1.10 noro 3666: else { tab[i] = 0; }
1.5 noro 3667: }
3668: t = 0;
3669: for ( i = ls-1; i >= 0; i-- )
3670: if ( as[i] ) {
3671: p = dpm_compress(as[i],tab);
3672: MKNODE(t1,(pointer)p,t); t = t1;
3673: }
3674: MKLIST(*s1,t);
1.10 noro 3675:
1.19 noro 3676: if ( dp_current_spec && dp_current_spec->module_top_weight ) {
1.10 noro 3677: new_w = (int *)MALLOC(j*sizeof(int));
3678: for ( j = 0, i = 1; i <= lm; i++ )
3679: if ( tab[i] ) { new_w[j++] = dp_current_spec->module_top_weight[i-1]; }
3680: t = 0;
3681: for ( i = j-1; i >= 0; i-- ) {
3682: STOZ(new_w[i],q);
3683: MKNODE(t1,q,t); t = t1;
3684: }
3685: } else
3686: t = 0;
3687: MKLIST(*w1,t);
3688:
1.5 noro 3689: t = 0;
3690: for ( i = lm; i >= 1; i-- )
3691: if ( am[i] ) {
3692: MKNODE(t1,(pointer)am[i],t); t = t1;
3693: }
3694: MKLIST(*m1,t);
3695: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[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 / asir2018 / engine