[BACK]Return to mp.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / pari-2.2 / src / kernel / none

Diff for /OpenXM_contrib/pari-2.2/src/kernel/none/Attic/mp.c between version 1.1 and 1.2

version 1.1, 2001/10/02 11:17:09 version 1.2, 2002/09/11 07:27:00
Line 30  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
Line 30  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  */   */
   
 /* z2 := z1[imin..imax].f shifted left sh bits (feeding f from the right) */  /* z2 := z1[imin..imax].f shifted left sh bits (feeding f from the right) */
   /* These macros work only for sh != 0 !!! */
 #define shift_left2(z2,z1,imin,imax,f, sh,m) {\  #define shift_left2(z2,z1,imin,imax,f, sh,m) {\
   register ulong _l,_k = ((ulong)f)>>m;\    register ulong _l,_k = ((ulong)f)>>m;\
   GEN t1 = z1 + imax, t2 = z2 + imax, T = z1 + imin;\    GEN t1 = z1 + imax, t2 = z2 + imax, T = z1 + imin;\
Line 45  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
Line 46  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
   shift_left2((z2),(z1),(imin),(imax),(f),(sh),(_m));\    shift_left2((z2),(z1),(imin),(imax),(f),(sh),(_m));\
 }  }
   
 /* z2 := f.z1[imin..imax-1] shifted right sh bits (feeding f from the left) */  #define shift_words_r(target,source,source_end,prepend, sh, sh_complement) {\
 #define shift_right2(z2,z1,imin,imax,f, sh,m) {\    register ulong _k,_l = *source++;\
   register GEN t1 = z1 + imin, t2 = z2 + imin, T = z1 + imax;\    *target++ = (_l>>(ulong)sh) | ((ulong)prepend<<(ulong)sh_complement);\
   register ulong _k,_l = *t1++;\    while (source < source_end) {\
   *t2++ = (_l>>(ulong)sh) | ((ulong)f<<(ulong)m);\      _k = _l<<(ulong)sh_complement; _l = *source++;\
   while (t1 < T) {\      *target++ = (_l>>(ulong)sh) | _k;\
     _k = _l<<(ulong)m; _l = *t1++;\  
     *t2++ = (_l>>(ulong)sh) | _k;\  
   }\    }\
 }  }
   #define shift_right2(z2,z1,imin,imax,f, sh,m) {\
     register GEN s = (z1) + (imin), ta = (z2) + (imin), se = (z1) + (imax);\
     shift_words_r(ta,s,se,(f),(sh),(m));                          \
   }
   /* z2 := f.z1[imin..imax-1] shifted right sh bits (feeding f from the left) */
 #define shift_right(z2,z1,imin,imax,f, sh) {\  #define shift_right(z2,z1,imin,imax,f, sh) {\
   register const ulong _m = BITS_IN_LONG - sh;\    register const ulong _m = BITS_IN_LONG - (sh);\
   shift_right2((z2),(z1),(imin),(imax),(f),(sh),(_m));\    shift_right2((z2),(z1),(imin),(imax),(f),(sh),(_m));\
 }  }
   
Line 95  addsispec(long s, GEN x, long nx)
Line 99  addsispec(long s, GEN x, long nx)
   while (xd > x) *--zd = *--xd;    while (xd > x) *--zd = *--xd;
   *--zd = evalsigne(1) | evallgefint(lz);    *--zd = evalsigne(1) | evallgefint(lz);
   *--zd = evaltyp(t_INT) | evallg(lz);    *--zd = evaltyp(t_INT) | evallg(lz);
   avma=(long)zd; return zd;    avma=(gpmem_t)zd; return zd;
 }  }
   
 #define swapspec(x,y, nx,ny) {long _a=nx;GEN _z=x; nx=ny; ny=_a; x=y; y=_z;}  #define swapspec(x,y, nx,ny) {long _a=nx;GEN _z=x; nx=ny; ny=_a; x=y; y=_z;}
Line 129  addiispec(GEN x, GEN y, long nx, long ny)
Line 133  addiispec(GEN x, GEN y, long nx, long ny)
   while (xd > x) *--zd = *--xd;    while (xd > x) *--zd = *--xd;
   *--zd = evalsigne(1) | evallgefint(lz);    *--zd = evalsigne(1) | evallgefint(lz);
   *--zd = evaltyp(t_INT) | evallg(lz);    *--zd = evaltyp(t_INT) | evallg(lz);
   avma=(long)zd; return zd;    avma=(gpmem_t)zd; return zd;
 }  }
   
 /* assume x >= y */  /* assume x >= y */
Line 158  subisspec(GEN x, long s, long nx)
Line 162  subisspec(GEN x, long s, long nx)
     do  *--zd = *--xd; while (xd > x);      do  *--zd = *--xd; while (xd > x);
   *--zd = evalsigne(1) | evallgefint(lz);    *--zd = evalsigne(1) | evallgefint(lz);
   *--zd = evaltyp(t_INT) | evallg(lz);    *--zd = evaltyp(t_INT) | evallg(lz);
   avma=(long)zd; return zd;    avma=(gpmem_t)zd; return zd;
 }  }
   
 /* assume x > y */  /* assume x > y */
Line 191  subiispec(GEN x, GEN y, long nx, long ny)
Line 195  subiispec(GEN x, GEN y, long nx, long ny)
     do  *--zd = *--xd; while (xd > x);      do  *--zd = *--xd; while (xd > x);
   *--zd = evalsigne(1) | evallgefint(lz);    *--zd = evalsigne(1) | evallgefint(lz);
   *--zd = evaltyp(t_INT) | evallg(lz);    *--zd = evaltyp(t_INT) | evallg(lz);
   avma=(long)zd; return zd;    avma=(gpmem_t)zd; return zd;
 }  }
   
 #ifndef __M68K__  #ifndef __M68K__
   
 /* prototype of positive small ints */  /* prototype of positive small ints */
 static long pos_s[] = {  static long pos_s[] = {
   evaltyp(t_INT) | m_evallg(3), evalsigne(1) | evallgefint(3), 0 };    evaltyp(t_INT) | _evallg(3), evalsigne(1) | evallgefint(3), 0 };
   
 /* prototype of negative small ints */  /* prototype of negative small ints */
 static long neg_s[] = {  static long neg_s[] = {
   evaltyp(t_INT) | m_evallg(3), evalsigne(-1) | evallgefint(3), 0 };    evaltyp(t_INT) | _evallg(3), evalsigne(-1) | evallgefint(3), 0 };
   
 void  void
 affir(GEN x, GEN y)  affir(GEN x, GEN y)
Line 213  affir(GEN x, GEN y)
Line 217  affir(GEN x, GEN y)
   if (!s)    if (!s)
   {    {
     y[1] = evalexpo(-bit_accuracy(ly));      y[1] = evalexpo(-bit_accuracy(ly));
     y[2] = 0; return;      return;
   }    }
   
   lx=lgefint(x); sh=bfffo(x[2]);    lx = lgefint(x); sh = bfffo(x[2]);
   y[1] = evalsigne(s) | evalexpo(bit_accuracy(lx)-sh-1);    y[1] = evalsigne(s) | evalexpo(bit_accuracy(lx)-sh-1);
   if (sh)    if (sh) {
   {  
     if (lx>ly)      if (lx>ly)
     { shift_left(y,x,2,ly-1, x[ly],sh); }      { shift_left(y,x,2,ly-1, x[ly],sh); }
     else      else
Line 227  affir(GEN x, GEN y)
Line 230  affir(GEN x, GEN y)
       for (i=lx; i<ly; i++) y[i]=0;        for (i=lx; i<ly; i++) y[i]=0;
       shift_left(y,x,2,lx-1, 0,sh);        shift_left(y,x,2,lx-1, 0,sh);
     }      }
     return;  
   }    }
     else {
   if (lx>=ly)      if (lx>=ly)
     for (i=2; i<ly; i++) y[i]=x[i];        for (i=2; i<ly; i++) y[i]=x[i];
   else      else
   {      {
     for (i=2; i<lx; i++) y[i]=x[i];        for (i=2; i<lx; i++) y[i]=x[i];
     for (   ; i<ly; i++) y[i]=0;        for (   ; i<ly; i++) y[i]=0;
       }
   }    }
 }  }
   
Line 244  affrr(GEN x, GEN y)
Line 247  affrr(GEN x, GEN y)
 {  {
   long lx,ly,i;    long lx,ly,i;
   
   y[1] = x[1]; if (!signe(x)) { y[2]=0; return; }    y[1] = x[1]; if (!signe(x)) return;
   
   lx=lg(x); ly=lg(y);    lx=lg(x); ly=lg(y);
   if (lx>=ly)    if (lx>=ly)
Line 259  affrr(GEN x, GEN y)
Line 262  affrr(GEN x, GEN y)
 GEN  GEN
 shifti(GEN x, long n)  shifti(GEN x, long n)
 {  {
     return shifti3(x, n, 0);
   }
   
   GEN
   shifti3(GEN x, long n, long flag)
   {
   const long s=signe(x);    const long s=signe(x);
   long lx,ly,i,m;    long lx,ly,i,m;
   GEN y;    GEN y;
Line 285  shifti(GEN x, long n)
Line 294  shifti(GEN x, long n)
   }    }
   else    else
   {    {
       long lyorig;
   
       if (s > 0)
           flag = 0;
     n = -n;      n = -n;
     ly = lx - (n>>TWOPOTBITS_IN_LONG);      ly = lyorig = lx - (n>>TWOPOTBITS_IN_LONG);
     if (ly<3) return gzero;      if (ly<3)
           return flag ? stoi(-1) : gzero;
     y = new_chunk(ly);      y = new_chunk(ly);
     m = n & (BITS_IN_LONG-1);      m = n & (BITS_IN_LONG-1);
     if (!m) for (i=2; i<ly; i++) y[i]=x[i];      if (m) {
     else  
     {  
       shift_right(y,x, 2,ly, 0,m);        shift_right(y,x, 2,ly, 0,m);
       if (y[2] == 0)        if (y[2] == 0)
       {        {
         if (ly==3) { avma = (long)(y+3); return gzero; }          if (ly==3) { avma = (gpmem_t)(y+3); return flag ? stoi(-1) : gzero; }
         ly--; avma = (long)(++y);          ly--; avma = (gpmem_t)(++y);
       }        }
       } else {
         for (i=2; i<ly; i++) y[i]=x[i];
     }      }
       /* With FLAG: round up instead of rounding down */
       if (flag) {                         /* Check low bits of x */
         i = lx; flag = 0;
         while (--i >= lyorig)
           if (x[i]) { flag = 1; break; }  /* Need to increment y by 1 */
         if (!flag && m)
           flag = x[lyorig - 1] & ((1<<m) - 1);
       }
       if (flag) {                         /* Increment y */
         for (i = ly;;)
         {
           if (--i < 2)
           { /* Extend y on the left */
             if (avma <= bot) err(errpile);
             avma = (gpmem_t)(--y); ly++;
             y[2] = 1; break;
           }
           if (++y[i]) break;
           /* Now propagate the bit into the next longword */
         }
       }
   }    }
   y[1] = evalsigne(s)|evallgefint(ly);    y[1] = evalsigne(s)|evallgefint(ly);
   y[0] = evaltyp(t_INT)|evallg(ly); return y;    y[0] = evaltyp(t_INT)|evallg(ly); return y;
Line 524  addir(GEN x, GEN y)
Line 559  addir(GEN x, GEN y)
 #else /* design issue: make 0.0 "absorbing" */  #else /* design issue: make 0.0 "absorbing" */
     if (e>0) return rcopy(y);      if (e>0) return rcopy(y);
 #endif  #endif
     z = cgetr(3 + ((-e)>>TWOPOTBITS_IN_LONG));      return itor(x, 3 + ((-e)>>TWOPOTBITS_IN_LONG));
     affir(x,z); return z;  
   }    }
   
   ly=lg(y);    ly=lg(y);
   if (e>0)    if (e > 0)
   {    {
     l = ly - (e>>TWOPOTBITS_IN_LONG);      l = ly - (e>>TWOPOTBITS_IN_LONG);
     if (l<3) return rcopy(y);      if (l<3) return rcopy(y);
   }    }
   else l = ly + ((-e)>>TWOPOTBITS_IN_LONG)+1;    else l = ly + ((-e)>>TWOPOTBITS_IN_LONG)+1;
   z=cgetr(l); affir(x,z); y=addrr(z,y);    y = addrr(itor(x,l), y);
   z = y+l; ly = lg(y); while (ly--) z[ly] = y[ly];    z = y+l; ly = lg(y); while (ly--) z[ly] = y[ly];
   avma=(long)z; return z;    avma=(gpmem_t)z; return z;
 }  }
   
 GEN  GEN
Line 554  addrr(GEN x, GEN y)
Line 588  addrr(GEN x, GEN y)
     if (!sx)      if (!sx)
     {      {
       if (e > 0) ex=ey;        if (e > 0) ex=ey;
       z=cgetr(3); z[1]=evalexpo(ex); z[2]=0; return z;        return realzero_bit(ex);
     }      }
     if (e > 0) { z=cgetr(3); z[1]=evalexpo(ey); z[2]=0; return z; }      if (e > 0) return realzero_bit(ey);
     lz = 3 + ((-e)>>TWOPOTBITS_IN_LONG);      lz = 3 + ((-e)>>TWOPOTBITS_IN_LONG);
     lx = lg(x); if (lz>lx) lz=lx;      lx = lg(x); if (lz>lx) lz=lx;
     z = cgetr(lz); while(--lz) z[lz]=x[lz];      z = cgetr(lz); while(--lz) z[lz]=x[lz];
Line 564  addrr(GEN x, GEN y)
Line 598  addrr(GEN x, GEN y)
   }    }
   if (!sx)    if (!sx)
   {    {
     if (e < 0) { z=cgetr(3); z[1]=evalexpo(ex); z[2]=0; return z; }      if (e < 0) return realzero_bit(ex);
     lz = 3 + (e>>TWOPOTBITS_IN_LONG);      lz = 3 + (e>>TWOPOTBITS_IN_LONG);
     ly = lg(y); if (lz>ly) lz=ly;      ly = lg(y); if (lz>ly) lz=ly;
     z = cgetr(lz); while (--lz) z[lz]=y[lz];      z = cgetr(lz); while (--lz) z[lz]=y[lz];
Line 602  addrr(GEN x, GEN y)
Line 636  addrr(GEN x, GEN y)
   
   if (sx==sy)    if (sx==sy)
   { /* addition */    { /* addition */
     i=lz-1; avma = (long)z;      i=lz-1; avma = (gpmem_t)z;
     if (flag==0) { z[i] = y[i]; i--; }      if (flag==0) { z[i] = y[i]; i--; }
     overflow=0;      overflow=0;
     for (j=lx-1; j>=2; i--,j--) z[i] = addllx(x[j],y[i]);      for (j=lx-1; j>=2; i--,j--) z[i] = addllx(x[j],y[i]);
Line 613  addrr(GEN x, GEN y)
Line 647  addrr(GEN x, GEN y)
         if (i==1)          if (i==1)
         {          {
           shift_right(z,z, 2,lz, 1,1);            shift_right(z,z, 2,lz, 1,1);
           ey=evalexpo(ey+1);            z[1] = evalsigne(sx) | evalexpo(ey+1); return z;
           z[1] = evalsigne(sx) | ey; return z;  
         }          }
         z[i] = y[i]+1; if (z[i--]) break;          z[i] = y[i]+1; if (z[i--]) break;
       }        }
Line 629  addrr(GEN x, GEN y)
Line 662  addrr(GEN x, GEN y)
     i=2; while (i<lx && x[i]==y[i]) i++;      i=2; while (i<lx && x[i]==y[i]) i++;
     if (i==lx)      if (i==lx)
     {      {
       avma = (long)(z+lz);        avma = (gpmem_t)(z+lz);
       e = evalexpo(ey - bit_accuracy(lx));        return realzero_bit(ey - bit_accuracy(lx));
       z=cgetr(3); z[1]=e; z[2]=0; return z;  
     }      }
     f2 = ((ulong)y[i] > (ulong)x[i]);      f2 = ((ulong)y[i] > (ulong)x[i]);
   }    }
Line 656  addrr(GEN x, GEN y)
Line 688  addrr(GEN x, GEN y)
   
   x = z+2; i=0; while (!x[i]) i++;    x = z+2; i=0; while (!x[i]) i++;
   lz -= i; z += i; m = bfffo(z[2]);    lz -= i; z += i; m = bfffo(z[2]);
   e = evalexpo(ey - (m | (i<<TWOPOTBITS_IN_LONG)));  
   if (m) shift_left(z,z,2,lz-1, 0,m);    if (m) shift_left(z,z,2,lz-1, 0,m);
   z[1] = evalsigne(sx) | e;    z[1] = evalsigne(sx) | evalexpo(ey - (m | (i<<TWOPOTBITS_IN_LONG)));
   z[0] = evaltyp(t_REAL) | evallg(lz);    z[0] = evaltyp(t_REAL) | evallg(lz);
   avma = (long)z; return z;    avma = (gpmem_t)z; return z;
 }  }
   
 GEN  GEN
Line 719  mulsispec(long x, GEN y, long ny)
Line 750  mulsispec(long x, GEN y, long ny)
   if (hiremainder) *--z = hiremainder; else lz--;    if (hiremainder) *--z = hiremainder; else lz--;
   *--z = evalsigne(1) | evallgefint(lz);    *--z = evalsigne(1) | evallgefint(lz);
   *--z = evaltyp(t_INT) | evallg(lz);    *--z = evaltyp(t_INT) | evallg(lz);
   avma=(long)z; return z;    avma=(gpmem_t)z; return z;
 }  }
   
 GEN  GEN
Line 755  addsmulsi(long a, long b, GEN Y)
Line 786  addsmulsi(long a, long b, GEN Y)
   if (hiremainder) *--z = hiremainder; else lz--;    if (hiremainder) *--z = hiremainder; else lz--;
   *--z = evalsigne(1) | evallgefint(lz);    *--z = evalsigne(1) | evallgefint(lz);
   *--z = evaltyp(t_INT) | evallg(lz);    *--z = evaltyp(t_INT) | evallg(lz);
   avma=(long)z; return z;    avma=(gpmem_t)z; return z;
 }  }
   
 #ifndef __M68K__  #ifndef __M68K__
Line 784  mulsr(long x, GEN y)
Line 815  mulsr(long x, GEN y)
   if (!s)    if (!s)
   {    {
     if (x<0) x = -x;      if (x<0) x = -x;
     e = y[1] + (BITS_IN_LONG-1)-bfffo(x);      e = expo(y) + (BITS_IN_LONG-1)-bfffo(x);
     if (e & ~EXPOBITS) err(muler2);      return realzero_bit(e);
     z=cgetr(3); z[1]=e; z[2]=0; return z;  
   }    }
   if (x<0) { s = -s; x = -x; }    if (x<0) { s = -s; x = -x; }
   if (x==1) { z=rcopy(y); setsigne(z,s); return z; }    if (x==1) { z=rcopy(y); setsigne(z,s); return z; }
Line 798  mulsr(long x, GEN y)
Line 828  mulsr(long x, GEN y)
   
   sh = bfffo(hiremainder); m = BITS_IN_LONG-sh;    sh = bfffo(hiremainder); m = BITS_IN_LONG-sh;
   if (sh) shift_left2(z,z, 2,lx-1, garde,sh,m);    if (sh) shift_left2(z,z, 2,lx-1, garde,sh,m);
   e = evalexpo(m+e);    z[1] = evalsigne(s) | evalexpo(m+e); return z;
   z[1] = evalsigne(s) | e; return z;  
 }  }
   
 GEN  GEN
Line 812  mulrr(GEN x, GEN y)
Line 841  mulrr(GEN x, GEN y)
   LOCAL_HIREMAINDER;    LOCAL_HIREMAINDER;
   LOCAL_OVERFLOW;    LOCAL_OVERFLOW;
   
   e = evalexpo(expo(x)+expo(y));    e = expo(x)+expo(y);
   if (!sx || !sy) { z=cgetr(3); z[1]=e; z[2]=0; return z; }    if (!sx || !sy) return realzero_bit(e);
   if (sy<0) sx = -sx;    if (sy<0) sx = -sx;
   lz=lg(x); ly=lg(y);    lz=lg(x); ly=lg(y);
   if (lz>ly) { lz=ly; z=x; x=y; y=z; flag=1; } else flag = (lz!=ly);    if (lz>ly) { lz=ly; z=x; x=y; y=z; flag=1; } else flag = (lz!=ly);
   z=cgetr(lz); z[1] = evalsigne(sx) | e;    z=cgetr(lz);
   if (lz==3)    if (lz==3)
   {    {
     if (flag)      if (flag)
Line 827  mulrr(GEN x, GEN y)
Line 856  mulrr(GEN x, GEN y)
     }      }
     else      else
       garde = mulll(x[2],y[2]);        garde = mulll(x[2],y[2]);
     if ((long)hiremainder<0) { z[2]=hiremainder; z[1]++; }      if ((long)hiremainder<0) { z[2]=hiremainder; e++; }
     else z[2]=(hiremainder<<1) | (garde>>(BITS_IN_LONG-1));      else z[2]=(hiremainder<<1) | (garde>>(BITS_IN_LONG-1));
       z[1] = evalsigne(sx)|evalexpo(e);
     return z;      return z;
   }    }
   
Line 865  mulrr(GEN x, GEN y)
Line 895  mulrr(GEN x, GEN y)
     z[i] = addll(addmul(p1,y1[i]), z[i]);      z[i] = addll(addmul(p1,y1[i]), z[i]);
   }    }
   z[2] = hiremainder+overflow;    z[2] = hiremainder+overflow;
   if (z[2] < 0) z[1]++;    if (z[2] < 0) e++; else shift_left(z,z,2,lzz,garde, 1);
   else    z[1] = evalsigne(sx) | evalexpo(e);
     shift_left(z,z,2,lzz,garde, 1);  
   return z;    return z;
 }  }
   
Line 883  mulir(GEN x, GEN y)
Line 912  mulir(GEN x, GEN y)
   if (!sx) return gzero;    if (!sx) return gzero;
   if (!is_bigint(x)) return mulsr(itos(x),y);    if (!is_bigint(x)) return mulsr(itos(x),y);
   sy=signe(y); ey=expo(y);    sy=signe(y); ey=expo(y);
   if (!sy)    if (!sy) return realzero_bit(expi(x)+ey);
   {  
     e = evalexpo(expi(x)+ey);  
     z=cgetr(3); z[1]=e; z[2]=0; return z;  
   }  
   
   if (sy<0) sx = -sx;    if (sy<0) sx = -sx;
   lz=lg(y); z=cgetr(lz);    lz=lg(y); z=cgetr(lz);
   y1=cgetr(lz+1);    y1 = itor(x, lz+1); x = y; y = y1;
   affir(x,y1); x=y; y=y1;    e = expo(y)+ey;
   e = evalexpo(expo(y)+ey);  
   z[1] = evalsigne(sx) | e;  
   if (lz==3)    if (lz==3)
   {    {
     (void)mulll(x[2],y[3]);      (void)mulll(x[2],y[3]);
     garde=addmul(x[2],y[2]);      garde=addmul(x[2],y[2]);
     if ((long)hiremainder < 0) { z[2]=hiremainder; z[1]++; }      if ((long)hiremainder < 0) { z[2]=hiremainder; e++; }
     else z[2]=(hiremainder<<1) | (garde>>(BITS_IN_LONG-1));      else z[2]=(hiremainder<<1) | (garde>>(BITS_IN_LONG-1));
     avma=(long)z; return z;      z[1] = evalsigne(sx) | evalexpo(e);
       avma=(gpmem_t)z; return z;
   }    }
   
   (void)mulll(x[2],y[lz]); garde=hiremainder;    (void)mulll(x[2],y[lz]); garde=hiremainder;
Line 937  mulir(GEN x, GEN y)
Line 961  mulir(GEN x, GEN y)
     z[i] = addll(addmul(p1,y1[i]), z[i]);      z[i] = addll(addmul(p1,y1[i]), z[i]);
   }    }
   z[2] = hiremainder+overflow;    z[2] = hiremainder+overflow;
   if (z[2] < 0) z[1]++;    if (z[2] < 0) e++;
   else    else
     shift_left(z,z,2,lzz,garde, 1);      shift_left(z,z,2,lzz,garde, 1);
   avma=(long)z; return z;    z[1] = evalsigne(sx) | evalexpo(e);
     avma=(gpmem_t)z; return z;
 }  }
   
 /* written by Bruno Haible following an idea of Robert Harley */  /* written by Bruno Haible following an idea of Robert Harley */
Line 974  modss(long x, long y)
Line 999  modss(long x, long y)
   
   if (!y) err(moder1);    if (!y) err(moder1);
   if (y<0) y=-y;    if (y<0) y=-y;
   hiremainder=0; divll(labs(x),y);    hiremainder=0; (void)divll(labs(x),y);
   if (!hiremainder) return gzero;    if (!hiremainder) return gzero;
   return (x < 0) ? stoi(y-hiremainder) : stoi(hiremainder);    return (x < 0) ? stoi(y-hiremainder) : stoi(hiremainder);
 }  }
Line 985  resss(long x, long y)
Line 1010  resss(long x, long y)
   LOCAL_HIREMAINDER;    LOCAL_HIREMAINDER;
   
   if (!y) err(reser1);    if (!y) err(reser1);
   hiremainder=0; divll(labs(x),labs(y));    hiremainder=0; (void)divll(labs(x),labs(y));
   return (x < 0) ? stoi(-((long)hiremainder)) : stoi(hiremainder);    return (x < 0) ? stoi(-((long)hiremainder)) : stoi(hiremainder);
 }  }
   
Line 1044  umodiu(GEN y, ulong x)
Line 1069  umodiu(GEN y, ulong x)
 GEN  GEN
 modiu(GEN y, ulong x) { return utoi(umodiu(y,x)); }  modiu(GEN y, ulong x) { return utoi(umodiu(y,x)); }
   
   /* return |y| \/ x */
   GEN
   diviu(GEN y, ulong x)
   {
     long sy=signe(y),ly,i;
     GEN z;
     LOCAL_HIREMAINDER;
   
     if (!x) err(diver4);
     if (!sy) { hiremainder=0; SAVE_HIREMAINDER; return gzero; }
   
     ly = lgefint(y);
     if (x <= (ulong)y[2]) hiremainder=0;
     else
     {
       if (ly==3) { hiremainder=itou(y); SAVE_HIREMAINDER; return gzero; }
       hiremainder=y[2]; ly--; y++;
     }
     z = cgeti(ly); z[1] = evallgefint(ly) | evalsigne(1);
     for (i=2; i<ly; i++) z[i]=divll(y[i],x);
     SAVE_HIREMAINDER; return z;
   }
   
 #ifndef __M68K__  #ifndef __M68K__
   
 GEN  GEN
Line 1096  divis(GEN y, long x)
Line 1144  divis(GEN y, long x)
 GEN  GEN
 divir(GEN x, GEN y)  divir(GEN x, GEN y)
 {  {
   GEN xr,z;    GEN z;
   long av,ly;    long ly;
     gpmem_t av;
   
   if (!signe(y)) err(diver5);    if (!signe(y)) err(diver5);
   if (!signe(x)) return gzero;    if (!signe(x)) return gzero;
   ly=lg(y); z=cgetr(ly); av=avma; xr=cgetr(ly+1); affir(x,xr);    ly = lg(y); z = cgetr(ly); av = avma;
   affrr(divrr(xr,y),z); avma=av; return z;    affrr(divrr(itor(x, ly+1), y), z);
     avma = av; return z;
 }  }
   
 GEN  GEN
 divri(GEN x, GEN y)  divri(GEN x, GEN y)
 {  {
   GEN yr,z;    long lx, s = signe(y);
   long av,lx,s=signe(y);    gpmem_t av;
     GEN z;
   
   if (!s) err(diver8);    if (!s) err(diver8);
   if (!signe(x))    if (!signe(x)) return realzero_bit(expo(x) - expi(y));
   {  
     const long ex = evalexpo(expo(x) - expi(y));  
   
     if (ex<0) err(diver12);  
     z=cgetr(3); z[1]=ex; z[2]=0; return z;  
   }  
   if (!is_bigint(y)) return divrs(x, s>0? y[2]: -y[2]);    if (!is_bigint(y)) return divrs(x, s>0? y[2]: -y[2]);
   
   lx=lg(x); z=cgetr(lx);    lx = lg(x); z = cgetr(lx); av = avma;
   av=avma; yr=cgetr(lx+1); affir(y,yr);    affrr(divrr(x, itor(y, lx+1)), z);
   affrr(divrr(x,yr),z); avma=av; return z;    avma = av; return z;
 }  }
   
 void  void
 diviiz(GEN x, GEN y, GEN z)  diviiz(GEN x, GEN y, GEN z)
 {  {
   long av=avma,lz;    gpmem_t av = avma;
   GEN xr,yr;    if (typ(z) == t_INT) affii(divii(x,y), z);
     else {
   if (typ(z)==t_INT) { affii(divii(x,y),z); avma=av; return; }      long lz = lg(z);
   lz=lg(z); xr=cgetr(lz); affir(x,xr); yr=cgetr(lz); affir(y,yr);      affrr(divrr(itor(x,lz), itor(y,lz)), z);
   affrr(divrr(xr,yr),z); avma=av;    }
     avma = av;
 }  }
   
 void  void
 mpdivz(GEN x, GEN y, GEN z)  mpdivz(GEN x, GEN y, GEN z)
 {  {
   long av=avma;    gpmem_t av = avma;
     GEN r;
   
   if (typ(z)==t_INT)    if (typ(z)==t_INT)
   {    {
     if (typ(x)==t_REAL || typ(y)==t_REAL) err(divzer1);      if (typ(x) == t_REAL || typ(y) == t_REAL) err(divzer1);
     affii(divii(x,y),z); avma=av; return;      affii(divii(x,y), z);
       avma = av; return;
   }    }
   if (typ(x)==t_INT)  
   {  
     GEN xr,yr;  
     long lz;  
   
     if (typ(y)==t_REAL) { affrr(divir(x,y),z); avma=av; return; }    if (typ(x) == t_INT)
     lz=lg(z); xr=cgetr(lz); affir(x,xr); yr=cgetr(lz); affir(y,yr);    {
     affrr(divrr(xr,yr),z); avma=av; return;      if (typ(y) == t_REAL)
         r = divir(x,y);
       else
       {
         long lz = lg(z);
         r = divrr(itor(x,lz), itor(y,lz));
       }
   }    }
   if (typ(y)==t_REAL) { affrr(divrr(x,y),z); avma=av; return; }    else if (typ(y) == t_REAL)
   affrr(divri(x,y),z); avma=av;      r = divrr(x,y);
     else
       r = divri(x,y);
     affrr(r, z); avma = av;
 }  }
   
 GEN  GEN
 divsr(long x, GEN y)  divsr(long x, GEN y)
 {  {
   long av,ly;    gpmem_t av;
   GEN p1,z;    long ly;
     GEN z;
   
   if (!signe(y)) err(diver3);    if (!signe(y)) err(diver3);
   if (!x) return gzero;    if (!x) return gzero;
   ly=lg(y); z=cgetr(ly); av=avma;    ly = lg(y); z = cgetr(ly); av = avma;
   p1=cgetr(ly+1); affsr(x,p1); affrr(divrr(p1,y),z);    affrr(divrr(stor(x,ly+1), y), z);
   avma=av; return z;    avma = av; return z;
 }  }
   
 GEN  GEN
Line 1182  modii(GEN x, GEN y)
Line 1236  modii(GEN x, GEN y)
     case 1: return resii(x,y);      case 1: return resii(x,y);
     default:      default:
     {      {
       long av = avma;        gpmem_t av = avma;
       (void)new_chunk(lgefint(y));        (void)new_chunk(lgefint(y));
       x = resii(x,y); avma=av;        x = resii(x,y); avma=av;
       if (x==gzero) return x;        if (x==gzero) return x;
Line 1194  modii(GEN x, GEN y)
Line 1248  modii(GEN x, GEN y)
 void  void
 modiiz(GEN x, GEN y, GEN z)  modiiz(GEN x, GEN y, GEN z)
 {  {
   const long av = avma;    const gpmem_t av = avma;
   affii(modii(x,y),z); avma=av;    affii(modii(x,y),z); avma=av;
 }  }
   
 GEN  GEN
 divrs(GEN x, long y)  divrs(GEN x, long y)
 {  {
   long i,lx,ex,garde,sh,s=signe(x);    long i,lx,garde,sh,s=signe(x);
   GEN z;    GEN z;
   LOCAL_HIREMAINDER;    LOCAL_HIREMAINDER;
   
   if (!y) err(diver6);    if (!y) err(diver6);
   if (!s)    if (!s) return realzero_bit(expo(x) - (BITS_IN_LONG-1)+bfffo(y));
   {  
     z=cgetr(3); z[1] = x[1] - (BITS_IN_LONG-1) + bfffo(y);  
     if (z[1]<0) err(diver7);  
     z[2]=0; return z;  
   }  
   if (y<0) { s = -s; y = -y; }    if (y<0) { s = -s; y = -y; }
   if (y==1) { z=rcopy(x); setsigne(z,s); return z; }    if (y==1) { z=rcopy(x); setsigne(z,s); return z; }
   
Line 1219  divrs(GEN x, long y)
Line 1268  divrs(GEN x, long y)
   for (i=2; i<lx; i++) z[i]=divll(x[i],y);    for (i=2; i<lx; i++) z[i]=divll(x[i],y);
   
   /* we may have hiremainder != 0 ==> garde */    /* we may have hiremainder != 0 ==> garde */
   garde=divll(0,y); sh=bfffo(z[2]); ex=evalexpo(expo(x)-sh);    garde=divll(0,y); sh=bfffo(z[2]);
   if (sh) shift_left(z,z, 2,lx-1, garde,sh);    if (sh) shift_left(z,z, 2,lx-1, garde,sh);
   z[1] = evalsigne(s) | ex;    z[1] = evalsigne(s) | evalexpo(expo(x)-sh);
   return z;    return z;
 }  }
   
Line 1233  divrr(GEN x, GEN y)
Line 1282  divrr(GEN x, GEN y)
   GEN z,x1;    GEN z,x1;
   
   if (!sy) err(diver9);    if (!sy) err(diver9);
   e = evalexpo(expo(x) - expo(y));    e = expo(x) - expo(y);
   if (!sx) { z=cgetr(3); z[1]=e; z[2]=0; return z; }    if (!sx) return realzero_bit(e);
   if (sy<0) sx = -sx;    if (sy<0) sx = -sx;
   e = evalsigne(sx) | e;  
   lx=lg(x); ly=lg(y);    lx=lg(x); ly=lg(y);
   if (ly==3)    if (ly==3)
   {    {
Line 1248  divrr(GEN x, GEN y)
Line 1296  divrr(GEN x, GEN y)
       l >>= 1; if (k&1) l |= HIGHBIT;        l >>= 1; if (k&1) l |= HIGHBIT;
       k >>= 1;        k >>= 1;
     }      }
     z=cgetr(3); z[1]=e;      z = cgetr(3); z[1] = evalsigne(sx) | evalexpo(e);
     hiremainder=k; z[2]=divll(l,y[2]); return z;      hiremainder=k; z[2]=divll(l,y[2]); return z;
   }    }
   
   lz = (lx<=ly)? lx: ly;    lz = min(lx,ly); z = new_chunk(lz);
   x1 = (z=new_chunk(lz))-1;    x1 = z-1;
   x1[1]=0; for (i=2; i<lz; i++) x1[i]=x[i];    x1[1]=0; for (i=2; i<lz; i++) x1[i]=x[i];
   x1[lz] = (lx>lz)? x[lz]: 0;    x1[lz] = (lx>lz)? x[lz]: 0;
   x=z; si=y[2]; saux=y[3];    x=z; si=y[2]; saux=y[3];
Line 1313  divrr(GEN x, GEN y)
Line 1361  divrr(GEN x, GEN y)
   }    }
   x1 = x-1; for (j=lz-1; j>=2; j--) x[j]=x1[j];    x1 = x-1; for (j=lz-1; j>=2; j--) x[j]=x1[j];
   if (*x) { shift_right(x,x, 2,lz, 1,1); } else e--;    if (*x) { shift_right(x,x, 2,lz, 1,1); } else e--;
   x[0]=evaltyp(t_REAL)|evallg(lz);    x[0] = evaltyp(t_REAL)|evallg(lz);
   x[1]=e; return x;    x[1] = evalsigne(sx) | evalexpo(e);
     return x;
 }  }
 #endif /* !defined(__M68K__) */  #endif /* !defined(__M68K__) */
   
 /* The following ones are not in mp.s (mulii is, with a different algorithm) */  /* The following ones are not in mp.s (mulii is, with a different algorithm) */
   
 /* Integer division x / y:  /* Integer division x / y: such that sign(r) = sign(x)
  *   if z = ONLY_REM return remainder, otherwise return quotient   *   if z = ONLY_REM return remainder, otherwise return quotient
  *   if z != NULL set *z to remainder   *   if z != NULL set *z to remainder
  *   *z is the last object on stack (and thus can be disposed of with cgiv   *   *z is the last object on stack (and thus can be disposed of with cgiv
Line 1331  GEN
Line 1381  GEN
 dvmdii(GEN x, GEN y, GEN *z)  dvmdii(GEN x, GEN y, GEN *z)
 {  {
   long sx=signe(x),sy=signe(y);    long sx=signe(x),sy=signe(y);
   long av,lx,ly,lz,i,j,sh,k,lq,lr;    long lx, ly, lz, i, j, sh, k, lq, lr;
     gpmem_t av;
   ulong si,qp,saux, *xd,*rd,*qd;    ulong si,qp,saux, *xd,*rd,*qd;
   GEN r,q,x1;    GEN r,q,x1;
   
Line 1447  DIVIDE: /* quotient is non-zero */
Line 1498  DIVIDE: /* quotient is non-zero */
     while (lz--) *--qd = *--xd;      while (lz--) *--qd = *--xd;
     *--qd = evalsigne(sy) | evallgefint(lq);      *--qd = evalsigne(sy) | evallgefint(lq);
     *--qd = evaltyp(t_INT) | evallg(lq);      *--qd = evaltyp(t_INT) | evallg(lq);
     avma = (long)qd; return (GEN)qd;      avma = (gpmem_t)qd; return (GEN)qd;
   }    }
   
   j=lq; while (j<lx && !x[j]) j++;    j=lq; while (j<lx && !x[j]) j++;
Line 1473  DIVIDE: /* quotient is non-zero */
Line 1524  DIVIDE: /* quotient is non-zero */
     }      }
     *--rd = evalsigne(sx) | evallgefint(lr);      *--rd = evalsigne(sx) | evallgefint(lr);
     *--rd = evaltyp(t_INT) | evallg(lr);      *--rd = evaltyp(t_INT) | evallg(lr);
     avma = (long)rd; return (GEN)rd;      avma = (gpmem_t)rd; return (GEN)rd;
   }    }
   
   lr = lz+2;    lr = lz+2;
Line 1517  DIVIDE: /* quotient is non-zero */
Line 1568  DIVIDE: /* quotient is non-zero */
     while (lr--) *--qd = *--rd;      while (lr--) *--qd = *--rd;
     *z = (GEN)qd;      *z = (GEN)qd;
   }    }
   avma = (long)qd; return q;    avma = (gpmem_t)qd; return q;
 }  }
   
 /* assume y > x > 0. return y mod x */  /* assume y > x > 0. return y mod x */
Line 1536  mppgcd_resiu(GEN y, ulong x)
Line 1587  mppgcd_resiu(GEN y, ulong x)
 void  void
 mppgcd_plus_minus(GEN x, GEN y, GEN res)  mppgcd_plus_minus(GEN x, GEN y, GEN res)
 {  {
   long av = avma;    gpmem_t av = avma;
   long lx = lgefint(x)-1;    long lx = lgefint(x)-1;
   long ly = lgefint(y)-1, lt,m,i;    long ly = lgefint(y)-1, lt,m,i;
   GEN t;    GEN t;
Line 1579  smulss(ulong x, ulong y, ulong *rem)
Line 1630  smulss(ulong x, ulong y, ulong *rem)
 #  define DIVCONVI 14  #  define DIVCONVI 14
 #endif  #endif
   
 /* convert integer --> base 10^9 [assume x > 0] */  /* convert integer --> base 10^9 */
 GEN  GEN
 convi(GEN x)  convi(GEN x)
 {  {
   ulong av=avma, lim;    gpmem_t av = avma, lim;
   long lz;    long lz;
   GEN z,p1;    GEN z,p1;
   
   lz = 3 + ((lgefint(x)-2)*15)/DIVCONVI;    lz = 3 + ((lgefint(x)-2)*15)/DIVCONVI;
   z=new_chunk(lz); z[1] = -1; p1 = z+2;    z = new_chunk(lz); z[1] = -1; p1 = z+2;
   lim = stack_lim(av,1);    lim = stack_lim(av,1);
   for(;;)    for(;;)
   {    {
     x = divis(x,1000000000); *p1++ = hiremainder;      x = diviu(x,1000000000); *p1++ = hiremainder;
     if (!signe(x)) { avma=av; return p1; }      if (!signe(x)) { avma = av; return p1; }
     if (low_stack(lim, stack_lim(av,1))) x = gerepileuptoint((long)z,x);      if (low_stack(lim, stack_lim(av,1))) x = gerepileuptoint((gpmem_t)z,x);
   }    }
 }  }
 #undef DIVCONVI  #undef DIVCONVI
Line 1641  confrac(GEN x)
Line 1692  confrac(GEN x)
 GEN  GEN
 truedvmdii(GEN x, GEN y, GEN *z)  truedvmdii(GEN x, GEN y, GEN *z)
 {  {
   long av = avma;    gpmem_t av = avma;
   GEN r, q = dvmdii(x,y,&r); /* assume that r is last on stack */    GEN r, q = dvmdii(x,y,&r); /* assume that r is last on stack */
   GEN *gptr[2];    GEN *gptr[2];
   
Line 1653  truedvmdii(GEN x, GEN y, GEN *z)
Line 1704  truedvmdii(GEN x, GEN y, GEN *z)
   }    }
   
   if (z == ONLY_REM)    if (z == ONLY_REM)
   {    { /* r += sign(y) * y, that is |y| */
     r = subiispec(y+2,r+2, lgefint(y)-2,lgefint(r)-2);      r = subiispec(y+2,r+2, lgefint(y)-2,lgefint(r)-2);
     return gerepileuptoint(av, r);      return gerepileuptoint(av, r);
   }    }
Line 1661  truedvmdii(GEN x, GEN y, GEN *z)
Line 1712  truedvmdii(GEN x, GEN y, GEN *z)
   if (!z) return gerepileuptoint(av, q);    if (!z) return gerepileuptoint(av, q);
   
   *z = subiispec(y+2,r+2, lgefint(y)-2,lgefint(r)-2);    *z = subiispec(y+2,r+2, lgefint(y)-2,lgefint(r)-2);
   gptr[0]=&q; gptr[1]=z; gerepilemanysp(av,(long)r,gptr,2);    gptr[0]=&q; gptr[1]=z; gerepilemanysp(av,(gpmem_t)r,gptr,2);
   return q;    return q;
 }  }
   
   /* Montgomery reduction.
    * N has k words, assume T >= 0 has less than 2k.
    * Return res := T / B^k mod N, where B = 2^BIL
    * such that 0 <= res < T/B^k + N  and  res has less than k words */
   GEN
   red_montgomery(GEN T, GEN N, ulong inv)
   {
     gpmem_t av;
     GEN Te,Td,Ne,Nd, scratch;
     ulong m,t,d,k = lgefint(N)-2;
     int carry;
     long i,j;
     LOCAL_HIREMAINDER;
     LOCAL_OVERFLOW;
   
     if (k == 0) return gzero;
     d = lgefint(T)-2; /* <= 2*k */
   #ifdef DEBUG
     if (d > 2*k) err(bugparier,"red_montgomery");
   #endif
     if (k == 1)
     { /* as below, special cased for efficiency */
       ulong n = (ulong)N[2];
       t = (ulong)T[d+1];
       m = t * inv;
       (void)addll(mulll(m, n), t); /* = 0 */
       t = hiremainder + overflow;
       if (d == 2)
       {
         t = addll((ulong)T[2], t);
         if (overflow) t -= n; /* t > n doesn't fit in 1 word */
       }
       return utoi(t);
     }
     /* assume k >= 2 */
     av = avma; scratch = new_chunk(k<<1); /* >= k + 2: result fits */
   
     /* copy T to scratch space (pad with zeroes to 2k words) */
     Td = (GEN)av;
     Te = T + (d+2);
     for (i=0; i < d     ; i++) *--Td = *--Te;
     for (   ; i < (k<<1); i++) *--Td = 0;
   
     Te = (GEN)av; /* 1 beyond end of T mantissa */
     Ne = N + k+2; /* 1 beyond end of N mantissa */
   
     carry = 0;
     for (i=0; i<k; i++) /* set T := T/B nod N, k times */
     {
       Td = Te; /* one beyond end of (new) T mantissa */
       Nd = Ne;
       m = *--Td * inv; /* solve T + m N = O(B) */
   
       /* set T := (T + mN) / B */
       Te = Td;
       (void)addll(mulll(m, *--Nd), *Td); /* = 0 */
       for (j=1; j<k; j++)
       {
         hiremainder += overflow;
         t = addll(addmul(m, *--Nd), *--Td); *Td = t;
       }
       overflow += hiremainder;
       t = addll(overflow, *--Td); *Td = t + carry;
       carry = (overflow || (carry && *Td == 0));
     }
     if (carry)
     { /* Td > N overflows (k+1 words), set Td := Td - N */
       Td = Te;
       Nd = Ne;
       t = subll(*--Td, *--Nd); *Td = t;
       while (Td > scratch) { t = subllx(*--Td, *--Nd); *Td = t; }
     }
   
     /* copy result */
     Td = (GEN)av;
     while (! *scratch) scratch++; /* strip leading zeroes */
     while (Te > scratch) *--Td = *--Te;
     k = ((GEN)av - Td) + 2;
     *--Td = evalsigne(1) | evallgefint(k);
     *--Td = evaltyp(t_INT) | evallg(k);
   #ifdef DEBUG
   {
     long l = lgefint(N)-2, s = BITS_IN_LONG*l;
     GEN R = shifti(gun, s);
     GEN res = resii(mulii(T, mpinvmod(R, N)), N);
     if (k > lgefint(N)
       || !egalii(resii(Td,N),res)
       || cmpii(Td, addii(shifti(T, -s), N)) >= 0) err(bugparier,"red_montgomery");
   }
   #endif
     avma = (gpmem_t)Td; return Td;
   }
   
 /* EXACT INTEGER DIVISION */  /* EXACT INTEGER DIVISION */
   
 /* Find c such that 1=c*b mod 2^BITS_IN_LONG, assuming b odd (unchecked) */  /* Find c such that 1=c*b mod 2^BITS_IN_LONG, assuming b odd (unchecked) */
 static ulong  ulong
 invrev(ulong b)  invrev(ulong b)
 {  {
   ulong x;    ulong x;
Line 1688  invrev(ulong b)
Line 1832  invrev(ulong b)
 }  }
   
 /* assume xy>0, y odd */  /* assume xy>0, y odd */
 GEN  GEN
 diviuexact(GEN x, ulong y)  diviuexact(GEN x, ulong y)
 {  {
   long i,lz,lx;    long i,lz,lx;
Line 1697  diviuexact(GEN x, ulong y)
Line 1841  diviuexact(GEN x, ulong y)
   
   if (y == 1) return icopy(x);    if (y == 1) return icopy(x);
   lx = lgefint(x);    lx = lgefint(x);
   if (lx == 3) return stoi((ulong)x[2] / y);    if (lx == 3) return utoi((ulong)x[2] / y);
   yinv = invrev(y);    yinv = invrev(y);
   lz = (y <= (ulong)x[2]) ? lx : lx-1;    lz = (y <= (ulong)x[2]) ? lx : lx-1;
   z = new_chunk(lz);    z = new_chunk(lz);
Line 1729  diviuexact(GEN x, ulong y)
Line 1873  diviuexact(GEN x, ulong y)
   z += i-2; lz -= i-2;    z += i-2; lz -= i-2;
   z[0] = evaltyp(t_INT)|evallg(lz);    z[0] = evaltyp(t_INT)|evallg(lz);
   z[1] = evalsigne(1)|evallg(lz);    z[1] = evalsigne(1)|evallg(lz);
   avma = (ulong)z; return z;    avma = (gpmem_t)z; return z;
 }  }
   
 /* Find z such that x=y*z, knowing that y | x (unchecked)  /* Find z such that x=y*z, knowing that y | x (unchecked)
Line 1738  diviuexact(GEN x, ulong y)
Line 1882  diviuexact(GEN x, ulong y)
 GEN  GEN
 diviiexact(GEN x, GEN y)  diviiexact(GEN x, GEN y)
 {  {
   long av,lx,ly,lz,vy,i,ii,sx = signe(x), sy = signe(y);    long lx, ly, lz, vy, i, ii, sx = signe(x), sy = signe(y);
     gpmem_t av;
   ulong y0inv,q;    ulong y0inv,q;
   GEN z;    GEN z;
   
Line 1758  diviiexact(GEN x, GEN y)
Line 1903  diviiexact(GEN x, GEN y)
   avma = av; /* will erase our x,y when exiting */    avma = av; /* will erase our x,y when exiting */
   /* now y is odd */    /* now y is odd */
   ly = lgefint(y);    ly = lgefint(y);
   if (ly == 3)    if (ly == 3)
   {    {
     x = diviuexact(x,(ulong)y[2]);      x = diviuexact(x,(ulong)y[2]);
     if (signe(x)) setsigne(x,sx*sy); /* should have x != 0 at this point */      if (signe(x)) setsigne(x,sx*sy); /* should have x != 0 at this point */
Line 1783  diviiexact(GEN x, GEN y)
Line 1928  diviiexact(GEN x, GEN y)
     /* x := x - q * y */      /* x := x - q * y */
     (void)mulll(q,y[0]); limj = max(lx - lz, ii+3-ly);      (void)mulll(q,y[0]); limj = max(lx - lz, ii+3-ly);
     { /* update neither lowest word (could set it to 0) nor highest ones */      { /* update neither lowest word (could set it to 0) nor highest ones */
       register GEN x0 = x + (ii - 1), y0 = y - 1, xlim = x + limj;        register GEN x0 = x + (ii - 1), y0 = y - 1, xlim = x + limj;
       for (; x0 >= xlim; x0--, y0--)        for (; x0 >= xlim; x0--, y0--)
       {        {
         *x0 = subll(*x0, addmul(q,*y0));          *x0 = subll(*x0, addmul(q,*y0));
Line 1810  diviiexact(GEN x, GEN y)
Line 1955  diviiexact(GEN x, GEN y)
   z += i-2; lz -= (i-2);    z += i-2; lz -= (i-2);
   z[0] = evaltyp(t_INT)|evallg(lz);    z[0] = evaltyp(t_INT)|evallg(lz);
   z[1] = evalsigne(sx*sy)|evallg(lz);    z[1] = evalsigne(sx*sy)|evallg(lz);
   avma = (ulong)z; return z;    avma = (gpmem_t)z; return z;
 }  }
   
 long  long
Line 1884  absr_cmp(GEN x, GEN y)
Line 2029  absr_cmp(GEN x, GEN y)
 #define _sqri_l 47  #define _sqri_l 47
 #define _muli_l 25 /* optimal on PII 350MHz + gcc 2.7.2.1 (gp-dyn) */  #define _muli_l 25 /* optimal on PII 350MHz + gcc 2.7.2.1 (gp-dyn) */
   
 #if 0 /* for tunings */  #if 1 /* for tunings */
 long KARATSUBA_SQRI_LIMIT = _sqri_l;  long KARATSUBA_SQRI_LIMIT = _sqri_l;
 long KARATSUBA_MULI_LIMIT = _muli_l;  long KARATSUBA_MULI_LIMIT = _muli_l;
   
Line 1945  muliispec(GEN x, GEN y, long nx, long ny)
Line 2090  muliispec(GEN x, GEN y, long nx, long ny)
   if (*zd == 0) { zd++; lz--; } /* normalize */    if (*zd == 0) { zd++; lz--; } /* normalize */
   *--zd = evalsigne(1) | evallgefint(lz);    *--zd = evalsigne(1) | evallgefint(lz);
   *--zd = evaltyp(t_INT) | evallg(lz);    *--zd = evaltyp(t_INT) | evallg(lz);
   avma=(long)zd; return zd;    avma=(gpmem_t)zd; return zd;
 }  }
   
 #ifdef INLINE  #ifdef INLINE
Line 1967  sqrispec(GEN x, long nx)
Line 2112  sqrispec(GEN x, long nx)
     *--zd = hiremainder; goto END;      *--zd = hiremainder; goto END;
   }    }
   xd = x + nx;    xd = x + nx;
   
   /* compute double products --> zd */    /* compute double products --> zd */
   p1 = *--xd; yd = xd; --zd;    p1 = *--xd; yd = xd; --zd;
   *--zd = mulll(p1, *--yd); z2e = zd;    *--zd = mulll(p1, *--yd); z2e = zd;
Line 2009  END:
Line 2154  END:
   if (*zd == 0) { zd++; lz--; } /* normalize */    if (*zd == 0) { zd++; lz--; } /* normalize */
   *--zd = evalsigne(1) | evallgefint(lz);    *--zd = evalsigne(1) | evallgefint(lz);
   *--zd = evaltyp(t_INT) | evallg(lz);    *--zd = evaltyp(t_INT) | evallg(lz);
   avma=(long)zd; return zd;    avma=(gpmem_t)zd; return zd;
 }  }
   
 /* return (x shifted left d words) + y. Assume d > 0, x > 0 and y >= 0 */  /* return (x shifted left d words) + y. Assume d > 0, x > 0 and y >= 0 */
Line 2048  static GEN
Line 2193  static GEN
 quickmulii(GEN a, GEN b, long na, long nb)  quickmulii(GEN a, GEN b, long na, long nb)
 {  {
   GEN a0,c,c0;    GEN a0,c,c0;
   long av,n0,n0a,i;    long n0, n0a, i;
     gpmem_t av;
   
   if (na < nb) swapspec(a,b, na,nb);    if (na < nb) swapspec(a,b, na,nb);
   if (nb == 1) return mulsispec(*b, a, na);    if (nb == 1) return mulsispec(*b, a, na);
Line 2111  GEN
Line 2257  GEN
 resiimul(GEN x, GEN sy)  resiimul(GEN x, GEN sy)
 {  {
   GEN r, q, y = (GEN)sy[1], invy;    GEN r, q, y = (GEN)sy[1], invy;
   long av = avma, k;    long k;
     gpmem_t av = avma;
   
   k = cmpii(x, y);    k = cmpii(x, y);
   if (k <= 0) return k? icopy(x): gzero;    if (k <= 0) return k? icopy(x): gzero;
Line 2140  resmod2n(GEN x, long n)
Line 2287  resmod2n(GEN x, long n)
 {  {
   long hi,l,k,lx,ly;    long hi,l,k,lx,ly;
   GEN z, xd, zd;    GEN z, xd, zd;
   
   if (!signe(x) || !n) return gzero;    if (!signe(x) || !n) return gzero;
   
   l = n & (BITS_IN_LONG-1);    /* n % BITS_IN_LONG */    l = n & (BITS_IN_LONG-1);    /* n % BITS_IN_LONG */
   k = n >> TWOPOTBITS_IN_LONG; /* n / BITS_IN_LONG */    k = n >> TWOPOTBITS_IN_LONG; /* n / BITS_IN_LONG */
   lx = lgefint(x);    lx = lgefint(x);
   if (lx < k+3) return icopy(x);    if (lx < k+3) return icopy(x);
Line 2172  static GEN
Line 2319  static GEN
 quicksqri(GEN a, long na)  quicksqri(GEN a, long na)
 {  {
   GEN a0,c;    GEN a0,c;
   long av,n0,n0a,i;    long n0, n0a, i;
     gpmem_t av;
   
   if (na < KARATSUBA_SQRI_LIMIT) return sqrispec(a,na);    if (na < KARATSUBA_SQRI_LIMIT) return sqrispec(a,na);
   i=(na>>1); n0=na-i; na=i;    i=(na>>1); n0=na-i; na=i;
Line 2216  karamulrr1(GEN x, GEN y)
Line 2364  karamulrr1(GEN x, GEN y)
   /* ensure that lg(y) >= lg(x) */    /* ensure that lg(y) >= lg(x) */
   if (lx>ly) { lx=ly; z=x; x=y; y=z; flag=1; } else flag = (lx!=ly);    if (lx>ly) { lx=ly; z=x; x=y; y=z; flag=1; } else flag = (lx!=ly);
   if (lx < MULRR_LIMIT) return mulrr(x,y);    if (lx < MULRR_LIMIT) return mulrr(x,y);
   sx=signe(x); sy=signe(y);    sx=signe(x); sy=signe(y); e = expo(x)+expo(y);
   e = evalexpo(expo(x)+expo(y));    if (!sx || !sy) return realzero_bit(e);
   if (!sx || !sy) { z=cgetr(3); z[2]=0; z[1]=e; return z; }  
   if (sy<0) sx = -sx;    if (sy<0) sx = -sx;
   ly=lx+flag; z=cgetr(lx);    ly=lx+flag; z=cgetr(lx);
   lz2 = (lx>>1); lz3 = lx-lz2;    lz2 = (lx>>1); lz3 = lx-lz2;
Line 2246  karamulrr1(GEN x, GEN y)
Line 2393  karamulrr1(GEN x, GEN y)
     garde = (hi[lx+2] << 1);      garde = (hi[lx+2] << 1);
     shift_left(z,hi,2,lx+1, garde, 1);      shift_left(z,hi,2,lx+1, garde, 1);
   }    }
   z[1]=evalsigne(sx) | e;    z[1]=evalsigne(sx) | evalexpo(e);
   if (garde < 0)    if (garde < 0)
   { /* round to nearest */    { /* round to nearest */
     i=lx+2; do z[--i]++; while (z[i]==0);      i=lx+2; do z[--i]++; while (z[i]==0);
     if (i==1) z[2]=HIGHBIT;      if (i==1) z[2]=HIGHBIT;
   }    }
   avma=(long)z; return z;    avma=(gpmem_t)z; return z;
 }  }
   
 GEN  GEN
Line 2264  karamulrr2(GEN x, GEN y)
Line 2411  karamulrr2(GEN x, GEN y)
   if (lx>ly) { lx=ly; z=x; x=y; y=z; flag=1; } else flag = (lx!=ly);    if (lx>ly) { lx=ly; z=x; x=y; y=z; flag=1; } else flag = (lx!=ly);
   if (lx < MULRR2_LIMIT) return mulrr(x,y);    if (lx < MULRR2_LIMIT) return mulrr(x,y);
   ly=lx+flag; sx=signe(x); sy=signe(y);    ly=lx+flag; sx=signe(x); sy=signe(y);
   e = evalexpo(expo(x)+expo(y));    e = expo(x)+expo(y);
   if (!sx || !sy) { z=cgetr(3); z[2]=0; z[1]=e; return z; }    if (!sx || !sy) return realzero_bit(e);
   if (sy<0) sx = -sx;    if (sy<0) sx = -sx;
   z=cgetr(lx);    z=cgetr(lx);
   hi=quickmulii(y+2,x+2,ly-2,lx-2);    hi=quickmulii(y+2,x+2,ly-2,lx-2);
Line 2279  karamulrr2(GEN x, GEN y)
Line 2426  karamulrr2(GEN x, GEN y)
     garde = (hi[lx] << 1);      garde = (hi[lx] << 1);
     shift_left(z,hi,2,lx-1, garde, 1);      shift_left(z,hi,2,lx-1, garde, 1);
   }    }
   z[1]=evalsigne(sx) | e;    z[1]=evalsigne(sx) | evalexpo(e);
   if (garde < 0)    if (garde < 0)
   { /* round to nearest */    { /* round to nearest */
     i=lx; do z[--i]++; while (z[i]==0);      i=lx; do z[--i]++; while (z[i]==0);
     if (i==1) z[2]=HIGHBIT;      if (i==1) z[2]=HIGHBIT;
   }    }
   avma=(long)z; return z;    avma=(gpmem_t)z; return z;
 }  }
   
 GEN  GEN
Line 2304  GEN
Line 2451  GEN
 karamulir(GEN x, GEN y, long flag)  karamulir(GEN x, GEN y, long flag)
 {  {
   long sx=signe(x),lz,i;    long sx=signe(x),lz,i;
   GEN z,temp,z1;    GEN z, z1;
   
   if (!sx) return gzero;    if (!sx) return gzero;
   lz=lg(y); z=cgetr(lz);    lz = lg(y); z = cgetr(lz);
   temp=cgetr(lz+1); affir(x,temp);    z1 = karamulrr(itor(x, lz+1), y, flag);
   z1=karamulrr(temp,y,flag);  
   for (i=1; i<lz; i++) z[i]=z1[i];    for (i=1; i<lz; i++) z[i]=z1[i];
   avma=(long)z; return z;    avma=(gpmem_t)z; return z;
 }  }
 #endif  #endif
   
 #ifdef LONG_IS_64BIT  #ifdef LONG_IS_64BIT
   long
   expodb(double x)
   {
     union { double f; ulong i; } fi;
     const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */
     const int exp_mid = 0x3ff;/* exponent bias */
   
 #if PARI_BYTE_ORDER == LITTLE_ENDIAN_64 || PARI_BYTE_ORDER == BIG_ENDIAN_64    if (x==0.) return -1023;
 #else    fi.f = x;
   error... unknown machine    return ((fi.i & (HIGHBIT-1)) >> mant_len) - exp_mid;
 #endif  }
   
 GEN  GEN
 dbltor(double x)  dbltor(double x)
 {  {
   GEN z = cgetr(3);    GEN z;
   long e;    long e;
   union { double f; ulong i; } fi;    union { double f; ulong i; } fi;
   const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */    const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */
Line 2333  dbltor(double x)
Line 2485  dbltor(double x)
   const int expo_len = 11; /* number of bits of exponent */    const int expo_len = 11; /* number of bits of exponent */
   LOCAL_HIREMAINDER;    LOCAL_HIREMAINDER;
   
   if (x==0) { z[1]=evalexpo(-308); z[2]=0; return z; }    if (x==0.) return realzero_bit(-1023);
   fi.f = x;    fi.f = x; z = cgetr(DEFAULTPREC);
   e = ((fi.i & (HIGHBIT-1)) >> mant_len) - exp_mid;    e = ((fi.i & (HIGHBIT-1)) >> mant_len) - exp_mid;
   z[1] = evalexpo(e) | evalsigne(x<0? -1: 1);    z[1] = evalexpo(e) | evalsigne(x<0? -1: 1);
   z[2] = (fi.i << expo_len) | HIGHBIT;    z[2] = (fi.i << expo_len) | HIGHBIT;
Line 2365  rtodbl(GEN x)
Line 2517  rtodbl(GEN x)
   return fi.f;    return fi.f;
 }  }
   
 #else  #else /* LONG_IS_64BIT */
   
 #if   PARI_BYTE_ORDER == LITTLE_ENDIAN  #if   PARI_DOUBLE_FORMAT == 1
 #  define INDEX0 1  #  define INDEX0 1
 #  define INDEX1 0  #  define INDEX1 0
 #elif PARI_BYTE_ORDER == BIG_ENDIAN  #elif PARI_DOUBLE_FORMAT == 0
 #  define INDEX0 0  #  define INDEX0 0
 #  define INDEX1 1  #  define INDEX1 1
 #else  
    error... unknown machine  
 #endif  #endif
   
   long
   expodb(double x)
   {
     union { double f; ulong i[2]; } fi;
     const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */
     const int exp_mid = 0x3ff;/* exponent bias */
     const int shift = mant_len-32;
   
     if (x==0.) return -1023;
     fi.f = x;
     {
       const ulong a = fi.i[INDEX0];
       return ((a & (HIGHBIT-1)) >> shift) - exp_mid;
     }
   }
   
 GEN  GEN
 dbltor(double x)  dbltor(double x)
 {  {
Line 2385  dbltor(double x)
Line 2551  dbltor(double x)
   union { double f; ulong i[2]; } fi;    union { double f; ulong i[2]; } fi;
   const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */    const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */
   const int exp_mid = 0x3ff;/* exponent bias */    const int exp_mid = 0x3ff;/* exponent bias */
   const int shift = mant_len-32;  
   const int expo_len = 11; /* number of bits of exponent */    const int expo_len = 11; /* number of bits of exponent */
     const int shift = mant_len-32;
   
   if (x==0) { z=cgetr(3); z[1]=evalexpo(-308); z[2]=0; return z; }    if (x==0.) return realzero_bit(-1023);
   fi.f = x; z=cgetr(4);    fi.f = x; z=cgetr(DEFAULTPREC);
   {    {
     const ulong a = fi.i[INDEX0];      const ulong a = fi.i[INDEX0];
     const ulong b = fi.i[INDEX1];      const ulong b = fi.i[INDEX1];
Line 2409  rtodbl(GEN x)
Line 2575  rtodbl(GEN x)
   union { double f; ulong i[2]; } fi;    union { double f; ulong i[2]; } fi;
   const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */    const int mant_len = 52;  /* mantissa bits (excl. hidden bit) */
   const int exp_mid = 0x3ff;/* exponent bias */    const int exp_mid = 0x3ff;/* exponent bias */
   const int shift = mant_len-32;  
   const int expo_len = 11; /* number of bits of exponent */    const int expo_len = 11; /* number of bits of exponent */
     const int shift = mant_len-32;
   
   if (typ(x)==t_INT && !s) return 0.0;    if (typ(x)==t_INT && !s) return 0.0;
   if (typ(x)!=t_REAL) err(typeer,"rtodbl");    if (typ(x)!=t_REAL) err(typeer,"rtodbl");
Line 2432  rtodbl(GEN x)
Line 2598  rtodbl(GEN x)
   fi.i[INDEX1] = (a << (BITS_IN_LONG-expo_len)) | (b >> expo_len);    fi.i[INDEX1] = (a << (BITS_IN_LONG-expo_len)) | (b >> expo_len);
   return fi.f;    return fi.f;
 }  }
 #endif  #endif /* LONG_IS_64BIT */
   
 /* Old cgiv without reference count (which was not used anyway)  /* Old cgiv without reference count (which was not used anyway)
  * Should be a macro.   * Should be a macro.
Line 2441  void
Line 2607  void
 cgiv(GEN x)  cgiv(GEN x)
 {  {
   if (x == (GEN) avma)    if (x == (GEN) avma)
     avma = (long) (x+lg(x));      avma = (gpmem_t) (x+lg(x));
 }  }
   
 /********************************************************************/  /********************************************************************/
Line 2526  cgiv(GEN x)
Line 2692  cgiv(GEN x)
  *   *
  * Note that these routines do not know and do not need to know about the   * Note that these routines do not know and do not need to know about the
  * PARI stack.   * PARI stack.
  *   *
  * Added 2001Jan15:   * Added 2001Jan15:
  * rgcduu() is a variant of xxgcduu() which does not have f  (the effect is   * rgcduu() is a variant of xxgcduu() which does not have f  (the effect is
  * that of f=0),  but instead has a ulong vmax parameter, for use in rational   * that of f=0),  but instead has a ulong vmax parameter, for use in rational
Line 2747  rgcduu(ulong d, ulong d1, ulong vmax,
Line 2913  rgcduu(ulong d, ulong d1, ulong vmax,
  * as many partial quotients of d/d1 as possible, and assigns to [u,u1;v,v1]   * as many partial quotients of d/d1 as possible, and assigns to [u,u1;v,v1]
  * the product of all the [0, 1; 1, q_j] thus obtained, where the leftmost   * the product of all the [0, 1; 1, q_j] thus obtained, where the leftmost
  * factor arises from the quotient of the first division step.   * factor arises from the quotient of the first division step.
  *   *
  * For use in rational reconstruction, input param vmax can be given a   * For use in rational reconstruction, input param vmax can be given a
  * nonzero value.  In this case, we will return early as soon as v1 > vmax   * nonzero value.  In this case, we will return early as soon as v1 > vmax
  * (i.e. it is guaranteed that v <= vmax). --2001Jan15   * (i.e. it is guaranteed that v <= vmax). --2001Jan15
Line 3282  lgcdii(ulong* d, ulong* d1,
Line 3448  lgcdii(ulong* d, ulong* d1,
       *u = xu1; *u1 = xu; *v = xv1; *v1 = xv;        *u = xu1; *u1 = xu; *v = xv1; *v1 = xv;
       break;        break;
     }      }
   
     res++; /* commit dd, xu1, xv1 */      res++; /* commit dd, xu1, xv1 */
     dd = tmpd; xu1 = tmpu; xv1 = tmpv;      dd = tmpd; xu1 = tmpu; xv1 = tmpv;
 #ifdef DEBUG_LEHMER  #ifdef DEBUG_LEHMER
Line 3303  lgcdii(ulong* d, ulong* d1,
Line 3469  lgcdii(ulong* d, ulong* d1,
  *==================================   *==================================
  *    If a is invertible, return 1, and set res  = a^{ -1 }   *    If a is invertible, return 1, and set res  = a^{ -1 }
  *    Otherwise, return 0, and set res = gcd(a,b)   *    Otherwise, return 0, and set res = gcd(a,b)
  *   *
  * This is sufficiently different from bezout() to be implemented separately   * This is sufficiently different from bezout() to be implemented separately
  * instead of having a bunch of extra conditionals in a single function body   * instead of having a bunch of extra conditionals in a single function body
  * to meet both purposes.   * to meet both purposes.
Line 3313  int
Line 3479  int
 invmod(GEN a, GEN b, GEN *res)  invmod(GEN a, GEN b, GEN *res)
 {  {
   GEN v,v1,d,d1,q,r;    GEN v,v1,d,d1,q,r;
   long av,av1,lim;    gpmem_t av, av1, lim;
   long s;    long s;
   ulong g;    ulong g;
   ulong xu,xu1,xv,xv1;          /* Lehmer stage recurrence matrix */    ulong xu,xu1,xv,xv1;          /* Lehmer stage recurrence matrix */
Line 3481  bezout(GEN a, GEN b, GEN *pu, GEN *pv)
Line 3647  bezout(GEN a, GEN b, GEN *pu, GEN *pv)
 {  {
   GEN t,u,u1,v,v1,d,d1,q,r;    GEN t,u,u1,v,v1,d,d1,q,r;
   GEN *pt;    GEN *pt;
   long av,av1,lim;    gpmem_t av, av1, lim;
   long s;    long s;
   int sa, sb;    int sa, sb;
   ulong g;    ulong g;
Line 3802  cbezout(long a,long b,long *uu,long *vv)
Line 3968  cbezout(long a,long b,long *uu,long *vv)
      AUTHOR = {Collins, George E. and Encarnaci{\'o}n, Mark J.},       AUTHOR = {Collins, George E. and Encarnaci{\'o}n, Mark J.},
       TITLE = {Efficient rational number reconstruction},        TITLE = {Efficient rational number reconstruction},
     JOURNAL = {J. Symbolic Comput.},      JOURNAL = {J. Symbolic Comput.},
    FJOURNAL = {Journal of Symbolic Computation},  
      VOLUME = {20},       VOLUME = {20},
        YEAR = {1995},         YEAR = {1995},
      NUMBER = {3},       NUMBER = {3},
       PAGES = {287--297},        PAGES = {287--297},
        ISSN = {0747-7171},  
     MRCLASS = {11Y16 (68Q40)},  
    MRNUMBER = {97c:11116},  
  MRREVIEWER = {Maurice Mignotte},  
  }   }
  *    Preprint available from:   * Preprint available from:
  *    ftp://ftp.risc.uni-linz.ac.at/pub/techreports/1994/94-64.ps.gz   * ftp://ftp.risc.uni-linz.ac.at/pub/techreports/1994/94-64.ps.gz
  */   */
   
 /* #define DEBUG_RATLIFT */  /* #define DEBUG_RATLIFT */
Line 3822  int
Line 3983  int
 ratlift(GEN x, GEN m, GEN *a, GEN *b, GEN amax, GEN bmax)  ratlift(GEN x, GEN m, GEN *a, GEN *b, GEN amax, GEN bmax)
 {  {
   GEN d,d1,v,v1,q,r;    GEN d,d1,v,v1,q,r;
   ulong av = avma, av1, lim;    gpmem_t av = avma, av1, lim;
   long lb,lr,lbb,lbr,s,s0;    long lb,lr,lbb,lbr,s,s0;
   ulong vmax;    ulong vmax;
   ulong xu,xu1,xv,xv1;          /* Lehmer stage recurrence matrix */    ulong xu,xu1,xv,xv1;          /* Lehmer stage recurrence matrix */
Line 3904  ratlift(GEN x, GEN m, GEN *a, GEN *b, GEN amax, GEN bm
Line 4065  ratlift(GEN x, GEN m, GEN *a, GEN *b, GEN amax, GEN bm
    * Furthermore, at the start of the loop body we have in fact     * Furthermore, at the start of the loop body we have in fact
    * (c') 0 <= v < v1 <= bmax, d > d1 > amax >= 0,     * (c') 0 <= v < v1 <= bmax, d > d1 > amax >= 0,
    * (and are never done already).     * (and are never done already).
    *     *
    * Main loop is similar to those of invmod() and bezout(), except for     * Main loop is similar to those of invmod() and bezout(), except for
    * having to determine appropriate vmax bounds, and checking termination     * having to determine appropriate vmax bounds, and checking termination
    * conditions.  The signe(d1) condition is only for paranoia     * conditions.  The signe(d1) condition is only for paranoia
Line 4209  mpsqrtl(GEN a)
Line 4370  mpsqrtl(GEN a)
   long l = lgefint(a);    long l = lgefint(a);
   ulong x,y,z,k,m;    ulong x,y,z,k,m;
   int L, s;    int L, s;
   
   if (l <= 3) return l==2? 0: usqrtsafe((ulong)a[2]);    if (l <= 3) return l==2? 0: usqrtsafe((ulong)a[2]);
   s = bfffo(a[2]);    s = bfffo(a[2]);
   if (s > 1)    if (s > 1)
Line 4242  mpsqrtl(GEN a)
Line 4403  mpsqrtl(GEN a)
   }    }
   while (x < y);    while (x < y);
   return y;    return y;
   }
   
   /* target should point to a buffer of source_end - source + 1 ulongs.
   
      fills this buffer by bits of ulongs in source..source_end-1 shifted
      right sh units; the "most significant" sh bits of the result are
      set to be the least significant sh bits of prepend.
   
      The ordering of bits in this bitmap is the same as for t_INT.
   
      sh should not exceed BITS_IN_LONG.
    */
   void
   shift_r(ulong *target, ulong *source, ulong *source_end, ulong prepend, ulong sh)
   {
     if (sh) {                             /* shift_words_r() works */
       register ulong sh_complement = BITS_IN_LONG - sh;
   
       shift_words_r(target, source, source_end, prepend, sh, sh_complement);
     } else {
       int i;
   
       for (i=0; i < source_end - source; i++)
         target[i] = source[i];
     }
 }  }
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>