[BACK]Return to arith1.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/arith1.c between version 1.1.1.1 and 1.2

version 1.1.1.1, 2001/10/02 11:17:01 version 1.2, 2002/09/11 07:26:48
Line 115  garith_proto2gs(GEN f(GEN,long), GEN x, long y)
Line 115  garith_proto2gs(GEN f(GEN,long), GEN x, long y)
 }  }
   
 GEN  GEN
   garith_proto3ggs(GEN f(GEN,GEN,long), GEN x, GEN y, long z)
   {
     long l,i,tx = typ(x);
     GEN t;
   
     if (is_matvec_t(tx))
     {
       l=lg(x); t=cgetg(l,tx);
       for (i=1; i<l; i++) t[i]= (long) garith_proto3ggs(f,(GEN) x[i],y,z);
       return t;
     }
     if (tx != t_INT) err(arither1);
     tx = typ(y);
     if (is_matvec_t(tx))
     {
       l=lg(y); t=cgetg(l,tx);
       for (i=1; i<l; i++) t[i]= (long) garith_proto3ggs(f,x,(GEN) y[i],z);
       return t;
     }
     if (tx != t_INT) err(arither1);
     return f(x,y,z);
   }
   
   GEN
 gassoc_proto(GEN f(GEN,GEN), GEN x, GEN y)  gassoc_proto(GEN f(GEN,GEN), GEN x, GEN y)
 {  {
   int tx=typ(x);    int tx=typ(x);
   if (!y)    if (!y)
   {    {
     ulong av=avma;      gpmem_t av = avma;
     if (tx!=t_VEC && tx!=t_COL)      if (tx!=t_VEC && tx!=t_COL)
       err(typeer,"association");        err(typeer,"association");
     return gerepileupto(av,divide_conquer_prod(x,f));      return gerepileupto(av,divide_conquer_prod(x,f));
Line 137  gassoc_proto(GEN f(GEN,GEN), GEN x, GEN y)
Line 161  gassoc_proto(GEN f(GEN,GEN), GEN x, GEN y)
 GEN  GEN
 order(GEN x)  order(GEN x)
 {  {
   long av=avma,av1,i,e;    gpmem_t av = avma,av1;
     long i,e;
   GEN o,m,p;    GEN o,m,p;
   
   if (typ(x) != t_INTMOD || !gcmp1(mppgcd((GEN) x[1],(GEN) x[2])))    if (typ(x) != t_INTMOD || !gcmp1(mppgcd((GEN) x[1],(GEN) x[2])))
Line 172  ggener(GEN m)
Line 197  ggener(GEN m)
 GEN  GEN
 gener(GEN m)  gener(GEN m)
 {  {
   long av=avma,av1,k,i,e;    gpmem_t av=avma,av1;
     long k,i,e;
   GEN x,t,q,p;    GEN x,t,q,p;
   
   if (typ(m) != t_INT) err(arither1);    if (typ(m) != t_INT) err(arither1);
   e = signe(m);    e = signe(m);
   if (!e) err(talker,"zero modulus in znprimroot");    if (!e) err(talker,"zero modulus in znprimroot");
   if (is_pm1(m)) { avma=av; return gmodulss(0,1); }    if (is_pm1(m)) return gmodulss(0,1);
   if (e<0) m = absi(m);    if (e < 0) m = absi(m);
   
   e = mod4(m);    e = mod4(m);
   if (e == 0) /* m = 0 mod 4 */    if (e == 0) /* m = 0 mod 4 */
Line 196  gener(GEN m)
Line 222  gener(GEN m)
   
   t=decomp(m); if (lg(t[1]) != 2) err(generer);    t=decomp(m); if (lg(t[1]) != 2) err(generer);
   p=gcoeff(t,1,1); e=itos(gcoeff(t,1,2)); q=subis(p,1);    p=gcoeff(t,1,1); e=itos(gcoeff(t,1,2)); q=subis(p,1);
   if (e>=2)    if (e >= 2)
   {    {
     x = (GEN) gener(p)[2];      x = (GEN)gener(p)[2];
     if (gcmp1(powmodulo(x,q, sqri(p)))) x = addii(x,p);      if (gcmp1(powmodulo(x,q, sqri(p)))) x = addii(x,p);
     av1=avma; return gerepile(av,av1,gmodulcp(x,m));      av1=avma; return gerepile(av,av1,gmodulcp(x,m));
   }    }
   
   p=(GEN)decomp(q)[1]; k=lg(p)-1; x=stoi(1);    p=(GEN)decomp(q)[1]; k=lg(p)-1; x=stoi(1);
     for (i=1; i<=k; i++) p[i] = (long)diviiexact(q, (GEN)p[i]);
   for(;;)    for(;;)
   {    {
     x[2]++;      x[2]++;
     if (gcmp1(mppgcd(m,x)))      if (gcmp1(mppgcd(m,x)))
     {      {
       for (i=k; i; i--)        for (i=k; i; i--)
         if (gcmp1(powmodulo(x, divii(q,(GEN)p[i]), m))) break;          if (gcmp1(powmodulo(x, (GEN)p[i], m))) break;
       if (!i) break;        if (!i) break;
     }      }
   }    }
   av1=avma; return gerepile(av,av1,gmodulcp(x,m));    av1=avma; return gerepile(av,av1,gmodulcp(x,m));
 }  }
   
   /* assume p odd prime */
   ulong
   u_gener(ulong p)
   {
     const gpmem_t av = avma;
     const long q = p - 1;
     const GEN L = (GEN)decomp(utoi(q))[1];
     const int k = lg(L) - 1;
     int i,x;
   
     for (x=2;;x++)
       if (x % p)
       {
         for (i=k; i; i--)
           if (powuumod(x, q/itos((GEN)L[i]), p) == 1) break;
         if (!i) break;
       }
     avma = av; return x;
   }
   
 GEN  GEN
 znstar(GEN n)  znstar(GEN n)
 {  {
   GEN p1,z,q,u,v,d,list,ep,h,gen,moduli,p,a;    GEN p1,z,q,u,v,d,list,ep,h,gen,moduli,p,a;
   long i,j,nbp,sizeh;    long i,j,nbp,sizeh;
   ulong av;    gpmem_t av;
   
   if (typ(n) != t_INT) err(arither1);    if (typ(n) != t_INT) err(arither1);
   if (!signe(n))    if (!signe(n))
Line 315  gracine(GEN a)
Line 362  gracine(GEN a)
 static GEN  static GEN
 racine_r(GEN a, long l)  racine_r(GEN a, long l)
 {  {
   long av,s;    gpmem_t av;
     long s;
   GEN x,y,z;    GEN x,y,z;
   
   if (l <= 4) return utoi(mpsqrtl(a));    if (l <= 4) return utoi(mpsqrtl(a));
Line 331  racine_r(GEN a, long l)
Line 379  racine_r(GEN a, long l)
     y = x; x = z;      y = x; x = z;
   }    }
   while (cmpii(x,y) < 0);    while (cmpii(x,y) < 0);
   avma = (long)y;    avma = (gpmem_t)y;
   return gerepileuptoint(av,y);    return gerepileuptoint(av,y);
 }  }
   
 GEN  GEN
 racine_i(GEN a, int roundup)  racine_i(GEN a, int roundup)
 {  {
   long k,m,l = lgefint(a), av = avma;    gpmem_t av = avma;
     long k,m,l = lgefint(a);
   GEN x = racine_r(a, l);    GEN x = racine_r(a, l);
   if (roundup && signe(x))    if (roundup && signe(x))
   {    {
     m = modBIL(x);      m = modBIL(x);
     k = (m * m != a[l-1] || !egalii(sqri(x),a));      k = (m * m != a[l-1] || !egalii(sqri(x),a));
     avma = (long)x;      avma = (gpmem_t)x;
     if (k) x = gerepileuptoint(av, addis(x,1));      if (k) x = gerepileuptoint(av, addis(x,1));
   }    }
   return x;    return x;
Line 409  ucarrecomplet(ulong A)
Line 458  ucarrecomplet(ulong A)
 long  long
 carrecomplet(GEN x, GEN *pt)  carrecomplet(GEN x, GEN *pt)
 {  {
   long av;    gpmem_t av;
   GEN y;    GEN y;
   
   switch(signe(x))    switch(signe(x))
Line 427  carrecomplet(GEN x, GEN *pt)
Line 476  carrecomplet(GEN x, GEN *pt)
   if (!carremod((ulong)smodis(x, 64*63*65*11))) return 0;    if (!carremod((ulong)smodis(x, 64*63*65*11))) return 0;
   av=avma; y = racine(x);    av=avma; y = racine(x);
   if (!egalii(sqri(y),x)) { avma=av; return 0; }    if (!egalii(sqri(y),x)) { avma=av; return 0; }
   if (pt) { *pt = y; avma=(long)y; } else avma=av;    if (pt) { *pt = y; avma=(gpmem_t)y; } else avma=av;
   return 1;    return 1;
 }  }
   
 static GEN  static GEN
 polcarrecomplet(GEN x, GEN *pt)  polcarrecomplet(GEN x, GEN *pt)
 {  {
   long i,l,av,av2;    gpmem_t av,av2;
     long i,l;
   GEN y,a,b;    GEN y,a,b;
   
   if (!signe(x)) return gun;    if (!signe(x)) return gun;
Line 493  gcarrecomplet(GEN x, GEN *pt)
Line 543  gcarrecomplet(GEN x, GEN *pt)
 GEN  GEN
 gcarreparfait(GEN x)  gcarreparfait(GEN x)
 {  {
     gpmem_t av;
   GEN p1,a,p;    GEN p1,a,p;
   long tx=typ(x),l,i,av,v;    long tx=typ(x),l,i,v;
   
   switch(tx)    switch(tx)
   {    {
Line 506  gcarreparfait(GEN x)
Line 557  gcarreparfait(GEN x)
   
     case t_INTMOD:      case t_INTMOD:
     {      {
       GEN b, q, e;        GEN b, q;
       long w;        long w;
       a = (GEN)x[2]; if (!signe(a)) return gun;        a = (GEN)x[2]; if (!signe(a)) return gun;
       av = avma;        av = avma;
Line 533  gcarreparfait(GEN x)
Line 584  gcarreparfait(GEN x)
       if (i==0)        if (i==0)
       {        {
         GEN d = mppgcd(a,q);          GEN d = mppgcd(a,q);
         p1 = factor(d);          p = (GEN)factor(d)[1]; l = lg(p);
         p = (GEN)p1[1];  
         e = (GEN)p1[2]; l = lg(e);  
         for (i=1; i<l; i++)          for (i=1; i<l; i++)
         {          {
           v = pvaluation(a,(GEN)p[i],&p1);            v = pvaluation(a,(GEN)p[i],&p1);
           w = itos((GEN)e[i]);            w = pvaluation(q,(GEN)p[i], &q);
           if (v < w && (v&1 || kronecker(p1,(GEN)p[i]) == -1))            if (v < w && (v&1 || kronecker(p1,(GEN)p[i]) == -1))
             { avma = av; return gzero; }              { avma = av; return gzero; }
           q = diviiexact(q, gpowgs((GEN)p[i], w));  
         }          }
         if (kronecker(a,q) == -1) { avma = av; return gzero; }          if (kronecker(a,q) == -1) { avma = av; return gzero; }
       }        }
       /* kro(a,q) = 1, q odd: need to factor q */        /* kro(a,q) = 1, q odd: need to factor q */
       p1 = factor(q);        p = (GEN)factor(q)[1]; l = lg(p) - 1;
       p = (GEN)p1[1];  
       e = (GEN)p1[2]; l = lg(e) - 1;  
       /* kro(a,q) = 1, check all p|q but the last (product formula) */        /* kro(a,q) = 1, check all p|q but the last (product formula) */
       for (i=1; i<l; i++)        for (i=1; i<l; i++)
         if (kronecker(a,(GEN)p[i]) == -1) { avma = av; return gzero; }          if (kronecker(a,(GEN)p[i]) == -1) { avma = av; return gzero; }
Line 618  gkronecker(GEN x, GEN y)
Line 664  gkronecker(GEN x, GEN y)
 long  long
 kronecker(GEN x, GEN y)  kronecker(GEN x, GEN y)
 {  {
     const gpmem_t av = avma;
   GEN z;    GEN z;
   long av,r,s=1;    long r,s=1;
   
   av=avma;  
   switch (signe(y))    switch (signe(y))
   {    {
     case -1: y=negi(y); if (signe(x)<0) s= -1; break;      case -1: y=negi(y); if (signe(x)<0) s= -1; break;
Line 662  gkrogs(GEN x, long y)
Line 708  gkrogs(GEN x, long y)
 long  long
 krogs(GEN x, long y)  krogs(GEN x, long y)
 {  {
   long av,r,s=1,x1,z;    const gpmem_t av = avma;
     long r,s=1,x1,z;
   
   av=avma;  
   if (y<=0)    if (y<=0)
   {    {
     if (y) { y = -y; if (signe(x)<0) s = -1; }      if (y) { y = -y; if (signe(x)<0) s = -1; }
Line 698  krogs(GEN x, long y)
Line 744  krogs(GEN x, long y)
 long  long
 krosg(long s, GEN x)  krosg(long s, GEN x)
 {  {
   long av = avma, y = kronecker(stoi(s),x);    gpmem_t av = avma;
     long y = kronecker(stoi(s),x);
   avma = av; return y;    avma = av; return y;
 }  }
   
Line 747  hil0(GEN x, GEN y, GEN p)
Line 794  hil0(GEN x, GEN y, GEN p)
 long  long
 hil(GEN x, GEN y, GEN p)  hil(GEN x, GEN y, GEN p)
 {  {
   long a,b,av,tx,ty,z;    gpmem_t av;
     long a,b,tx,ty,z;
   GEN p1,p2,u,v;    GEN p1,p2,u,v;
   
   if (gcmp0(x) || gcmp0(y)) return 0;    if (gcmp0(x) || gcmp0(y)) return 0;
Line 779  hil(GEN x, GEN y, GEN p)
Line 827  hil(GEN x, GEN y, GEN p)
         case t_REAL:          case t_REAL:
           return (signe(x)<0 && signe(y)<0)? -1: 1;            return (signe(x)<0 && signe(y)<0)? -1: 1;
         case t_INTMOD:          case t_INTMOD:
           if (egalii(gdeux,(GEN)y[1])) err(hiler1);            p = (GEN)y[1];
           return hil(x,(GEN)y[2],(GEN)y[1]);            if (egalii(gdeux,p)) err(hiler1);
             return hil(x,(GEN)y[2],p);
         case t_FRAC: case t_FRACN:          case t_FRAC: case t_FRACN:
           p1=mulii((GEN)y[1],(GEN)y[2]); z=hil(x,p1,p);            p1=mulii((GEN)y[1],(GEN)y[2]); z=hil(x,p1,p);
           avma=av; return z;            avma=av; return z;
         case t_PADIC:          case t_PADIC:
           if (egalii(gdeux,(GEN)y[2]) && precp(y) <= 2) err(hiler1);            p = (GEN)y[2];
           p1 = odd(valp(y))? mulii((GEN)y[2],(GEN)y[4]): (GEN)y[4];            if (egalii(gdeux,p) && precp(y) <= 1) err(hiler1);
           z=hil(x,p1,(GEN)y[2]); avma=av; return z;            p1 = odd(valp(y))? mulii(p,(GEN)y[4]): (GEN)y[4];
             z=hil(x,p1,p); avma=av; return z;
       }        }
       break;        break;
   
Line 797  hil(GEN x, GEN y, GEN p)
Line 847  hil(GEN x, GEN y, GEN p)
       return signe(y[1])*signe(y[2]);        return signe(y[1])*signe(y[2]);
   
     case t_INTMOD:      case t_INTMOD:
       if (egalii(gdeux,(GEN)y[1])) err(hiler1);        p = (GEN)x[1];
         if (egalii(gdeux,p)) err(hiler1);
       switch(ty)        switch(ty)
       {        {
         case t_INTMOD:          case t_INTMOD:
           if (!egalii((GEN)x[1],(GEN)y[1])) break;            if (!egalii(p, (GEN)y[1])) break;
           return hil((GEN)x[2],(GEN)y[2],(GEN)x[1]);            return hil((GEN)x[2],(GEN)y[2],p);
         case t_FRAC: case t_FRACN:          case t_FRAC: case t_FRACN:
           return hil((GEN)x[2],y,(GEN)x[1]);            return hil((GEN)x[2],y,p);
         case t_PADIC:          case t_PADIC:
           if (!egalii((GEN)x[1],(GEN)y[2])) break;            if (!egalii(p, (GEN)y[2])) break;
           return hil((GEN)x[2],y,(GEN)x[1]);            return hil((GEN)x[2],y,p);
       }        }
       break;        break;
   
Line 824  hil(GEN x, GEN y, GEN p)
Line 875  hil(GEN x, GEN y, GEN p)
       break;        break;
   
     case t_PADIC:      case t_PADIC:
       if (ty != t_PADIC || !egalii((GEN)x[2],(GEN)y[2])) break;        p = (GEN)x[2];
       p1 = odd(valp(x))? mulii((GEN)x[2],(GEN)x[4]): (GEN)x[4];        if (ty != t_PADIC || !egalii(p,(GEN)y[2])) break;
       p2 = odd(valp(y))? mulii((GEN)y[2],(GEN)y[4]): (GEN)y[4];        if (egalii(p, gdeux) && (precp(x) <= 1 || precp(y) <= 1)) err(hiler1);
       z=hil(p1,p2,(GEN)x[2]); avma=av; return z;        p1 = odd(valp(x))? mulii(p,(GEN)x[4]): (GEN)x[4];
         p2 = odd(valp(y))? mulii(p,(GEN)y[4]): (GEN)y[4];
         z=hil(p1,p2,p); avma=av; return z;
   }    }
   err(talker,"forbidden or incompatible types in hil");    err(talker,"forbidden or incompatible types in hil");
   return 0; /* not reached */    return 0; /* not reached */
Line 844  hil(GEN x, GEN y, GEN p)
Line 897  hil(GEN x, GEN y, GEN p)
   
 #define sqrmod(x,p) (resii(sqri(x),p))  #define sqrmod(x,p) (resii(sqri(x),p))
   
 /* Assume p is prime. return NULL if (a,p) = -1 */  static GEN ffsqrtmod(GEN a, GEN p);
   
   /* Tonelli-Shanks. Assume p is prime and return NULL if (a,p) = -1. */
 GEN  GEN
 mpsqrtmod(GEN a, GEN p)  mpsqrtmod(GEN a, GEN p)
 {  {
   long av = avma, av1,lim,i,k,e;    gpmem_t av = avma, av1,lim;
     long i,k,e;
   GEN p1,q,v,y,w,m;    GEN p1,q,v,y,w,m;
   
   if (typ(a) != t_INT || typ(p) != t_INT) err(arither1);    if (typ(a) != t_INT || typ(p) != t_INT) err(arither1);
   if (signe(p) <= 0 || is_pm1(p)) err(talker,"not a prime in mpsqrtmod");    if (signe(p) <= 0 || is_pm1(p)) err(talker,"not a prime in mpsqrtmod");
   p1 = addsi(-1,p); e = vali(p1);    p1 = addsi(-1,p); e = vali(p1);
   
     /* If `e' is "too big", use Cipolla algorithm ! [GTL] */
     if (e*(e-1) > 20 + 8 * bit_accuracy(lgefint(p)))
     {
       v = ffsqrtmod(a,p);
       if (!v) { avma = av; return NULL; }
       return gerepileuptoint(av,v);
     }
   
   if (e == 0) /* p = 2 */    if (e == 0) /* p = 2 */
   {    {
     avma = av;      avma = av;
Line 911  mpsqrtmod(GEN a, GEN p)
Line 976  mpsqrtmod(GEN a, GEN p)
   return gerepileuptoint(av, v);    return gerepileuptoint(av, v);
 }  }
   
   /* Cipolla's algorithm is better when e = v_2(p-1) is "too big".
    * Otherwise, is a constant times worse than the above one.
    * For p = 3 (mod 4), is about 3 times worse, and in average
    * is about 2 or 2.5 times worse.
    *
    * But try both algorithms for e.g. S(n)=(2^n+3)^2-8 with
    * n = 750, 771, 779, 790, 874, 1176, 1728, 2604, etc.
    *
    * If X^2 = t^2 - a  is not a square in F_p, then
    *
    *   (t+X)^(p+1) = (t+X)(t-X) = a
    *
    * so we get sqrt(a) in F_p^2 by
    *
    *   sqrt(a) = (t+X)^((p+1)/2)
    *
    * If (a|p)=1, then sqrt(a) is in F_p.
    *
    * cf: LNCS 2286, pp 430-434 (2002)  [Gonzalo Tornaria] */
   static GEN
   ffsqrtmod(GEN a, GEN p)
   {
     gpmem_t av = avma, av1, lim;
     long e, t, man, k, nb;
     GEN q, n, u, v, d, d2, b, u2, v2, qptr;
   
     if (kronecker(a, p) < 0) return NULL;
   
     av1 = avma;
     for(t=1; ; t++)
     {
       n = subsi(t*t, a);
       if (kronecker(n, p) < 0) break;
       avma = av1;
     }
   
     u = stoi(t); v = gun; /* u+vX = t+X */
     e = vali(addsi(-1,p)); q = shifti(p, -e);
     /* p = 2^e q + 1  and  (p+1)/2 = 2^(e-1)q + 1 */
   
     /* Compute u+vX := (t+X)^q */
     av1 = avma; lim = stack_lim(av1, 1);
     qptr = q+2; man = *qptr;
     k = 1+bfffo(man); man<<=k; k=BITS_IN_LONG-k;
     for (nb=lgefint(q)-2;nb; nb--)
     {
       for (; k; man<<=1, k--)
       {
         if (man < 0)
         { /* u+vX := (u+vX)^2 (t+X) */
           d = addii(u, mulsi(t,v));
           d2 = sqri(d);
           b = resii(mulii(a,v), p);
           u = modii(subii(mulsi(t,d2), mulii(b,addii(u,d))), p);
           v = modii(subii(d2, mulii(b,v)), p);
         }
         else
         { /* u+vX := (u+vX)^2 */
           u2 = sqri(u);
           v2 = sqri(v);
           v = modii(subii(sqri(addii(v,u)), addii(u2,v2)), p);
           u = modii(addii(u2, mulii(v2,n)), p);
         }
       }
   
       if (low_stack(lim, stack_lim(av1, 1)))
       {
          GEN *gptr[2]; gptr[0]=&u; gptr[1]=&v;
          if (DEBUGMEM>1) err(warnmem, "ffsqrtmod");
          gerepilemany(av1,gptr,2);
       }
   
       man = *++qptr; k = BITS_IN_LONG;
     }
   
     while (--e)
     { /* u+vX := (u+vX)^2 */
       u2 = sqri(u);
       v2 = sqri(v);
       v = modii(subii(sqri(addii(v,u)), addii(u2,v2)), p);
       u = modii(addii(u2, mulii(v2,n)), p);
   
       if (low_stack(lim, stack_lim(av1, 1)))
       {
          GEN *gptr[2]; gptr[0]=&u; gptr[1]=&v;
          if (DEBUGMEM>1) err(warnmem, "ffsqrtmod");
          gerepilemany(av1,gptr,2);
       }
     }
   
     /* Now u+vX = (t+X)^(2^(e-1)q); thus
      *
      *   (u+vX)(t+X) = sqrt(a) + 0 X
      *
      * Whence,
      *
      *   sqrt(a) = (u+vt)t - v*a
      *   0       = (u+vt)
      *
      * Thus a square root is v*a */
   
     v = modii(mulii(v,a), p);
   
     u = subii(p,v); if (cmpii(v,u) > 0) v = u;
     return gerepileuptoint(av,v);
   }
   
 /*******************************************************************/  /*******************************************************************/
 /*                                                                 */  /*                                                                 */
 /*                       n-th ROOT MODULO p                        */  /*                       n-th ROOT MODULO p                        */
Line 930  mpsqrtmod(GEN a, GEN p)
Line 1102  mpsqrtmod(GEN a, GEN p)
 static GEN  static GEN
 mplgenmod(GEN l, long e, GEN r,GEN p,GEN *zeta)  mplgenmod(GEN l, long e, GEN r,GEN p,GEN *zeta)
 {  {
   ulong av1;    const gpmem_t av1 = avma;
   GEN m,m1;    GEN m,m1;
   long k,i;    long k,i;
   av1 = avma;  
   for (k=1; ; k++)    for (k=1; ; k++)
   {    {
     m1 = m = powmodulo(stoi(k+1),r,p);      m1 = m = powmodulo(stoi(k+1),r,p);
       if (gcmp1(m)) {avma=av1; continue;}
     for (i=1; i<e; i++)      for (i=1; i<e; i++)
       if (gcmp1(m=powmodulo(m,l,p))) break;        if (gcmp1(m=powmodulo(m,l,p))) break;
     if (i==e) break;      if (i==e) break;
Line 960  GEN Fp_shanks(GEN x,GEN g0,GEN p, GEN q);
Line 1132  GEN Fp_shanks(GEN x,GEN g0,GEN p, GEN q);
 static GEN  static GEN
 mpsqrtlmod(GEN a, GEN l, GEN p, GEN q,long e, GEN r, GEN y, GEN m)  mpsqrtlmod(GEN a, GEN l, GEN p, GEN q,long e, GEN r, GEN y, GEN m)
 {  {
   ulong av = avma, tetpil,lim;    gpmem_t av = avma, tetpil,lim;
   long k;    long k;
   GEN p1,u1,u2,v,w,z,dl;    GEN p1,u1,u2,v,w,z,dl;
   /* y contient un generateur de (Z/pZ)^* eleve a la puis (p-1)/(l^e) */    /* y contient un generateur de (Z/pZ)^* eleve a la puis (p-1)/(l^e) */
   bezout(r,l,&u1,&u2);    (void)bezout(r,l,&u1,&u2);
   v=powmodulo(a,u2,p);    v=powmodulo(a,u2,p);
   w=powmodulo(a,modii(mulii(negi(u1),r),q),p);    w=powmodulo(a,modii(mulii(negi(u1),r),q),p);
   lim = stack_lim(av,1);    lim = stack_lim(av,1);
Line 1015  If a=0 ,return 0 and if zetan is not NULL zetan is set
Line 1187  If a=0 ,return 0 and if zetan is not NULL zetan is set
 GEN  GEN
 mpsqrtnmod(GEN a, GEN n, GEN p, GEN *zetan)  mpsqrtnmod(GEN a, GEN n, GEN p, GEN *zetan)
 {  {
   ulong ltop=avma,lbot=0,av1,lim;    gpmem_t ltop=avma,lbot=0,av1,lim;
   GEN m,u1,u2;    GEN m,u1,u2;
   GEN q,z;    GEN q,z;
   GEN *gptr[2];    GEN *gptr[2];
Line 1183  ulong mppgcd_resiu(GEN y, ulong x);
Line 1355  ulong mppgcd_resiu(GEN y, ulong x);
 GEN  GEN
 mppgcd(GEN a, GEN b)  mppgcd(GEN a, GEN b)
 {  {
   long av,v,w;    long v, w;
     gpmem_t av;
   GEN t, p1;    GEN t, p1;
   
   if (typ(a) != t_INT || typ(b) != t_INT) err(arither1);    if (typ(a) != t_INT || typ(b) != t_INT) err(arither1);
Line 1204  mppgcd(GEN a, GEN b)
Line 1377  mppgcd(GEN a, GEN b)
   }    }
   
   /* larger than gcd: "avma=av" gerepile (erasing t) is valid */    /* larger than gcd: "avma=av" gerepile (erasing t) is valid */
   av = avma; new_chunk(lgefint(b)); /* HACK */    av = avma; (void)new_chunk(lgefint(b)); /* HACK */
   t = resii(a,b);    t = resii(a,b);
   if (!signe(t)) { avma=av; return absi(b); }    if (!signe(t)) { avma=av; return absi(b); }
   
Line 1237  mppgcd(GEN a, GEN b)
Line 1410  mppgcd(GEN a, GEN b)
     }      }
   }    }
   {    {
     long r[] = {evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3), 0};      long r[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 0};
     r[2] = (long) ugcd((ulong)b[2], (ulong)a[2]);      r[2] = (long) ugcd((ulong)b[2], (ulong)a[2]);
     avma = av; return shifti(r,v);      avma = av; return shifti(r,v);
   }    }
Line 1246  mppgcd(GEN a, GEN b)
Line 1419  mppgcd(GEN a, GEN b)
 GEN  GEN
 mpppcm(GEN x, GEN y)  mpppcm(GEN x, GEN y)
 {  {
   ulong av;    gpmem_t av;
   GEN p1,p2;    GEN p1,p2;
   if (typ(x) != t_INT || typ(y) != t_INT) err(arither1);    if (typ(x) != t_INT || typ(y) != t_INT) err(arither1);
   if (!signe(x)) return gzero;    if (!signe(x)) return gzero;
Line 1292  GEN chinese(GEN x, GEN y)
Line 1465  GEN chinese(GEN x, GEN y)
 GEN  GEN
 chinois(GEN x, GEN y)  chinois(GEN x, GEN y)
 {  {
   ulong av,tetpil, tx = typ(x);    gpmem_t av,tetpil;
   long i,lx,vx;    long i,lx,vx, tx = typ(x);
   GEN z,p1,p2,d,u,v;    GEN z,p1,p2,d,u,v;
   
   if (gegal(x,y)) return gcopy(x);    if (gegal(x,y)) return gcopy(x);
Line 1339  chinois(GEN x, GEN y)
Line 1512  chinois(GEN x, GEN y)
 GEN  GEN
 chinois_int_coprime(GEN x2, GEN y2, GEN x1, GEN y1, GEN z1)  chinois_int_coprime(GEN x2, GEN y2, GEN x1, GEN y1, GEN z1)
 {  {
   ulong av = avma;    gpmem_t av = avma;
   GEN ax,p1;    GEN ax,p1;
   (void)new_chunk((lgefint(z1)<<1)+lgefint(x1)+lgefint(y1)); /* HACK */    (void)new_chunk((lgefint(z1)<<1)+lgefint(x1)+lgefint(y1)); /* HACK */
   ax = mulii(mpinvmod(x1,y1), x1);    ax = mulii(mpinvmod(x1,y1), x1);
Line 1364  mpinvmod(GEN a, GEN m)
Line 1537  mpinvmod(GEN a, GEN m)
   
 /*********************************************************************/  /*********************************************************************/
 /**                                                                 **/  /**                                                                 **/
 /**                   POWERING MODULO (A^N mod M)                   **/  /**                    MODULAR EXPONENTIATION                       **/
 /**                                                                 **/  /**                                                                 **/
 /*********************************************************************/  /*********************************************************************/
 GEN resiimul(GEN x, GEN y);  extern ulong invrev(ulong b);
 GEN resmod2n(GEN x, long y);  extern GEN resiimul(GEN x, GEN y);
 GEN _resii(GEN x, GEN y) { return resii(x,y); }  
   
   static GEN _resii(GEN x, GEN y) { return resii(x,y); }
   
   /* Montgomery reduction */
   
   typedef struct {
     GEN N;
     ulong inv; /* inv = -N^(-1) mod B, */
   } montdata;
   
   void
   init_montdata(GEN N, montdata *s)
   {
     s->N = N;
     s->inv = (ulong) -invrev(modBIL(N));
   }
   
 GEN  GEN
 init_remainder(GEN M)  init_remainder(GEN M)
 {  {
   GEN sM = cgetg(3, t_VEC);    GEN sM = cgetg(3, t_VEC);
   GEN Mr = cgetr(lgefint(M) + 1);  
   affir(M, Mr);  
   sM[1] = (long)M;    sM[1] = (long)M;
   sM[2] = (long)linv(Mr);    sM[2] = (long)linv( itor(M, lgefint(M) + 1) ); /* 1. / M */
   return sM;    return sM;
 }  }
   
 /* optimal on UltraSPARC */  /* optimal on UltraSPARC */
 static long RESIIMUL_LIMIT = 150;  static long RESIIMUL_LIMIT   = 150;
   static long MONTGOMERY_LIMIT = 32;
   
 void  void
 setresiilimit(long n) { RESIIMUL_LIMIT = n; }  setresiilimit(long n) { RESIIMUL_LIMIT = n; }
   void
   setmontgomerylimit(long n) { MONTGOMERY_LIMIT = n; }
   
   typedef struct {
     GEN N;
     GEN (*res)(GEN,GEN);
     GEN (*mulred)(GEN,GEN,GEN);
   } muldata;
   
   /* reduction for multiplication by 2 */
   static GEN
   _redsub(GEN x, GEN N)
   {
     return (cmpii(x,N) >= 0)? subii(x,N): x;
   }
   /* Montgomery reduction */
   extern GEN red_montgomery(GEN T, GEN N, ulong inv);
   static GEN
   montred(GEN x, GEN N)
   {
     return red_montgomery(x, ((montdata*)N)->N, ((montdata*)N)->inv);
   }
   /* 2x mod N */
   static GEN
   _muli2red(GEN x, GEN y/* ignored */, GEN N) {
     (void)y; return _redsub(shifti(x,1), N);
   }
   static GEN
   _muli2montred(GEN x, GEN y/* ignored */, GEN N) {
     GEN n = ((montdata*)N)->N;
     GEN z = _muli2red(x,y, n);
     long l = lgefint(n);
     while (lgefint(z) > l) z = subii(z,n);
     return z;
   }
   static GEN
   _muli2invred(GEN x, GEN y/* ignored */, GEN N) {
     return _muli2red(x,y, (GEN)N[1]);
   }
   /* xy mod N */
   static GEN
   _muliired(GEN x, GEN y, GEN N) { return resii(mulii(x,y), N); }
   static GEN
   _muliimontred(GEN x, GEN y, GEN N) { return montred(mulii(x,y), N); }
   static GEN
   _muliiinvred(GEN x, GEN y, GEN N) { return resiimul(mulii(x,y), N); }
   
   static GEN
   _mul(void *data, GEN x, GEN y)
   {
     muldata *D = (muldata *)data;
     return D->mulred(x,y,D->N);
   }
   static GEN
   _sqr(void *data, GEN x)
   {
     muldata *D = (muldata *)data;
     return D->res(sqri(x), D->N);
   }
   
   /* A^k mod N */
 GEN  GEN
 powmodulo(GEN A, GEN N, GEN M)  powmodulo(GEN A, GEN k, GEN N)
 {  {
     gpmem_t av = avma;
     long t,s, lN;
     int base_is_2, use_montgomery;
   GEN y;    GEN y;
   long av = avma,av1,lim,man,k,nb,s, *p;    muldata  D;
   GEN (*mul)(GEN,GEN) = mulii;    montdata S;
   GEN (*res)(GEN,GEN) = _resii;  
   
   if (typ(A) != t_INT || typ(N) != t_INT || typ(M) != t_INT) err(arither1);    if (typ(A) != t_INT || typ(k) != t_INT || typ(N) != t_INT) err(arither1);
   s = signe(N);    s = signe(k);
   if (!s)    if (!s)
   {    {
     k = signe(resii(A,M)); avma=av;      t = signe(resii(A,N)); avma = av;
     return k? gun: gzero;      return t? gun: gzero;
   }    }
   if (s < 0) { A = mpinvmod(A,M); N = absi(N); }    if (s < 0) y = mpinvmod(A,N);
   else    else
   {    {
     A = modii(A,M);      y = modii(A,N);
     if (!signe(A)) { avma=av; return gzero; }      if (!signe(y)) { avma = av; return gzero; }
   }    }
   y = A;  
   if (lgefint(y)==3) switch(y[2])    base_is_2 = 0;
     if (lgefint(y) == 3) switch(y[2])
   {    {
     case 1: /* y = 1 */      case 1: avma = av; return gun;
       avma=av; return gun;      case 2: base_is_2 = 1; break;
     case 2: /* y = 2, use shifti not mulii */  
     mul = (GEN (*)(GEN,GEN))shifti; A = (GEN)1L;  
   }    }
   
   /* TODO: Move this out of here and use for general modular computations */    /* TODO: Move this out of here and use for general modular computations */
   if ((k = vali(M)) && k == expi(M))    lN = lgefint(N);
   { /* M is a power of 2 */    use_montgomery = mod2(N) && lN < MONTGOMERY_LIMIT;
     M = (GEN)k;    if (use_montgomery)
     res = (GEN(*)(GEN,GEN))resmod2n;    {
       init_montdata(N, &S);
       y = resii(shifti(y, bit_accuracy(lN)), N);
       D.mulred = base_is_2? &_muli2montred: &_muliimontred;
       D.res = &montred;
       D.N = (GEN)&S;
   }    }
   else if (lgefint(M) > RESIIMUL_LIMIT && (lgefint(N) > 3 || N[2] > 4))    else if (lN > RESIIMUL_LIMIT
   { /* compute x % M using multiplication by 1./M */         && (lgefint(k) > 3 || (((double)k[2])*expi(A)) > 2 + expi(N)))
     M = init_remainder(M);    {
     res = (GEN(*)(GEN,GEN))resiimul;      D.mulred = base_is_2? &_muli2invred: &_muliiinvred;
       D.res = &resiimul;
       D.N = init_remainder(N);
   }    }
     else
   p = N+2; man = *p; /* see puissii */  
   k = 1+bfffo(man); man<<=k; k = BITS_IN_LONG-k;  
   av1=avma; lim=stack_lim(av1,1);  
   for (nb=lgefint(N)-2;;)  
   {    {
     for (; k; man<<=1,k--)      D.mulred = base_is_2? &_muli2red: &_muliired;
     {      D.res = &_resii;
       y = res(sqri(y), M);      D.N = N;
       if (man < 0) y = res(mul(y,A), M);  
       if (low_stack(lim, stack_lim(av1,1)))  
       {  
         if (DEBUGMEM>1) err(warnmem,"powmodulo");  
         y = gerepileuptoint(av1,y);  
       }  
     }  
     if (--nb == 0) break;  
     man = *++p, k = BITS_IN_LONG;  
   }    }
   return gerepileupto(av,y);    y = leftright_pow(y, k, (void*)&D, &_sqr, &_mul);
     if (use_montgomery)
     {
       y = montred(y, (GEN)&S);
       if (cmpii(y,N) >= 0) y = subii(y,N);
     }
     return gerepileuptoint(av,y);
 }  }
   
 /*********************************************************************/  /*********************************************************************/
Line 1457  powmodulo(GEN A, GEN N, GEN M)
Line 1705  powmodulo(GEN A, GEN N, GEN M)
 /**                                                                 **/  /**                                                                 **/
 /*********************************************************************/  /*********************************************************************/
 GEN  GEN
 gnextprime(GEN n)  gnextprime(GEN n) { return garith_proto(nextprime,n,0); }
 {  
   return garith_proto(nextprime,n,0);  
 }  
   
 GEN  GEN
 gprecprime(GEN n)  gprecprime(GEN n) { return garith_proto(precprime,n,0); }
 {  
   return garith_proto(precprime,n,0);  
 }  
   
 GEN  GEN
 gisprime(GEN x, long flag)  gisprime(GEN x, long flag)
 {  {
   switch (flag)    switch (flag)
   {    {
   case 0:      case 0: return arith_proto(isprime,x,1);
     return arith_proto(isprime,x,1);      case 1: return garith_proto2gs(plisprime,x,1);
   case 1:      case 2: return arith_proto(isprimeAPRCL,x,1);
     return garith_proto2gs(plisprime,x,0);  
   case 2:  
     return garith_proto2gs(plisprime,x,1);  
   }    }
   err(flagerr,"gisprime");    err(flagerr,"gisprime");
   return 0;    return 0;
 }  }
   
 long  long
 isprime(GEN x)  isprimeSelfridge(GEN x) { return (plisprime(x,0)==gun); }
   
   /* assume x BSW pseudoprime. Check whether it's small enough to be certified
    * prime */
   int
   BSW_isprime_small(GEN x)
 {  {
   return millerrabin(x,10);    long l = lgefint(x);
     if (l < 4) return 1;
     if (l == 4)
     {
       gpmem_t av = avma;
       long t = cmpii(x, u2toi(0x918UL, 0x4e72a000UL)); /* 10^13 */
       avma = av;
       if (t < 0) return 1;
     }
     return 0;
 }  }
   
 GEN  /* assume x a BSW pseudoprime */
 gispsp(GEN x)  int
   BSW_isprime(GEN x)
 {  {
   return arith_proto(ispsp,x,1);    gpmem_t av = avma;
     long l, res;
     GEN F, p;
   
     if (BSW_isprime_small(x)) return 1;
     F = (GEN)auxdecomp(subis(x,1), 0)[1];
     l = lg(F); p = (GEN)F[l-1];
     if (BSW_psp(p))
     { /* smooth */
       GEN z = cgetg(3, t_VEC);
       z[1] = (long)x;
       z[2] = (long)F; res = isprimeSelfridge(z);
     }
     else
       res = isprimeAPRCL(x);
     avma = av; return res;
 }  }
   
 long  long
 ispsp(GEN x)  isprime(GEN x)
 {  {
   return millerrabin(x,1);    return BSW_psp(x) && BSW_isprime(x);
 }  }
   
 GEN  GEN
 gmillerrabin(GEN n, long k)  gispseudoprime(GEN x, long flag)
 {  {
   return arith_proto2gs(millerrabin,n,k);    if (flag == 0) return arith_proto(BSW_psp,x,1);
     return gmillerrabin(x, flag);
 }  }
   
   long
   ispseudoprime(GEN x, long flag)
   {
     if (flag == 0) return BSW_psp(x);
     return millerrabin(x, flag);
   }
   
   GEN
   gispsp(GEN x) { return arith_proto(ispsp,x,1); }
   
   long
   ispsp(GEN x) { return millerrabin(x,1); }
   
   GEN
   gmillerrabin(GEN n, long k) { return arith_proto2gs(millerrabin,n,k); }
   
 /*********************************************************************/  /*********************************************************************/
 /**                                                                 **/  /**                                                                 **/
 /**                    FUNDAMENTAL DISCRIMINANTS                    **/  /**                    FUNDAMENTAL DISCRIMINANTS                    **/
 /**                                                                 **/  /**                                                                 **/
 /*********************************************************************/  /*********************************************************************/
   
 /* gissquarefree, issquarefree moved into arith2.c with the other  
    basic arithmetic functions (issquarefree is the square of the  
    Moebius mu function, after all...) --GN */  
   
 GEN  GEN
 gisfundamental(GEN x)  gisfundamental(GEN x) { return arith_proto(isfundamental,x,1); }
 {  
   return arith_proto(isfundamental,x,1);  
 }  
   
 long  long
 isfundamental(GEN x)  isfundamental(GEN x)
 {  {
   long av,r;    long r;
   GEN p1;    GEN p1;
   
   if (gcmp0(x)) return 0;    if (gcmp0(x)) return 0;
   r=mod4(x);    r=mod4(x);
   if (!r)    if (!r)
   {    {
     av=avma; p1=shifti(x,-2);      const gpmem_t av = avma;
       p1=shifti(x,-2);
     r=mod4(p1); if (!r) return 0;      r=mod4(p1); if (!r) return 0;
     if (signe(x)<0) r=4-r;      if (signe(x)<0) r=4-r;
     r = (r==1)? 0: issquarefree(p1);      r = (r==1)? 0: issquarefree(p1);
Line 1547  isfundamental(GEN x)
Line 1826  isfundamental(GEN x)
 GEN  GEN
 quaddisc(GEN x)  quaddisc(GEN x)
 {  {
   long av=avma,tetpil=avma,i,r,tx=typ(x);    const gpmem_t av = avma;
     long i,r,tx=typ(x);
   GEN p1,p2,f,s;    GEN p1,p2,f,s;
   
   if (tx != t_INT && ! is_frac_t(tx)) err(arither1);    if (tx != t_INT && ! is_frac_t(tx)) err(arither1);
   f=factor(x); p1=(GEN)f[1]; p2=(GEN)f[2];    f=factor(x); p1=(GEN)f[1]; p2=(GEN)f[2];
   s=gun;    s = gun;
   for (i=1; i<lg(p1); i++)    for (i=1; i<lg(p1); i++)
     if ( odd(mael(p2,i,2)) ) { tetpil=avma; s=gmul(s,(GEN)p1[i]); }      if (odd(mael(p2,i,2))) s = gmul(s,(GEN)p1[i]);
   r=mod4(s); if (gsigne(x)<0) r=4-r;    r=mod4(s); if (gsigne(x)<0) r=4-r;
   if (r>1) { tetpil=avma; s=shifti(s,2); }    if (r>1) s = shifti(s,2);
   return gerepile(av,tetpil,s);    return gerepileuptoint(av, s);
 }  }
   
 /*********************************************************************/  /*********************************************************************/
Line 1565  quaddisc(GEN x)
Line 1845  quaddisc(GEN x)
 /**                              FACTORIAL                          **/  /**                              FACTORIAL                          **/
 /**                                                                 **/  /**                                                                 **/
 /*********************************************************************/  /*********************************************************************/
 GEN muluu(ulong x, ulong y);  
   
 GEN  GEN
 mpfact(long n)  mpfact(long n)
 {  {
   long lx,k,l, av = avma;    const gpmem_t av = avma;
     long lx,k,l;
   GEN x;    GEN x;
   
   if (n<2)    if (n<2)
Line 1598  mpfact(long n)
Line 1877  mpfact(long n)
 GEN  GEN
 mpfactr(long n, long prec)  mpfactr(long n, long prec)
 {  {
   long av = avma, lim,k;    gpmem_t av = avma, lim;
   GEN f = cgetr(prec);    long k;
     GEN f = realun(prec);
   
   affsr(1,f);  
   if (n<2)    if (n<2)
   {    {
     if (n<0) err(facter);      if (n<0) err(facter);
Line 1629  mpfactr(long n, long prec)
Line 1908  mpfactr(long n, long prec)
 void  void
 lucas(long n, GEN *ln, GEN *ln1)  lucas(long n, GEN *ln, GEN *ln1)
 {  {
   long taille,av;    gpmem_t av;
     long taille;
   GEN z,t;    GEN z,t;
   
   if (!n) { *ln=stoi(2); *ln1=stoi(1); return; }    if (!n) { *ln = stoi(2); *ln1 = stoi(1); return; }
   
   taille=(long)(pariC3*(1+labs(n))+3);    taille = 3 + (long)(pariC3 * (1+labs(n)));
   *ln=cgeti(taille); *ln1=cgeti(taille);    *ln = cgeti(taille);
   av=avma; lucas(n/2,&z,&t);    *ln1= cgeti(taille);
     av = avma; lucas(n/2, &z, &t);
   switch(n % 4)    switch(n % 4)
   {    {
     case -3:      case -3:
       addsiz(2,sqri(z),*ln1);        addsiz(2,sqri(z), *ln1);
       subiiz(addsi(1,mulii(z,t)),*ln1,*ln); break;        subiiz(addsi(1,mulii(z,t)),*ln1, *ln); break;
     case -2:  
       addsiz(2,sqri(z),*ln); addsiz(1,mulii(z,t),*ln1); break;  
     case -1:      case -1:
       subisz(sqri(z),2,*ln1);        subisz(sqri(z),2, *ln1);
       subiiz(subis(mulii(z,t),1),*ln1,*ln); break;        subiiz(subis(mulii(z,t),1),*ln1, *ln); break;
     case  0: subisz(sqri(z),2,*ln); subisz(mulii(z,t),1,*ln1); break;      case  0: subisz(sqri(z),2,    *ln); subisz(mulii(z,t),1, *ln1); break;
     case  1: subisz(mulii(z,t),1,*ln); addsiz(2,sqri(t),*ln1); break;      case  1: subisz(mulii(z,t),1, *ln); addsiz(2,sqri(t),    *ln1); break;
     case  2: addsiz(2,sqri(z),*ln); addsiz(1,mulii(z,t),*ln1); break;      case -2:
     case  3: addsiz(1,mulii(z,t),*ln); subisz(sqri(t),2,*ln1);      case  2: addsiz(2,sqri(z),    *ln); addsiz(1,mulii(z,t), *ln1); break;
       case  3: addsiz(1,mulii(z,t), *ln); subisz(sqri(t),2,    *ln1);
   }    }
   avma=av;    avma = av;
 }  }
   
 GEN  GEN
 fibo(long n)  fibo(long n)
 {  {
   long av = avma;    gpmem_t av = avma;
   GEN ln,ln1;    GEN ln,ln1;
   
   lucas(n-1,&ln,&ln1);    lucas(n-1,&ln,&ln1);
Line 1670  fibo(long n)
Line 1950  fibo(long n)
 /*                      CONTINUED FRACTIONS                        */  /*                      CONTINUED FRACTIONS                        */
 /*                                                                 */  /*                                                                 */
 /*******************************************************************/  /*******************************************************************/
 static GEN sfcont2(GEN b, GEN x, long k);  static GEN
   icopy_lg(GEN x, long l)
 GEN  
 gcf(GEN x)  
 {  {
   return sfcont(x,x,0);    long lx = lgefint(x);
     GEN y;
   
     if (lx >= l) return icopy(x);
     y = cgeti(l); affii(x, y); return y;
 }  }
   
 GEN  /* if y != NULL, stop as soon as partial quotients differ from y */
 gcf2(GEN b, GEN x)  static GEN
   Qsfcont(GEN x, GEN y, long k)
 {  {
   return contfrac0(x,b,0);    GEN  z, p1, p2, p3;
 }    long i, l, ly;
   
 GEN    ly = lgefint(x[2]);
 gboundcf(GEN x, long k)    l = 3 + (long) ((ly-2) / pariC3);
 {    if (k > 0 && ++k > 0 && l > k) l = k; /* beware overflow */
   return sfcont(x,x,k);    if ((ulong)l > LGBITS) l = LGBITS;
 }  
   
 GEN    p1 = (GEN)x[1];
 contfrac0(GEN x, GEN b, long flag)    p2 = (GEN)x[2];
 {    z = cgetg(l,t_VEC);
   long lb,tb,i;    l--;
   GEN y;    if (y) {
       gpmem_t av = avma;
   if (!b || gcmp0(b)) return sfcont(x,x,flag);      if (l >= lg(y)) l = lg(y)-1;
   tb = typ(b);      for (i = 1; i <= l; i++)
   if (tb == t_INT) return sfcont(x,x,itos(b));      {
   if (! is_matvec_t(tb)) err(typeer,"contfrac0");        z[i] = y[i];
         p3 = p2; if (!gcmp1((GEN)z[i])) p3 = mulii((GEN)z[i], p2);
   lb=lg(b); if (lb==1) return cgetg(1,t_VEC);        p3 = subii(p1, p3);
   if (tb != t_MAT) return sfcont2(b,x,flag);        if (signe(p3) < 0)
   if (lg(b[1])==1) return sfcont(x,x,flag);        { /* partial quotient too large */
   y = (GEN) gpmalloc(lb*sizeof(long));          p3 = addii(p3, p2);
   for (i=1; i<lb; i++) y[i]=coeff(b,1,i);          if (signe(p3) >= 0) i++; /* by 1 */
   x = sfcont2(y,x,flag); free(y); return x;          break;
         }
         if (cmpii(p3, p2) >= 0)
         { /* partial quotient too small */
           p3 = subii(p3, p2);
           if (cmpii(p3, p2) < 0) i++; /* by 1 */
           break;
         }
         if ((i & 0xff) == 0) gerepileall(av, 2, &p2, &p3);
         p1 = p2; p2 = p3;
       }
     } else {
       p1 = icopy_lg(p1, ly);
       p2 = icopy(p2);
       for (i = 1; i <= l; i++)
       {
         z[i] = (long)truedvmdii(p1,p2,&p3);
         if (p3 == gzero) { i++; break; }
         affii(p3, p1); cgiv(p3); p3 = p1;
         p1 = p2; p2 = p3;
       }
     }
     i--;
     if (i > 2 && gcmp1((GEN)z[i]))
     {
       cgiv((GEN)z[i]); --i;
       addsiz(1,(GEN)z[i], (GEN)z[i]);
     }
     setlg(z,i+1); return z;
 }  }
   
 GEN  static GEN
 sfcont(GEN x, GEN x1, long k)  sfcont(GEN x, long k)
 {  {
   long  av,lx=lg(x),tx=typ(x),e,i,j,l,l1,lx1,tetpil,f;    gpmem_t av;
   GEN   y,p1,p2,p3,p4,yp;    long lx,tx=typ(x),e,i,l;
     GEN  y,p1,p2,p3;
   
   if (is_scalar_t(tx))    if (is_scalar_t(tx))
   {    {
     if (gcmp0(x)) { y=cgetg(2,t_VEC); y[1]=zero; return y; }      if (gcmp0(x)) return _vec(gzero);
     switch(tx)      switch(tx)
     {      {
       case t_INT:        case t_INT: y = cgetg(2,t_VEC); y[1] = (long)icopy(x); return y;
         y=cgetg(2,t_VEC); y[1]=lcopy(x); return y;  
       case t_REAL:        case t_REAL:
         l=avma; p1=cgetg(3,t_FRACN);          av = avma; lx = lg(x);
         p2 = rcopy(x); settyp(p2,t_INT); setlgefint(p2,lx);          p2 = rcopy(x); settyp(p2,t_INT); setlgefint(p2,lx);
         p1[1] = (long) p2; e = bit_accuracy(lx)-1-expo(x);          e = bit_accuracy(lx)-1-expo(x);
         if (e<0) err(talker,"integral part not significant in scfont");          if (e < 0) err(talker,"integral part not significant in scfont");
           p1 = cgetg(3, t_FRACN);
           p1[1] = (long)p2;
           p1[2] = lshifti(gun, e);
   
           p3 = cgetg(3, t_FRACN);
           p3[1] = laddsi(signe(x), p2);
           p3[2] = p1[2];
           p1 = Qsfcont(p1,NULL,k);
           return gerepilecopy(av, Qsfcont(p3,p1,k));
   
         l1 = (e>>TWOPOTBITS_IN_LONG)+3; p2=cgeti(l1);  
         p2[1]= evalsigne(1)|evallgefint(l1);  
         p2[2]= (1L<<(e&(BITS_IN_LONG-1)));  
         for (i=3; i<l1; i++) p2[i]=0;  
         p1[2]=(long) p2; p3=cgetg(3,t_FRACN);  
         p3[2]=lcopy(p2);  
         p3[1]=laddsi(signe(x),(GEN)p1[1]);  
         p1=sfcont(p1,p1,k); tetpil=avma;  
         return gerepile(l,tetpil,sfcont(p3,p1,k));  
   
       case t_FRAC: case t_FRACN:        case t_FRAC: case t_FRACN:
         l=avma; lx1=lgefint(x[2]);          av = avma;
         l1 = (long) ((double) BYTES_IN_LONG/4.0*46.093443*(lx1-2)+3);          return gerepileupto(av, Qsfcont(x, NULL, k));
         if (k>0 && ++k > 0 && l1 > k) l1 = k; /* beware overflow */  
         if ((ulong)l1>LGBITS) l1=LGBITS;  
         if (lgefint(x[1]) >= lx1) p1=icopy((GEN)x[1]);  
         else affii((GEN)x[1],p1=cgeti(lx1));  
         p2=icopy((GEN)x[2]); lx1=lg(x1);  
         y=cgetg(l1,t_VEC); f=(x!=x1); i=0;  
         while (!gcmp0(p2) && i<=l1-2)  
         {  
           i++; y[i]=ldvmdii(p1,p2,&p3);  
           if (signe(p3)>=0) affii(p3,p1);  
           else  
           {  
             p4=addii(p3,p2); affii(p4,p1); cgiv(p4);  
             y[1]=laddsi(-1,(GEN)y[1]);  
           }  
           cgiv(p3); p4=p1; p1=p2; p2=p4;  
           if (f)  
           {  
             if (i>=lx1) { i--; break; }  
             if (!egalii((GEN)y[i],(GEN)x1[i]))  
             {  
               av=avma;  
               if (gcmp1(absi(subii((GEN)x1[i],(GEN)y[i]))))  
               {  
                 if (i>=lx1-1 || !gcmp1((GEN)x1[i+1]))  
                   affii((GEN)x1[i],(GEN)y[i]);  
               }  
               else i--;  
               avma=av; break;  
             }  
           }  
         }  
         if (i>=2 && gcmp1((GEN)y[i]))  
         {  
           cgiv((GEN)y[i]); --i; cgiv((GEN)y[i]);  
           y[i]=laddsi(1,(GEN)y[i]);  
         }  
         setlg(y,i+1); return gerepileupto(l, y);  
     }      }
     err(typeer,"sfcont");      err(typeer,"sfcont");
   }    }
   
   switch(tx)    switch(tx)
   {    {
     case t_POL:      case t_POL: return _veccopy(x);
       y=cgetg(2,t_VEC); y[1]=lcopy(x); break;  
     case t_SER:      case t_SER:
       av=avma; p1=gtrunc(x); tetpil=avma;        av = avma; p1 = gtrunc(x);
       y=gerepile(av,tetpil,sfcont(p1,p1,k)); break;        return gerepileupto(av,sfcont(p1,k));
     case t_RFRAC:      case t_RFRAC:
     case t_RFRACN:      case t_RFRACN:
       av=avma; l1=lgef(x[1]); if (lgef(x[2])>l1) l1=lgef(x[2]);        av = avma;
       if (k>0 && l1>k+1) l1=k+1;        l = typ(x[1]) == t_POL? lgef(x[1]): 3;
       yp=cgetg(l1,t_VEC); p1=(GEN)x[1]; i=0; p2=(GEN)x[2];        if (lgef(x[2]) > l) l = lgef(x[2]);
       while (!gcmp0(p2) && i<=l1-2)        if (k > 0 && l > k+1) l = k+1;
         y = cgetg(l,t_VEC);
         p1 = (GEN)x[1];
         p2 = (GEN)x[2];
         for (i=1; i<l; i++)
       {        {
         i++; yp[i]=ldivres(p1,p2,&p3);          y[i] = ldivres(p1,p2,&p3);
         p1=p2; p2=p3;          if (gcmp0(p3)) { i++; break; }
           p1 = p2; p2 = p3;
       }        }
       tetpil=avma; y=cgetg(i+1,t_VEC);        setlg(y, i); return gerepilecopy(av, y);
       for (j=1; j<=i; j++) y[j]=lcopy((GEN)yp[j]);  
       y=gerepile(av,tetpil,y); break;  
     default: err(typeer,"sfcont");  
       return NULL; /* not reached */  
   }    }
   return y;    err(typeer,"sfcont");
     return NULL; /* not reached */
 }  }
   
 static GEN  static GEN
 sfcont2(GEN b, GEN x, long k)  sfcont2(GEN b, GEN x, long k)
 {  {
   ulong av = avma;    gpmem_t av = avma;
   long l1 = lg(b), tx = typ(x), i;    long l1 = lg(b), tx = typ(x), i;
   GEN y,p1;    GEN y,p1;
   
Line 1846  sfcont2(GEN b, GEN x, long k)
Line 2118  sfcont2(GEN b, GEN x, long k)
   return gerepilecopy(av,y);    return gerepilecopy(av,y);
 }  }
   
   
 GEN  GEN
   gcf(GEN x)
   {
     return sfcont(x,0);
   }
   
   GEN
   gcf2(GEN b, GEN x)
   {
     return contfrac0(x,b,0);
   }
   
   GEN
   gboundcf(GEN x, long k)
   {
     return sfcont(x,k);
   }
   
   GEN
   contfrac0(GEN x, GEN b, long flag)
   {
     long lb,tb,i;
     GEN y;
   
     if (!b || gcmp0(b)) return sfcont(x,flag);
     tb = typ(b);
     if (tb == t_INT) return sfcont(x,itos(b));
     if (! is_matvec_t(tb)) err(typeer,"contfrac0");
   
     lb=lg(b); if (lb==1) return cgetg(1,t_VEC);
     if (tb != t_MAT) return sfcont2(b,x,flag);
     if (lg(b[1])==1) return sfcont(x,flag);
     y = (GEN) gpmalloc(lb*sizeof(long));
     for (i=1; i<lb; i++) y[i]=coeff(b,1,i);
     x = sfcont2(y,x,flag); free(y); return x;
   }
   
   GEN
 pnqn(GEN x)  pnqn(GEN x)
 {  {
   long av=avma,tetpil,lx,ly,tx=typ(x),i;    gpmem_t av=avma,tetpil;
     long lx,ly,tx=typ(x),i;
   GEN y,p0,p1,q0,q1,a,b,p2,q2;    GEN y,p0,p1,q0,q1,a,b,p2,q2;
   
   if (! is_matvec_t(tx)) err(typeer,"pnqn");    if (! is_matvec_t(tx)) err(typeer,"pnqn");
Line 1899  bestappr_mod(GEN x, GEN A, GEN B)
Line 2210  bestappr_mod(GEN x, GEN A, GEN B)
   {    {
     case t_INTMOD:      case t_INTMOD:
     {      {
       ulong av = avma;        gpmem_t av = avma;
       GEN a,b,d, t = cgetg(3, t_FRAC);        GEN a,b,d, t = cgetg(3, t_FRAC);
       if (! ratlift((GEN)x[2], (GEN)x[1], &a,&b,A,B)) return NULL;        if (! ratlift((GEN)x[2], (GEN)x[1], &a,&b,A,B)) return NULL;
       if (is_pm1(b)) return icopy_av(a, (GEN)av);        if (is_pm1(b)) return icopy_av(a, (GEN)av);
Line 1928  bestappr_mod(GEN x, GEN A, GEN B)
Line 2239  bestappr_mod(GEN x, GEN A, GEN B)
 GEN  GEN
 bestappr(GEN x, GEN k)  bestappr(GEN x, GEN k)
 {  {
   ulong av = avma, tetpil;    gpmem_t av = avma, tetpil;
   long tx,tk=typ(k),lx,i,e;    long tx,tk=typ(k),lx,i,e;
   GEN p0,p1,p,q0,q1,q,a,y;    GEN p0,p1,p,q0,q1,q,a,y;
   
Line 1939  bestappr(GEN x, GEN k)
Line 2250  bestappr(GEN x, GEN k)
     k = gcvtoi(k,&e);      k = gcvtoi(k,&e);
   }    }
   if (signe(k) <= 0) k=gun;    if (signe(k) <= 0) k=gun;
   tx=typ(x); if (tx == t_FRACN) x = gred(x);    tx=typ(x); if (tx == t_FRACN) { x = gred(x); tx = typ(x); }
   switch(tx)    switch(tx)
   {    {
     case t_INT:      case t_INT:
Line 1954  bestappr(GEN x, GEN k)
Line 2265  bestappr(GEN x, GEN k)
       }        }
   
     case t_REAL:      case t_REAL:
       p1=gun; p0=gfloor(x); q1=gzero; q0=gun; a=p0;        p1=gun; a=p0=gfloor(x); q1=gzero; q0=gun;
       while (gcmp(q0,k)<=0)        while (cmpii(q0,k)<=0)
       {        {
         x = gsub(x,a);          x = gsub(x,a);
         if (gcmp0(x)) { p1=p0; q1=q0; break; }          if (gcmp0(x)) { p1=p0; q1=q0; break; }
Line 1980  bestappr(GEN x, GEN k)
Line 2291  bestappr(GEN x, GEN k)
 GEN  GEN
 bestappr0(GEN x, GEN a, GEN b)  bestappr0(GEN x, GEN a, GEN b)
 {  {
   ulong av;    gpmem_t av;
   GEN t;    GEN t;
   if (!b) return bestappr(x,a);    if (!b) return bestappr(x,a);
   av = avma;    av = avma;
Line 2027  update_f(GEN f, GEN a)
Line 2338  update_f(GEN f, GEN a)
 GEN  GEN
 fundunit(GEN x)  fundunit(GEN x)
 {  {
   long av = avma, av2,tetpil,lim,r,flp,flq;    gpmem_t av = avma, av2, lim;
     long r,flp,flq;
   GEN pol,y,a,u,v,u1,v1,sqd,f;    GEN pol,y,a,u,v,u1,v1,sqd,f;
   
   if (typ(x) != t_INT) err(arither1);    if (typ(x) != t_INT) err(arither1);
Line 2056  fundunit(GEN x)
Line 2368  fundunit(GEN x)
   pol = quadpoly(x); y = get_quad(f,pol,r);    pol = quadpoly(x); y = get_quad(f,pol,r);
   if (!flq) v1 = y; else { update_f(f,a); v1 = get_quad(f,pol,r); }    if (!flq) v1 = y; else { update_f(f,a); v1 = get_quad(f,pol,r); }
   
   u1=gconj(y); tetpil=avma; y=gdiv(v1,u1);    u1=gconj(y); y=gdiv(v1,u1);
   if (signe(y[3]) < 0) { tetpil=avma; y=gneg(y); }    if (signe(y[3]) < 0) y = gneg(y);
   return gerepile(av,tetpil,y);    return gerepileupto(av, y);
 }  }
   
 GEN  GEN
Line 2070  gregula(GEN x, long prec)
Line 2382  gregula(GEN x, long prec)
 GEN  GEN
 regula(GEN x, long prec)  regula(GEN x, long prec)
 {  {
   long av = avma,av2,lim,r,fl,rexp;    gpmem_t av = avma,av2,lim;
     long r,fl,rexp;
   GEN reg,rsqd,y,u,v,u1,v1, sqd = racine(x);    GEN reg,rsqd,y,u,v,u1,v1, sqd = racine(x);
   
   if (signe(x)<=0) err(arither2);    if (signe(x)<=0) err(arither2);
Line 2079  regula(GEN x, long prec)
Line 2392  regula(GEN x, long prec)
   rsqd = gsqrt(x,prec);    rsqd = gsqrt(x,prec);
   if (egalii(sqri(sqd),x)) err(talker,"square argument in regula");    if (egalii(sqri(sqd),x)) err(talker,"square argument in regula");
   
   rexp=0; reg=cgetr(prec); affsr(2,reg);    rexp=0; reg = stor(2, prec);
   av2=avma; lim = stack_lim(av2,2);    av2=avma; lim = stack_lim(av2,2);
   u = stoi(r); v = gdeux;    u = stoi(r); v = gdeux;
   for(;;)    for(;;)
Line 2103  regula(GEN x, long prec)
Line 2416  regula(GEN x, long prec)
   y = mplog(divri(reg,v));    y = mplog(divri(reg,v));
   if (rexp)    if (rexp)
   {    {
     u1 = mulsr(rexp, glog(gdeux, prec));      u1 = mulsr(rexp, mplog2(prec));
     setexpo(u1, expo(u1)+1);      setexpo(u1, expo(u1)+1);
     y = addrr(y,u1);      y = addrr(y,u1);
   }    }
Line 2174  end_classno(GEN h, GEN hin, GEN forms, long lform)
Line 2487  end_classno(GEN h, GEN hin, GEN forms, long lform)
   q = ground(gdiv(hin,h)); /* approximate order of G/H */    q = ground(gdiv(hin,h)); /* approximate order of G/H */
   for (i=1; i < lform; i++)    for (i=1; i < lform; i++)
   {    {
     long av = avma;      gpmem_t av = avma;
     fg = powgi((GEN)forms[i], h);      fg = powgi((GEN)forms[i], h);
     fh = powgi(fg, q);      fh = powgi(fg, q);
     a = (GEN)fh[1];      a = (GEN)fh[1];
Line 2198  to_form(GEN x, long f)
Line 2511  to_form(GEN x, long f)
   return redimag(primeform(x, stoi(f), 0));    return redimag(primeform(x, stoi(f), 0));
 }  }
   
   static GEN
   conductor_part(GEN x, GEN *ptD, GEN *ptreg, GEN *ptfa)
   {
     long n,i,k,s=signe(x),fl2;
     GEN e,p,H,d,D,fa,reg;
   
     fa = auxdecomp(absi(x),1);
     e = gtovecsmall((GEN)fa[2]);
     fa = (GEN)fa[1];
     n = lg(fa); d = gun;
     for (i=1; i<n; i++)
       if (e[i] & 1) d = mulii(d,(GEN)fa[i]);
     if (mod4(d) == 2-s) fl2 = 0;
     else
     {
       fl2 = (mod4(x)==0);
       if (!fl2) err(funder2,"classno");
       d = shifti(d,2);
     }
     H = gun; D = (s<0)? negi(d): d; /* d = abs(D) */
     /* f \prod_{p|f}  [ 1 - (D/p) p^-1 ] */
     for (i=1; i<n; i++)
     {
       p = (GEN)fa[i];
       k = e[i];
       if (fl2 && i==1) k -= 2; /* p = 2 */
       if (k >= 2)
       {
         H = mulii(H, subis(p, kronecker(D,p)));
         if (k>=4) H = mulii(H,gpowgs(p,(k>>1)-1));
       }
     }
   
     /* divide by [ O^* : O_K^* ] */
     if (s < 0)
     {
       reg = NULL;
       if (!is_bigint(d))
         switch(itos(d))
         {
           case 4: H = divis(H,2); break;
           case 3: H = divis(H,3); break;
         }
     } else {
       reg = regula(D,DEFAULTPREC);
       if (!egalii(x,D))
         H = divii(H, ground(gdiv(regula(x,DEFAULTPREC), reg)));
     }
     *ptfa = fa; *ptD = D; if (ptreg) *ptreg = reg;
     return H;
   }
   
   
 static long  static long
 two_rank(GEN x)  two_rank(GEN x)
 {  {
Line 2215  two_rank(GEN x)
Line 2581  two_rank(GEN x)
 }  }
   
 #define MAXFORM 11  #define MAXFORM 11
 #define _low(x) (__x=(GEN)x, modBIL(__x))  #if 0 /* gcc-ism */
   #  define _low(x) ({GEN __x=(GEN)x; modBIL(__x);})
   #else
   #  define _low(x) (__x = (GEN)(x), modBIL(__x))
   #endif
   
 /* h(x) for x<0 using Baby Step/Giant Step.  /* h(x) for x<0 using Baby Step/Giant Step.
  * Assumes G is not too far from being cyclic.   * Assumes G is not too far from being cyclic.
Line 2224  two_rank(GEN x)
Line 2594  two_rank(GEN x)
 GEN  GEN
 classno(GEN x)  classno(GEN x)
 {  {
   long av = avma, av2,c,lforms,k,l,i,j,j1,j2,lim,com,s, forms[MAXFORM];    gpmem_t av = avma, av2, lim;
   long r2;    long r2,c,lforms,k,l,i,j,j1,j2,com,s, forms[MAXFORM];
   GEN a,b,count,index,tabla,tablb,hash,p1,p2,hin,h,f,fh,fg,ftest;    GEN a,b,count,index,tabla,tablb,hash,p1,p2,hin,h,f,fh,fg,ftest;
   GEN __x;    GEN Hf,D,fa;
   byteptr p = diffptr;    byteptr p = diffptr;
     GEN __x;
   
   if (typ(x) != t_INT) err(arither1);    if (typ(x) != t_INT) err(arither1);
   s=signe(x); if (s>=0) return classno2(x);    s=signe(x); if (s>=0) return classno2(x);
   
   k=mod4(x); if (k==1 || k==2) err(funder2,"classno");    k=mod4(x); if (k==1 || k==2) err(funder2,"classno");
   if (gcmpgs(x,-12) >= 0) return gun;    if (cmpis(x,-12) >= 0) return gun;
   
   p2 = gsqrt(absi(x),DEFAULTPREC);    Hf = conductor_part(x,&D,NULL,&fa);
     if (cmpis(D,-12) >= 0) return gerepilecopy(av, Hf);
   
     p2 = gsqrt(absi(D),DEFAULTPREC);
   p1 = divrr(p2,mppi(DEFAULTPREC));    p1 = divrr(p2,mppi(DEFAULTPREC));
   p2 = gtrunc(shiftr(gsqrt(p2,DEFAULTPREC),1));    p2 = gtrunc(shiftr(gsqrt(p2,DEFAULTPREC),1));
   s = 1000;    s = 1000;
Line 2245  classno(GEN x)
Line 2619  classno(GEN x)
     if (is_bigint(p2)) err(talker,"discriminant too big in classno");      if (is_bigint(p2)) err(talker,"discriminant too big in classno");
     s = itos(p2); if (s < 1000) s = 1000;      s = itos(p2); if (s < 1000) s = 1000;
   }    }
   r2 = two_rank(x);    r2 = two_rank(D);
   c=lforms=0;    c=lforms=0;
   while (c <= s && *p)    while (c <= s && *p)
   {    {
     c += *p++; k = krogs(x,c);      NEXT_PRIME_VIADIFF(c,p);
       k = krogs(D,c);
     if (!k) continue;      if (!k) continue;
   
     av2 = avma;      av2 = avma;
Line 2270  classno(GEN x)
Line 2645  classno(GEN x)
   tabla = new_chunk(10000);    tabla = new_chunk(10000);
   tablb = new_chunk(10000);    tablb = new_chunk(10000);
   hash  = new_chunk(10000);    hash  = new_chunk(10000);
   f = gsqr(to_form(x, forms[0]));    f = gsqr(to_form(D, forms[0]));
   p1 = fh = powgi(f, h);    p1 = fh = powgi(f, h);
   for (i=0; i<s; i++)    for (i=0; i<s; i++)
   {    {
Line 2300  classno(GEN x)
Line 2675  classno(GEN x)
           h = addii(addis(h,j2), mulss(s,com));            h = addii(addis(h,j2), mulss(s,com));
           forms[0] = (long)f;            forms[0] = (long)f;
           for (i=1; i<lforms; i++)            for (i=1; i<lforms; i++)
             forms[i] = (long)gsqr(to_form(x, forms[i]));              forms[i] = (long)gsqr(to_form(D, forms[i]));
           h = end_classno(h,hin,forms,lforms);            h = end_classno(h,hin,forms,lforms);
             h = mulii(h,Hf);
           return gerepileuptoint(av, shifti(h, r2));            return gerepileuptoint(av, shifti(h, r2));
         }          }
       }        }
Line 2316  classno(GEN x)
Line 2692  classno(GEN x)
 GEN  GEN
 classno2(GEN x)  classno2(GEN x)
 {  {
   long av=avma,tetpil,n,i,k,s=signe(x),fl2;    gpmem_t av = avma;
   GEN p1,p2,p3,p4,p5,p7,p8,pi4,reg,logd,d,fd;    long n,i,k,s = signe(x);
     GEN p1,p2,p3,p4,p5,p7,Hf,Pi,reg,logd,d,D;
   
   if (typ(x) != t_INT) err(arither1);    if (typ(x) != t_INT) err(arither1);
   if (!s) err(arither2);    if (!s) err(arither2);
   if (s<0 && gcmpgs(x,-12) >= 0) return gun;    if (s < 0 && cmpis(x,-12) >= 0) return gun;
   
   p1=auxdecomp(absi(x),1); p2=(GEN)p1[2]; p1=(GEN)p1[1];    Hf = conductor_part(x, &D, &reg, &p1);
   n=lg(p1); d=gun;    if (s < 0 && cmpis(D,-12) >= 0)
   for (i=1; i<n; i++)      return gerepilecopy(av, Hf); /* |D| < 12*/
     if (mod2((GEN)p2[i])) d=mulii(d,(GEN)p1[i]);  
   fl2 = (mod4(x)==0); /* 4 | x */  
   if (mod4(d) != 2-s)  
   {  
     if (!fl2) err(funder2,"classno2");  
     d = shifti(d,2);  
   }  
   else fl2=0;  
   p8=gun; fd = (s<0)? negi(d): d; /* d = abs(fd) */  
   for (i=1; i<n; i++)  
   {  
     k=itos((GEN)p2[i]); p4=(GEN)p1[i];  
     if (fl2 && i==1) k -= 2; /* p4 = 2 */  
     if (k>=2)  
     {  
       p8=mulii(p8,subis(p4,kronecker(fd,p4)));  
       if (k>=4) p8=mulii(p8,gpuigs(p4,(k>>1)-1));  
     }  
   }  
   if (s<0 && lgefint(d)==3)  
   {  
     switch(d[2])  
     {  
       case 4: p8=gdivgs(p8,2); break;  
       case 3: p8=gdivgs(p8,3); break;  
     }  
     if (d[2] < 12) /* |fd| < 12*/  
       { tetpil=avma; return gerepile(av,tetpil,icopy(p8)); }  
   }  
   
   pi4 = mppi(DEFAULTPREC); logd = glog(d,DEFAULTPREC);    Pi = mppi(DEFAULTPREC);
   p1 = gsqrt(gdiv(gmul(d,logd),gmul2n(pi4,1)),DEFAULTPREC);    d = absi(D); logd = glog(d,DEFAULTPREC);
   p4 = divri(pi4,d); p7 = ginv(mpsqrt(pi4));    p1 = mpsqrt(gdiv(gmul(d,logd), gmul2n(Pi,1)));
   if (s > 0)    if (s > 0)
   {    {
     reg = regula(d,DEFAULTPREC);      p2 = subsr(1, gmul2n(divrr(mplog(reg),logd),1));
     p2 = gsubsg(1, gmul2n(gdiv(glog(reg,DEFAULTPREC),logd),1));      if (gcmp(gsqr(p2), divsr(2,logd)) >= 0) p1 = gmul(p2,p1);
     p3 = gsqrt(gdivsg(2,logd),DEFAULTPREC);  
     if (gcmp(p2,p3)>=0) p1 = gmul(p2,p1);  
   }    }
   else reg = NULL; /* for gcc -Wall */    p1 = gtrunc(p1);
   p1 = gtrunc(p1); n=p1[2];    if (is_bigint(p1)) err(talker,"discriminant too large in classno");
   if (lgefint(p1)!=3 || n<0)    n = itos(p1);
     err(talker,"discriminant too large in classno");  
   
     p4 = divri(Pi,d); p7 = ginv(mpsqrt(Pi));
   p1 = gsqrt(d,DEFAULTPREC); p3 = gzero;    p1 = gsqrt(d,DEFAULTPREC); p3 = gzero;
   if (s > 0)    if (s > 0)
   {    {
     for (i=1; i<=n; i++)      for (i=1; i<=n; i++)
     {      {
       k=krogs(fd,i);        k = krogs(D,i); if (!k) continue;
       if (k)        p2 = mulir(mulss(i,i),p4);
       {        p5 = subsr(1,mulrr(p7,incgam3(ghalf,p2,DEFAULTPREC)));
         p2 = mulir(mulss(i,i),p4);        p5 = addrr(divrs(mulrr(p1,p5),i), eint1(p2,DEFAULTPREC));
         p5 = subsr(1,mulrr(p7,incgam3(ghalf,p2,DEFAULTPREC)));        p3 = (k>0)? addrr(p3,p5): subrr(p3,p5);
         p5 = addrr(divrs(mulrr(p1,p5),i),eint1(p2,DEFAULTPREC));  
         p3 = (k>0)? addrr(p3,p5): subrr(p3,p5);  
       }  
     }      }
     p3 = shiftr(divrr(p3,reg),-1);      p3 = shiftr(divrr(p3,reg),-1);
     if (!egalii(x,fd)) /* x != p3 */  
       p8 = gdiv(p8,ground(gdiv(regula(x,DEFAULTPREC),reg)));  
   }    }
   else    else
   {    {
     p1 = gdiv(p1,pi4);      p1 = gdiv(p1,Pi);
     for (i=1; i<=n; i++)      for (i=1; i<=n; i++)
     {      {
       k=krogs(fd,i);        k = krogs(D,i); if (!k) continue;
       if (k)        p2 = mulir(mulss(i,i),p4);
       {        p5 = subsr(1,mulrr(p7,incgam3(ghalf,p2,DEFAULTPREC)));
         p2=mulir(mulss(i,i),p4);        p5 = addrr(p5, divrr(divrs(p1,i),mpexp(p2)));
         p5=subsr(1,mulrr(p7,incgam3(ghalf,p2,DEFAULTPREC)));        p3 = (k>0)? addrr(p3,p5): subrr(p3,p5);
         p5=addrr(p5,divrr(divrs(p1,i),mpexp(p2)));  
         p3 = (k>0)? addrr(p3,p5): subrr(p3,p5);  
       }  
     }      }
   }    }
   p3=ground(p3); tetpil=avma;    return gerepileuptoint(av, mulii(Hf,ground(p3)));
   return gerepile(av,tetpil,mulii(p8,p3));  
 }  }
   
 GEN  GEN
Legend:
Removed from v.1.1.1.1  
changed lines
  Added in v.1.2
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>