=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/EZ.c,v retrieving revision 1.1 retrieving revision 1.4 diff -u -p -r1.1 -r1.4 --- OpenXM_contrib2/asir2000/engine/EZ.c 1999/12/03 07:39:07 1.1 +++ OpenXM_contrib2/asir2000/engine/EZ.c 2018/03/29 01:32:51 1.4 @@ -1,4 +1,52 @@ -/* $OpenXM: OpenXM/src/asir99/engine/EZ.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/EZ.c,v 1.3 2000/08/22 05:04:03 noro Exp $ +*/ #include "ca.h" void ezgcdp(vl,p1,p2,pr) @@ -6,23 +54,23 @@ VL vl; P p1,p2; P *pr; { - P f1,f2; - Q c; + P f1,f2; + Q c; - if ( !p1 ) - if ( !p2 ) - *pr = 0; - else - *pr = p2; - else if ( !p2 ) - *pr = p1; - else { - if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) ) - error("ezgcdp : invalid argument"); - ptozp(p1,1,&c,&f1); - ptozp(p2,1,&c,&f2); - ezgcdpz(vl,f1,f2,pr); - } + if ( !p1 ) + if ( !p2 ) + *pr = 0; + else + *pr = p2; + else if ( !p2 ) + *pr = p1; + else { + if ( !qpcheck((Obj)p1) || !qpcheck((Obj)p2) ) + error("ezgcdp : invalid argument"); + ptozp(p1,1,&c,&f1); + ptozp(p2,1,&c,&f2); + ezgcdpz(vl,f1,f2,pr); + } } /* @@ -34,137 +82,137 @@ void ezgcdpz(vl,p1,p2,pr) VL vl; P p1,p2,*pr; { - P f1,f2,t,s,g,gcd; - P fc1,fc2,cont; - VL tvl,svl; - V v1,v2; - DCP dc,dc1,dc2; - N cn; - Q cq; - Q c1,c2,c; - P g1,g2; - P mgcd; - extern int nez; - P pm[2]; + P f1,f2,t,s,g,gcd; + P fc1,fc2,cont; + VL tvl,svl; + V v1,v2; + DCP dc,dc1,dc2; + N cn; + Q cq; + Q c1,c2,c; + P g1,g2; + P mgcd; + extern int nez; + P pm[2]; - if ( nez ) { - pm[0] = p1; pm[1] = p2; nezgcdnpz(vl,pm,2,pr); return; - } - monomialfctr(vl,p1,&f1,&dc1); p1 = f1; - monomialfctr(vl,p2,&f2,&dc2); p2 = f2; - for ( mgcd = (P)ONE; dc1; dc1 = NEXT(dc1) ) - for ( v1 = VR(dc1->c), dc = dc2; dc; dc = NEXT(dc) ) - if ( v1 == VR(dc->c) ) { - pwrp(vl,dc->c,cmpq(dc1->d,dc->d)>=0?dc->d:dc1->d,&t); - mulp(vl,mgcd,t,&s); mgcd = s; break; - } + if ( nez ) { + pm[0] = p1; pm[1] = p2; nezgcdnpz(vl,pm,2,pr); return; + } + monomialfctr(vl,p1,&f1,&dc1); p1 = f1; + monomialfctr(vl,p2,&f2,&dc2); p2 = f2; + for ( mgcd = (P)ONE; dc1; dc1 = NEXT(dc1) ) + for ( v1 = VR(dc1->c), dc = dc2; dc; dc = NEXT(dc) ) + if ( v1 == VR(dc->c) ) { + pwrp(vl,dc->c,cmpq(dc1->d,dc->d)>=0?dc->d:dc1->d,&t); + mulp(vl,mgcd,t,&s); mgcd = s; break; + } - if ( NUM(p1) ) { - if ( NUM(p2) ) { - gcdn(NM((Q)p1),NM((Q)p2),&cn); NTOQ(cn,1,c); - } else { - ptozp(p2,1,&c2,&t); - gcdn(NM((Q)p1),NM(c2),&cn); NTOQ(cn,1,c); - } - mulp(vl,(P)c,mgcd,pr); - return; - } else if ( NUM(p2) ) { - ptozp(p1,1,&c1,&t); - gcdn(NM(c1),NM((Q)p2),&cn); NTOQ(cn,1,c); - mulp(vl,(P)c,mgcd,pr); - return; - } + if ( NUM(p1) ) { + if ( NUM(p2) ) { + gcdn(NM((Q)p1),NM((Q)p2),&cn); NTOQ(cn,1,c); + } else { + ptozp(p2,1,&c2,&t); + gcdn(NM((Q)p1),NM(c2),&cn); NTOQ(cn,1,c); + } + mulp(vl,(P)c,mgcd,pr); + return; + } else if ( NUM(p2) ) { + ptozp(p1,1,&c1,&t); + gcdn(NM(c1),NM((Q)p2),&cn); NTOQ(cn,1,c); + mulp(vl,(P)c,mgcd,pr); + return; + } - ptozp(p1,1,&c1,&g1); ptozp(p2,1,&c2,&g2); - gcdn(NM(c1),NM(c2),&cn); NTOQ(cn,1,cq); + ptozp(p1,1,&c1,&g1); ptozp(p2,1,&c2,&g2); + gcdn(NM(c1),NM(c2),&cn); NTOQ(cn,1,cq); - if ( ( v1 = VR(f1) ) != ( v2 = VR(f2) ) ) { - while ( v1 != VR(vl) && v2 != VR(vl) ) - vl = NEXT(vl); - if ( v1 == VR(vl) ) { - pcp(vl,f1,&g,&fc1); - ezgcdpz(vl,fc1,f2,&t); - } else { - pcp(vl,f2,&g,&fc2); - ezgcdpz(vl,fc2,f1,&t); - } - mulp(vl,t,(P)cq,&s); mulp(vl,s,mgcd,pr); - return; - } + if ( ( v1 = VR(f1) ) != ( v2 = VR(f2) ) ) { + while ( v1 != VR(vl) && v2 != VR(vl) ) + vl = NEXT(vl); + if ( v1 == VR(vl) ) { + pcp(vl,f1,&g,&fc1); + ezgcdpz(vl,fc1,f2,&t); + } else { + pcp(vl,f2,&g,&fc2); + ezgcdpz(vl,fc2,f1,&t); + } + mulp(vl,t,(P)cq,&s); mulp(vl,s,mgcd,pr); + return; + } - pcp(vl,f1,&g1,&fc1); pcp(vl,f2,&g2,&fc2); - ezgcdpz(vl,fc1,fc2,&cont); - clctv(vl,g1,&tvl); clctv(vl,g2,&svl); - if ( !NEXT(tvl) && !NEXT(svl) ) { - uezgcdpz(vl,g1,g2,&t); - mulp(vl,t,cont,&s); mulp(vl,s,(P)cq,&t); mulp(vl,t,mgcd,pr); - return; - } + pcp(vl,f1,&g1,&fc1); pcp(vl,f2,&g2,&fc2); + ezgcdpz(vl,fc1,fc2,&cont); + clctv(vl,g1,&tvl); clctv(vl,g2,&svl); + if ( !NEXT(tvl) && !NEXT(svl) ) { + uezgcdpz(vl,g1,g2,&t); + mulp(vl,t,cont,&s); mulp(vl,s,(P)cq,&t); mulp(vl,t,mgcd,pr); + return; + } - if ( homdeg(g1) > homdeg(g2) ) { - t = g1; g1 = g2; g2 = t; - } - sqfrp(vl,g1,&dc); - for ( g = g2, gcd = (P)ONE; dc && !NUM(g); dc = NEXT(dc) ) { - if ( NUM(COEF(dc)) ) - continue; - ezgcdpp(vl,dc,g,&t); - if ( NUM(t) ) - continue; - mulp(vl,gcd,t,&s); gcd = s; - divsp(vl,g,t,&s); g = s; - } - mulp(vl,gcd,cont,&t); mulp(vl,t,(P)cq,&s); mulp(vl,s,mgcd,pr); + if ( homdeg(g1) > homdeg(g2) ) { + t = g1; g1 = g2; g2 = t; + } + sqfrp(vl,g1,&dc); + for ( g = g2, gcd = (P)ONE; dc && !NUM(g); dc = NEXT(dc) ) { + if ( NUM(COEF(dc)) ) + continue; + ezgcdpp(vl,dc,g,&t); + if ( NUM(t) ) + continue; + mulp(vl,gcd,t,&s); gcd = s; + divsp(vl,g,t,&s); g = s; + } + mulp(vl,gcd,cont,&t); mulp(vl,t,(P)cq,&s); mulp(vl,s,mgcd,pr); } void uezgcdpz(vl,p1,p2,pr) VL vl; P p1,p2,*pr; { - P g1,g2,gcd,t,s,g; - DCP dc; - N cn; - Q c1,c2,cq; - P f1,f2; + P g1,g2,gcd,t,s,g; + DCP dc; + N cn; + Q c1,c2,cq; + P f1,f2; - if ( NUM(p1) ) { - if ( NUM(p2) ) { - gcdn(NM((Q)p1),NM((Q)p2),&cn); NTOQ(cn,1,cq),*pr = (P)cq; - return; - } else { - ptozp(p2,1,&c2,&t); - gcdn(NM((Q)p1),NM(c2),&cn); NTOQ(cn,1,cq),*pr = (P)cq; - return; - } - } else if ( NUM(p2) ) { - ptozp(p1,1,&c1,&t); - gcdn(NM(c1),NM((Q)p2),&cn); NTOQ(cn,1,cq),*pr = (P)cq; - return; - } + if ( NUM(p1) ) { + if ( NUM(p2) ) { + gcdn(NM((Q)p1),NM((Q)p2),&cn); NTOQ(cn,1,cq),*pr = (P)cq; + return; + } else { + ptozp(p2,1,&c2,&t); + gcdn(NM((Q)p1),NM(c2),&cn); NTOQ(cn,1,cq),*pr = (P)cq; + return; + } + } else if ( NUM(p2) ) { + ptozp(p1,1,&c1,&t); + gcdn(NM(c1),NM((Q)p2),&cn); NTOQ(cn,1,cq),*pr = (P)cq; + return; + } - ptozp(p1,1,&c1,&f1); ptozp(p2,1,&c2,&f2); - gcdn(NM(c1),NM(c2),&cn); NTOQ(cn,1,cq); - if ( UDEG(f1) > UDEG(f2) ) { - g1 = f2; g2 = f1; - } else { - g1 = f1; g2 = f2; - } - usqp(g1,&dc); - for ( g = g2, gcd = (P)ONE; dc && !NUM(g); dc = NEXT(dc) ) { - if ( NUM(COEF(dc)) ) - continue; - uezgcdpp(dc,g,&t); - if ( NUM(t) ) - continue; - mulp(CO,gcd,t,&s); gcd = s; - divsp(CO,g,t,&s); g = s; - } - mulp(vl,gcd,(P)cq,pr); + ptozp(p1,1,&c1,&f1); ptozp(p2,1,&c2,&f2); + gcdn(NM(c1),NM(c2),&cn); NTOQ(cn,1,cq); + if ( UDEG(f1) > UDEG(f2) ) { + g1 = f2; g2 = f1; + } else { + g1 = f1; g2 = f2; + } + usqp(g1,&dc); + for ( g = g2, gcd = (P)ONE; dc && !NUM(g); dc = NEXT(dc) ) { + if ( NUM(COEF(dc)) ) + continue; + uezgcdpp(dc,g,&t); + if ( NUM(t) ) + continue; + mulp(CO,gcd,t,&s); gcd = s; + divsp(CO,g,t,&s); g = s; + } + mulp(vl,gcd,(P)cq,pr); } /* - dc : dc pair c^d - r = gcd(c^d,f) + dc : dc pair c^d + r = gcd(c^d,f) */ void ezgcdpp(vl,dc,f,r) @@ -173,26 +221,26 @@ DCP dc; P f; P *r; { - int k; - P g,h,t,s,gcd; + int k; + P g,h,t,s,gcd; - if ( NUM(f) ) { - *r = (P)ONE; - return; - } - k = QTOS(DEG(dc)) - 1; -/* ezgcd1p(vl,COEF(dc),f,&gcd); */ - ezgcdnp(vl,COEF(dc),&f,1,&gcd); - g = gcd; divsp(vl,f,gcd,&h); - for ( ; k > 0; k-- ) { -/* ezgcd1p(vl,g,h,&t); */ - ezgcdnp(vl,g,&h,1,&t); - if ( NUM(t) ) - break; - mulp(vl,gcd,t,&s); gcd = s; - divsp(vl,h,t,&s); h = s; - } - *r = gcd; + if ( NUM(f) ) { + *r = (P)ONE; + return; + } + k = QTOS(DEG(dc)) - 1; +/* ezgcd1p(vl,COEF(dc),f,&gcd); */ + ezgcdnp(vl,COEF(dc),&f,1,&gcd); + g = gcd; divsp(vl,f,gcd,&h); + for ( ; k > 0; k-- ) { +/* ezgcd1p(vl,g,h,&t); */ + ezgcdnp(vl,g,&h,1,&t); + if ( NUM(t) ) + break; + mulp(vl,gcd,t,&s); gcd = s; + divsp(vl,h,t,&s); h = s; + } + *r = gcd; } void ezgcd1p(vl,p0,p1,gcdp) @@ -200,14 +248,14 @@ VL vl; P p0,p1; P *gcdp; { - VL nvl,tvl,vl0,vl1,vl2; - P f0,f1,t; + VL nvl,tvl,vl0,vl1,vl2; + P f0,f1,t; - clctv(vl,p0,&vl0); clctv(vl,p1,&vl1); mergev(vl,vl0,vl1,&vl2); - minchdegp(vl,p0,&tvl,&t); reordvar(vl2,tvl->v,&nvl); - reorderp(nvl,vl,p0,&f0); reorderp(nvl,vl,p1,&f1); - ezgcdnp(nvl,f0,&f1,1,&t); - reorderp(vl,nvl,t,gcdp); + clctv(vl,p0,&vl0); clctv(vl,p1,&vl1); mergev(vl,vl0,vl1,&vl2); + minchdegp(vl,p0,&tvl,&t); reordvar(vl2,tvl->v,&nvl); + reorderp(nvl,vl,p0,&f0); reorderp(nvl,vl,p1,&f1); + ezgcdnp(nvl,f0,&f1,1,&t); + reorderp(vl,nvl,t,gcdp); } void uezgcdpp(dc,f,r) @@ -215,25 +263,25 @@ DCP dc; P f; P *r; { - int k; - P g,h,t,s,gcd; + int k; + P g,h,t,s,gcd; - if ( NUM(f) ) { - *r = (P)ONE; - return; - } + if ( NUM(f) ) { + *r = (P)ONE; + return; + } - k = QTOS(DEG(dc)) - 1; - uezgcd1p(COEF(dc),f,&gcd); - g = gcd; divsp(CO,f,gcd,&h); - for ( ; k > 0; k-- ) { - uezgcd1p(g,h,&t); - if ( NUM(t) ) - break; - mulp(CO,gcd,t,&s); gcd = s; - divsp(CO,h,t,&s); h = s; - } - *r = gcd; + k = QTOS(DEG(dc)) - 1; + uezgcd1p(COEF(dc),f,&gcd); + g = gcd; divsp(CO,f,gcd,&h); + for ( ; k > 0; k-- ) { + uezgcd1p(g,h,&t); + if ( NUM(t) ) + break; + mulp(CO,gcd,t,&s); gcd = s; + divsp(CO,h,t,&s); h = s; + } + *r = gcd; } /* @@ -249,176 +297,176 @@ 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 *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 *tps; + VL nvl1,nvl2,nvl,tvl; + V v; + int i,j,k,d0,dg,dg0,dmin; + VN vn0,vn1,vnt; + int nv,flag; - /* 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); lg = lp0; dmin = d0 = deg(v,p0); - clctv(vl,p0,&nvl); - for ( i = 0; i < m; i++ ) { - clctv(vl,tps[i],&nvl1); mergev(vl,nvl,nvl1,&nvl2); - nvl = nvl2; - ezgcdpz(nvl,lg,LC(tps[i]),&t); lg = t; - } + lp0 = LC(p0); lg = lp0; dmin = d0 = deg(v,p0); + clctv(vl,p0,&nvl); + for ( i = 0; i < m; i++ ) { + clctv(vl,tps[i],&nvl1); mergev(vl,nvl,nvl1,&nvl2); + nvl = nvl2; + ezgcdpz(nvl,lg,LC(tps[i]),&t); lg = t; + } - 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,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; - /* 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 */ + 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; - /* 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 */ + 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; - /* 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; + } + } + } + } while ( !nextbin(vnt,j) ); } /* - f : sqfr + f : sqfr */ void uezgcd1p(f,g,r) P f,g; P *r; { - UM wf,wf1,wf2,wg,wg1,wgcd,wcof; - int d1,d2,dg,m,index,n; - ML list; - DCP dc; - P t; - Q c; + UM wf,wf1,wf2,wg,wg1,wgcd,wcof; + int d1,d2,dg,m,index,n; + ML list; + DCP dc; + P t; + Q c; - if ( NUM(f) || NUM(g) ) { - *r = (P)ONE; - return; - } - ptozp(f,1,&c,&t); f = t; ptozp(g,1,&c,&t); g = t; - d1 = UDEG(f); d2 = UDEG(g); - n = MAX(d1,d2)+1; - wf = W_UMALLOC(n); wf1 = W_UMALLOC(n); - wf2 = W_UMALLOC(n); wg = W_UMALLOC(n); - wg1 = W_UMALLOC(n); wgcd = W_UMALLOC(n); - wcof = W_UMALLOC(n); + if ( NUM(f) || NUM(g) ) { + *r = (P)ONE; + return; + } + ptozp(f,1,&c,&t); f = t; ptozp(g,1,&c,&t); g = t; + d1 = UDEG(f); d2 = UDEG(g); + n = MAX(d1,d2)+1; + wf = W_UMALLOC(n); wf1 = W_UMALLOC(n); + wf2 = W_UMALLOC(n); wg = W_UMALLOC(n); + wg1 = W_UMALLOC(n); wgcd = W_UMALLOC(n); + wcof = W_UMALLOC(n); - for ( index = INDEX, dg = MIN(d1,d2)+1; ; ) { - m = sprime[index++]; - if ( !rem(NM((Q)UCOEF(f)),m) || !rem(NM((Q)UCOEF(g)),m)) - continue; - ptoum(m,f,wf); cpyum(wf,wf1); - diffum(m,wf1,wf2); gcdum(m,wf1,wf2,wgcd); - if ( DEG(wgcd) > 0 ) - continue; - ptoum(m,g,wg); cpyum(wg,wg1); cpyum(wf,wf1); - gcdum(m,wf1,wg1,wgcd); - if ( dg <= DEG(wgcd) ) - continue; - else - dg = DEG(wgcd); - if ( dg == 0 ) { - *r = (P)ONE; - return; - } else if ( dg == d1 ) { - if ( divtpz(CO,g,f,&t) ) { - *r = f; - return; - } - } else if ( dg == d2 ) { - if ( divtpz(CO,f,g,&t) ) { - *r = g; - return; - } - } else { - divum(m,wf,wgcd,wcof); - ezgcdhensel(f,m,wcof,wgcd,&list); - dtest(f,list,1,&dc); - if ( NEXT(dc) ) { - if ( divtpz(CO,g,COEF(dc),&t) ) { - *r = COEF(dc); - return; - } - } - } - } + for ( index = INDEX, dg = MIN(d1,d2)+1; ; ) { + m = sprime[index++]; + if ( !rem(NM((Q)UCOEF(f)),m) || !rem(NM((Q)UCOEF(g)),m)) + continue; + ptoum(m,f,wf); cpyum(wf,wf1); + diffum(m,wf1,wf2); gcdum(m,wf1,wf2,wgcd); + if ( DEG(wgcd) > 0 ) + continue; + ptoum(m,g,wg); cpyum(wg,wg1); cpyum(wf,wf1); + gcdum(m,wf1,wg1,wgcd); + if ( dg <= DEG(wgcd) ) + continue; + else + dg = DEG(wgcd); + if ( dg == 0 ) { + *r = (P)ONE; + return; + } else if ( dg == d1 ) { + if ( divtpz(CO,g,f,&t) ) { + *r = f; + return; + } + } else if ( dg == d2 ) { + if ( divtpz(CO,f,g,&t) ) { + *r = g; + return; + } + } else { + divum(m,wf,wgcd,wcof); + ezgcdhensel(f,m,wcof,wgcd,&list); + dtest(f,list,1,&dc); + if ( NEXT(dc) ) { + if ( divtpz(CO,g,COEF(dc),&t) ) { + *r = COEF(dc); + return; + } + } + } + } } - + void uexgcdnp(vl,p1,ps,n,vn,cbd,g,h,a,b,q) VL vl; P p1; @@ -429,49 +477,49 @@ Q cbd; P *g,*h,*a,*b; Q *q; { - UM wg,wh,wg1,wh1,wgcd,wt; - P t,s,gcd,cof,gk,hk,ak,bk; - Q c,tq,bb,tq1,qk; - int mod,k,i,index,d; + UM wg,wh,wg1,wh1,wgcd,wt; + P t,s,gcd,cof,gk,hk,ak,bk; + Q c,tq,bb,tq1,qk; + int mod,k,i,index,d; - ptozp(p1,1,&c,&gcd); - for ( i = 0; i < n; i++ ) { - substvp(vl,ps[i],vn,&t); - uezgcd1p(gcd,t,&s); - if ( NUM(s) ) { - *g = (P)ONE; *h = p1; *a = 0; *b = (P)ONE; - return; - } else - gcd = s; - } + ptozp(p1,1,&c,&gcd); + for ( i = 0; i < n; i++ ) { + substvp(vl,ps[i],vn,&t); + uezgcd1p(gcd,t,&s); + if ( NUM(s) ) { + *g = (P)ONE; *h = p1; *a = 0; *b = (P)ONE; + return; + } else + gcd = s; + } - if ( !dcomp(p1,gcd) ) { - *g = p1; *h = (P)ONE; *a = 0; *b = (P)ONE; - return; - } - - divsp(CO,p1,gcd,&cof); - d = UDEG(p1)+1; - wg = W_UMALLOC(d); wh = W_UMALLOC(d); - wg1 = W_UMALLOC(d); wh1 = W_UMALLOC(d); - wgcd = W_UMALLOC(d); wt = W_UMALLOC(d); - for ( index = INDEX; ; ) { - mod = sprime[index++]; - if ( !rem(NM((Q)LC(p1)),mod) ) - continue; - ptoum(mod,gcd,wg); ptoum(mod,cof,wh); - cpyum(wg,wg1); cpyum(wh,wh1); - gcdum(mod,wg,wh,wt); cpyum(wt,wgcd); - if ( DEG(wgcd) > 0 ) - continue; - STOQ(mod,tq); addq(cbd,cbd,&bb); - for ( k = 1; cmpq(tq,bb) < 0; k++ ) { - mulq(tq,tq,&tq1); tq = tq1; - } - henzq(p1,gcd,wg1,cof,wh1,mod,k,&gk,&hk,&bk,&ak,&qk); - break; - } - *g = gk; *h = hk; *a = ak; *b = bk; *q = tq; + if ( !dcomp(p1,gcd) ) { + *g = p1; *h = (P)ONE; *a = 0; *b = (P)ONE; + return; + } + + divsp(CO,p1,gcd,&cof); + d = UDEG(p1)+1; + wg = W_UMALLOC(d); wh = W_UMALLOC(d); + wg1 = W_UMALLOC(d); wh1 = W_UMALLOC(d); + wgcd = W_UMALLOC(d); wt = W_UMALLOC(d); + for ( index = INDEX; ; ) { + mod = sprime[index++]; + if ( !rem(NM((Q)LC(p1)),mod) ) + continue; + ptoum(mod,gcd,wg); ptoum(mod,cof,wh); + cpyum(wg,wg1); cpyum(wh,wh1); + gcdum(mod,wg,wh,wt); cpyum(wt,wgcd); + if ( DEG(wgcd) > 0 ) + continue; + STOQ(mod,tq); addq(cbd,cbd,&bb); + for ( k = 1; cmpq(tq,bb) < 0; k++ ) { + mulq(tq,tq,&tq1); tq = tq1; + } + henzq(p1,gcd,wg1,cof,wh1,mod,k,&gk,&hk,&bk,&ak,&qk); + break; + } + *g = gk; *h = hk; *a = ak; *b = bk; *q = tq; } void ezgcdhensel(f,mod,cof,gcd,listp) @@ -480,45 +528,45 @@ int mod; UM gcd,cof; ML *listp; { - register int i,j; - int q,n,bound; - int *p; - int **pp; - ML blist,clist,bqlist,cqlist,rlist; - UM *b; - LUM fl,tl; - LUM *l; - LUM LUMALLOC(); - - blist = W_MLALLOC(2); - blist->n = 2; blist->mod = mod; - blist->c[0] = (pointer)cof; blist->c[1] = (pointer)gcd; - gcdgen(f,blist,&clist); - henprep(f,blist,clist,&bqlist,&cqlist); - n = bqlist->n; q = bqlist->mod; - bqlist->bound = cqlist->bound = bound = mignotte(q,f); + register int i,j; + int q,n,bound; + int *p; + int **pp; + ML blist,clist,bqlist,cqlist,rlist; + UM *b; + LUM fl,tl; + LUM *l; + LUM LUMALLOC(); + + blist = W_MLALLOC(2); + blist->n = 2; blist->mod = mod; + blist->c[0] = (pointer)cof; blist->c[1] = (pointer)gcd; + gcdgen(f,blist,&clist); + henprep(f,blist,clist,&bqlist,&cqlist); + n = bqlist->n; q = bqlist->mod; + bqlist->bound = cqlist->bound = bound = mignotte(q,f); - if ( bound == 1 ) { - *listp = rlist = MLALLOC(n); - rlist->n = n; - rlist->mod = q; - rlist->bound = bound; + if ( bound == 1 ) { + *listp = rlist = MLALLOC(n); + rlist->n = n; + rlist->mod = q; + rlist->bound = bound; - for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; - i < n; i++ ) { - tl = LUMALLOC(DEG(b[i]),1); - l[i] = tl; - p = COEF(b[i]); + for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; + i < n; i++ ) { + tl = LUMALLOC(DEG(b[i]),1); + l[i] = tl; + p = COEF(b[i]); - for ( j = 0, pp = COEF(tl); - j <= DEG(tl); j++ ) - pp[j][0] = p[j]; - } - } else { - W_LUMALLOC((int)UDEG(f),bound,fl); - ptolum(q,bound,f,fl); - henmain(fl,bqlist,cqlist,listp); - } + for ( j = 0, pp = COEF(tl); + j <= DEG(tl); j++ ) + pp[j][0] = p[j]; + } + } else { + W_LUMALLOC((int)UDEG(f),bound,fl); + ptolum(q,bound,f,fl); + henmain(fl,bqlist,cqlist,listp); + } } void ezgcdnpz(vl,ps,m,pr) @@ -526,60 +574,60 @@ VL vl; P *ps,*pr; int m; { - P t,s,gcd; - P cont; - VL tvl,svl,nvl; - int i; - DCP dc; - N cn; - Q cq; - P *pl,*ppl,*pcl; - Q *cl; - V v; + P t,s,gcd; + P cont; + VL tvl,svl,nvl; + int i; + DCP dc; + N cn; + Q cq; + P *pl,*ppl,*pcl; + Q *cl; + V v; - pl = (P *)ALLOCA((m+1)*sizeof(P)); - cl = (Q *)ALLOCA((m+1)*sizeof(Q)); - for ( i = 0; i < m; i++ ) - ptozp(ps[i],1,&cl[i],&pl[i]); - for ( i = 1, cq = cl[0]; i < m; i++ ) { - gcdn(NM(cl[i]),NM(cq),&cn); NTOQ(cn,1,cq); - } - for ( i = 0; i < m; i++ ) - if ( NUM(pl[i]) ) { - *pr = (P)cq; - return; - } + pl = (P *)ALLOCA((m+1)*sizeof(P)); + cl = (Q *)ALLOCA((m+1)*sizeof(Q)); + for ( i = 0; i < m; i++ ) + ptozp(ps[i],1,&cl[i],&pl[i]); + for ( i = 1, cq = cl[0]; i < m; i++ ) { + gcdn(NM(cl[i]),NM(cq),&cn); NTOQ(cn,1,cq); + } + for ( i = 0; i < m; i++ ) + if ( NUM(pl[i]) ) { + *pr = (P)cq; + return; + } - for ( i = 0, nvl = 0; i < m; i++ ) { - clctv(vl,ps[i],&tvl); mergev(vl,nvl,tvl,&svl); nvl = svl; - } + for ( i = 0, nvl = 0; i < m; i++ ) { + clctv(vl,ps[i],&tvl); mergev(vl,nvl,tvl,&svl); nvl = svl; + } - ppl = (P *)ALLOCA((m+1)*sizeof(P)); - pcl = (P *)ALLOCA((m+1)*sizeof(P)); - for ( i = 0, v = nvl->v; i < m; i++ ) - if ( v == VR(pl[i]) ) - pcp(nvl,pl[i],&ppl[i],&pcl[i]); - else { - ppl[i] = (P)ONE; pcl[i] = pl[i]; - } - ezgcdnpz(nvl,pcl,m,&cont); + ppl = (P *)ALLOCA((m+1)*sizeof(P)); + pcl = (P *)ALLOCA((m+1)*sizeof(P)); + for ( i = 0, v = nvl->v; i < m; i++ ) + if ( v == VR(pl[i]) ) + pcp(nvl,pl[i],&ppl[i],&pcl[i]); + else { + ppl[i] = (P)ONE; pcl[i] = pl[i]; + } + ezgcdnpz(nvl,pcl,m,&cont); - for ( i = 0; i < m; i++ ) - if ( NUM(ppl[i]) ) { - mulp(nvl,cont,(P)cq,pr); - return; - } - sortplist(ppl,m); - sqfrp(nvl,ppl[0],&dc); - for ( dc = NEXT(dc), gcd = (P)ONE; dc; dc = NEXT(dc) ) { - if ( NUM(COEF(dc)) ) - continue; - ezgcdnpp(vl,dc,ppl+1,m-1,&t); - if ( NUM(t) ) - continue; - mulp(vl,gcd,t,&s); gcd = s; - } - mulp(vl,gcd,cont,&t); mulp(vl,t,(P)cq,pr); + for ( i = 0; i < m; i++ ) + if ( NUM(ppl[i]) ) { + mulp(nvl,cont,(P)cq,pr); + return; + } + sortplist(ppl,m); + sqfrp(nvl,ppl[0],&dc); + for ( dc = NEXT(dc), gcd = (P)ONE; dc; dc = NEXT(dc) ) { + if ( NUM(COEF(dc)) ) + continue; + ezgcdnpp(vl,dc,ppl+1,m-1,&t); + if ( NUM(t) ) + continue; + mulp(vl,gcd,t,&s); gcd = s; + } + mulp(vl,gcd,cont,&t); mulp(vl,t,(P)cq,pr); } void ezgcdnpp(vl,dc,pl,m,r) @@ -589,38 +637,38 @@ 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; - ezgcdnp(vl,COEF(dc),pl,m,&gcd); - if ( NUM(gcd) ) { - *r = (P)ONE; - return; - } - pm = (P *) ALLOCA((m+1)*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-- ) { - ezgcdnp(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; + ezgcdnp(vl,COEF(dc),pl,m,&gcd); + if ( NUM(gcd) ) { + *r = (P)ONE; + return; + } + pm = (P *) ALLOCA((m+1)*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-- ) { + ezgcdnp(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; } void pcp(vl,p,pp,c) @@ -628,32 +676,32 @@ VL vl; P p; P *pp,*c; { - P f,g,n; - DCP dc; - P *pl; - int i,m; - extern int nez; + P f,g,n; + DCP dc; + P *pl; + int i,m; + extern int nez; - if ( NUM(p) ) { - *c = p; - *pp = (P)ONE; - } else { - for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ ); - if ( m == 1 ) { - *c = COEF(DC(p)); - divsp(vl,p,*c,pp); - return; - } - ptozp(p,1,(Q *)&n,&f); - pl = (P *)ALLOCA((m+1)*sizeof(P)); - for ( i = 0, dc = DC(f); i < m; dc = NEXT(dc), i++ ) - pl[i] = COEF(dc); - if ( nez ) - nezgcdnpz(vl,pl,m,&g); - else - ezgcdnpz(vl,pl,m,&g); - mulp(vl,g,n,c); divsp(vl,f,g,pp); - } + if ( NUM(p) ) { + *c = p; + *pp = (P)ONE; + } else { + for ( m = 0, dc = DC(p); dc; dc = NEXT(dc), m++ ); + if ( m == 1 ) { + *c = COEF(DC(p)); + divsp(vl,p,*c,pp); + return; + } + ptozp(p,1,(Q *)&n,&f); + pl = (P *)ALLOCA((m+1)*sizeof(P)); + for ( i = 0, dc = DC(f); i < m; dc = NEXT(dc), i++ ) + pl[i] = COEF(dc); + if ( nez ) + nezgcdnpz(vl,pl,m,&g); + else + ezgcdnpz(vl,pl,m,&g); + mulp(vl,g,n,c); divsp(vl,f,g,pp); + } } #if 0 @@ -662,20 +710,20 @@ VL vl; P p; P *pp,*c; { - P f,g,n,t; - DCP dc; + P f,g,n,t; + DCP dc; - if ( NUM(p) ) { - *c = p; - *pp = (P)ONE; - } else { - ptozp(p,1,&n,&f); - for ( g = LC(f), dc = NEXT(DC(f)); dc; dc = NEXT(dc) ) { - ezgcdpz(vl,g,COEF(dc),&t); g = t; - } - mulp(vl,g,n,c); - divsp(vl,f,g,pp); - } + if ( NUM(p) ) { + *c = p; + *pp = (P)ONE; + } else { + ptozp(p,1,&n,&f); + for ( g = LC(f), dc = NEXT(DC(f)); dc; dc = NEXT(dc) ) { + ezgcdpz(vl,g,COEF(dc),&t); g = t; + } + mulp(vl,g,n,c); + divsp(vl,f,g,pp); + } } #endif @@ -684,33 +732,33 @@ VL vl; int m; P *ps,l,p; { - int i; - P r,t; + int i; + P r,t; - for ( i = 0; i < m; i++ ) { - mulp(vl,ps[i],l,&t); - if ( !divtpz(vl,t,p,&r) ) - return ( 0 ); - } - return ( 1 ); + for ( i = 0; i < m; i++ ) { + mulp(vl,ps[i],l,&t); + if ( !divtpz(vl,t,p,&r) ) + return ( 0 ); + } + return ( 1 ); } void sortplist(plist,n) P *plist; int n; { - int i,j; - P t; + int i,j; + P t; - for ( i = 0; i < n; i++ ) - for ( j = i + 1; j < n; j++ ) - if ( cmpq(DEG(DC(plist[i])),DEG(DC(plist[j]))) > 0 ) { - t = plist[i]; plist[i] = plist[j]; plist[j] = t; - } + for ( i = 0; i < n; i++ ) + for ( j = i + 1; j < n; j++ ) + if ( cmpq(DEG(DC(plist[i])),DEG(DC(plist[j]))) > 0 ) { + t = plist[i]; plist[i] = plist[j]; plist[j] = t; + } } int dcomp(p1,p2) P p1,p2; { - return ( cmpq(DEG(DC(p1)),DEG(DC(p2))) ); + return ( cmpq(DEG(DC(p1)),DEG(DC(p2))) ); }