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