![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to PU.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 PU.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 ! noro 1: /* $OpenXM: OpenXM/src/asir99/engine/PU.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
! 2: #include "ca.h"
! 3:
! 4: void reorderp(nvl,ovl,p,pr)
! 5: VL nvl,ovl;
! 6: P p;
! 7: P *pr;
! 8: {
! 9: DCP dc;
! 10: P x,m,s,t,c;
! 11:
! 12: if ( !p )
! 13: *pr = 0;
! 14: else if ( NUM(p) )
! 15: *pr = p;
! 16: else {
! 17: MKV(VR(p),x);
! 18: for ( s = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
! 19: reorderp(nvl,ovl,COEF(dc),&c);
! 20: if ( DEG(dc) ) {
! 21: pwrp(nvl,x,DEG(dc),&t); mulp(nvl,c,t,&m);
! 22: addp(nvl,m,s,&t);
! 23: } else
! 24: addp(nvl,s,c,&t);
! 25: s = t;
! 26: }
! 27: *pr = s;
! 28: }
! 29: }
! 30:
! 31: void substp(vl,p,v0,p0,pr)
! 32: VL vl;
! 33: V v0;
! 34: P p,p0;
! 35: P *pr;
! 36: {
! 37: P x,t,m,c,s,a;
! 38: DCP dc;
! 39: Q d;
! 40:
! 41: if ( !p )
! 42: *pr = 0;
! 43: else if ( NUM(p) )
! 44: *pr = p;
! 45: else if ( VR(p) != v0 ) {
! 46: MKV(VR(p),x);
! 47: for ( c = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
! 48: substp(vl,COEF(dc),v0,p0,&t);
! 49: if ( DEG(dc) ) {
! 50: pwrp(vl,x,DEG(dc),&s); mulp(vl,s,t,&m);
! 51: addp(vl,m,c,&a);
! 52: c = a;
! 53: } else {
! 54: addp(vl,t,c,&a);
! 55: c = a;
! 56: }
! 57: }
! 58: *pr = c;
! 59: } else {
! 60: dc = DC(p);
! 61: c = COEF(dc);
! 62: for ( d = DEG(dc), dc = NEXT(dc);
! 63: dc; d = DEG(dc), dc = NEXT(dc) ) {
! 64: subq(d,DEG(dc),(Q *)&t); pwrp(vl,p0,(Q)t,&s);
! 65: mulp(vl,s,c,&m);
! 66: addp(vl,m,COEF(dc),&c);
! 67: }
! 68: if ( d ) {
! 69: pwrp(vl,p0,d,&t); mulp(vl,t,c,&m);
! 70: c = m;
! 71: }
! 72: *pr = c;
! 73: }
! 74: }
! 75:
! 76: void detp(vl,rmat,n,dp)
! 77: VL vl;
! 78: P **rmat;
! 79: int n;
! 80: P *dp;
! 81: {
! 82: int i,j,k,sgn;
! 83: P mjj,mij,t,s,u,d;
! 84: P **mat;
! 85: P *mi,*mj;
! 86:
! 87: mat = (P **)almat_pointer(n,n);
! 88: for ( i = 0; i < n; i++ )
! 89: for ( j = 0; j < n; j++ )
! 90: mat[i][j] = rmat[i][j];
! 91: for ( j = 0, d = (P)ONE, sgn = 1; j < n; j++ ) {
! 92: for ( i = j; (i < n) && !mat[i][j]; i++ );
! 93: if ( i == n ) {
! 94: *dp = 0; return;
! 95: }
! 96: for ( k = i; k < n; k++ )
! 97: if ( mat[k][j] && (nmonop(mat[k][j]) < nmonop(mat[i][j]) ) )
! 98: i = k;
! 99: if ( j != i ) {
! 100: mj = mat[j]; mat[j] = mat[i]; mat[i] = mj; sgn = -sgn;
! 101: }
! 102: for ( i = j + 1, mj = mat[j], mjj = mj[j]; i < n; i++ )
! 103: for ( k = j + 1, mi = mat[i], mij = mi[j]; k < n; k++ ) {
! 104: mulp(vl,mi[k],mjj,&t); mulp(vl,mj[k],mij,&s);
! 105: subp(vl,t,s,&u); divsp(vl,u,d,&mi[k]);
! 106: }
! 107: d = mjj;
! 108: }
! 109: if ( sgn < 0 )
! 110: chsgnp(d,dp);
! 111: else
! 112: *dp = d;
! 113: }
! 114:
! 115: void reordvar(vl,v,nvlp)
! 116: VL vl;
! 117: V v;
! 118: VL *nvlp;
! 119: {
! 120: VL nvl,nvl0;
! 121:
! 122: for ( NEWVL(nvl0), nvl0->v = v, nvl = nvl0;
! 123: vl; vl = NEXT(vl) )
! 124: if ( vl->v == v )
! 125: continue;
! 126: else {
! 127: NEWVL(NEXT(nvl));
! 128: nvl = NEXT(nvl);
! 129: nvl->v = vl->v;
! 130: }
! 131: NEXT(nvl) = 0;
! 132: *nvlp = nvl0;
! 133: }
! 134:
! 135: void gcdprsp(vl,p1,p2,pr)
! 136: VL vl;
! 137: P p1,p2,*pr;
! 138: {
! 139: P g1,g2,gc1,gc2,gp1,gp2,g,gc,gp,gcr;
! 140: V v1,v2;
! 141:
! 142: if ( !p1 )
! 143: ptozp0(p2,pr);
! 144: else if ( !p2 )
! 145: ptozp0(p1,pr);
! 146: else if ( NUM(p1) || NUM(p2) )
! 147: *pr = (P)ONE;
! 148: else {
! 149: ptozp0(p1,&g1); ptozp0(p2,&g2);
! 150: if ( ( v1 = VR(g1) ) == ( v2 = VR(g2) ) ) {
! 151: gcdcp(vl,g1,&gc1); divsp(vl,g1,gc1,&gp1);
! 152: gcdcp(vl,g2,&gc2); divsp(vl,g2,gc2,&gp2);
! 153: gcdprsp(vl,gc1,gc2,&gcr);
! 154: sprs(vl,v1,gp1,gp2,&g);
! 155:
! 156: if ( VR(g) == v1 ) {
! 157: ptozp0(g,&gp);
! 158: gcdcp(vl,gp,&gc); divsp(vl,gp,gc,&gp1);
! 159: mulp(vl,gp1,gcr,pr);
! 160:
! 161: } else
! 162: *pr = gcr;
! 163: } else {
! 164: while ( v1 != vl->v && v2 != vl->v )
! 165: vl = NEXT(vl);
! 166: if ( v1 == vl->v ) {
! 167: gcdcp(vl,g1,&gc1); gcdprsp(vl,gc1,g2,pr);
! 168: } else {
! 169: gcdcp(vl,g2,&gc2); gcdprsp(vl,gc2,g1,pr);
! 170: }
! 171: }
! 172: }
! 173: }
! 174:
! 175: void gcdcp(vl,p,pr)
! 176: VL vl;
! 177: P p,*pr;
! 178: {
! 179: P g,g1;
! 180: DCP dc;
! 181:
! 182: if ( NUM(p) )
! 183: *pr = (P)ONE;
! 184: else {
! 185: dc = DC(p);
! 186: g = COEF(dc);
! 187: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
! 188: gcdprsp(vl,g,COEF(dc),&g1);
! 189: g = g1;
! 190: }
! 191: *pr = g;
! 192: }
! 193: }
! 194:
! 195: void sprs(vl,v,p1,p2,pr)
! 196: VL vl;
! 197: V v;
! 198: P p1,p2,*pr;
! 199: {
! 200: P q1,q2,m,m1,m2,x,h,r,g1,g2;
! 201: int d;
! 202: Q dq;
! 203: VL nvl;
! 204:
! 205: reordvar(vl,v,&nvl);
! 206: reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
! 207:
! 208: if ( ( VR(q1) != v ) || ( VR(q2) != v ) ) {
! 209: *pr = 0;
! 210: return;
! 211: }
! 212:
! 213: if ( deg(v,q1) >= deg(v,q2) ) {
! 214: g1 = q1; g2 = q2;
! 215: } else {
! 216: g2 = q1; g1 = q2;
! 217: }
! 218:
! 219: for ( h = (P)ONE, x = (P)ONE; ; ) {
! 220: if ( !deg(v,g2) )
! 221: break;
! 222:
! 223: premp(nvl,g1,g2,&r);
! 224: if ( !r )
! 225: break;
! 226:
! 227: d = deg(v,g1) - deg(v,g2); STOQ(d,dq);
! 228: pwrp(nvl,h,dq,&m); mulp(nvl,m,x,&m1); g1 = g2;
! 229: divsp(nvl,r,m1,&g2); x = LC(g1); /* g1 is not const w.r.t v */
! 230: pwrp(nvl,x,dq,&m1); mulp(nvl,m1,h,&m2);
! 231: divsp(nvl,m2,m,&h);
! 232: }
! 233: *pr = g2;
! 234: }
! 235:
! 236: void resultp(vl,v,p1,p2,pr)
! 237: VL vl;
! 238: V v;
! 239: P p1,p2,*pr;
! 240: {
! 241: P q1,q2,m,m1,m2,lc,q,r,t,g1,g2,adj;
! 242: int d,d1,d2,j,k;
! 243: VL nvl;
! 244: Q dq;
! 245:
! 246: if ( !p1 || !p2 ) {
! 247: *pr = 0; return;
! 248: }
! 249: reordvar(vl,v,&nvl);
! 250: reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
! 251:
! 252: if ( VR(q1) != v )
! 253: if ( VR(q2) != v ) {
! 254: *pr = 0;
! 255: return;
! 256: } else {
! 257: d = deg(v,q2); STOQ(d,dq);
! 258: pwrp(vl,q1,dq,pr);
! 259: return;
! 260: }
! 261: else if ( VR(q2) != v ) {
! 262: d = deg(v,q1); STOQ(d,dq);
! 263: pwrp(vl,q2,dq,pr);
! 264: return;
! 265: }
! 266:
! 267: if ( ( VR(q1) != v ) || ( VR(q2) != v ) ) {
! 268: *pr = 0;
! 269: return;
! 270: }
! 271:
! 272: d1 = deg(v,q1); d2 = deg(v,q2);
! 273: if ( d1 > d2 ) {
! 274: g1 = q1; g2 = q2; adj = (P)ONE;
! 275: } else if ( d1 < d2 ) {
! 276: g2 = q1; g1 = q2;
! 277: if ( (d1 % 2) && (d2 % 2) ) {
! 278: STOQ(-1,dq); adj = (P)dq;
! 279: } else
! 280: adj = (P)ONE;
! 281: } else {
! 282: premp(nvl,q1,q2,&t);
! 283: d = deg(v,t); STOQ(d,dq); pwrp(nvl,LC(q2),dq,&adj);
! 284: g1 = q2; g2 = t;
! 285: if ( d1 % 2 ) {
! 286: chsgnp(adj,&t); adj = t;
! 287: }
! 288: }
! 289: d1 = deg(v,g1); j = d1 - 1;
! 290:
! 291: for ( lc = (P)ONE; ; ) {
! 292: if ( ( k = deg(v,g2) ) < 0 ) {
! 293: *pr = 0;
! 294: return;
! 295: }
! 296:
! 297: if ( k == j )
! 298: if ( !k ) {
! 299: divsp(nvl,g2,adj,pr);
! 300: return;
! 301: } else {
! 302: premp(nvl,g1,g2,&r); mulp(nvl,lc,lc,&m);
! 303: divsp(nvl,r,m,&q);
! 304: g1 = g2; g2 = q;
! 305: lc = LC(g1); /* g1 is not const */
! 306: j = k - 1;
! 307: }
! 308: else {
! 309: d = j - k; STOQ(d,dq);
! 310: pwrp(nvl,(VR(g2)==v?LC(g2):g2),dq,&m);
! 311: mulp(nvl,g2,m,&m1);
! 312: pwrp(nvl,lc,dq,&m); divsp(nvl,m1,m,&t);
! 313: if ( k == 0 ) {
! 314: divsp(nvl,t,adj,pr);
! 315: return;
! 316: } else {
! 317: premp(nvl,g1,g2,&r);
! 318: mulp(nvl,lc,lc,&m1); mulp(nvl,m,m1,&m2);
! 319: divsp(nvl,r,m2,&q);
! 320: g1 = t; g2 = q;
! 321: if ( d % 2 ) {
! 322: chsgnp(g2,&t); g2 = t;
! 323: }
! 324: lc = LC(g1); /* g1 is not const */
! 325: j = k - 1;
! 326: }
! 327: }
! 328: }
! 329: }
! 330:
! 331: void srch2(vl,v,p1,p2,pr)
! 332: VL vl;
! 333: V v;
! 334: P p1,p2,*pr;
! 335: {
! 336: P q1,q2,m,m1,m2,lc,q,r,t,s,g1,g2,adj;
! 337: int d,d1,d2,j,k;
! 338: VL nvl;
! 339: Q dq;
! 340:
! 341: reordvar(vl,v,&nvl);
! 342: reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
! 343:
! 344: if ( VR(q1) != v )
! 345: if ( VR(q2) != v ) {
! 346: *pr = 0;
! 347: return;
! 348: } else {
! 349: d = deg(v,q2); STOQ(d,dq);
! 350: pwrp(vl,q1,dq,pr);
! 351: return;
! 352: }
! 353: else if ( VR(q2) != v ) {
! 354: d = deg(v,q1); STOQ(d,dq);
! 355: pwrp(vl,q2,dq,pr);
! 356: return;
! 357: }
! 358:
! 359: if ( ( VR(q1) != v ) || ( VR(q2) != v ) ) {
! 360: *pr = 0;
! 361: return;
! 362: }
! 363:
! 364: if ( deg(v,q1) >= deg(v,q2) ) {
! 365: g1 = q1; g2 = q2;
! 366: } else {
! 367: g2 = q1; g1 = q2;
! 368: }
! 369:
! 370:
! 371: if ( ( d1 = deg(v,g1) ) > ( d2 = deg(v,g2) ) ) {
! 372: j = d1 - 1;
! 373: adj = (P)ONE;
! 374: } else {
! 375: premp(nvl,g1,g2,&t);
! 376: d = deg(v,t); STOQ(d,dq);
! 377: pwrp(nvl,LC(g2),dq,&adj);
! 378: g1 = g2; g2 = t;
! 379: j = deg(v,g1) - 1;
! 380: }
! 381:
! 382: for ( lc = (P)ONE; ; ) {
! 383: if ( ( k = deg(v,g2) ) < 0 ) {
! 384: *pr = 0;
! 385: return;
! 386: }
! 387:
! 388: ptozp(g1,1,(Q *)&t,&s); g1 = s; ptozp(g2,1,(Q *)&t,&s); g2 = s;
! 389: if ( k == j )
! 390: if ( !k ) {
! 391: divsp(nvl,g2,adj,pr);
! 392: return;
! 393: } else {
! 394: premp(nvl,g1,g2,&r); mulp(nvl,lc,lc,&m);
! 395: divsp(nvl,r,m,&q);
! 396: g1 = g2; g2 = q;
! 397: lc = LC(g1); /* g1 is not const */
! 398: j = k - 1;
! 399: }
! 400: else {
! 401: d = j - k; STOQ(d,dq);
! 402: pwrp(nvl,(VR(g2)==v?LC(g2):g2),dq,&m);
! 403: mulp(nvl,g2,m,&m1);
! 404: pwrp(nvl,lc,dq,&m); divsp(nvl,m1,m,&t);
! 405: if ( k == 0 ) {
! 406: divsp(nvl,t,adj,pr);
! 407: return;
! 408: } else {
! 409: premp(nvl,g1,g2,&r);
! 410: mulp(nvl,lc,lc,&m1); mulp(nvl,m,m1,&m2);
! 411: divsp(nvl,r,m2,&q);
! 412: g1 = t; g2 = q;
! 413: if ( d % 2 ) {
! 414: chsgnp(g2,&t); g2 = t;
! 415: }
! 416: lc = LC(g1); /* g1 is not const */
! 417: j = k - 1;
! 418: }
! 419: }
! 420: }
! 421: }
! 422:
! 423: void srcr(vl,v,p1,p2,pr)
! 424: VL vl;
! 425: V v;
! 426: P p1,p2,*pr;
! 427: {
! 428: P q1,q2,c,c1;
! 429: P tg,tg1,tg2,resg;
! 430: int index,mod;
! 431: Q a,b,f,q,t,s,u,n,m;
! 432: VL nvl;
! 433:
! 434: reordvar(vl,v,&nvl);
! 435: reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
! 436:
! 437: if ( ( VR(q1) != v ) && ( VR(q2) != v ) ) {
! 438: *pr = 0;
! 439: return;
! 440: }
! 441: if ( VR(q1) != v ) {
! 442: pwrp(vl,q1,DEG(DC(q2)),pr);
! 443: return;
! 444: }
! 445: if ( VR(q2) != v ) {
! 446: pwrp(vl,q2,DEG(DC(q1)),pr);
! 447: return;
! 448: }
! 449: norm1c(q1,&a); norm1c(q2,&b);
! 450: n = DEG(DC(q1)); m = DEG(DC(q2));
! 451: pwrq(a,m,&t); pwrq(b,n,&s); mulq(t,s,&u);
! 452: factorial(QTOS(n)+QTOS(m),&t);
! 453: mulq(u,t,&s); addq(s,s,&f);
! 454: for ( index = 0, q = ONE, c = 0; cmpq(f,q) >= 0; ) {
! 455: mod = lprime[index++];
! 456: if ( !mod )
! 457: error("sqfrum : lprime[] exhausted.");
! 458: ptomp(mod,LC(q1),&tg);
! 459: if ( !tg )
! 460: continue;
! 461: ptomp(mod,LC(q2),&tg);
! 462: if ( !tg )
! 463: continue;
! 464: ptomp(mod,q1,&tg1); ptomp(mod,q2,&tg2);
! 465: srchmp(nvl,mod,v,tg1,tg2,&resg);
! 466: chnremp(nvl,mod,c,q,resg,&c1); c = c1;
! 467: STOQ(mod,t); mulq(q,t,&s); q = s;
! 468: }
! 469: *pr = c;
! 470: }
! 471:
! 472: void res_ch_det(vl,v,p1,p2,pr)
! 473: VL vl;
! 474: V v;
! 475: P p1,p2,*pr;
! 476: {
! 477: P q1,q2,c,c1;
! 478: P tg,tg1,tg2,resg;
! 479: int index,mod;
! 480: Q a,b,f,q,t,s,u,n,m;
! 481: VL nvl;
! 482:
! 483: reordvar(vl,v,&nvl);
! 484: reorderp(nvl,vl,p1,&q1); reorderp(nvl,vl,p2,&q2);
! 485:
! 486: if ( ( VR(q1) != v ) && ( VR(q2) != v ) ) {
! 487: *pr = 0;
! 488: return;
! 489: }
! 490: if ( VR(q1) != v ) {
! 491: pwrp(vl,q1,DEG(DC(q2)),pr);
! 492: return;
! 493: }
! 494: if ( VR(q2) != v ) {
! 495: pwrp(vl,q2,DEG(DC(q1)),pr);
! 496: return;
! 497: }
! 498: norm1c(q1,&a); norm1c(q2,&b);
! 499: n = DEG(DC(q1)); m = DEG(DC(q2));
! 500: pwrq(a,m,&t); pwrq(b,n,&s); mulq(t,s,&u);
! 501: factorial(QTOS(n)+QTOS(m),&t);
! 502: mulq(u,t,&s); addq(s,s,&f);
! 503: for ( index = 0, q = ONE, c = 0; cmpq(f,q) >= 0; ) {
! 504: mod = lprime[index++];
! 505: if ( !mod )
! 506: error("sqfrum : lprime[] exhausted.");
! 507: ptomp(mod,LC(q1),&tg);
! 508: if ( !tg )
! 509: continue;
! 510: ptomp(mod,LC(q2),&tg);
! 511: if ( !tg )
! 512: continue;
! 513: ptomp(mod,q1,&tg1); ptomp(mod,q2,&tg2);
! 514: res_detmp(nvl,mod,v,tg1,tg2,&resg);
! 515: chnremp(nvl,mod,c,q,resg,&c1); c = c1;
! 516: STOQ(mod,t); mulq(q,t,&s); q = s;
! 517: }
! 518: *pr = c;
! 519: }
! 520:
! 521: void res_detmp(vl,mod,v,p1,p2,dp)
! 522: VL vl;
! 523: int mod;
! 524: V v;
! 525: P p1,p2;
! 526: P *dp;
! 527: {
! 528: int n1,n2,n,sgn;
! 529: int i,j,k;
! 530: P mjj,mij,t,s,u,d;
! 531: P **mat;
! 532: P *mi,*mj;
! 533: DCP dc;
! 534:
! 535: n1 = UDEG(p1); n2 = UDEG(p2); n = n1+n2;
! 536: mat = (P **)almat_pointer(n,n);
! 537: for ( dc = DC(p1); dc; dc = NEXT(dc) )
! 538: mat[0][n1-QTOS(DEG(dc))] = COEF(dc);
! 539: for ( i = 1; i < n2; i++ )
! 540: for ( j = 0; j <= n1; j++ )
! 541: mat[i][i+j] = mat[0][j];
! 542: for ( dc = DC(p2); dc; dc = NEXT(dc) )
! 543: mat[n2][n2-QTOS(DEG(dc))] = COEF(dc);
! 544: for ( i = 1; i < n1; i++ )
! 545: for ( j = 0; j <= n2; j++ )
! 546: mat[i+n2][i+j] = mat[n2][j];
! 547: for ( j = 0, d = (P)ONEM, sgn = 1; j < n; j++ ) {
! 548: for ( i = j; (i < n) && !mat[i][j]; i++ );
! 549: if ( i == n ) {
! 550: *dp = 0; return;
! 551: }
! 552: for ( k = i; k < n; k++ )
! 553: if ( mat[k][j] && (nmonop(mat[k][j]) < nmonop(mat[i][j]) ) )
! 554: i = k;
! 555: if ( j != i ) {
! 556: mj = mat[j]; mat[j] = mat[i]; mat[i] = mj; sgn = -sgn;
! 557: }
! 558: for ( i = j + 1, mj = mat[j], mjj = mj[j]; i < n; i++ )
! 559: for ( k = j + 1, mi = mat[i], mij = mi[j]; k < n; k++ ) {
! 560: mulmp(vl,mod,mi[k],mjj,&t); mulmp(vl,mod,mj[k],mij,&s);
! 561: submp(vl,mod,t,s,&u); divsmp(vl,mod,u,d,&mi[k]);
! 562: }
! 563: d = mjj;
! 564: }
! 565: if ( sgn > 0 )
! 566: *dp = d;
! 567: else
! 568: chsgnmp(mod,d,dp);
! 569: }
! 570:
! 571: #if 0
! 572: showmat(vl,mat,n)
! 573: VL vl;
! 574: P **mat;
! 575: int n;
! 576: {
! 577: int i,j;
! 578: P t;
! 579:
! 580: for ( i = 0; i < n; i++ ) {
! 581: for ( j = 0; j < n; j++ ) {
! 582: mptop(mat[i][j],&t); printp(vl,t); fprintf(out," ");
! 583: }
! 584: fprintf(out,"\n");
! 585: }
! 586: fflush(out);
! 587: }
! 588:
! 589: showmp(vl,p)
! 590: VL vl;
! 591: P p;
! 592: {
! 593: P t;
! 594:
! 595: mptop(p,&t); printp(vl,t); fprintf(out,"\n");
! 596: }
! 597: #endif
! 598:
! 599: void premp(vl,p1,p2,pr)
! 600: VL vl;
! 601: P p1,p2,*pr;
! 602: {
! 603: P m,m1,m2;
! 604: P *pw;
! 605: DCP dc;
! 606: V v1,v2;
! 607: register int i,j;
! 608: int n1,n2,d;
! 609:
! 610: if ( NUM(p1) )
! 611: if ( NUM(p2) )
! 612: *pr = 0;
! 613: else
! 614: *pr = p1;
! 615: else if ( NUM(p2) )
! 616: *pr = 0;
! 617: else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
! 618: n1 = deg(v1,p1); n2 = deg(v1,p2);
! 619: pw = (P *)ALLOCA((n1+1)*sizeof(P));
! 620: bzero((char *)pw,(int)((n1+1)*sizeof(P)));
! 621:
! 622: for ( dc = DC(p1); dc; dc = NEXT(dc) )
! 623: pw[QTOS(DEG(dc))] = COEF(dc);
! 624:
! 625: for ( i = n1; i >= n2; i-- ) {
! 626: if ( pw[i] ) {
! 627: m = pw[i];
! 628: for ( j = i; j >= 0; j-- ) {
! 629: mulp(vl,pw[j],LC(p2),&m1); pw[j] = m1;
! 630: }
! 631:
! 632: for ( dc = DC(p2), d = i - n2; dc; dc = NEXT(dc) ) {
! 633: mulp(vl,COEF(dc),m,&m1);
! 634: subp(vl,pw[QTOS(DEG(dc))+d],m1,&m2);
! 635: pw[QTOS(DEG(dc))+d] = m2;
! 636: }
! 637: } else
! 638: for ( j = i; j >= 0; j-- )
! 639: if ( pw[j] ) {
! 640: mulp(vl,pw[j],LC(p2),&m1); pw[j] = m1;
! 641: }
! 642: }
! 643: plisttop(pw,v1,n2-1,pr);
! 644: } else {
! 645: while ( v1 != vl->v && v2 != vl->v )
! 646: vl = NEXT(vl);
! 647: if ( v1 == vl->v )
! 648: *pr = 0;
! 649: else
! 650: *pr = p1;
! 651: }
! 652: }
! 653:
! 654: void ptozp0(p,pr)
! 655: P p;
! 656: P *pr;
! 657: {
! 658: Q c;
! 659:
! 660: ptozp(p,1,&c,pr);
! 661: }
! 662:
! 663: void mindegp(vl,p,mvlp,pr)
! 664: VL vl,*mvlp;
! 665: P p,*pr;
! 666: {
! 667: P t;
! 668: VL nvl,tvl,avl;
! 669: V v;
! 670: int n,d;
! 671:
! 672: clctv(vl,p,&nvl); v = nvl->v; d = getdeg(nvl->v,p);
! 673: for ( avl = NEXT(nvl); avl; avl = NEXT(avl) ) {
! 674: n = getdeg(avl->v,p);
! 675: if ( n < d ) {
! 676: v = avl->v; d = n;
! 677: }
! 678: }
! 679: if ( v != nvl->v ) {
! 680: reordvar(nvl,v,&tvl); reorderp(tvl,nvl,p,&t);
! 681: *pr = t; *mvlp = tvl;
! 682: } else {
! 683: *pr = p; *mvlp = nvl;
! 684: }
! 685: }
! 686:
! 687: void maxdegp(vl,p,mvlp,pr)
! 688: VL vl,*mvlp;
! 689: P p,*pr;
! 690: {
! 691: P t;
! 692: VL nvl,tvl,avl;
! 693: V v;
! 694: int n,d;
! 695:
! 696: clctv(vl,p,&nvl); v = nvl->v; d = getdeg(nvl->v,p);
! 697: for ( avl = NEXT(nvl); avl; avl = NEXT(avl) ) {
! 698: n = getdeg(avl->v,p);
! 699: if ( n > d ) {
! 700: v = avl->v; d = n;
! 701: }
! 702: }
! 703: if ( v != nvl->v ) {
! 704: reordvar(nvl,v,&tvl); reorderp(tvl,nvl,p,&t);
! 705: *pr = t; *mvlp = tvl;
! 706: } else {
! 707: *pr = p; *mvlp = nvl;
! 708: }
! 709: }
! 710:
! 711: void min_common_vars_in_coefp(vl,p,mvlp,pr)
! 712: VL vl,*mvlp;
! 713: P p,*pr;
! 714: {
! 715: P u,p0;
! 716: VL tvl,cvl,svl,uvl,avl,vl0;
! 717: int n,n0,d,d0;
! 718: DCP dc;
! 719:
! 720: clctv(vl,p,&tvl); vl = tvl;
! 721: for ( tvl = vl, n0 = 0; tvl; tvl = NEXT(tvl), n0++ );
! 722: for ( avl = vl; avl; avl = NEXT(avl) ) {
! 723: if ( VR(p) != avl->v ) {
! 724: reordvar(vl,avl->v,&uvl); reorderp(uvl,vl,p,&u);
! 725: } else {
! 726: uvl = vl; u = p;
! 727: }
! 728: for ( cvl = NEXT(uvl), dc = DC(u); dc && cvl; dc = NEXT(dc) ) {
! 729: clctv(uvl,COEF(dc),&tvl); intersectv(cvl,tvl,&svl); cvl = svl;
! 730: }
! 731: if ( !cvl ) {
! 732: *mvlp = uvl; *pr = u; return;
! 733: } else {
! 734: for ( tvl = cvl, n = 0; tvl; tvl = NEXT(tvl), n++ );
! 735: if ( n < n0 ) {
! 736: vl0 = uvl; p0 = u; n0 = n; d0 = getdeg(uvl->v,u);
! 737: } else if ( (n == n0) && ((d = getdeg(uvl->v,u)) < d0) ) {
! 738: vl0 = uvl; p0 = u; n0 = n; d0 = d;
! 739: }
! 740: }
! 741: }
! 742: *pr = p0; *mvlp = vl0;
! 743: }
! 744:
! 745: void minlcdegp(vl,p,mvlp,pr)
! 746: VL vl,*mvlp;
! 747: P p,*pr;
! 748: {
! 749: P u,p0;
! 750: VL tvl,uvl,avl,vl0;
! 751: int d,d0;
! 752:
! 753: clctv(vl,p,&tvl); vl = tvl;
! 754: vl0 = vl; p0 = p; d0 = homdeg(COEF(DC(p)));
! 755: for ( avl = NEXT(vl); avl; avl = NEXT(avl) ) {
! 756: reordvar(vl,avl->v,&uvl); reorderp(uvl,vl,p,&u);
! 757: d = homdeg(COEF(DC(u)));
! 758: if ( d < d0 ) {
! 759: vl0 = uvl; p0 = u; d0 = d;
! 760: }
! 761: }
! 762: *pr = p0; *mvlp = vl0;
! 763: }
! 764:
! 765: void sort_by_deg(n,p,pr)
! 766: int n;
! 767: P *p,*pr;
! 768: {
! 769: int j,k,d,k0;
! 770: V v;
! 771:
! 772: for ( j = 0; j < n; j++ ) {
! 773: for ( k0 = k = j, d = deg(v = VR(p[0]),p[j]);
! 774: k < n; k++ )
! 775: if ( deg(v,p[k]) < d ) {
! 776: k0 = k;
! 777: d = deg(v,p[k]);
! 778: }
! 779: pr[j] = p[k0];
! 780: if ( j != k0 )
! 781: p[k0] = p[j];
! 782: }
! 783: }
! 784:
! 785: void sort_by_deg_rev(n,p,pr)
! 786: int n;
! 787: P *p,*pr;
! 788: {
! 789: int j,k,d,k0;
! 790: V v;
! 791:
! 792: for ( j = 0; j < n; j++ ) {
! 793: for ( k0 = k = j, d = deg(v = VR(p[0]),p[j]);
! 794: k < n; k++ )
! 795: if ( deg(v,p[k]) > d ) {
! 796: k0 = k;
! 797: d = deg(v,p[k]);
! 798: }
! 799: pr[j] = p[k0];
! 800: if ( j != k0 )
! 801: p[k0] = p[j];
! 802: }
! 803: }
! 804:
! 805:
! 806: void getmindeg(v,p,dp)
! 807: V v;
! 808: P p;
! 809: Q *dp;
! 810: {
! 811: Q dt,d;
! 812: DCP dc;
! 813:
! 814: if ( !p || NUM(p) )
! 815: *dp = 0;
! 816: else if ( v == VR(p) ) {
! 817: for ( dc = DC(p); NEXT(dc); dc = NEXT(dc) );
! 818: *dp = DEG(dc);
! 819: } else {
! 820: dc = DC(p);
! 821: getmindeg(v,COEF(dc),&d);
! 822: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
! 823: getmindeg(v,COEF(dc),&dt);
! 824: if ( cmpq(dt,d) < 0 )
! 825: d = dt;
! 826: }
! 827: *dp = d;
! 828: }
! 829: }
! 830:
! 831: void minchdegp(vl,p,mvlp,pr)
! 832: VL vl,*mvlp;
! 833: P p,*pr;
! 834: {
! 835: P t;
! 836: VL tvl,nvl,avl;
! 837: int n,d,z;
! 838: V v;
! 839:
! 840: if ( NUM(p) ) {
! 841: *mvlp = vl;
! 842: *pr = p;
! 843: return;
! 844: }
! 845: clctv(vl,p,&nvl);
! 846: v = nvl->v;
! 847: d = getchomdeg(v,p) + getlchomdeg(v,p,&z);
! 848: for ( avl = NEXT(nvl); avl; avl = NEXT(avl) ) {
! 849: n = getchomdeg(avl->v,p) + getlchomdeg(avl->v,p,&z);
! 850: if ( n < d ) {
! 851: v = avl->v; d = n;
! 852: }
! 853: }
! 854: if ( v != nvl->v ) {
! 855: reordvar(nvl,v,&tvl); reorderp(tvl,nvl,p,&t);
! 856: *pr = t; *mvlp = tvl;
! 857: } else {
! 858: *pr = p; *mvlp = nvl;
! 859: }
! 860: }
! 861:
! 862: int getchomdeg(v,p)
! 863: V v;
! 864: P p;
! 865: {
! 866: int m,m1;
! 867: DCP dc;
! 868:
! 869: if ( !p || NUM(p) )
! 870: return ( 0 );
! 871: else if ( VR(p) == v )
! 872: return ( dbound(v,p) );
! 873: else {
! 874: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
! 875: m1 = getchomdeg(v,COEF(dc))+QTOS(DEG(dc));
! 876: m = MAX(m,m1);
! 877: }
! 878: return ( m );
! 879: }
! 880: }
! 881:
! 882: int getlchomdeg(v,p,d)
! 883: V v;
! 884: P p;
! 885: int *d;
! 886: {
! 887: int m0,m1,d0,d1;
! 888: DCP dc;
! 889:
! 890: if ( !p || NUM(p) ) {
! 891: *d = 0;
! 892: return ( 0 );
! 893: } else if ( VR(p) == v ) {
! 894: *d = QTOS(DEG(DC(p)));
! 895: return ( homdeg(LC(p)) );
! 896: } else {
! 897: for ( dc = DC(p), m0 = 0, d0 = 0; dc; dc = NEXT(dc) ) {
! 898: m1 = getlchomdeg(v,COEF(dc),&d1)+QTOS(DEG(dc));
! 899: if ( d1 > d0 ) {
! 900: m0 = m1;
! 901: d0 = d1;
! 902: } else if ( d1 == d0 )
! 903: m0 = MAX(m0,m1);
! 904: }
! 905: *d = d0;
! 906: return ( m0 );
! 907: }
! 908: }
! 909:
! 910: int nmonop(p)
! 911: P p;
! 912: {
! 913: int s;
! 914: DCP dc;
! 915:
! 916: if ( !p )
! 917: return 0;
! 918: else
! 919: switch ( OID(p) ) {
! 920: case O_N:
! 921: return 1; break;
! 922: case O_P:
! 923: for ( dc = DC((P)p), s = 0; dc; dc = NEXT(dc) )
! 924: s += nmonop(COEF(dc));
! 925: return s; break;
! 926: default:
! 927: return 0;
! 928: }
! 929: }
! 930:
! 931: int qpcheck(p)
! 932: Obj p;
! 933: {
! 934: DCP dc;
! 935:
! 936: if ( !p )
! 937: return 1;
! 938: else
! 939: switch ( OID(p) ) {
! 940: case O_N:
! 941: return RATN((Num)p)?1:0;
! 942: case O_P:
! 943: for ( dc = DC((P)p); dc; dc = NEXT(dc) )
! 944: if ( !qpcheck((Obj)COEF(dc)) )
! 945: return 0;
! 946: return 1;
! 947: default:
! 948: return 0;
! 949: }
! 950: }
! 951:
! 952: /* check if p is univariate and all coeffs are INT or LM */
! 953:
! 954: int uzpcheck(p)
! 955: Obj p;
! 956: {
! 957: DCP dc;
! 958: P c;
! 959:
! 960: if ( !p )
! 961: return 1;
! 962: else
! 963: switch ( OID(p) ) {
! 964: case O_N:
! 965: return (RATN((Num)p)&&INT(p))||(NID((Num)p)==N_LM);
! 966: case O_P:
! 967: for ( dc = DC((P)p); dc; dc = NEXT(dc) ) {
! 968: c = COEF(dc);
! 969: if ( !NUM(c) || !uzpcheck(c) )
! 970: return 0;
! 971: }
! 972: return 1;
! 973: default:
! 974: return 0;
! 975: }
! 976: }
! 977:
! 978: int p_mag(p)
! 979: P p;
! 980: {
! 981: int s;
! 982: DCP dc;
! 983:
! 984: if ( !p )
! 985: return 0;
! 986: else if ( OID(p) == O_N )
! 987: return n_bits(NM((Q)p))+(INT((Q)p)?0:n_bits(DN((Q)p)));
! 988: else {
! 989: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) )
! 990: s += p_mag(COEF(dc));
! 991: return s;
! 992: }
! 993: }