![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to PUM.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 PUM.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 ! noro 1: /* $OpenXM: OpenXM/src/asir99/engine/PUM.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
! 2: #include "ca.h"
! 3:
! 4: void gcdcmp(), sprsm();
! 5:
! 6: void detmp(vl,mod,rmat,n,dp)
! 7: VL vl;
! 8: int mod;
! 9: P **rmat;
! 10: int n;
! 11: P *dp;
! 12: {
! 13: int i,j,k,sgn;
! 14: P mjj,mij,t,s,u,d;
! 15: P **mat;
! 16: P *mi,*mj;
! 17:
! 18: mat = (P **)almat_pointer(n,n);
! 19: for ( i = 0; i < n; i++ )
! 20: for ( j = 0; j < n; j++ )
! 21: mat[i][j] = rmat[i][j];
! 22: for ( j = 0, d = (P)ONEM, sgn = 1; j < n; j++ ) {
! 23: for ( i = j; (i < n) && !mat[i][j]; i++ );
! 24: if ( i == n ) {
! 25: *dp = 0; return;
! 26: }
! 27: for ( k = i; k < n; k++ )
! 28: if ( mat[k][j] && (nmonop(mat[k][j]) < nmonop(mat[i][j]) ) )
! 29: i = k;
! 30: if ( j != i ) {
! 31: mj = mat[j]; mat[j] = mat[i]; mat[i] = mj; sgn = -sgn;
! 32: }
! 33: for ( i = j + 1, mj = mat[j], mjj = mj[j]; i < n; i++ )
! 34: for ( k = j + 1, mi = mat[i], mij = mi[j]; k < n; k++ ) {
! 35: mulmp(vl,mod,mi[k],mjj,&t); mulmp(vl,mod,mj[k],mij,&s);
! 36: submp(vl,mod,t,s,&u); divsmp(vl,mod,u,d,&mi[k]);
! 37: }
! 38: d = mjj;
! 39: }
! 40: if ( sgn < 0 )
! 41: chsgnmp(mod,d,dp);
! 42: else
! 43: *dp = d;
! 44: }
! 45:
! 46: void resultmp(vl,mod,v,p1,p2,pr)
! 47: VL vl;
! 48: int mod;
! 49: V v;
! 50: P p1,p2,*pr;
! 51: {
! 52: P q1,q2,m,m1,m2,lc,q,r,t,g1,g2,adj;
! 53: int d,d1,d2,j,k;
! 54: VL nvl;
! 55: Q dq;
! 56: MQ mq;
! 57:
! 58: if ( !p1 || !p2 ) {
! 59: *pr = 0; return;
! 60: }
! 61: reordvar(vl,v,&nvl);
! 62: reordermp(nvl,mod,vl,p1,&q1); reordermp(nvl,mod,vl,p2,&q2);
! 63:
! 64: if ( VR(q1) != v )
! 65: if ( VR(q2) != v ) {
! 66: *pr = 0;
! 67: return;
! 68: } else {
! 69: d = deg(v,q2); STOQ(d,dq);
! 70: pwrmp(vl,mod,q1,dq,pr);
! 71: return;
! 72: }
! 73: else if ( VR(q2) != v ) {
! 74: d = deg(v,q1); STOQ(d,dq);
! 75: pwrmp(vl,mod,q2,dq,pr);
! 76: return;
! 77: }
! 78:
! 79: if ( ( VR(q1) != v ) || ( VR(q2) != v ) ) {
! 80: *pr = 0;
! 81: return;
! 82: }
! 83:
! 84: d1 = deg(v,q1); d2 = deg(v,q2);
! 85: if ( d1 > d2 ) {
! 86: g1 = q1; g2 = q2; adj = (P)ONEM;
! 87: } else if ( d1 < d2 ) {
! 88: g2 = q1; g1 = q2;
! 89: if ( (d1 % 2) && (d2 % 2) ) {
! 90: STOMQ(-1,mq); adj = (P)mq;
! 91: } else
! 92: adj = (P)ONEM;
! 93: } else {
! 94: premmp(nvl,mod,q1,q2,&t);
! 95: d = deg(v,t); STOQ(d,dq); pwrmp(nvl,mod,LC(q2),dq,&adj);
! 96: g1 = q2; g2 = t;
! 97: if ( d1 % 2 ) {
! 98: chsgnmp(mod,adj,&t); adj = t;
! 99: }
! 100: }
! 101: d1 = deg(v,g1); j = d1 - 1;
! 102:
! 103: for ( lc = (P)ONEM; ; ) {
! 104: if ( ( k = deg(v,g2) ) < 0 ) {
! 105: *pr = 0;
! 106: return;
! 107: }
! 108:
! 109: if ( k == j )
! 110: if ( !k ) {
! 111: divsmp(nvl,mod,g2,adj,pr);
! 112: return;
! 113: } else {
! 114: premmp(nvl,mod,g1,g2,&r); mulmp(nvl,mod,lc,lc,&m);
! 115: divsmp(nvl,mod,r,m,&q);
! 116: g1 = g2; g2 = q;
! 117: lc = LC(g1); /* g1 is not const */
! 118: j = k - 1;
! 119: }
! 120: else {
! 121: d = j - k; STOQ(d,dq);
! 122: pwrmp(nvl,mod,(VR(g2)==v?LC(g2):g2),dq,&m);
! 123: mulmp(nvl,mod,g2,m,&m1);
! 124: pwrmp(nvl,mod,lc,dq,&m); divsmp(nvl,mod,m1,m,&t);
! 125: if ( k == 0 ) {
! 126: divsmp(nvl,mod,t,adj,pr);
! 127: return;
! 128: } else {
! 129: premmp(nvl,mod,g1,g2,&r);
! 130: mulmp(nvl,mod,lc,lc,&m1); mulmp(nvl,mod,m,m1,&m2);
! 131: divsmp(nvl,mod,r,m2,&q);
! 132: g1 = t; g2 = q;
! 133: if ( d % 2 ) {
! 134: chsgnmp(mod,g2,&t); g2 = t;
! 135: }
! 136: lc = LC(g1); /* g1 is not const */
! 137: j = k - 1;
! 138: }
! 139: }
! 140: }
! 141: }
! 142:
! 143: void premmp(vl,mod,p1,p2,pr)
! 144: VL vl;
! 145: int mod;
! 146: P p1,p2,*pr;
! 147: {
! 148: P m,m1,m2;
! 149: P *pw;
! 150: DCP dc;
! 151: V v1,v2;
! 152: register int i,j;
! 153: int n1,n2,d;
! 154:
! 155: if ( NUM(p1) )
! 156: if ( NUM(p2) )
! 157: *pr = 0;
! 158: else
! 159: *pr = p1;
! 160: else if ( NUM(p2) )
! 161: *pr = 0;
! 162: else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
! 163: n1 = deg(v1,p1); n2 = deg(v1,p2);
! 164: pw = (P *)ALLOCA((n1+1)*sizeof(P));
! 165: bzero((char *)pw,(int)((n1+1)*sizeof(P)));
! 166:
! 167: for ( dc = DC(p1); dc; dc = NEXT(dc) )
! 168: pw[QTOS(DEG(dc))] = COEF(dc);
! 169:
! 170: for ( i = n1; i >= n2; i-- ) {
! 171: if ( pw[i] ) {
! 172: m = pw[i];
! 173: for ( j = i; j >= 0; j-- ) {
! 174: mulmp(vl,mod,pw[j],LC(p2),&m1); pw[j] = m1;
! 175: }
! 176:
! 177: for ( dc = DC(p2), d = i - n2; dc; dc = NEXT(dc) ) {
! 178: mulmp(vl,mod,COEF(dc),m,&m1);
! 179: submp(vl,mod,pw[QTOS(DEG(dc))+d],m1,&m2);
! 180: pw[QTOS(DEG(dc))+d] = m2;
! 181: }
! 182: } else
! 183: for ( j = i; j >= 0; j-- )
! 184: if ( pw[j] ) {
! 185: mulmp(vl,mod,pw[j],LC(p2),&m1); pw[j] = m1;
! 186: }
! 187: }
! 188: plisttop(pw,v1,n2-1,pr);
! 189: } else {
! 190: while ( v1 != vl->v && v2 != vl->v )
! 191: vl = NEXT(vl);
! 192: if ( v1 == vl->v )
! 193: *pr = 0;
! 194: else
! 195: *pr = p1;
! 196: }
! 197: }
! 198:
! 199: void srchmp(vl,mod,v,p1,p2,pr)
! 200: VL vl;
! 201: int mod;
! 202: V v;
! 203: P p1,p2,*pr;
! 204: {
! 205: P a,b,q1,q2,x,t,s,d,bg,c,c0,db;
! 206: int i,m,k;
! 207: V v0;
! 208: VL nvl,tvl,nvl1,nvl2;
! 209: Q dq;
! 210: MQ q;
! 211:
! 212: if ( vl->v != v ) {
! 213: reordvar(vl,v,&tvl);
! 214: reordermp(tvl,mod,vl,p1,&q1); reordermp(tvl,mod,vl,p2,&q2);
! 215: } else {
! 216: q1 = p1; q2 = p2; tvl = vl;
! 217: }
! 218: clctv(tvl,q1,&nvl1); clctv(tvl,q2,&nvl2); mergev(tvl,nvl1,nvl2,&nvl);
! 219: if ( VR(q1) != v )
! 220: if ( VR(q2) != v )
! 221: *pr = 0;
! 222: else {
! 223: m = getdeg(v,q2); STOQ(m,dq); pwrmp(vl,mod,q1,dq,pr);
! 224: }
! 225: else if ( VR(q2) != v ) {
! 226: m = getdeg(v,q1); STOQ(m,dq); pwrmp(vl,mod,q2,dq,pr);
! 227: } else if ( !NEXT(nvl) )
! 228: srchump(mod,p1,p2,pr);
! 229: else {
! 230: v0 = NEXT(nvl)->v;
! 231: k = getdeg(v,q1)*getdeg(v0,q2)+getdeg(v,q2)*getdeg(v0,q1)+1;
! 232: for ( i = 0, c = 0, d = (P)ONEM, MKMV(v0,x);
! 233: ( i < mod ) && (getdeg(v0,d) < k) ; i++ ) {
! 234: STOMQ(i,q),bg = (P)q; substmp(nvl,mod,LC(q1),v0,bg,&t);
! 235: if ( !t )
! 236: continue;
! 237: substmp(nvl,mod,LC(q2),v0,bg,&t);
! 238: if ( !t )
! 239: continue;
! 240: substmp(nvl,mod,q1,v0,bg,&a); substmp(nvl,mod,q2,v0,bg,&b);
! 241: srchmp(nvl,mod,v,a,b,&c0); substmp(nvl,mod,c,v0,bg,&t);
! 242: submp(nvl,mod,c0,t,&s); mulmp(nvl,mod,s,d,&t);
! 243: substmp(nvl,mod,d,v0,bg,&db);
! 244: divsmp(nvl,mod,t,db,&s); addmp(nvl,mod,s,c,&t); c = t;
! 245: submp(nvl,mod,x,bg,&t); mulmp(nvl,mod,d,t,&s); d = s;
! 246: }
! 247: if ( i == mod )
! 248: error("srchmp : ???");
! 249: *pr = c;
! 250: }
! 251: }
! 252:
! 253: int ucmpp(p,q)
! 254: P p,q;
! 255: {
! 256: DCP dcp,dcq;
! 257:
! 258: if ( !p )
! 259: if ( !q )
! 260: return ( 0 );
! 261: else
! 262: return ( 1 );
! 263: else if ( !q )
! 264: return ( 1 );
! 265: else if ( NUM(p) )
! 266: if ( !NUM(q) )
! 267: return ( 1 );
! 268: else
! 269: return ( cmpq((Q)p,(Q)q) );
! 270: else if ( NUM(q) )
! 271: return ( 1 );
! 272: else {
! 273: for ( dcp = DC(p), dcq = DC(q); dcp && dcq;
! 274: dcp = NEXT(dcp), dcq = NEXT(dcq) )
! 275: if ( cmpq(DEG(dcp),DEG(dcq) ) )
! 276: return ( 1 );
! 277: else if ( cmpq((Q)COEF(dcp),(Q)COEF(dcq) ) )
! 278: return ( 1 );
! 279: if ( dcp || dcq )
! 280: return ( 1 );
! 281: else
! 282: return ( 0 );
! 283: }
! 284: }
! 285:
! 286: #if 0
! 287: srchump(mod,p1,p2,pr)
! 288: int mod;
! 289: P p1,p2,*pr;
! 290: {
! 291: int r;
! 292: MQ mq;
! 293:
! 294: r = eucresum(mod,p1,p2);
! 295: STOMQ(r,mq); *pr = (P)mq;
! 296: }
! 297: #endif
! 298:
! 299: void srchump(mod,p1,p2,pr)
! 300: int mod;
! 301: P p1,p2,*pr;
! 302: {
! 303: UM m,m1,q,r,t,g1,g2;
! 304: int lc,d,d1,d2,i,j,k,l,l1,l2,tmp,adj;
! 305: V v;
! 306:
! 307: v = VR(p1); d = MAX(UDEG(p1),UDEG(p2));
! 308: g1 = W_UMALLOC(d); g2 = W_UMALLOC(d);
! 309: bzero((char *)g1,(int)((d+2)*sizeof(int))); bzero((char *)g2,(int)((d+2)*sizeof(int)));
! 310: if ( d == (int)UDEG(p1) ) {
! 311: mptoum(p1,g1); mptoum(p2,g2);
! 312: } else {
! 313: mptoum(p2,g1); mptoum(p1,g2);
! 314: }
! 315: if ( ( d1 = DEG(g1) ) > ( d2 = DEG(g2) ) ) {
! 316: j = d1 - 1; adj = 1;
! 317: } else
! 318: j = d2;
! 319: lc = 1;
! 320: r = W_UMALLOC(d1+d2); q = W_UMALLOC(d1+d2);
! 321: m1 = W_UMALLOC(d1+d2); t = W_UMALLOC(d1+d2);
! 322: bzero((char *)r,(int)((d1+d2+2)*sizeof(int))); bzero((char *)q,(int)((d1+d2+2)*sizeof(int)));
! 323: bzero((char *)m1,(int)((d1+d2+2)*sizeof(int))); bzero((char *)t,(int)((d1+d2+2)*sizeof(int)));
! 324: m = W_UMALLOC(0); bzero((char *)m,(int)(2*sizeof(int)));
! 325: adj = pwrm(mod,COEF(g2)[DEG(g2)],DEG(g1));
! 326: DEG(m) = 0; COEF(m)[0] = invm(COEF(g2)[DEG(g2)],mod);
! 327: mulum(mod,g2,m,r); cpyum(r,g2);
! 328: while ( 1 ) {
! 329: if ( ( k = DEG(g2) ) < 0 ) {
! 330: *pr = 0;
! 331: return;
! 332: }
! 333: if ( k == j ) {
! 334: if ( k == 0 ) {
! 335: DEG(m) = 0; COEF(m)[0] = adj;
! 336: mulum(mod,g2,m,r); umtomp(v,r,pr);
! 337: return;
! 338: } else {
! 339: DEG(m) = 0;
! 340: COEF(m)[0] = pwrm(mod,COEF(g2)[k],DEG(g1)-k+1);
! 341: mulum(mod,g1,m,r); DEG(r) = divum(mod,r,g2,t);
! 342: DEG(m) = 0; COEF(m)[0] = dmb(mod,lc,lc,&tmp);
! 343: divum(mod,r,m,q); cpyum(g2,g1); cpyum(q,g2);
! 344: lc = COEF(g1)[DEG(g1)]; j = k - 1;
! 345: }
! 346: } else {
! 347: d = j - k;
! 348: DEG(m) = 0; COEF(m)[0] = pwrm(mod,COEF(g2)[DEG(g2)],d);
! 349: mulum(mod,g2,m,m1); l = pwrm(mod,lc,d);
! 350: DEG(m) = 0; COEF(m)[0] = l; divum(mod,m1,m,t);
! 351: if ( k == 0 ) {
! 352: DEG(m) = 0; COEF(m)[0] = adj;
! 353: mulum(mod,t,m,r); umtomp(v,r,pr);
! 354: return;
! 355: } else {
! 356: DEG(m) = 0;
! 357: COEF(m)[0] = pwrm(mod,COEF(g2)[k],DEG(g1)-k+1);
! 358: mulum(mod,g1,m,r); DEG(r) = divum(mod,r,g2,q);
! 359: l1 = dmb(mod,lc,lc,&tmp); l2 = dmb(mod,l,l1,&tmp);
! 360: DEG(m) = 0; COEF(m)[0] = l2;
! 361: divum(mod,r,m,q); cpyum(t,g1); cpyum(q,g2);
! 362: if ( d % 2 )
! 363: for ( i = DEG(g2); i >= 0; i-- )
! 364: COEF(g2)[i] = ( mod - COEF(g2)[i] ) % mod;
! 365: lc = COEF(g1)[DEG(g1)]; j = k - 1;
! 366: }
! 367: }
! 368: }
! 369: }
! 370:
! 371: void substmp(vl,mod,p,v0,p0,pr)
! 372: VL vl;
! 373: int mod;
! 374: V v0;
! 375: P p,p0;
! 376: P *pr;
! 377: {
! 378: P x,t,m,c,s,a;
! 379: DCP dc;
! 380: Q d;
! 381:
! 382: if ( !p )
! 383: *pr = 0;
! 384: else if ( NUM(p) )
! 385: *pr = p;
! 386: else if ( VR(p) != v0 ) {
! 387: MKMV(VR(p),x);
! 388: for ( c = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
! 389: substmp(vl,mod,COEF(dc),v0,p0,&t);
! 390: if ( DEG(dc) ) {
! 391: pwrmp(vl,mod,x,DEG(dc),&s); mulmp(vl,mod,s,t,&m);
! 392: addmp(vl,mod,m,c,&a);
! 393: c = a;
! 394: } else {
! 395: addmp(vl,mod,t,c,&a);
! 396: c = a;
! 397: }
! 398: }
! 399: *pr = c;
! 400: } else {
! 401: dc = DC(p);
! 402: c = COEF(dc);
! 403: for ( d = DEG(dc), dc = NEXT(dc);
! 404: dc; d = DEG(dc), dc = NEXT(dc) ) {
! 405: subq(d,DEG(dc),(Q *)&t); pwrmp(vl,mod,p0,(Q)t,&s);
! 406: mulmp(vl,mod,s,c,&m);
! 407: addmp(vl,mod,m,COEF(dc),&c);
! 408: }
! 409: if ( d ) {
! 410: pwrmp(vl,mod,p0,d,&t); mulmp(vl,mod,t,c,&m);
! 411: c = m;
! 412: }
! 413: *pr = c;
! 414: }
! 415: }
! 416:
! 417: void reordermp(nvl,mod,ovl,p,pr)
! 418: VL nvl,ovl;
! 419: int mod;
! 420: P p;
! 421: P *pr;
! 422: {
! 423: DCP dc;
! 424: P x,m,s,t,c;
! 425:
! 426: if ( !p )
! 427: *pr = 0;
! 428: else if ( NUM(p) )
! 429: *pr = p;
! 430: else {
! 431: MKMV(VR(p),x);
! 432: for ( s = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
! 433: reordermp(nvl,mod,ovl,COEF(dc),&c);
! 434: if ( DEG(dc) ) {
! 435: pwrmp(nvl,mod,x,DEG(dc),&t); mulmp(nvl,mod,c,t,&m);
! 436: addmp(nvl,mod,m,s,&t);
! 437: } else
! 438: addmp(nvl,mod,s,c,&t);
! 439: s = t;
! 440: }
! 441: *pr = s;
! 442: }
! 443: }
! 444:
! 445: void chnremp(vl,mod,p,q,c,r)
! 446: VL vl;
! 447: int mod;
! 448: P p;
! 449: Q q;
! 450: P c;
! 451: P *r;
! 452: {
! 453: P tg,sg,ug;
! 454: P t,u;
! 455: MQ mq;
! 456:
! 457: ptomp(mod,p,&tg); submp(vl,mod,c,tg,&sg);
! 458: UTOMQ(rem(NM(q),mod),mq),tg = (P)mq; divsmp(vl,mod,sg,tg,&ug);
! 459: normalizemp(mod,ug);
! 460: mptop(ug,&u); mulp(vl,u,(P)q,&t); addp(vl,t,p,r);
! 461: }
! 462:
! 463: /* XXX strange behavior of invm() on SPARC */
! 464:
! 465: void chnrem(mod,v,c,q,t,cr,qr)
! 466: int mod;
! 467: V v;
! 468: UM t;
! 469: Q q,*qr;
! 470: P c,*cr;
! 471: {
! 472: int n,m,i,d,a,sd,tmp;
! 473: Q b,s,z;
! 474: Q *pc,*pcr;
! 475: DCP dc;
! 476:
! 477: if ( !c || NUM(c) )
! 478: n = 0;
! 479: else
! 480: n = UDEG(c);
! 481: m = DEG(t); d = MAX(n,m); W_CALLOC(n,Q,pc); W_CALLOC(d,Q,pcr);
! 482: if ( !c )
! 483: pc[0] = 0;
! 484: else if ( NUM(c) )
! 485: pc[0] = (Q)c;
! 486: else
! 487: for ( dc = DC(c); dc; dc = NEXT(dc) )
! 488: pc[QTOS(DEG(dc))] = (Q)COEF(dc);
! 489: for ( i = 0; i <= d; i++ ) {
! 490: b = (i>n?0:pc[i]); a = (i>m?0:COEF(t)[i]);
! 491: if ( b )
! 492: a = (a-SGN(pc[i])*((int)rem(NM(pc[i]),mod)))%mod;
! 493: sd = dmb(mod,(a>=0?a:a+mod),invm(rem(NM(q),mod),mod),&tmp);
! 494: if ( ( 2 * sd ) > mod )
! 495: sd -= mod;
! 496: STOQ(sd,z); mulq(z,q,&s); addq(s,b,&pcr[i]);
! 497: }
! 498: STOQ(mod,z); mulq(q,z,qr); plisttop((P *)pcr,v,d,cr);
! 499: }
! 500:
! 501: void normalizemp(mod,g)
! 502: int mod;
! 503: P g;
! 504: {
! 505: DCP dc;
! 506:
! 507: if ( !g )
! 508: return;
! 509: else if ( NUM(g) ) {
! 510: if ( 2 * CONT((MQ)g) > mod )
! 511: CONT((MQ)g) -= mod;
! 512: return;
! 513: } else
! 514: for ( dc = DC(g); dc; dc = NEXT(dc) )
! 515: normalizemp(mod,COEF(dc));
! 516: }
! 517:
! 518: void norm(p,r)
! 519: P p;
! 520: Q *r;
! 521: {
! 522: N t;
! 523: DCP dc;
! 524:
! 525: for ( dc = DC(p), t = ONEN; dc; dc = NEXT(dc) )
! 526: if ( cmpn(NM((Q)COEF(dc)),t) > 0 ) {
! 527: t = NM((Q)COEF(dc));
! 528: }
! 529: NTOQ(t,1,*r);
! 530: }
! 531:
! 532: void norm1(p,r)
! 533: P p,*r;
! 534: {
! 535: DCP dc;
! 536: P t,s,u;
! 537: Q q;
! 538:
! 539: if ( NUM(p) )
! 540: NTOQ(NM((Q)p),1,q),*r = (P)q;
! 541: else {
! 542: for ( t = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
! 543: norm1(COEF(dc),&s); addq((Q)t,(Q)s,(Q *)&u); t = u;
! 544: }
! 545: *r = t;
! 546: }
! 547: }
! 548:
! 549: void norm1c(p,r)
! 550: P p;
! 551: Q *r;
! 552: {
! 553: N t;
! 554: Q s;
! 555: DCP dc;
! 556:
! 557: if ( NUM(p) )
! 558: norm1(p,(P *)r);
! 559: else {
! 560: for ( dc = DC(p), t = ONEN; dc; dc = NEXT(dc) ) {
! 561: norm1(COEF(dc),(P *)&s);
! 562: if ( cmpn(NM(s),t) > 0 )
! 563: t = NM(s);
! 564: }
! 565: NTOQ(t,1,*r);
! 566: }
! 567: }
! 568:
! 569: void gcdprsmp(vl,mod,p1,p2,pr)
! 570: VL vl;
! 571: int mod;
! 572: P p1,p2,*pr;
! 573: {
! 574: P g1,g2,gc1,gc2,gp1,gp2,g,gc,gp,gcr;
! 575: V v1,v2;
! 576:
! 577: if ( !p1 )
! 578: *pr = p2;
! 579: else if ( !p2 )
! 580: *pr = p1;
! 581: else if ( NUM(p1) || NUM(p2) )
! 582: *pr = (P)ONEM;
! 583: else {
! 584: g1 = p1; g2 = p2;
! 585: if ( ( v1 = VR(g1) ) == ( v2 = VR(g2) ) ) {
! 586: gcdcmp(vl,mod,g1,&gc1); divsmp(vl,mod,g1,gc1,&gp1);
! 587: gcdcmp(vl,mod,g2,&gc2); divsmp(vl,mod,g2,gc2,&gp2);
! 588: gcdprsmp(vl,mod,gc1,gc2,&gcr);
! 589: sprsm(vl,mod,v1,gp1,gp2,&g);
! 590:
! 591: if ( VR(g) == v1 ) {
! 592: gp = g;
! 593: gcdcmp(vl,mod,gp,&gc); divsmp(vl,mod,gp,gc,&gp1);
! 594: mulmp(vl,mod,gp1,gcr,pr);
! 595: } else
! 596: *pr = gcr;
! 597: } else {
! 598: while ( v1 != vl->v && v2 != vl->v )
! 599: vl = NEXT(vl);
! 600: if ( v1 == vl->v ) {
! 601: gcdcmp(vl,mod,g1,&gc1); gcdprsmp(vl,mod,gc1,g2,pr);
! 602: } else {
! 603: gcdcmp(vl,mod,g2,&gc2); gcdprsmp(vl,mod,gc2,g1,pr);
! 604: }
! 605: }
! 606: }
! 607: }
! 608:
! 609: void gcdcmp(vl,mod,p,pr)
! 610: VL vl;
! 611: int mod;
! 612: P p,*pr;
! 613: {
! 614: P g,g1;
! 615: DCP dc;
! 616:
! 617: if ( NUM(p) )
! 618: *pr = (P)ONEM;
! 619: else {
! 620: dc = DC(p);
! 621: g = COEF(dc);
! 622: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
! 623: gcdprsmp(vl,mod,g,COEF(dc),&g1);
! 624: g = g1;
! 625: }
! 626: *pr = g;
! 627: }
! 628: }
! 629:
! 630: void sprsm(vl,mod,v,p1,p2,pr)
! 631: VL vl;
! 632: int mod;
! 633: V v;
! 634: P p1,p2,*pr;
! 635: {
! 636: P q1,q2,m,m1,m2,x,h,r,g1,g2;
! 637: int d;
! 638: Q dq;
! 639: VL nvl;
! 640:
! 641: reordvar(vl,v,&nvl);
! 642: reordermp(nvl,mod,vl,p1,&q1); reordermp(nvl,mod,vl,p2,&q2);
! 643:
! 644: if ( ( VR(q1) != v ) || ( VR(q2) != v ) ) {
! 645: *pr = 0;
! 646: return;
! 647: }
! 648:
! 649: if ( deg(v,q1) >= deg(v,q2) ) {
! 650: g1 = q1; g2 = q2;
! 651: } else {
! 652: g2 = q1; g1 = q2;
! 653: }
! 654:
! 655: for ( h = (P)ONEM, x = (P)ONEM; ; ) {
! 656: if ( !deg(v,g2) )
! 657: break;
! 658:
! 659: premmp(nvl,mod,g1,g2,&r);
! 660: if ( !r )
! 661: break;
! 662:
! 663: d = deg(v,g1) - deg(v,g2); STOQ(d,dq);
! 664: pwrmp(nvl,mod,h,dq,&m); mulmp(nvl,mod,m,x,&m1); g1 = g2;
! 665: divsmp(nvl,mod,r,m1,&g2); x = LC(g1); /* g1 is not const w.r.t v */
! 666: pwrmp(nvl,mod,x,dq,&m1); mulmp(nvl,mod,m1,h,&m2);
! 667: divsmp(nvl,mod,m2,m,&h);
! 668: }
! 669: *pr = g2;
! 670: }