=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/A.c,v retrieving revision 1.2 retrieving revision 1.8 diff -u -p -r1.2 -r1.8 --- OpenXM_contrib2/asir2000/engine/A.c 2000/08/21 08:31:24 1.2 +++ OpenXM_contrib2/asir2000/engine/A.c 2018/03/29 01:32:51 1.8 @@ -23,7 +23,7 @@ * 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@flab.fujitsu.co.jp of the detailed specification + * 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. * @@ -45,83 +45,86 @@ * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. * - * $OpenXM: OpenXM_contrib2/asir2000/engine/A.c,v 1.1.1.1 1999/12/03 07:39:07 noro Exp $ + * $OpenXM: OpenXM_contrib2/asir2000/engine/A.c,v 1.7 2013/06/13 18:40:31 ohara Exp $ */ #include "ca.h" +#include "parse.h" +int get_lprime(); + void pdiva(vl,p0,p1,p2,pr) VL vl; P p1,p2,p0; P *pr; { - Q m,d; - P t,q,r; + Q m,d; + P t,q,r; - if ( VR(p2) == VR(p0) ) { - *pr = p1; - return; - } - if ( VR(p1) == VR(p0) ) { - error("pdiva : ???"); - return; - } - subq(DEG(DC(p1)),DEG(DC(p2)),&m); addq(m,ONE,&d); - pwrq((Q)COEF(DC(p2)),d,&m); mulp(vl,p1,(P)m,&t); - divsrp(vl,t,p2,&q,&r); remsdcp(vl,q,p0,pr); + if ( VR(p2) == VR(p0) ) { + *pr = p1; + return; + } + if ( VR(p1) == VR(p0) ) { + error("pdiva : ???"); + return; + } + subq(DEG(DC(p1)),DEG(DC(p2)),&m); addq(m,ONE,&d); + pwrq((Q)COEF(DC(p2)),d,&m); mulp(vl,p1,(P)m,&t); + divsrp(vl,t,p2,&q,&r); remsdcp(vl,q,p0,pr); } void rema(vl,p0,p1,p2,pr) VL vl; P p1,p2,p0,*pr; { - P m,m1,m2,g,lc,x,t; - Q q; - V v; - int d1,d2; + P m,m1,m2,g,lc,x,t; + Q q; + V v; + int d1,d2; - if ( !p1 || !p2 || NUM(p2) || ( VR(p2) == VR(p0) ) ) - *pr = 0; - else if ( NUM(p1) ) - *pr = p1; - else if ( ( v = VR(p1) ) == VR(p0) ) - remsdcp(vl,p1,p0,pr); - else if ( ( d1 = deg(v,p1) ) < ( d2 = deg(v,p2) ) ) - *pr = p1; - else if ( d1 == d2 ) { - mulp(vl,p1,LC(p2),&m); mulp(vl,p2,LC(p1),&m1); - subp(vl,m,m1,&m2); - remsdcp(vl,m2,p0,pr); - } else { - g = p1; lc = LC(p2); MKV(v,x); - while ( ( d1 = deg(v,g) ) >= d2 ) { - mulp(vl,g,lc,&m); mulp(vl,p2,LC(g),&m1); - STOQ(d1-d2,q); pwrp(vl,x,q,&t); - mulp(vl,m1,t,&m2); subp(vl,m,m2,&m1); - remsdcp(vl,m1,p0,&g); - } - *pr = g; - } + if ( !p1 || !p2 || NUM(p2) || ( VR(p2) == VR(p0) ) ) + *pr = 0; + else if ( NUM(p1) ) + *pr = p1; + else if ( ( v = VR(p1) ) == VR(p0) ) + remsdcp(vl,p1,p0,pr); + else if ( ( d1 = deg(v,p1) ) < ( d2 = deg(v,p2) ) ) + *pr = p1; + else if ( d1 == d2 ) { + mulp(vl,p1,LC(p2),&m); mulp(vl,p2,LC(p1),&m1); + subp(vl,m,m1,&m2); + remsdcp(vl,m2,p0,pr); + } else { + g = p1; lc = LC(p2); MKV(v,x); + while ( ( d1 = deg(v,g) ) >= d2 ) { + mulp(vl,g,lc,&m); mulp(vl,p2,LC(g),&m1); + STOQ(d1-d2,q); pwrp(vl,x,q,&t); + mulp(vl,m1,t,&m2); subp(vl,m,m2,&m1); + remsdcp(vl,m1,p0,&g); + } + *pr = g; + } } void gcda(vl,p0,p1,p2,pr) VL vl; P p1,p2,p0,*pr; { - P g1,g2,r,t,c; + P g1,g2,r,t,c; - remsdcp(vl,p1,p0,&t); pcp(vl,t,&g1,&c); remsdcp(vl,p2,p0,&g2); - while ( 1 ) { - if ( NUM(g2) || (VR(g2) == VR(p0)) ) { - *pr = (P)ONE; - return; - } - pcp(vl,g2,&t,&c); pmonic(vl,p0,t,&g2); rema(vl,p0,g1,g2,&r); - if ( !r ) { - pmonic(vl,p0,g2,pr); - return; - } - g1 = g2; g2 = r; - } + remsdcp(vl,p1,p0,&t); pcp(vl,t,&g1,&c); remsdcp(vl,p2,p0,&g2); + while ( 1 ) { + if ( NUM(g2) || (VR(g2) == VR(p0)) ) { + *pr = (P)ONE; + return; + } + pcp(vl,g2,&t,&c); pmonic(vl,p0,t,&g2); rema(vl,p0,g1,g2,&r); + if ( !r ) { + pmonic(vl,p0,g2,pr); + return; + } + g1 = g2; g2 = r; + } } void sqa(vl,p0,p,dcp) @@ -129,52 +132,52 @@ VL vl; P p,p0; DCP *dcp; { - V v; - int i,j,n; - P f,f1,f2,gcd,flat,q,r,t; - P *a; - int *b; - DCP dc,dc0; + V v; + int i,j,n; + P f,f1,f2,gcd,flat,q,r,t; + P *a; + int *b; + DCP dc,dc0; - if ( ( v = VR(p) ) == VR(p0) ) { - NEWDC(dc); *dcp = dc; COEF(dc) = (P)ONE; DEG(dc) = ONE; - NEXT(dc) = 0; - } else { - n = deg(v,p); W_CALLOC(n,P,a); W_CALLOC(n,int,b); - for ( i = 0, f = p; ;i++ ) { - if ( VR(f) != v ) - break; - remsdcp(vl,f,p0,&f1); diffp(vl,f1,VR(f1),&f2); - while ( 1 ) { - pmonic(vl,p0,f2,&t); f2 = t; rema(vl,p0,f1,f2,&r); - if ( !r ) { - pmonic(vl,p0,f2,&gcd); pqra(vl,p0,f,gcd,&flat,&r); - break; - } - f1 = f2; f2 = r; - } - if ( VR(gcd) != v ) { - a[i] = f; b[i] = 1; - break; - } - for ( j = 1, f = gcd; ; j++ ) { - pqra(vl,p0,f,flat,&q,&r); - if ( r ) - break; - else - f = q; - } - a[i] = flat; b[i] = j; - } - for ( i = 0, j = 0, dc0 = 0; a[i]; i++ ) { - NEXTDC(dc0,dc); j += b[i]; STOQ(j,DEG(dc)); - if ( a[i+1] ) - pqra(vl,p0,a[i],a[i+1],&COEF(dc),&t); - else - COEF(dc) = a[i]; - } - NEXT(dc) = 0; *dcp = dc0; - } + if ( ( v = VR(p) ) == VR(p0) ) { + NEWDC(dc); *dcp = dc; COEF(dc) = (P)ONE; DEG(dc) = ONE; + NEXT(dc) = 0; + } else { + n = deg(v,p); W_CALLOC(n,P,a); W_CALLOC(n,int,b); + for ( i = 0, f = p; ;i++ ) { + if ( VR(f) != v ) + break; + remsdcp(vl,f,p0,&f1); diffp(vl,f1,VR(f1),&f2); + while ( 1 ) { + pmonic(vl,p0,f2,&t); f2 = t; rema(vl,p0,f1,f2,&r); + if ( !r ) { + pmonic(vl,p0,f2,&gcd); pqra(vl,p0,f,gcd,&flat,&r); + break; + } + f1 = f2; f2 = r; + } + if ( VR(gcd) != v ) { + a[i] = f; b[i] = 1; + break; + } + for ( j = 1, f = gcd; ; j++ ) { + pqra(vl,p0,f,flat,&q,&r); + if ( r ) + break; + else + f = q; + } + a[i] = flat; b[i] = j; + } + for ( i = 0, j = 0, dc0 = 0; a[i]; i++ ) { + NEXTDC(dc0,dc); j += b[i]; STOQ(j,DEG(dc)); + if ( a[i+1] ) + pqra(vl,p0,a[i],a[i+1],&COEF(dc),&t); + else + COEF(dc) = a[i]; + } + NEXT(dc) = 0; *dcp = dc0; + } } void pqra(vl, p0, p1, p2, p, pr) @@ -182,120 +185,179 @@ VL vl; P p0,p1,p2; P *p,*pr; { - V v,v1; - P m,m1,q,r; - Q tq; - int n; + V v,v1; + P m,m1,q,r; + Q tq; + int n; - if ( ( v = VR(p1) ) != ( v1 = VR(p2) ) ) { - while ( ( VR(vl) != v ) && ( VR(vl) != v1 ) ) - vl = NEXT(vl); - if ( VR(vl) == v ) { - *p = p1; *pr = 0; - } else { - *p = 0; *pr = p1; - } - } else if ( (n = deg(v,p1) - deg(v,p2)) < 0 ) { - *pr = p1; *p = 0; - } else { - n++; STOQ(n,tq); - pwrp(vl,LC(p2),tq,&m); mulp(vl,m,p1,&m1); divsrp(vl,m1,p2,&q,&r); - remsdcp(vl,q,p0,p); remsdcp(vl,r,p0,pr); - } + if ( ( v = VR(p1) ) != ( v1 = VR(p2) ) ) { + while ( ( VR(vl) != v ) && ( VR(vl) != v1 ) ) + vl = NEXT(vl); + if ( VR(vl) == v ) { + *p = p1; *pr = 0; + } else { + *p = 0; *pr = p1; + } + } else if ( (n = deg(v,p1) - deg(v,p2)) < 0 ) { + *pr = p1; *p = 0; + } else { + n++; STOQ(n,tq); + pwrp(vl,LC(p2),tq,&m); mulp(vl,m,p1,&m1); divsrp(vl,m1,p2,&q,&r); + remsdcp(vl,q,p0,p); remsdcp(vl,r,p0,pr); + } } void pmonic(vl, p0, p, pr) VL vl; P p,p0,*pr; { - P d,m,pp,c,t; + P d,m,pp,c,t; - if ( NUM(p) || ( VR(p) == VR(p0) ) ) - *pr = (P)ONE; - else if ( NUM(COEF(DC(p))) ) - *pr = p; - else { - ptozp(COEF(DC(p)),1,(Q *)&c,&pp); pinva(p0,pp,&d); - mulp(vl,p,d,&m); remsdcp(vl,m,p0,&t); ptozp(t,1,(Q *)&c,pr); - } + if ( NUM(p) || ( VR(p) == VR(p0) ) ) + *pr = (P)ONE; + else if ( NUM(COEF(DC(p))) ) + *pr = p; + else { + ptozp(COEF(DC(p)),1,(Q *)&c,&pp); pinva(p0,pp,&d); + mulp(vl,p,d,&m); remsdcp(vl,m,p0,&t); ptozp(t,1,(Q *)&c,pr); + } } void remsdcp(vl, p, p0, pr) VL vl; P p,p0,*pr; { - P q,r,c; - DCP dc,dcr,dcr0; + P q,r,c; + DCP dc,dcr,dcr0; - if ( !p ) - *pr = 0; - else if ( NUM(p) ) - *pr = (P)ONE; - else if ( VR(p) == VR(p0) ) { - divsrp(vl,p,p0,&q,&r); ptozp(r,1,(Q *)&c,pr); - } else { - for ( dcr = dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) { - divsrp(vl,COEF(dc),p0,&q,&r); - if ( r ) { - NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = r; - } - } - if ( !dcr0 ) - *pr = 0; - else { - NEXT(dcr) = 0; MKP(VR(p),dcr0,r); ptozp(r,1,(Q *)&c,pr); - } - } + if ( !p ) + *pr = 0; + else if ( NUM(p) ) + *pr = (P)ONE; + else if ( VR(p) == VR(p0) ) { + divsrp(vl,p,p0,&q,&r); ptozp(r,1,(Q *)&c,pr); + } else { + for ( dcr = dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) { + divsrp(vl,COEF(dc),p0,&q,&r); + if ( r ) { + NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = r; + } + } + if ( !dcr0 ) + *pr = 0; + else { + NEXT(dcr) = 0; MKP(VR(p),dcr0,r); ptozp(r,1,(Q *)&c,pr); + } + } } void pinva(p0, p, pr) P p,p0,*pr; { - V v; - Q f,q,q1; - Q u,t,a,b,s; - P c,c1; - P tg,th,tr; - P z,zz,zzz; - UM wg,wh,ws,wt,wm; - int n,m,modres,mod,index,i; + V v; + Q f,q,q1; + Q u,t,a,b,s; + P c,c1; + P tg,th,tr; + P z,zz,zzz; + UM wg,wh,ws,wt,wm; + int n,m,modres,mod,index,i; - v = VR(p); n=deg(v,p); m=deg(v,p0); - norm(p,&a); norm(p0,&b); - STOQ(m,u); pwrq(a,u,&t); STOQ(n,u); pwrq(b,u,&s); - mulq(t,s,&u); - factorial(n+m,&t); mulq(u,t,&s); addq(s,s,&f); - wg = W_UMALLOC(m+n); wh = W_UMALLOC(m+n); - wt = W_UMALLOC(m+n); ws = W_UMALLOC(m+n); - wm = W_UMALLOC(m+n); - for ( q = ONE, t = 0, c = 0, index = 0; ; ) { - mod = lprime[index++]; - if ( !mod ) - error("pinva : lprime[] exhausted."); - if ( !rem(NM((Q)LC(p)),mod) || !rem(NM((Q)LC(p0)),mod) ) - continue; - ptomp(mod,p,&tg); ptomp(mod,p0,&th); srchump(mod,tg,th,&tr); - if ( !tr ) - continue; - else - modres = CONT((MQ)tr); - mptoum(tg,wg); mptoum(th,wh); - eucum(mod,wg,wh,ws,wt); /* wg * ws + wh * wt = 1 */ - DEG(wm) = 0; COEF(wm)[0] = modres; - mulum(mod,ws,wm,wt); - for ( i = DEG(wt); i >= 0; i-- ) - if ( ( COEF(wt)[i] * 2 ) > mod ) - COEF(wt)[i] -= mod; - chnrem(mod,v,c,q,wt,&c1,&q1); - if ( !ucmpp(c,c1) ) { - mulp(CO,c,p,&z); divsrp(CO,z,p0,&zz,&zzz); - if ( NUM(zzz) ) { - q = q1; c = c1; break; - } - } - q = q1; c = c1; - if ( cmpq(f,q) < 0 ) - break; - } - ptozp(c,1,&s,pr); + v = VR(p); n=deg(v,p); m=deg(v,p0); + norm(p,&a); norm(p0,&b); + STOQ(m,u); pwrq(a,u,&t); STOQ(n,u); pwrq(b,u,&s); + mulq(t,s,&u); + factorial(n+m,&t); mulq(u,t,&s); addq(s,s,&f); + wg = W_UMALLOC(m+n); wh = W_UMALLOC(m+n); + wt = W_UMALLOC(m+n); ws = W_UMALLOC(m+n); + wm = W_UMALLOC(m+n); + for ( q = ONE, t = 0, c = 0, index = 0; ; ) { + mod = get_lprime(index++); + if ( !rem(NM((Q)LC(p)),mod) || !rem(NM((Q)LC(p0)),mod) ) + continue; + ptomp(mod,p,&tg); ptomp(mod,p0,&th); srchump(mod,tg,th,&tr); + if ( !tr ) + continue; + else + modres = CONT((MQ)tr); + mptoum(tg,wg); mptoum(th,wh); + eucum(mod,wg,wh,ws,wt); /* wg * ws + wh * wt = 1 */ + DEG(wm) = 0; COEF(wm)[0] = modres; + mulum(mod,ws,wm,wt); + for ( i = DEG(wt); i >= 0; i-- ) + if ( ( COEF(wt)[i] * 2 ) > mod ) + COEF(wt)[i] -= mod; + chnrem(mod,v,c,q,wt,&c1,&q1); + if ( !ucmpp(c,c1) ) { + mulp(CO,c,p,&z); divsrp(CO,z,p0,&zz,&zzz); + if ( NUM(zzz) ) { + q = q1; c = c1; break; + } + } + q = q1; c = c1; + if ( cmpq(f,q) < 0 ) + break; + } + ptozp(c,1,&s,pr); +} + +/* + * compute g s.t. h*p0+g*p = res(v,p0,p) + */ + +void inva_chrem(P p0, P p, NODE *pr) +{ + V v; + Q f,q,q1; + Q u,t,a,b,s,dn; + P c,c0,c1; + P tg,th,tr; + P z,zz,zzz; + UM wg,wh,ws,wt,wm; + int n,m,modres,mod,index,i; + + v = VR(p); n=deg(v,p); m=deg(v,p0); + norm(p,&a); norm(p0,&b); + STOQ(m,u); pwrq(a,u,&t); STOQ(n,u); pwrq(b,u,&s); + mulq(t,s,&u); + factorial(n+m,&t); mulq(u,t,&s); addq(s,s,&f); + wg = W_UMALLOC(m+n); wh = W_UMALLOC(m+n); + wt = W_UMALLOC(m+n); ws = W_UMALLOC(m+n); + wm = W_UMALLOC(m+n); + for ( q = ONE, t = 0, c = 0, index = 0; ; ) { + mod = get_lprime(index++); + if ( !rem(NM((Q)LC(p)),mod) || !rem(NM((Q)LC(p0)),mod) ) + continue; + ptomp(mod,p,&tg); ptomp(mod,p0,&th); srchump(mod,tg,th,&tr); + if ( !tr ) + continue; + else + modres = CONT((MQ)tr); + mptoum(tg,wg); mptoum(th,wh); + eucum(mod,wg,wh,ws,wt); /* wg * ws + wh * wt = 1 */ + DEG(wm) = 0; COEF(wm)[0] = modres; + mulum(mod,ws,wm,wt); + for ( i = DEG(wt); i >= 0; i-- ) + if ( ( COEF(wt)[i] * 2 ) > mod ) + COEF(wt)[i] -= mod; + chnrem(mod,v,c,q,wt,&c1,&q1); + if ( !ucmpp(c,c1) ) { + mulp(CO,c,p,&z); divsrp(CO,z,p0,&zz,&zzz); + if ( NUM(zzz) ) { + q = (Q)zzz; break; + } + } + q = q1; c = c1; + if ( cmpq(f,q) < 0 ) { + mulp(CO,c,p,&z); divsrp(CO,z,p0,&zz,&zzz); + q = (Q)zzz; + break; + } + } + /* c*p = q mod p0 */ + /* c = s*c0 */ + ptozp(c,1,&s,&c0); + /* c/q = c0/(q/s) */ + divq(q,s,&dn); + *pr = (NODE)mknode(2,c0,dn); }