![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to PM.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 PM.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 ! noro 1: /* $OpenXM: OpenXM/src/asir99/engine/PM.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
! 2: #ifndef MODULAR
! 3: #define MODULAR
! 4: #endif
! 5:
! 6: #include "b.h"
! 7: #include "ca.h"
! 8:
! 9: void D_ADDP(vl,p1,p2,pr)
! 10: VL vl;
! 11: P p1,p2,*pr;
! 12: {
! 13: register DCP dc1,dc2,dcr0,dcr;
! 14: V v1,v2;
! 15: P t;
! 16:
! 17: if ( !p1 )
! 18: *pr = p2;
! 19: else if ( !p2 )
! 20: *pr = p1;
! 21: else if ( NUM(p1) )
! 22: if ( NUM(p2) )
! 23: ADDNUM(p1,p2,pr);
! 24: else
! 25: ADDPQ(p2,p1,pr);
! 26: else if ( NUM(p2) )
! 27: ADDPQ(p1,p2,pr);
! 28: else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
! 29: for ( dc1 = DC(p1), dc2 = DC(p2), dcr0 = 0; dc1 && dc2; )
! 30: switch ( cmpq(DEG(dc1),DEG(dc2)) ) {
! 31: case 0:
! 32: ADDP(vl,COEF(dc1),COEF(dc2),&t);
! 33: if ( t ) {
! 34: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc1); COEF(dcr) = t;
! 35: }
! 36: dc1 = NEXT(dc1); dc2 = NEXT(dc2); break;
! 37: case 1:
! 38: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc1); COEF(dcr) = COEF(dc1);
! 39: dc1 = NEXT(dc1); break;
! 40: case -1:
! 41: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc2); COEF(dcr) = COEF(dc2);
! 42: dc2 = NEXT(dc2); break;
! 43: }
! 44: if ( !dcr0 )
! 45: if ( dc1 )
! 46: dcr0 = dc1;
! 47: else if ( dc2 )
! 48: dcr0 = dc2;
! 49: else {
! 50: *pr = 0;
! 51: return;
! 52: }
! 53: else
! 54: if ( dc1 )
! 55: NEXT(dcr) = dc1;
! 56: else if ( dc2 )
! 57: NEXT(dcr) = dc2;
! 58: else
! 59: NEXT(dcr) = 0;
! 60: MKP(v1,dcr0,*pr);
! 61: } else {
! 62: while ( v1 != VR(vl) && v2 != VR(vl) )
! 63: vl = NEXT(vl);
! 64: if ( v1 == VR(vl) )
! 65: ADDPTOC(vl,p1,p2,pr);
! 66: else
! 67: ADDPTOC(vl,p2,p1,pr);
! 68: }
! 69: }
! 70:
! 71: void D_ADDPQ(p,q,pr)
! 72: P p,q,*pr;
! 73: {
! 74: DCP dc,dcr,dcr0;
! 75: P t;
! 76:
! 77: if ( NUM(p) )
! 78: ADDNUM(p,q,pr);
! 79: else {
! 80: for ( dcr0 = 0, dc = DC(p); dc && DEG(dc); dc = NEXT(dc) ) {
! 81: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = COEF(dc);
! 82: }
! 83: if ( !dc ) {
! 84: NEWDC(NEXT(dcr)); dcr = NEXT(dcr); DEG(dcr) = 0; COEF(dcr) = q;
! 85: } else {
! 86: ADDPQ(COEF(dc),q,&t);
! 87: if ( t ) {
! 88: NEWDC(NEXT(dcr)); dcr = NEXT(dcr); DEG(dcr) = 0; COEF(dcr) = t;
! 89: }
! 90: }
! 91: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
! 92: }
! 93: }
! 94:
! 95: void D_ADDPTOC(vl,p,c,pr)
! 96: VL vl;
! 97: P p,c,*pr;
! 98: {
! 99: DCP dc,dcr,dcr0;
! 100: P t;
! 101:
! 102: for ( dcr0 = 0, dc = DC(p); dc && DEG(dc); dc = NEXT(dc) ) {
! 103: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = COEF(dc);
! 104: }
! 105: if ( !dc ) {
! 106: NEWDC(NEXT(dcr)); dcr = NEXT(dcr); DEG(dcr) = 0; COEF(dcr) = c;
! 107: } else {
! 108: ADDP(vl,COEF(dc),c,&t);
! 109: if ( t ) {
! 110: NEWDC(NEXT(dcr)); dcr = NEXT(dcr); DEG(dcr) = 0; COEF(dcr) = t;
! 111: }
! 112: }
! 113: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
! 114: }
! 115:
! 116: void D_SUBP(vl,p1,p2,pr)
! 117: VL vl;
! 118: P p1,p2,*pr;
! 119: {
! 120: P t;
! 121:
! 122: if ( !p2 )
! 123: *pr = p1;
! 124: else {
! 125: CHSGNP(p2,&t); ADDP(vl,p1,t,pr);
! 126: }
! 127: }
! 128:
! 129: void D_MULP(vl,p1,p2,pr)
! 130: VL vl;
! 131: P p1,p2,*pr;
! 132: {
! 133: register DCP dc,dct,dcr,dcr0;
! 134: V v1,v2;
! 135: P t,s,u;
! 136: int n1,n2;
! 137:
! 138: if ( !p1 || !p2 ) *pr = 0;
! 139: else if ( NUM(p1) )
! 140: MULPQ(p2,p1,pr);
! 141: else if ( NUM(p2) )
! 142: MULPQ(p1,p2,pr);
! 143: else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
! 144: for ( dc = DC(p1), n1 = 0; dc; dc = NEXT(dc), n1++ );
! 145: for ( dc = DC(p2), n2 = 0; dc; dc = NEXT(dc), n2++ );
! 146: if ( n1 > n2 )
! 147: for ( dc = DC(p2), s = 0; dc; dc = NEXT(dc) ) {
! 148: for ( dcr0 = 0, dct = DC(p1); dct; dct = NEXT(dct) ) {
! 149: NEXTDC(dcr0,dcr); MULP(vl,COEF(dct),COEF(dc),&COEF(dcr));
! 150: addq(DEG(dct),DEG(dc),&DEG(dcr));
! 151: }
! 152: NEXT(dcr) = 0; MKP(v1,dcr0,t);
! 153: ADDP(vl,s,t,&u); s = u; t = u = 0;
! 154: }
! 155: else
! 156: for ( dc = DC(p1), s = 0; dc; dc = NEXT(dc) ) {
! 157: for ( dcr0 = 0, dct = DC(p2); dct; dct = NEXT(dct) ) {
! 158: NEXTDC(dcr0,dcr); MULP(vl,COEF(dct),COEF(dc),&COEF(dcr));
! 159: addq(DEG(dct),DEG(dc),&DEG(dcr));
! 160: }
! 161: NEXT(dcr) = 0; MKP(v1,dcr0,t);
! 162: ADDP(vl,s,t,&u); s = u; t = u = 0;
! 163: }
! 164: *pr = s;
! 165: } else {
! 166: while ( v1 != VR(vl) && v2 != VR(vl) )
! 167: vl = NEXT(vl);
! 168: if ( v1 == VR(vl) )
! 169: MULPC(vl,p1,p2,pr);
! 170: else
! 171: MULPC(vl,p2,p1,pr);
! 172: }
! 173: }
! 174:
! 175: void D_MULPQ(p,q,pr)
! 176: P p,q,*pr;
! 177: {
! 178: DCP dc,dcr,dcr0;
! 179:
! 180: if (!p || !q)
! 181: *pr = 0;
! 182: else if ( Uniq(q) )
! 183: *pr = p;
! 184: else if ( NUM(p) )
! 185: MULNUM(p,q,pr);
! 186: else {
! 187: for ( dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
! 188: NEXTDC(dcr0,dcr); MULPQ(COEF(dc),q,&COEF(dcr)); DEG(dcr) = DEG(dc);
! 189: }
! 190: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
! 191: }
! 192: }
! 193:
! 194: void D_MULPC(vl,p,c,pr)
! 195: VL vl;
! 196: P p,c,*pr;
! 197: {
! 198: DCP dc,dcr,dcr0;
! 199:
! 200: if ( NUM(c) )
! 201: MULPQ(p,c,pr);
! 202: else {
! 203: for ( dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
! 204: NEXTDC(dcr0,dcr); MULP(vl,COEF(dc),c,&COEF(dcr)); DEG(dcr) = DEG(dc);
! 205: }
! 206: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
! 207: }
! 208: }
! 209:
! 210: void D_PWRP(vl,p,q,pr)
! 211: VL vl;
! 212: P p,*pr;
! 213: Q q;
! 214: {
! 215: register DCP dc,dcr;
! 216: int n,i;
! 217: P *x,*y;
! 218: P t,s,u;
! 219: DCP dct;
! 220: P *pt;
! 221:
! 222: if ( !q ) {
! 223: *pr = (P)One;
! 224: } else if ( !p )
! 225: *pr = 0;
! 226: else if ( UNIQ(q) )
! 227: *pr = p;
! 228: else if ( NUM(p) )
! 229: PWRNUM(p,q,pr);
! 230: else {
! 231: dc = DC(p);
! 232: if ( !NEXT(dc) ) {
! 233: NEWDC(dcr);
! 234: PWRP(vl,COEF(dc),q,&COEF(dcr)); mulq(DEG(dc),q,&DEG(dcr));
! 235: NEXT(dcr) = 0; MKP(VR(p),dcr,*pr);
! 236: } else if ( !INT(q) ) {
! 237: error("pwrp: can't calculate fractional power."); *pr = 0;
! 238: } else if ( PL(NM(q)) == 1 ) {
! 239: n = QTOS(q); x = (P *)ALLOCA((n+1)*sizeof(pointer));
! 240: NEWDC(dct); DEG(dct) = DEG(dc); COEF(dct) = COEF(dc);
! 241: NEXT(dct) = 0; MKP(VR(p),dct,t);
! 242: for ( i = 0, u = (P)One; i < n; i++ ) {
! 243: x[i] = u; MULP(vl,u,t,&s); u = s;
! 244: }
! 245: x[n] = u; y = (P *)ALLOCA((n+1)*sizeof(pointer));
! 246: if ( DEG(NEXT(dc)) ) {
! 247: dct = NEXT(dc); MKP(VR(p),dct,t);
! 248: } else
! 249: t = COEF(NEXT(dc));
! 250: for ( i = 0, u = (P)One; i < n; i++ ) {
! 251: y[i] = u; MULP(vl,u,t,&s); u = s;
! 252: }
! 253: y[n] = u;
! 254: pt = (P *)ALLOCA((n+1)*sizeof(pointer)); MKBC(n,pt);
! 255: for ( i = 0, u = 0; i <= n; i++ ) {
! 256: MULP(vl,x[i],y[n-i],&t); MULP(vl,t,pt[i],&s);
! 257: ADDP(vl,u,s,&t); u = t;
! 258: }
! 259: *pr = u;
! 260: } else {
! 261: error("exponent too big");
! 262: *pr = 0;
! 263: }
! 264: }
! 265: }
! 266:
! 267: void D_CHSGNP(p,pr)
! 268: P p,*pr;
! 269: {
! 270: register DCP dc,dcr,dcr0;
! 271: P t;
! 272:
! 273: if ( !p )
! 274: *pr = NULL;
! 275: else if ( NUM(p) ) {
! 276: #if defined(THINK_C) || defined(_PA_RISC1_1) || defined(__alpha) || defined(mips)
! 277: #ifdef FBASE
! 278: chsgnnum((Num)p,(Num *)pr);
! 279: #else
! 280: MQ mq;
! 281:
! 282: NEWMQ(mq); CONT(mq)=mod-CONT((MQ)p); *pr = (P)mq;
! 283: #endif
! 284: #else
! 285: CHSGNNUM(p,t); *pr = t;
! 286: #endif
! 287: } else {
! 288: for ( dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
! 289: NEXTDC(dcr0,dcr); CHSGNP(COEF(dc),&COEF(dcr)); DEG(dcr) = DEG(dc);
! 290: }
! 291: NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);
! 292: }
! 293: }
! 294:
! 295:
! 296: #ifdef FBASE
! 297: void ptozp(p,sgn,c,pr)
! 298: P p;
! 299: int sgn;
! 300: Q *c;
! 301: P *pr;
! 302: {
! 303: N nm,dn;
! 304:
! 305: if ( !p ) {
! 306: *c = 0; *pr = 0;
! 307: } else {
! 308: lgp(p,&nm,&dn);
! 309: if ( UNIN(dn) )
! 310: NTOQ(nm,sgn,*c);
! 311: else
! 312: NDTOQ(nm,dn,sgn,*c);
! 313: divcp(p,*c,pr);
! 314: }
! 315: }
! 316:
! 317: void lgp(p,g,l)
! 318: P p;
! 319: N *g,*l;
! 320: {
! 321: N g1,g2,l1,l2,l3,l4,tmp;
! 322: DCP dc;
! 323:
! 324: if ( NUM(p) ) {
! 325: *g = NM((Q)p);
! 326: if ( INT((Q)p) )
! 327: *l = ONEN;
! 328: else
! 329: *l = DN((Q)p);
! 330: } else {
! 331: dc = DC(p); lgp(COEF(dc),g,l);
! 332: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
! 333: lgp(COEF(dc),&g1,&l1); gcdn(*g,g1,&g2); *g = g2;
! 334: gcdn(*l,l1,&l2); kmuln(*l,l1,&l3); divn(l3,l2,&l4,&tmp); *l = l4;
! 335: }
! 336: }
! 337: }
! 338:
! 339: void divcp(f,q,rp)
! 340: P f;
! 341: Q q;
! 342: P *rp;
! 343: {
! 344: DCP dc,dcr,dcr0;
! 345:
! 346: if ( !f )
! 347: *rp = 0;
! 348: else if ( NUM(f) )
! 349: DIVNUM(f,q,rp);
! 350: else {
! 351: for ( dc = DC(f), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 352: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc);
! 353: divcp(COEF(dc),q,&COEF(dcr));
! 354: }
! 355: NEXT(dcr) = 0; MKP(VR(f),dcr0,*rp);
! 356: }
! 357: }
! 358:
! 359: void diffp(vl,p,v,r)
! 360: VL vl;
! 361: P p;
! 362: V v;
! 363: P *r;
! 364: {
! 365: P t;
! 366: DCP dc,dcr,dcr0;
! 367:
! 368: if ( !p || NUM(p) )
! 369: *r = 0;
! 370: else {
! 371: if ( v == VR(p) ) {
! 372: for ( dc = DC(p), dcr0 = 0;
! 373: dc && DEG(dc); dc = NEXT(dc) ) {
! 374: NEXTDC(dcr0,dcr); SUBQ(DEG(dc),ONE,&DEG(dcr));
! 375: MULPQ(COEF(dc),(P)DEG(dc),&COEF(dcr));
! 376: }
! 377: NEXT(dcr) = 0; MKP(v,dcr0,*r);
! 378: } else {
! 379: for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 380: diffp(vl,COEF(dc),v,&t);
! 381: if ( t ) {
! 382: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t;
! 383: }
! 384: }
! 385: if ( !dcr0 )
! 386: *r = 0;
! 387: else {
! 388: NEXT(dcr) = 0; MKP(VR(p),dcr0,*r);
! 389: }
! 390: }
! 391: }
! 392: }
! 393:
! 394: void coefp(p,d,pr)
! 395: P p;
! 396: int d;
! 397: P *pr;
! 398: {
! 399: DCP dc;
! 400: int sgn;
! 401: Q dq;
! 402:
! 403: if ( NUM(p) )
! 404: if ( d == 0 )
! 405: *pr = p;
! 406: else
! 407: *pr = 0;
! 408: else {
! 409: for ( STOQ(d,dq), dc = DC(p); dc; dc = NEXT(dc) )
! 410: if ( (sgn = cmpq(DEG(dc),dq)) > 0 )
! 411: continue;
! 412: else if ( sgn == 0 ) {
! 413: *pr = COEF(dc);
! 414: return;
! 415: } else {
! 416: *pr = 0;
! 417: break;
! 418: }
! 419: *pr = 0;
! 420: }
! 421: }
! 422:
! 423: int compp(vl,p1,p2)
! 424: VL vl;
! 425: P p1,p2;
! 426: {
! 427: DCP dc1,dc2;
! 428: V v1,v2;
! 429:
! 430: if ( !p1 )
! 431: return p2 ? -1 : 0;
! 432: else if ( !p2 )
! 433: return 1;
! 434: else if ( NUM(p1) )
! 435: return NUM(p2) ? (*cmpnumt[MAX(NID(p1),NID(p2))])(p1,p2) : -1;
! 436: else if ( NUM(p2) )
! 437: return 1;
! 438: if ( (v1 = VR(p1)) == (v2 = VR(p2)) ) {
! 439: for ( dc1 = DC(p1), dc2 = DC(p2);
! 440: dc1 && dc2; dc1 = NEXT(dc1), dc2 = NEXT(dc2) )
! 441: switch ( cmpq(DEG(dc1),DEG(dc2)) ) {
! 442: case 1:
! 443: return 1; break;
! 444: case -1:
! 445: return -1; break;
! 446: default:
! 447: switch ( compp(vl,COEF(dc1),COEF(dc2)) ) {
! 448: case 1:
! 449: return 1; break;
! 450: case -1:
! 451: return -1; break;
! 452: default:
! 453: break;
! 454: }
! 455: break;
! 456: }
! 457: return dc1 ? 1 : (dc2 ? -1 : 0);
! 458: } else {
! 459: for ( ; v1 != VR(vl) && v2 != VR(vl); vl = NEXT(vl) );
! 460: return v1 == VR(vl) ? 1 : -1;
! 461: }
! 462: }
! 463:
! 464: void csump(vl,p,c)
! 465: VL vl;
! 466: P p;
! 467: Q *c;
! 468: {
! 469: DCP dc;
! 470: Q s,s1,s2;
! 471:
! 472: if ( !p || NUM(p) )
! 473: *c = (Q)p;
! 474: else {
! 475: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {
! 476: csump(vl,COEF(dc),&s1); addq(s,s1,&s2); s = s2;
! 477: }
! 478: *c = s;
! 479: }
! 480: }
! 481:
! 482: void degp(v,p,d)
! 483: V v;
! 484: P p;
! 485: Q *d;
! 486: {
! 487: DCP dc;
! 488: Q m,m1;
! 489:
! 490: if ( !p || NUM(p) )
! 491: *d = 0;
! 492: else if ( v == VR(p) )
! 493: *d = DEG(DC(p));
! 494: else {
! 495: for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {
! 496: degp(v,COEF(dc),&m1);
! 497: if ( cmpq(m,m1) < 0 )
! 498: m = m1;
! 499: }
! 500: *d = m;
! 501: }
! 502: }
! 503: #endif