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

Diff for /OpenXM_contrib2/asir2000/engine/PUM.c between version 1.4 and 1.5

version 1.4, 2001/02/21 07:10:18 version 1.5, 2018/03/29 01:32:51
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/PUM.c,v 1.3 2000/08/22 05:04:04 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/engine/PUM.c,v 1.4 2001/02/21 07:10:18 noro Exp $
 */  */
 #include "ca.h"  #include "ca.h"
   
Line 58  P **rmat;
Line 58  P **rmat;
 int n;  int n;
 P *dp;  P *dp;
 {  {
         int i,j,k,sgn;    int i,j,k,sgn;
         P mjj,mij,t,s,u,d;    P mjj,mij,t,s,u,d;
         P **mat;    P **mat;
         P *mi,*mj;    P *mi,*mj;
   
         mat = (P **)almat_pointer(n,n);    mat = (P **)almat_pointer(n,n);
         for ( i = 0; i < n; i++ )    for ( i = 0; i < n; i++ )
                 for ( j = 0; j < n; j++ )      for ( j = 0; j < n; j++ )
                         mat[i][j] = rmat[i][j];        mat[i][j] = rmat[i][j];
         for ( j = 0, d = (P)ONEM, sgn = 1; j < n; j++ ) {    for ( j = 0, d = (P)ONEM, sgn = 1; j < n; j++ ) {
                 for ( i = j; (i < n) && !mat[i][j]; i++ );      for ( i = j; (i < n) && !mat[i][j]; i++ );
                 if ( i == n ) {      if ( i == n ) {
                         *dp = 0; return;        *dp = 0; return;
                 }      }
                 for ( k = i; k < n; k++ )      for ( k = i; k < n; k++ )
                         if ( mat[k][j] && (nmonop(mat[k][j]) < nmonop(mat[i][j]) ) )        if ( mat[k][j] && (nmonop(mat[k][j]) < nmonop(mat[i][j]) ) )
                                 i = k;          i = k;
                 if ( j != i ) {      if ( j != i ) {
                         mj = mat[j]; mat[j] = mat[i]; mat[i] = mj; sgn = -sgn;        mj = mat[j]; mat[j] = mat[i]; mat[i] = mj; sgn = -sgn;
                 }      }
                 for ( i = j + 1, mj = mat[j], mjj = mj[j]; i < n; i++ )      for ( i = j + 1, mj = mat[j], mjj = mj[j]; i < n; i++ )
                         for ( k = j + 1, mi = mat[i], mij = mi[j]; k < n; k++ ) {        for ( k = j + 1, mi = mat[i], mij = mi[j]; k < n; k++ ) {
                                 mulmp(vl,mod,mi[k],mjj,&t); mulmp(vl,mod,mj[k],mij,&s);          mulmp(vl,mod,mi[k],mjj,&t); mulmp(vl,mod,mj[k],mij,&s);
                                 submp(vl,mod,t,s,&u); divsmp(vl,mod,u,d,&mi[k]);          submp(vl,mod,t,s,&u); divsmp(vl,mod,u,d,&mi[k]);
                         }        }
                 d = mjj;      d = mjj;
         }    }
         if ( sgn < 0 )    if ( sgn < 0 )
                 chsgnmp(mod,d,dp);      chsgnmp(mod,d,dp);
         else    else
                 *dp = d;      *dp = d;
 }  }
   
 void resultmp(vl,mod,v,p1,p2,pr)  void resultmp(vl,mod,v,p1,p2,pr)
Line 97  int mod;
Line 97  int mod;
 V v;  V v;
 P p1,p2,*pr;  P p1,p2,*pr;
 {  {
         P q1,q2,m,m1,m2,lc,q,r,t,g1,g2,adj;    P q1,q2,m,m1,m2,lc,q,r,t,g1,g2,adj;
         int d,d1,d2,j,k;    int d,d1,d2,j,k;
         VL nvl;    VL nvl;
         Q dq;    Q dq;
         MQ mq;    MQ mq;
   
         if ( !p1 || !p2 ) {    if ( !p1 || !p2 ) {
                 *pr = 0; return;      *pr = 0; return;
         }    }
         reordvar(vl,v,&nvl);    reordvar(vl,v,&nvl);
         reordermp(nvl,mod,vl,p1,&q1); reordermp(nvl,mod,vl,p2,&q2);    reordermp(nvl,mod,vl,p1,&q1); reordermp(nvl,mod,vl,p2,&q2);
   
         if ( VR(q1) != v )    if ( VR(q1) != v )
                 if ( VR(q2) != v ) {      if ( VR(q2) != v ) {
                         *pr = 0;        *pr = 0;
                         return;        return;
                 } else {      } else {
                         d = deg(v,q2); STOQ(d,dq);        d = deg(v,q2); STOQ(d,dq);
                         pwrmp(vl,mod,q1,dq,pr);        pwrmp(vl,mod,q1,dq,pr);
                         return;        return;
                 }      }
         else if ( VR(q2) != v ) {    else if ( VR(q2) != v ) {
                 d = deg(v,q1); STOQ(d,dq);      d = deg(v,q1); STOQ(d,dq);
                 pwrmp(vl,mod,q2,dq,pr);      pwrmp(vl,mod,q2,dq,pr);
                 return;      return;
         }    }
   
         if ( ( VR(q1) != v ) || ( VR(q2) != v ) ) {    if ( ( VR(q1) != v ) || ( VR(q2) != v ) ) {
                 *pr = 0;      *pr = 0;
                 return;      return;
         }    }
   
         d1 = deg(v,q1); d2 = deg(v,q2);    d1 = deg(v,q1); d2 = deg(v,q2);
         if ( d1 > d2 ) {    if ( d1 > d2 ) {
                 g1 = q1; g2 = q2; adj = (P)ONEM;      g1 = q1; g2 = q2; adj = (P)ONEM;
         } else if ( d1 < d2 ) {    } else if ( d1 < d2 ) {
                 g2 = q1; g1 = q2;      g2 = q1; g1 = q2;
                 if ( (d1 % 2) && (d2 % 2) ) {      if ( (d1 % 2) && (d2 % 2) ) {
                         STOMQ(-1,mq); adj = (P)mq;        STOMQ(-1,mq); adj = (P)mq;
                 } else      } else
                         adj = (P)ONEM;        adj = (P)ONEM;
         } else {    } else {
                 premmp(nvl,mod,q1,q2,&t);      premmp(nvl,mod,q1,q2,&t);
                 d = deg(v,t); STOQ(d,dq); pwrmp(nvl,mod,LC(q2),dq,&adj);      d = deg(v,t); STOQ(d,dq); pwrmp(nvl,mod,LC(q2),dq,&adj);
                 g1 = q2; g2 = t;      g1 = q2; g2 = t;
                 if ( d1 % 2 ) {      if ( d1 % 2 ) {
                         chsgnmp(mod,adj,&t); adj = t;        chsgnmp(mod,adj,&t); adj = t;
                 }      }
         }    }
         d1 = deg(v,g1); j = d1 - 1;    d1 = deg(v,g1); j = d1 - 1;
   
         for ( lc = (P)ONEM; ; ) {    for ( lc = (P)ONEM; ; ) {
                 if ( ( k = deg(v,g2) ) < 0 ) {      if ( ( k = deg(v,g2) ) < 0 ) {
                         *pr = 0;        *pr = 0;
                         return;        return;
                 }      }
   
                 if ( k == j )      if ( k == j )
                         if ( !k ) {        if ( !k ) {
                                 divsmp(nvl,mod,g2,adj,pr);          divsmp(nvl,mod,g2,adj,pr);
                                 return;          return;
                         } else {        } else {
                                 premmp(nvl,mod,g1,g2,&r); mulmp(nvl,mod,lc,lc,&m);          premmp(nvl,mod,g1,g2,&r); mulmp(nvl,mod,lc,lc,&m);
                                 divsmp(nvl,mod,r,m,&q);          divsmp(nvl,mod,r,m,&q);
                                 g1 = g2; g2 = q;          g1 = g2; g2 = q;
                                 lc = LC(g1); /* g1 is not const */          lc = LC(g1); /* g1 is not const */
                                 j = k - 1;          j = k - 1;
                         }        }
                 else {      else {
                         d = j - k; STOQ(d,dq);        d = j - k; STOQ(d,dq);
                         pwrmp(nvl,mod,(VR(g2)==v?LC(g2):g2),dq,&m);        pwrmp(nvl,mod,(VR(g2)==v?LC(g2):g2),dq,&m);
                         mulmp(nvl,mod,g2,m,&m1);        mulmp(nvl,mod,g2,m,&m1);
                         pwrmp(nvl,mod,lc,dq,&m); divsmp(nvl,mod,m1,m,&t);        pwrmp(nvl,mod,lc,dq,&m); divsmp(nvl,mod,m1,m,&t);
                         if ( k == 0 ) {        if ( k == 0 ) {
                                 divsmp(nvl,mod,t,adj,pr);          divsmp(nvl,mod,t,adj,pr);
                                 return;          return;
                         } else {        } else {
                                 premmp(nvl,mod,g1,g2,&r);          premmp(nvl,mod,g1,g2,&r);
                                 mulmp(nvl,mod,lc,lc,&m1); mulmp(nvl,mod,m,m1,&m2);          mulmp(nvl,mod,lc,lc,&m1); mulmp(nvl,mod,m,m1,&m2);
                                 divsmp(nvl,mod,r,m2,&q);          divsmp(nvl,mod,r,m2,&q);
                                 g1 = t; g2 = q;          g1 = t; g2 = q;
                                 if ( d % 2 ) {          if ( d % 2 ) {
                                         chsgnmp(mod,g2,&t); g2 = t;            chsgnmp(mod,g2,&t); g2 = t;
                                 }          }
                                 lc = LC(g1); /* g1 is not const */          lc = LC(g1); /* g1 is not const */
                                 j = k - 1;          j = k - 1;
                         }        }
                 }      }
         }    }
 }  }
   
 void premmp(vl,mod,p1,p2,pr)  void premmp(vl,mod,p1,p2,pr)
Line 193  VL vl;
Line 193  VL vl;
 int mod;  int mod;
 P p1,p2,*pr;  P p1,p2,*pr;
 {  {
         P m,m1,m2;    P m,m1,m2;
         P *pw;    P *pw;
         DCP dc;    DCP dc;
         V v1,v2;    V v1,v2;
         register int i,j;    register int i,j;
         int n1,n2,d;    int n1,n2,d;
   
         if ( NUM(p1) )    if ( NUM(p1) )
                 if ( NUM(p2) )      if ( NUM(p2) )
                         *pr = 0;        *pr = 0;
                 else      else
                         *pr = p1;        *pr = p1;
         else if ( NUM(p2) )    else if ( NUM(p2) )
                 *pr = 0;      *pr = 0;
         else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {    else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
                 n1 = deg(v1,p1); n2 = deg(v1,p2);      n1 = deg(v1,p1); n2 = deg(v1,p2);
                 pw = (P *)ALLOCA((n1+1)*sizeof(P));      pw = (P *)ALLOCA((n1+1)*sizeof(P));
                 bzero((char *)pw,(int)((n1+1)*sizeof(P)));      bzero((char *)pw,(int)((n1+1)*sizeof(P)));
   
                 for ( dc = DC(p1); dc; dc = NEXT(dc) )      for ( dc = DC(p1); dc; dc = NEXT(dc) )
                         pw[QTOS(DEG(dc))] = COEF(dc);        pw[QTOS(DEG(dc))] = COEF(dc);
   
                 for ( i = n1; i >= n2; i-- ) {      for ( i = n1; i >= n2; i-- ) {
                         if ( pw[i] ) {        if ( pw[i] ) {
                                 chsgnmp(mod,pw[i],&m);          chsgnmp(mod,pw[i],&m);
                                 for ( j = i; j >= 0; j-- ) {          for ( j = i; j >= 0; j-- ) {
                                         mulmp(vl,mod,pw[j],LC(p2),&m1); pw[j] = m1;            mulmp(vl,mod,pw[j],LC(p2),&m1); pw[j] = m1;
                                 }          }
   
                                 for ( dc = DC(p2), d = i - n2; dc; dc = NEXT(dc) ) {          for ( dc = DC(p2), d = i - n2; dc; dc = NEXT(dc) ) {
                                         mulmp(vl,mod,COEF(dc),m,&m1);            mulmp(vl,mod,COEF(dc),m,&m1);
                                         addmp(vl,mod,pw[QTOS(DEG(dc))+d],m1,&m2);            addmp(vl,mod,pw[QTOS(DEG(dc))+d],m1,&m2);
                                         pw[QTOS(DEG(dc))+d] = m2;            pw[QTOS(DEG(dc))+d] = m2;
                                 }          }
                         } else        } else
                                 for ( j = i; j >= 0; j-- )          for ( j = i; j >= 0; j-- )
                                         if ( pw[j] ) {            if ( pw[j] ) {
                                                 mulmp(vl,mod,pw[j],LC(p2),&m1); pw[j] = m1;              mulmp(vl,mod,pw[j],LC(p2),&m1); pw[j] = m1;
                                         }            }
                 }      }
                 plisttop(pw,v1,n2-1,pr);      plisttop(pw,v1,n2-1,pr);
         } else {    } else {
                 while ( v1 != vl->v && v2 != vl->v )      while ( v1 != vl->v && v2 != vl->v )
                         vl = NEXT(vl);        vl = NEXT(vl);
                 if ( v1 == vl->v )      if ( v1 == vl->v )
                         *pr = 0;        *pr = 0;
                 else      else
                         *pr = p1;        *pr = p1;
         }    }
 }  }
   
 void srchmp(vl,mod,v,p1,p2,pr)  void srchmp(vl,mod,v,p1,p2,pr)
Line 250  int mod;
Line 250  int mod;
 V v;  V v;
 P p1,p2,*pr;  P p1,p2,*pr;
 {  {
         P a,b,q1,q2,x,t,s,d,bg,c,c0,db;    P a,b,q1,q2,x,t,s,d,bg,c,c0,db;
         int i,m,k;    int i,m,k;
         V v0;    V v0;
         VL nvl,tvl,nvl1,nvl2;    VL nvl,tvl,nvl1,nvl2;
         Q dq;    Q dq;
         MQ q;    MQ q;
   
         if ( vl->v != v ) {    if ( vl->v != v ) {
                 reordvar(vl,v,&tvl);      reordvar(vl,v,&tvl);
                 reordermp(tvl,mod,vl,p1,&q1); reordermp(tvl,mod,vl,p2,&q2);      reordermp(tvl,mod,vl,p1,&q1); reordermp(tvl,mod,vl,p2,&q2);
         } else {    } else {
                 q1 = p1; q2 = p2; tvl = vl;      q1 = p1; q2 = p2; tvl = vl;
         }    }
         clctv(tvl,q1,&nvl1); clctv(tvl,q2,&nvl2); mergev(tvl,nvl1,nvl2,&nvl);    clctv(tvl,q1,&nvl1); clctv(tvl,q2,&nvl2); mergev(tvl,nvl1,nvl2,&nvl);
         if ( VR(q1) != v )    if ( VR(q1) != v )
                 if ( VR(q2) != v )      if ( VR(q2) != v )
                         *pr = 0;        *pr = 0;
                 else {      else {
                         m = getdeg(v,q2); STOQ(m,dq); pwrmp(vl,mod,q1,dq,pr);        m = getdeg(v,q2); STOQ(m,dq); pwrmp(vl,mod,q1,dq,pr);
                 }      }
         else if ( VR(q2) != v ) {    else if ( VR(q2) != v ) {
                 m = getdeg(v,q1); STOQ(m,dq); pwrmp(vl,mod,q2,dq,pr);      m = getdeg(v,q1); STOQ(m,dq); pwrmp(vl,mod,q2,dq,pr);
         } else if ( !NEXT(nvl) )    } else if ( !NEXT(nvl) )
                 srchump(mod,p1,p2,pr);      srchump(mod,p1,p2,pr);
         else {    else {
                 v0 = NEXT(nvl)->v;      v0 = NEXT(nvl)->v;
                 k = getdeg(v,q1)*getdeg(v0,q2)+getdeg(v,q2)*getdeg(v0,q1)+1;      k = getdeg(v,q1)*getdeg(v0,q2)+getdeg(v,q2)*getdeg(v0,q1)+1;
                 for ( i = 0, c = 0, d = (P)ONEM, MKMV(v0,x);      for ( i = 0, c = 0, d = (P)ONEM, MKMV(v0,x);
                         ( i < mod ) && (getdeg(v0,d) < k) ; i++ ) {        ( i < mod ) && (getdeg(v0,d) < k) ; i++ ) {
                         STOMQ(i,q),bg = (P)q; substmp(nvl,mod,LC(q1),v0,bg,&t);        STOMQ(i,q),bg = (P)q; substmp(nvl,mod,LC(q1),v0,bg,&t);
                         if ( !t )        if ( !t )
                                 continue;          continue;
                         substmp(nvl,mod,LC(q2),v0,bg,&t);        substmp(nvl,mod,LC(q2),v0,bg,&t);
                         if ( !t )        if ( !t )
                                 continue;          continue;
                         substmp(nvl,mod,q1,v0,bg,&a); substmp(nvl,mod,q2,v0,bg,&b);        substmp(nvl,mod,q1,v0,bg,&a); substmp(nvl,mod,q2,v0,bg,&b);
                         srchmp(nvl,mod,v,a,b,&c0); substmp(nvl,mod,c,v0,bg,&t);        srchmp(nvl,mod,v,a,b,&c0); substmp(nvl,mod,c,v0,bg,&t);
                         submp(nvl,mod,c0,t,&s); mulmp(nvl,mod,s,d,&t);        submp(nvl,mod,c0,t,&s); mulmp(nvl,mod,s,d,&t);
                         substmp(nvl,mod,d,v0,bg,&db);        substmp(nvl,mod,d,v0,bg,&db);
                         divsmp(nvl,mod,t,db,&s); addmp(nvl,mod,s,c,&t); c = t;        divsmp(nvl,mod,t,db,&s); addmp(nvl,mod,s,c,&t); c = t;
                         submp(nvl,mod,x,bg,&t); mulmp(nvl,mod,d,t,&s); d = s;        submp(nvl,mod,x,bg,&t); mulmp(nvl,mod,d,t,&s); d = s;
                 }      }
                 if ( i == mod )      if ( i == mod )
                         error("srchmp : ???");        error("srchmp : ???");
                 *pr = c;      *pr = c;
         }    }
 }  }
   
 int ucmpp(p,q)  int ucmpp(p,q)
 P p,q;  P p,q;
 {  {
         DCP dcp,dcq;    DCP dcp,dcq;
   
         if ( !p )    if ( !p )
                 if ( !q )      if ( !q )
                         return ( 0 );        return ( 0 );
                 else      else
                         return ( 1 );        return ( 1 );
         else if ( !q )    else if ( !q )
                 return ( 1 );      return ( 1 );
         else if ( NUM(p) )    else if ( NUM(p) )
                 if ( !NUM(q) )      if ( !NUM(q) )
                         return ( 1 );        return ( 1 );
                 else      else
                         return ( cmpq((Q)p,(Q)q) );        return ( cmpq((Q)p,(Q)q) );
         else if ( NUM(q) )    else if ( NUM(q) )
                         return ( 1 );        return ( 1 );
         else {    else {
                 for ( dcp = DC(p), dcq = DC(q); dcp && dcq;      for ( dcp = DC(p), dcq = DC(q); dcp && dcq;
                         dcp = NEXT(dcp), dcq = NEXT(dcq) )        dcp = NEXT(dcp), dcq = NEXT(dcq) )
                         if ( cmpq(DEG(dcp),DEG(dcq) ) )        if ( cmpq(DEG(dcp),DEG(dcq) ) )
                                 return ( 1 );          return ( 1 );
                         else if ( cmpq((Q)COEF(dcp),(Q)COEF(dcq) ) )        else if ( cmpq((Q)COEF(dcp),(Q)COEF(dcq) ) )
                                 return ( 1 );          return ( 1 );
                 if ( dcp || dcq )      if ( dcp || dcq )
                         return ( 1 );        return ( 1 );
                 else      else
                         return ( 0 );        return ( 0 );
         }    }
 }  }
   
 #if 0  #if 0
Line 336  srchump(mod,p1,p2,pr)
Line 336  srchump(mod,p1,p2,pr)
 int mod;  int mod;
 P p1,p2,*pr;  P p1,p2,*pr;
 {  {
         int r;    int r;
         MQ mq;    MQ mq;
   
         r = eucresum(mod,p1,p2);    r = eucresum(mod,p1,p2);
         STOMQ(r,mq); *pr = (P)mq;    STOMQ(r,mq); *pr = (P)mq;
 }  }
 #endif  #endif
   
Line 348  void srchump(mod,p1,p2,pr)
Line 348  void srchump(mod,p1,p2,pr)
 int mod;  int mod;
 P p1,p2,*pr;  P p1,p2,*pr;
 {  {
         UM m,m1,q,r,t,g1,g2;    UM m,m1,q,r,t,g1,g2;
         int lc,d,d1,d2,i,j,k,l,l1,l2,tmp,adj;    int lc,d,d1,d2,i,j,k,l,l1,l2,tmp,adj;
         V v;    V v;
   
         v = VR(p1); d = MAX(UDEG(p1),UDEG(p2));    v = VR(p1); d = MAX(UDEG(p1),UDEG(p2));
         g1 = W_UMALLOC(d); g2 = W_UMALLOC(d);    g1 = W_UMALLOC(d); g2 = W_UMALLOC(d);
         bzero((char *)g1,(int)((d+2)*sizeof(int))); bzero((char *)g2,(int)((d+2)*sizeof(int)));    bzero((char *)g1,(int)((d+2)*sizeof(int))); bzero((char *)g2,(int)((d+2)*sizeof(int)));
         if ( d == (int)UDEG(p1) ) {    if ( d == (int)UDEG(p1) ) {
                 mptoum(p1,g1); mptoum(p2,g2);      mptoum(p1,g1); mptoum(p2,g2);
         } else {    } else {
                 mptoum(p2,g1); mptoum(p1,g2);      mptoum(p2,g1); mptoum(p1,g2);
         }    }
         if ( ( d1 = DEG(g1) ) > ( d2 = DEG(g2) ) ) {    if ( ( d1 = DEG(g1) ) > ( d2 = DEG(g2) ) ) {
                 j = d1 - 1; adj = 1;      j = d1 - 1; adj = 1;
         } else    } else
                 j = d2;      j = d2;
         lc = 1;    lc = 1;
         r = W_UMALLOC(d1+d2); q = W_UMALLOC(d1+d2);    r = W_UMALLOC(d1+d2); q = W_UMALLOC(d1+d2);
         m1 = W_UMALLOC(d1+d2); t = W_UMALLOC(d1+d2);    m1 = W_UMALLOC(d1+d2); t = W_UMALLOC(d1+d2);
         bzero((char *)r,(int)((d1+d2+2)*sizeof(int))); bzero((char *)q,(int)((d1+d2+2)*sizeof(int)));    bzero((char *)r,(int)((d1+d2+2)*sizeof(int))); bzero((char *)q,(int)((d1+d2+2)*sizeof(int)));
         bzero((char *)m1,(int)((d1+d2+2)*sizeof(int))); bzero((char *)t,(int)((d1+d2+2)*sizeof(int)));    bzero((char *)m1,(int)((d1+d2+2)*sizeof(int))); bzero((char *)t,(int)((d1+d2+2)*sizeof(int)));
         m = W_UMALLOC(0); bzero((char *)m,(int)(2*sizeof(int)));    m = W_UMALLOC(0); bzero((char *)m,(int)(2*sizeof(int)));
         adj = pwrm(mod,COEF(g2)[DEG(g2)],DEG(g1));    adj = pwrm(mod,COEF(g2)[DEG(g2)],DEG(g1));
         DEG(m) = 0; COEF(m)[0] = invm(COEF(g2)[DEG(g2)],mod);    DEG(m) = 0; COEF(m)[0] = invm(COEF(g2)[DEG(g2)],mod);
         mulum(mod,g2,m,r); cpyum(r,g2);    mulum(mod,g2,m,r); cpyum(r,g2);
         while ( 1 ) {    while ( 1 ) {
                 if ( ( k = DEG(g2) ) < 0 ) {      if ( ( k = DEG(g2) ) < 0 ) {
                         *pr = 0;        *pr = 0;
                         return;        return;
                 }      }
                 if ( k == j ) {      if ( k == j ) {
                         if ( k == 0 ) {        if ( k == 0 ) {
                                 DEG(m) = 0; COEF(m)[0] = adj;          DEG(m) = 0; COEF(m)[0] = adj;
                                 mulum(mod,g2,m,r); umtomp(v,r,pr);          mulum(mod,g2,m,r); umtomp(v,r,pr);
                                 return;          return;
                         } else {        } else {
                                 DEG(m) = 0;          DEG(m) = 0;
                                 COEF(m)[0] = pwrm(mod,COEF(g2)[k],DEG(g1)-k+1);          COEF(m)[0] = pwrm(mod,COEF(g2)[k],DEG(g1)-k+1);
                                 mulum(mod,g1,m,r); DEG(r) = divum(mod,r,g2,t);          mulum(mod,g1,m,r); DEG(r) = divum(mod,r,g2,t);
                                 DEG(m) = 0; COEF(m)[0] = dmb(mod,lc,lc,&tmp);          DEG(m) = 0; COEF(m)[0] = dmb(mod,lc,lc,&tmp);
                                 divum(mod,r,m,q); cpyum(g2,g1); cpyum(q,g2);          divum(mod,r,m,q); cpyum(g2,g1); cpyum(q,g2);
                                 lc = COEF(g1)[DEG(g1)]; j = k - 1;          lc = COEF(g1)[DEG(g1)]; j = k - 1;
                         }        }
                 } else {      } else {
                         d = j - k;        d = j - k;
                         DEG(m) = 0; COEF(m)[0] = pwrm(mod,COEF(g2)[DEG(g2)],d);        DEG(m) = 0; COEF(m)[0] = pwrm(mod,COEF(g2)[DEG(g2)],d);
                         mulum(mod,g2,m,m1); l = pwrm(mod,lc,d);        mulum(mod,g2,m,m1); l = pwrm(mod,lc,d);
                         DEG(m) = 0; COEF(m)[0] = l; divum(mod,m1,m,t);        DEG(m) = 0; COEF(m)[0] = l; divum(mod,m1,m,t);
                         if ( k == 0 ) {        if ( k == 0 ) {
                                 DEG(m) = 0; COEF(m)[0] = adj;          DEG(m) = 0; COEF(m)[0] = adj;
                                 mulum(mod,t,m,r); umtomp(v,r,pr);          mulum(mod,t,m,r); umtomp(v,r,pr);
                                 return;          return;
                         } else {        } else {
                                 DEG(m) = 0;          DEG(m) = 0;
                                 COEF(m)[0] = pwrm(mod,COEF(g2)[k],DEG(g1)-k+1);          COEF(m)[0] = pwrm(mod,COEF(g2)[k],DEG(g1)-k+1);
                                 mulum(mod,g1,m,r); DEG(r) = divum(mod,r,g2,q);          mulum(mod,g1,m,r); DEG(r) = divum(mod,r,g2,q);
                                 l1 = dmb(mod,lc,lc,&tmp); l2 = dmb(mod,l,l1,&tmp);          l1 = dmb(mod,lc,lc,&tmp); l2 = dmb(mod,l,l1,&tmp);
                                 DEG(m) = 0; COEF(m)[0] = l2;          DEG(m) = 0; COEF(m)[0] = l2;
                                 divum(mod,r,m,q); cpyum(t,g1); cpyum(q,g2);          divum(mod,r,m,q); cpyum(t,g1); cpyum(q,g2);
                                 if ( d % 2 )          if ( d % 2 )
                                         for ( i = DEG(g2); i >= 0; i-- )            for ( i = DEG(g2); i >= 0; i-- )
                                                 COEF(g2)[i] = ( mod - COEF(g2)[i] ) % mod;              COEF(g2)[i] = ( mod - COEF(g2)[i] ) % mod;
                                 lc = COEF(g1)[DEG(g1)]; j = k - 1;          lc = COEF(g1)[DEG(g1)]; j = k - 1;
                         }        }
                 }      }
         }    }
 }  }
   
 void substmp(vl,mod,p,v0,p0,pr)  void substmp(vl,mod,p,v0,p0,pr)
Line 423  V v0;
Line 423  V v0;
 P p,p0;  P p,p0;
 P *pr;  P *pr;
 {  {
         P x,t,m,c,s,a;    P x,t,m,c,s,a;
         DCP dc;    DCP dc;
         Q d;    Q d;
   
         if ( !p )    if ( !p )
                 *pr = 0;      *pr = 0;
         else if ( NUM(p) )    else if ( NUM(p) )
                 *pr = p;      *pr = p;
         else if ( VR(p) != v0 ) {    else if ( VR(p) != v0 ) {
                 MKMV(VR(p),x);      MKMV(VR(p),x);
                 for ( c = 0, dc = DC(p); dc; dc = NEXT(dc) ) {      for ( c = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
                         substmp(vl,mod,COEF(dc),v0,p0,&t);        substmp(vl,mod,COEF(dc),v0,p0,&t);
                         if ( DEG(dc) ) {        if ( DEG(dc) ) {
                                 pwrmp(vl,mod,x,DEG(dc),&s); mulmp(vl,mod,s,t,&m);          pwrmp(vl,mod,x,DEG(dc),&s); mulmp(vl,mod,s,t,&m);
                                 addmp(vl,mod,m,c,&a);          addmp(vl,mod,m,c,&a);
                                 c = a;          c = a;
                         } else {        } else {
                                 addmp(vl,mod,t,c,&a);          addmp(vl,mod,t,c,&a);
                                 c = a;          c = a;
                         }        }
                 }      }
                 *pr = c;      *pr = c;
         } else {    } else {
                 dc = DC(p);      dc = DC(p);
                 c = COEF(dc);      c = COEF(dc);
                 for ( d = DEG(dc), dc = NEXT(dc);      for ( d = DEG(dc), dc = NEXT(dc);
                         dc; d = DEG(dc), dc = NEXT(dc) ) {        dc; d = DEG(dc), dc = NEXT(dc) ) {
                                 subq(d,DEG(dc),(Q *)&t); pwrmp(vl,mod,p0,(Q)t,&s);          subq(d,DEG(dc),(Q *)&t); pwrmp(vl,mod,p0,(Q)t,&s);
                                 mulmp(vl,mod,s,c,&m);          mulmp(vl,mod,s,c,&m);
                                 addmp(vl,mod,m,COEF(dc),&c);          addmp(vl,mod,m,COEF(dc),&c);
                 }      }
                 if ( d ) {      if ( d ) {
                         pwrmp(vl,mod,p0,d,&t); mulmp(vl,mod,t,c,&m);        pwrmp(vl,mod,p0,d,&t); mulmp(vl,mod,t,c,&m);
                         c = m;        c = m;
                 }      }
                 *pr = c;      *pr = c;
         }    }
 }  }
   
 void reordermp(nvl,mod,ovl,p,pr)  void reordermp(nvl,mod,ovl,p,pr)
Line 468  int mod;
Line 468  int mod;
 P p;  P p;
 P *pr;  P *pr;
 {  {
         DCP dc;    DCP dc;
         P x,m,s,t,c;    P x,m,s,t,c;
   
         if ( !p )    if ( !p )
                 *pr = 0;      *pr = 0;
         else if ( NUM(p) )    else if ( NUM(p) )
                 *pr = p;      *pr = p;
         else {    else {
                 MKMV(VR(p),x);      MKMV(VR(p),x);
                 for ( s = 0, dc = DC(p); dc; dc = NEXT(dc) ) {      for ( s = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
                         reordermp(nvl,mod,ovl,COEF(dc),&c);        reordermp(nvl,mod,ovl,COEF(dc),&c);
                         if ( DEG(dc) ) {        if ( DEG(dc) ) {
                                 pwrmp(nvl,mod,x,DEG(dc),&t); mulmp(nvl,mod,c,t,&m);          pwrmp(nvl,mod,x,DEG(dc),&t); mulmp(nvl,mod,c,t,&m);
                                 addmp(nvl,mod,m,s,&t);          addmp(nvl,mod,m,s,&t);
                         } else        } else
                                 addmp(nvl,mod,s,c,&t);          addmp(nvl,mod,s,c,&t);
                         s = t;        s = t;
                 }      }
                 *pr = s;      *pr = s;
         }    }
 }  }
   
 void chnremp(vl,mod,p,q,c,r)  void chnremp(vl,mod,p,q,c,r)
 VL vl;  VL vl;
 int mod;  int mod;
Line 498  Q q;
Line 498  Q q;
 P c;  P c;
 P *r;  P *r;
 {  {
         P tg,sg,ug;    P tg,sg,ug;
         P t,u;    P t,u;
         MQ mq;    MQ mq;
   
         ptomp(mod,p,&tg); submp(vl,mod,c,tg,&sg);    ptomp(mod,p,&tg); submp(vl,mod,c,tg,&sg);
         UTOMQ(rem(NM(q),mod),mq),tg = (P)mq; divsmp(vl,mod,sg,tg,&ug);    UTOMQ(rem(NM(q),mod),mq),tg = (P)mq; divsmp(vl,mod,sg,tg,&ug);
         normalizemp(mod,ug);    normalizemp(mod,ug);
         mptop(ug,&u); mulp(vl,u,(P)q,&t); addp(vl,t,p,r);    mptop(ug,&u); mulp(vl,u,(P)q,&t); addp(vl,t,p,r);
 }  }
   
 /* XXX  strange behavior of invm() on SPARC */  /* XXX  strange behavior of invm() on SPARC */
Line 517  UM t;
Line 517  UM t;
 Q q,*qr;  Q q,*qr;
 P c,*cr;  P c,*cr;
 {  {
         int n,m,i,d,a,sd,tmp;    int n,m,i,d,a,sd,tmp;
         Q b,s,z;    Q b,s,z;
         Q *pc,*pcr;    Q *pc,*pcr;
         DCP dc;    DCP dc;
   
         if ( !c || NUM(c) )    if ( !c || NUM(c) )
                 n = 0;      n = 0;
         else    else
                 n = UDEG(c);      n = UDEG(c);
         m = DEG(t); d = MAX(n,m); W_CALLOC(n,Q,pc); W_CALLOC(d,Q,pcr);    m = DEG(t); d = MAX(n,m); W_CALLOC(n,Q,pc); W_CALLOC(d,Q,pcr);
         if ( !c )    if ( !c )
                 pc[0] = 0;      pc[0] = 0;
         else if ( NUM(c) )    else if ( NUM(c) )
                 pc[0] = (Q)c;      pc[0] = (Q)c;
         else    else
                 for ( dc = DC(c); dc; dc = NEXT(dc) )      for ( dc = DC(c); dc; dc = NEXT(dc) )
                         pc[QTOS(DEG(dc))] = (Q)COEF(dc);        pc[QTOS(DEG(dc))] = (Q)COEF(dc);
         for ( i = 0; i <= d; i++ ) {    for ( i = 0; i <= d; i++ ) {
                 b = (i>n?0:pc[i]); a = (i>m?0:COEF(t)[i]);      b = (i>n?0:pc[i]); a = (i>m?0:COEF(t)[i]);
                 if ( b )      if ( b )
                         a = (a-SGN(pc[i])*((int)rem(NM(pc[i]),mod)))%mod;        a = (a-SGN(pc[i])*((int)rem(NM(pc[i]),mod)))%mod;
                 sd = dmb(mod,(a>=0?a:a+mod),invm(rem(NM(q),mod),mod),&tmp);      sd = dmb(mod,(a>=0?a:a+mod),invm(rem(NM(q),mod),mod),&tmp);
                 if ( ( 2 * sd ) > mod )      if ( ( 2 * sd ) > mod )
                         sd -= mod;        sd -= mod;
                 STOQ(sd,z); mulq(z,q,&s); addq(s,b,&pcr[i]);      STOQ(sd,z); mulq(z,q,&s); addq(s,b,&pcr[i]);
         }    }
         STOQ(mod,z); mulq(q,z,qr); plisttop((P *)pcr,v,d,cr);    STOQ(mod,z); mulq(q,z,qr); plisttop((P *)pcr,v,d,cr);
 }  }
   
 void normalizemp(mod,g)  void normalizemp(mod,g)
 int mod;  int mod;
 P g;  P g;
 {  {
         DCP dc;    DCP dc;
   
         if ( !g )    if ( !g )
                 return;      return;
         else if ( NUM(g) ) {    else if ( NUM(g) ) {
                 if ( 2 * CONT((MQ)g) > mod )      if ( 2 * CONT((MQ)g) > mod )
                         CONT((MQ)g) -= mod;        CONT((MQ)g) -= mod;
                 return;      return;
         } else    } else
                 for ( dc = DC(g); dc; dc = NEXT(dc) )      for ( dc = DC(g); dc; dc = NEXT(dc) )
                         normalizemp(mod,COEF(dc));        normalizemp(mod,COEF(dc));
 }  }
   
 void norm(p,r)  void norm(p,r)
 P p;  P p;
 Q *r;  Q *r;
 {  {
         N t;    N t;
         DCP dc;    DCP dc;
   
         for ( dc = DC(p), t = ONEN; dc; dc = NEXT(dc) )    for ( dc = DC(p), t = ONEN; dc; dc = NEXT(dc) )
                 if ( cmpn(NM((Q)COEF(dc)),t) > 0 ) {      if ( cmpn(NM((Q)COEF(dc)),t) > 0 ) {
                         t = NM((Q)COEF(dc));        t = NM((Q)COEF(dc));
                 }      }
         NTOQ(t,1,*r);    NTOQ(t,1,*r);
 }  }
   
 void norm1(p,r)  void norm1(p,r)
 P p,*r;  P p,*r;
 {  {
         DCP dc;    DCP dc;
         P t,s,u;    P t,s,u;
         Q q;    Q q;
   
         if ( NUM(p) )    if ( NUM(p) )
                 NTOQ(NM((Q)p),1,q),*r = (P)q;      NTOQ(NM((Q)p),1,q),*r = (P)q;
         else {    else {
                 for ( t = 0, dc = DC(p); dc; dc = NEXT(dc) ) {      for ( t = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
                         norm1(COEF(dc),&s); addq((Q)t,(Q)s,(Q *)&u); t = u;        norm1(COEF(dc),&s); addq((Q)t,(Q)s,(Q *)&u); t = u;
                 }      }
                 *r = t;      *r = t;
         }    }
 }  }
   
 void norm1c(p,r)  void norm1c(p,r)
 P p;  P p;
 Q *r;  Q *r;
 {  {
         N t;    N t;
         Q s;    Q s;
         DCP dc;    DCP dc;
   
         if ( NUM(p) )    if ( NUM(p) )
                 norm1(p,(P *)r);      norm1(p,(P *)r);
         else {    else {
                 for ( dc = DC(p), t = ONEN; dc; dc = NEXT(dc) ) {      for ( dc = DC(p), t = ONEN; dc; dc = NEXT(dc) ) {
                         norm1(COEF(dc),(P *)&s);        norm1(COEF(dc),(P *)&s);
                         if ( cmpn(NM(s),t) > 0 )        if ( cmpn(NM(s),t) > 0 )
                                 t = NM(s);          t = NM(s);
                 }      }
                 NTOQ(t,1,*r);      NTOQ(t,1,*r);
         }    }
 }  }
   
 void gcdprsmp(vl,mod,p1,p2,pr)  void gcdprsmp(vl,mod,p1,p2,pr)
Line 619  VL vl;
Line 619  VL vl;
 int mod;  int mod;
 P p1,p2,*pr;  P p1,p2,*pr;
 {  {
         P g1,g2,gc1,gc2,gp1,gp2,g,gc,gp,gcr;    P g1,g2,gc1,gc2,gp1,gp2,g,gc,gp,gcr;
         V v1,v2;    V v1,v2;
   
         if ( !p1 )    if ( !p1 )
                 *pr = p2;      *pr = p2;
         else if ( !p2 )    else if ( !p2 )
                 *pr = p1;      *pr = p1;
         else if ( NUM(p1) || NUM(p2) )    else if ( NUM(p1) || NUM(p2) )
                 *pr = (P)ONEM;      *pr = (P)ONEM;
         else {    else {
                 g1 = p1; g2 = p2;      g1 = p1; g2 = p2;
                 if ( ( v1 = VR(g1) ) == ( v2 = VR(g2) ) ) {      if ( ( v1 = VR(g1) ) == ( v2 = VR(g2) ) ) {
                         gcdcmp(vl,mod,g1,&gc1); divsmp(vl,mod,g1,gc1,&gp1);        gcdcmp(vl,mod,g1,&gc1); divsmp(vl,mod,g1,gc1,&gp1);
                         gcdcmp(vl,mod,g2,&gc2); divsmp(vl,mod,g2,gc2,&gp2);        gcdcmp(vl,mod,g2,&gc2); divsmp(vl,mod,g2,gc2,&gp2);
                         gcdprsmp(vl,mod,gc1,gc2,&gcr);        gcdprsmp(vl,mod,gc1,gc2,&gcr);
                         sprsm(vl,mod,v1,gp1,gp2,&g);        sprsm(vl,mod,v1,gp1,gp2,&g);
   
                         if ( VR(g) == v1 ) {        if ( VR(g) == v1 ) {
                                 gp = g;          gp = g;
                                 gcdcmp(vl,mod,gp,&gc); divsmp(vl,mod,gp,gc,&gp1);          gcdcmp(vl,mod,gp,&gc); divsmp(vl,mod,gp,gc,&gp1);
                                 mulmp(vl,mod,gp1,gcr,pr);          mulmp(vl,mod,gp1,gcr,pr);
                         } else        } else
                                 *pr = gcr;          *pr = gcr;
                 } else {      } else {
                         while ( v1 != vl->v && v2 != vl->v )        while ( v1 != vl->v && v2 != vl->v )
                                 vl = NEXT(vl);          vl = NEXT(vl);
                         if ( v1 == vl->v ) {        if ( v1 == vl->v ) {
                                 gcdcmp(vl,mod,g1,&gc1); gcdprsmp(vl,mod,gc1,g2,pr);          gcdcmp(vl,mod,g1,&gc1); gcdprsmp(vl,mod,gc1,g2,pr);
                         } else {        } else {
                                 gcdcmp(vl,mod,g2,&gc2); gcdprsmp(vl,mod,gc2,g1,pr);          gcdcmp(vl,mod,g2,&gc2); gcdprsmp(vl,mod,gc2,g1,pr);
                         }        }
                 }      }
         }    }
 }  }
   
 void gcdcmp(vl,mod,p,pr)  void gcdcmp(vl,mod,p,pr)
Line 659  VL vl;
Line 659  VL vl;
 int mod;  int mod;
 P p,*pr;  P p,*pr;
 {  {
         P g,g1;    P g,g1;
         DCP dc;    DCP dc;
   
         if ( NUM(p) )    if ( NUM(p) )
                 *pr = (P)ONEM;      *pr = (P)ONEM;
         else {    else {
                 dc = DC(p);      dc = DC(p);
                 g = COEF(dc);      g = COEF(dc);
                 for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {      for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
                         gcdprsmp(vl,mod,g,COEF(dc),&g1);        gcdprsmp(vl,mod,g,COEF(dc),&g1);
                         g = g1;        g = g1;
                 }      }
                 *pr = g;      *pr = g;
         }    }
 }  }
   
 void sprsm(vl,mod,v,p1,p2,pr)  void sprsm(vl,mod,v,p1,p2,pr)
Line 681  int mod;
Line 681  int mod;
 V v;  V v;
 P p1,p2,*pr;  P p1,p2,*pr;
 {  {
         P q1,q2,m,m1,m2,x,h,r,g1,g2;    P q1,q2,m,m1,m2,x,h,r,g1,g2;
         int d;    int d;
         Q dq;    Q dq;
         VL nvl;    VL nvl;
   
         reordvar(vl,v,&nvl);    reordvar(vl,v,&nvl);
         reordermp(nvl,mod,vl,p1,&q1); reordermp(nvl,mod,vl,p2,&q2);    reordermp(nvl,mod,vl,p1,&q1); reordermp(nvl,mod,vl,p2,&q2);
   
         if ( ( VR(q1) != v ) || ( VR(q2) != v ) ) {    if ( ( VR(q1) != v ) || ( VR(q2) != v ) ) {
                 *pr = 0;      *pr = 0;
                 return;      return;
         }    }
   
         if ( deg(v,q1) >= deg(v,q2) ) {    if ( deg(v,q1) >= deg(v,q2) ) {
                 g1 = q1; g2 = q2;      g1 = q1; g2 = q2;
         } else {    } else {
                 g2 = q1; g1 = q2;      g2 = q1; g1 = q2;
         }    }
   
         for ( h = (P)ONEM, x = (P)ONEM; ; ) {    for ( h = (P)ONEM, x = (P)ONEM; ; ) {
                 if ( !deg(v,g2) )      if ( !deg(v,g2) )
                         break;        break;
   
                 premmp(nvl,mod,g1,g2,&r);      premmp(nvl,mod,g1,g2,&r);
                 if ( !r )      if ( !r )
                         break;        break;
   
                 d = deg(v,g1) - deg(v,g2); STOQ(d,dq);      d = deg(v,g1) - deg(v,g2); STOQ(d,dq);
                 pwrmp(nvl,mod,h,dq,&m); mulmp(nvl,mod,m,x,&m1); g1 = g2;      pwrmp(nvl,mod,h,dq,&m); mulmp(nvl,mod,m,x,&m1); g1 = g2;
                 divsmp(nvl,mod,r,m1,&g2); x = LC(g1); /* g1 is not const w.r.t v */      divsmp(nvl,mod,r,m1,&g2); x = LC(g1); /* g1 is not const w.r.t v */
                 pwrmp(nvl,mod,x,dq,&m1); mulmp(nvl,mod,m1,h,&m2);      pwrmp(nvl,mod,x,dq,&m1); mulmp(nvl,mod,m1,h,&m2);
                 divsmp(nvl,mod,m2,m,&h);      divsmp(nvl,mod,m2,m,&h);
         }    }
         *pr = g2;    *pr = g2;
 }  }
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.5
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>