![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to A.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine|
Return to A.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 ! noro 1: /* $OpenXM: OpenXM/src/asir99/engine/A.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
! 2: #include "ca.h"
! 3:
! 4: void pdiva(vl,p0,p1,p2,pr)
! 5: VL vl;
! 6: P p1,p2,p0;
! 7: P *pr;
! 8: {
! 9: Q m,d;
! 10: P t,q,r;
! 11:
! 12: if ( VR(p2) == VR(p0) ) {
! 13: *pr = p1;
! 14: return;
! 15: }
! 16: if ( VR(p1) == VR(p0) ) {
! 17: error("pdiva : ???");
! 18: return;
! 19: }
! 20: subq(DEG(DC(p1)),DEG(DC(p2)),&m); addq(m,ONE,&d);
! 21: pwrq((Q)COEF(DC(p2)),d,&m); mulp(vl,p1,(P)m,&t);
! 22: divsrp(vl,t,p2,&q,&r); remsdcp(vl,q,p0,pr);
! 23: }
! 24:
! 25: void rema(vl,p0,p1,p2,pr)
! 26: VL vl;
! 27: P p1,p2,p0,*pr;
! 28: {
! 29: P m,m1,m2,g,lc,x,t;
! 30: Q q;
! 31: V v;
! 32: int d1,d2;
! 33:
! 34: if ( !p1 || !p2 || NUM(p2) || ( VR(p2) == VR(p0) ) )
! 35: *pr = 0;
! 36: else if ( NUM(p1) )
! 37: *pr = p1;
! 38: else if ( ( v = VR(p1) ) == VR(p0) )
! 39: remsdcp(vl,p1,p0,pr);
! 40: else if ( ( d1 = deg(v,p1) ) < ( d2 = deg(v,p2) ) )
! 41: *pr = p1;
! 42: else if ( d1 == d2 ) {
! 43: mulp(vl,p1,LC(p2),&m); mulp(vl,p2,LC(p1),&m1);
! 44: subp(vl,m,m1,&m2);
! 45: remsdcp(vl,m2,p0,pr);
! 46: } else {
! 47: g = p1; lc = LC(p2); MKV(v,x);
! 48: while ( ( d1 = deg(v,g) ) >= d2 ) {
! 49: mulp(vl,g,lc,&m); mulp(vl,p2,LC(g),&m1);
! 50: STOQ(d1-d2,q); pwrp(vl,x,q,&t);
! 51: mulp(vl,m1,t,&m2); subp(vl,m,m2,&m1);
! 52: remsdcp(vl,m1,p0,&g);
! 53: }
! 54: *pr = g;
! 55: }
! 56: }
! 57:
! 58: void gcda(vl,p0,p1,p2,pr)
! 59: VL vl;
! 60: P p1,p2,p0,*pr;
! 61: {
! 62: P g1,g2,r,t,c;
! 63:
! 64: remsdcp(vl,p1,p0,&t); pcp(vl,t,&g1,&c); remsdcp(vl,p2,p0,&g2);
! 65: while ( 1 ) {
! 66: if ( NUM(g2) || (VR(g2) == VR(p0)) ) {
! 67: *pr = (P)ONE;
! 68: return;
! 69: }
! 70: pcp(vl,g2,&t,&c); pmonic(vl,p0,t,&g2); rema(vl,p0,g1,g2,&r);
! 71: if ( !r ) {
! 72: pmonic(vl,p0,g2,pr);
! 73: return;
! 74: }
! 75: g1 = g2; g2 = r;
! 76: }
! 77: }
! 78:
! 79: void sqa(vl,p0,p,dcp)
! 80: VL vl;
! 81: P p,p0;
! 82: DCP *dcp;
! 83: {
! 84: V v;
! 85: int i,j,n;
! 86: P f,f1,f2,gcd,flat,q,r,t;
! 87: P *a;
! 88: int *b;
! 89: DCP dc,dc0;
! 90:
! 91: if ( ( v = VR(p) ) == VR(p0) ) {
! 92: NEWDC(dc); *dcp = dc; COEF(dc) = (P)ONE; DEG(dc) = ONE;
! 93: NEXT(dc) = 0;
! 94: } else {
! 95: n = deg(v,p); W_CALLOC(n,P,a); W_CALLOC(n,int,b);
! 96: for ( i = 0, f = p; ;i++ ) {
! 97: if ( VR(f) != v )
! 98: break;
! 99: remsdcp(vl,f,p0,&f1); diffp(vl,f1,VR(f1),&f2);
! 100: while ( 1 ) {
! 101: pmonic(vl,p0,f2,&t); f2 = t; rema(vl,p0,f1,f2,&r);
! 102: if ( !r ) {
! 103: pmonic(vl,p0,f2,&gcd); pqra(vl,p0,f,gcd,&flat,&r);
! 104: break;
! 105: }
! 106: f1 = f2; f2 = r;
! 107: }
! 108: if ( VR(gcd) != v ) {
! 109: a[i] = f; b[i] = 1;
! 110: break;
! 111: }
! 112: for ( j = 1, f = gcd; ; j++ ) {
! 113: pqra(vl,p0,f,flat,&q,&r);
! 114: if ( r )
! 115: break;
! 116: else
! 117: f = q;
! 118: }
! 119: a[i] = flat; b[i] = j;
! 120: }
! 121: for ( i = 0, j = 0, dc0 = 0; a[i]; i++ ) {
! 122: NEXTDC(dc0,dc); j += b[i]; STOQ(j,DEG(dc));
! 123: if ( a[i+1] )
! 124: pqra(vl,p0,a[i],a[i+1],&COEF(dc),&t);
! 125: else
! 126: COEF(dc) = a[i];
! 127: }
! 128: NEXT(dc) = 0; *dcp = dc0;
! 129: }
! 130: }
! 131:
! 132: void pqra(vl, p0, p1, p2, p, pr)
! 133: VL vl;
! 134: P p0,p1,p2;
! 135: P *p,*pr;
! 136: {
! 137: V v,v1;
! 138: P m,m1,q,r;
! 139: Q tq;
! 140: int n;
! 141:
! 142: if ( ( v = VR(p1) ) != ( v1 = VR(p2) ) ) {
! 143: while ( ( VR(vl) != v ) && ( VR(vl) != v1 ) )
! 144: vl = NEXT(vl);
! 145: if ( VR(vl) == v ) {
! 146: *p = p1; *pr = 0;
! 147: } else {
! 148: *p = 0; *pr = p1;
! 149: }
! 150: } else if ( (n = deg(v,p1) - deg(v,p2)) < 0 ) {
! 151: *pr = p1; *p = 0;
! 152: } else {
! 153: n++; STOQ(n,tq);
! 154: pwrp(vl,LC(p2),tq,&m); mulp(vl,m,p1,&m1); divsrp(vl,m1,p2,&q,&r);
! 155: remsdcp(vl,q,p0,p); remsdcp(vl,r,p0,pr);
! 156: }
! 157: }
! 158:
! 159: void pmonic(vl, p0, p, pr)
! 160: VL vl;
! 161: P p,p0,*pr;
! 162: {
! 163: P d,m,pp,c,t;
! 164:
! 165: if ( NUM(p) || ( VR(p) == VR(p0) ) )
! 166: *pr = (P)ONE;
! 167: else if ( NUM(COEF(DC(p))) )
! 168: *pr = p;
! 169: else {
! 170: ptozp(COEF(DC(p)),1,(Q *)&c,&pp); pinva(p0,pp,&d);
! 171: mulp(vl,p,d,&m); remsdcp(vl,m,p0,&t); ptozp(t,1,(Q *)&c,pr);
! 172: }
! 173: }
! 174:
! 175: void remsdcp(vl, p, p0, pr)
! 176: VL vl;
! 177: P p,p0,*pr;
! 178: {
! 179: P q,r,c;
! 180: DCP dc,dcr,dcr0;
! 181:
! 182: if ( !p )
! 183: *pr = 0;
! 184: else if ( NUM(p) )
! 185: *pr = (P)ONE;
! 186: else if ( VR(p) == VR(p0) ) {
! 187: divsrp(vl,p,p0,&q,&r); ptozp(r,1,(Q *)&c,pr);
! 188: } else {
! 189: for ( dcr = dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
! 190: divsrp(vl,COEF(dc),p0,&q,&r);
! 191: if ( r ) {
! 192: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = r;
! 193: }
! 194: }
! 195: if ( !dcr0 )
! 196: *pr = 0;
! 197: else {
! 198: NEXT(dcr) = 0; MKP(VR(p),dcr0,r); ptozp(r,1,(Q *)&c,pr);
! 199: }
! 200: }
! 201: }
! 202:
! 203: void pinva(p0, p, pr)
! 204: P p,p0,*pr;
! 205: {
! 206: V v;
! 207: Q f,q,q1;
! 208: Q u,t,a,b,s;
! 209: P c,c1;
! 210: P tg,th,tr;
! 211: P z,zz,zzz;
! 212: UM wg,wh,ws,wt,wm;
! 213: int n,m,modres,mod,index,i;
! 214:
! 215: v = VR(p); n=deg(v,p); m=deg(v,p0);
! 216: norm(p,&a); norm(p0,&b);
! 217: STOQ(m,u); pwrq(a,u,&t); STOQ(n,u); pwrq(b,u,&s);
! 218: mulq(t,s,&u);
! 219: factorial(n+m,&t); mulq(u,t,&s); addq(s,s,&f);
! 220: wg = W_UMALLOC(m+n); wh = W_UMALLOC(m+n);
! 221: wt = W_UMALLOC(m+n); ws = W_UMALLOC(m+n);
! 222: wm = W_UMALLOC(m+n);
! 223: for ( q = ONE, t = 0, c = 0, index = 0; ; ) {
! 224: mod = lprime[index++];
! 225: if ( !mod )
! 226: error("pinva : lprime[] exhausted.");
! 227: if ( !rem(NM((Q)LC(p)),mod) || !rem(NM((Q)LC(p0)),mod) )
! 228: continue;
! 229: ptomp(mod,p,&tg); ptomp(mod,p0,&th); srchump(mod,tg,th,&tr);
! 230: if ( !tr )
! 231: continue;
! 232: else
! 233: modres = CONT((MQ)tr);
! 234: mptoum(tg,wg); mptoum(th,wh);
! 235: eucum(mod,wg,wh,ws,wt); /* wg * ws + wh * wt = 1 */
! 236: DEG(wm) = 0; COEF(wm)[0] = modres;
! 237: mulum(mod,ws,wm,wt);
! 238: for ( i = DEG(wt); i >= 0; i-- )
! 239: if ( ( COEF(wt)[i] * 2 ) > mod )
! 240: COEF(wt)[i] -= mod;
! 241: chnrem(mod,v,c,q,wt,&c1,&q1);
! 242: if ( !ucmpp(c,c1) ) {
! 243: mulp(CO,c,p,&z); divsrp(CO,z,p0,&zz,&zzz);
! 244: if ( NUM(zzz) ) {
! 245: q = q1; c = c1; break;
! 246: }
! 247: }
! 248: q = q1; c = c1;
! 249: if ( cmpq(f,q) < 0 )
! 250: break;
! 251: }
! 252: ptozp(c,1,&s,pr);
! 253: }