[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.12

version 1.1, 2003/09/02 07:00:51 version 1.12, 2013/11/05 11:36:58
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) ) {
                           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;                  return r;
         }          }
 }  }
Line 52  Q ztoq(Z n)
Line 152  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) ) {
                   c = ZTOS(n);
                   STOQ(c,r);
                   return r;
           } else {
                 nm = dupz(n);                  nm = dupz(n);
                 if ( SL(nm) < 0 ) {                  if ( SL(nm) < 0 ) {
                         sgn = -1;                          sgn = -1;
Line 69  Q ztoq(Z n)
Line 173  Q ztoq(Z n)
   
 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) ) {
                   c = -ZTOS(n);
                   STOZ(c,r);
                   return r;
           } else {
                 r = dupz(n);                  r = dupz(n);
                 SL(r) = -SL(r);                  SL(r) = -SL(r);
                 return 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 ) {  
                                                 nr=ZALLOC(2); SL(nr) = 2; BD(nr)[1] = 1;  
                                         } else {  
                                                 nr=ZALLOC(1); SL(nr) = 1;  
                                         }  
                                         BD(nr)[0] = t;  
                                 } else {  
                                         if ( u1 == u2 ) nr = 0;  
                                         else if ( u1 > u2 ) {  
                                                 nr=ZALLOC(1); SL(nr) = 1; BD(nr)[0] = u1-u2;  
                                         } else {  
                                                 nr=ZALLOC(1); SL(nr) = -1; BD(nr)[0] = u2-u1;  
                                         }  
                                 }  
                         } else {                          } else {
                                 if ( sd2 > 0 ) {                                  t.p = 1; t.b[0] = c;
                                         if ( u2 == u1 ) nr = 0;  
                                         else if ( u2 > u1 ) {  
                                                 nr=ZALLOC(1); SL(nr) = 1; BD(nr)[0] = u2-u1;  
                                         } else {  
                                                 nr=ZALLOC(1); SL(nr) = -1; BD(nr)[0] = u1-u2;  
                                         }  
                                 } else {  
                                         t = u1+u2;  
                                         if ( t < u1 ) {  
                                                 nr=ZALLOC(2); SL(nr) = -2; BD(nr)[1] = 1;  
                                         } else {  
                                                 nr=ZALLOC(1); SL(nr) = -1;  
                                         }  
                                         BD(nr)[0] = t;  
                                 }  
                         }                          }
                           if ( (d2 = SL(n2)) < 0 ) d2 = -d2;
                           r = ZALLOC(d2+1);
                           _addz(&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 {                  } else {
                         d = MAX(d1,d2)+1;                          t.p = 1; t.b[0] = c;
                         nr = ZALLOC(d);  
                         _addz(n1,n2,nr);  
                         if ( !SL(nr) ) nr = 0;  
                 }                  }
                 return nr;                  if ( (d1 = SL(n1)) < 0 ) d1 = -d1;
                   r = ZALLOC(d1+1);
                   _addz(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);
                   _addz(n1,n2,r);
                   if ( !SL(r) ) r = 0;
         }          }
           if ( r && !((int)r&1) && IS_SZ(r) ) {
                   SZTOZ(r,r1); r = r1;
           }
           return r;
 }  }
   
 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);
                           STOZ(c,r);
                   } else {
                           c = ZTOS(n1);
                           if ( c < 0 ) {
                                   t.p = -1; t.b[0] = -c;
                           } else {
                                   t.p = 1; t.b[0] = c;
                           }
                           if ( (d2 = SL(n2)) < 0 ) d2 = -d2;
                           r = ZALLOC(d2+1);
                           _subz(&t,n2,r);
                           if ( !SL(r) ) r = 0;
                 }                  }
         else if ( !n2 ) return dupz(n1);          } else if ( IS_IMM(n2) ) {
         else {                  c = ZTOS(n2);
                 if ( (sd1 = SL(n1)) < 0 ) sd1 = -sd1;                  if ( c < 0 ) {
                 if ( (sd2 = SL(n2)) < 0 ) sd2 = -sd2;                          t.p = -1; t.b[0] = -c;
                 d = MAX(sd1,sd2)+1;                  } else {
                 nr = ZALLOC(d);                          t.p = 1; t.b[0] = c;
                 _subz(n1,n2,nr);                  }
                 if ( !SL(nr) ) nr = 0;                  if ( (d1 = SL(n1)) < 0 ) d1 = -d1;
                 return nr;                  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;
           int c1,c2;
         unsigned int u1,u2,u,l;          unsigned int u1,u2,u,l;
         Z nr;          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;
                           else if ( c2 == -1 )
                                   return chsgnz(n1);
                           else {
                                   sgn = 1;
                                   if ( c1 < 0 ) { c1 = -c1; sgn = -sgn; }
                                   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;
                                   }
                         }                          }
                         BD(nr)[0] = l;  
                 } else {                  } else {
                         nr = ZALLOC(d1+d2);                          sgn = 1;
                         _mulz(n1,n2,nr);                          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;
                 }                  }
                 return nr;          } 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;
           if ( IS_IMM(n) ) {
                   c = ZTOS(n)%m;
                   if ( c < 0 ) c += m;
                   return c;
           }
   
         i = SL(n);          i = SL(n);
         if ( i < 0 ) i = -i;          if ( i < 0 ) i = -i;
         for ( i--, x = BD(n)+i, r = 0; i >= 0; i--, x-- ) {          for ( i--, x = BD(n)+i, r = 0; i >= 0; i--, x-- ) {
Line 281  int remzi(Z n,int m)
Line 539  int remzi(Z n,int m)
                 DSAB(m,r,*x,t,r);                  DSAB(m,r,*x,t,r);
 #endif  #endif
         }          }
           if ( r && SL(n) < 0 )
                   r = m-r;
         return r;          return r;
 }  }
   
Line 294  Z gcdz_cofactor(Z n1,Z n2,Z *c1,Z *c2)
Line 554  Z gcdz_cofactor(Z n1,Z n2,Z *c1,Z *c2)
         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;
Line 334  Z array_gcdz(Z *b,int n)
Line 595  Z array_gcdz(Z *b,int n)
                 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)
 {  {
Line 349  void _copyz(Z n1,Z n2)
Line 611  void _copyz(Z n1,Z n2)
   
 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;
Line 393  void _mulz(Z n1,Z n2,Z nr)
Line 655  void _mulz(Z n1,Z n2,Z nr)
         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);
         }          }
 }  }
Line 415  int _addz_main(unsigned int *m1,int d1,unsigned int *m
Line 677  int _addz_main(unsigned int *m1,int d1,unsigned int *m
                 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)
         __asm {          __asm {
         push    esi          push    esi
         push    edi          push    edi
Line 442  int _addz_main(unsigned int *m1,int d1,unsigned int *m
Line 704  int _addz_main(unsigned int *m1,int d1,unsigned int *m
         }          }
 #elif defined(i386)  #elif defined(i386)
         asm volatile("\          asm volatile("\
           pushl   %%ebx;\
         movl    %1,%%esi;\          movl    %1,%%esi;\
         movl    %2,%%edi;\          movl    %2,%%edi;\
         movl    %3,%%ebx;\          movl    %3,%%ebx;\
Line 459  int _addz_main(unsigned int *m1,int d1,unsigned int *m
Line 722  int _addz_main(unsigned int *m1,int d1,unsigned int *m
         jnz Lstart__addz;\          jnz Lstart__addz;\
         movl    $0,%%eax;\          movl    $0,%%eax;\
         adcl    %%eax,%%eax;\          adcl    %%eax,%%eax;\
         movl    %%eax,%0"\          movl    %%eax,%0;\
           popl    %%ebx"\
         :"=m"(c)\          :"=m"(c)\
         :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\          :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
         :"eax","ebx","ecx","edx","esi","edi");          :"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;
Line 505  int _subz_main(unsigned int *m1,int d1,unsigned int *m
Line 769  int _subz_main(unsigned int *m1,int d1,unsigned int *m
                 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)
         __asm {          __asm {
         push    esi          push    esi
         push    edi          push    edi
Line 532  int _subz_main(unsigned int *m1,int d1,unsigned int *m
Line 796  int _subz_main(unsigned int *m1,int d1,unsigned int *m
         }          }
 #elif defined(i386)  #elif defined(i386)
         asm volatile("\          asm volatile("\
           pushl   %%ebx;\
         movl    %1,%%esi;\          movl    %1,%%esi;\
         movl    %2,%%edi;\          movl    %2,%%edi;\
         movl    %3,%%ebx;\          movl    %3,%%ebx;\
Line 549  int _subz_main(unsigned int *m1,int d1,unsigned int *m
Line 814  int _subz_main(unsigned int *m1,int d1,unsigned int *m
         jnz Lstart__subz;\          jnz Lstart__subz;\
         movl    $0,%%eax;\          movl    $0,%%eax;\
         adcl    %%eax,%%eax;\          adcl    %%eax,%%eax;\
         movl    %%eax,%0"\          movl    %%eax,%0;\
           popl    %%ebx"\
         :"=m"(br)\          :"=m"(br)\
         :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\          :"m"(m1),"m"(m2),"m"(mr),"m"(d2)\
         :"eax","ebx","ecx","edx","esi","edi");          :"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 ) {
Line 582  int _subz_main(unsigned int *m1,int d1,unsigned int *m
Line 848  int _subz_main(unsigned int *m1,int d1,unsigned int *m
   
 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) ) {
                   u = ZTOS(n);
                   fprintf(asir_out,"%d",u);
           } else {
                 if ( (sd = SL(n)) < 0 ) { SL(n) = -SL(n); fprintf(asir_out,"-"); }                  if ( (sd = SL(n)) < 0 ) { SL(n) = -SL(n); fprintf(asir_out,"-"); }
                 printn((N)n);                  printn((N)n);
                 if ( sd < 0 ) SL(n) = -SL(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.12
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>