=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/distm.c,v retrieving revision 1.1.1.1 retrieving revision 1.13 diff -u -p -r1.1.1.1 -r1.13 --- OpenXM_contrib2/asir2000/engine/distm.c 1999/12/03 07:39:08 1.1.1.1 +++ OpenXM_contrib2/asir2000/engine/distm.c 2003/07/21 01:57:51 1.13 @@ -1,24 +1,59 @@ -/* $OpenXM: OpenXM/src/asir99/engine/distm.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */ +/* + * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED + * All rights reserved. + * + * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited, + * non-exclusive and royalty-free license to use, copy, modify and + * redistribute, solely for non-commercial and non-profit purposes, the + * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and + * conditions of this Agreement. For the avoidance of doubt, you acquire + * only a limited right to use the SOFTWARE hereunder, and FLL or any + * third party developer retains all rights, including but not limited to + * copyrights, in and to the SOFTWARE. + * + * (1) FLL does not grant you a license in any way for commercial + * purposes. You may use the SOFTWARE only for non-commercial and + * non-profit purposes only, such as academic, research and internal + * business use. + * (2) The SOFTWARE is protected by the Copyright Law of Japan and + * international copyright treaties. If you make copies of the SOFTWARE, + * with or without modification, as permitted hereunder, you shall affix + * to all such copies of the SOFTWARE the above copyright notice. + * (3) An explicit reference to this SOFTWARE and its copyright owner + * shall be made on your publication or presentation in any form of the + * results obtained by use of the SOFTWARE. + * (4) In the event that you modify the SOFTWARE, you shall notify FLL by + * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification + * for such modification or the source code of the modified part of the + * SOFTWARE. + * + * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL + * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND + * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS + * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES' + * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY + * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY. + * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT, + * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY + * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL + * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES + * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES + * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY + * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF + * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART + * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY + * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, + * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. + * + * $OpenXM: OpenXM_contrib2/asir2000/engine/distm.c,v 1.12 2003/07/20 08:55:23 noro Exp $ +*/ #include "ca.h" #include "inline.h" -#define NV(p) ((p)->nv) -#define C(p) ((p)->c) -#if 0 -#define ITOS(p) (((unsigned int)(p))>>1) -#define STOI(i) ((P)((((unsigned int)(i))<<1)|1)) -#else -#define ITOS(p) (((unsigned int)(p))&0x7fffffff) -#define STOI(i) ((P)((unsigned int)(i)|0x80000000)) -#endif - extern int (*cmpdl)(); +extern int do_weyl; -void ptomd(vl,mod,dvl,p,pr) -VL vl,dvl; -int mod; -P p; -DP *pr; +void ptomd(VL vl,int mod,VL dvl,P p,DP *pr) { P t; @@ -26,11 +61,7 @@ DP *pr; mptomd(vl,mod,dvl,t,pr); } -void mptomd(vl,mod,dvl,p,pr) -VL vl,dvl; -int mod; -P p; -DP *pr; +void mptomd(VL vl,int mod,VL dvl,P p,DP *pr) { int n,i; VL tvl; @@ -60,9 +91,10 @@ DP *pr; } else { for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) { mptomd(vl,mod,dvl,COEF(dc),&t); - NEWDL(d,n); d->td = QTOS(DEG(dc)); d->d[i] = d->td; + NEWDL(d,n); d->d[i] = QTOS(DEG(dc)); + d->td = MUL_WEIGHT(d->d[i],i); NEWMP(m); m->dl = d; C(m) = (P)ONEM; NEXT(m) = 0; MKDP(n,m,u); - mulmd(vl,mod,t,u,&r); addmd(vl,mod,r,s,&t); s = t; + comm_mulmd(vl,mod,t,u,&r); addmd(vl,mod,r,s,&t); s = t; } *pr = s; } @@ -70,11 +102,7 @@ DP *pr; } } -void mdtop(vl,mod,dvl,p,pr) -VL vl,dvl; -int mod; -DP p; -P *pr; +void mdtop(VL vl,int mod,VL dvl,DP p,P *pr) { int n,i; DL d; @@ -98,10 +126,7 @@ P *pr; } } -void addmd(vl,mod,p1,p2,pr) -VL vl; -int mod; -DP p1,p2,*pr; +void addmd(VL vl,int mod,DP p1,DP p2,DP *pr) { int n; MP m1,m2,mr,mr0; @@ -159,10 +184,7 @@ DP p1,p2,*pr; } } -void submd(vl,mod,p1,p2,pr) -VL vl; -int mod; -DP p1,p2,*pr; +void submd(VL vl,int mod,DP p1,DP p2,DP *pr) { DP t; @@ -173,9 +195,7 @@ DP p1,p2,*pr; } } -void chsgnmd(mod,p,pr) -int mod; -DP p,*pr; +void chsgnmd(int mod,DP p,DP *pr) { MP m,mr,mr0; @@ -191,13 +211,21 @@ DP p,*pr; } } -void mulmd(vl,mod,p1,p2,pr) -VL vl; -int mod; -DP p1,p2,*pr; +void mulmd(VL vl,int mod,DP p1,DP p2,DP *pr) { + if ( !do_weyl ) + comm_mulmd(vl,mod,p1,p2,pr); + else + weyl_mulmd(vl,mod,p1,p2,pr); +} + +void comm_mulmd(VL vl,int mod,DP p1,DP p2,DP *pr) +{ MP m; DP s,t,u; + int i,l,l1; + static MP *w; + static int wlen; if ( !p1 || !p2 ) *pr = 0; @@ -206,20 +234,60 @@ DP p1,p2,*pr; else if ( OID(p2) <= O_P ) mulmdc(vl,mod,p1,(P)p2,pr); else { - for ( m = BDY(p2), s = 0; m; m = NEXT(m) ) { - mulmdm(vl,mod,p1,m,&t); addmd(vl,mod,s,t,&u); s = u; + for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ ); + for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ ); + if ( l1 < l ) { + t = p1; p1 = p2; p2 = t; + l = l1; } + if ( l > wlen ) { + if ( w ) GC_free(w); + w = (MP *)MALLOC(l*sizeof(MP)); + wlen = l; + } + for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ ) + w[i] = m; + for ( s = 0, i = l-1; i >= 0; i-- ) { + mulmdm(vl,mod,p1,w[i],&t); addmd(vl,mod,s,t,&u); s = u; + } + bzero(w,l*sizeof(MP)); *pr = s; } } -void mulmdm(vl,mod,p,m0,pr) -VL vl; -int mod; -DP p; -MP m0; -DP *pr; +void weyl_mulmd(VL vl,int mod,DP p1,DP p2,DP *pr) { + MP m; + DP s,t,u; + int i,l; + static MP *w; + static int wlen; + + if ( !p1 || !p2 ) + *pr = 0; + else if ( OID(p1) <= O_P ) + mulmdc(vl,mod,p2,(P)p1,pr); + else if ( OID(p2) <= O_P ) + mulmdc(vl,mod,p1,(P)p2,pr); + else { + for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ ); + if ( l > wlen ) { + if ( w ) GC_free(w); + w = (MP *)MALLOC(l*sizeof(MP)); + wlen = l; + } + for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ ) + w[i] = m; + for ( s = 0, i = l-1; i >= 0; i-- ) { + weyl_mulmdm(vl,mod,p1,w[i],&t); addmd(vl,mod,s,t,&u); s = u; + } + bzero(w,l*sizeof(MP)); + *pr = s; + } +} + +void mulmdm(VL vl,int mod,DP p,MP m0,DP *pr) +{ MP m,mr,mr0; P c; DL d; @@ -245,13 +313,113 @@ DP *pr; } } -void mulmdc(vl,mod,p,c,pr) -VL vl; -int mod; -DP p; -P c; -DP *pr; +void weyl_mulmdm(VL vl,int mod,DP p,MP m0,DP *pr) { + DP r,t,t1; + MP m; + int n,l,i; + static MP *w; + static int wlen; + + if ( !p ) + *pr = 0; + else { + for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ ); + if ( l > wlen ) { + if ( w ) GC_free(w); + w = (MP *)MALLOC(l*sizeof(MP)); + wlen = l; + } + for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ ) + w[i] = m; + for ( r = 0, i = l-1, n = NV(p); i >= 0; i-- ) { + weyl_mulmmm(vl,mod,w[i],m0,n,&t); + addmd(vl,mod,r,t,&t1); r = t1; + } + bzero(w,l*sizeof(MP)); + if ( r ) + r->sugar = p->sugar + m0->dl->td; + *pr = r; + } +} + +/* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */ + +void weyl_mulmmm(VL vl,int mod,MP m0,MP m1,int n,DP *pr) +{ + MP mr,mr0; + MQ mq; + DP r,t,t1; + P c,c0,c1; + DL d,d0,d1; + int i,j,a,b,k,l,n2,s,min; + static int *tab; + static int tablen; + + if ( !m0 || !m1 ) + *pr = 0; + else { + c0 = C(m0); c1 = C(m1); + mulmp(vl,mod,c0,c1,&c); + d0 = m0->dl; d1 = m1->dl; + n2 = n>>1; + if ( n & 1 ) { + /* homogenized computation; dx-xd=h^2 */ + /* offset of h-degree */ + NEWDL(d,n); + d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1]; + NEWMP(mr); mr->c = (P)ONEM; mr->dl = d; NEXT(mr) = 0; + MKDP(n,mr,r); r->sugar = d->td; + } else + r = (DP)ONEM; + for ( i = 0; i < n2; i++ ) { + a = d0->d[i]; b = d1->d[n2+i]; + k = d0->d[n2+i]; l = d1->d[i]; + + /* degree of xi^a*(Di^k*xi^l)*Di^b */ + a += l; + b += k; + s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i); + + /* compute xi^a*(Di^k*xi^l)*Di^b */ + min = MIN(k,l); + + if ( min+1 > tablen ) { + if ( tab ) GC_free(tab); + tab = (int *)MALLOC((min+1)*sizeof(int)); + tablen = min+1; + } + mkwcm(k,l,mod,tab); + if ( n & 1 ) + for ( mr0 = 0, j = 0; j <= min; j++ ) { + NEXTMP(mr0,mr); NEWDL(d,n); + d->d[i] = a-j; d->d[n2+i] = b-j; + d->td = s; + d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i)); + STOMQ(tab[j],mq); mr->c = (P)mq; mr->dl = d; + } + else + for ( mr0 = 0, s = 0, j = 0; j <= min; j++ ) { + NEXTMP(mr0,mr); NEWDL(d,n); + d->d[i] = a-j; d->d[n2+i] = b-j; + d->td = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i); /* XXX */ + s = MAX(s,d->td); /* XXX */ + STOMQ(tab[j],mq); mr->c = (P)mq; mr->dl = d; + } + bzero(tab,(min+1)*sizeof(int)); + if ( mr0 ) + NEXT(mr) = 0; + MKDP(n,mr0,t); + if ( t ) + t->sugar = s; + comm_mulmd(vl,mod,r,t,&t1); r = t1; + } + mulmdc(vl,mod,r,c,pr); + } +} + +void mulmdc(VL vl,int mod,DP p,P c,DP *pr) +{ MP m,mr,mr0; int t; MQ q; @@ -274,12 +442,7 @@ DP *pr; } } -void divsmdc(vl,mod,p,c,pr) -VL vl; -int mod; -DP p; -P c; -DP *pr; +void divsmdc(VL vl,int mod,DP p,P c,DP *pr) { MP m,mr,mr0; @@ -297,45 +460,100 @@ DP *pr; } } -#define MKDPM(n,m,d) (NEWDP(d),(d)->nv=(n),BDY(d)=(m)) - -void _mdtop(vl,mod,dvl,p,pr) -VL vl,dvl; -int mod; -DP p; -P *pr; +void _dtop_mod(VL vl,VL dvl,DP p,P *pr) { int n,i; DL d; MP m; P r,s,t,u,w; Q q; - MQ tq; VL tvl; if ( !p ) *pr = 0; else { for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) { - STOMQ(ITOS(C(m)),tq); t = (P)tq; + i = ITOS(m->c); STOQ(i,q); t = (P)q; for ( i = 0, d = m->dl, tvl = dvl; i < n; tvl = NEXT(tvl), i++ ) { - MKV(tvl->v,r); STOQ(d->d[i],q); pwrmp(vl,mod,r,q,&u); - mulmp(vl,mod,t,u,&w); t = w; + MKV(tvl->v,r); STOQ(d->d[i],q); pwrp(vl,r,q,&u); + mulp(vl,t,u,&w); t = w; } - addmp(vl,mod,s,t,&u); s = u; + addp(vl,s,t,&u); s = u; } - mptop(s,pr); + *pr = s; } } -void _addmd(vl,mod,p1,p2,pr) -VL vl; -int mod; -DP p1,p2,*pr; +void _dp_mod(DP p,int mod,NODE subst,DP *rp) { + MP m,mr,mr0; + P t,s; + NODE tn; + + if ( !p ) + *rp = 0; + else { + for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) { + for ( tn = subst, s = m->c; tn; tn = NEXT(NEXT(tn)) ) { + substp(CO,s,BDY(tn),BDY(NEXT(tn)),&t); + s = t; + } + ptomp(mod,s,&t); + if ( t ) { + NEXTMP(mr0,mr); mr->c = STOI(CONT((MQ)t)); mr->dl = m->dl; + } + } + if ( mr0 ) { + NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar; + } else + *rp = 0; + } +} + +void _dp_monic(DP p,int mod,DP *rp) +{ + MP m,mr,mr0; + int c,c1; + + if ( !p ) + *rp = 0; + else { + c = invm(ITOS(BDY(p)->c),mod); + for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) { + c1 = dmar(ITOS(m->c),c,0,mod); + NEXTMP(mr0,mr); mr->c = STOI(c1); mr->dl = m->dl; + } + NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar; + } +} + +void _printdp(DP d) +{ + int n,i; + MP m; + DL dl; + + if ( !d ) { + printf("0"); return; + } + for ( n = d->nv, m = BDY(d); m; m = NEXT(m) ) { + printf("%d*<<",ITOS(m->c)); + for ( i = 0, dl = m->dl; i < n-1; i++ ) + printf("%d,",dl->d[i]); + printf("%d",dl->d[i]); + printf(">>"); + if ( NEXT(m) ) + printf("+"); + } +} + +/* merge p1 and p2 into pr */ + +void addmd_destructive(int mod,DP p1,DP p2,DP *pr) +{ int n; - MP m1,m2,mr,mr0; + MP m1,m2,mr,mr0,s; int t; if ( !p1 ) @@ -349,16 +567,18 @@ DP p1,p2,*pr; t = (ITOS(C(m1))+ITOS(C(m2))) - mod; if ( t < 0 ) t += mod; + s = m1; m1 = NEXT(m1); if ( t ) { - NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = STOI(t); + NEXTMP2(mr0,mr,s); C(mr) = STOI(t); } - m1 = NEXT(m1); m2 = NEXT(m2); break; + s = m2; m2 = NEXT(m2); + break; case 1: - NEXTMP(mr0,mr); mr->dl = m1->dl; C(mr) = C(m1); - m1 = NEXT(m1); break; + s = m1; m1 = NEXT(m1); NEXTMP2(mr0,mr,s); + break; case -1: - NEXTMP(mr0,mr); mr->dl = m2->dl; C(mr) = C(m2); - m2 = NEXT(m2); break; + s = m2; m2 = NEXT(m2); NEXTMP2(mr0,mr,s); + break; } if ( !mr0 ) if ( m1 ) @@ -375,72 +595,85 @@ DP p1,p2,*pr; NEXT(mr) = m2; else NEXT(mr) = 0; - MKDPM(NV(p1),mr0,*pr); + MKDP(NV(p1),mr0,*pr); if ( *pr ) (*pr)->sugar = MAX(p1->sugar,p2->sugar); } } -void _submd(vl,mod,p1,p2,pr) -VL vl; -int mod; -DP p1,p2,*pr; +void mulmd_dup(int mod,DP p1,DP p2,DP *pr) { - DP t; - - if ( !p2 ) - *pr = p1; - else { - _chsgnmd(mod,p2,&t); _addmd(vl,mod,p1,t,pr); - } + if ( !do_weyl ) + comm_mulmd_dup(mod,p1,p2,pr); + else + weyl_mulmd_dup(mod,p1,p2,pr); } -void _chsgnmd(mod,p,pr) -int mod; -DP p,*pr; +void comm_mulmd_dup(int mod,DP p1,DP p2,DP *pr) { - MP m,mr,mr0; + MP m; + DP s,t,u; + int i,l,l1; + static MP *w; + static int wlen; - if ( !p ) + if ( !p1 || !p2 ) *pr = 0; else { - for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) { - NEXTMP(mr0,mr); C(mr) = STOI(mod - ITOS(C(m))); mr->dl = m->dl; + for ( m = BDY(p1), l1 = 0; m; m = NEXT(m), l1++ ); + for ( m = BDY(p2), l = 0; m; m = NEXT(m), l++ ); + if ( l1 < l ) { + t = p1; p1 = p2; p2 = t; + l = l1; } - NEXT(mr) = 0; MKDPM(NV(p),mr0,*pr); - if ( *pr ) - (*pr)->sugar = p->sugar; + if ( l > wlen ) { + if ( w ) GC_free(w); + w = (MP *)MALLOC(l*sizeof(MP)); + wlen = l; + } + for ( m = BDY(p2), i = 0; i < l; m = NEXT(m), i++ ) + w[i] = m; + for ( s = 0, i = l-1; i >= 0; i-- ) { + mulmdm_dup(mod,p1,w[i],&t); addmd_destructive(mod,s,t,&u); s = u; + } + bzero(w,l*sizeof(MP)); + *pr = s; } } -void _mulmd(vl,mod,p1,p2,pr) -VL vl; -int mod; -DP p1,p2,*pr; +void weyl_mulmd_dup(int mod,DP p1,DP p2,DP *pr) { MP m; DP s,t,u; + int i,l; + static MP *w; + static int wlen; if ( !p1 || !p2 ) *pr = 0; else { - for ( m = BDY(p2), s = 0; m; m = NEXT(m) ) { - _mulmdm(vl,mod,p1,m,&t); _addmd(vl,mod,s,t,&u); s = u; + for ( m = BDY(p1), l = 0; m; m = NEXT(m), l++ ); + if ( l > wlen ) { + if ( w ) GC_free(w); + w = (MP *)MALLOC(l*sizeof(MP)); + wlen = l; } + for ( m = BDY(p1), i = 0; i < l; m = NEXT(m), i++ ) + w[i] = m; + for ( s = 0, i = l-1; i >= 0; i-- ) { + weyl_mulmdm_dup(mod,w[i],p2,&t); addmd_destructive(mod,s,t,&u); s = u; + } + bzero(w,l*sizeof(MP)); *pr = s; } } -void _mulmdm(vl,mod,p,m0,pr) -VL vl; -int mod; -DP p; -MP m0; -DP *pr; +void mulmdm_dup(int mod,DP p,MP m0,DP *pr) { MP m,mr,mr0; - DL d; - int c,n,r; + DL d,dt,dm; + int c,n,i; + int *pt,*p1,*p2; if ( !p ) *pr = 0; @@ -449,167 +682,1631 @@ DP *pr; m; m = NEXT(m) ) { NEXTMP(mr0,mr); C(mr) = STOI(dmar(ITOS(C(m)),c,0,mod)); - adddl(n,m->dl,d,&mr->dl); + NEWDL_NOINIT(dt,n); mr->dl = dt; + dm = m->dl; + dt->td = d->td + dm->td; + for ( i = 0, pt = dt->d, p1=d->d, p2 = dm->d; i < n; i++ ) + *pt++ = *p1++ + *p2++; } - NEXT(mr) = 0; MKDPM(NV(p),mr0,*pr); + NEXT(mr) = 0; MKDP(NV(p),mr0,*pr); if ( *pr ) (*pr)->sugar = p->sugar + m0->dl->td; } } -void _dtop_mod(vl,dvl,p,pr) -VL vl,dvl; -DP p; -P *pr; +void weyl_mulmdm_dup(int mod,MP m0,DP p,DP *pr) { - int n,i; - DL d; + DP r,t,t1; MP m; - P r,s,t,u,w; - Q q; - VL tvl; + DL d0; + int n,n2,l,i,j,tlen; + static MP *w,*psum; + static struct cdlm *tab; + static int wlen; + static int rtlen; if ( !p ) *pr = 0; else { - for ( n = p->nv, m = BDY(p), s = 0; m; m = NEXT(m) ) { - i = ITOS(m->c); STOQ(i,q); t = (P)q; - for ( i = 0, d = m->dl, tvl = dvl; - i < n; tvl = NEXT(tvl), i++ ) { - MKV(tvl->v,r); STOQ(d->d[i],q); pwrp(vl,r,q,&u); - mulp(vl,t,u,&w); t = w; + for ( m = BDY(p), l = 0; m; m = NEXT(m), l++ ); + if ( l > wlen ) { + if ( w ) GC_free(w); + w = (MP *)MALLOC(l*sizeof(MP)); + wlen = l; + } + for ( m = BDY(p), i = 0; i < l; m = NEXT(m), i++ ) + w[i] = m; + n = NV(p); n2 = n>>1; + d0 = m0->dl; + + for ( i = 0, tlen = 1; i < n2; i++ ) + tlen *= d0->d[n2+i]+1; + if ( tlen > rtlen ) { + if ( tab ) GC_free(tab); + if ( psum ) GC_free(psum); + rtlen = tlen; + tab = (struct cdlm *)MALLOC(rtlen*sizeof(struct cdlm)); + psum = (MP *)MALLOC(rtlen*sizeof(MP)); + } + bzero(psum,tlen*sizeof(MP)); + for ( i = l-1; i >= 0; i-- ) { + bzero(tab,tlen*sizeof(struct cdlm)); + weyl_mulmmm_dup(mod,m0,w[i],n,tab,tlen); + for ( j = 0; j < tlen; j++ ) { + if ( tab[j].c ) { + NEWMP(m); m->dl = tab[j].d; + C(m) = STOI(tab[j].c); NEXT(m) = psum[j]; + psum[j] = m; + } } - addp(vl,s,t,&u); s = u; } - *pr = s; + for ( j = tlen-1, r = 0; j >= 0; j-- ) + if ( psum[j] ) { + MKDP(n,psum[j],t); addmd_destructive(mod,r,t,&t1); r = t1; + } + if ( r ) + r->sugar = p->sugar + m0->dl->td; + *pr = r; } } -void _dp_red_mod(p1,p2,mod,rp) -DP p1,p2; -int mod; -DP *rp; +/* m0 = x0^d0*x1^d1*... * dx0^d(n/2)*dx1^d(n/2+1)*... */ + +void weyl_mulmmm_dup(int mod,MP m0,MP m1,int n,struct cdlm *rtab,int rtablen) { - int i,n; - DL d1,d2,d; - MP m; - DP t,s; - int c,c1; + int c,c0,c1; + DL d,d0,d1,dt; + int i,j,a,b,k,l,n2,s,min,curlen; + struct cdlm *p; + static int *ctab; + static struct cdlm *tab; + static int tablen; + static struct cdlm *tmptab; + static int tmptablen; - n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl; - NEWDL(d,n); d->td = d1->td - d2->td; + if ( !m0 || !m1 ) { + rtab[0].c = 0; + rtab[0].d = 0; + return; + } + c0 = ITOS(C(m0)); c1 = ITOS(C(m1)); + c = dmar(c0,c1,0,mod); + d0 = m0->dl; d1 = m1->dl; + n2 = n>>1; + curlen = 1; + + NEWDL(d,n); + if ( n & 1 ) + /* offset of h-degree */ + d->td = d->d[n-1] = d0->d[n-1]+d1->d[n-1]; + else + d->td = 0; + rtab[0].c = c; + rtab[0].d = d; + + if ( rtablen > tmptablen ) { + if ( tmptab ) GC_free(tmptab); + tmptab = (struct cdlm *)MALLOC(rtablen*sizeof(struct cdlm)); + tmptablen = rtablen; + } + + for ( i = 0; i < n2; i++ ) { + a = d0->d[i]; b = d1->d[n2+i]; + k = d0->d[n2+i]; l = d1->d[i]; + + /* degree of xi^a*(Di^k*xi^l)*Di^b */ + a += l; + b += k; + s = MUL_WEIGHT(a,i)+MUL_WEIGHT(b,n2+i); + + if ( !k || !l ) { + for ( j = 0, p = rtab; j < curlen; j++, p++ ) { + if ( p->c ) { + dt = p->d; + dt->d[i] = a; + dt->d[n2+i] = b; + dt->td += s; + } + } + curlen *= k+1; + continue; + } + if ( k+1 > tablen ) { + if ( tab ) GC_free(tab); + if ( ctab ) GC_free(ctab); + tablen = k+1; + tab = (struct cdlm *)MALLOC(tablen*sizeof(struct cdlm)); + ctab = (int *)MALLOC(tablen*sizeof(int)); + } + /* compute xi^a*(Di^k*xi^l)*Di^b */ + min = MIN(k,l); + mkwcm(k,l,mod,ctab); + bzero(tab,(k+1)*sizeof(struct cdlm)); + /* n&1 != 0 => homogenized computation; dx-xd=h^2 */ + if ( n & 1 ) + for ( j = 0; j <= min; j++ ) { + NEWDL(d,n); + d->d[i] = a-j; d->d[n2+i] = b-j; + d->td = s; + d->d[n-1] = s-(MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i)); + tab[j].d = d; + tab[j].c = ctab[j]; + } + else + for ( j = 0; j <= min; j++ ) { + NEWDL(d,n); + d->d[i] = a-j; d->d[n2+i] = b-j; + d->td = MUL_WEIGHT(a-j,i)+MUL_WEIGHT(b-j,n2+i); /* XXX */ + tab[j].d = d; + tab[j].c = ctab[j]; + } + comm_mulmd_tab_destructive(mod,n,rtab,curlen,tab,k+1); + curlen *= k+1; + } +} + +void comm_mulmd_tab_destructive(int mod,int nv,struct cdlm *t,int n,struct cdlm *t1,int n1) +{ + int i,j; + struct cdlm *p; + int c; + DL d; + + for ( j = 1, p = t+n; j < n1; j++ ) { + c = t1[j].c; + d = t1[j].d; + if ( !c ) + break; + for ( i = 0; i < n; i++, p++ ) { + if ( t[i].c ) { + p->c = dmar(t[i].c,c,0,mod); + adddl_dup(nv,t[i].d,d,&p->d); + } + } + } + c = t1[0].c; + d = t1[0].d; + for ( i = 0, p = t; i < n; i++, p++ ) + if ( t[i].c ) { + p->c = dmar(t[i].c,c,0,mod); + /* t[i].d += d */ + adddl_destructive(nv,t[i].d,d); + } +} + +void adddl_dup(int n,DL d1,DL d2,DL *dr) +{ + DL dt; + int i; + + NEWDL(dt,n); + *dr = dt; + dt->td = d1->td + d2->td; for ( i = 0; i < n; i++ ) - d->d[i] = d1->d[i]-d2->d[i]; - c = invm(ITOS(BDY(p2)->c),mod); c1 = dmar(c,ITOS(BDY(p1)->c),0,mod); - NEWMP(m); m->dl = d; m->c = STOI(mod-c1); NEXT(m) = 0; - MKDP(n,m,s); s->sugar = d->td; - _mulmd(CO,mod,p2,s,&t); _addmd(CO,mod,p1,t,rp); + dt->d[i] = d1->d[i]+d2->d[i]; } -void _dp_mod(p,mod,subst,rp) -DP p; -int mod; -NODE subst; -DP *rp; +/* ------------------------------------------------------ */ + +#if defined(__GNUC__) +#define INLINE inline +#elif defined(VISUAL) +#define INLINE __inline +#else +#define INLINE +#endif + +typedef struct oPGeoBucket { + int m; + struct oND *body[32]; +} *PGeoBucket; + +typedef struct oND { + struct oNM *body; + int nv; + int sugar; +} *ND; + +typedef struct oNM { + struct oNM *next; + int td; + int c; + unsigned int dl[1]; +} *NM; + +typedef struct oND_pairs { + struct oND_pairs *next; + int i1,i2; + int td,sugar; + unsigned int lcm[1]; +} *ND_pairs; + +static ND *nps; +int nd_mod,nd_nvar; +int is_rlex; +int nd_epw,nd_bpe,nd_wpd; +unsigned int nd_mask[32]; +unsigned int nd_mask0; +NM _nm_free_list; +ND _nd_free_list; +ND_pairs _ndp_free_list; + +extern int Top,Reverse; +int nd_psn,nd_pslen; + +void GC_gcollect(); +NODE append_one(NODE,int); + +#define HTD(d) ((d)->body->td) +#define HDL(d) ((d)->body->dl) +#define HC(d) ((d)->body->c) + +#define NEWND_pairs(m) if(!_ndp_free_list)_NDP_alloc(); (m)=_ndp_free_list; _ndp_free_list = NEXT(_ndp_free_list) +#define NEWNM(m) if(!_nm_free_list)_NM_alloc(); (m)=_nm_free_list; _nm_free_list = NEXT(_nm_free_list) +#define MKND(n,m,d) if(!_nd_free_list)_ND_alloc(); (d)=_nd_free_list; _nd_free_list = (ND)BDY(_nd_free_list); (d)->nv=(n); BDY(d)=(m) + +#define NEXTNM(r,c) \ +if(!(r)){NEWNM(r);(c)=(r);}else{NEWNM(NEXT(c));(c)=NEXT(c);} +#define NEXTNM2(r,c,s) \ +if(!(r)){(c)=(r)=(s);}else{NEXT(c)=(s);(c)=(s);} +#define FREENM(m) NEXT(m)=_nm_free_list; _nm_free_list=(m) +#define FREEND(m) BDY(m)=(NM)_nd_free_list; _nd_free_list=(m) + +#define NEXTND_pairs(r,c) \ +if(!(r)){NEWND_pairs(r);(c)=(r);}else{NEWND_pairs(NEXT(c));(c)=NEXT(c);} + +ND_pairs crit_B( ND_pairs d, int s ); +void nd_gr(LIST f,LIST v,int m,struct order_spec *ord,LIST *rp); +NODE nd_setup(NODE f); +int nd_newps(ND a); +ND_pairs nd_minp( ND_pairs d, ND_pairs *prest ); +NODE update_base(NODE nd,int ndp); +static ND_pairs equivalent_pairs( ND_pairs d1, ND_pairs *prest ); +int crit_2( int dp1, int dp2 ); +ND_pairs crit_F( ND_pairs d1 ); +ND_pairs crit_M( ND_pairs d1 ); +ND_pairs nd_newpairs( NODE g, int t ); +ND_pairs update_pairs( ND_pairs d, NODE /* of index */ g, int t); +NODE nd_gb(NODE f); +void nd_free_private_storage(); +void _NM_alloc(); +void _ND_alloc(); +int ndl_td(unsigned int *d); +ND nd_add(ND p1,ND p2); +ND nd_mul_nm(ND p,NM m0); +ND nd_mul_term(ND p,int td,unsigned int *d); +ND nd_sp(ND_pairs p); +ND nd_find_reducer(ND g); +ND nd_nf(ND g,int full); +void ndl_print(unsigned int *dl); +void nd_print(ND p); +void ndp_print(ND_pairs d); +int nd_length(ND p); +void nd_monic(ND p); +void nd_mul_c(ND p,int mul); + +void nd_free_private_storage() { - MP m,mr,mr0; - P t,s; - NODE tn; + _nd_free_list = 0; + _nm_free_list = 0; + GC_gcollect(); +} +void _NM_alloc() +{ + NM p; + int i; + + for ( i = 0; i < 10240; i++ ) { + p = (NM)GC_malloc(sizeof(struct oNM)+(nd_wpd-1)*sizeof(unsigned int)); + p->next = _nm_free_list; _nm_free_list = p; + } +} + +void _ND_alloc() +{ + ND p; + int i; + + for ( i = 0; i < 1024; i++ ) { + p = (ND)GC_malloc(sizeof(struct oND)); + p->body = (NM)_nd_free_list; _nd_free_list = p; + } +} + +void _NDP_alloc() +{ + ND_pairs p; + int i; + + for ( i = 0; i < 10240; i++ ) { + p = (ND_pairs)GC_malloc(sizeof(struct oND_pairs) + +(nd_wpd-1)*sizeof(unsigned int)); + p->next = _ndp_free_list; _ndp_free_list = p; + } +} + +INLINE nd_length(ND p) +{ + NM m; + int i; + if ( !p ) - *rp = 0; + return 0; else { - for ( mr0 = 0, m = BDY(p); m; m = NEXT(m) ) { - for ( tn = subst, s = m->c; tn; tn = NEXT(NEXT(tn)) ) { - substp(CO,s,BDY(tn),BDY(NEXT(tn)),&t); - s = t; + for ( i = 0, m = BDY(p); m; m = NEXT(m), i++ ); + return i; + } +} + +INLINE int ndl_reducible(unsigned int *d1,unsigned int *d2) +{ + unsigned int u1,u2; + int i,j; + + switch ( nd_bpe ) { + case 4: + for ( i = 0; i < nd_wpd; i++ ) { + u1 = d1[i]; u2 = d2[i]; + if ( (u1&0xf0000000) < (u2&0xf0000000) ) return 0; + if ( (u1&0xf000000) < (u2&0xf000000) ) return 0; + if ( (u1&0xf00000) < (u2&0xf00000) ) return 0; + if ( (u1&0xf0000) < (u2&0xf0000) ) return 0; + if ( (u1&0xf000) < (u2&0xf000) ) return 0; + if ( (u1&0xf00) < (u2&0xf00) ) return 0; + if ( (u1&0xf0) < (u2&0xf0) ) return 0; + if ( (u1&0xf) < (u2&0xf) ) return 0; } - ptomp(mod,s,&t); - if ( t ) { - NEXTMP(mr0,mr); mr->c = STOI(CONT((MQ)t)); mr->dl = m->dl; + return 1; + break; + case 6: + for ( i = 0; i < nd_wpd; i++ ) { + u1 = d1[i]; u2 = d2[i]; + if ( (u1&0x3f000000) < (u2&0x3f000000) ) return 0; + if ( (u1&0xfc0000) < (u2&0xfc0000) ) return 0; + if ( (u1&0x3f000) < (u2&0x3f000) ) return 0; + if ( (u1&0xfc0) < (u2&0xfc0) ) return 0; + if ( (u1&0x3f) < (u2&0x3f) ) return 0; } + return 1; + break; + case 8: + for ( i = 0; i < nd_wpd; i++ ) { + u1 = d1[i]; u2 = d2[i]; + if ( (u1&0xff000000) < (u2&0xff000000) ) return 0; + if ( (u1&0xff0000) < (u2&0xff0000) ) return 0; + if ( (u1&0xff00) < (u2&0xff00) ) return 0; + if ( (u1&0xff) < (u2&0xff) ) return 0; + } + return 1; + break; + case 16: + for ( i = 0; i < nd_wpd; i++ ) { + u1 = d1[i]; u2 = d2[i]; + if ( (u1&0xffff0000) < (u2&0xffff0000) ) return 0; + if ( (u1&0xffff) < (u2&0xffff) ) return 0; + } + return 1; + break; + case 32: + for ( i = 0; i < nd_wpd; i++ ) + if ( d1[i] < d2[i] ) return 0; + return 1; + break; + default: + for ( i = 0; i < nd_wpd; i++ ) { + u1 = d1[i]; u2 = d2[i]; + for ( j = 0; j < nd_epw; j++ ) + if ( (u1&nd_mask[j]) < (u2&nd_mask[j]) ) return 0; + } + return 1; + } +} + +void ndl_lcm(unsigned int *d1,unsigned *d2,unsigned int *d) +{ + unsigned int t1,t2,u,u1,u2; + int i,j; + + switch ( nd_bpe ) { + case 4: + for ( i = 0; i < nd_wpd; i++ ) { + u1 = d1[i]; u2 = d2[i]; + t1 = (u1&0xf0000000); t2 = (u2&0xf0000000); u = t1>t2?t1:t2; + t1 = (u1&0xf000000); t2 = (u2&0xf000000); u |= t1>t2?t1:t2; + t1 = (u1&0xf00000); t2 = (u2&0xf00000); u |= t1>t2?t1:t2; + t1 = (u1&0xf0000); t2 = (u2&0xf0000); u |= t1>t2?t1:t2; + t1 = (u1&0xf000); t2 = (u2&0xf000); u |= t1>t2?t1:t2; + t1 = (u1&0xf00); t2 = (u2&0xf00); u |= t1>t2?t1:t2; + t1 = (u1&0xf0); t2 = (u2&0xf0); u |= t1>t2?t1:t2; + t1 = (u1&0xf); t2 = (u2&0xf); u |= t1>t2?t1:t2; + d[i] = u; + } + break; + case 6: + for ( i = 0; i < nd_wpd; i++ ) { + u1 = d1[i]; u2 = d2[i]; + t1 = (u1&0x3f000000); t2 = (u2&0x3f000000); u = t1>t2?t1:t2; + t1 = (u1&0xfc0000); t2 = (u2&0xfc0000); u |= t1>t2?t1:t2; + t1 = (u1&0x3f000); t2 = (u2&0x3f000); u |= t1>t2?t1:t2; + t1 = (u1&0xfc0); t2 = (u2&0xfc0); u |= t1>t2?t1:t2; + t1 = (u1&0x3f); t2 = (u2&0x3f); u |= t1>t2?t1:t2; + d[i] = u; + } + break; + case 8: + for ( i = 0; i < nd_wpd; i++ ) { + u1 = d1[i]; u2 = d2[i]; + t1 = (u1&0xff000000); t2 = (u2&0xff000000); u = t1>t2?t1:t2; + t1 = (u1&0xff0000); t2 = (u2&0xff0000); u |= t1>t2?t1:t2; + t1 = (u1&0xff00); t2 = (u2&0xff00); u |= t1>t2?t1:t2; + t1 = (u1&0xff); t2 = (u2&0xff); u |= t1>t2?t1:t2; + d[i] = u; + } + break; + case 16: + for ( i = 0; i < nd_wpd; i++ ) { + u1 = d1[i]; u2 = d2[i]; + t1 = (u1&0xffff0000); t2 = (u2&0xffff0000); u = t1>t2?t1:t2; + t1 = (u1&0xffff); t2 = (u2&0xffff); u |= t1>t2?t1:t2; + d[i] = u; + } + break; + case 32: + for ( i = 0; i < nd_wpd; i++ ) { + u1 = d1[i]; u2 = d2[i]; + d[i] = u1>u2?u1:u2; + } + break; + default: + for ( i = 0; i < nd_wpd; i++ ) { + u1 = d1[i]; u2 = d2[i]; + for ( j = 0, u = 0; j < nd_epw; j++ ) { + t1 = (u1&nd_mask[j]); t2 = (u2&nd_mask[j]); u |= t1>t2?t1:t2; + } + d[i] = u; + } + break; + } +} + +int ndl_td(unsigned int *d) +{ + unsigned int t,u; + int i,j; + + for ( t = 0, i = 0; i < nd_wpd; i++ ) { + u = d[i]; + for ( j = 0; j < nd_epw; j++, u>>=nd_bpe ) + t += (u&nd_mask0); + } + return t; +} + +INLINE int ndl_compare(unsigned int *d1,unsigned int *d2) +{ + int i; + + for ( i = 0; i < nd_wpd; i++, d1++, d2++ ) + if ( *d1 > *d2 ) + return is_rlex ? -1 : 1; + else if ( *d1 < *d2 ) + return is_rlex ? 1 : -1; + return 0; +} + +INLINE int ndl_equal(unsigned int *d1,unsigned int *d2) +{ + int i; + + for ( i = 0; i < nd_wpd; i++ ) + if ( d1[i] != d2[i] ) + return 0; + return 1; +} + +INLINE void ndl_add(unsigned int *d1,unsigned int *d2,unsigned int *d) +{ + int i; + + for ( i = 0; i < nd_wpd; i++ ) + d[i] = d1[i]+d2[i]; +} + +void ndl_sub(unsigned int *d1,unsigned int *d2,unsigned int *d) +{ + int i; + + for ( i = 0; i < nd_wpd; i++ ) + d[i] = d1[i]-d2[i]; +} + +int ndl_disjoint(unsigned int *d1,unsigned int *d2) +{ + unsigned int t1,t2,u,u1,u2; + int i,j; + + switch ( nd_bpe ) { + case 4: + for ( i = 0; i < nd_wpd; i++ ) { + u1 = d1[i]; u2 = d2[i]; + t1 = u1&0xf0000000; t2 = u2&0xf0000000; if ( t1&&t2 ) return 0; + t1 = u1&0xf000000; t2 = u2&0xf000000; if ( t1&&t2 ) return 0; + t1 = u1&0xf00000; t2 = u2&0xf00000; if ( t1&&t2 ) return 0; + t1 = u1&0xf0000; t2 = u2&0xf0000; if ( t1&&t2 ) return 0; + t1 = u1&0xf000; t2 = u2&0xf000; if ( t1&&t2 ) return 0; + t1 = u1&0xf00; t2 = u2&0xf00; if ( t1&&t2 ) return 0; + t1 = u1&0xf0; t2 = u2&0xf0; if ( t1&&t2 ) return 0; + t1 = u1&0xf; t2 = u2&0xf; if ( t1&&t2 ) return 0; + } + return 1; + break; + case 6: + for ( i = 0; i < nd_wpd; i++ ) { + u1 = d1[i]; u2 = d2[i]; + t1 = u1&0x3f000000; t2 = u2&0x3f000000; if ( t1&&t2 ) return 0; + t1 = u1&0xfc0000; t2 = u2&0xfc0000; if ( t1&&t2 ) return 0; + t1 = u1&0x3f000; t2 = u2&0x3f000; if ( t1&&t2 ) return 0; + t1 = u1&0xfc0; t2 = u2&0xfc0; if ( t1&&t2 ) return 0; + t1 = u1&0x3f; t2 = u2&0x3f; if ( t1&&t2 ) return 0; + } + return 1; + break; + case 8: + for ( i = 0; i < nd_wpd; i++ ) { + u1 = d1[i]; u2 = d2[i]; + t1 = u1&0xff000000; t2 = u2&0xff000000; if ( t1&&t2 ) return 0; + t1 = u1&0xff0000; t2 = u2&0xff0000; if ( t1&&t2 ) return 0; + t1 = u1&0xff00; t2 = u2&0xff00; if ( t1&&t2 ) return 0; + t1 = u1&0xff; t2 = u2&0xff; if ( t1&&t2 ) return 0; + } + return 1; + break; + case 16: + for ( i = 0; i < nd_wpd; i++ ) { + u1 = d1[i]; u2 = d2[i]; + t1 = u1&0xffff0000; t2 = u2&0xffff0000; if ( t1&&t2 ) return 0; + t1 = u1&0xffff; t2 = u2&0xffff; if ( t1&&t2 ) return 0; + } + return 1; + break; + case 32: + for ( i = 0; i < nd_wpd; i++ ) + if ( d1[i] && d2[i] ) return 0; + return 1; + break; + default: + for ( i = 0; i < nd_wpd; i++ ) { + u1 = d1[i]; u2 = d2[i]; + for ( j = 0; j < nd_epw; j++ ) { + if ( (u1&nd_mask0) && (u2&nd_mask0) ) return 0; + u1 >>= nd_bpe; u2 >>= nd_bpe; + } + } + return 1; + break; + } +} + +ND nd_reduce(ND p1,ND p2) +{ + int c,c1,c2,t,td,td2,mul; + NM m2,prev,head,cur,new; + unsigned int *d; + + if ( !p1 ) + return 0; + else { + c2 = invm(HC(p2),nd_mod); + c1 = nd_mod-HC(p1); + DMAR(c1,c2,0,nd_mod,mul); + td = HTD(p1)-HTD(p2); + d = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int)); + ndl_sub(HDL(p1),HDL(p2),d); + prev = 0; head = cur = BDY(p1); + NEWNM(new); + for ( m2 = BDY(p2); m2; ) { + td2 = new->td = m2->td+td; + ndl_add(m2->dl,d,new->dl); + if ( !cur ) { + c1 = C(m2); + DMAR(c1,mul,0,nd_mod,c2); + C(new) = c2; + if ( !prev ) { + prev = new; + NEXT(prev) = 0; + head = prev; + } else { + NEXT(prev) = new; + NEXT(new) = 0; + prev = new; + } + m2 = NEXT(m2); + NEWNM(new); + continue; + } + if ( cur->td > td2 ) + c = 1; + else if ( cur->td < td2 ) + c = -1; + else + c = ndl_compare(cur->dl,new->dl); + switch ( c ) { + case 0: + c2 = C(m2); + c1 = C(cur); + DMAR(c2,mul,c1,nd_mod,t); + if ( t ) + C(cur) = t; + else if ( !prev ) { + head = NEXT(cur); + FREENM(cur); + cur = head; + } else { + NEXT(prev) = NEXT(cur); + FREENM(cur); + cur = NEXT(prev); + } + m2 = NEXT(m2); + break; + case 1: + prev = cur; + cur = NEXT(cur); + break; + case -1: + if ( !prev ) { + /* cur = head */ + prev = new; + c2 = C(m2); + DMAR(c2,mul,0,nd_mod,c1); + C(prev) = c1; + NEXT(prev) = head; + head = prev; + } else { + c2 = C(m2); + DMAR(c2,mul,0,nd_mod,c1); + C(new) = c1; + NEXT(prev) = new; + NEXT(new) = cur; + prev = new; + } + NEWNM(new); + m2 = NEXT(m2); + break; + } } - if ( mr0 ) { - NEXT(mr) = 0; MKDP(p->nv,mr0,*rp); (*rp)->sugar = p->sugar; + FREENM(new); + if ( head ) { + BDY(p1) = head; + p1->sugar = MAX(p1->sugar,p2->sugar+td); + return p1; + } else { + FREEND(p1); + return 0; + } + + } +} + +ND nd_sp(ND_pairs p) +{ + NM m; + ND p1,p2,t1,t2; + unsigned int *lcm; + int td; + + p1 = nps[p->i1]; + p2 = nps[p->i2]; + lcm = p->lcm; + td = p->td; + NEWNM(m); + C(m) = HC(p2); m->td = td-HTD(p1); ndl_sub(lcm,HDL(p1),m->dl); NEXT(m) = 0; + t1 = nd_mul_nm(p1,m); + C(m) = nd_mod-HC(p1); m->td = td-HTD(p2); ndl_sub(lcm,HDL(p2),m->dl); + t2 = nd_mul_nm(p2,m); + FREENM(m); + return nd_add(t1,t2); +} + +ND nd_find_reducer(ND g) +{ + NM m; + ND r,p; + int i,c1,c2,c; + + for ( i = 0; i < nd_psn; i++ ) { + p = nps[i]; + if ( HTD(g) >= HTD(p) && ndl_reducible(HDL(g),HDL(p)) ) { +#if 1 + NEWNM(m); + c1 = invm(HC(p),nd_mod); + c2 = nd_mod-HC(g); + DMAR(c1,c2,0,nd_mod,c); + C(m) = c; + m->td = HTD(g)-HTD(p); + ndl_sub(HDL(g),HDL(p),m->dl); + NEXT(m) = 0; + r = nd_mul_nm(p,m); + FREENM(m); + r->sugar = m->td + p->sugar; + return r; +#else + return p; +#endif + } + } + return 0; +} + +ND nd_find_monic_reducer(ND g) +{ + int *d; + ND p,r; + int i; + + for ( i = 0; i < nd_psn; i++ ) { + p = nps[i]; + if ( HTD(g) >= HTD(p) && ndl_reducible(HDL(g),HDL(p)) ) { + d = (int *)ALLOCA(nd_wpd*sizeof(int)); + ndl_sub(HDL(g),HDL(p),d); + r = nd_mul_term(p,HTD(g)-HTD(p),d); + return r; + } + } + return 0; +} + +ND nd_add(ND p1,ND p2) +{ + int n,c; + int t; + ND r; + NM m1,m2,mr0,mr,s; + + if ( !p1 ) + return p2; + else if ( !p2 ) + return p1; + else { + for ( n = NV(p1), m1 = BDY(p1), m2 = BDY(p2), mr0 = 0; m1 && m2; ) { + if ( m1->td > m2->td ) + c = 1; + else if ( m1->td < m2->td ) + c = -1; + else + c = ndl_compare(m1->dl,m2->dl); + switch ( c ) { + case 0: + t = ((C(m1))+(C(m2))) - nd_mod; + if ( t < 0 ) + t += nd_mod; + s = m1; m1 = NEXT(m1); + if ( t ) { + NEXTNM2(mr0,mr,s); C(mr) = (t); + } else { + FREENM(s); + } + s = m2; m2 = NEXT(m2); FREENM(s); + break; + case 1: + s = m1; m1 = NEXT(m1); NEXTNM2(mr0,mr,s); + break; + case -1: + s = m2; m2 = NEXT(m2); NEXTNM2(mr0,mr,s); + break; + } + } + if ( !mr0 ) + if ( m1 ) + mr0 = m1; + else if ( m2 ) + mr0 = m2; + else + return 0; + else if ( m1 ) + NEXT(mr) = m1; + else if ( m2 ) + NEXT(mr) = m2; + else + NEXT(mr) = 0; + BDY(p1) = mr0; + p1->sugar = MAX(p1->sugar,p2->sugar); + FREEND(p2); + return p1; + } +} + +ND nd_mul_nm(ND p,NM m0) +{ + NM m,mr,mr0; + unsigned int *d,*dt,*dm; + int c,n,td,i,c1,c2; + int *pt,*p1,*p2; + ND r; + + if ( !p ) + return 0; + else { + n = NV(p); m = BDY(p); + d = m0->dl; td = m0->td; c = C(m0); + mr0 = 0; + for ( ; m; m = NEXT(m) ) { + NEXTNM(mr0,mr); + c1 = C(m); + DMAR(c1,c,0,nd_mod,c2); + C(mr) = c2; + mr->td = m->td+td; + ndl_add(m->dl,d,mr->dl); + } + NEXT(mr) = 0; + MKND(NV(p),mr0,r); + r->sugar = p->sugar + td; + return r; + } +} + +ND nd_mul_term(ND p,int td,unsigned int *d) +{ + NM m,mr,mr0; + int c,n; + ND r; + + if ( !p ) + return 0; + else { + n = NV(p); m = BDY(p); + mr0 = 0; + for ( ; m; m = NEXT(m) ) { + NEXTNM(mr0,mr); + C(mr) = C(m); + mr->td = m->td+td; + ndl_add(m->dl,d,mr->dl); + } + NEXT(mr) = 0; + MKND(NV(p),mr0,r); + r->sugar = p->sugar + td; + return r; + } +} + +#if 1 +ND nd_nf(ND g,int full) +{ + ND p,d,red; + NM m,mrd,tail; + int n,sugar,psugar; + + if ( !g ) + return 0; + sugar = g->sugar; + n = NV(g); + for ( d = 0; g; ) { + red = nd_find_reducer(g); + if ( red ) { +#if 1 + g = nd_add(g,red); + sugar = MAX(sugar,red->sugar); +#else + psugar = (HTD(g)-HTD(red))+red->sugar; + g = nd_reduce(g,red); + sugar = MAX(sugar,psugar); +#endif + } else if ( !full ) + return g; + else { + m = BDY(g); + if ( NEXT(m) ) { + BDY(g) = NEXT(m); NEXT(m) = 0; + } else { + FREEND(g); g = 0; + } + if ( d ) { + NEXT(tail)=m; + tail=m; + } else { + MKND(n,m,d); + tail = BDY(d); + } + } + } + if ( d ) + d->sugar = sugar; + return d; +} +#else + +ND nd_remove_head(ND p) +{ + NM m; + + m = BDY(p); + if ( !NEXT(m) ) { + FREEND(p); + p = 0; + } else + BDY(p) = NEXT(m); + FREENM(m); + return p; +} + +PGeoBucket create_pbucket() +{ + PGeoBucket g; + + g = CALLOC(1,sizeof(struct oPGeoBucket)); + g->m = -1; + return g; +} + +void add_pbucket(PGeoBucket g,ND d) +{ + int l,k,m; + + l = nd_length(d); + for ( k = 0, m = 1; l > m; k++, m <<= 2 ); + /* 4^(k-1) < l <= 4^k */ + d = nd_add(g->body[k],d); + for ( ; d && nd_length(d) > 1<<(2*k); k++ ) { + g->body[k] = 0; + d = nd_add(g->body[k+1],d); + } + g->body[k] = d; + g->m = MAX(g->m,k); +} + +int head_pbucket(PGeoBucket g) +{ + int j,i,c,k,nv,sum; + unsigned int *di,*dj; + ND gi,gj; + + k = g->m; + while ( 1 ) { + j = -1; + for ( i = 0; i <= k; i++ ) { + if ( !(gi = g->body[i]) ) + continue; + if ( j < 0 ) { + j = i; + gj = g->body[j]; + dj = HDL(gj); + sum = HC(gj); + } else { + di = HDL(gi); + nv = NV(gi); + if ( HTD(gi) > HTD(gj) ) + c = 1; + else if ( HTD(gi) < HTD(gj) ) + c = -1; + else + c = ndl_compare(di,dj); + if ( c > 0 ) { + if ( sum ) + HC(gj) = sum; + else + g->body[j] = nd_remove_head(gj); + j = i; + gj = g->body[j]; + dj = HDL(gj); + sum = HC(gj); + } else if ( c == 0 ) { + sum = sum+HC(gi)-nd_mod; + if ( sum < 0 ) + sum += nd_mod; + g->body[i] = nd_remove_head(gi); + } + } + } + if ( j < 0 ) + return -1; + else if ( sum ) { + HC(gj) = sum; + return j; } else - *rp = 0; + g->body[j] = nd_remove_head(gj); } } -void _dp_sp_mod(p1,p2,mod,rp) -DP p1,p2; -int mod; -DP *rp; +ND normalize_pbucket(PGeoBucket g) { - int i,n,td; - int *w; - DL d1,d2,d; - MP m; - DP t,s,u; + int i; + ND r,t; - n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl; - w = (int *)ALLOCA(n*sizeof(int)); - for ( i = 0, td = 0; i < n; i++ ) { - w[i] = MAX(d1->d[i],d2->d[i]); td += w[i]; + r = 0; + for ( i = 0; i <= g->m; i++ ) + r = nd_add(r,g->body[i]); + return r; +} + +ND nd_nf(ND g,int full) +{ + ND u,p,d,red; + NODE l; + NM m,mrd; + int sugar,psugar,n,h_reducible,h; + PGeoBucket bucket; + + if ( !g ) { + return 0; } - NEWDL(d,n); d->td = td - d1->td; - for ( i = 0; i < n; i++ ) - d->d[i] = w[i] - d1->d[i]; - NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0; - MKDP(n,m,s); s->sugar = d->td; _mulmd(CO,mod,p1,s,&t); - NEWDL(d,n); d->td = td - d2->td; - for ( i = 0; i < n; i++ ) - d->d[i] = w[i] - d2->d[i]; - NEWMP(m); m->dl = d; m->c = STOI(mod - ITOS(BDY(p1)->c)); NEXT(m) = 0; - MKDP(n,m,s); s->sugar = d->td; _mulmd(CO,mod,p2,s,&u); - _addmd(CO,mod,t,u,rp); + sugar = g->sugar; + n = g->nv; + bucket = create_pbucket(); + add_pbucket(bucket,g); + d = 0; + while ( 1 ) { + h = head_pbucket(bucket); + if ( h < 0 ) { + if ( d ) + d->sugar = sugar; + return d; + } + g = bucket->body[h]; + red = nd_find_reducer(g); + if ( red ) { + bucket->body[h] = nd_remove_head(g); + red = nd_remove_head(red); + add_pbucket(bucket,red); + sugar = MAX(sugar,red->sugar); + } else if ( !full ) { + g = normalize_pbucket(bucket); + if ( g ) + g->sugar = sugar; + return g; + } else { + m = BDY(g); + if ( NEXT(m) ) { + BDY(g) = NEXT(m); NEXT(m) = 0; + } else { + FREEND(g); g = 0; + } + bucket->body[h] = g; + NEXT(m) = 0; + if ( d ) { + for ( mrd = BDY(d); NEXT(mrd); mrd = NEXT(mrd) ); + NEXT(mrd) = m; + } else { + MKND(n,m,d); + } + } + } } +#endif -void _dp_sp_component_mod(p1,p2,mod,f1,f2) -DP p1,p2; -int mod; -DP *f1,*f2; +NODE nd_gb(NODE f) { - int i,n,td; - int *w; - DL d1,d2,d; - MP m; - DP t,s,u; + int i,nh; + NODE r,g,gall; + ND_pairs d; + ND_pairs l; + ND h,nf; - n = p1->nv; d1 = BDY(p1)->dl; d2 = BDY(p2)->dl; - w = (int *)ALLOCA(n*sizeof(int)); - for ( i = 0, td = 0; i < n; i++ ) { - w[i] = MAX(d1->d[i],d2->d[i]); td += w[i]; + for ( gall = g = 0, d = 0, r = f; r; r = NEXT(r) ) { + i = (int)BDY(r); + d = update_pairs(d,g,i); + g = update_base(g,i); + gall = append_one(gall,i); } - NEWDL(d,n); d->td = td - d1->td; - for ( i = 0; i < n; i++ ) - d->d[i] = w[i] - d1->d[i]; - NEWMP(m); m->dl = d; m->c = BDY(p2)->c; NEXT(m) = 0; - MKDP(n,m,s); s->sugar = d->td; _mulmd(CO,mod,p1,s,f1); - NEWDL(d,n); d->td = td - d2->td; - for ( i = 0; i < n; i++ ) - d->d[i] = w[i] - d2->d[i]; - NEWMP(m); m->dl = d; m->c = BDY(p1)->c; NEXT(m) = 0; - MKDP(n,m,s); s->sugar = d->td; _mulmd(CO,mod,p2,s,f2); + while ( d ) { +#if 0 + ndp_print(d); +#endif + l = nd_minp(d,&d); + h = nd_sp(l); + nf = nd_nf(h,!Top); + if ( nf ) { + printf("+"); fflush(stdout); +#if 0 + ndl_print(HDL(nf)); fflush(stdout); +#endif + nh = nd_newps(nf); + d = update_pairs(d,g,nh); + g = update_base(g,nh); + gall = append_one(gall,nh); + } else { + printf("."); fflush(stdout); + } + } + return g; } -void _printdp(d) -DP d; +ND_pairs update_pairs( ND_pairs d, NODE /* of index */ g, int t) { - int n,i; - MP m; + ND_pairs d1,nd,cur,head,prev; + + if ( !g ) return d; + d = crit_B(d,t); + d1 = nd_newpairs(g,t); + d1 = crit_M(d1); + d1 = crit_F(d1); + prev = 0; cur = head = d1; + while ( cur ) { + if ( crit_2( cur->i1,cur->i2 ) ) { + if ( !prev ) { + head = cur = NEXT(cur); + } else { + cur = NEXT(prev) = NEXT(cur); + } + } else { + prev = cur; + cur = NEXT(cur); + } + } + if ( !d ) + return head; + else { + nd = d; + while ( NEXT(nd) ) + nd = NEXT(nd); + NEXT(nd) = head; + return d; + } +} + +ND_pairs nd_newpairs( NODE g, int t ) +{ + NODE h; + unsigned int *dl; + int td,ts,s; + ND_pairs r,r0; + + dl = HDL(nps[t]); + td = HTD(nps[t]); + ts = nps[t]->sugar - td; + for ( r0 = 0, h = g; h; h = NEXT(h) ) { + NEXTND_pairs(r0,r); + r->i1 = (int)BDY(h); + r->i2 = t; + ndl_lcm(HDL(nps[r->i1]),dl,r->lcm); + r->td = ndl_td(r->lcm); + s = nps[r->i1]->sugar-HTD(nps[r->i1]); + r->sugar = MAX(s,ts) + r->td; + } + NEXT(r) = 0; + return r0; +} + +ND_pairs crit_B( ND_pairs d, int s ) +{ + ND_pairs cur,head,prev; + unsigned int *t,*tl,*lcm; + int td,tdl; + + if ( !d ) return 0; + t = HDL(nps[s]); + prev = 0; + head = cur = d; + lcm = (unsigned int *)ALLOCA(nd_wpd*sizeof(unsigned int)); + while ( cur ) { + tl = cur->lcm; + if ( ndl_reducible(tl,t) + && (ndl_lcm(HDL(nps[cur->i1]),t,lcm),!ndl_equal(lcm,tl)) + && (ndl_lcm(HDL(nps[cur->i2]),t,lcm),!ndl_equal(lcm,tl)) ) { + if ( !prev ) { + head = cur = NEXT(cur); + } else { + cur = NEXT(prev) = NEXT(cur); + } + } else { + prev = cur; + cur = NEXT(cur); + } + } + return head; +} + +ND_pairs crit_M( ND_pairs d1 ) +{ + ND_pairs e,d2,d3,dd,p; + unsigned int *id,*jd; + int itd,jtd; + + for ( dd = 0, e = d1; e; e = d3 ) { + if ( !(d2 = NEXT(e)) ) { + NEXT(e) = dd; + return e; + } + id = e->lcm; + itd = e->td; + for ( d3 = 0; d2; d2 = p ) { + p = NEXT(d2), + jd = d2->lcm; + jtd = d2->td; + if ( jtd == itd ) + if ( id == jd ); + else if ( ndl_reducible(jd,id) ) continue; + else if ( ndl_reducible(id,jd) ) goto delit; + else ; + else if ( jtd > itd ) + if ( ndl_reducible(jd,id) ) continue; + else ; + else if ( ndl_reducible(id,jd ) ) goto delit; + NEXT(d2) = d3; + d3 = d2; + } + NEXT(e) = dd; + dd = e; + continue; + /**/ + delit: NEXT(d2) = d3; + d3 = d2; + for ( ; p; p = d2 ) { + d2 = NEXT(p); + NEXT(p) = d3; + d3 = p; + } + } + return dd; +} + +ND_pairs crit_F( ND_pairs d1 ) +{ + ND_pairs rest, head; + ND_pairs last, p, r, w; + int s; + + for ( head = last = 0, p = d1; NEXT(p); ) { + r = w = equivalent_pairs(p,&rest); + s = r->sugar; + while ( w = NEXT(w) ) + if ( crit_2(w->i1,w->i2) ) { + r = w; + break; + } else if ( w->sugar < s ) { + r = w; + s = r->sugar; + } + if ( last ) NEXT(last) = r; + else head = r; + NEXT(last = r) = 0; + p = rest; + if ( !p ) return head; + } + if ( !last ) return p; + NEXT(last) = p; + return head; +} + +int crit_2( int dp1, int dp2 ) +{ + return ndl_disjoint(HDL(nps[dp1]),HDL(nps[dp2])); +} + +static ND_pairs equivalent_pairs( ND_pairs d1, ND_pairs *prest ) +{ + ND_pairs w,p,r,s; + unsigned int *d; + int td; + + w = d1; + d = w->lcm; + td = w->td; + s = NEXT(w); + NEXT(w) = 0; + for ( r = 0; s; s = p ) { + p = NEXT(s); + if ( td == s->td && ndl_equal(d,s->lcm) ) { + NEXT(s) = w; + w = s; + } else { + NEXT(s) = r; + r = s; + } + } + *prest = r; + return w; +} + +NODE update_base(NODE nd,int ndp) +{ + unsigned int *dl, *dln; + NODE last, p, head; + int td,tdn; + + dl = HDL(nps[ndp]); + td = HTD(nps[ndp]); + for ( head = last = 0, p = nd; p; ) { + dln = HDL(nps[(int)BDY(p)]); + tdn = HTD(nps[(int)BDY(p)]); + if ( tdn >= td && ndl_reducible( dln, dl ) ) { + p = NEXT(p); + if ( last ) NEXT(last) = p; + } else { + if ( !last ) head = p; + p = NEXT(last = p); + } + } + head = append_one(head,ndp); + return head; +} + +ND_pairs nd_minp( ND_pairs d, ND_pairs *prest ) +{ + ND_pairs m,ml,p,l; + unsigned int *lcm; + int s,td,len,tlen,c; + + if ( !(p = NEXT(m = d)) ) { + *prest = p; + NEXT(m) = 0; + return m; + } + lcm = m->lcm; + s = m->sugar; + td = m->td; + len = nd_length(nps[m->i1])+nd_length(nps[m->i2]); + for ( ml = 0, l = m; p; p = NEXT(l = p) ) { + if (p->sugar < s) + goto find; + else if ( p->sugar == s ) { + if ( p->td < td ) + goto find; + else if ( p->td == td ) { + c = ndl_compare(p->lcm,lcm); + if ( c < 0 ) + goto find; + else if ( c == 0 ) { + tlen = nd_length(nps[p->i1])+nd_length(nps[p->i2]); + if ( tlen < len ) + goto find; + } + } + } + continue; +find: + ml = l; + m = p; + lcm = m->lcm; + s = m->sugar; + td = m->td; + len = tlen; + } + if ( !ml ) *prest = NEXT(m); + else { + NEXT(ml) = NEXT(m); + *prest = d; + } + NEXT(m) = 0; + return m; +} + +int nd_newps(ND a) +{ + if ( nd_psn == nd_pslen ) { + nd_pslen *= 2; + nps = (ND *)REALLOC((char *)nps,nd_pslen*sizeof(ND)); + } + nd_monic(a); + nps[nd_psn] = a; + return nd_psn++; +} + +NODE NODE_sortb(NODE f,int); +ND dptond(DP); +DP ndtodp(ND); + +NODE nd_setup(NODE f) +{ + int i; + NODE s,s0,f0; + +#if 0 + f0 = f = NODE_sortb(f,1); +#endif + nd_psn = length(f); nd_pslen = 2*nd_psn; + nps = (ND *)MALLOC(nd_pslen*sizeof(ND)); + nd_bpe = 4; + nd_epw = (sizeof(unsigned int)*8)/nd_bpe; + nd_wpd = nd_nvar/nd_epw+(nd_nvar%nd_epw?1:0); + if ( nd_bpe < 32 ) { + nd_mask0 = (1<id ) + error("nd_gr : unsupported order"); + switch ( ord->ord.simple ) { + case 0: + is_rlex = 1; + break; + case 1: + is_rlex = 0; + break; + default: + error("nd_gr : unsupported order"); + } + initd(ord); + nd_mod = m; + for ( fd0 = 0, t = BDY(f); t; t = NEXT(t) ) { + ptod(CO,vv,(P)BDY(t),&b); + _dp_mod(b,m,0,&c); + if ( c ) { + NEXTNODE(fd0,fd); BDY(fd) = (pointer)c; + } + } + if ( fd0 ) NEXT(fd) = 0; + s = nd_setup(fd0); + x = nd_gb(s); +#if 0 + x = nd_reduceall(x,m); +#endif + for ( r0 = 0; x; x = NEXT(x) ) { + NEXTNODE(r0,r); + a = ndtodp(nps[(int)BDY(x)]); + _dtop_mod(CO,vv,a,(P *)&BDY(r)); + } + if ( r0 ) NEXT(r) = 0; + MKLIST(*rp,r0); +} + +void dltondl(int n,DL dl,unsigned int *r) +{ + unsigned int *d; + int i; + + d = dl->d; + bzero(r,nd_wpd*sizeof(unsigned int)); + if ( is_rlex ) + for ( i = 0; i < n; i++ ) + r[(n-1-i)/nd_epw] |= (d[i]<<((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe)); + else + for ( i = 0; i < n; i++ ) + r[i/nd_epw] |= d[i]<<((nd_epw-(i%nd_epw)-1)*nd_bpe); +} + +DL ndltodl(int n,int td,unsigned int *ndl) +{ DL dl; + int *d; + int i; - if ( !d ) { - printf("0"); return; + NEWDL(dl,n); + dl->td = td; + d = dl->d; + if ( is_rlex ) + for ( i = 0; i < n; i++ ) + d[i] = (ndl[(n-1-i)/nd_epw]>>((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe)) + &((1<>((nd_epw-(i%nd_epw)-1)*nd_bpe)) + &((1<c = ITOS(t->c); + m->td = t->dl->td; + dltondl(n,t->dl,m->dl); } - for ( n = d->nv, m = BDY(d); m; m = NEXT(m) ) { - printf("%d*<<",ITOS(m->c)); - for ( i = 0, dl = m->dl; i < n-1; i++ ) - printf("%d,",dl->d[i]); - printf("%d",dl->d[i]); - printf(">>"); - if ( NEXT(m) ) - printf("+"); + NEXT(m) = 0; + MKND(n,m0,d); + d->nv = n; + d->sugar = p->sugar; + return d; +} + +DP ndtodp(ND p) +{ + DP d; + MP m0,m; + NM t; + int n; + + if ( !p ) + return 0; + n = NV(p); + m0 = 0; + for ( t = BDY(p); t; t = NEXT(t) ) { + NEXTMP(m0,m); + m->c = STOI(t->c); + m->dl = ndltodl(n,t->td,t->dl); + } + NEXT(m) = 0; + MKDP(n,m0,d); + d->sugar = p->sugar; + return d; +} + +void ndl_print(unsigned int *dl) +{ + int n; + int i; + + n = nd_nvar; + printf("<<"); + if ( is_rlex ) + for ( i = 0; i < n; i++ ) + printf(i==n-1?"%d":"%d,", + (dl[(n-1-i)/nd_epw]>>((nd_epw-((n-1-i)%nd_epw)-1)*nd_bpe)) + &((1<>((nd_epw-(i%nd_epw)-1)*nd_bpe)) + &((1<>"); +} + +void nd_print(ND p) +{ + NM m; + + if ( !p ) + printf("0\n"); + else { + for ( m = BDY(p); m; m = NEXT(m) ) { + printf("+%d*",m->c); + ndl_print(m->dl); + } + printf("\n"); + } +} + +void ndp_print(ND_pairs d) +{ + ND_pairs t; + + for ( t = d; t; t = NEXT(t) ) { + printf("%d,%d ",t->i1,t->i2); + } + printf("\n"); +} + +void nd_monic(ND p) +{ + if ( !p ) + return; + else + nd_mul_c(p,invm(HC(p),nd_mod)); +} + +void nd_mul_c(ND p,int mul) +{ + NM m; + int c,c1; + + if ( !p ) + return; + for ( m = BDY(p); m; m = NEXT(m) ) { + c1 = C(m); + DMAR(c1,mul,0,nd_mod,c); + C(m) = c; } }