[BACK]Return to PM.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / engine

Diff for /OpenXM_contrib2/asir2000/engine/PM.c between version 1.5 and 1.6

version 1.5, 2001/06/15 07:56:03 version 1.6, 2004/08/18 00:17:02
Line 45 
Line 45 
  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,   * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.   * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
  *   *
  * $OpenXM: OpenXM_contrib2/asir2000/engine/PM.c,v 1.4 2000/12/05 01:24:52 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/engine/PM.c,v 1.5 2001/06/15 07:56:03 noro Exp $
 */  */
 #ifndef MODULAR  #ifndef MODULAR
 #define MODULAR  #define MODULAR
Line 54 
Line 54 
 #include "b.h"  #include "b.h"
 #include "ca.h"  #include "ca.h"
   
 void D_ADDP(vl,p1,p2,pr)  #include "poly.c"
 VL vl;  
 P p1,p2,*pr;  
 {  
         register DCP dc1,dc2,dcr0,dcr;  
         V v1,v2;  
         P t;  
   
         if ( !p1 )  
                 *pr = p2;  
         else if ( !p2 )  
                 *pr = p1;  
         else if ( NUM(p1) )  
                 if ( NUM(p2) )  
                         ADDNUM(p1,p2,pr);  
                 else  
                         ADDPQ(p2,p1,pr);  
         else if ( NUM(p2) )  
                 ADDPQ(p1,p2,pr);  
         else if ( ( v1 = VR(p1) ) ==  ( v2 = VR(p2) ) ) {  
                 for ( dc1 = DC(p1), dc2 = DC(p2), dcr0 = 0; dc1 && dc2; )  
                         switch ( cmpq(DEG(dc1),DEG(dc2)) ) {  
                                 case 0:  
                                         ADDP(vl,COEF(dc1),COEF(dc2),&t);  
                                         if ( t )  {  
                                                 NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc1); COEF(dcr) = t;  
                                         }  
                                         dc1 = NEXT(dc1); dc2 = NEXT(dc2); break;  
                                 case 1:  
                                         NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc1); COEF(dcr) = COEF(dc1);  
                                         dc1 = NEXT(dc1); break;  
                                 case -1:  
                                         NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc2); COEF(dcr) = COEF(dc2);  
                                         dc2 = NEXT(dc2); break;  
                         }  
                 if ( !dcr0 )  
                         if ( dc1 )  
                                 dcr0 = dc1;  
                         else if ( dc2 )  
                                 dcr0 = dc2;  
                         else {  
                                 *pr = 0;  
                                 return;  
                         }  
                 else  
                         if ( dc1 )  
                                 NEXT(dcr) = dc1;  
                         else if ( dc2 )  
                                 NEXT(dcr) = dc2;  
                         else  
                                 NEXT(dcr) = 0;  
                 MKP(v1,dcr0,*pr);  
         } else {  
                 while ( v1 != VR(vl) && v2 != VR(vl) )  
                         vl = NEXT(vl);  
                 if ( v1 == VR(vl) )  
                         ADDPTOC(vl,p1,p2,pr);  
                 else  
                         ADDPTOC(vl,p2,p1,pr);  
         }  
 }  
   
 void D_ADDPQ(p,q,pr)  
 P p,q,*pr;  
 {  
         DCP dc,dcr,dcr0;  
         P t;  
   
         if ( NUM(p) )  
                 ADDNUM(p,q,pr);  
         else {  
                 for ( dcr0 = 0, dc = DC(p); dc && DEG(dc); dc = NEXT(dc) ) {  
                         NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = COEF(dc);  
                 }  
                 if ( !dc ) {  
                         NEWDC(NEXT(dcr)); dcr = NEXT(dcr); DEG(dcr) = 0; COEF(dcr) = q;  
                 } else {  
                         ADDPQ(COEF(dc),q,&t);  
                         if ( t ) {  
                                 NEWDC(NEXT(dcr)); dcr = NEXT(dcr); DEG(dcr) = 0; COEF(dcr) = t;  
                         }  
                 }  
                 NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);  
         }  
 }  
   
 void D_ADDPTOC(vl,p,c,pr)  
 VL vl;  
 P p,c,*pr;  
 {  
         DCP dc,dcr,dcr0;  
         P t;  
   
         for ( dcr0 = 0, dc = DC(p); dc && DEG(dc); dc = NEXT(dc) ) {  
                 NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = COEF(dc);  
         }  
         if ( !dc ) {  
                 NEWDC(NEXT(dcr)); dcr = NEXT(dcr); DEG(dcr) = 0; COEF(dcr) = c;  
         } else {  
                 ADDP(vl,COEF(dc),c,&t);  
                 if ( t ) {  
                         NEWDC(NEXT(dcr)); dcr = NEXT(dcr); DEG(dcr) = 0; COEF(dcr) = t;  
                 }  
         }  
         NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);  
 }  
   
 void D_SUBP(vl,p1,p2,pr)  
 VL vl;  
 P p1,p2,*pr;  
 {  
         P t;  
   
         if ( !p2 )  
                 *pr = p1;  
         else {  
                 CHSGNP(p2,&t); ADDP(vl,p1,t,pr);  
         }  
 }  
   
 void D_MULP(vl,p1,p2,pr)  
 VL vl;  
 P p1,p2,*pr;  
 {  
         register DCP dc,dct,dcr,dcr0;  
         V v1,v2;  
         P t,s,u;  
         int n1,n2;  
   
         if ( !p1 || !p2 ) *pr = 0;  
         else if ( NUM(p1) )  
                 MULPQ(p2,p1,pr);  
         else if ( NUM(p2) )  
                 MULPQ(p1,p2,pr);  
         else if ( ( v1 = VR(p1) ) ==  ( v2 = VR(p2) ) ) {  
                 for ( dc = DC(p1), n1 = 0; dc; dc = NEXT(dc), n1++ );  
                 for ( dc = DC(p2), n2 = 0; dc; dc = NEXT(dc), n2++ );  
                 if ( n1 > n2 )  
                         for ( dc = DC(p2), s = 0; dc; dc = NEXT(dc) ) {  
                                 for ( dcr0 = 0, dct = DC(p1); dct; dct = NEXT(dct) ) {  
                                         NEXTDC(dcr0,dcr); MULP(vl,COEF(dct),COEF(dc),&COEF(dcr));  
                                         addq(DEG(dct),DEG(dc),&DEG(dcr));  
                                 }  
                                 NEXT(dcr) = 0; MKP(v1,dcr0,t);  
                                 ADDP(vl,s,t,&u); s = u; t = u = 0;  
                         }  
                 else  
                         for ( dc = DC(p1), s = 0; dc; dc = NEXT(dc) ) {  
                                 for ( dcr0 = 0, dct = DC(p2); dct; dct = NEXT(dct) ) {  
                                         NEXTDC(dcr0,dcr); MULP(vl,COEF(dct),COEF(dc),&COEF(dcr));  
                                         addq(DEG(dct),DEG(dc),&DEG(dcr));  
                                 }  
                                 NEXT(dcr) = 0; MKP(v1,dcr0,t);  
                                 ADDP(vl,s,t,&u); s = u; t = u = 0;  
                         }  
                 *pr = s;  
         } else {  
                 while ( v1 != VR(vl) && v2 != VR(vl) )  
                         vl = NEXT(vl);  
                 if ( v1 == VR(vl) )  
                         MULPC(vl,p1,p2,pr);  
                 else  
                         MULPC(vl,p2,p1,pr);  
         }  
 }  
   
 void D_MULPQ(p,q,pr)  
 P p,q,*pr;  
 {  
         DCP dc,dcr,dcr0;  
   
         if (!p || !q)  
                 *pr = 0;  
         else if ( Uniq(q) )  
                 *pr = p;  
         else if ( NUM(p) )  
                 MULNUM(p,q,pr);  
         else {  
                 for ( dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {  
                         NEXTDC(dcr0,dcr); MULPQ(COEF(dc),q,&COEF(dcr)); DEG(dcr) = DEG(dc);  
                 }  
                 NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);  
         }  
 }  
   
 void D_MULPC(vl,p,c,pr)  
 VL vl;  
 P p,c,*pr;  
 {  
         DCP dc,dcr,dcr0;  
   
         if ( NUM(c) )  
                 MULPQ(p,c,pr);  
         else {  
                 for ( dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {  
                         NEXTDC(dcr0,dcr); MULP(vl,COEF(dc),c,&COEF(dcr)); DEG(dcr) = DEG(dc);  
                 }  
                 NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);  
         }  
 }  
   
 void D_PWRP(vl,p,q,pr)  
 VL vl;  
 P p,*pr;  
 Q q;  
 {  
         register DCP dc,dcr;  
         int n,i;  
         P *x,*y;  
         P t,s,u;  
         DCP dct;  
         P *pt;  
   
         if ( !q ) {  
                 *pr = (P)One;  
         } else if ( !p )  
                 *pr = 0;  
         else if ( UNIQ(q) )  
                 *pr = p;  
         else if ( NUM(p) )  
                 PWRNUM(p,q,pr);  
         else {  
                 dc = DC(p);  
                 if ( !NEXT(dc) ) {  
                         NEWDC(dcr);  
                         PWRP(vl,COEF(dc),q,&COEF(dcr)); mulq(DEG(dc),q,&DEG(dcr));  
                         NEXT(dcr) = 0; MKP(VR(p),dcr,*pr);  
                 } else if ( !INT(q) ) {  
                         error("pwrp: can't calculate fractional power."); *pr = 0;  
                 } else if ( PL(NM(q)) == 1 ) {  
                         n = QTOS(q); x = (P *)ALLOCA((n+1)*sizeof(pointer));  
                         NEWDC(dct); DEG(dct) = DEG(dc); COEF(dct) = COEF(dc);  
                         NEXT(dct) = 0; MKP(VR(p),dct,t);  
                         for ( i = 0, u = (P)One; i < n; i++ ) {  
                                 x[i] = u; MULP(vl,u,t,&s); u = s;  
                         }  
                         x[n] = u; y = (P *)ALLOCA((n+1)*sizeof(pointer));  
                         if ( DEG(NEXT(dc)) ) {  
                                 dct = NEXT(dc); MKP(VR(p),dct,t);  
                         } else  
                                 t = COEF(NEXT(dc));  
                         for ( i = 0, u = (P)One; i < n; i++ ) {  
                                 y[i] = u; MULP(vl,u,t,&s); u = s;  
                         }  
                         y[n] = u;  
                         pt = (P *)ALLOCA((n+1)*sizeof(pointer)); MKBC(n,pt);  
                         for ( i = 0, u = 0; i <= n; i++ ) {  
                                 MULP(vl,x[i],y[n-i],&t); MULP(vl,t,pt[i],&s);  
                                 ADDP(vl,u,s,&t); u = t;  
                         }  
                         *pr = u;  
                 } else {  
                         error("exponent too big");  
                         *pr = 0;  
                 }  
         }  
 }  
   
 void D_CHSGNP(p,pr)  
 P p,*pr;  
 {  
         register DCP dc,dcr,dcr0;  
         P t;  
   
         if ( !p )  
                 *pr = NULL;  
         else if ( NUM(p) ) {  
 #if defined(_PA_RISC1_1) || defined(__alpha) || defined(mips) || defined(_IBMR2)  
 #ifdef FBASE  
                 chsgnnum((Num)p,(Num *)pr);  
 #else  
                 MQ mq;  
   
                 NEWMQ(mq); CONT(mq)=mod-CONT((MQ)p); *pr = (P)mq;  
 #endif  
 #else  
                 CHSGNNUM(p,t); *pr = t;  
 #endif  
         } else {  
                 for ( dcr0 = 0, dc = DC(p); dc; dc = NEXT(dc) ) {  
                         NEXTDC(dcr0,dcr); CHSGNP(COEF(dc),&COEF(dcr)); DEG(dcr) = DEG(dc);  
                 }  
                 NEXT(dcr) = 0; MKP(VR(p),dcr0,*pr);  
         }  
 }  
   
   
 #ifdef FBASE  
 void ptozp(p,sgn,c,pr)  
 P p;  
 int sgn;  
 Q *c;  
 P *pr;  
 {  
         N nm,dn;  
   
         if ( !p ) {  
                 *c = 0; *pr = 0;  
         } else {  
                 lgp(p,&nm,&dn);  
                 if ( UNIN(dn) )  
                         NTOQ(nm,sgn,*c);  
                 else  
                         NDTOQ(nm,dn,sgn,*c);  
                 divcp(p,*c,pr);  
         }  
 }  
   
 void lgp(p,g,l)  
 P p;  
 N *g,*l;  
 {  
         N g1,g2,l1,l2,l3,l4,tmp;  
         DCP dc;  
   
         if ( NUM(p) ) {  
                 *g = NM((Q)p);  
                 if ( INT((Q)p) )  
                         *l = ONEN;  
                 else  
                         *l = DN((Q)p);  
         } else {  
                 dc = DC(p); lgp(COEF(dc),g,l);  
                 for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {  
                         lgp(COEF(dc),&g1,&l1); gcdn(*g,g1,&g2); *g = g2;  
                         gcdn(*l,l1,&l2); kmuln(*l,l1,&l3); divn(l3,l2,&l4,&tmp); *l = l4;  
                 }  
         }  
 }  
   
 void divcp(f,q,rp)  
 P f;  
 Q q;  
 P *rp;  
 {  
         DCP dc,dcr,dcr0;  
   
         if ( !f )  
                 *rp = 0;  
         else if ( NUM(f) )  
                 DIVNUM(f,q,rp);  
         else {  
                 for ( dc = DC(f), dcr0 = 0; dc; dc = NEXT(dc) ) {  
                         NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc);  
                         divcp(COEF(dc),q,&COEF(dcr));  
                 }  
                 NEXT(dcr) = 0; MKP(VR(f),dcr0,*rp);  
         }  
 }  
   
 void diffp(vl,p,v,r)  
 VL vl;  
 P p;  
 V v;  
 P *r;  
 {  
         P t;  
         DCP dc,dcr,dcr0;  
   
         if ( !p || NUM(p) )  
                 *r = 0;  
         else {  
                 if ( v == VR(p) ) {  
                         for ( dc = DC(p), dcr0 = 0;  
                                 dc && DEG(dc); dc = NEXT(dc) ) {  
                                 NEXTDC(dcr0,dcr); SUBQ(DEG(dc),ONE,&DEG(dcr));  
                                 MULPQ(COEF(dc),(P)DEG(dc),&COEF(dcr));  
                         }  
                         NEXT(dcr) = 0; MKP(v,dcr0,*r);  
                 } else {  
                         for ( dc = DC(p), dcr0 = 0; dc; dc = NEXT(dc) ) {  
                                 diffp(vl,COEF(dc),v,&t);  
                                 if ( t ) {  
                                         NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = t;  
                                 }  
                         }  
                         if ( !dcr0 )  
                                 *r = 0;  
                         else {  
                                 NEXT(dcr) = 0; MKP(VR(p),dcr0,*r);  
                         }  
                 }  
         }  
 }  
   
 void coefp(p,d,pr)  
 P p;  
 int d;  
 P *pr;  
 {  
         DCP dc;  
         int sgn;  
         Q dq;  
   
         if ( NUM(p) )  
                 if ( d == 0 )  
                         *pr = p;  
                 else  
                         *pr = 0;  
         else {  
                 for ( STOQ(d,dq), dc = DC(p); dc; dc = NEXT(dc) )  
                         if ( (sgn = cmpq(DEG(dc),dq)) > 0 )  
                                 continue;  
                         else if ( sgn == 0 ) {  
                                 *pr = COEF(dc);  
                                 return;  
                         } else {  
                                 *pr = 0;  
                                 break;  
                         }  
                 *pr = 0;  
         }  
 }  
   
 int compp(vl,p1,p2)  
 VL vl;  
 P p1,p2;  
 {  
         DCP dc1,dc2;  
         V v1,v2;  
   
         if ( !p1 )  
                 return p2 ? -1 : 0;  
         else if ( !p2 )  
                 return 1;  
         else if ( NUM(p1) )  
                 return NUM(p2) ? (*cmpnumt[MAX(NID(p1),NID(p2))])(p1,p2) : -1;  
         else if ( NUM(p2) )  
                 return 1;  
         if ( (v1 = VR(p1)) == (v2 = VR(p2)) ) {  
                 for ( dc1 = DC(p1), dc2 = DC(p2);  
                         dc1 && dc2; dc1 = NEXT(dc1), dc2 = NEXT(dc2) )  
                         switch ( cmpq(DEG(dc1),DEG(dc2)) ) {  
                                 case 1:  
                                         return 1; break;  
                                 case -1:  
                                         return -1; break;  
                                 default:  
                                         switch ( compp(vl,COEF(dc1),COEF(dc2)) ) {  
                                                 case 1:  
                                                         return 1; break;  
                                                 case -1:  
                                                         return -1; break;  
                                                 default:  
                                                         break;  
                                         }  
                                         break;  
                         }  
                 return dc1 ? 1 : (dc2 ? -1 : 0);  
         } else {  
                 for ( ; v1 != VR(vl) && v2 != VR(vl); vl = NEXT(vl) );  
                 return v1 == VR(vl) ? 1 : -1;  
         }  
 }  
   
 void csump(vl,p,c)  
 VL vl;  
 P p;  
 Q *c;  
 {  
         DCP dc;  
         Q s,s1,s2;  
   
         if ( !p || NUM(p) )  
                 *c = (Q)p;  
         else {  
                 for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) {  
                         csump(vl,COEF(dc),&s1); addq(s,s1,&s2); s = s2;  
                 }  
                 *c = s;  
         }  
 }  
   
 void degp(v,p,d)  
 V v;  
 P p;  
 Q *d;  
 {  
         DCP dc;  
         Q m,m1;  
   
         if ( !p || NUM(p) )  
                 *d = 0;  
         else if ( v == VR(p) )  
                 *d = DEG(DC(p));  
         else {  
                 for ( dc = DC(p), m = 0; dc; dc = NEXT(dc) ) {  
                         degp(v,COEF(dc),&m1);  
                         if ( cmpq(m,m1) < 0 )  
                                 m = m1;  
                 }  
                 *d = m;  
         }  
 }  
 #endif  
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.6
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>