![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to _distm.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine|
Return to _distm.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.2 noro 1: /*
2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
3: * All rights reserved.
4: *
5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
6: * non-exclusive and royalty-free license to use, copy, modify and
7: * redistribute, solely for non-commercial and non-profit purposes, the
8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
9: * conditions of this Agreement. For the avoidance of doubt, you acquire
10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
11: * third party developer retains all rights, including but not limited to
12: * copyrights, in and to the SOFTWARE.
13: *
14: * (1) FLL does not grant you a license in any way for commercial
15: * purposes. You may use the SOFTWARE only for non-commercial and
16: * non-profit purposes only, such as academic, research and internal
17: * business use.
18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
19: * international copyright treaties. If you make copies of the SOFTWARE,
20: * with or without modification, as permitted hereunder, you shall affix
21: * to all such copies of the SOFTWARE the above copyright notice.
22: * (3) An explicit reference to this SOFTWARE and its copyright owner
23: * shall be made on your publication or presentation in any form of the
24: * results obtained by use of the SOFTWARE.
25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.3 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2 noro 27: * for such modification or the source code of the modified part of the
28: * SOFTWARE.
29: *
30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
47: *
1.4 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/engine/_distm.c,v 1.3 2000/08/22 05:04:05 noro Exp $
1.2 noro 49: */
1.1 noro 50: #include "ca.h"
51: #include "inline.h"
52:
53: #define NV(p) ((p)->nv)
54: #define C(p) ((p)->c)
55: #if 0
56: #define ITOS(p) (((unsigned int)(p))>>1)
57: #define STOI(i) ((P)((((unsigned int)(i))<<1)|1))
58: #else
59: #define ITOS(p) (((unsigned int)(p))&0x7fffffff)
60: #define STOI(i) ((P)((unsigned int)(i)|0x80000000))
61: #endif
62:
1.4 ! noro 63: struct cdlm {
! 64: int c;
! 65: DL d;
! 66: };
! 67:
1.1 noro 68: extern int (*cmpdl)();
69: extern int do_weyl;
70:
71: void _dptodp();
72: void adddl_dup();
1.4 ! noro 73: void adddl_destructive();
1.1 noro 74: void _mulmdm_dup();
75: void _free_dp();
76: void _comm_mulmd_dup();
77: void _weyl_mulmd_dup();
78: void _weyl_mulmmm_dup();
79: void _weyl_mulmdm_dup();
1.4 ! noro 80: void _comm_mulmd_tab();
! 81: void _comm_mulmd_tab_destructive();
1.1 noro 82:
83: MP _mp_free_list;
84: DP _dp_free_list;
85: DL _dl_free_list;
86: int current_dl_length;
87:
88: #define _NEWDL_NOINIT(d,n) if ((n)!= current_dl_length){_dl_free_list=0; current_dl_length=(n);} if(!_dl_free_list)_DL_alloc(); (d)=_dl_free_list; _dl_free_list = *((DL *)_dl_free_list)
89: #define _NEWDL(d,n) if ((n)!= current_dl_length){_dl_free_list=0; current_dl_length=(n);} if(!_dl_free_list)_DL_alloc(); (d)=_dl_free_list; _dl_free_list = *((DL *)_dl_free_list); bzero((d),(((n)+1)*sizeof(int)))
90: #define _NEWMP(m) if(!_mp_free_list)_MP_alloc(); (m)=_mp_free_list; _mp_free_list = NEXT(_mp_free_list)
91: #define _MKDP(n,m,d) if(!_dp_free_list)_DP_alloc(); (d)=_dp_free_list; _dp_free_list = (DP)BDY(_dp_free_list); (d)->nv=(n); BDY(d)=(m)
92:
93: #define _NEXTMP(r,c) \
94: if(!(r)){_NEWMP(r);(c)=(r);}else{_NEWMP(NEXT(c));(c)=NEXT(c);}
95:
96: #define _NEXTMP2(r,c,s) \
97: if(!(r)){(c)=(r)=(s);}else{NEXT(c)=(s);(c)=(s);}
98:
99: #define FREEDL(m) *((DL *)m)=_dl_free_list; _dl_free_list=(m)
100: #define FREEMP(m) NEXT(m)=_mp_free_list; _mp_free_list=(m)
101: #define FREEDP(m) BDY(m)=(MP)_dp_free_list; _dp_free_list=(m)
102:
103:
104: void _DL_alloc()
105: {
106: int *p;
107: int i,dl_len;
108: static int DL_alloc_count;
109:
110: /* fprintf(stderr,"DL_alloc : %d \n",++DL_alloc_count); */
111: dl_len = (current_dl_length+1);
112: p = (int *)GC_malloc(128*dl_len*sizeof(int));
113: for ( i = 0; i < 128; i++, p += dl_len ) {
114: *(DL *)p = _dl_free_list;
115: _dl_free_list = (DL)p;
116: }
117: }
118:
119: void _MP_alloc()
120: {
121: MP p;
122: int i;
123: static int MP_alloc_count;
124:
125: /* fprintf(stderr,"MP_alloc : %d \n",++MP_alloc_count); */
126: p = (MP)GC_malloc(1024*sizeof(struct oMP));
127: for ( i = 0; i < 1024; i++ ) {
128: p[i].next = _mp_free_list; _mp_free_list = &p[i];
129: }
130: }
131:
132: void _DP_alloc()
133: {
134: DP p;
135: int i;
136: static int DP_alloc_count;
137:
138: /* fprintf(stderr,"DP_alloc : %d \n",++DP_alloc_count); */
139: p = (DP)GC_malloc(1024*sizeof(struct oDP));
140: for ( i = 0; i < 1024; i++ ) {
141: p[i].body = (MP)_dp_free_list; _dp_free_list = &p[i];
142: }
143: }
144:
145: /* merge p1 and p2 into pr */
146:
147: void _addmd_destructive(mod,p1,p2,pr)
148: int mod;
149: DP p1,p2,*pr;
150: {
151: int n;
152: MP m1,m2,mr,mr0,s;
153: int t;
154:
155: if ( !p1 )
156: *pr = p2;
157: else if ( !p2 )
158: *pr = p1;
159: else {
160: for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; )
161: switch ( (*cmpdl)(n,m1->dl,m2->dl) ) {
162: case 0:
163: t = (ITOS(C(m1))+ITOS(C(m2))) - mod;
164: if ( t < 0 )
165: t += mod;
166: s = m1; m1 = NEXT(m1);
167: if ( t ) {
168: _NEXTMP2(mr0,mr,s); C(mr) = STOI(t);
169: } else {
170: FREEDL(s->dl); FREEMP(s);
171: }
172: s = m2; m2 = NEXT(m2); FREEDL(s->dl); FREEMP(s);
173: break;
174: case 1:
175: s = m1; m1 = NEXT(m1); _NEXTMP2(mr0,mr,s);
176: break;
177: case -1:
178: s = m2; m2 = NEXT(m2); _NEXTMP2(mr0,mr,s);
179: break;
180: }
181: if ( !mr0 )
182: if ( m1 )
183: mr0 = m1;
184: else if ( m2 )
185: mr0 = m2;
186: else {
187: *pr = 0;
188: return;
189: }
190: else if ( m1 )
191: NEXT(mr) = m1;
192: else if ( m2 )
193: NEXT(mr) = m2;
194: else
195: NEXT(mr) = 0;
196: _MKDP(NV(p1),mr0,*pr);
197: if ( *pr )
198: (*pr)->sugar = MAX(p1->sugar,p2->sugar);
199: FREEDP(p1); FREEDP(p2);
200: }
201: }
202:
203: void _mulmd_dup(mod,p1,p2,pr)
204: int mod;
205: DP p1,p2,*pr;
206: {
207: if ( !do_weyl )
208: _comm_mulmd_dup(mod,p1,p2,pr);
209: else
210: _weyl_mulmd_dup(mod,p1,p2,pr);
211: }
212:
213: void _comm_mulmd_dup(mod,p1,p2,pr)
214: int mod;
215: DP p1,p2,*pr;
216: {
217: MP m;
218: DP s,t,u;
219: int i,l,l1;
220: static MP *w;
221: static int wlen;
222:
223: if ( !p1 || !p2 )
224: *pr = 0;
225: else {
226: for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ );
227: for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ );
228: if ( l1 < l ) {
229: t = p1; p1 = p2; p2 = t;
230: l = l1;
231: }
232: if ( l > wlen ) {
233: if ( w ) GC_free(w);
234: w = (MP *)MALLOC(l*sizeof(MP));
235: wlen = l;
236: }
237: for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ )
238: w[i] = m;
239: for ( s = 0, i = l-1; i >= 0; i-- ) {
240: _mulmdm_dup(mod,p1,w[i],&t); _addmd_destructive(mod,s,t,&u); s = u;
241: }
242: bzero(w,l*sizeof(MP));
243: *pr = s;
244: }
245: }
246:
247: void _weyl_mulmd_dup(mod,p1,p2,pr)
248: int mod;
249: DP p1,p2,*pr;
250: {
251: MP m;
252: DP s,t,u;
253: int i,l,l1;
254: static MP *w;
255: static int wlen;
256:
257: if ( !p1 || !p2 )
258: *pr = 0;
259: else {
1.4 ! noro 260: for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ );
1.1 noro 261: if ( l > wlen ) {
262: if ( w ) GC_free(w);
263: w = (MP *)MALLOC(l*sizeof(MP));
264: wlen = l;
265: }
1.4 ! noro 266: for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ )
1.1 noro 267: w[i] = m;
268: for ( s = 0, i = l-1; i >= 0; i-- ) {
1.4 ! noro 269: _weyl_mulmdm_dup(mod,w[i],p2,&t); _addmd_destructive(mod,s,t,&u); s = u;
1.1 noro 270: }
271: bzero(w,l*sizeof(MP));
272: *pr = s;
273: }
274: }
275:
276: void _mulmdm_dup(mod,p,m0,pr)
277: int mod;
278: DP p;
279: MP m0;
280: DP *pr;
281: {
282: MP m,mr,mr0;
283: DL d,dt,dm;
284: int c,n,r,i;
285: int *pt,*p1,*p2;
286:
287: if ( !p )
288: *pr = 0;
289: else {
290: for ( mr0 = 0, m = BDY(p), c = ITOS(C(m0)), d = m0->dl, n = NV(p);
291: m; m = NEXT(m) ) {
292: _NEXTMP(mr0,mr);
293: C(mr) = STOI(dmar(ITOS(C(m)),c,0,mod));
294: _NEWDL_NOINIT(dt,n); mr->dl = dt;
295: dm = m->dl;
296: dt->td = d->td + dm->td;
297: for ( i = 0, pt = dt->d, p1=d->d, p2 = dm->d; i < n; i++ )
298: *pt++ = *p1++ + *p2++;
299: }
300: NEXT(mr) = 0; _MKDP(NV(p),mr0,*pr);
301: if ( *pr )
302: (*pr)->sugar = p->sugar + m0->dl->td;
303: }
304: }
305:
1.4 ! noro 306: void _weyl_mulmmm_dup_bug();
! 307:
! 308: void _weyl_mulmdm_dup(mod,m0,p,pr)
1.1 noro 309: int mod;
1.4 ! noro 310: MP m0;
1.1 noro 311: DP p;
312: DP *pr;
313: {
314: DP r,t,t1;
315: MP m;
1.4 ! noro 316: DL d0;
! 317: int n,n2,l,i,j,tlen;
! 318: static MP *w,*psum;
! 319: static struct cdlm *tab;
1.1 noro 320: static int wlen;
1.4 ! noro 321: static int rtlen;
1.1 noro 322:
323: if ( !p )
324: *pr = 0;
325: else {
326: for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ );
327: if ( l > wlen ) {
328: if ( w ) GC_free(w);
329: w = (MP *)MALLOC(l*sizeof(MP));
330: wlen = l;
331: }
332: for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ )
333: w[i] = m;
1.4 ! noro 334: n = NV(p); n2 = n>>1;
! 335: d0 = m0->dl;
! 336:
! 337: for ( i = 0, tlen = 1; i < n2; i++ )
! 338: tlen *= d0->d[n2+i]+1;
! 339: if ( tlen > rtlen ) {
! 340: if ( tab ) GC_free(tab);
! 341: if ( psum ) GC_free(psum);
! 342: rtlen = tlen;
! 343: tab = (struct cdlm *)MALLOC(rtlen*sizeof(struct cdlm));
! 344: psum = (MP *)MALLOC(rtlen*sizeof(MP));
1.1 noro 345: }
1.4 ! noro 346: bzero(psum,tlen*sizeof(MP));
! 347: for ( i = l-1; i >= 0; i-- ) {
! 348: bzero(tab,tlen*sizeof(struct cdlm));
! 349: _weyl_mulmmm_dup(mod,m0,w[i],n,tab,tlen);
! 350: for ( j = 0; j < tlen; j++ ) {
! 351: if ( tab[j].c ) {
! 352: _NEWMP(m); m->dl = tab[j].d;
! 353: C(m) = STOI(tab[j].c); NEXT(m) = psum[j];
! 354: psum[j] = m;
! 355: }
! 356: }
! 357: }
! 358: for ( j = tlen-1, r = 0; j >= 0; j-- )
! 359: if ( psum[j] ) {
! 360: _MKDP(n,psum[j],t); _addmd_destructive(mod,r,t,&t1); r = t1;
! 361: }
1.1 noro 362: if ( r )
363: r->sugar = p->sugar + m0->dl->td;
364: *pr = r;
365: }
366: }
1.4 ! noro 367:
1.1 noro 368: /* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */
369:
1.4 ! noro 370: void _weyl_mulmmm_dup(mod,m0,m1,n,rtab,rtablen)
1.1 noro 371: int mod;
372: MP m0,m1;
373: int n;
1.4 ! noro 374: struct cdlm *rtab;
! 375: int rtablen;
1.1 noro 376: {
377: MP m,mr,mr0;
378: DP r,t,t1;
379: int c,c0,c1,cc;
1.4 ! noro 380: DL d,d0,d1,dt;
! 381: int i,j,a,b,k,l,n2,s,min,h,curlen;
! 382: struct cdlm *p;
! 383: static int *ctab;
! 384: static struct cdlm *tab;
1.1 noro 385: static int tablen;
1.4 ! noro 386: static struct cdlm *tmptab;
! 387: static int tmptablen;
1.1 noro 388:
1.4 ! noro 389: if ( !m0 || !m1 ) {
! 390: rtab[0].c = 0;
! 391: rtab[0].d = 0;
! 392: return;
! 393: }
! 394: c0 = ITOS(C(m0)); c1 = ITOS(C(m1));
! 395: c = dmar(c0,c1,0,mod);
! 396: d0 = m0->dl; d1 = m1->dl;
! 397: n2 = n>>1;
! 398: curlen = 1;
! 399:
! 400: _NEWDL(d,n);
! 401: if ( n & 1 )
! 402: /* offset of h-degree */
! 403: d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1];
! 404: else
! 405: d->td = 0;
! 406: rtab[0].c = c;
! 407: rtab[0].d = d;
! 408:
! 409: if ( rtablen > tmptablen ) {
! 410: if ( tmptab ) GC_free(tmptab);
! 411: tmptab = (struct cdlm *)MALLOC(rtablen*sizeof(struct cdlm));
! 412: tmptablen = rtablen;
! 413: }
1.1 noro 414:
1.4 ! noro 415: for ( i = 0; i < n2; i++ ) {
! 416: a = d0->d[i]; b = d1->d[n2+i];
! 417: k = d0->d[n2+i]; l = d1->d[i];
! 418: if ( !k || !l ) {
! 419: a += l;
! 420: b += k;
! 421: s = a+b;
! 422: for ( j = 0, p = rtab; j < curlen; j++, p++ ) {
! 423: if ( p->c ) {
! 424: dt = p->d;
! 425: dt->d[i] = a;
! 426: dt->d[n2+i] = b;
! 427: dt->td += s;
! 428: }
! 429: }
! 430: curlen *= k+1;
! 431: continue;
! 432: }
! 433: if ( k+1 > tablen ) {
! 434: if ( tab ) GC_free(tab);
! 435: if ( ctab ) GC_free(ctab);
! 436: tablen = k+1;
! 437: tab = (struct cdlm *)MALLOC(tablen*sizeof(struct cdlm));
! 438: ctab = (int *)MALLOC(tablen*sizeof(int));
! 439: }
! 440: /* degree of xi^a*(Di^k*xi^l)*Di^b */
! 441: s = a+k+l+b;
! 442: /* compute xi^a*(Di^k*xi^l)*Di^b */
! 443: min = MIN(k,l);
! 444: mkwcm(k,l,mod,ctab);
! 445: bzero(tab,(k+1)*sizeof(struct cdlm));
! 446: /* n&1 != 0 => homogenized computation; dx-xd=h^2 */
1.1 noro 447: if ( n & 1 )
1.4 ! noro 448: for ( j = 0; j <= min; j++ ) {
! 449: _NEWDL(d,n);
! 450: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
! 451: d->td = s;
! 452: d->d[n-1] = s-(d->d[i]+d->d[n2+i]);
! 453: tab[j].d = d;
! 454: tab[j].c = ctab[j];
! 455: }
1.1 noro 456: else
1.4 ! noro 457: for ( j = 0; j <= min; j++ ) {
! 458: _NEWDL(d,n);
! 459: d->d[i] = l-j+a; d->d[n2+i] = k-j+b;
! 460: d->td = d->d[i]+d->d[n2+i]; /* XXX */
! 461: tab[j].d = d;
! 462: tab[j].c = ctab[j];
! 463: }
! 464: #if 0
! 465: _comm_mulmd_tab(mod,n,rtab,curlen,tab,k+1,tmptab);
! 466: for ( j = 0; j < curlen; j++ )
! 467: if ( rtab[j].d ) { FREEDL(rtab[j].d); }
! 468: for ( j = 0; j <= min; j++ )
! 469: if ( tab[j].d ) { FREEDL(tab[j].d); }
! 470: curlen *= k+1;
! 471: bcopy(tmptab,rtab,curlen*sizeof(struct cdlm));
! 472: #else
! 473: _comm_mulmd_tab_destructive(mod,n,rtab,curlen,tab,k+1);
! 474: for ( j = 0; j <= min; j++ )
! 475: if ( tab[j].d ) { FREEDL(tab[j].d); }
! 476: curlen *= k+1;
! 477: #endif
! 478: }
! 479: }
! 480:
! 481: /* direct product of two cdlm tables
! 482: rt[] = [
! 483: t[0]*t1[0],...,t[n-1]*t1[0],
! 484: t[0]*t1[1],...,t[n-1]*t1[1],
! 485: ...
! 486: t[0]*t1[n1-1],...,t[n-1]*t1[n1-1]
! 487: ]
! 488: */
! 489:
! 490: void _comm_mulmd_tab(mod,nv,t,n,t1,n1,rt)
! 491: int mod;
! 492: int nv;
! 493: struct cdlm *t;
! 494: int n;
! 495: struct cdlm *t1;
! 496: int n1;
! 497: struct cdlm *rt;
! 498: {
! 499: int i,j;
! 500: struct cdlm *p;
! 501: int c;
! 502: DL d;
! 503:
! 504: bzero(rt,n*n1*sizeof(struct cdlm));
! 505: for ( j = 0, p = rt; j < n1; j++ ) {
! 506: c = t1[j].c;
! 507: d = t1[j].d;
! 508: if ( !c )
! 509: break;
! 510: for ( i = 0; i < n; i++, p++ ) {
! 511: if ( t[i].c ) {
! 512: p->c = dmar(t[i].c,c,0,mod);
! 513: adddl_dup(nv,t[i].d,d,&p->d);
! 514: }
! 515: }
! 516: }
! 517: }
! 518:
! 519: void _comm_mulmd_tab_destructive(mod,nv,t,n,t1,n1)
! 520: int mod;
! 521: int nv;
! 522: struct cdlm *t;
! 523: int n;
! 524: struct cdlm *t1;
! 525: int n1;
! 526: {
! 527: int i,j;
! 528: struct cdlm *p;
! 529: int c;
! 530: DL d;
! 531:
! 532: for ( j = 1, p = t+n; j < n1; j++ ) {
! 533: c = t1[j].c;
! 534: d = t1[j].d;
! 535: if ( !c )
! 536: break;
! 537: for ( i = 0; i < n; i++, p++ ) {
! 538: if ( t[i].c ) {
! 539: p->c = dmar(t[i].c,c,0,mod);
! 540: adddl_dup(nv,t[i].d,d,&p->d);
1.1 noro 541: }
542: }
543: }
1.4 ! noro 544: c = t1[0].c;
! 545: d = t1[0].d;
! 546: for ( i = 0, p = t; i < n; i++, p++ )
! 547: if ( t[i].c ) {
! 548: p->c = dmar(t[i].c,c,0,mod);
! 549: /* t[i].d += d */
! 550: adddl_destructive(nv,t[i].d,d);
! 551: }
1.1 noro 552: }
553:
554: void _dp_red_mod_destructive(p1,p2,mod,rp)
555: DP p1,p2;
556: int mod;
557: DP *rp;
558: {
559: int i,n;
560: DL d1,d2,d;
561: MP m;
562: DP t,s;
563: int c,c1;
564:
565: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
566: _NEWDL(d,n); d->td = d1->td - d2->td;
567: for ( i = 0; i < n; i++ )
568: d->d[i] = d1->d[i]-d2->d[i];
569: c = invm(ITOS(BDY(p2)->c),mod); c1 = dmar(c,ITOS(BDY(p1)->c),0,mod);
570: _NEWMP(m); m->dl = d; m->c = STOI(mod-c1); NEXT(m) = 0;
571: _MKDP(n,m,s); s->sugar = d->td;
572: _mulmd_dup(mod,s,p2,&t); _free_dp(s);
573: _addmd_destructive(mod,p1,t,rp);
574: }
575:
576: void _dp_sp_mod_dup(p1,p2,mod,rp)
577: DP p1,p2;
578: int mod;
579: DP *rp;
580: {
581: int i,n,td;
582: int *w;
583: DL d1,d2,d;
584: MP m;
585: DP t,s,u;
586:
587: n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl;
588: w = (int *)ALLOCA(n*sizeof(int));
589: for ( i = 0, td = 0; i < n; i++ ) {
590: w[i] = MAX(d1->d[i],d2->d[i]); td += w[i];
591: }
592: _NEWDL(d,n); d->td = td - d1->td;
593: for ( i = 0; i < n; i++ )
594: d->d[i] = w[i] - d1->d[i];
595: _NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0;
596: _MKDP(n,m,s); s->sugar = d->td; _mulmd_dup(mod,s,p1,&t); _free_dp(s);
597: _NEWDL(d,n); d->td = td - d2->td;
598: for ( i = 0; i < n; i++ )
599: d->d[i] = w[i] - d2->d[i];
600: _NEWMP(m); m->dl = d; m->c = STOI(mod - ITOS(BDY(p1)->c)); NEXT(m) = 0;
601: _MKDP(n,m,s); s->sugar = d->td; _mulmd_dup(mod,s,p2,&u); _free_dp(s);
602: _addmd_destructive(mod,t,u,rp);
603: }
604:
605: void _dltodl(d,dr)
606: DL d;
607: DL *dr;
608: {
609: int i,n;
610: DL t;
611:
612: n = current_dl_length;
613: NEWDL(t,n); *dr = t;
614: t->td = d->td;
615: for ( i = 0; i < n; i++ )
616: t->d[i] = d->d[i];
617: }
618:
619: void adddl_dup(n,d1,d2,dr)
620: int n;
621: DL d1,d2;
622: DL *dr;
623: {
624: DL dt;
625: int i;
626:
1.4 ! noro 627: _NEWDL(dt,n);
! 628: *dr = dt;
1.1 noro 629: dt->td = d1->td + d2->td;
630: for ( i = 0; i < n; i++ )
631: dt->d[i] = d1->d[i]+d2->d[i];
1.4 ! noro 632: }
! 633:
! 634: /* d1 += d2 */
! 635:
! 636: void adddl_destructive(n,d1,d2)
! 637: int n;
! 638: DL d1,d2;
! 639: {
! 640: DL dt;
! 641: int i;
! 642:
! 643: d1->td += d2->td;
! 644: for ( i = 0; i < n; i++ )
! 645: d1->d[i] += d2->d[i];
! 646: }
! 647:
! 648: void _free_dlarray(a,n)
! 649: DL *a;
! 650: int n;
! 651: {
! 652: int i;
! 653:
! 654: for ( i = 0; i < n; i++ ) { FREEDL(a[i]); }
1.1 noro 655: }
656:
657: void _free_dp(f)
658: DP f;
659: {
660: MP m,s;
661:
662: if ( !f )
663: return;
664: m = f->body;
665: while ( m ) {
666: s = m; m = NEXT(m); FREEDL(s->dl); FREEMP(s);
667: }
668: FREEDP(f);
669: }
670:
671: void _dptodp(p,r)
672: DP p;
673: DP *r;
674: {
675: MP m,mr0,mr;
676:
677: if ( !p )
678: *r = 0;
679: else {
680: for ( m = BDY(p), mr0 = 0; m; m = NEXT(m) ) {
681: NEXTMP(mr0,mr);
682: _dltodl(m->dl,&mr->dl);
683: mr->c = m->c;
684: }
685: NEXT(mr) = 0;
686: MKDP(p->nv,mr0,*r);
687: (*r)->sugar = p->sugar;
688: }
689: }
690:
691: void _dp_nf_mod_destructive(b,g,ps,mod,full,rp)
692: NODE b;
693: DP g;
694: DP *ps;
695: int mod,full;
696: DP *rp;
697: {
698: DP u,p,d,s,t;
699: NODE l;
700: MP m,mr,mrd;
701: int sugar,psugar,n,h_reducible,i;
702:
703: if ( !g ) {
704: *rp = 0; return;
705: }
706: sugar = g->sugar;
707: n = g->nv;
708: for ( d = 0; g; ) {
709: for ( h_reducible = 0, l = b; l; l = NEXT(l) ) {
710: if ( dp_redble(g,p = ps[(int)BDY(l)]) ) {
711: h_reducible = 1;
712: psugar = (BDY(g)->dl->td - BDY(p)->dl->td) + p->sugar;
713: _dp_red_mod_destructive(g,p,mod,&u); g = u;
714: sugar = MAX(sugar,psugar);
715: if ( !g ) {
716: if ( d )
717: d->sugar = sugar;
718: _dptodp(d,rp); _free_dp(d); return;
719: }
720: break;
721: }
722: }
723: if ( !h_reducible ) {
724: /* head term is not reducible */
725: if ( !full ) {
726: if ( g )
727: g->sugar = sugar;
728: _dptodp(g,rp); _free_dp(g); return;
729: } else {
730: m = BDY(g);
731: if ( NEXT(m) ) {
732: BDY(g) = NEXT(m); NEXT(m) = 0;
733: } else {
734: FREEDP(g); g = 0;
735: }
736: if ( d ) {
737: for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) );
738: NEXT(mrd) = m;
739: } else {
740: _MKDP(n,m,d);
741: }
742: }
743: }
744: }
745: if ( d )
746: d->sugar = sugar;
747: _dptodp(d,rp); _free_dp(d);
748: }