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

Diff for /OpenXM_contrib2/asir2000/engine/Z.c between version 1.1 and 1.16

version 1.1, 2003/09/02 07:00:51 version 1.16, 2018/03/29 01:32:51
Line 1 
Line 1 
 #include "ca.h"  #include "ca.h"
   #include "base.h"
 #include "inline.h"  #include "inline.h"
   
 typedef struct oZ {  #if defined(__GNUC__)
         int p;  #define INLINE static inline
         unsigned int b[1];  #elif defined(VISUAL)
 } *Z;  #define INLINE __inline
   #else
   #define INLINE
   #endif
   
 #define ZALLOC(d) ((Z)MALLOC_ATOMIC(TRUESIZE(oZ,(d)-1,int)))  INLINE void _addz(Z n1,Z n2,Z nr);
 #define SL(n) ((n)->p)  INLINE void _subz(Z n1,Z n2,Z nr);
   INLINE void _mulz(Z n1,Z n2,Z nr);
   int _addz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr);
   int _subz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr);
   
 Z qtoz(Q n);  /* immediate int -> Z */
 Q ztoq(Z n);  #define UTOZ(c,n) (n)=(!((unsigned int)(c))?0:(((unsigned int)(c))<=IMM_MAX?((Z)((((unsigned int)(c))<<1)|1)):utoz((unsigned int)(c))))
 Z chsgnz(Z n);  #define STOZ(c,n) (n)=(!((int)(c))?0:(((int)(c))>=IMM_MIN&&((int)(c))<=IMM_MAX?((Z)((((int)(c))<<1)|1)):stoz((int)(c))))
 Z dupz(Z n);  /* immediate Z ? */
 Z addz(Z n1,Z n2);  #define IS_IMM(c) (((unsigned int)c)&1)
 Z subz(Z n1,Z n2);  /* Z can be conver to immediate ? */
 Z mulz(Z n1,Z n2);  #define IS_SZ(n) (((SL(n) == 1)||(SL(n)==-1))&&BD(n)[0]<=IMM_MAX)
 Z divsz(Z n1,Z n2);  /* Z -> immediate Z */
 Z divz(Z n1,Z n2,Z *rem);  #define SZTOZ(n,z) (z)=(Z)(SL(n)<0?(((-BD(n)[0])<<1)|1):(((BD(n)[0])<<1)|1))
 Z gcdz(Z n1,Z n2);  /* Z -> immediate int */
 Z gcdz_cofactor(Z n1,Z n2,Z *c1,Z *c2);  #define ZTOS(c) (((int)(c))>>1)
 Z estimate_array_gcdz(Z *a,int n);  
 Z array_gcdz(Z *a,int n);  
 int remzi(Z n,int m);  
 inline void _addz(Z n1,Z n2,Z nr);  
 inline void _subz(Z n1,Z n2,Z nr);  
 inline void _mulz(Z n1,Z n2,Z nr);  
 inline int _addz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr);  
 inline int _subz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr);  
   
   int uniz(Z a)
   {
     if ( IS_IMM(a) && ZTOS(a) == 1 ) return 1;
     else return 0;
   }
   
   int cmpz(Z a,Z b)
   {
     int *ma,*mb;
     int sa,sb,da,db,ca,cb,i;
   
     if ( !a )
       return -sgnz(b);
     else if ( !b )
       return sgnz(a);
     else {
       sa = sgnz(a); sb = sgnz(b);
       if ( sa > sb ) return 1;
       else if ( sa < sb ) return -1;
       else if ( IS_IMM(a) )
         if ( IS_IMM(b) ) {
           ca = ZTOS(a); cb = ZTOS(b);
           if ( ca > cb ) return sa;
           else if ( ca < cb ) return -sa;
           else return 0;
         } else
           return -sa;
       else if ( IS_IMM(b) )
         return sa;
       else {
         da = SL(a)*sa; db = SL(b)*sa;
         if ( da > db ) return sa;
         else if ( da < db ) return -sa;
         else {
           for ( i = da-1, ma = BD(a)+i, mb = BD(b)+i; i >= 0; i-- )
             if ( *ma > *mb ) return sa;
             else if ( *ma < *mb ) return -sa;
           return 0;
         }
       }
     }
   }
   
   Z stoz(int c)
   {
     Z z;
   
     z = ZALLOC(1);
     if ( c < 0 ) {
       SL(z) = -1; BD(z)[0] = -c;
     } else {
       SL(z) = 1; BD(z)[0] = c;
     }
     return z;
   }
   
   Z utoz(unsigned int c)
   {
     Z z;
   
     z = ZALLOC(1);
     SL(z) = 1; BD(z)[0] = c;
     return z;
   }
   
   Z simpz(Z n)
   {
     Z n1;
   
     if ( !n ) return 0;
     else if ( IS_IMM(n) ) return n;
     else if ( IS_SZ(n) ) {
       SZTOZ(n,n1); return n1;
     } else
       return n;
   }
   
   int sgnz(Z n)
   {
     if ( !n ) return 0;
     else if ( IS_IMM(n) ) return ZTOS(n)>0?1:-1;
     else if ( SL(n) < 0 ) return -1;
     else return 1;
   }
   
 z_mag(Z n)  z_mag(Z n)
 {  {
         return n_bits((N)n);    int c,i;
   
     if ( !n ) return 0;
     else if ( IS_IMM(n) ) {
       c = ZTOS(n);
       if ( c < 0 ) c = -c;
       for ( i = 0; c; c >>= 1, i++ );
       return i;
     }
     else return n_bits((N)n);
 }  }
   
 Z qtoz(Q n)  Z qtoz(Q n)
 {  {
         Z r;    Z r,t;
     int c;
   
         if ( !n ) return 0;    if ( !n ) return 0;
         else if ( !INT(n) )    else if ( !INT(n) )
                 error("qtoz : invalid input");      error("qtoz : invalid input");
         else {    else {
                 r = dupz((Z)NM(n));      t = (Z)NM(n);
                 if ( SGN(n) < 0 ) SL(r) = -SL(r);      if ( IS_SZ(t) ) {
                 return r;        c = SGN(n) < 0 ? -BD(t)[0] : BD(t)[0];
         }        STOZ(c,r);
       } else {
         r = dupz((Z)t);
         if ( SGN(n) < 0 ) SL(r) = -SL(r);
       }
       return r;
     }
 }  }
   
 Q ztoq(Z n)  Q ztoq(Z n)
 {  {
         Q r;    Q r;
         Z nm;    Z nm;
         int sgn;    int sgn,c;
   
         if ( !n ) return 0;    if ( !n ) return 0;
         else {    else if ( IS_IMM(n) ) {
                 nm = dupz(n);      c = ZTOS(n);
                 if ( SL(nm) < 0 ) {      STOQ(c,r);
                         sgn = -1;      return r;
                         SL(nm) = -SL(nm);    } else {
                 } else      nm = dupz(n);
                         sgn = 1;      if ( SL(nm) < 0 ) {
                 NTOQ((N)nm,sgn,r);        sgn = -1;
                 return r;        SL(nm) = -SL(nm);
         }      } else
         sgn = 1;
       NTOQ((N)nm,sgn,r);
       return r;
     }
 }  }
   
 Z dupz(Z n)  Z dupz(Z n)
 {  {
         Z nr;    Z r;
         int sd,i;    int sd,i;
   
         if ( !n ) return 0;    if ( !n ) return 0;
         if ( (sd = SL(n)) < 0 ) sd = -sd;    else if ( IS_IMM(n) ) return n;
         nr = ZALLOC(sd);    else {
         SL(nr) = SL(n);      if ( (sd = SL(n)) < 0 ) sd = -sd;
         for ( i = 0; i < sd; i++ ) BD(nr)[i] = BD(n)[i];      r = ZALLOC(sd);
         return nr;      SL(r) = SL(n);
       for ( i = 0; i < sd; i++ ) BD(r)[i] = BD(n)[i];
       return r;
     }
 }  }
   
 Z chsgnz(Z n)  Z chsgnz(Z n)
 {  {
         Z r;    Z r;
     int c;
   
         if ( !n ) return 0;    if ( !n ) return 0;
         else {    else if ( IS_IMM(n) ) {
                 r = dupz(n);      c = -ZTOS(n);
                 SL(r) = -SL(r);      STOZ(c,r);
                 return r;      return r;
         }    } else {
       r = dupz(n);
       SL(r) = -SL(r);
       return r;
     }
 }  }
   
   Z absz(Z n)
   {
     Z r;
     int c;
   
     if ( !n ) return 0;
     else if ( sgnz(n) > 0 )
       return n;
     else
       return chsgnz(n);
   }
   
 Z addz(Z n1,Z n2)  Z addz(Z n1,Z n2)
 {  {
         unsigned int u1,u2,t;    int d1,d2,d,c;
         int sd1,sd2,d,d1,d2;    Z r,r1;
         Z nr;    struct oZ t;
   
         if ( !n1 ) return dupz(n2);    if ( !n1 ) return dupz(n2);
         else if ( !n2 ) return dupz(n1);    else if ( !n2 ) return dupz(n1);
         else {    else if ( IS_IMM(n1) ) {
                 d1 = ((sd1 = SL(n1))< 0)?-sd1:sd1;      if ( IS_IMM(n2) ) {
                 d2 = ((sd2 = SL(n2))< 0)?-sd2:sd2;        c = ZTOS(n1)+ZTOS(n2);
                 if ( d1 == 1 && d2 == 1 ) {        STOZ(c,r);
                         u1 = BD(n1)[0]; u2 = BD(n2)[0];      } else {
                         if ( sd1 > 0 ) {        c = ZTOS(n1);
                                 if ( sd2 > 0 ) {        if ( c < 0 ) {
                                         t = u1+u2;          t.p = -1; t.b[0] = -c;
                                         if ( t < u1 ) {        } else {
                                                 nr=ZALLOC(2); SL(nr) = 2; BD(nr)[1] = 1;          t.p = 1; t.b[0] = c;
                                         } else {        }
                                                 nr=ZALLOC(1); SL(nr) = 1;        if ( (d2 = SL(n2)) < 0 ) d2 = -d2;
                                         }        r = ZALLOC(d2+1);
                                         BD(nr)[0] = t;        _addz(&t,n2,r);
                                 } else {        if ( !SL(r) ) r = 0;
                                         if ( u1 == u2 ) nr = 0;      }
                                         else if ( u1 > u2 ) {    } else if ( IS_IMM(n2) ) {
                                                 nr=ZALLOC(1); SL(nr) = 1; BD(nr)[0] = u1-u2;      c = ZTOS(n2);
                                         } else {      if ( c < 0 ) {
                                                 nr=ZALLOC(1); SL(nr) = -1; BD(nr)[0] = u2-u1;        t.p = -1; t.b[0] = -c;
                                         }      } else {
                                 }        t.p = 1; t.b[0] = c;
                         } else {      }
                                 if ( sd2 > 0 ) {      if ( (d1 = SL(n1)) < 0 ) d1 = -d1;
                                         if ( u2 == u1 ) nr = 0;      r = ZALLOC(d1+1);
                                         else if ( u2 > u1 ) {      _addz(n1,&t,r);
                                                 nr=ZALLOC(1); SL(nr) = 1; BD(nr)[0] = u2-u1;      if ( !SL(r) ) r = 0;
                                         } else {    } else {
                                                 nr=ZALLOC(1); SL(nr) = -1; BD(nr)[0] = u1-u2;      if ( (d1 = SL(n1)) < 0 ) d1 = -d1;
                                         }      if ( (d2 = SL(n2)) < 0 ) d2 = -d2;
                                 } else {      d = MAX(d1,d2)+1;
                                         t = u1+u2;      r = ZALLOC(d);
                                         if ( t < u1 ) {      _addz(n1,n2,r);
                                                 nr=ZALLOC(2); SL(nr) = -2; BD(nr)[1] = 1;      if ( !SL(r) ) r = 0;
                                         } else {    }
                                                 nr=ZALLOC(1); SL(nr) = -1;    if ( r && !((int)r&1) && IS_SZ(r) ) {
                                         }      SZTOZ(r,r1); r = r1;
                                         BD(nr)[0] = t;    }
                                 }    return r;
                         }  
                 } else {  
                         d = MAX(d1,d2)+1;  
                         nr = ZALLOC(d);  
                         _addz(n1,n2,nr);  
                         if ( !SL(nr) ) nr = 0;  
                 }  
                 return nr;  
         }  
 }  }
   
 Z subz(Z n1,Z n2)  Z subz(Z n1,Z n2)
 {  {
         int sd1,sd2,d;    int d1,d2,d,c;
         Z nr;    Z r,r1;
     struct oZ t;
   
         if ( !n1 )    if ( !n1 )
                 if ( !n2 ) return 0;      return chsgnz(n2);
                 else {    else if ( !n2 ) return n1;
                         nr = dupz(n2);    else if ( IS_IMM(n1) ) {
                         SL(nr) = -SL(nr);      if ( IS_IMM(n2) ) {
                 }        c = ZTOS(n1)-ZTOS(n2);
         else if ( !n2 ) return dupz(n1);        STOZ(c,r);
         else {      } else {
                 if ( (sd1 = SL(n1)) < 0 ) sd1 = -sd1;        c = ZTOS(n1);
                 if ( (sd2 = SL(n2)) < 0 ) sd2 = -sd2;        if ( c < 0 ) {
                 d = MAX(sd1,sd2)+1;          t.p = -1; t.b[0] = -c;
                 nr = ZALLOC(d);        } else {
                 _subz(n1,n2,nr);          t.p = 1; t.b[0] = c;
                 if ( !SL(nr) ) nr = 0;        }
                 return nr;        if ( (d2 = SL(n2)) < 0 ) d2 = -d2;
         }        r = ZALLOC(d2+1);
         _subz(&t,n2,r);
         if ( !SL(r) ) r = 0;
       }
     } else if ( IS_IMM(n2) ) {
       c = ZTOS(n2);
       if ( c < 0 ) {
         t.p = -1; t.b[0] = -c;
       } else {
         t.p = 1; t.b[0] = c;
       }
       if ( (d1 = SL(n1)) < 0 ) d1 = -d1;
       r = ZALLOC(d1+1);
       _subz(n1,&t,r);
       if ( !SL(r) ) r = 0;
     } else {
       if ( (d1 = SL(n1)) < 0 ) d1 = -d1;
       if ( (d2 = SL(n2)) < 0 ) d2 = -d2;
       d = MAX(d1,d2)+1;
       r = ZALLOC(d);
       _subz(n1,n2,r);
       if ( !SL(r) ) r = 0;
     }
     if ( r && !((int)r&1) && IS_SZ(r) ) {
       SZTOZ(r,r1); r = r1;
     }
     return r;
 }  }
   
 Z mulz(Z n1,Z n2)  Z mulz(Z n1,Z n2)
 {  {
         int sd1,sd2,d1,d2;    int d1,d2,sgn,i;
         unsigned int u1,u2,u,l;    int c1,c2;
         Z nr;    unsigned int u1,u2,u,l;
     Z r;
   
         if ( !n1 || !n2 ) return 0;    if ( !n1 || !n2 ) return 0;
         else {  
                 d1 = ((sd1 = SL(n1))< 0)?-sd1:sd1;    if ( IS_IMM(n1) ) {
                 d2 = ((sd2 = SL(n2))< 0)?-sd2:sd2;      c1 = ZTOS(n1);
                 if ( d1 == 1 && d2 == 1 ) {      if ( IS_IMM(n2) ) {
                         u1 = BD(n1)[0]; u2 = BD(n2)[0];        c2 = ZTOS(n2);
                         DM(u1,u2,u,l);        if ( c1 == 1 )
                         if ( u ) {          return n2;
                                 nr = ZALLOC(2); SL(nr) = sd1*sd2*2; BD(nr)[1] = u;        else if ( c1 == -1 )
                         } else {          return chsgnz(n2);
                                 nr = ZALLOC(1); SL(nr) = sd1*sd2;        else if ( c2 == 1 )
                         }          return n1;
                         BD(nr)[0] = l;        else if ( c2 == -1 )
                 } else {          return chsgnz(n1);
                         nr = ZALLOC(d1+d2);        else {
                         _mulz(n1,n2,nr);          sgn = 1;
                 }          if ( c1 < 0 ) { c1 = -c1; sgn = -sgn; }
                 return nr;          if ( c2 < 0 ) { c2 = -c2; sgn = -sgn; }
         }          u1 = (unsigned int)c1; u2 = (unsigned int)c2;
           DM(u1,u2,u,l);
           if ( !u ) {
             UTOZ(l,r);
           } else {
             r = ZALLOC(2); SL(r) = 2; BD(r)[1] = u; BD(r)[0] = l;
           }
         }
       } else {
         sgn = 1;
         if ( c1 < 0 ) { c1 = -c1; sgn = -sgn; }
         u1 = (unsigned int)c1;
         if ( (d2 = SL(n2)) < 0 ) { sgn = -sgn; d2 = -d2; }
         r = ZALLOC(d2+1);
         for ( i = d2; i >= 0; i-- ) BD(r)[i] = 0;
         muln_1(BD(n2),d2,u1,BD(r));
         SL(r) = BD(r)[d2]?d2+1:d2;
       }
     } else if ( IS_IMM(n2) ) {
       c2 = ZTOS(n2);
       if ( c2 == 1 )
         return n1;
       else if ( c2 == -1 )
         return chsgnz(n1);
   
       sgn = 1;
       if ( c2 < 0 ) { sgn = -sgn; c2 = -c2; }
       u2 = (unsigned int)c2;
       if ( (d1 = SL(n1)) < 0 ) { sgn = -sgn; d1 = -d1; }
       r = ZALLOC(d1+1);
       for ( i = d1; i >= 0; i-- ) BD(r)[i] = 0;
       muln_1(BD(n1),d1,u2,BD(r));
       SL(r) = BD(r)[d1]?d1+1:d1;
     } else {
       sgn = 1;
       if ( (d1 = SL(n1)) < 0 ) d1 = -d1;
       if ( (d2 = SL(n2)) < 0 ) d2 = -d2;
       r = ZALLOC(d1+d2);
       _mulz(n1,n2,r);
     }
     if ( sgn < 0 ) r = chsgnz(r);
     return r;
 }  }
   
 /* XXX */  /* kokokara */
 Z divz(Z n1,Z n2,Z *rem)  #if 0
   Z divsz(Z n1,Z n2)
 {  {
         int sgn,sd1,sd2;    int sgn,d1,d2;
         N q,r;    Z q;
   
         if ( !n2 ) error("divz : division by 0");    if ( !n2 ) error("divsz : division by 0");
         else if ( !n1 ) {    if ( !n1 ) return 0;
                 *rem = 0; return 0;  
         } else {  
                 sd2 = SL(n2);  
                 if ( sd2 == 1 && BD(n2)[0] == 1 ) {  
                         *rem = 0; return n1;  
                 } else if ( sd2 == -1 && BD(n2)[0] == 1 ) {  
                         *rem = 0; return chsgnz(n1);  
                 }  
   
                 sd1 = SL(n1);    if ( IS_IMM(n1) ) {
                 n1 = dupz(n1);      if ( !IS_IMM(n2) )
                 if ( SL(n1) < 0 ) SL(n1) = -SL(n1);        error("divsz : cannot happen");
                 divn((N)n1,(N)n2,&q,&r);      c = ZTOS(n1)/ZTOS(n2);
                 if ( q && ((sd1>0&&sd2<0)||(sd1<0&&sd2>0)) ) SL((Z)q) = -SL((Z)q);      STOZ(c,q);
                 if ( r && sd1 < 0 ) SL((Z)r) = -SL((Z)r);      return q;
                 *rem = (Z)r;    }
                 return (Z)q;    if ( IS_IMM(n2) ) {
         }      sgn = 1;
       u2 = ZTOS(n2); if ( u2 < 0 ) { sgn = -sgn; u2 = -u2; }
       diviz(n1,u2,&q);
       if ( sgn < 0 ) SL(q) = -SL(q);
       return q;
     }
   
     sgn = 1;
     if ( (d2 = SL(n2)) < 0 ) { sgn = -sgn; d2 = -d2; }
     if ( d2 == 1 ) {
       diviz(n1,BD(u2)[0],&q);
       if ( sgn < 0 ) SL(q) = -SL(q);
       return q;
     }
     if ( (d1 = SL(n1)) < 0 ) { sgn = -sgn; d1 = -d1; }
     if ( d1 < d2 ) error("divsz : cannot happen");
     return q;
 }  }
   #endif
   
   /* XXX */
   Z divz(Z n1,Z n2,Z *r)
   {
     int s1,s2;
     Q t1,t2,qq,rq;
     N qn,rn;
   
     if ( !n1 ) {
       *r = 0; return 0;
     }
     if ( !n2 )
       error("divz : division by 0");
     t1 = ztoq(n1); t2 = ztoq(n2);
     s1 = sgnz(n1); s2 = sgnz(n2);
     /* n1 = qn*SGN(n1)*SGN(n2)*n2+SGN(n1)*rn */
     divn(NM(t1),NM(t2),&qn,&rn);
     NTOQ(qn,s1*s2,qq);
     NTOQ(rn,s1,rq);
     *r = qtoz(rq);
     return qtoz(qq);
   }
   
 Z divsz(Z n1,Z n2)  Z divsz(Z n1,Z n2)
 {  {
         int sgn,sd1,sd2;    int s1,s2;
         N q;    Q t1,t2,qq;
     N qn;
   
         if ( !n2 ) error("divsz : division by 0");    if ( !n1 )
         else if ( !n1 ) return 0;      return 0;
         else {    if ( !n2 )
                 sd2 = SL(n2);      error("divsz : division by 0");
                 if ( sd2 == 1 && BD(n2)[0] == 1 ) return n1;    t1 = ztoq(n1); t2 = ztoq(n2);
                 else if ( sd2 == -1 && BD(n2)[0] == 1 ) return chsgnz(n1);    s1 = sgnz(n1); s2 = sgnz(n2);
     /* n1 = qn*SGN(n1)*SGN(n2)*n2 */
                 sd1 = SL(n1);    divsn(NM(t1),NM(t2),&qn);
                 n1 = dupz(n1);    NTOQ(qn,s1*s2,qq);
                 if ( SL(n1) < 0 ) SL(n1) = -SL(n1);    return qtoz(qq);
                 divsn((N)n1,(N)n2,&q);  
                 if ( q && ((sd1>0&&sd2<0)||(sd1<0&&sd2>0)) ) SL((Z)q) = -SL((Z)q);  
                 return (Z)q;  
         }  
 }  }
   
   int gcdimm(int c1,int c2)
   {
     int r;
   
     if ( !c1 ) return c2;
     else if ( !c2 ) return c1;
     while ( 1 ) {
       r = c1%c2;
       if ( !r ) return c2;
       c1 = c2; c2 = r;
     }
   }
   
 Z gcdz(Z n1,Z n2)  Z gcdz(Z n1,Z n2)
 {  {
         int sd1,sd2;    int c1,c2,g;
         N gcd;    Z gcd,r;
     N gn;
   
         if ( !n1 ) return n2;    if ( !n1 ) return n2;
         else if ( !n2 ) return n1;    else if ( !n2 ) return n1;
   
         n1 = dupz(n1);    if ( IS_IMM(n1) ) {
         if ( SL(n1) < 0 ) SL(n1) = -SL(n1);      c1 = ZTOS(n1);
         n2 = dupz(n2);      if ( c1 < 0 ) c1 = -c1;
         if ( SL(n2) < 0 ) SL(n2) = -SL(n2);      if ( IS_IMM(n2) )
         gcdn((N)n1,(N)n2,&gcd);        c2 = ZTOS(n2);
         return (Z)gcd;      else
         c2 = remzi(n2,c1>0?c1:-c1);
       if ( c2 < 0 ) c2 = -c2;
       g = gcdimm(c1,c2);
       STOZ(g,gcd);
       return gcd;
     } else if ( IS_IMM(n2) ) {
       c2 = ZTOS(n2);
       if ( c2 < 0 ) c2 = -c2;
       c1 = remzi(n1,c2>0?c2:-c2);
       if ( c1 < 0 ) c1 = -c1;
       g = gcdimm(c1,c2);
       STOZ(g,gcd);
       return gcd;
     } else {
       n1 = dupz(n1);
       if ( SL(n1) < 0 ) SL(n1) = -SL(n1);
       n2 = dupz(n2);
       if ( SL(n2) < 0 ) SL(n2) = -SL(n2);
       gcdn((N)n1,(N)n2,&gn);
       gcd = (Z)gn;
       if ( IS_SZ(gcd) ) {
         SZTOZ(gcd,r); gcd = r;
       }
       return gcd;
     }
 }  }
   
 int remzi(Z n,int m)  int remzi(Z n,int m)
 {  {
         unsigned int *x;    unsigned int *x;
         unsigned int t,r;    unsigned int t,r;
         int i;    int i,c;
   
         if ( !n ) return 0;    if ( !n ) return 0;
         i = SL(n);    if ( IS_IMM(n) ) {
         if ( i < 0 ) i = -i;      c = ZTOS(n)%m;
         for ( i--, x = BD(n)+i, r = 0; i >= 0; i--, x-- ) {      if ( c < 0 ) c += m;
       return c;
     }
   
     i = SL(n);
     if ( i < 0 ) i = -i;
     for ( i--, x = BD(n)+i, r = 0; i >= 0; i--, x-- ) {
 #if defined(sparc)  #if defined(sparc)
                 r = dsa(m,r,*x);      r = dsa(m,r,*x);
 #else  #else
                 DSAB(m,r,*x,t,r);      DSAB(m,r,*x,t,r);
 #endif  #endif
         }    }
         return r;    if ( r && SL(n) < 0 )
       r = m-r;
     return r;
 }  }
   
 Z gcdz_cofactor(Z n1,Z n2,Z *c1,Z *c2)  Z gcdz_cofactor(Z n1,Z n2,Z *c1,Z *c2)
 {  {
         Z gcd;    Z gcd;
   
         gcd = gcdz(n1,n2);    gcd = gcdz(n1,n2);
         *c1 = divsz(n1,gcd);    *c1 = divsz(n1,gcd);
         *c2 = divsz(n2,gcd);    *c2 = divsz(n2,gcd);
         return gcd;    return gcd;
 }  }
   
   #if 0
 Z estimate_array_gcdz(Z *b,int n)  Z estimate_array_gcdz(Z *b,int n)
 {  {
         int m,i,j,sd;    int m,i,j,sd;
         Z *a;    Z *a;
         Z s,t;    Z s,t;
   
         a = (Z *)ALLOCA(n*sizeof(Z));    a = (Z *)ALLOCA(n*sizeof(Z));
         for ( i = j = 0; i < n; i++ ) if ( b[i] ) a[j++] = b[i];    for ( i = j = 0; i < n; i++ ) if ( b[i] ) a[j++] = b[i];
         n = j;    n = j;
         if ( !n ) return 0;    if ( !n ) return 0;
         if ( n == 1 ) return a[0];    if ( n == 1 ) return a[0];
   
         m = n/2;    m = n/2;
         for ( i = 0, s = 0; i < m; i++ ) {    for ( i = 0, s = 0; i < m; i++ ) {
                 if ( !a[i] ) continue;      if ( !a[i] ) continue;
                 else s = (SL(a[i])<0)?subz(s,a[i]):addz(s,a[i]);      else s = (SL(a[i])<0)?subz(s,a[i]):addz(s,a[i]);
         }    }
         for ( t = 0; i < n; i++ ) {    for ( t = 0; i < n; i++ ) {
                 if ( !a[i] ) continue;      if ( !a[i] ) continue;
                 else t = (SL(a[i])<0)?subz(t,a[i]):addz(t,a[i]);      else t = (SL(a[i])<0)?subz(t,a[i]):addz(t,a[i]);
         }    }
         return gcdz(s,t);    return gcdz(s,t);
 }  }
   
 Z array_gcdz(Z *b,int n)  Z array_gcdz(Z *b,int n)
 {  {
         int m,i,j,sd;    int m,i,j,sd;
         Z *a;    Z *a;
         Z gcd;    Z gcd;
   
         a = (Z *)ALLOCA(n*sizeof(Z));    a = (Z *)ALLOCA(n*sizeof(Z));
         for ( i = j = 0; i < n; i++ ) if ( b[i] ) a[j++] = b[i];    for ( i = j = 0; i < n; i++ ) if ( b[i] ) a[j++] = b[i];
         n = j;    n = j;
         if ( !n ) return 0;    if ( !n ) return 0;
         if ( n == 1 ) return a[0];    if ( n == 1 ) return a[0];
         gcd = a[0];    gcd = a[0];
         for ( i = 1; i < n; i++ )    for ( i = 1; i < n; i++ )
                 gcd = gcdz(gcd,a[i]);      gcd = gcdz(gcd,a[i]);
         return gcd;    return gcd;
 }  }
   #endif
   
 void _copyz(Z n1,Z n2)  void _copyz(Z n1,Z n2)
 {  {
         int n,i;    int n,i;
   
         if ( !n1 || !SL(n1) ) SL(n2) = 0;    if ( !n1 || !SL(n1) ) SL(n2) = 0;
         else {    else {
                 n = SL(n2) = SL(n1);      n = SL(n2) = SL(n1);
                 if ( n < 0 ) n = -n;      if ( n < 0 ) n = -n;
                 for ( i = 0; i < n; i++ ) BD(n2)[i] = BD(n1)[i];      for ( i = 0; i < n; i++ ) BD(n2)[i] = BD(n1)[i];
         }    }
 }  }
   
 void _addz(Z n1,Z n2,Z nr)  void _addz(Z n1,Z n2,Z nr)
 {  {
         int sd1,sd2;    int d1,d2;
   
         if ( !n1 || !SL(n1) ) _copyz(n2,nr);    if ( !n1 || !SL(n1) ) _copyz(n2,nr);
         else if ( !n2 || !SL(n2) ) _copyz(n1,nr);    else if ( !n2 || !SL(n2) ) _copyz(n1,nr);
         else if ( (sd1=SL(n1)) > 0 )    else if ( (d1=SL(n1)) > 0 )
                 if ( (sd2=SL(n2)) > 0 )      if ( (d2=SL(n2)) > 0 )
                         SL(nr) = _addz_main(BD(n1),sd1,BD(n2),sd2,BD(nr));        SL(nr) = _addz_main(BD(n1),d1,BD(n2),d2,BD(nr));
                 else      else
                         SL(nr) = _subz_main(BD(n1),sd1,BD(n2),-sd2,BD(nr));        SL(nr) = _subz_main(BD(n1),d1,BD(n2),-d2,BD(nr));
         else if ( (sd2=SL(n2)) > 0 )    else if ( (d2=SL(n2)) > 0 )
                 SL(nr) = _subz_main(BD(n2),sd2,BD(n1),-sd1,BD(nr));      SL(nr) = _subz_main(BD(n2),d2,BD(n1),-d1,BD(nr));
         else    else
                 SL(nr) = -_addz_main(BD(n1),-sd1,BD(n2),-sd2,BD(nr));      SL(nr) = -_addz_main(BD(n1),-d1,BD(n2),-d2,BD(nr));
 }  }
   
 void _subz(Z n1,Z n2,Z nr)  void _subz(Z n1,Z n2,Z nr)
 {  {
         int sd1,sd2;    int d1,d2;
   
         if ( !n1 || !SL(n1) ) _copyz(n2,nr);    if ( !n1 || !SL(n1) ) _copyz(n2,nr);
         else if ( !n2 || !SL(n2) ) {    else if ( !n2 || !SL(n2) ) {
                 _copyz(n1,nr);      _copyz(n1,nr);
                 SL(nr) = -SL(nr);      SL(nr) = -SL(nr);
         } else if ( (sd1=SL(n1)) > 0 )    } else if ( (d1=SL(n1)) > 0 )
                 if ( (sd2=SL(n2)) > 0 )      if ( (d2=SL(n2)) > 0 )
                         SL(nr) = _subz_main(BD(n1),sd1,BD(n2),sd2,BD(nr));        SL(nr) = _subz_main(BD(n1),d1,BD(n2),d2,BD(nr));
                 else      else
                         SL(nr) = _addz_main(BD(n1),sd1,BD(n2),-sd2,BD(nr));        SL(nr) = _addz_main(BD(n1),d1,BD(n2),-d2,BD(nr));
         else if ( (sd2=SL(n2)) > 0 )    else if ( (d2=SL(n2)) > 0 )
                 SL(nr) = -_addz_main(BD(n1),-sd1,BD(n2),sd2,BD(nr));      SL(nr) = -_addz_main(BD(n1),-d1,BD(n2),d2,BD(nr));
         else    else
                 SL(nr) = -_subz_main(BD(n1),-sd1,BD(n2),-sd2,BD(nr));      SL(nr) = -_subz_main(BD(n1),-d1,BD(n2),-d2,BD(nr));
 }  }
   
 void _mulz(Z n1,Z n2,Z nr)  void _mulz(Z n1,Z n2,Z nr)
 {  {
         int sd1,sd2;    int d1,d2;
         int i,d,sgn;    int i,d,sgn;
         unsigned int mul;    unsigned int mul;
         unsigned int *m1,*m2;    unsigned int *m1,*m2;
   
         if ( !n1 || !SL(n1) || !n2 || !SL(n2) )    if ( !n1 || !SL(n1) || !n2 || !SL(n2) )
                 SL(nr) = 0;      SL(nr) = 0;
         else {    else {
                 sd1 = SL(n1); sd2 = SL(n2);      d1 = SL(n1); d2 = SL(n2);
                 sgn = 1;      sgn = 1;
                 if ( sd1 < 0 ) { sgn = -sgn; sd1 = -sd1; }      if ( d1 < 0 ) { sgn = -sgn; d1 = -d1; }
                 if ( sd2 < 0 ) { sgn = -sgn; sd2 = -sd2; }      if ( d2 < 0 ) { sgn = -sgn; d2 = -d2; }
                 d = sd1+sd2;      d = d1+d2;
                 for ( i = d-1, m1 = BD(nr); i >= 0; i-- ) *m1++ = 0;      for ( i = d-1, m1 = BD(nr); i >= 0; i-- ) *m1++ = 0;
                 for ( i = 0, m1 = BD(n1), m2 = BD(n2); i < sd2; i++, m2++ )      for ( i = 0, m1 = BD(n1), m2 = BD(n2); i < d2; i++, m2++ )
                         if ( mul = *m2 ) muln_1(m1,sd1,mul,BD(nr)+i);        if ( mul = *m2 ) muln_1(m1,d1,mul,BD(nr)+i);
                 SL(nr) = sgn*(BD(nr)[d-1]?d:d-1);      SL(nr) = sgn*(BD(nr)[d-1]?d:d-1);
         }    }
 }  }
   
 int _addz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr)  int _addz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr)
 {  {
         int d,i;    int d,i;
         unsigned int tmp,c;    unsigned int tmp,c;
         unsigned int *t;    unsigned int *t;
   
         if ( d2 > d1 ) {    if ( d2 > d1 ) {
                 t = m1; m1 = m2; m2 = t;      t = m1; m1 = m2; m2 = t;
                 d = d1; d1 = d2; d2 = d;      d = d1; d1 = d2; d2 = d;
         }    }
 #if defined(VISUAL)  #if defined(_M_IX86) && !defined(__MINGW32__)
         __asm {    __asm {
         push    esi    push  esi
         push    edi    push  edi
         mov esi,m1    mov esi,m1
         mov edi,m2    mov edi,m2
         mov ebx,mr    mov ebx,mr
         mov ecx,d2    mov ecx,d2
         xor     eax,eax    xor  eax,eax
         Lstart__addz:    Lstart__addz:
         mov eax,DWORD PTR [esi]    mov eax,DWORD PTR [esi]
         mov edx,DWORD PTR [edi]    mov edx,DWORD PTR [edi]
         adc eax,edx    adc eax,edx
         mov DWORD PTR [ebx],eax    mov DWORD PTR [ebx],eax
         lea esi,DWORD PTR [esi+4]    lea esi,DWORD PTR [esi+4]
         lea edi,DWORD PTR [edi+4]    lea edi,DWORD PTR [edi+4]
         lea ebx,DWORD PTR [ebx+4]    lea ebx,DWORD PTR [ebx+4]
         dec ecx    dec ecx
         jnz Lstart__addz    jnz Lstart__addz
         pop     edi    pop  edi
         pop     esi    pop  esi
         mov eax,0    mov eax,0
         adc eax,eax    adc eax,eax
         mov c,eax    mov c,eax
         }    }
 #elif defined(i386)  #elif defined(i386) && !defined(__MINGW32__)
         asm volatile("\    asm volatile("\
         movl    %1,%%esi;\    pushl  %%ebx;\
         movl    %2,%%edi;\    movl  %1,%%esi;\
         movl    %3,%%ebx;\    movl  %2,%%edi;\
         movl    %4,%%ecx;\    movl  %3,%%ebx;\
         testl   %%eax,%%eax;\    movl  %4,%%ecx;\
         Lstart__addz:;\    testl  %%eax,%%eax;\
         movl    (%%esi),%%eax;\    Lstart__addz:;\
         movl    (%%edi),%%edx;\    movl  (%%esi),%%eax;\
         adcl    %%edx,%%eax;\    movl  (%%edi),%%edx;\
         movl    %%eax,(%%ebx);\    adcl  %%edx,%%eax;\
         leal    4(%%esi),%%esi;\    movl  %%eax,(%%ebx);\
         leal    4(%%edi),%%edi;\    leal  4(%%esi),%%esi;\
         leal    4(%%ebx),%%ebx;\    leal  4(%%edi),%%edi;\
         decl    %%ecx;\    leal  4(%%ebx),%%ebx;\
         jnz Lstart__addz;\    decl  %%ecx;\
         movl    $0,%%eax;\    jnz Lstart__addz;\
         adcl    %%eax,%%eax;\    movl  $0,%%eax;\
         movl    %%eax,%0"\    adcl  %%eax,%%eax;\
         :"=m"(c)\    movl  %%eax,%0;\
         :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\    popl  %%ebx"\
         :"eax","ebx","ecx","edx","esi","edi");    :"=m"(c)\
     :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
     :"eax","ecx","edx","esi","edi");
 #else  #else
         for ( i = 0, c = 0, mr = BD(nr); i < d2; i++, m1++, m2++, mr++ ) {    for ( i = 0, c = 0; i < d2; i++, m1++, m2++, mr++ ) {
                 tmp = *m1 + *m2;      tmp = *m1 + *m2;
                 if ( tmp < *m1 ) {      if ( tmp < *m1 ) {
                         tmp += c;        tmp += c;
                         c = 1;        c = 1;
                 } else {      } else {
                         tmp += c;        tmp += c;
                         c = tmp < c ? 1 : 0;        c = tmp < c ? 1 : 0;
                 }      }
                 *mr = tmp;      *mr = tmp;
         }    }
 #endif  #endif
         for ( i = d2, m1 += d2, mr += d2; (i < d1) && c ; i++ ) {    for ( i = d2, m1 += d2, mr += d2; (i < d1) && c ; i++ ) {
                 tmp = *m1++ + c;      tmp = *m1++ + c;
                 c = tmp < c ? 1 : 0;      c = tmp < c ? 1 : 0;
                 *mr++ = tmp;      *mr++ = tmp;
         }    }
         for ( ; i < d1; i++ )    for ( ; i < d1; i++ )
                         *mr++ = *m1++;        *mr++ = *m1++;
         *mr = c;    *mr = c;
         return (c?d1+1:d1);    return (c?d1+1:d1);
 }  }
   
 int _subz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr)  int _subz_main(unsigned int *m1,int d1,unsigned int *m2,int d2,unsigned int *mr)
 {  {
         int d,i,sgn;    int d,i,sgn;
         unsigned int t,tmp,br;    unsigned int t,tmp,br;
         unsigned int *m;    unsigned int *m;
   
         if ( d1 > d2 ) sgn = 1;    if ( d1 > d2 ) sgn = 1;
         else if ( d1 < d2 ) sgn = -1;    else if ( d1 < d2 ) sgn = -1;
         else {    else {
                 for ( i = d1-1; i >= 0 && m1[i] == m2[i]; i-- );      for ( i = d1-1; i >= 0 && m1[i] == m2[i]; i-- );
                 if ( i < 0 ) return 0;      if ( i < 0 ) return 0;
                 if ( m1[i] > m2[i] ) sgn = 1;      if ( m1[i] > m2[i] ) sgn = 1;
                 else if ( m1[i] < m2[i] ) sgn = -1;      else if ( m1[i] < m2[i] ) sgn = -1;
         }    }
         if ( sgn < 0 ) {    if ( sgn < 0 ) {
                 m = m1; m1 = m2; m2 = m;      m = m1; m1 = m2; m2 = m;
                 d = d1; d1 = d2; d2 = d;      d = d1; d1 = d2; d2 = d;
         }    }
 #if defined(VISUAL)  #if defined(_M_IX86) && !defined(__MINGW32__)
         __asm {    __asm {
         push    esi    push  esi
         push    edi    push  edi
         mov esi,m1    mov esi,m1
         mov edi,m2    mov edi,m2
         mov ebx,mr    mov ebx,mr
         mov ecx,d2    mov ecx,d2
         xor     eax,eax    xor  eax,eax
         Lstart__subz:    Lstart__subz:
         mov eax,DWORD PTR [esi]    mov eax,DWORD PTR [esi]
         mov edx,DWORD PTR [edi]    mov edx,DWORD PTR [edi]
         sbb eax,edx    sbb eax,edx
         mov DWORD PTR [ebx],eax    mov DWORD PTR [ebx],eax
         lea esi,DWORD PTR [esi+4]    lea esi,DWORD PTR [esi+4]
         lea edi,DWORD PTR [edi+4]    lea edi,DWORD PTR [edi+4]
         lea ebx,DWORD PTR [ebx+4]    lea ebx,DWORD PTR [ebx+4]
         dec ecx    dec ecx
         jnz Lstart__subz    jnz Lstart__subz
         pop     edi    pop  edi
         pop     esi    pop  esi
         mov eax,0    mov eax,0
         adc eax,eax    adc eax,eax
         mov br,eax    mov br,eax
         }    }
 #elif defined(i386)  #elif defined(i386) && !defined(__MINGW32__)
         asm volatile("\    asm volatile("\
         movl    %1,%%esi;\    pushl  %%ebx;\
         movl    %2,%%edi;\    movl  %1,%%esi;\
         movl    %3,%%ebx;\    movl  %2,%%edi;\
         movl    %4,%%ecx;\    movl  %3,%%ebx;\
         testl   %%eax,%%eax;\    movl  %4,%%ecx;\
         Lstart__subz:;\    testl  %%eax,%%eax;\
         movl    (%%esi),%%eax;\    Lstart__subz:;\
         movl    (%%edi),%%edx;\    movl  (%%esi),%%eax;\
         sbbl    %%edx,%%eax;\    movl  (%%edi),%%edx;\
         movl    %%eax,(%%ebx);\    sbbl  %%edx,%%eax;\
         leal    4(%%esi),%%esi;\    movl  %%eax,(%%ebx);\
         leal    4(%%edi),%%edi;\    leal  4(%%esi),%%esi;\
         leal    4(%%ebx),%%ebx;\    leal  4(%%edi),%%edi;\
         decl    %%ecx;\    leal  4(%%ebx),%%ebx;\
         jnz Lstart__subz;\    decl  %%ecx;\
         movl    $0,%%eax;\    jnz Lstart__subz;\
         adcl    %%eax,%%eax;\    movl  $0,%%eax;\
         movl    %%eax,%0"\    adcl  %%eax,%%eax;\
         :"=m"(br)\    movl  %%eax,%0;\
         :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\    popl  %%ebx"\
         :"eax","ebx","ecx","edx","esi","edi");    :"=m"(br)\
     :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
     :"eax","ecx","edx","esi","edi");
 #else  #else
         for ( i = 0, br = 0, mr = BD(nr); i < d2; i++, mr++ ) {    for ( i = 0, br = 0; i < d2; i++, mr++ ) {
                 t = *m1++;      t = *m1++;
                 tmp = *m2++ + br;      tmp = *m2++ + br;
                 if ( br > 0 && !tmp ) {      if ( br > 0 && !tmp ) {
                         /* tmp = 2^32 => br = 1 */        /* tmp = 2^32 => br = 1 */
                 }else {      }else {
                         tmp = t-tmp;        tmp = t-tmp;
                         br = tmp > t ? 1 : 0;        br = tmp > t ? 1 : 0;
                         *mr = tmp;        *mr = tmp;
                 }      }
         }    }
 #endif  #endif
         for ( i = d2, m1 += d2, mr += d2; (i < d1) && br; i++ ) {    for ( i = d2, m1 += d2, mr += d2; (i < d1) && br; i++ ) {
                 t = *m1++;      t = *m1++;
                 tmp = t - br;      tmp = t - br;
                 br = tmp > t ? 1 : 0;      br = tmp > t ? 1 : 0;
                 *mr++ = tmp;      *mr++ = tmp;
         }    }
         for ( ; i < d1; i++ )    for ( ; i < d1; i++ )
                 *mr++ = *m1++;      *mr++ = *m1++;
         for ( i = d1-1, mr--; i >= 0 && !*mr--; i-- );    for ( i = d1-1, mr--; i >= 0 && !*mr--; i-- );
         return sgn*(i+1);    return sgn*(i+1);
 }  }
   
 /* XXX */  /* XXX */
   
 void printz(Z n)  void printz(Z n)
 {  {
         int sd;    int sd,u;
   
         if ( !n )    if ( !n )
                 fprintf(asir_out,"0");      fprintf(asir_out,"0");
         else {    else if ( IS_IMM(n) ) {
                 if ( (sd = SL(n)) < 0 ) { SL(n) = -SL(n); fprintf(asir_out,"-"); }      u = ZTOS(n);
                 printn((N)n);      fprintf(asir_out,"%d",u);
                 if ( sd < 0 ) SL(n) = -SL(n);    } else {
         }      if ( (sd = SL(n)) < 0 ) { SL(n) = -SL(n); fprintf(asir_out,"-"); }
       printn((N)n);
       if ( sd < 0 ) SL(n) = -SL(n);
     }
   }
   
   /*
    *  Dx^k*x^l = W(k,l,0)*x^l*Dx^k+W(k,l,1)*x^(l-1)*x^(k-1)*+...
    *
    *  t = [W(k,l,0) W(k,l,1) ... W(k,l,min(k,l)]
    *  where W(k,l,i) = i! * kCi * lCi
    */
   
   void mkwcz(int k,int l,Z *t)
   {
     int i,n,up,low;
     N nm,d,c;
   
     n = MIN(k,l);
     for ( t[0] = (Z)ONEN, i = 1; i <= n; i++ ) {
       DM(k-i+1,l-i+1,up,low);
       if ( up ) {
         nm = NALLOC(2); PL(nm) = 2; BD(nm)[0] = low; BD(nm)[1] = up;
       } else {
         nm = NALLOC(1); PL(nm) = 1; BD(nm)[0] = low;
       }
       kmuln((N)t[i-1],nm,&d); divin(d,i,&c); t[i] = (Z)c;
     }
 }  }
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.16
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>