[BACK]Return to arith2.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / pari-2.2 / src / basemath

Diff for /OpenXM_contrib/pari-2.2/src/basemath/Attic/arith2.c between version 1.1 and 1.2

version 1.1, 2001/10/02 11:17:01 version 1.2, 2002/09/11 07:26:48
Line 25  extern GEN arith_proto(long f(GEN), GEN x, int do_erro
Line 25  extern GEN arith_proto(long f(GEN), GEN x, int do_erro
 extern GEN arith_proto2(long f(GEN,GEN), GEN x, GEN n);  extern GEN arith_proto2(long f(GEN,GEN), GEN x, GEN n);
 extern GEN garith_proto(GEN f(GEN), GEN x, int do_error);  extern GEN garith_proto(GEN f(GEN), GEN x, int do_error);
 extern GEN garith_proto2gs(GEN f(GEN,long), GEN x, long y);  extern GEN garith_proto2gs(GEN f(GEN,long), GEN x, long y);
   extern GEN garith_proto3ggs(GEN f(GEN,GEN,long), GEN x, GEN y, long z);
   
 #define sqru(i) (muluu((i),(i)))  #define sqru(i) (muluu((i),(i)))
   
Line 38  GEN
Line 39  GEN
 prime(long n)  prime(long n)
 {  {
   byteptr p = diffptr;    byteptr p = diffptr;
   long c, prime = 0;    long prime = 0;
   
   if (n <= 0) err(talker, "n-th prime meaningless if n = %ld",n);    if (n <= 0) err(talker, "n-th prime meaningless if n = %ld",n);
   while (n--)    while (n--)
   {    {
     c = *p++; if (!c) err(primer1);      NEXT_PRIME_VIADIFF_CHECK(prime,p);
     prime += c;  
   }    }
   return stoi(prime);    return stoi(prime);
 }  }
Line 56  pith(long n)
Line 56  pith(long n)
   ulong prime = 0, res = 0;    ulong prime = 0, res = 0;
   
   if (n <= 0) err(talker, "pith meaningless if n = %ld",n);    if (n <= 0) err(talker, "pith meaningless if n = %ld",n);
   if (maxprime() <= n) err(primer1);    if (maxprime() <= (ulong)n) err(primer1);
   while (prime<=n) { prime += *p++; res++; }    while (prime<=(ulong)n) {
         NEXT_PRIME_VIADIFF(prime,p);
         res++;
     }
   return stoi(res-1);    return stoi(res-1);
 }  }
   
Line 65  GEN
Line 68  GEN
 primes(long n)  primes(long n)
 {  {
   byteptr p = diffptr;    byteptr p = diffptr;
   long c, prime = 0;    long prime = 0;
   GEN y,z;    GEN y,z;
   
   if (n < 0) n = 0;    if (n < 0) n = 0;
   z = y = cgetg(n+1,t_VEC);    z = y = cgetg(n+1,t_VEC);
   while (n--)    while (n--)
   {    {
     c = *p++; if (!c) err(primer1);      NEXT_PRIME_VIADIFF_CHECK(prime,p);
     prime += c; *++z = lstoi(prime);      *++z = lstoi(prime);
   }    }
   return y;    return y;
 }  }
Line 85  primes(long n)
Line 88  primes(long n)
  * when maxnum (size) is moderate.   * when maxnum (size) is moderate.
  */   */
 static byteptr  static byteptr
 initprimes1(long size, long *lenp, long *lastp)  initprimes1(ulong size, long *lenp, long *lastp)
 {  {
   long k;    long k;
   byteptr q,r,s,fin, p = (byteptr) gpmalloc(size+2);    byteptr q,r,s,fin, p = (byteptr) gpmalloc(size+2);
Line 110  initprimes1(long size, long *lenp, long *lastp)
Line 113  initprimes1(long size, long *lenp, long *lastp)
   fprintferr("initprimes1: q = %ld, len = %ld, last = %ld\n",    fprintferr("initprimes1: q = %ld, len = %ld, last = %ld\n",
              q, *lenp, *lastp);               q, *lenp, *lastp);
 #endif  #endif
   return (byteptr) gprealloc(p,r-p,size + 2);    return (byteptr) gprealloc(p,r-p);
 }  }
   
 #define PRIME_ARENA (200 * 1024) /* No slowdown even with 256K level-2 cache */  #define PRIME_ARENA (200 * 1024) /* No slowdown even with 256K level-2 cache */
Line 119  initprimes1(long size, long *lenp, long *lastp)
Line 122  initprimes1(long size, long *lenp, long *lastp)
    recursive call will bottom out and invoke initprimes1() at once.     recursive call will bottom out and invoke initprimes1() at once.
    (Not static;  might conceivably be useful to someone in library mode) */     (Not static;  might conceivably be useful to someone in library mode) */
 byteptr  byteptr
 initprimes0(ulong maxnum, long *lenp, long *lastp)  initprimes0(ulong maxnum, long *lenp, ulong *lastp)
 {  {
   long k, size, alloced, asize, psize, rootnum, curlow, last, remains;    long k, size, alloced, asize, psize, rootnum;
   byteptr q,r,s,fin, p, p1, fin1, plast, curdiff;    byteptr q,r,s,fin, p, p1, fin1, plast, curdiff;
     ulong last, remains, curlow;
   
 #if 0  #if 0
   fprintferr("initprimes0: called for maxnum = %lu\n", maxnum);    fprintferr("initprimes0: called for maxnum = %lu\n", maxnum);
Line 160  initprimes0(ulong maxnum, long *lenp, long *lastp)
Line 164  initprimes0(ulong maxnum, long *lenp, long *lastp)
   asize = ARENA_IN_ROOTS * rootnum;     /* Make % overhead negligeable. */    asize = ARENA_IN_ROOTS * rootnum;     /* Make % overhead negligeable. */
   if (asize < PRIME_ARENA) asize = PRIME_ARENA;    if (asize < PRIME_ARENA) asize = PRIME_ARENA;
   if (asize > remains) asize = remains + 1;/* + room for a sentinel byte */    if (asize > remains) asize = remains + 1;/* + room for a sentinel byte */
   alloced = (avma - bot < asize>>1); /* enough room on the stack ? */    alloced = (avma - bot < (size_t)asize>>1); /* enough room on the stack ? */
   if (alloced)    if (alloced)
     p = (byteptr) gpmalloc(asize);      p = (byteptr) gpmalloc(asize);
   else    else
Line 182  initprimes0(ulong maxnum, long *lenp, long *lastp)
Line 186  initprimes0(ulong maxnum, long *lenp, long *lastp)
     }      }
     memset(p, 0, asize);      memset(p, 0, asize);
     /* p corresponds to curlow.  q runs over primediffs.  */      /* p corresponds to curlow.  q runs over primediffs.  */
       /* Don't care about DIFFPTR_SKIP: false positives provide no problem */
     for (q = p1+2, k = 3; q <= fin1; k += *q++)      for (q = p1+2, k = 3; q <= fin1; k += *q++)
     {      {
       /* The first odd number which is >= curlow and divisible by p        /* The first odd number which is >= curlow and divisible by p
Line 207  initprimes0(ulong maxnum, long *lenp, long *lastp)
Line 212  initprimes0(ulong maxnum, long *lenp, long *lastp)
     /* now q runs over addresses corresponding to primes */      /* now q runs over addresses corresponding to primes */
     for (q = p; ; plast = q++)      for (q = p; ; plast = q++)
     {      {
         long d;
   
       while (*q) q++;           /* use 0-sentinel at end */        while (*q) q++;           /* use 0-sentinel at end */
       if (q >= fin) break;        if (q >= fin) break;
       *curdiff++ = (q - plast) << 1;        d = (q - plast) << 1;
         while (d >= DIFFPTR_SKIP)
           *curdiff++ = DIFFPTR_SKIP, d -= DIFFPTR_SKIP;
         *curdiff++ = d;
     }      }
     plast -= asize - 1;      plast -= asize - 1;
     remains -= asize - 1;      remains -= asize - 1;
Line 220  initprimes0(ulong maxnum, long *lenp, long *lastp)
Line 230  initprimes0(ulong maxnum, long *lenp, long *lastp)
   *lenp = curdiff - p1;    *lenp = curdiff - p1;
   *lastp = last;    *lastp = last;
   if (alloced) free(p);    if (alloced) free(p);
   return (byteptr) gprealloc(p1, *lenp, size);    return (byteptr) gprealloc(p1, *lenp);
 }  }
 #if 0 /* not yet... GN */  #if 0 /* not yet... GN */
 /* The diffptr table will contain at least 6548 entries (up to and including  /* The diffptr table will contain at least 6548 entries (up to and including
Line 235  init_tinyprimes_tridiv(byteptr p); /* in ifactor2.c */
Line 245  init_tinyprimes_tridiv(byteptr p); /* in ifactor2.c */
 static ulong _maxprime = 0;  static ulong _maxprime = 0;
   
 ulong  ulong
 maxprime() { return _maxprime; }  maxprime(void) { return _maxprime; }
   
 byteptr  byteptr
 initprimes(long maxnum)  initprimes(ulong maxnum)
 {  {
   long len, last;    long len;
     ulong last;
   byteptr p;    byteptr p;
   /* The algorithm must see the next prime beyond maxnum, whence the +257. */    /* The algorithm must see the next prime beyond maxnum, whence the +512. */
   /* switch to unsigned: adding the 257 _could_ wrap to a negative number. */    /* switch to unsigned: adding the 512 _could_ wrap to a negative number. */
   ulong maxnum1 = ((maxnum<65302?65302:maxnum)+257ul);    ulong maxnum1 = ((maxnum<65302?65302:maxnum)+512ul);
   
   if (maxnum1 > 436273000)    if ((maxnum>>1) > VERYBIGINT - 1024)
     err(talker,"impossible to have prestored primes > 436273009");        err(talker, "Too large primelimit");
   
   p = initprimes0(maxnum1, &len, &last);    p = initprimes0(maxnum1, &len, &last);
 #if 0 /* not yet... GN */  #if 0 /* not yet... GN */
   static int build_the_tables = 1;    static int build_the_tables = 1;
Line 273  cleanprimetab(void)
Line 283  cleanprimetab(void)
 GEN  GEN
 addprimes(GEN p)  addprimes(GEN p)
 {  {
   ulong av;    gpmem_t av;
   long i,k,tx,lp;    long i,k,tx,lp;
   GEN L;    GEN L;
   
Line 281  addprimes(GEN p)
Line 291  addprimes(GEN p)
   tx = typ(p);    tx = typ(p);
   if (is_vec_t(tx))    if (is_vec_t(tx))
   {    {
     for (i=1; i < lg(p); i++) addprimes((GEN) p[i]);      for (i=1; i < lg(p); i++) (void)addprimes((GEN) p[i]);
     return primetab;      return primetab;
   }    }
   if (tx != t_INT) err(typeer,"addprime");    if (tx != t_INT) err(typeer,"addprime");
Line 302  addprimes(GEN p)
Line 312  addprimes(GEN p)
       gunclone(n); primetab[i] = 0;        gunclone(n); primetab[i] = 0;
     }      }
   }    }
   primetab = (GEN) gprealloc(primetab, (lp+1)*sizeof(long), lp*sizeof(long));    primetab = (GEN) gprealloc(primetab, (lp+1)*sizeof(long));
   primetab[i] = lclone(p); setlg(primetab, lp+1);    primetab[i] = lclone(p); setlg(primetab, lp+1);
   if (k > 1) { cleanprimetab(); setlg(L,k); addprimes(L); }    if (k > 1) { cleanprimetab(); setlg(L,k); (void)addprimes(L); }
   avma = av; return primetab;    avma = av; return primetab;
 }  }
   
Line 340  removeprimes(GEN prime)
Line 350  removeprimes(GEN prime)
     }      }
     else      else
     {      {
       for (i=1; i < lg(prime); i++) removeprime((GEN) prime[i]);        for (i=1; i < lg(prime); i++) (void)removeprime((GEN) prime[i]);
     }      }
     return primetab;      return primetab;
   }    }
Line 361  removeprimes(GEN prime)
Line 371  removeprimes(GEN prime)
 /* Some overkill removed from this  (15 spsp for an integer < 2^32 ?!).  /* Some overkill removed from this  (15 spsp for an integer < 2^32 ?!).
    Should really revert to isprime() once the new primality tester is ready     Should really revert to isprime() once the new primality tester is ready
    --GN */     --GN */
   #if 0
 #define pseudoprime(p) millerrabin(p,3*lgefint(p))  #define pseudoprime(p) millerrabin(p,3*lgefint(p))
   #else /* remove further overkill :-)  --KB */
   #define pseudoprime(p) BSW_psp(p)
   #endif
   
 /* where to stop trial dividing in factorization */  /* where to stop trial dividing in factorization */
   
Line 408  GEN 
Line 422  GEN 
 auxdecomp1(GEN n, long (*ifac_break)(GEN n, GEN pairs, GEN here, GEN state),  auxdecomp1(GEN n, long (*ifac_break)(GEN n, GEN pairs, GEN here, GEN state),
                   GEN state, long all, long hint)                    GEN state, long all, long hint)
 {  {
   long pp[] = { evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3),2,0 };    gpmem_t av;
   long nb = 0,i,k,lim1,av,lp;    long pp[] = { evaltyp(t_INT)|_evallg(4), 0,0,0 };
     long nb = 0,i,k,lp,p,lim1;
   byteptr d=diffptr+1;    byteptr d=diffptr+1;
   
   if (typ(n) != t_INT) err(arither1);    if (typ(n) != t_INT) err(arither1);
   i=signe(n); if (!i) err(arither2);    i=signe(n); if (!i) err(arither2);
   cgetg(3,t_MAT);    (void)cgetg(3,t_MAT);
   if (i<0) { stoi(-1); stoi(1); nb++; }    if (i<0) { (void)stoi(-1); (void)stoi(1); nb++; }
   if (is_pm1(n)) return aux_end(NULL,nb);    if (is_pm1(n)) return aux_end(NULL,nb);
   
   n = gclone(n);  setsigne(n,1);    n = gclone(n);  setsigne(n,1);
   i = vali(n);    i = vali(n);
   if (i)    if (i)
   {    {
     stoi(2); stoi(i); nb++;      (void)stoi(2); (void)stoi(i); nb++;
     av=avma; affii(shifti(n,-i), n); avma=av;      av=avma; affii(shifti(n,-i), n); avma=av;
   }    }
   if (is_pm1(n)) return aux_end(n,nb);    if (is_pm1(n)) return aux_end(n,nb);
   lim1 = tridiv_bound(n, all);    lim1 = tridiv_bound(n, all);
   
   /* trial division */    /* trial division */
   while (*d && pp[2] <= lim1)    p = 2;
     while (*d && p <= lim1)
   {    {
     pp[2] += *d++;      NEXT_PRIME_VIADIFF(p,d);
     if (mpdivis(n,pp,n))      if (mpdivisis(n,p,n))
     {      {
       nb++; k=1; while (mpdivis(n,pp,n)) k++;        nb++; k=1; while (mpdivisis(n,p,n)) k++;
       icopy(pp); stoi(k);        (void)utoi(p); (void)stoi(k);
       if (is_pm1(n)) return aux_end(n,nb);        if (is_pm1(n)) return aux_end(n,nb);
     }      }
   }    }
   
   /* pp = square of biggest p tried so far */    /* pp = square of biggest p tried so far */
   av=avma; setlg(pp,4); affii(sqri(pp),pp); avma=av;    av=avma; affii(muluu(p,p),pp); avma=av;
   if (cmpii(pp,n) > 0)    if (cmpii(pp,n) > 0)
   {    {
     nb++;      nb++;
     icopy(n); stoi(1); return aux_end(n,nb);      (void)icopy(n); (void)stoi(1); return aux_end(n,nb);
   }    }
   
   /* trial divide by the "special primes" (usually huge composites...) */    /* trial divide by the "special primes" (usually huge composites...) */
Line 454  auxdecomp1(GEN n, long (*ifac_break)(GEN n, GEN pairs,
Line 470  auxdecomp1(GEN n, long (*ifac_break)(GEN n, GEN pairs,
     if (mpdivis(n,(GEN) primetab[i],n))      if (mpdivis(n,(GEN) primetab[i],n))
     {      {
       nb++; k=1; while (mpdivis(n,(GEN) primetab[i],n)) k++;        nb++; k=1; while (mpdivis(n,(GEN) primetab[i],n)) k++;
       icopy((GEN) primetab[i]); stoi(k);        (void)icopy((GEN) primetab[i]); (void)stoi(k);
       if (is_pm1(n)) return aux_end(n,nb);        if (is_pm1(n)) return aux_end(n,nb);
     }      }
   
Line 462  auxdecomp1(GEN n, long (*ifac_break)(GEN n, GEN pairs,
Line 478  auxdecomp1(GEN n, long (*ifac_break)(GEN n, GEN pairs,
   if ((k && cmpii(pp,n) > 0) || (all==1 && pseudoprime(n)))    if ((k && cmpii(pp,n) > 0) || (all==1 && pseudoprime(n)))
   {    {
     nb++;      nb++;
     icopy(n); stoi(1); return aux_end(n,nb);      (void)icopy(n); (void)stoi(1); return aux_end(n,nb);
   }    }
   
   /* now we have a large composite.  Use hint as is if all==1 */    /* now we have a large composite.  Use hint as is if all==1 */
Line 484  auxdecomp1(GEN n, long (*ifac_break)(GEN n, GEN pairs,
Line 500  auxdecomp1(GEN n, long (*ifac_break)(GEN n, GEN pairs,
 static long  static long
 ifac_break_limit(GEN n, GEN pairs/*unused*/, GEN here, GEN state)  ifac_break_limit(GEN n, GEN pairs/*unused*/, GEN here, GEN state)
 {  {
   ulong ltop = avma;    gpmem_t ltop = avma;
   int res;    int res;
   if (!here) /* initial call */    if (!here) /* initial call */
    /* Small prime have been removed since start, n is the new unfactored part.     /* Small prime have been removed since start, n is the new unfactored part.
Line 557  GEN
Line 573  GEN
 boundfact(GEN n, long lim)  boundfact(GEN n, long lim)
 {  {
   GEN p1,p2,p3,p4,p5,y;    GEN p1,p2,p3,p4,p5,y;
   long av = avma,tetpil;    gpmem_t av = avma,tetpil;
   
   if (lim<=1) lim=0;    if (lim<=1) lim=0;
   switch(typ(n))    switch(typ(n))
Line 565  boundfact(GEN n, long lim)
Line 581  boundfact(GEN n, long lim)
     case t_INT:      case t_INT:
       return auxdecomp(n,lim);        return auxdecomp(n,lim);
     case t_FRACN:      case t_FRACN:
       n=gred(n); /* fall through */        n=gred(n);
         if (typ(n) == t_INT)
           return gerepileupto(av, boundfact(n,lim));
         /* fall through */
     case t_FRAC:      case t_FRAC:
       p1=auxdecomp((GEN)n[1],lim);        p1=auxdecomp((GEN)n[1],lim);
       p2=auxdecomp((GEN)n[2],lim);        p2=auxdecomp((GEN)n[2],lim);
Line 581  boundfact(GEN n, long lim)
Line 600  boundfact(GEN n, long lim)
   err(arither1);    err(arither1);
   return NULL; /* not reached */    return NULL; /* not reached */
 }  }
   /***********************************************************************/
   /**                                                                   **/
   /**                    SIMPLE FACTORISATIONS                          **/
   /**                                                                   **/
   /***********************************************************************/
   
   /*Factorize n and output [fp,fe] where
    * fp and fe are vecsmall and n=prod{fp[i]^fe[i]}
    */
   
   GEN
   decomp_small(long n)
   {
     gpmem_t ltop=avma;
     GEN F = factor(stoi(n));
     long i, l=lg(F[1]);
     GEN f=cgetg(3,t_VEC);
     GEN fp=cgetg(l,t_VECSMALL);
     GEN fe=cgetg(l,t_VECSMALL);
     f[1] = (long) fp;
     f[2] = (long) fe;
     for(i = 1; i < l; i++)
     {
       fp[i]=itos(gcoeff(F,i,1));
       fe[i]=itos(gcoeff(F,i,2));
     }
     return gerepileupto(ltop,f);
   }
   
   /*Return the primary factors of a small integer as a vecsmall*/
   GEN
   decomp_primary_small(long n)
   {
     gpmem_t ltop=avma;
     GEN F = factor(stoi(n));
     GEN fc = cgetg(lg(F[1]), t_VECSMALL);
     gpmem_t av=avma;
     long i;
     for (i = 1; i < lg(fc); i++)
       fc[i] = itos(powgi(gmael(F,1,i), gmael(F,2,i)));
     avma=av;
     return gerepileupto(ltop,fc);
   }
   
   
 /***********************************************************************/  /***********************************************************************/
 /**                                                                   **/  /**                                                                   **/
 /**                    BASIC ARITHMETIC FUNCTIONS                     **/  /**                    BASIC ARITHMETIC FUNCTIONS                     **/
Line 623  long
Line 686  long
 mu(GEN n)  mu(GEN n)
 {  {
   byteptr d = diffptr+1;        /* point at 3 - 2 */    byteptr d = diffptr+1;        /* point at 3 - 2 */
   ulong p,lim1, av = avma;    gpmem_t av = avma;
   long s, v;    ulong p;
     long s, v, lim1;
   
   if (typ(n) != t_INT) err(arither1);    if (typ(n) != t_INT) err(arither1);
   if (!signe(n)) err(arither2);    if (!signe(n)) err(arither2);
Line 638  mu(GEN n)
Line 702  mu(GEN n)
   
   while (*d && p < lim1)    while (*d && p < lim1)
   {    {
     p += *d++;      NEXT_PRIME_VIADIFF(p,d);
     if (mpdivisis(n,p,n))      if (mpdivisis(n,p,n))
     {      {
       if (smodis(n,p) == 0) { avma=av; return 0; }        if (smodis(n,p) == 0) { avma=av; return 0; }
Line 664  gissquarefree(GEN x)
Line 728  gissquarefree(GEN x)
 long  long
 issquarefree(GEN x)  issquarefree(GEN x)
 {  {
   ulong av = avma,lim1,p;    ulong p;
   long tx;    gpmem_t av = avma;
     long tx, lim1;
   GEN d;    GEN d;
   
   if (gcmp0(x)) return 0;    if (gcmp0(x)) return 0;
Line 683  issquarefree(GEN x)
Line 748  issquarefree(GEN x)
     p = 2;      p = 2;
     while (*d && p < lim1)      while (*d && p < lim1)
     {      {
       p += *d++;        NEXT_PRIME_VIADIFF(p,d);
       if (mpdivisis(x,p,x))        if (mpdivisis(x,p,x))
       {        {
         if (smodis(x,p) == 0) { avma = av; return 0; }          if (smodis(x,p) == 0) { avma = av; return 0; }
Line 709  long
Line 774  long
 omega(GEN n)  omega(GEN n)
 {  {
   byteptr d=diffptr+1;    byteptr d=diffptr+1;
   ulong p, lim1, av = avma;    gpmem_t av = avma;
   long nb,v;    long p,nb,v, lim1;
   
   if (typ(n) != t_INT) err(arither1);    if (typ(n) != t_INT) err(arither1);
   if (!signe(n)) err(arither2);    if (!signe(n)) err(arither2);
Line 723  omega(GEN n)
Line 788  omega(GEN n)
   
   while (*d && p < lim1)    while (*d && p < lim1)
   {    {
     p += *d++;      NEXT_PRIME_VIADIFF(p,d);
     if (mpdivisis(n,p,n))      if (mpdivisis(n,p,n))
     {      {
       nb++; while (mpdivisis(n,p,n)); /* empty */        nb++; while (mpdivisis(n,p,n)); /* empty */
Line 746  long
Line 811  long
 bigomega(GEN n)  bigomega(GEN n)
 {  {
   byteptr d=diffptr+1;    byteptr d=diffptr+1;
   ulong av = avma, p,lim1;    ulong p;
   long nb,v;    gpmem_t av = avma;
     long nb,v, lim1;
   
   if (typ(n) != t_INT) err(arither1);    if (typ(n) != t_INT) err(arither1);
   if (!signe(n)) err(arither2);    if (!signe(n)) err(arither2);
Line 759  bigomega(GEN n)
Line 825  bigomega(GEN n)
   
   while (*d && p < lim1)    while (*d && p < lim1)
   {    {
     p += *d++;      NEXT_PRIME_VIADIFF(p,d);
     if (mpdivisis(n,p,n))      if (mpdivisis(n,p,n))
     {      {
       do nb++; while (mpdivisis(n,p,n));        do nb++; while (mpdivisis(n,p,n));
Line 782  phi(GEN n)
Line 848  phi(GEN n)
 {  {
   byteptr d = diffptr+1;    byteptr d = diffptr+1;
   GEN m;    GEN m;
   ulong av = avma, lim1, p;    ulong p;
   long v;    gpmem_t av = avma;
     long v, lim1;
   
   if (typ(n) != t_INT) err(arither1);    if (typ(n) != t_INT) err(arither1);
   if (!signe(n)) err(arither2);    if (!signe(n)) err(arither2);
Line 796  phi(GEN n)
Line 863  phi(GEN n)
   
   while (*d && p < lim1)    while (*d && p < lim1)
   {    {
     p += *d++;      NEXT_PRIME_VIADIFF(p,d);
     if (mpdivisis(n,p,n))      if (mpdivisis(n,p,n))
     {      {
       m = mulis(m, p-1);        m = mulis(m, p-1);
Line 824  numbdiv(GEN n)
Line 891  numbdiv(GEN n)
 {  {
   byteptr d=diffptr+1;    byteptr d=diffptr+1;
   GEN m;    GEN m;
   long l, v;    long l, v, lim1;
   ulong av = avma, p,lim1;    ulong p;
     gpmem_t av = avma;
   
   if (typ(n) != t_INT) err(arither1);    if (typ(n) != t_INT) err(arither1);
   if (!signe(n)) err(arither2);    if (!signe(n)) err(arither2);
Line 838  numbdiv(GEN n)
Line 906  numbdiv(GEN n)
   
   while (*d && p < lim1)    while (*d && p < lim1)
   {    {
     p += *d++;      NEXT_PRIME_VIADIFF(p,d);
     l=1; while (mpdivisis(n,p,n)) l++;      l=1; while (mpdivisis(n,p,n)) l++;
     m = mulsi(l,m); if (is_pm1(n)) { return gerepileupto(av,m); }      m = mulsi(l,m); if (is_pm1(n)) { return gerepileupto(av,m); }
   }    }
Line 862  sumdiv(GEN n)
Line 930  sumdiv(GEN n)
 {  {
   byteptr d=diffptr+1;    byteptr d=diffptr+1;
   GEN m,m1;    GEN m,m1;
   ulong av=avma,lim1,p;    ulong p;
   long v;    gpmem_t av=avma;
     long v, lim1;
   
   if (typ(n) != t_INT) err(arither1);    if (typ(n) != t_INT) err(arither1);
   if (!signe(n)) err(arither2);    if (!signe(n)) err(arither2);
Line 876  sumdiv(GEN n)
Line 945  sumdiv(GEN n)
   
   while (*d && p < lim1)    while (*d && p < lim1)
   {    {
     p += *d++;      NEXT_PRIME_VIADIFF(p,d);
     if (mpdivisis(n,p,n))      if (mpdivisis(n,p,n))
     {      {
       m1 = utoi(p+1);        m1 = utoi(p+1);
Line 904  sumdivk(GEN n, long k)
Line 973  sumdivk(GEN n, long k)
 {  {
   byteptr d=diffptr+1;    byteptr d=diffptr+1;
   GEN n1,m,m1,pk;    GEN n1,m,m1,pk;
   ulong av = avma, p, lim1;    ulong p;
   long k1,v;    gpmem_t av = avma;
     long k1,v, lim1;
   
   if (!k) return numbdiv(n);    if (!k) return numbdiv(n);
   if (k==1) return sumdiv(n);    if (k==1) return sumdiv(n);
Line 924  sumdivk(GEN n, long k)
Line 994  sumdivk(GEN n, long k)
   
   while (*d && p < lim1)    while (*d && p < lim1)
   {    {
     p += *d++;      NEXT_PRIME_VIADIFF(p,d);
     if (mpdivisis(n,p,n))      if (mpdivisis(n,p,n))
     {      {
       pk = gpowgs(utoi(p),k); m1 = addsi(1,pk);        pk = gpowgs(utoi(p),k); m1 = addsi(1,pk);
Line 954  sumdivk(GEN n, long k)
Line 1024  sumdivk(GEN n, long k)
 GEN  GEN
 divisors(GEN n)  divisors(GEN n)
 {  {
   long tetpil,av=avma,i,j,l;    gpmem_t tetpil,av=avma;
     long i,j,l;
   GEN *d,*t,*t1,*t2,*t3, nbdiv,e;    GEN *d,*t,*t1,*t2,*t3, nbdiv,e;
   
   if (typ(n) != t_MAT || lg(n) != 3) n = auxdecomp(n,1);    if (typ(n) != t_MAT || lg(n) != 3) n = auxdecomp(n,1);
Line 978  divisors(GEN n)
Line 1049  divisors(GEN n)
   tetpil=avma; return gerepile(av,tetpil,sort((GEN)t));    tetpil=avma; return gerepile(av,tetpil,sort((GEN)t));
 }  }
   
 GEN  static GEN
 core(GEN n)  _core(GEN n, int all)
 {  {
   long av=avma,tetpil,i;    gpmem_t av = avma;
   GEN fa,p1,p2,res=gun;    long i;
     GEN fa,p1,p2,c = gun;
   
   fa=auxdecomp(n,1); p1=(GEN)fa[1]; p2=(GEN)fa[2];    fa = auxdecomp(n,all);
     p1 = (GEN)fa[1];
     p2 = (GEN)fa[2];
   for (i=1; i<lg(p1); i++)    for (i=1; i<lg(p1); i++)
     if (mod2((GEN)p2[i])) res = mulii(res,(GEN)p1[i]);      if (mod2((GEN)p2[i])) c = mulii(c,(GEN)p1[i]);
   tetpil=avma; return gerepile(av,tetpil,icopy(res));    return gerepileuptoint(av, c);
 }  }
   
 GEN  static GEN
 core2(GEN n)  _core2(GEN n, int all)
 {  {
   long av=avma,tetpil,i;    gpmem_t av = avma;
     long i;
     GEN fa,p1,p2,e,c=gun,f=gun,y;
   
   GEN fa,p1,p2,p3,res=gun,res2=gun,y;    fa = auxdecomp(n,all);
     p1 = (GEN)fa[1];
   fa=auxdecomp(n,1); p1=(GEN)fa[1]; p2=(GEN)fa[2];    p2 = (GEN)fa[2];
   for (i=1; i<lg(p1); i++)    for (i=1; i<lg(p1); i++)
   {    {
     p3=(GEN)p2[i];      e = (GEN)p2[i];
     if (mod2(p3)) res=mulii(res,(GEN)p1[i]);      if (mod2(e))   c = mulii(c, (GEN)p1[i]);
     if (!gcmp1(p3)) res2=mulii(res2,gpui((GEN)p1[i],shifti(p3,-1),0));      if (!gcmp1(e)) f = mulii(f, powgi((GEN)p1[i], shifti(e,-1)));
   }    }
   tetpil=avma; y=cgetg(3,t_VEC);    y = cgetg(3,t_VEC);
   y[1]=(long)icopy(res); y[2]=(long)icopy(res2);    y[1] = (long)c;
   return gerepile(av,tetpil,y);    y[2] = (long)f; return gerepilecopy(av, y);
 }  }
   
   GEN core(GEN n)        { return _core(n,1); }
   GEN corepartial(GEN n) { return _core(n,0); }
   GEN core2(GEN n)       { return _core2(n,1); }
   GEN core2partial(GEN n){ return _core2(n,0); }
   
 GEN  GEN
 core0(GEN n,long flag)  core0(GEN n,long flag) { return flag? core2(n): core(n); }
 {  
   return flag? core2(n): core(n);  
 }  
   
   static long
   _mod4(GEN c) { long r = mod4(c); if (signe(c) < 0) r = 4-r; return r; }
   
 GEN  GEN
 coredisc(GEN n)  coredisc(GEN n)
 {  {
   long av=avma,tetpil,r;    gpmem_t av = avma;
   GEN p1=core(n);    GEN c = core(n);
     long r = _mod4(c);
   r=mod4(p1); if (signe(p1)<0) r = 4-r;    if (r==1 || r==4) return c;
   if (r==1 || r==4) return p1;    return gerepileuptoint(av, shifti(c,2));
   tetpil=avma; return gerepile(av,tetpil,shifti(p1,2));  
 }  }
   
 GEN  GEN
 coredisc2(GEN n)  coredisc2(GEN n)
 {  {
   long av=avma,tetpil,r;    gpmem_t av = avma;
   GEN y,p1,p2=core2(n);    GEN y = core2(n);
     GEN c = (GEN)y[1], f = (GEN)y[2];
   p1=(GEN)p2[1]; r=mod4(p1); if (signe(p1)<0) r=4-r;    long r = _mod4(c);
   if (r==1 || r==4) return p2;    if (r==1 || r==4) return y;
   tetpil=avma; y=cgetg(3,t_VEC);    y = cgetg(3,t_VEC);
   y[1]=lshifti(p1,2); y[2]=lmul2n((GEN)p2[2],-1);    y[1] = lshifti(c,2);
   return gerepile(av,tetpil,y);    y[2] = lmul2n(f,-1); return gerepileupto(av, y);
 }  }
   
 GEN  GEN
 coredisc0(GEN n,long flag)  coredisc0(GEN n,long flag) { return flag? coredisc2(n): coredisc(n); }
 {  
   return flag? coredisc2(n): coredisc(n);  
 }  
   
 /***********************************************************************/  /***********************************************************************/
 /**                                                                   **/  /**                                                                   **/
Line 1065  qf_create(GEN x, GEN y, GEN z, long s)
Line 1142  qf_create(GEN x, GEN y, GEN z, long s)
   if (typ(x)!=t_INT || typ(y)!=t_INT || typ(z)!=t_INT) err(typeer,"Qfb");    if (typ(x)!=t_INT || typ(y)!=t_INT || typ(z)!=t_INT) err(typeer,"Qfb");
   if (!s)    if (!s)
   {    {
     long av=avma; s = signe(qf_disc(x,y,z)); avma=av;      gpmem_t av=avma; s = signe(qf_disc(x,y,z)); avma=av;
     if (!s) err(talker,"zero discriminant in Qfb");      if (!s) err(talker,"zero discriminant in Qfb");
   }    }
   t = (s>0)? cgetg(5,t_QFR): cgetg(4,t_QFI);    t = (s>0)? cgetg(5,t_QFR): cgetg(4,t_QFI);
Line 1155  comp_gen(GEN z,GEN x,GEN y)
Line 1232  comp_gen(GEN z,GEN x,GEN y)
 static GEN  static GEN
 compimag0(GEN x, GEN y, int raw)  compimag0(GEN x, GEN y, int raw)
 {  {
   ulong tx = typ(x), av = avma;    gpmem_t av = avma;
     long tx = typ(x);
   GEN z;    GEN z;
   
   if (typ(y) != tx || tx!=t_QFI) err(typeer,"composition");    if (typ(y) != tx || tx!=t_QFI) err(typeer,"composition");
Line 1168  compimag0(GEN x, GEN y, int raw)
Line 1246  compimag0(GEN x, GEN y, int raw)
 static GEN  static GEN
 compreal0(GEN x, GEN y, int raw)  compreal0(GEN x, GEN y, int raw)
 {  {
   ulong tx = typ(x), av = avma;    gpmem_t av = avma;
     long tx = typ(x);
   GEN z;    GEN z;
   
   if (typ(y) != tx || tx!=t_QFR) err(typeer,"composition");    if (typ(y) != tx || tx!=t_QFR) err(typeer,"composition");
Line 1199  compraw(GEN x, GEN y)
Line 1278  compraw(GEN x, GEN y)
 GEN  GEN
 sqcompimag0(GEN x, long raw)  sqcompimag0(GEN x, long raw)
 {  {
   long av = avma;    gpmem_t av = avma;
   GEN z = cgetg(4,t_QFI);    GEN z = cgetg(4,t_QFI);
   
   if (typ(x)!=t_QFI) err(typeer,"composition");    if (typ(x)!=t_QFI) err(typeer,"composition");
Line 1211  sqcompimag0(GEN x, long raw)
Line 1290  sqcompimag0(GEN x, long raw)
 static GEN  static GEN
 sqcompreal0(GEN x, long raw)  sqcompreal0(GEN x, long raw)
 {  {
   long av = avma;    gpmem_t av = avma;
   GEN z = cgetg(5,t_QFR);    GEN z = cgetg(5,t_QFR);
   
   if (typ(x)!=t_QFR) err(typeer,"composition");    if (typ(x)!=t_QFR) err(typeer,"composition");
Line 1236  static GEN
Line 1315  static GEN
 real_unit_form_by_disc(GEN D, long prec)  real_unit_form_by_disc(GEN D, long prec)
 {  {
   GEN y = cgetg(5,t_QFR), isqrtD;    GEN y = cgetg(5,t_QFR), isqrtD;
   long av = avma;    gpmem_t av = avma;
   
   if (typ(D) != t_INT || signe(D) <= 0) err(typeer,"real_unit_form_by_disc");    if (typ(D) != t_INT || signe(D) <= 0) err(typeer,"real_unit_form_by_disc");
   switch(mod4(D))    switch(mod4(D))
Line 1249  real_unit_form_by_disc(GEN D, long prec)
Line 1328  real_unit_form_by_disc(GEN D, long prec)
   if (mod2(D) != mod2(isqrtD))    if (mod2(D) != mod2(isqrtD))
     isqrtD = gerepileuptoint(av, addsi(-1,isqrtD));      isqrtD = gerepileuptoint(av, addsi(-1,isqrtD));
   y[2] = (long)isqrtD; av = avma;    y[2] = (long)isqrtD; av = avma;
   y[3] = (long)gerepileuptoint(av, shifti(subii(sqri(isqrtD),D),-2));    y[3] = lpileuptoint(av, shifti(subii(sqri(isqrtD),D),-2));
   y[4] = (long)realzero(prec); return y;    y[4] = (long)realzero(prec); return y;
 }  }
   
 GEN  GEN
 real_unit_form(GEN x)  real_unit_form(GEN x)
 {  {
   long av = avma, prec;    gpmem_t av = avma;
     long prec;
   GEN D;    GEN D;
   if (typ(x) != t_QFR) err(typeer,"real_unit_form");    if (typ(x) != t_QFR) err(typeer,"real_unit_form");
   prec = precision((GEN)x[4]);    prec = precision((GEN)x[4]);
Line 1283  imag_unit_form_by_disc(GEN D)
Line 1363  imag_unit_form_by_disc(GEN D)
   y[3] = lshifti(D,-2); setsigne(y[3],1);    y[3] = lshifti(D,-2); setsigne(y[3],1);
   if (isodd)    if (isodd)
   {    {
     long av = avma;      gpmem_t av = avma;
     y[3] = (long)gerepileuptoint(av, addis((GEN)y[3],1));      y[3] = lpileuptoint(av, addis((GEN)y[3],1));
   }    }
   return y;    return y;
 }  }
Line 1293  GEN
Line 1373  GEN
 imag_unit_form(GEN x)  imag_unit_form(GEN x)
 {  {
   GEN p1,p2, y = cgetg(4,t_QFI);    GEN p1,p2, y = cgetg(4,t_QFI);
   long av;    gpmem_t av;
   if (typ(x) != t_QFI) err(typeer,"imag_unit_form");    if (typ(x) != t_QFI) err(typeer,"imag_unit_form");
   y[1] = un;    y[1] = un;
   y[2] = mpodd((GEN)x[2])? un: zero;    y[2] = mpodd((GEN)x[2])? un: zero;
   av = avma; p1 = mulii((GEN)x[1],(GEN)x[3]);    av = avma; p1 = mulii((GEN)x[1],(GEN)x[3]);
   p2 = shifti(sqri((GEN)x[2]),-2);    p2 = shifti(sqri((GEN)x[2]),-2);
   y[3] = (long)gerepileuptoint(av, subii(p1,p2));    y[3] = lpileuptoint(av, subii(p1,p2));
   return y;    return y;
 }  }
   
 GEN  GEN
 powrealraw(GEN x, long n)  powrealraw(GEN x, long n)
 {  {
   long av = avma, m;    gpmem_t av = avma;
     long m;
   GEN y;    GEN y;
   
   if (typ(x) != t_QFR)    if (typ(x) != t_QFR)
Line 1329  powrealraw(GEN x, long n)
Line 1410  powrealraw(GEN x, long n)
 GEN  GEN
 powimagraw(GEN x, long n)  powimagraw(GEN x, long n)
 {  {
   long av = avma, m;    gpmem_t av = avma;
     long m;
   GEN y;    GEN y;
   
   if (typ(x) != t_QFI)    if (typ(x) != t_QFI)
Line 1360  powraw(GEN x, long n)
Line 1442  powraw(GEN x, long n)
 GEN  GEN
 nucomp(GEN x, GEN y, GEN l)  nucomp(GEN x, GEN y, GEN l)
 {  {
   long av=avma,tetpil,cz;    gpmem_t av=avma,tetpil;
     long cz;
   GEN a,a1,a2,b,b2,d,d1,e,g,n,p1,p2,p3,q1,q2,q3,q4,s,t2,t3,u,u1,v,v1,v2,v3,z;    GEN a,a1,a2,b,b2,d,d1,e,g,n,p1,p2,p3,q1,q2,q3,q4,s,t2,t3,u,u1,v,v1,v2,v3,z;
   
   if (x==y) return nudupl(x,l);    if (x==y) return nudupl(x,l);
Line 1424  nucomp(GEN x, GEN y, GEN l)
Line 1507  nucomp(GEN x, GEN y, GEN l)
 GEN  GEN
 nudupl(GEN x, GEN l)  nudupl(GEN x, GEN l)
 {  {
   long av=avma,tetpil,cz;    gpmem_t av=avma,tetpil;
     long cz;
   GEN u,v,d,d1,p1,a,b,c,b2,z,v2,v3,t2,t3,e,g;    GEN u,v,d,d1,p1,a,b,c,b2,z,v2,v3,t2,t3,e,g;
   
   if (typ(x) != t_QFI)    if (typ(x) != t_QFI)
Line 1464  nudupl(GEN x, GEN l)
Line 1548  nudupl(GEN x, GEN l)
 GEN  GEN
 nupow(GEN x, GEN n)  nupow(GEN x, GEN n)
 {  {
   long av,tetpil,i,j;    gpmem_t av,tetpil;
     long i,j;
   unsigned long m;    unsigned long m;
   GEN y,l;    GEN y,l;
   
Line 1551  rhoreal_aux(GEN x, GEN D, GEN sqrtD, GEN isqrtD)
Line 1636  rhoreal_aux(GEN x, GEN D, GEN sqrtD, GEN isqrtD)
   GEN p1,p2, y = cgetg(6,t_VEC);    GEN p1,p2, y = cgetg(6,t_VEC);
   GEN b = (GEN)x[2];    GEN b = (GEN)x[2];
   GEN c = (GEN)x[3];    GEN c = (GEN)x[3];
   long s = signe(c);  
   
   y[1] = (long)c;    y[1] = (long)c;
   p2 = (absi_cmp(isqrtD,c) >= 0)? isqrtD: c;    p2 = (absi_cmp(isqrtD,c) >= 0)? isqrtD: absi(c);
   setsigne(c,1);  
   p1 = shifti(c,1);    p1 = shifti(c,1);
     if (p1 == gzero) err(talker, "reducible form in rhoreal");
     setsigne(p1,1); /* |2c| */
   p2 = divii(addii(p2,b), p1);    p2 = divii(addii(p2,b), p1);
   setsigne(c,s);  
   y[2] = lsubii(mulii(p2,p1), b);    y[2] = lsubii(mulii(p2,p1), b);
   
   p1 = shifti(subii(sqri((GEN)y[2]),D),-2);    p1 = shifti(subii(sqri((GEN)y[2]),D),-2);
Line 1586  GEN
Line 1670  GEN
 codeform5(GEN x, long prec)  codeform5(GEN x, long prec)
 {  {
   GEN y = cgetg(6,t_VEC);    GEN y = cgetg(6,t_VEC);
   y[1]=x[1];    y[1] = x[1];
   y[2]=x[2];    y[2] = x[2];
   y[3]=x[3];    y[3] = x[3];
   y[4]=zero;    y[4] = zero;
   y[5]=(long)realun(prec); return y;    y[5] = (long)realun(prec); return y;
 }  }
   
 static GEN  static GEN
Line 1619  decodeform(GEN x, GEN d0)
Line 1703  decodeform(GEN x, GEN d0)
     p1 = shiftr(p1,-e);      p1 = shiftr(p1,-e);
     p2 = addis(mulsi(EXP220,p2), e);      p2 = addis(mulsi(EXP220,p2), e);
     p1 = mplog(p1);      p1 = mplog(p1);
     p1 = mpadd(p1, mulir(p2, glog(gdeux, lg(d0))));      p1 = mpadd(p1, mulir(p2, mplog2(lg(d0))));
   }    }
   else    else
   { /* to avoid loss of precision */    { /* to avoid loss of precision */
Line 1665  redrealform5(GEN x, GEN D, GEN sqrtD, GEN isqrtD)
Line 1749  redrealform5(GEN x, GEN D, GEN sqrtD, GEN isqrtD)
 static GEN  static GEN
 redreal0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD)  redreal0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD)
 {  {
   long av = avma, prec;    gpmem_t av = avma;
     long prec;
   GEN d0;    GEN d0;
   
   if (typ(x) != t_QFR) err(talker,"not a real quadratic form in redreal");    if (typ(x) != t_QFR) err(talker,"not a real quadratic form in redreal");
Line 1703  redreal0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrt
Line 1788  redreal0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrt
 GEN  GEN
 comprealform5(GEN x, GEN y, GEN D, GEN sqrtD, GEN isqrtD)  comprealform5(GEN x, GEN y, GEN D, GEN sqrtD, GEN isqrtD)
 {  {
   ulong av = avma;    gpmem_t av = avma;
   GEN z = cgetg(6,t_VEC); comp_gen(z,x,y);    GEN z = cgetg(6,t_VEC); comp_gen(z,x,y);
   if (x == y)    if (x == y)
   {    {
Line 1723  comprealform5(GEN x, GEN y, GEN D, GEN sqrtD, GEN isqr
Line 1808  comprealform5(GEN x, GEN y, GEN D, GEN sqrtD, GEN isqr
 GEN  GEN
 powrealform(GEN x, GEN n)  powrealform(GEN x, GEN n)
 {  {
   long av = avma, i,m;    gpmem_t av = avma;
     long i,m;
   GEN y,D,sqrtD,isqrtD,d0;    GEN y,D,sqrtD,isqrtD,d0;
   
   if (typ(x) != t_QFR)    if (typ(x) != t_QFR)
Line 1755  powrealform(GEN x, GEN n)
Line 1841  powrealform(GEN x, GEN n)
 GEN  GEN
 redimag(GEN x)  redimag(GEN x)
 {  {
   long av=avma, fl;    gpmem_t av=avma;
     long fl;
   do x = rhoimag0(x, &fl); while (fl == 0);    do x = rhoimag0(x, &fl); while (fl == 0);
   x = gerepilecopy(av,x);    x = gerepilecopy(av,x);
   if (fl == 2) setsigne(x[2], -signe(x[2]));    if (fl == 2) setsigne(x[2], -signe(x[2]));
Line 1790  GEN
Line 1877  GEN
 qfbred0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD)  qfbred0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD)
 {  {
   long tx=typ(x),fl;    long tx=typ(x),fl;
   ulong av;    gpmem_t av;
   
   if (!is_qf_t(tx)) err(talker,"not a quadratic form in qfbred");    if (!is_qf_t(tx)) err(talker,"not a quadratic form in qfbred");
   if (!D) D = qf_disc(x,NULL,NULL);    if (!D) D = qf_disc(x,NULL,NULL);
Line 1815  qfbred0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD
Line 1902  qfbred0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD
 GEN  GEN
 primeform(GEN x, GEN p, long prec)  primeform(GEN x, GEN p, long prec)
 {  {
   long av,tetpil,s=signe(x);    gpmem_t av;
   GEN y,b;    long s, sx = signe(x);
     GEN y, b, c;
   
   if (typ(x) != t_INT || !s) err(arither1);    if (typ(x) != t_INT || !sx) err(arither1);
   if (typ(p) != t_INT || signe(p) <= 0) err(arither1);    if (typ(p) != t_INT || signe(p) <= 0) err(arither1);
   if (is_pm1(p))    if (is_pm1(p))
     return s<0? imag_unit_form_by_disc(x)      return sx<0? imag_unit_form_by_disc(x)
               : real_unit_form_by_disc(x,prec);                 : real_unit_form_by_disc(x,prec);
   if (s < 0)    s = mod8(x);
     if (sx < 0)
   {    {
     y = cgetg(4,t_QFI); s = 8-mod8(x);      if (s) s = 8-s;
       y = cgetg(4, t_QFI);
   }    }
   else    else
   {    {
     y = cgetg(5,t_QFR); s = mod8(x);      y = cgetg(5, t_QFR);
     y[4]=(long)realzero(prec);      y[4] = (long)realzero(prec);
   }    }
   switch(s&3)    switch(s&3)
   {    {
     case 2: case 3: err(funder2,"primeform");      case 2: case 3: err(funder2,"primeform");
   }    }
   y[1] = (long)icopy(p); av = avma;    av = avma;
   if (egalii(p,gdeux))    if (egalii(p, gdeux))
   {    {
     switch(s)      switch(s)
     {      {
       case 0: y[2]=zero; break;        case 0: b = gzero; break;
       case 8: s=0; y[2]=zero; break;        case 1: b = gun;   break;
       case 1: s=1; y[2]=un; break;        case 4: b = gdeux; break;
       case 4: s=4; y[2]=deux; break;        default: err(sqrter5); b = NULL; /* -Wall */
       default: err(sqrter5);  
     }      }
     b=subsi(s,x); tetpil=avma; b=shifti(b,-3);      c = shifti(subsi(s,x), -3);
   }    }
   else    else
   {    {
     b = mpsqrtmod(x,p); if (!b) err(sqrter5);      b = mpsqrtmod(x,p); if (!b) err(sqrter5);
     if (mod2(b) == mod2(x)) y[2] = (long)b;      s &= 1; /* s = x mod 2 */
     else      /* mod(b) != mod2(x) ? [Warning: we may have b == 0] */
     { tetpil = avma; y[2] = lpile(av,tetpil,subii(p,b)); }      if ((!signe(b) && s) || mod2(b) != s) b = gerepileuptoint(av, subii(p,b));
   
     av=avma; b=shifti(subii(sqri((GEN)y[2]),x),-2);      av = avma;
     tetpil=avma; b=divii(b,p);      c = diviiexact(shifti(subii(sqri(b), x), -2), p);
   }    }
   y[3]=lpile(av,tetpil,b); return y;    y[3] = lpileuptoint(av, c);
     y[2] = (long)b;
     y[1] = licopy(p); return y;
 }  }
   
 /*********************************************************************/  /*********************************************************************/
Line 1955  bittest(GEN x, long n)
Line 2046  bittest(GEN x, long n)
   return n? 1: 0;    return n? 1: 0;
 }  }
   
   GEN
   bittest_many(GEN x, GEN gn, long c)
   {
     long lx = lgefint(x), l1, l2, b1, b2, two_adic = 0, l_res, skip;
     ulong *p, *pnew;
     long n = itos(gn), extra_words = 0, partial_bits;
     GEN res;
   
     if (c == 1)                           /* Shortcut */
         return bittest(x,n) ? gun : gzero;
     /* Negative n with n+c>0 are too hairy to implement... */
     if (!signe(x) || c == 0)
         return gzero;
     if (c < 0) {                          /* Negative x means 2s complemenent */
         c = -c;
         if (signe(x) < 0)
             two_adic = 1;                 /* treat x 2-adically... */
     }
     if (n < 0) {
         gpmem_t oa, na;
   
         if (n + c <= 0)
             return gzero;
         oa = avma;
         res = bittest_many(x, gzero, two_adic? -(n+c) : n+c);
         na = avma;
         res = shifti3(res,-n,0);
         return gerepile(oa,na,res);
     }
     partial_bits = (c & (BITS_IN_LONG-1));
     l2 = lx-1 - (n>>TWOPOTBITS_IN_LONG); /* Last=least-significant word */
     l1 = lx-1 - ((n + c - 1)>>TWOPOTBITS_IN_LONG); /* First word */
     b2 = (n & (BITS_IN_LONG-1));          /* Last bit, skip less-signifant */
     b1 = ((n + c - 1) & (BITS_IN_LONG-1)); /* First bit, skip more-significant */
     if (l2 < 2) {                         /* always: l1 <= l2 */
         if (!two_adic)
             return gzero;
         /* x is non-zero, so it behaves as -1.  */
         return gbitneg(gzero,c);
     }
     if (l1 < 2) {
         if (two_adic)     /* If b1 < b2, bits are set by prepend in shift_r */
             extra_words = 2 - l1 - (b1 < b2);
         else
             partial_bits = b2 ? BITS_IN_LONG - b2 : 0;
         l1 = 2;
         b1 = (BITS_IN_LONG-1);            /* Include all bits in this word */
     }
     skip = (b1 < b2);                     /* Skip the first word of the shift */
     l_res = l2 - l1 + 1 + 2 - skip + extra_words; /* A coarse estimate */
     p = (ulong*) (new_chunk(l_res) + 2 + extra_words);
     shift_r(p - skip, (ulong*)x + l1, (ulong*)x + l2 + 1, 0, b2);
     if (two_adic) {                       /* Check the low bits of x */
           int i = lx, nonzero = 0;
   
           /* Flip the bits */
           pnew = p + l_res - 2 - extra_words;
           while (--pnew >= p)
               *pnew = MAXULONG - *pnew;
           /* Fill the extra words */
           while (extra_words--)
               *--p = MAXULONG;
           /* The result is one less than 2s-complement if the lower-bits
              of x were all 0. */
           while (--i > l2) {
               if (x[i]) {
                   nonzero = 1;
                   break;
               }
           }
           if (!nonzero && b2)             /* Check the partial word too */
               nonzero = x[l2] & ((1<<b2) - 1);
           if (!nonzero) { /* Increment res.  Do not underflow to before p */
               pnew = p + l_res - 2;
               while (--pnew >= p)
                   if (++*pnew)
                       break;
           }
     }
     if (partial_bits)
         *p &= (1<<partial_bits) - 1;
     pnew = p;
     while ((pnew < p + l_res - 2) && !*pnew) /* Skip 0 words */
         pnew++;
     avma = (gpmem_t)(pnew - 2);
     if (pnew >= p - 2 + l_res)
         return gzero;
     l_res -= (pnew - p);
     res = (GEN)(pnew - 2);
     res[1] = evalsigne(1)|evallgefint(l_res);
     res[0] = evaltyp(t_INT)|evallg(l_res);
     return res;
   }
   
 static long  static long
 bittestg(GEN x, GEN n)  bittestg(GEN x, GEN n)
 {  {
Line 1967  gbittest(GEN x, GEN n)
Line 2152  gbittest(GEN x, GEN n)
   return arith_proto2(bittestg,x,n);    return arith_proto2(bittestg,x,n);
 }  }
   
   GEN
   gbittest3(GEN x, GEN n, long c)
   {
     return garith_proto3ggs(bittest_many,x,n,c);
   }
   
 /***********************************************************************/  /***********************************************************************/
 /**                                                                   **/  /**                                                                   **/
 /**                          BITMAP OPS                               **/  /**                          BITMAP OPS                               **/
Line 2226  ibitnegimply(GEN x, GEN y)
Line 2417  ibitnegimply(GEN x, GEN y)
 }  }
   
 static GEN  static GEN
 inegate_inplace(GEN z, long ltop)  inegate_inplace(GEN z, gpmem_t ltop)
 {  {
   long lbot, o;    long o;
   
   /* Negate z */    /* Negate z */
   o = incdec(z, 1);                     /* Can overflow z... */    o = incdec(z, 1);                     /* Can overflow z... */
Line 2238  inegate_inplace(GEN z, long ltop)
Line 2429  inegate_inplace(GEN z, long ltop)
   else if (lgefint(z) == 2)    else if (lgefint(z) == 2)
       setsigne(z, 0);        setsigne(z, 0);
   incdec(z,-1);                 /* Restore a normalized value */    incdec(z,-1);                 /* Restore a normalized value */
   lbot = avma;    return gerepileupto(ltop, subis(z,1));
   z = gsub(z,gun);  
   return gerepile(ltop,lbot,z);  
 }  }
   
 GEN  GEN
 gbitor(GEN x, GEN y)  gbitor(GEN x, GEN y)
 {  {
     gpmem_t ltop;
   long sx,sy;    long sx,sy;
   long ltop;  
   GEN z;    GEN z;
   
   if (typ(x) != t_INT || typ(y) != t_INT)    if (typ(x) != t_INT || typ(y) != t_INT)
Line 2278  gbitor(GEN x, GEN y)
Line 2467  gbitor(GEN x, GEN y)
 GEN  GEN
 gbitand(GEN x, GEN y)  gbitand(GEN x, GEN y)
 {  {
     gpmem_t ltop;
   long sx,sy;    long sx,sy;
   long ltop;  
   GEN z;    GEN z;
   
   if (typ(x) != t_INT || typ(y) != t_INT)    if (typ(x) != t_INT || typ(y) != t_INT)
Line 2312  gbitand(GEN x, GEN y)
Line 2501  gbitand(GEN x, GEN y)
 GEN  GEN
 gbitxor(GEN x, GEN y)  gbitxor(GEN x, GEN y)
 {  {
     gpmem_t ltop;
   long sx,sy;    long sx,sy;
   long ltop;  
   GEN z;    GEN z;
   
   if (typ(x) != t_INT || typ(y) != t_INT)    if (typ(x) != t_INT || typ(y) != t_INT)
Line 2345  gbitxor(GEN x, GEN y)
Line 2534  gbitxor(GEN x, GEN y)
 GEN  GEN
 gbitnegimply(GEN x, GEN y)              /* x & ~y */  gbitnegimply(GEN x, GEN y)              /* x & ~y */
 {  {
     gpmem_t ltop;
   long sx,sy;    long sx,sy;
   long ltop;  
   GEN z;    GEN z;
   
   if (typ(x) != t_INT || typ(y) != t_INT)    if (typ(x) != t_INT || typ(y) != t_INT)
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>