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

Diff for /OpenXM_contrib2/asir2000/engine/Mgfs.c between version 1.15 and 1.20

version 1.15, 2002/12/18 06:15:40 version 1.20, 2018/03/29 01:32:51
Line 1 
Line 1 
 /* $OpenXM: OpenXM_contrib2/asir2000/engine/Mgfs.c,v 1.14 2002/09/27 08:40:49 noro Exp $ */  /* $OpenXM: OpenXM_contrib2/asir2000/engine/Mgfs.c,v 1.19 2015/08/14 13:51:54 fujimoto Exp $ */
   
 #include "ca.h"  #include "ca.h"
 #include "inline.h"  #include "inline.h"
Line 8  extern int up_kara_mag, current_gfs_q1;
Line 8  extern int up_kara_mag, current_gfs_q1;
 extern int *current_gfs_plus1;  extern int *current_gfs_plus1;
   
 #if defined(__GNUC__)  #if defined(__GNUC__)
 #define INLINE inline  #define INLINE static inline
 #elif defined(VISUAL)  #elif defined(VISUAL) || defined(__MINGW32__)
 #define INLINE __inline  #define INLINE __inline
 #else  #else
 #define INLINE  #define INLINE
Line 17  extern int *current_gfs_plus1;
Line 17  extern int *current_gfs_plus1;
   
 INLINE int _ADDSF(int a,int b)  INLINE int _ADDSF(int a,int b)
 {  {
         if ( !a )    if ( !a )
                 return b;      return b;
         else if ( !b )    else if ( !b )
                 return a;      return a;
   
         a = IFTOF(a); b = IFTOF(b);    a = IFTOF(a); b = IFTOF(b);
   
         if ( !current_gfs_ntoi ) {    if ( !current_gfs_ntoi ) {
                 a = a+b-current_gfs_q;      a = a+b-current_gfs_q;
                 if ( a == 0 )      if ( a == 0 )
                         return 0;        return 0;
                 else {      else {
                         if ( a < 0 )        if ( a < 0 )
                                 a += current_gfs_q;          a += current_gfs_q;
                         return FTOIF(a);        return FTOIF(a);
                 }      }
         }    }
   
         if ( a > b ) {    if ( a > b ) {
                 /* tab[a]+tab[b] = tab[b](tab[a-b]+1) */      /* tab[a]+tab[b] = tab[b](tab[a-b]+1) */
                 a = current_gfs_plus1[a-b];      a = current_gfs_plus1[a-b];
                 if ( a < 0 )      if ( a < 0 )
                         return 0;        return 0;
                 else {      else {
                         a += b;        a += b;
                         if ( a >= current_gfs_q1 )        if ( a >= current_gfs_q1 )
                                 a -= current_gfs_q1;          a -= current_gfs_q1;
                         return FTOIF(a);        return FTOIF(a);
                 }      }
         } else {    } else {
                 /* tab[a]+tab[b] = tab[a](tab[b-a]+1) */      /* tab[a]+tab[b] = tab[a](tab[b-a]+1) */
                 b = current_gfs_plus1[b-a];      b = current_gfs_plus1[b-a];
                 if ( b < 0 )      if ( b < 0 )
                         return 0;        return 0;
                 else {      else {
                         b += a;        b += a;
                         if ( b >= current_gfs_q1 )        if ( b >= current_gfs_q1 )
                                 b -= current_gfs_q1;          b -= current_gfs_q1;
                         return FTOIF(b);        return FTOIF(b);
                 }      }
         }    }
 }  }
   
 INLINE int _MULSF(int a,int b)  INLINE int _MULSF(int a,int b)
 {  {
         int c;    int c;
   
         if ( !a || !b )    if ( !a || !b )
                 return 0;      return 0;
         else if ( !current_gfs_ntoi ) {    else if ( !current_gfs_ntoi ) {
                 a = IFTOF(a); b = IFTOF(b);      a = IFTOF(a); b = IFTOF(b);
                 DMAR(a,b,0,current_gfs_q,c);      DMAR(a,b,0,current_gfs_q,c);
                 return FTOIF(c);      return FTOIF(c);
         } else {    } else {
                 a = IFTOF(a) + IFTOF(b);      a = IFTOF(a) + IFTOF(b);
                 if ( a >= current_gfs_q1 )      if ( a >= current_gfs_q1 )
                         a -= current_gfs_q1;        a -= current_gfs_q1;
                 return FTOIF(a);      return FTOIF(a);
         }    }
 }  }
   
 void addsfum(UM p1,UM p2,UM pr)  void addsfum(UM p1,UM p2,UM pr)
 {  {
         int *c1,*c2,*cr,i,dmax,dmin;    int *c1,*c2,*cr,i,dmax,dmin;
   
         if ( DEG(p1) == -1 ) {    if ( DEG(p1) == -1 ) {
                 cpyum(p2,pr);      cpyum(p2,pr);
                 return;      return;
         }    }
         if ( DEG(p2) == -1 ) {    if ( DEG(p2) == -1 ) {
                 cpyum(p1,pr);      cpyum(p1,pr);
                 return;      return;
         }    }
         if ( DEG(p1) >= DEG(p2) ) {    if ( DEG(p1) >= DEG(p2) ) {
                 c1 = COEF(p1); c2 = COEF(p2); dmax = DEG(p1); dmin = DEG(p2);      c1 = COEF(p1); c2 = COEF(p2); dmax = DEG(p1); dmin = DEG(p2);
         } else {    } else {
                 c1 = COEF(p2); c2 = COEF(p1); dmax = DEG(p2); dmin = DEG(p1);      c1 = COEF(p2); c2 = COEF(p1); dmax = DEG(p2); dmin = DEG(p1);
         }    }
         for ( i = 0, cr = COEF(pr); i <= dmin; i++ )    for ( i = 0, cr = COEF(pr); i <= dmin; i++ )
                 cr[i] = _ADDSF(c1[i],c2[i]);      cr[i] = _ADDSF(c1[i],c2[i]);
         for ( ; i <= dmax; i++ )    for ( ; i <= dmax; i++ )
                 cr[i] = c1[i];      cr[i] = c1[i];
         if ( dmax == dmin )    if ( dmax == dmin )
                 degum(pr,dmax);      degum(pr,dmax);
         else    else
                 DEG(pr) = dmax;      DEG(pr) = dmax;
 }  }
   
 void subsfum(UM p1,UM p2,UM pr)  void subsfum(UM p1,UM p2,UM pr)
 {  {
         int *c1,*c2,*cr,i;    int *c1,*c2,*cr,i;
         int dmax,dmin;    int dmax,dmin;
   
         if ( DEG(p1) == -1 ) {    if ( DEG(p1) == -1 ) {
                 for ( i = DEG(pr) = DEG(p2), c2 = COEF(p2), cr = COEF(pr);      for ( i = DEG(pr) = DEG(p2), c2 = COEF(p2), cr = COEF(pr);
                           i >= 0; i-- )          i >= 0; i-- )
                         cr[i] = _chsgnsf(c2[i]);        cr[i] = _chsgnsf(c2[i]);
                 return;      return;
         }    }
         if ( DEG(p2) == -1 ) {    if ( DEG(p2) == -1 ) {
                 cpyum(p1,pr);      cpyum(p1,pr);
                 return;      return;
         }    }
         c1 = COEF(p1); c2 = COEF(p2); cr = COEF(pr);    c1 = COEF(p1); c2 = COEF(p2); cr = COEF(pr);
         if ( DEG(p1) >= DEG(p2) ) {    if ( DEG(p1) >= DEG(p2) ) {
                 dmax = DEG(p1); dmin = DEG(p2);      dmax = DEG(p1); dmin = DEG(p2);
                 for ( i = 0; i <= dmin; i++ )      for ( i = 0; i <= dmin; i++ )
                         cr[i] = _subsf(c1[i],c2[i]);        cr[i] = _subsf(c1[i],c2[i]);
                 for ( ; i <= dmax; i++ )      for ( ; i <= dmax; i++ )
                         cr[i] = c1[i];        cr[i] = c1[i];
         } else {    } else {
                 dmax = DEG(p2); dmin = DEG(p1);      dmax = DEG(p2); dmin = DEG(p1);
                 for ( i = 0; i <= dmin; i++ )      for ( i = 0; i <= dmin; i++ )
                         cr[i] = _subsf(c1[i],c2[i]);        cr[i] = _subsf(c1[i],c2[i]);
                 for ( ; i <= dmax; i++ )      for ( ; i <= dmax; i++ )
                         cr[i] = _chsgnsf(c2[i]);        cr[i] = _chsgnsf(c2[i]);
         }    }
         if ( dmax == dmin )    if ( dmax == dmin )
                 degum(pr,dmax);      degum(pr,dmax);
         else    else
                 DEG(pr) = dmax;      DEG(pr) = dmax;
 }  }
   
 void gcdsfum(UM p1,UM p2,UM pr)  void gcdsfum(UM p1,UM p2,UM pr)
 {  {
         int inv;    int inv;
         UM t1,t2,q,tum;    UM t1,t2,q,tum;
         int drem;    int drem;
   
         if ( DEG(p1) < 0 )    if ( DEG(p1) < 0 )
                 cpyum(p2,pr);      cpyum(p2,pr);
         else if ( DEG(p2) < 0 )    else if ( DEG(p2) < 0 )
                 cpyum(p1,pr);      cpyum(p1,pr);
         else {    else {
                 if ( DEG(p1) >= DEG(p2) ) {      if ( DEG(p1) >= DEG(p2) ) {
                         t1 = p1; t2 = p2;        t1 = p1; t2 = p2;
                 } else {      } else {
                         t1 = p2; t2 = p1;        t1 = p2; t2 = p1;
                 }      }
                 q = W_UMALLOC(DEG(t1));      q = W_UMALLOC(DEG(t1));
                 while ( ( drem = divsfum(t1,t2,q) ) >= 0 ) {      while ( ( drem = divsfum(t1,t2,q) ) >= 0 ) {
                         tum = t1; t1 = t2; t2 = tum; DEG(t2) = drem;        tum = t1; t1 = t2; t2 = tum; DEG(t2) = drem;
                 }      }
                 inv = _invsf(COEF(t2)[DEG(t2)]);      inv = _invsf(COEF(t2)[DEG(t2)]);
                 mulssfum(t2,inv,pr);      mulssfum(t2,inv,pr);
         }    }
 }  }
   
 void mulsfum(UM p1,UM p2,UM pr)  void mulsfum(UM p1,UM p2,UM pr)
 {  {
         int *pc1,*pcr;    int *pc1,*pcr;
         int *c1,*c2,*cr;    int *c1,*c2,*cr;
         int mul;    int mul;
         int i,j,d1,d2;    int i,j,d1,d2;
   
         if ( ( (d1 = DEG(p1)) < 0) || ( (d2 = DEG(p2)) < 0 ) ) {    if ( ( (d1 = DEG(p1)) < 0) || ( (d2 = DEG(p2)) < 0 ) ) {
                 DEG(pr) = -1;      DEG(pr) = -1;
                 return;      return;
         }    }
         c1 = COEF(p1); c2 = COEF(p2); cr = COEF(pr);    c1 = COEF(p1); c2 = COEF(p2); cr = COEF(pr);
         bzero((char *)cr,(int)((d1+d2+1)*sizeof(int)));    bzero((char *)cr,(int)((d1+d2+1)*sizeof(int)));
         for ( i = 0; i <= d2; i++, cr++ )    for ( i = 0; i <= d2; i++, cr++ )
                 if ( mul = *c2++ )      if ( mul = *c2++ )
                         for ( j = 0, pc1 = c1, pcr = cr; j <= d1; j++, pc1++, pcr++ )        for ( j = 0, pc1 = c1, pcr = cr; j <= d1; j++, pc1++, pcr++ )
                                 if ( *pc1 )          if ( *pc1 )
                                         *pcr = _ADDSF(_MULSF(*pc1,mul),*pcr);            *pcr = _ADDSF(_MULSF(*pc1,mul),*pcr);
         DEG(pr) = d1 + d2;    DEG(pr) = d1 + d2;
 }  }
   
 void mulssfum(UM p,int n,UM pr)  void mulssfum(UM p,int n,UM pr)
 {  {
         int *sp,*dp;    int *sp,*dp;
         int i;    int i;
   
         for ( i = DEG(pr) = DEG(p), sp = COEF(p)+i, dp = COEF(pr)+i;    for ( i = DEG(pr) = DEG(p), sp = COEF(p)+i, dp = COEF(pr)+i;
                   i >= 0; i--, dp--, sp-- )        i >= 0; i--, dp--, sp-- )
                 *dp = _MULSF(*sp,n);      *dp = _MULSF(*sp,n);
 }  }
   
 void kmulsfum(UM n1,UM n2,UM nr)  void kmulsfum(UM n1,UM n2,UM nr)
 {  {
         UM n,t,s,m,carry;    UM n,t,s,m,carry;
         int d,d1,d2,len,i,l;    int d,d1,d2,len,i,l;
         unsigned int *r,*r0;    unsigned int *r,*r0;
   
         if ( !n1 || !n2 ) {    if ( !n1 || !n2 ) {
                 nr->d = -1; return;      nr->d = -1; return;
         }    }
         d1 = DEG(n1)+1; d2 = DEG(n2)+1;    d1 = DEG(n1)+1; d2 = DEG(n2)+1;
         if ( (d1 < up_kara_mag) || (d2 < up_kara_mag) ) {    if ( (d1 < up_kara_mag) || (d2 < up_kara_mag) ) {
                 mulsfum(n1,n2,nr); return;      mulsfum(n1,n2,nr); return;
         }    }
         if ( d1 < d2 ) {    if ( d1 < d2 ) {
                 n = n1; n1 = n2; n2 = n;      n = n1; n1 = n2; n2 = n;
                 d = d1; d1 = d2; d2 = d;      d = d1; d1 = d2; d2 = d;
         }    }
         if ( d2 > (d1+1)/2 ) {    if ( d2 > (d1+1)/2 ) {
                 kmulsfummain(n1,n2,nr); return;      kmulsfummain(n1,n2,nr); return;
         }    }
         d = (d1/d2)+((d1%d2)!=0);    d = (d1/d2)+((d1%d2)!=0);
         len = (d+1)*d2;    len = (d+1)*d2;
         r0 = (unsigned int *)ALLOCA(len*sizeof(int));    r0 = (unsigned int *)ALLOCA(len*sizeof(int));
         bzero((char *)r0,len*sizeof(int));    bzero((char *)r0,len*sizeof(int));
         m = W_UMALLOC(d2+1);    m = W_UMALLOC(d1+d2+1);
         carry = W_UMALLOC(d2+1);    carry = W_UMALLOC(d1+d2+1);
         t = W_UMALLOC(d1+d2+1);    t = W_UMALLOC(d1+d2+1);
         s = W_UMALLOC(d1+d2+1);    s = W_UMALLOC(d1+d2+1);
         for ( DEG(carry) = -1, i = 0, r = r0; i < d; i++, r += d2 ) {    for ( DEG(carry) = -1, i = 0, r = r0; i < d; i++, r += d2 ) {
                 extractum(n1,i*d2,d2,m);      extractum(n1,i*d2,d2,m);
                 if ( m ) {      if ( m ) {
                         kmulsfum(m,n2,t);        kmulsfum(m,n2,t);
                         addsfum(t,carry,s);        addsfum(t,carry,s);
                         c_copyum(s,d2,r);        c_copyum(s,d2,r);
                         extractum(s,d2,d2,carry);        extractum(s,d2,d2,carry);
                 } else {      } else {
                         c_copyum(carry,d2,r);        c_copyum(carry,d2,r);
                         carry = 0;        carry = 0;
                 }      }
         }    }
         c_copyum(carry,d2,r);    c_copyum(carry,d2,r);
         for ( l = len - 1; !r0[l]; l-- );    for ( l = len - 1; !r0[l]; l-- );
         l++;    l++;
         DEG(nr) = l-1;    DEG(nr) = l-1;
         bcopy((char *)r0,(char *)COEF(nr),l*sizeof(int));    bcopy((char *)r0,(char *)COEF(nr),l*sizeof(int));
 }  }
   
 void kmulsfummain(UM n1,UM n2,UM nr)  void kmulsfummain(UM n1,UM n2,UM nr)
 {  {
         int d1,d2,h,len;    int d1,d2,h,len,d;
         UM n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;    UM n1lo,n1hi,n2lo,n2hi,hi,lo,mid1,mid2,mid,s1,s2,t1,t2;
   
         d1 = DEG(n1)+1; d2 = DEG(n2)+1; h = (d1+1)/2;    d1 = DEG(n1)+1; d2 = DEG(n2)+1; h = (d1+1)/2;
         n1lo = W_UMALLOC(d1+1); n1hi = W_UMALLOC(d1+1);    d = d1+d2+1;
         n2lo = W_UMALLOC(d2+1); n2hi = W_UMALLOC(d2+1);    n1lo = W_UMALLOC(d); n1hi = W_UMALLOC(d);
         lo = W_UMALLOC(d1+d2+1); hi = W_UMALLOC(d1+d2+1);    n2lo = W_UMALLOC(d); n2hi = W_UMALLOC(d);
         mid1 = W_UMALLOC(d1+d2+1); mid2 = W_UMALLOC(d1+d2+1);    lo = W_UMALLOC(d); hi = W_UMALLOC(d);
         mid = W_UMALLOC(d1+d2+1);    mid1 = W_UMALLOC(d); mid2 = W_UMALLOC(d);
         s1 = W_UMALLOC(d1+d2+1); s2 = W_UMALLOC(d1+d2+1);    mid = W_UMALLOC(d);
         extractum(n1,0,h,n1lo); extractum(n1,h,d1-h,n1hi);    s1 = W_UMALLOC(d); s2 = W_UMALLOC(d);
         extractum(n2,0,h,n2lo); extractum(n2,h,d2-h,n2hi);    extractum(n1,0,h,n1lo); extractum(n1,h,d1-h,n1hi);
         kmulsfum(n1hi,n2hi,hi); kmulsfum(n1lo,n2lo,lo);    extractum(n2,0,h,n2lo); extractum(n2,h,d2-h,n2hi);
         len = DEG(hi)+1+2*h; t1 = W_UMALLOC(len-1); DEG(t1) = len-1;    kmulsfum(n1hi,n2hi,hi); kmulsfum(n1lo,n2lo,lo);
         bzero((char *)COEF(t1),len*sizeof(int));    len = DEG(hi)+1+2*h; t1 = W_UMALLOC(len-1); DEG(t1) = len-1;
         if ( lo )    bzero((char *)COEF(t1),len*sizeof(int));
                 bcopy((char *)COEF(lo),(char *)COEF(t1),(DEG(lo)+1)*sizeof(int));    if ( lo )
         if ( hi )      bcopy((char *)COEF(lo),(char *)COEF(t1),(DEG(lo)+1)*sizeof(int));
                 bcopy((char *)COEF(hi),(char *)(COEF(t1)+2*h),(DEG(hi)+1)*sizeof(int));    if ( hi )
       bcopy((char *)COEF(hi),(char *)(COEF(t1)+2*h),(DEG(hi)+1)*sizeof(int));
   
         addsfum(hi,lo,mid1);    addsfum(hi,lo,mid1);
         subsfum(n1hi,n1lo,s1); subsfum(n2lo,n2hi,s2);    subsfum(n1hi,n1lo,s1); subsfum(n2lo,n2hi,s2);
         kmulsfum(s1,s2,mid2); addsfum(mid1,mid2,mid);    kmulsfum(s1,s2,mid2); addsfum(mid1,mid2,mid);
         if ( mid ) {    if ( mid ) {
                 len = DEG(mid)+1+h; t2 = W_UMALLOC(len-1); DEG(t2) = len-1;      len = DEG(mid)+1+h; t2 = W_UMALLOC(len-1); DEG(t2) = len-1;
                 bzero((char *)COEF(t2),len*sizeof(int));      bzero((char *)COEF(t2),len*sizeof(int));
                 bcopy((char *)COEF(mid),(char *)(COEF(t2)+h),(DEG(mid)+1)*sizeof(int));      bcopy((char *)COEF(mid),(char *)(COEF(t2)+h),(DEG(mid)+1)*sizeof(int));
                 addsfum(t1,t2,nr);      addsfum(t1,t2,nr);
         } else    } else
                 copyum(t1,nr);      copyum(t1,nr);
 }  }
   
 int divsfum(UM p1,UM p2,UM pq)  int divsfum(UM p1,UM p2,UM pq)
 {  {
         int *pc1,*pct;    int *pc1,*pct;
         int *c1,*c2,*ct;    int *c1,*c2,*ct;
         int inv,hd,tmp;    int inv,hd,tmp;
         int i,j, d1,d2,dd;    int i,j, d1,d2,dd;
   
         if ( (d1 = DEG(p1)) < (d2 = DEG(p2)) ) {    if ( (d1 = DEG(p1)) < (d2 = DEG(p2)) ) {
                 DEG(pq) = -1;      DEG(pq) = -1;
                 return d1;      return d1;
         }    }
         c1 = COEF(p1); c2 = COEF(p2); dd = d1-d2;    c1 = COEF(p1); c2 = COEF(p2); dd = d1-d2;
         if ( ( hd = c2[d2] ) != _onesf() ) {    if ( ( hd = c2[d2] ) != _onesf() ) {
                 inv = _invsf(hd);      inv = _invsf(hd);
                 for ( pc1 = c2 + d2; pc1 >= c2; pc1-- )      for ( pc1 = c2 + d2; pc1 >= c2; pc1-- )
                         *pc1 = _MULSF(*pc1,inv);        *pc1 = _MULSF(*pc1,inv);
         } else    } else
                 inv = _onesf();      inv = _onesf();
         for ( i = dd, ct = c1+d1; i >= 0; i-- )    for ( i = dd, ct = c1+d1; i >= 0; i-- )
                 if ( tmp = *ct-- ) {      if ( tmp = *ct-- ) {
                         tmp = _chsgnsf(tmp);        tmp = _chsgnsf(tmp);
                         for ( j = d2-1, pct = ct, pc1 = c2+j; j >= 0; j--, pct--, pc1-- )        for ( j = d2-1, pct = ct, pc1 = c2+j; j >= 0; j--, pct--, pc1-- )
                                 *pct = _ADDSF(_MULSF(*pc1,tmp),*pct);          *pct = _ADDSF(_MULSF(*pc1,tmp),*pct);
                 }      }
         if ( inv != _onesf() ) {    if ( inv != _onesf() ) {
                 for ( pc1 = c1+d2, pct = c1+d1; pc1 <= pct; pc1++ )      for ( pc1 = c1+d2, pct = c1+d1; pc1 <= pct; pc1++ )
                         *pc1 = _MULSF(*pc1,inv);        *pc1 = _MULSF(*pc1,inv);
                 for ( pc1 = c2, pct = c2+d2; pc1 <= pct; pc1++ )      for ( pc1 = c2, pct = c2+d2; pc1 <= pct; pc1++ )
                         *pc1 = _MULSF(*pc1,hd);        *pc1 = _MULSF(*pc1,hd);
         }    }
         for ( i = d2-1, pc1 = c1+i; i >= 0 && !(*pc1); pc1--, i-- );    for ( i = d2-1, pc1 = c1+i; i >= 0 && !(*pc1); pc1--, i-- );
         for ( DEG(pq) = j = dd, pc1 = c1+d1, pct = COEF(pq)+j; j >= 0; j-- )    for ( DEG(pq) = j = dd, pc1 = c1+d1, pct = COEF(pq)+j; j >= 0; j-- )
                 *pct-- = *pc1--;      *pct-- = *pc1--;
         return i;    return i;
 }  }
   
 void diffsfum(UM f,UM fd)  void diffsfum(UM f,UM fd)
 {  {
         int *dp,*sp;    int *dp,*sp;
         int i;    int i;
   
         for ( i = DEG(f), dp = COEF(fd)+i-1, sp = COEF(f)+i;    for ( i = DEG(f), dp = COEF(fd)+i-1, sp = COEF(f)+i;
                 i >= 1; i--, dp--, sp-- ) {      i >= 1; i--, dp--, sp-- ) {
                 *dp = _MULSF(*sp,_itosf(i%current_gfs_p));      *dp = _MULSF(*sp,_itosf(i%current_gfs_p));
         }    }
         degum(fd,DEG(f) - 1);    degum(fd,DEG(f) - 1);
 }  }
   
 void monicsfum(UM f)  void monicsfum(UM f)
 {  {
         int *sp;    int *sp;
         int i,inv;    int i,inv;
   
         i = DEG(f); sp = COEF(f)+i;    i = DEG(f); sp = COEF(f)+i;
         inv = _invsf(*sp);    inv = _invsf(*sp);
         for ( ; i >= 0; i--, sp-- )    for ( ; i >= 0; i--, sp-- )
                 *sp = _MULSF(*sp,inv);      *sp = _MULSF(*sp,inv);
 }  }
   
 int compsfum(UM a,UM b)  int compsfum(UM a,UM b)
 {  {
         int i,da,db;    int i,da,db;
   
         if ( !a )    if ( !a )
                 return !b?0:1;      return !b?0:1;
         else if ( !b )    else if ( !b )
                 return 1;      return 1;
         else if ( (da = DEG(a)) > (db = DEG(b)) )    else if ( (da = DEG(a)) > (db = DEG(b)) )
                 return 1;      return 1;
         else if ( da < db )    else if ( da < db )
                 return -1;      return -1;
         else {    else {
                 for ( i = da; i >= 0 && COEF(a)[i] == COEF(b)[i]; i-- );      for ( i = da; i >= 0 && COEF(a)[i] == COEF(b)[i]; i-- );
                 if ( i < 0 )      if ( i < 0 )
                         return 0;        return 0;
                 else if ( (unsigned int)COEF(a)[i] > (unsigned int)COEF(b)[i] )      else if ( (unsigned int)COEF(a)[i] > (unsigned int)COEF(b)[i] )
                         return 1;        return 1;
                 else      else
                         return -1;        return -1;
         }    }
 }  }
   
 void addsfarray(int,int *,int *);  void addsfarray(int,int *,int *);
Line 366  void mulsfarray_trunc(int,int *,int *,int *);
Line 367  void mulsfarray_trunc(int,int *,int *,int *);
   
 void mulsfbm(BM f1,BM f2,BM fr)  void mulsfbm(BM f1,BM f2,BM fr)
 {  {
         UM mul,t,s;    UM mul,t,s;
         int i,j,d1,d2,dy;    int i,j,d1,d2,dy;
   
         dy = MIN(DEG(f1),DEG(f2));    dy = MIN(DEG(f1),DEG(f2));
   
         d1 = degbm(f1);    d1 = degbm(f1);
         d2 = degbm(f2);    d2 = degbm(f2);
         t = W_UMALLOC(d1+d2);    t = W_UMALLOC(d1+d2);
         s = W_UMALLOC(d1+d2);    s = W_UMALLOC(d1+d2);
   
         for ( i = 0; i <= dy; i++ ) {    for ( i = 0; i <= dy; i++ ) {
                 mul = COEF(f2)[i];      mul = COEF(f2)[i];
                 if ( DEG(mul) >= 0 )      if ( DEG(mul) >= 0 )
                 for ( j = 0; i+j <= dy; j++ ) {      for ( j = 0; i+j <= dy; j++ ) {
                         if ( COEF(f1)[j] ) {        if ( COEF(f1)[j] ) {
                                 kmulsfum(COEF(f1)[j],mul,t);          kmulsfum(COEF(f1)[j],mul,t);
                                 addsfum(t,COEF(fr)[i+j],s);          addsfum(t,COEF(fr)[i+j],s);
                                 cpyum(s,COEF(fr)[i+j]);          cpyum(s,COEF(fr)[i+j]);
                         }        }
                 }      }
         }    }
         DEG(fr) = dy;    DEG(fr) = dy;
 }  }
   
 int degbm(BM f)  int degbm(BM f)
 {  {
         int d,i,dy;    int d,i,dy;
   
         dy = DEG(f);    dy = DEG(f);
         d = DEG(COEF(f)[0]);    d = DEG(COEF(f)[0]);
         for ( i = 1; i <= dy; i++ )    for ( i = 1; i <= dy; i++ )
                 d = MAX(DEG(COEF(f)[i]),d);      d = MAX(DEG(COEF(f)[i]),d);
         return d;    return d;
 }  }
   
 /* g += f */  /* g += f */
   
 void addtosfbm(BM f,BM g)  void addtosfbm(BM f,BM g)
 {  {
         int i,d1,d2,dy;    int i,d1,d2,dy;
         UM t;    UM t;
   
         dy = DEG(g);    dy = DEG(g);
         if ( DEG(f) > dy )    if ( DEG(f) > dy )
                 error("addtosfbm : invalid input");      error("addtosfbm : invalid input");
         dy = MIN(DEG(f),dy);    dy = MIN(DEG(f),dy);
         d1 = degbm(f);    d1 = degbm(f);
         d2 = degbm(g);    d2 = degbm(g);
         t = W_UMALLOC(MAX(d1,d2));    t = W_UMALLOC(MAX(d1,d2));
         for ( i = 0; i <= dy; i++ ) {    for ( i = 0; i <= dy; i++ ) {
                 addsfum(COEF(f)[i],COEF(g)[i],t);      addsfum(COEF(f)[i],COEF(g)[i],t);
                 cpyum(t,COEF(g)[i]);      cpyum(t,COEF(g)[i]);
         }    }
 }  }
   
 void eucsfum(UM f1,UM f2,UM a,UM b)  void eucsfum(UM f1,UM f2,UM a,UM b)
 {  {
         UM g1,g2,a1,a2,a3,wm,q,tum;    UM g1,g2,a1,a2,a3,wm,q,tum;
         int d,dr;    int d,dr;
   
         /* XXX */    /* XXX */
         d = DEG(f1) + DEG(f2) + 10;    d = DEG(f1) + DEG(f2) + 10;
         g1 = W_UMALLOC(d); g2 = W_UMALLOC(d); a1 = W_UMALLOC(d);    g1 = W_UMALLOC(d); g2 = W_UMALLOC(d); a1 = W_UMALLOC(d);
         a2 = W_UMALLOC(d); a3 = W_UMALLOC(d); wm = W_UMALLOC(d);    a2 = W_UMALLOC(d); a3 = W_UMALLOC(d); wm = W_UMALLOC(d);
         q = W_UMALLOC(d);    q = W_UMALLOC(d);
         DEG(a1) = 0; COEF(a1)[0] = _onesf(); DEG(a2) = -1;    DEG(a1) = 0; COEF(a1)[0] = _onesf(); DEG(a2) = -1;
         cpyum(f1,g1); cpyum(f2,g2);    cpyum(f1,g1); cpyum(f2,g2);
         while ( 1 ) {    while ( 1 ) {
                 dr = divsfum(g1,g2,q); tum = g1; g1 = g2; g2 = tum;      dr = divsfum(g1,g2,q); tum = g1; g1 = g2; g2 = tum;
                 if ( ( DEG(g2) = dr ) == -1 )      if ( ( DEG(g2) = dr ) == -1 )
                         break;        break;
                 mulsfum(a2,q,wm); subsfum(a1,wm,a3); dr = divsfum(a3,f2,q);      mulsfum(a2,q,wm); subsfum(a1,wm,a3); dr = divsfum(a3,f2,q);
                 tum = a1; a1 = a2; a2 = a3; a3 = tum; DEG(a3) = dr;      tum = a1; a1 = a2; a2 = a3; a3 = tum; DEG(a3) = dr;
         }    }
         if ( COEF(g1)[0] != _onesf() )    if ( COEF(g1)[0] != _onesf() )
                 mulssfum(a2,_invsf(COEF(g1)[0]),a);      mulssfum(a2,_invsf(COEF(g1)[0]),a);
         else    else
                 cpyum(a2,a);      cpyum(a2,a);
         mulsfum(a,f1,wm);    mulsfum(a,f1,wm);
         if ( DEG(wm) >= 0 )    if ( DEG(wm) >= 0 )
                 COEF(wm)[0] = _subsf(COEF(wm)[0],_onesf());      COEF(wm)[0] = _subsf(COEF(wm)[0],_onesf());
         else {    else {
                 DEG(wm) = 0; COEF(wm)[0] = _chsgnsf(_onesf());      DEG(wm) = 0; COEF(wm)[0] = _chsgnsf(_onesf());
         }    }
         divsfum(wm,f2,q); mulssfum(q,_chsgnsf(_onesf()),b);    divsfum(wm,f2,q); mulssfum(q,_chsgnsf(_onesf()),b);
 }  }
   
 void shiftsfum(UM,int,UM);  void shiftsfum(UM,int,UM);
   
 void shiftsflum(int n,LUM f,int ev)  void shiftsflum(int n,LUM f,int ev)
 {  {
         int d,i,j;    int d,i,j;
         UM w,w1;    UM w,w1;
         int *coef;    int *coef;
         int **c;    int **c;
   
         d = f->d;    d = f->d;
         c = f->c;    c = f->c;
         w = W_UMALLOC(n);    w = W_UMALLOC(n);
         w1 = W_UMALLOC(n);    w1 = W_UMALLOC(n);
         for ( i = 0; i <= d; i++ ) {    for ( i = 0; i <= d; i++ ) {
                 coef = c[i];      coef = c[i];
                 for ( j = n-1; j >= 0 && coef[j] == 0; j-- );      for ( j = n-1; j >= 0 && coef[j] == 0; j-- );
                 if ( j <= 0 )      if ( j <= 0 )
                         continue;        continue;
                 DEG(w) = j; bcopy(coef,COEF(w),(j+1)*sizeof(int));      DEG(w) = j; bcopy(coef,COEF(w),(j+1)*sizeof(int));
                 shiftsfum(w,ev,w1);      shiftsfum(w,ev,w1);
                 bcopy(COEF(w1),coef,(j+1)*sizeof(int));      bcopy(COEF(w1),coef,(j+1)*sizeof(int));
         }    }
 }  }
   
 /* f(x) -> g(x) = f(x+a) */  /* f(x) -> g(x) = f(x+a) */
   
 void shiftsfum(UM f,int a,UM g)  void shiftsfum(UM f,int a,UM g)
 {  {
         int n,i;    int n,i;
         UM pwr,xa,w,t;    UM pwr,xa,w,t;
         int *c;    int *c;
   
         n = DEG(f);    n = DEG(f);
         if ( n <= 0 )    if ( n <= 0 )
                 cpyum(f,g);      cpyum(f,g);
         else {    else {
                 c = COEF(f);      c = COEF(f);
   
                 /* g = 0 */      /* g = 0 */
                 DEG(g) = -1;      DEG(g) = -1;
   
                 /* pwr = 1 */      /* pwr = 1 */
                 pwr = W_UMALLOC(n); DEG(pwr) = 0; COEF(pwr)[0] = _onesf();      pwr = W_UMALLOC(n); DEG(pwr) = 0; COEF(pwr)[0] = _onesf();
   
                 /* xa = x+a */      /* xa = x+a */
                 xa = W_UMALLOC(n); DEG(xa) = 1;      xa = W_UMALLOC(n); DEG(xa) = 1;
                 COEF(xa)[0] = a; COEF(xa)[1] = _onesf();      COEF(xa)[0] = a; COEF(xa)[1] = _onesf();
   
                 w = W_UMALLOC(n);      w = W_UMALLOC(n);
                 t = W_UMALLOC(n);      t = W_UMALLOC(n);
   
                 /* g += pwr*c[0] */      /* g += pwr*c[0] */
                 for ( i = 0; i <= n; i++ ) {      for ( i = 0; i <= n; i++ ) {
                         mulssfum(pwr,c[i],w); addsfum(g,w,t); cpyum(t,g);        mulssfum(pwr,c[i],w); addsfum(g,w,t); cpyum(t,g);
                         if ( i < n ) {        if ( i < n ) {
                                 mulsfum(pwr,xa,t); cpyum(t,pwr);          mulsfum(pwr,xa,t); cpyum(t,pwr);
                         }        }
                 }      }
         }    }
 }  }
   
 /* f(y) -> f(y+a) */  /* f(y) -> f(y+a) */
   
 void shiftsfbm(BM f,int a)  void shiftsfbm(BM f,int a)
 {  {
         int i,j,d,dy;    int i,j,d,dy;
         UM pwr,ya,w,t,s;    UM pwr,ya,w,t,s;
         UM *c;    UM *c;
   
         dy = DEG(f);    dy = DEG(f);
         c = COEF(f);    c = COEF(f);
         d = degbm(f);    d = degbm(f);
   
         w = W_UMALLOC(d);    w = W_UMALLOC(d);
         t = W_UMALLOC(d);    t = W_UMALLOC(d);
         s = W_UMALLOC(dy);    s = W_UMALLOC(dy);
   
         /* pwr = 1 */    /* pwr = 1 */
         pwr = W_UMALLOC(dy); DEG(pwr) = 0; COEF(pwr)[0] = _onesf();    pwr = W_UMALLOC(dy); DEG(pwr) = 0; COEF(pwr)[0] = _onesf();
   
         /* ya = y+a */    /* ya = y+a */
         ya = W_UMALLOC(1); DEG(ya) = 1;    ya = W_UMALLOC(1); DEG(ya) = 1;
         COEF(ya)[0] = a; COEF(ya)[1] = _onesf();    COEF(ya)[0] = a; COEF(ya)[1] = _onesf();
   
         for ( i = 0; i <= dy; i++ ) {    for ( i = 0; i <= dy; i++ ) {
                 /* c[i] does not change */      /* c[i] does not change */
                 for ( j = 0; j < i; j++ ) {      for ( j = 0; j < i; j++ ) {
                         mulssfum(c[i],COEF(pwr)[j],w);        mulssfum(c[i],COEF(pwr)[j],w);
                         addsfum(w,c[j],t); cpyum(t,c[j]);        addsfum(w,c[j],t); cpyum(t,c[j]);
                 }      }
                 if ( i < dy ) {      if ( i < dy ) {
                         mulsfum(pwr,ya,s); cpyum(s,pwr);        mulsfum(pwr,ya,s); cpyum(s,pwr);
                 }      }
         }    }
 }  }
   
 void clearbm(int n,BM f)  void clearbm(int n,BM f)
 {  {
         int i,dy;    int i,dy;
         UM *c;    UM *c;
   
         dy = DEG(f);    dy = DEG(f);
         for ( c = COEF(f), i = 0; i <= dy; i++ )    for ( c = COEF(f), i = 0; i <= dy; i++ )
                 clearum(c[i],n);      clearum(c[i],n);
 }  }
   
 static void create_bmono(P c,V x,int i,V y,int j,P *mono);  static void create_bmono(P c,V x,int i,V y,int j,P *mono);
   
 struct lb {  struct lb {
         int pos,len;    int pos,len;
         int *r;    int *r;
         int *hist;    int *hist;
 };  };
   
 static NODE insert_lb(NODE g,struct lb *a)  static NODE insert_lb(NODE g,struct lb *a)
 {  {
         NODE prev,cur,n;    NODE prev,cur,n;
   
         prev = 0; cur = g;    prev = 0; cur = g;
         while ( cur ) {    while ( cur ) {
                 if ( a->pos < ((struct lb *)BDY(cur))->pos ) {      if ( a->pos < ((struct lb *)BDY(cur))->pos ) {
                         MKNODE(n,a,cur);        MKNODE(n,a,cur);
                         if ( !prev )        if ( !prev )
                                 return n;          return n;
                         else {        else {
                                 NEXT(prev) = n;          NEXT(prev) = n;
                                 return g;          return g;
                         }        }
                 } else {      } else {
                         prev = cur;        prev = cur;
                         cur = NEXT(cur);        cur = NEXT(cur);
                 }      }
         }    }
         MKNODE(n,a,0);    MKNODE(n,a,0);
         NEXT(prev) = n;    NEXT(prev) = n;
         return g;    return g;
 }  }
   
 static void lnf(int *r,int *h,int n,int len,NODE g)  static void lnf(int *r,int *h,int n,int len,NODE g)
 {  {
         struct lb *t;    struct lb *t;
         int pos,i,j,len1,c;    int pos,i,j,len1,c;
         int *r1,*h1;    int *r1,*h1;
   
         for ( ; g; g = NEXT(g) ) {    for ( ; g; g = NEXT(g) ) {
                 t = (struct lb *)BDY(g);      t = (struct lb *)BDY(g);
                 pos = t->pos;      pos = t->pos;
                 if ( c = r[pos] ) {      if ( c = r[pos] ) {
                         c = _chsgnsf(c);        c = _chsgnsf(c);
                         r1 = t->r;        r1 = t->r;
                         h1 = t->hist;        h1 = t->hist;
                         len1 = t->len;        len1 = t->len;
                         for ( i = pos; i < n; i++ )        for ( i = pos; i < n; i++ )
                                 if ( r1[i] )          if ( r1[i] )
                                         r[i] = _ADDSF(r[i],_MULSF(r1[i],c));            r[i] = _ADDSF(r[i],_MULSF(r1[i],c));
                         for ( i = 0; i < len1; i++ )        for ( i = 0; i < len1; i++ )
                                 if ( h1[i] )          if ( h1[i] )
                                         h[i] = _ADDSF(h[i],_MULSF(h1[i],c));            h[i] = _ADDSF(h[i],_MULSF(h1[i],c));
                 }      }
         }    }
         for ( i = 0; i < n && !r[i]; i++ );    for ( i = 0; i < n && !r[i]; i++ );
         if ( i < n ) {    if ( i < n ) {
                 c = _invsf(r[i]);      c = _invsf(r[i]);
                 for ( j = i; j < n; j++ )      for ( j = i; j < n; j++ )
                         if ( r[j] )        if ( r[j] )
                                 r[j] = _MULSF(r[j],c);          r[j] = _MULSF(r[j],c);
                 for ( j = 0; j < len; j++ )      for ( j = 0; j < len; j++ )
                         if ( h[j] )        if ( h[j] )
                                 h[j] = _MULSF(h[j],c);          h[j] = _MULSF(h[j],c);
         }    }
 }  }
   
 void sfmintdeg(VL vl,P fx,int dy,int c,P *fr)  void sfmintdeg(VL vl,P fx,int dy,int c,P *fr)
 {  {
         V x,y;    V x,y;
         int dx,dxdy,i,j,k,l,d,len,len0,u,dyk;    int dx,dxdy,i,j,k,l,d,len,len0,u,dyk;
         UP *rx;    UP *rx;
         DCP dc;    DCP dc;
         P t,f,mono,f1;    P t,f,mono,f1;
         UP ut,h;    UP ut,h;
         int ***nf;    int ***nf;
         int *r,*hist,*prev,*r1;    int *r,*hist,*prev,*r1;
         struct lb *lb;    struct lb *lb;
         GFS s;    GFS s;
         NODE g;    NODE g;
   
         x = vl->v;    x = vl->v;
         y = NEXT(vl)->v;    y = NEXT(vl)->v;
         dx = getdeg(x,fx);    dx = getdeg(x,fx);
         dxdy = dx*dy;    dxdy = dx*dy;
         /* rx = -(fx-x^dx) */    /* rx = -(fx-x^dx) */
         rx = (UP *)CALLOC(dx,sizeof(UP));    rx = (UP *)CALLOC(dx,sizeof(UP));
         for ( dc = DC(fx); dc; dc = NEXT(dc)) {    for ( dc = DC(fx); dc; dc = NEXT(dc)) {
                 chsgnp(COEF(dc),&t);      chsgnp(COEF(dc),&t);
                 ptoup(t,&ut);      ptoup(t,&ut);
                 rx[QTOS(DEG(dc))] = ut;      rx[QTOS(DEG(dc))] = ut;
         }    }
         /* nf[d] = normal form table of monomials with total degree d */    /* nf[d] = normal form table of monomials with total degree d */
         nf = (int ***)CALLOC(dx+dy+1,sizeof(int **)); /* xxx */    nf = (int ***)CALLOC(dx+dy+1,sizeof(int **)); /* xxx */
         nf[0] = (int **)CALLOC(1,sizeof(int *));    nf[0] = (int **)CALLOC(1,sizeof(int *));
   
         /* nf[0][0] = 1 */    /* nf[0][0] = 1 */
         r = (int *)CALLOC(dxdy,sizeof(int));    r = (int *)CALLOC(dxdy,sizeof(int));
         r[0] = _onesf();    r[0] = _onesf();
         nf[0][0] = r;    nf[0][0] = r;
   
         hist = (int *)CALLOC(1,sizeof(int));    hist = (int *)CALLOC(1,sizeof(int));
         hist[0] = _onesf();    hist[0] = _onesf();
   
         lb = (struct lb *)CALLOC(1,sizeof(struct lb));    lb = (struct lb *)CALLOC(1,sizeof(struct lb));
         lb->pos = 0;    lb->pos = 0;
         lb->r = r;    lb->r = r;
         lb->hist = hist;    lb->hist = hist;
         lb->len = 1;    lb->len = 1;
   
         /* g : table of normal form as linear form */    /* g : table of normal form as linear form */
         MKNODE(g,lb,0);    MKNODE(g,lb,0);
   
         len = 1;    len = 1;
         h = UPALLOC(dy);    h = UPALLOC(dy);
         for ( d = 1; ; d++ ) {    for ( d = 1; ; d++ ) {
                 if ( d > c ){      if ( d > c ){
                         return;        return;
                 }      }
                 nf[d] = (int **)CALLOC(d+1,sizeof(int *));      nf[d] = (int **)CALLOC(d+1,sizeof(int *));
                 len0 = len;      len0 = len;
                 len += d+1;      len += d+1;
   
                 for ( i = d; i >= 0; i-- ) {      for ( i = d; i >= 0; i-- ) {
                         /* nf(x^(d-i)*y^i) = nf(y*nf(x^(d-i)*y^(i-1))) */        /* nf(x^(d-i)*y^i) = nf(y*nf(x^(d-i)*y^(i-1))) */
                         /* nf(x^d) = nf(nf(x^(d-1))*x) */        /* nf(x^d) = nf(nf(x^(d-1))*x) */
                         r = (int *)CALLOC(dxdy,sizeof(int));        r = (int *)CALLOC(dxdy,sizeof(int));
                         if ( i == 0 ) {        if ( i == 0 ) {
                                 prev = nf[d-1][0];          prev = nf[d-1][0];
                                 bcopy(prev,r+dy,(dxdy-dy)*sizeof(int));          bcopy(prev,r+dy,(dxdy-dy)*sizeof(int));
   
                                 /* create the head coeff */          /* create the head coeff */
                                 for ( l = 0, k = dxdy-dy; l < dy; l++, k++ ) {          for ( l = 0, k = dxdy-dy; l < dy; l++, k++ ) {
                                         iftogfs(prev[k],&s);            iftogfs(prev[k],&s);
                                         COEF(h)[l] = (Num)s;            COEF(h)[l] = (Num)s;
                                 }          }
                                 for ( l = dy-1; l >= 0 && !COEF(h)[l]; l--);          for ( l = dy-1; l >= 0 && !COEF(h)[l]; l--);
                                 DEG(h) = l;          DEG(h) = l;
   
                                 for ( k = 0, dyk = 0; k < dx; k++, dyk += dy ) {          for ( k = 0, dyk = 0; k < dx; k++, dyk += dy ) {
                                         tmulup(rx[k],h,dy,&ut);            tmulup(rx[k],h,dy,&ut);
                                         if ( ut )            if ( ut )
                                                 for ( l = 0; l < dy; l++ ) {              for ( l = 0; l < dy; l++ ) {
                                                         s = (GFS)COEF(ut)[l];                s = (GFS)COEF(ut)[l];
                                                         if ( s ) {                if ( s ) {
                                                                 u = CONT(s);                  u = CONT(s);
                                                                 r[dyk+l] = _ADDSF(r[dyk+l],FTOIF(u));                  r[dyk+l] = _ADDSF(r[dyk+l],FTOIF(u));
                                                         }                }
                                                 }              }
                                 }          }
                         } else {        } else {
                                 prev = nf[d-1][i-1];          prev = nf[d-1][i-1];
                                 for ( k = 0, dyk = 0; k < dx; k++, dyk += dy ) {          for ( k = 0, dyk = 0; k < dx; k++, dyk += dy ) {
                                         for ( l = 1; l < dy; l++ )            for ( l = 1; l < dy; l++ )
                                                 r[dyk+l] = prev[dyk+l-1];              r[dyk+l] = prev[dyk+l-1];
                                 }          }
                         }        }
                         nf[d][i] = r;        nf[d][i] = r;
                         hist = (int *)CALLOC(len,sizeof(int));        hist = (int *)CALLOC(len,sizeof(int));
                         hist[len0+i] = _onesf();        hist[len0+i] = _onesf();
                         r1 = (int *)CALLOC(dxdy,sizeof(int));        r1 = (int *)CALLOC(dxdy,sizeof(int));
                         bcopy(r,r1,dxdy*sizeof(int));        bcopy(r,r1,dxdy*sizeof(int));
                         lnf(r1,hist,dxdy,len,g);        lnf(r1,hist,dxdy,len,g);
                         for ( k = 0; k < dxdy && !r1[k]; k++ );        for ( k = 0; k < dxdy && !r1[k]; k++ );
                         if ( k == dxdy ) {        if ( k == dxdy ) {
                                 f = 0;          f = 0;
                                 for ( k = j = 0; k <= d; k++ )          for ( k = j = 0; k <= d; k++ )
                                         for ( i = 0; i <= k; i++, j++ )            for ( i = 0; i <= k; i++, j++ )
                                                 if ( hist[j] ) {              if ( hist[j] ) {
                                                         iftogfs(hist[j],&s);                iftogfs(hist[j],&s);
                                                         /* mono = s*x^(k-i)*y^i */                /* mono = s*x^(k-i)*y^i */
                                                         create_bmono((P)s,x,k-i,y,i,&mono);                create_bmono((P)s,x,k-i,y,i,&mono);
                                                         addp(vl,f,mono,&f1);                addp(vl,f,mono,&f1);
                                                         f = f1;                f = f1;
                                                 }              }
                                 *fr = f;          *fr = f;
                                 return;          return;
                         }       else {        }  else {
                                 lb = (struct lb *)CALLOC(1,sizeof(struct lb));          lb = (struct lb *)CALLOC(1,sizeof(struct lb));
                                 lb->pos = k;          lb->pos = k;
                                 lb->r = r1;          lb->r = r1;
                                 lb->hist = hist;          lb->hist = hist;
                                 lb->len = len;          lb->len = len;
                                 g = insert_lb(g,lb);          g = insert_lb(g,lb);
                         }        }
                 }      }
         }    }
 }  }
   
 static void create_bmono(P c,V x,int i,V y,int j,P *mono)  static void create_bmono(P c,V x,int i,V y,int j,P *mono)
 {  {
         P t,s;    P t,s;
   
         if ( !i )    if ( !i )
                 if ( !j )      if ( !j )
                         t = c;        t = c;
                 else {      else {
                         /* c*y^j */        /* c*y^j */
                         MKV(y,t);        MKV(y,t);
                         COEF(DC(t)) = c;        COEF(DC(t)) = c;
                         STOQ(j,DEG(DC(t)));        STOQ(j,DEG(DC(t)));
                 }      }
         else if ( !j ) {    else if ( !j ) {
                 /* c*x^i */      /* c*x^i */
                 MKV(x,t);      MKV(x,t);
                 COEF(DC(t)) = c;      COEF(DC(t)) = c;
                 STOQ(i,DEG(DC(t)));      STOQ(i,DEG(DC(t)));
         } else {    } else {
                 MKV(y,s);      MKV(y,s);
                 COEF(DC(s)) = c;      COEF(DC(s)) = c;
                 STOQ(j,DEG(DC(s)));      STOQ(j,DEG(DC(s)));
                 MKV(x,t);      MKV(x,t);
                 COEF(DC(t)) = s;      COEF(DC(t)) = s;
                 STOQ(i,DEG(DC(t)));      STOQ(i,DEG(DC(t)));
         }    }
         *mono = t;    *mono = t;
 }  }

Legend:
Removed from v.1.15  
changed lines
  Added in v.1.20

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>