=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/NEZ.c,v retrieving revision 1.1 retrieving revision 1.7 diff -u -p -r1.1 -r1.7 --- OpenXM_contrib2/asir2000/engine/NEZ.c 1999/12/03 07:39:08 1.1 +++ OpenXM_contrib2/asir2000/engine/NEZ.c 2018/03/29 01:32:51 1.7 @@ -1,4 +1,52 @@ -/* $OpenXM: OpenXM/src/asir99/engine/NEZ.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/NEZ.c,v 1.6 2010/02/19 08:27:19 noro Exp $ +*/ #include "ca.h" void nezgcdnpz(vl,ps,m,pr) @@ -6,108 +54,111 @@ VL vl; P *ps,*pr; int m; { - P t,s,mg; - VL tvl,svl,avl,nvl; - int i,j,k; - N c; - Q cq; - P *pl,*pm; - DCP *ml; - Q *cl; - P **cp; - int *cn; + P t,s,mg; + VL tvl,svl,avl,nvl; + int i,j,k; + N c; + Q cq; + P *pl,*pm; + DCP *ml; + Q *cl; + P **cp; + int *cn; - pl = (P *)ALLOCA(m*sizeof(P)); - ml = (DCP *)ALLOCA(m*sizeof(DCP)); - cl = (Q *)ALLOCA(m*sizeof(Q)); - for ( i = 0; i < m; i++ ) - monomialfctr(vl,ps[i],&pl[i],&ml[i]); - gcdmonomial(vl,ml,m,&mg); - for ( i = 0; i < m; i++ ) { - ptozp(pl[i],1,&cl[i],&t); pl[i] = t; - } - for ( i = 1, cq = cl[0]; i < m; i++ ) { - gcdn(NM(cl[i]),NM(cq),&c); NTOQ(c,1,cq); - } - for ( i = 0; i < m; i++ ) - if ( NUM(pl[i]) ) { - mulp(vl,(P)cq,mg,pr); return; - } - for ( i = 0, nvl = vl, avl = 0; nvl && i < m; i++ ) { - clctv(vl,pl[i],&tvl); - intersectv(nvl,tvl,&svl); nvl = svl; - mergev(vl,avl,tvl,&svl); avl = svl; - } - if ( !nvl ) { - mulp(vl,(P)cq,mg,pr); return; - } - if ( !NEXT(avl) ) { - nuezgcdnpzmain(vl,pl,m,&t); mulp(vl,t,(P)cq,&s); mulp(vl,s,mg,pr); - return; - } - for ( tvl = nvl, i = 0; tvl; tvl = NEXT(tvl), i++ ); - for ( tvl = avl, j = 0; tvl; tvl = NEXT(tvl), j++ ); - if ( i == j ) { - /* all the pl[i]'s have the same variables */ - sortplistbyhomdeg(pl,m); nezgcdnpzmain(nvl,pl,m,&t); + if ( m == 1 ) { + *pr = ps[0]; return; + } + pl = (P *)ALLOCA(m*sizeof(P)); + ml = (DCP *)ALLOCA(m*sizeof(DCP)); + cl = (Q *)ALLOCA(m*sizeof(Q)); + for ( i = 0; i < m; i++ ) + monomialfctr(vl,ps[i],&pl[i],&ml[i]); + gcdmonomial(vl,ml,m,&mg); + for ( i = 0; i < m; i++ ) { + ptozp(pl[i],1,&cl[i],&t); pl[i] = t; + } + for ( i = 1, cq = cl[0]; i < m; i++ ) { + gcdn(NM(cl[i]),NM(cq),&c); NTOQ(c,1,cq); + } + for ( i = 0; i < m; i++ ) + if ( NUM(pl[i]) ) { + mulp(vl,(P)cq,mg,pr); return; + } + for ( i = 0, nvl = vl, avl = 0; nvl && i < m; i++ ) { + clctv(vl,pl[i],&tvl); + intersectv(nvl,tvl,&svl); nvl = svl; + mergev(vl,avl,tvl,&svl); avl = svl; + } + if ( !nvl ) { + mulp(vl,(P)cq,mg,pr); return; + } + if ( !NEXT(avl) ) { + nuezgcdnpzmain(vl,pl,m,&t); mulp(vl,t,(P)cq,&s); mulp(vl,s,mg,pr); + return; + } + for ( tvl = nvl, i = 0; tvl; tvl = NEXT(tvl), i++ ); + for ( tvl = avl, j = 0; tvl; tvl = NEXT(tvl), j++ ); + if ( i == j ) { + /* all the pl[i]'s have the same variables */ + sortplistbyhomdeg(pl,m); nezgcdnpzmain(nvl,pl,m,&t); #if 0 - /* search the minimal degree poly */ - for ( i = 0; i < m; i++ ) { - for ( tvl = nvl; tvl; tvl = NEXT(tvl) ) { + /* search the minimal degree poly */ + for ( i = 0; i < m; i++ ) { + for ( tvl = nvl; tvl; tvl = NEXT(tvl) ) { dt = getdeg(tvl->v,pl[i]); if ( tvl == nvl || dt < d1 ) { v1 = tvl->v; d1 = dt; } - } - if ( i == 0 || d1 < d ) { - v = v1; d = d1; j = i; - } - } - t = pl[0]; pl[0] = pl[j]; pl[j] = t; - if ( v != nvl->v ) { - reordvar(nvl,v,&mvl); - for ( i = 0; i < m; i++ ) { - reorderp(mvl,nvl,pl[i],&t); pl[i] = t; - } - nezgcdnpzmain(mvl,pl,m,&s); reorderp(nvl,mvl,s,&t); - } else - nezgcdnpzmain(nvl,pl,m,&t); + } + if ( i == 0 || d1 < d ) { + v = v1; d = d1; j = i; + } + } + t = pl[0]; pl[0] = pl[j]; pl[j] = t; + if ( v != nvl->v ) { + reordvar(nvl,v,&mvl); + for ( i = 0; i < m; i++ ) { + reorderp(mvl,nvl,pl[i],&t); pl[i] = t; + } + nezgcdnpzmain(mvl,pl,m,&s); reorderp(nvl,mvl,s,&t); + } else + nezgcdnpzmain(nvl,pl,m,&t); #endif - } else { - cp = (P **)ALLOCA(m*sizeof(P *)); - cn = (int *)ALLOCA(m*sizeof(int)); - for ( i = 0; i < m; i++ ) { - cp[i] = (P *)ALLOCA(lengthp(pl[i])*sizeof(P)); - cn[i] = pcoef(vl,nvl,pl[i],cp[i]); - } - for ( i = j = 0; i < m; i++ ) - j += cn[i]; - pm = (P *)ALLOCA(j*sizeof(P)); - for ( i = k = 0; i < m; i++ ) - for ( j = 0; j < cn[i]; j++ ) - pm[k++] = cp[i][j]; - nezgcdnpz(vl,pm,k,&t); - } - mulp(vl,t,(P)cq,&s); mulp(vl,s,mg,pr); + } else { + cp = (P **)ALLOCA(m*sizeof(P *)); + cn = (int *)ALLOCA(m*sizeof(int)); + for ( i = 0; i < m; i++ ) { + cp[i] = (P *)ALLOCA(lengthp(pl[i])*sizeof(P)); + cn[i] = pcoef(vl,nvl,pl[i],cp[i]); + } + for ( i = j = 0; i < m; i++ ) + j += cn[i]; + pm = (P *)ALLOCA(j*sizeof(P)); + for ( i = k = 0; i < m; i++ ) + for ( j = 0; j < cn[i]; j++ ) + pm[k++] = cp[i][j]; + nezgcdnpz(vl,pm,k,&t); + } + mulp(vl,t,(P)cq,&s); mulp(vl,s,mg,pr); } void sortplistbyhomdeg(p,n) P *p; int n; { - int i,j,k; - int *l; - P t; + int i,j,k; + int *l; + P t; - l = (int *)ALLOCA(n*sizeof(int)); - for ( i = 0; i < n; i++ ) - l[i] = homdeg(p[i]); - for ( i = 0; i < n; i++ ) - for ( j = i + 1; j < n; j++ ) - if ( l[j] < l[i] ) { - t = p[i]; p[i] = p[j]; p[j] = t; - k = l[i]; l[i] = l[j]; l[j] = k; - } + l = (int *)ALLOCA(n*sizeof(int)); + for ( i = 0; i < n; i++ ) + l[i] = homdeg(p[i]); + for ( i = 0; i < n; i++ ) + for ( j = i + 1; j < n; j++ ) + if ( l[j] < l[i] ) { + t = p[i]; p[i] = p[j]; p[j] = t; + k = l[i]; l[i] = l[j]; l[j] = k; + } } void nuezgcdnpzmain(vl,ps,m,r) @@ -116,22 +167,22 @@ P *ps; int m; P *r; { - P *tps; - P f,t; - int i; + P *tps; + P f,t; + int i; - for ( i = 0; i < m; i++ ) - if ( NUM(ps[i]) ) { - *r = (P)ONE; return; - } - tps = (P *) ALLOCA(m*sizeof(P)); - for ( i = 0; i < m; i++ ) - tps[i] = ps[i]; - sortplist(tps,m); - for ( i = 1, f = tps[0]; i < m && !NUM(f); i++ ) { - uezgcdpz(vl,f,tps[i],&t); f = t; - } - *r = f; + for ( i = 0; i < m; i++ ) + if ( NUM(ps[i]) ) { + *r = (P)ONE; return; + } + tps = (P *) ALLOCA(m*sizeof(P)); + for ( i = 0; i < m; i++ ) + tps[i] = ps[i]; + sortplist(tps,m); + for ( i = 1, f = tps[0]; i < m && !NUM(f); i++ ) { + uezgcdpz(vl,f,tps[i],&t); f = t; + } + *r = f; } void gcdmonomial(vl,dcl,m,r) @@ -140,42 +191,42 @@ DCP *dcl; int m; P *r; { - int i,j,n; - P g,x,s,t; - Q d; - DCP dc; - VN vn; + int i,j,n; + P g,x,s,t; + Q d; + DCP dc; + VN vn; - for ( i = 0; i < m; i++ ) - if ( !dcl[i] ) { - *r = (P)ONE; return; - } - for ( n = 0, dc = dcl[0]; dc; dc = NEXT(dc), n++ ); - vn = (VN)ALLOCA(n*sizeof(struct oVN)); - for ( i = 0, dc = dcl[0]; i < n; dc = NEXT(dc), i++ ) { - vn[i].v = VR(COEF(dc)); vn[i].n = QTOS(DEG(dc)); - } - for ( i = 1; i < m; i++ ) { - for ( j = 0; j < n; j++ ) { - for ( dc = dcl[i]; dc; dc = NEXT(dc) ) - if ( VR(COEF(dc)) == vn[j].v ) { - vn[j].n = MIN(vn[j].n,QTOS(DEG(dc))); break; - } - if ( !dc ) - vn[j].n = 0; - } - for ( j = n-1; j >= 0 && !vn[j].n; j-- ); - if ( j < 0 ) { - *r = (P)ONE; return; - } else - n = j+1; - } - for ( j = 0, g = (P)ONE; j < n; j++ ) - if ( vn[j].n ) { - MKV(vn[j].v,x); STOQ(vn[j].n,d); - pwrp(vl,x,d,&t); mulp(vl,t,g,&s); g = s; - } - *r = g; + for ( i = 0; i < m; i++ ) + if ( !dcl[i] ) { + *r = (P)ONE; return; + } + for ( n = 0, dc = dcl[0]; dc; dc = NEXT(dc), n++ ); + vn = (VN)ALLOCA(n*sizeof(struct oVN)); + for ( i = 0, dc = dcl[0]; i < n; dc = NEXT(dc), i++ ) { + vn[i].v = VR(COEF(dc)); vn[i].n = QTOS(DEG(dc)); + } + for ( i = 1; i < m; i++ ) { + for ( j = 0; j < n; j++ ) { + for ( dc = dcl[i]; dc; dc = NEXT(dc) ) + if ( VR(COEF(dc)) == vn[j].v ) { + vn[j].n = MIN(vn[j].n,QTOS(DEG(dc))); break; + } + if ( !dc ) + vn[j].n = 0; + } + for ( j = n-1; j >= 0 && !vn[j].n; j-- ); + if ( j < 0 ) { + *r = (P)ONE; return; + } else + n = j+1; + } + for ( j = 0, g = (P)ONE; j < n; j++ ) + if ( vn[j].n ) { + MKV(vn[j].v,x); STOQ(vn[j].n,d); + pwrp(vl,x,d,&t); mulp(vl,t,g,&s); g = s; + } + *r = g; } void nezgcdnpzmain(vl,pl,m,pr) @@ -183,26 +234,26 @@ VL vl; P *pl,*pr; int m; { - P *ppl,*pcl; - int i; - P cont,gcd,t,s; - DCP dc; - - ppl = (P *)ALLOCA(m*sizeof(P)); - pcl = (P *)ALLOCA(m*sizeof(P)); - for ( i = 0; i < m; i++ ) - pcp(vl,pl[i],&ppl[i],&pcl[i]); - nezgcdnpz(vl,pcl,m,&cont); - sqfrp(vl,ppl[0],&dc); - for ( dc = NEXT(dc), gcd = (P)ONE; dc; dc = NEXT(dc) ) { - if ( NUM(COEF(dc)) ) - continue; - nezgcdnpp(vl,dc,ppl+1,m-1,&t); - if ( NUM(t) ) - continue; - mulp(vl,gcd,t,&s); gcd = s; - } - mulp(vl,gcd,cont,pr); + P *ppl,*pcl; + int i; + P cont,gcd,t,s; + DCP dc; + + ppl = (P *)ALLOCA(m*sizeof(P)); + pcl = (P *)ALLOCA(m*sizeof(P)); + for ( i = 0; i < m; i++ ) + pcp(vl,pl[i],&ppl[i],&pcl[i]); + nezgcdnpz(vl,pcl,m,&cont); + sqfrp(vl,ppl[0],&dc); + for ( dc = NEXT(dc), gcd = (P)ONE; dc; dc = NEXT(dc) ) { + if ( NUM(COEF(dc)) ) + continue; + nezgcdnpp(vl,dc,ppl+1,m-1,&t); + if ( NUM(t) ) + continue; + mulp(vl,gcd,t,&s); gcd = s; + } + mulp(vl,gcd,cont,pr); } void nezgcdnpp(vl,dc,pl,m,r) @@ -212,35 +263,35 @@ P *pl; int m; P *r; { - int i,k; - P g,t,s,gcd; - P *pm; + int i,k; + P g,t,s,gcd; + P *pm; - nezgcdnp_sqfr_primitive(vl,COEF(dc),pl,m,&gcd); - if ( NUM(gcd) ) { - *r = (P)ONE; return; - } - pm = (P *) ALLOCA(m*sizeof(P)); - for ( i = 0; i < m; i++ ) { - divsp(vl,pl[i],gcd,&pm[i]); - if ( NUM(pm[i]) ) { - *r = gcd; return; - } - } - for ( g = gcd, k = QTOS(DEG(dc))-1; k > 0; k-- ) { - nezgcdnp_sqfr_primitive(vl,g,pm,m,&t); - if ( NUM(t) ) - break; - mulp(vl,gcd,t,&s); gcd = s; - for ( i = 0; i < m; i++ ) { - divsp(vl,pm[i],t,&s); - if ( NUM(s) ) { - *r = gcd; return; - } - pm[i] = s; - } - } - *r = gcd; + nezgcdnp_sqfr_primitive(vl,COEF(dc),pl,m,&gcd); + if ( NUM(gcd) ) { + *r = (P)ONE; return; + } + pm = (P *) ALLOCA(m*sizeof(P)); + for ( i = 0; i < m; i++ ) { + divsp(vl,pl[i],gcd,&pm[i]); + if ( NUM(pm[i]) ) { + *r = gcd; return; + } + } + for ( g = gcd, k = QTOS(DEG(dc))-1; k > 0; k-- ) { + nezgcdnp_sqfr_primitive(vl,g,pm,m,&t); + if ( NUM(t) ) + break; + mulp(vl,gcd,t,&s); gcd = s; + for ( i = 0; i < m; i++ ) { + divsp(vl,pm[i],t,&s); + if ( NUM(s) ) { + *r = gcd; return; + } + pm[i] = s; + } + } + *r = gcd; } /* @@ -256,120 +307,120 @@ VL vl; int m; P p0,*ps,*pr; { - /* variables */ - P p00,g,h,g0,h0,a0,b0; - P lp0,lgp0,lp00,lg,lg0,cbd,tq,t; - P *lc; - P *tps; - VL nvl1,nvl2,nvl,tvl; - V v; - int i,j,k,d0,dg,dg0,dmin; - VN vn0,vn1,vnt; - int nv,flag; + /* variables */ + P p00,g,h,g0,h0,a0,b0; + P lp0,lgp0,lp00,lg,lg0,cbd,tq,t; + P *lc; + P *tps; + VL nvl1,nvl2,nvl,tvl; + V v; + int i,j,k,d0,dg,dg0,dmin,z; + VN vn1; + int nv,flag,max; - /* set-up */ - if ( NUM(p0) ) { - *pr = (P) ONE; return; - } - for ( v = VR(p0), i = 0; i < m; i++ ) - if ( NUM(ps[i]) || (v != VR(ps[i])) ) { - *pr = (P)ONE; return; - } - tps = (P *) ALLOCA(m*sizeof(P)); - for ( i = 0; i < m; i++ ) - tps[i] = ps[i]; - sortplist(tps,m); - /* deg(tps[0]) <= deg(tps[1]) <= ... */ - - if ( !cmpq(DEG(DC(p0)),ONE) ) { - if ( divcheck(vl,tps,m,(P)ONE,p0) ) - *pr = p0; - else - *pr = (P)ONE; - return; - } + /* set-up */ + if ( NUM(p0) ) { + *pr = (P) ONE; return; + } + for ( v = VR(p0), i = 0; i < m; i++ ) + if ( NUM(ps[i]) || (v != VR(ps[i])) ) { + *pr = (P)ONE; return; + } + tps = (P *) ALLOCA(m*sizeof(P)); + for ( i = 0; i < m; i++ ) + tps[i] = ps[i]; + sortplist(tps,m); + /* deg(tps[0]) <= deg(tps[1]) <= ... */ + + if ( !cmpq(DEG(DC(p0)),ONE) ) { + if ( divcheck(vl,tps,m,(P)ONE,p0) ) + *pr = p0; + else + *pr = (P)ONE; + return; + } - lp0 = LC(p0); dmin = d0 = deg(v,p0); lc = (P *)ALLOCA((m+1)*sizeof(P)); - for ( lc[0] = lp0, i = 0; i < m; i++ ) - lc[i+1] = LC(tps[i]); - clctv(vl,p0,&nvl); - for ( i = 0; i < m; i++ ) { - clctv(vl,tps[i],&nvl1); mergev(vl,nvl,nvl1,&nvl2); nvl = nvl2; - } - nezgcdnpz(nvl,lc,m+1,&lg); + lp0 = LC(p0); dmin = d0 = deg(v,p0); lc = (P *)ALLOCA((m+1)*sizeof(P)); + for ( lc[0] = lp0, i = 0; i < m; i++ ) + lc[i+1] = LC(tps[i]); + clctv(vl,p0,&nvl); + for ( i = 0; i < m; i++ ) { + clctv(vl,tps[i],&nvl1); mergev(vl,nvl,nvl1,&nvl2); nvl = nvl2; + } + nezgcdnpz(nvl,lc,m+1,&lg); - mulp(nvl,p0,lg,&lgp0); k = dbound(v,lgp0)+1; cbound(nvl,lgp0,(Q *)&cbd); - for ( nv = 0, tvl = nvl; tvl; tvl = NEXT(tvl), nv++ ); - W_CALLOC(nv,struct oVN,vn0); W_CALLOC(nv,struct oVN,vnt); - W_CALLOC(nv,struct oVN,vn1); - for ( i = 0, tvl = NEXT(nvl); tvl; tvl = NEXT(tvl), i++ ) - vn1[i].v = vn0[i].v = tvl->v; + mulp(nvl,p0,lg,&lgp0); k = dbound(v,lgp0)+1; cbound(nvl,lgp0,(Q *)&cbd); + for ( nv = 0, tvl = nvl; tvl; tvl = NEXT(tvl), nv++ ); + W_CALLOC(nv,struct oVN,vn1); + for ( i = 0, tvl = NEXT(nvl); tvl; tvl = NEXT(tvl), i++ ) + vn1[i].v = tvl->v; - /* main loop */ - for ( dg = deg(v,tps[0]) + 1; ; next(vn0) ) - do { - for ( i = 0, j = 0; vn0[i].v; i++ ) - if ( vn0[i].n ) vnt[j++].v = (V)i; - vnt[j].n = 0; + /* main loop */ + /* finally, 'max' random evaluations will be generated */ + for ( max = 1, dg = deg(v,tps[0]) + 1; ; max = 2*max ) + for ( z = 0; z < max; z++ ) { + for ( i = 0; vn1[i].v; i++ ) + vn1[i].n = mt_genrand()%max; - /* find b s.t. LC(p0)(b), LC(tps[i])(b) != 0 */ - mulsgn(vn0,vnt,j,vn1); substvp(nvl,p0,vn1,&p00); - flag = (!zerovpchk(nvl,lp0,vn1) && sqfrchk(p00)); - for ( i = 0; flag && (i < m); i++ ) - flag &= (!zerovpchk(nvl,LC(tps[i]),vn1)); - if ( !flag ) - continue; + /* find b s.t. LC(p0)(b), LC(tps[i])(b) != 0 */ + + substvp(nvl,p0,vn1,&p00); + flag = (!zerovpchk(nvl,lp0,vn1) && sqfrchk(p00)); + for ( i = 0; flag && (i < m); i++ ) + flag &= (!zerovpchk(nvl,LC(tps[i]),vn1)); + if ( !flag ) + continue; - /* substitute y -> b */ - substvp(nvl,lg,vn1,&lg0); lp00 = LC(p00); - /* extended-gcd in 1-variable */ - uexgcdnp(nvl,p00,tps,m,vn1,(Q)cbd,&g0,&h0,&a0,&b0,(Q *)&tq); - if ( NUM(g0) ) { - *pr = (P)ONE; return; - } else if ( dg > ( dg0 = deg(v,g0) ) ) { - dg = dg0; - if ( dg == d0 ) { - if ( divcheck(nvl,tps,m,lp0,p0) ) { - *pr = p0; return; - } - } else if ( dg == deg(v,tps[0]) ) { - if ( divtpz(nvl,p0,tps[0],&t) && - divcheck(nvl,tps+1,m-1,LC(tps[0]),tps[0]) ) { - *pr = tps[0]; return; - } else - break; - } else { - henmv(nvl,vn1,lgp0,g0,h0,a0,b0,lg,lp0,lg0,lp00,(Q)tq,k,&g,&h); - if ( divtpz(nvl,lgp0,g,&t) && - divcheck(nvl,tps,m,lg,g) ) { - pcp(nvl,g,pr,&t); return; - } - } - } - } while ( !nextbin(vnt,j) ); + /* substitute y -> b */ + substvp(nvl,lg,vn1,&lg0); lp00 = LC(p00); + /* extended-gcd in 1-variable */ + uexgcdnp(nvl,p00,tps,m,vn1,(Q)cbd,&g0,&h0,&a0,&b0,(Q *)&tq); + if ( NUM(g0) ) { + *pr = (P)ONE; return; + } else if ( dg > ( dg0 = deg(v,g0) ) ) { + dg = dg0; + if ( dg == d0 ) { + if ( divcheck(nvl,tps,m,lp0,p0) ) { + *pr = p0; return; + } + } else if ( dg == deg(v,tps[0]) ) { + if ( divtpz(nvl,p0,tps[0],&t) && + divcheck(nvl,tps+1,m-1,LC(tps[0]),tps[0]) ) { + *pr = tps[0]; return; + } else + break; + } else { + henmv(nvl,vn1,lgp0,g0,h0,a0,b0,lg,lp0,lg0,lp00,(Q)tq,k,&g,&h); + if ( divtpz(nvl,lgp0,g,&t) && + divcheck(nvl,tps,m,lg,g) ) { + pcp(nvl,g,pr,&t); return; + } + } + } + } } void intersectv(vl1,vl2,vlp) VL vl1,vl2,*vlp; { - int i,n; - VL tvl; - VN tvn; + int i,n; + VL tvl; + VN tvn; - if ( !vl1 || !vl2 ) { - *vlp = 0; return; - } - for ( n = 0, tvl = vl1; tvl; tvl = NEXT(tvl), n++ ); - tvn = (VN) ALLOCA(n*sizeof(struct oVN)); - for ( i = 0, tvl = vl1; i < n; tvl = NEXT(tvl), i++ ) { - tvn[i].v = tvl->v; tvn[i].n = 0; - } - for ( tvl = vl2; tvl; tvl = NEXT(tvl) ) - for ( i = 0; i < n; i++ ) - if ( tvn[i].v == tvl->v ) { - tvn[i].n = 1; break; - } - vntovl(tvn,n,vlp); + if ( !vl1 || !vl2 ) { + *vlp = 0; return; + } + for ( n = 0, tvl = vl1; tvl; tvl = NEXT(tvl), n++ ); + tvn = (VN) ALLOCA(n*sizeof(struct oVN)); + for ( i = 0, tvl = vl1; i < n; tvl = NEXT(tvl), i++ ) { + tvn[i].v = tvl->v; tvn[i].n = 0; + } + for ( tvl = vl2; tvl; tvl = NEXT(tvl) ) + for ( i = 0; i < n; i++ ) + if ( tvn[i].v == tvl->v ) { + tvn[i].n = 1; break; + } + vntovl(tvn,n,vlp); } int pcoef(vl,ivl,p,cp) @@ -377,33 +428,33 @@ VL vl,ivl; P p; P *cp; { - VL nvl,tvl,svl,mvl,mvl0; - P t; + VL nvl,tvl,svl,mvl,mvl0; + P t; - if ( NUM(p) ) { - cp[0] = p; return 1; - } else { - clctv(vl,p,&nvl); - for ( tvl = nvl, mvl0 = 0; tvl; tvl = NEXT(tvl) ) { - for ( svl = ivl; svl; svl = NEXT(svl) ) - if ( tvl->v == svl->v ) - break; - if ( !svl ) { - if ( !mvl0 ) { - NEWVL(mvl0); mvl = mvl0; - } else { - NEWVL(NEXT(mvl)); mvl = NEXT(mvl); - } - mvl->v = tvl->v; - } - } - if ( mvl0 ) - NEXT(mvl) = ivl; - else - mvl0 = ivl; - reorderp(mvl0,nvl,p,&t); - return pcoef0(mvl0,ivl,t,cp); - } + if ( NUM(p) ) { + cp[0] = p; return 1; + } else { + clctv(vl,p,&nvl); + for ( tvl = nvl, mvl0 = 0; tvl; tvl = NEXT(tvl) ) { + for ( svl = ivl; svl; svl = NEXT(svl) ) + if ( tvl->v == svl->v ) + break; + if ( !svl ) { + if ( !mvl0 ) { + NEWVL(mvl0); mvl = mvl0; + } else { + NEWVL(NEXT(mvl)); mvl = NEXT(mvl); + } + mvl->v = tvl->v; + } + } + if ( mvl0 ) + NEXT(mvl) = ivl; + else + mvl0 = ivl; + reorderp(mvl0,nvl,p,&t); + return pcoef0(mvl0,ivl,t,cp); + } } int pcoef0(vl,ivl,p,cp) @@ -411,39 +462,57 @@ VL vl,ivl; P p; P *cp; { - int cn,n; - DCP dc; - V v; - VL tvl; + int cn,n; + DCP dc; + V v; + VL tvl; - if ( NUM(p) ) { - cp[0] = p; return 1; - } else { - for ( v = VR(p), tvl = ivl; tvl; tvl = NEXT(tvl) ) - if ( v == tvl->v ) - break; - if ( tvl ) { - cp[0] = p; return 1; - } else { - for ( dc = DC(p), n = 0; dc; dc = NEXT(dc) ) { - cn = pcoef0(vl,ivl,COEF(dc),cp); cp += cn; n += cn; - } - return n; - } - } + if ( NUM(p) ) { + cp[0] = p; return 1; + } else { + for ( v = VR(p), tvl = ivl; tvl; tvl = NEXT(tvl) ) + if ( v == tvl->v ) + break; + if ( tvl ) { + cp[0] = p; return 1; + } else { + for ( dc = DC(p), n = 0; dc; dc = NEXT(dc) ) { + cn = pcoef0(vl,ivl,COEF(dc),cp); cp += cn; n += cn; + } + return n; + } + } } int lengthp(p) P p; { - int n; - DCP dc; + int n; + DCP dc; - if ( NUM(p) ) - return 1; - else { - for ( dc = DC(p), n = 0; dc; dc = NEXT(dc) ) - n += lengthp(COEF(dc)); - return n; - } + if ( NUM(p) ) + return 1; + else { + for ( dc = DC(p), n = 0; dc; dc = NEXT(dc) ) + n += lengthp(COEF(dc)); + return n; + } +} + +int nonzero_const_term(P p) +{ + DCP dc; + + if ( !p ) + return 0; + else if ( NUM(p) ) + return 1; + else { + dc = DC(p); + for ( ; NEXT(dc); dc = NEXT(dc) ); + if ( DEG(dc) ) + return 0; + else + return nonzero_const_term(COEF(dc)); + } }