===================================================================
RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/arith1.c,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -p -r1.1 -r1.2
--- OpenXM_contrib/pari-2.2/src/basemath/Attic/arith1.c	2001/10/02 11:17:01	1.1
+++ OpenXM_contrib/pari-2.2/src/basemath/Attic/arith1.c	2002/09/11 07:26:48	1.2
@@ -1,4 +1,4 @@
-/* $Id: arith1.c,v 1.1 2001/10/02 11:17:01 noro Exp $
+/* $Id: arith1.c,v 1.2 2002/09/11 07:26:48 noro Exp $
 
 Copyright (C) 2000  The PARI group.
 
@@ -115,12 +115,36 @@ garith_proto2gs(GEN f(GEN,long), GEN x, long y)
 }
 
 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)
 {
   int tx=typ(x);
   if (!y)
   {
-    ulong av=avma;
+    gpmem_t av = avma;
     if (tx!=t_VEC && tx!=t_COL)
       err(typeer,"association");
     return gerepileupto(av,divide_conquer_prod(x,f));
@@ -137,7 +161,8 @@ gassoc_proto(GEN f(GEN,GEN), GEN x, GEN y)
 GEN
 order(GEN x)
 {
-  long av=avma,av1,i,e;
+  gpmem_t av = avma,av1;
+  long i,e;
   GEN o,m,p;
 
   if (typ(x) != t_INTMOD || !gcmp1(mppgcd((GEN) x[1],(GEN) x[2])))
@@ -172,14 +197,15 @@ ggener(GEN m)
 GEN
 gener(GEN m)
 {
-  long av=avma,av1,k,i,e;
+  gpmem_t av=avma,av1;
+  long k,i,e;
   GEN x,t,q,p;
 
   if (typ(m) != t_INT) err(arither1);
   e = signe(m);
   if (!e) err(talker,"zero modulus in znprimroot");
-  if (is_pm1(m)) { avma=av; return gmodulss(0,1); }
-  if (e<0) m = absi(m);
+  if (is_pm1(m)) return gmodulss(0,1);
+  if (e < 0) m = absi(m);
 
   e = mod4(m);
   if (e == 0) /* m = 0 mod 4 */
@@ -196,33 +222,54 @@ gener(GEN m)
 
   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);
-  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);
     av1=avma; return gerepile(av,av1,gmodulcp(x,m));
   }
 
   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(;;)
   {
     x[2]++;
     if (gcmp1(mppgcd(m,x)))
     {
       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;
     }
   }
   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
 znstar(GEN n)
 {
   GEN p1,z,q,u,v,d,list,ep,h,gen,moduli,p,a;
   long i,j,nbp,sizeh;
-  ulong av;
+  gpmem_t av;
 
   if (typ(n) != t_INT) err(arither1);
   if (!signe(n))
@@ -315,7 +362,8 @@ gracine(GEN a)
 static GEN
 racine_r(GEN a, long l)
 {
-  long av,s;
+  gpmem_t av;
+  long s;
   GEN x,y,z;
 
   if (l <= 4) return utoi(mpsqrtl(a));
@@ -331,20 +379,21 @@ racine_r(GEN a, long l)
     y = x; x = z;
   }
   while (cmpii(x,y) < 0);
-  avma = (long)y;
+  avma = (gpmem_t)y;
   return gerepileuptoint(av,y);
 }
 
 GEN
 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);
   if (roundup && signe(x))
   {
     m = modBIL(x);
     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));
   }
   return x;
@@ -409,7 +458,7 @@ ucarrecomplet(ulong A)
 long
 carrecomplet(GEN x, GEN *pt)
 {
-  long av;
+  gpmem_t av;
   GEN y;
 
   switch(signe(x))
@@ -427,14 +476,15 @@ carrecomplet(GEN x, GEN *pt)
   if (!carremod((ulong)smodis(x, 64*63*65*11))) return 0;
   av=avma; y = racine(x);
   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;
 }
 
 static GEN
 polcarrecomplet(GEN x, GEN *pt)
 {
-  long i,l,av,av2;
+  gpmem_t av,av2;
+  long i,l;
   GEN y,a,b;
 
   if (!signe(x)) return gun;
@@ -493,8 +543,9 @@ gcarrecomplet(GEN x, GEN *pt)
 GEN
 gcarreparfait(GEN x)
 {
+  gpmem_t av;
   GEN p1,a,p;
-  long tx=typ(x),l,i,av,v;
+  long tx=typ(x),l,i,v;
 
   switch(tx)
   {
@@ -506,7 +557,7 @@ gcarreparfait(GEN x)
 
     case t_INTMOD:
     {
-      GEN b, q, e;
+      GEN b, q;
       long w;
       a = (GEN)x[2]; if (!signe(a)) return gun;
       av = avma;
@@ -533,23 +584,18 @@ gcarreparfait(GEN x)
       if (i==0)
       {
         GEN d = mppgcd(a,q);
-        p1 = factor(d);
-        p = (GEN)p1[1];
-        e = (GEN)p1[2]; l = lg(e);
+        p = (GEN)factor(d)[1]; l = lg(p);
         for (i=1; i<l; i++)
         {
           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)) 
             { avma = av; return gzero; }
-          q = diviiexact(q, gpowgs((GEN)p[i], w));
         }
         if (kronecker(a,q) == -1) { avma = av; return gzero; }
       }
       /* kro(a,q) = 1, q odd: need to factor q */
-      p1 = factor(q);
-      p = (GEN)p1[1]; 
-      e = (GEN)p1[2]; l = lg(e) - 1;
+      p = (GEN)factor(q)[1]; l = lg(p) - 1;
       /* kro(a,q) = 1, check all p|q but the last (product formula) */
       for (i=1; i<l; i++)
         if (kronecker(a,(GEN)p[i]) == -1) { avma = av; return gzero; }
@@ -618,10 +664,10 @@ gkronecker(GEN x, GEN y)
 long
 kronecker(GEN x, GEN y)
 {
+  const gpmem_t av = avma;
   GEN z;
-  long av,r,s=1;
+  long r,s=1;
 
-  av=avma;
   switch (signe(y))
   {
     case -1: y=negi(y); if (signe(x)<0) s= -1; break;
@@ -662,9 +708,9 @@ gkrogs(GEN x, long y)
 long
 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) { y = -y; if (signe(x)<0) s = -1; }
@@ -698,7 +744,8 @@ krogs(GEN x, long y)
 long
 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;
 }
 
@@ -747,7 +794,8 @@ hil0(GEN x, GEN y, GEN p)
 long
 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;
 
   if (gcmp0(x) || gcmp0(y)) return 0;
@@ -779,15 +827,17 @@ hil(GEN x, GEN y, GEN p)
 	case t_REAL:
 	  return (signe(x)<0 && signe(y)<0)? -1: 1;
 	case t_INTMOD:
-	  if (egalii(gdeux,(GEN)y[1])) err(hiler1);
-	  return hil(x,(GEN)y[2],(GEN)y[1]);
+          p = (GEN)y[1];
+	  if (egalii(gdeux,p)) err(hiler1);
+	  return hil(x,(GEN)y[2],p);
 	case t_FRAC: case t_FRACN:
 	  p1=mulii((GEN)y[1],(GEN)y[2]); z=hil(x,p1,p);
 	  avma=av; return z;
 	case t_PADIC:
-	  if (egalii(gdeux,(GEN)y[2]) && precp(y) <= 2) err(hiler1);
-	  p1 = odd(valp(y))? mulii((GEN)y[2],(GEN)y[4]): (GEN)y[4];
-	  z=hil(x,p1,(GEN)y[2]); avma=av; return z;
+          p = (GEN)y[2];
+	  if (egalii(gdeux,p) && precp(y) <= 1) err(hiler1);
+	  p1 = odd(valp(y))? mulii(p,(GEN)y[4]): (GEN)y[4];
+	  z=hil(x,p1,p); avma=av; return z;
       }
       break;
 
@@ -797,17 +847,18 @@ hil(GEN x, GEN y, GEN p)
       return signe(y[1])*signe(y[2]);
 
     case t_INTMOD:
-      if (egalii(gdeux,(GEN)y[1])) err(hiler1);
+      p = (GEN)x[1];
+      if (egalii(gdeux,p)) err(hiler1);
       switch(ty)
       {
         case t_INTMOD:
-          if (!egalii((GEN)x[1],(GEN)y[1])) break;
-          return hil((GEN)x[2],(GEN)y[2],(GEN)x[1]);
+          if (!egalii(p, (GEN)y[1])) break;
+          return hil((GEN)x[2],(GEN)y[2],p);
         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:
-          if (!egalii((GEN)x[1],(GEN)y[2])) break;
-          return hil((GEN)x[2],y,(GEN)x[1]);
+          if (!egalii(p, (GEN)y[2])) break;
+          return hil((GEN)x[2],y,p);
       }
       break;
 
@@ -824,10 +875,12 @@ hil(GEN x, GEN y, GEN p)
       break;
 
     case t_PADIC:
-      if (ty != t_PADIC || !egalii((GEN)x[2],(GEN)y[2])) break;
-      p1 = odd(valp(x))? mulii((GEN)x[2],(GEN)x[4]): (GEN)x[4];
-      p2 = odd(valp(y))? mulii((GEN)y[2],(GEN)y[4]): (GEN)y[4];
-      z=hil(p1,p2,(GEN)x[2]); avma=av; return z;
+      p = (GEN)x[2];
+      if (ty != t_PADIC || !egalii(p,(GEN)y[2])) break;
+      if (egalii(p, gdeux) && (precp(x) <= 1 || precp(y) <= 1)) err(hiler1);
+      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");
   return 0; /* not reached */
@@ -844,16 +897,28 @@ hil(GEN x, GEN y, GEN 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
 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;
 
   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");
   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 */
   {
     avma = av;
@@ -911,6 +976,113 @@ mpsqrtmod(GEN a, GEN p)
   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                        */
@@ -930,13 +1102,13 @@ mpsqrtmod(GEN a, GEN p)
 static GEN
 mplgenmod(GEN l, long e, GEN r,GEN p,GEN *zeta)
 {
-  ulong av1;
+  const gpmem_t av1 = avma;
   GEN m,m1;
   long k,i; 
-  av1 = avma;
   for (k=1; ; k++)
   {
     m1 = m = powmodulo(stoi(k+1),r,p);
+    if (gcmp1(m)) {avma=av1; continue;}
     for (i=1; i<e; i++)
       if (gcmp1(m=powmodulo(m,l,p))) break;
     if (i==e) break;
@@ -960,11 +1132,11 @@ GEN Fp_shanks(GEN x,GEN g0,GEN p, GEN q);
 static GEN
 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;
   GEN p1,u1,u2,v,w,z,dl;
   /* 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);
   w=powmodulo(a,modii(mulii(negi(u1),r),q),p);
   lim = stack_lim(av,1);
@@ -1015,7 +1187,7 @@ If a=0 ,return 0 and if zetan is not NULL zetan is set
 GEN 
 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 q,z;
   GEN *gptr[2];
@@ -1183,7 +1355,8 @@ ulong mppgcd_resiu(GEN y, ulong x);
 GEN
 mppgcd(GEN a, GEN b)
 {
-  long av,v,w;
+  long v, w;
+  gpmem_t av;
   GEN t, p1;
 
   if (typ(a) != t_INT || typ(b) != t_INT) err(arither1);
@@ -1204,7 +1377,7 @@ mppgcd(GEN a, GEN b)
   }
 
   /* 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);
   if (!signe(t)) { avma=av; return absi(b); }
 
@@ -1237,7 +1410,7 @@ 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]);
     avma = av; return shifti(r,v);
   }
@@ -1246,7 +1419,7 @@ mppgcd(GEN a, GEN b)
 GEN
 mpppcm(GEN x, GEN y)
 {
-  ulong av;
+  gpmem_t av;
   GEN p1,p2;
   if (typ(x) != t_INT || typ(y) != t_INT) err(arither1);
   if (!signe(x)) return gzero;
@@ -1292,8 +1465,8 @@ GEN chinese(GEN x, GEN y)
 GEN
 chinois(GEN x, GEN y)
 {
-  ulong av,tetpil, tx = typ(x);
-  long i,lx,vx;
+  gpmem_t av,tetpil;
+  long i,lx,vx, tx = typ(x);
   GEN z,p1,p2,d,u,v;
 
   if (gegal(x,y)) return gcopy(x);
@@ -1339,7 +1512,7 @@ chinois(GEN x, GEN y)
 GEN
 chinois_int_coprime(GEN x2, GEN y2, GEN x1, GEN y1, GEN z1)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   GEN ax,p1;
   (void)new_chunk((lgefint(z1)<<1)+lgefint(x1)+lgefint(y1)); /* HACK */
   ax = mulii(mpinvmod(x1,y1), x1);
@@ -1364,91 +1537,166 @@ mpinvmod(GEN a, GEN m)
 
 /*********************************************************************/
 /**                                                                 **/
-/**                   POWERING MODULO (A^N mod M)                   **/
+/**                    MODULAR EXPONENTIATION                       **/
 /**                                                                 **/
 /*********************************************************************/
-GEN resiimul(GEN x, GEN y);
-GEN resmod2n(GEN x, long y);
-GEN _resii(GEN x, GEN y) { return resii(x,y); }
+extern ulong invrev(ulong b);
+extern GEN resiimul(GEN x, GEN 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
 init_remainder(GEN M)
 {
   GEN sM = cgetg(3, t_VEC);
-  GEN Mr = cgetr(lgefint(M) + 1);
-  affir(M, Mr);
   sM[1] = (long)M;
-  sM[2] = (long)linv(Mr);
+  sM[2] = (long)linv( itor(M, lgefint(M) + 1) ); /* 1. / M */
   return sM;
 }
 
 /* optimal on UltraSPARC */
-static long RESIIMUL_LIMIT = 150;
+static long RESIIMUL_LIMIT   = 150;
+static long MONTGOMERY_LIMIT = 32;
 
 void
 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
-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;
-  long av = avma,av1,lim,man,k,nb,s, *p;
-  GEN (*mul)(GEN,GEN) = mulii;
-  GEN (*res)(GEN,GEN) = _resii;
+  muldata  D;
+  montdata S;
 
-  if (typ(A) != t_INT || typ(N) != t_INT || typ(M) != t_INT) err(arither1);
-  s = signe(N);
+  if (typ(A) != t_INT || typ(k) != t_INT || typ(N) != t_INT) err(arither1);
+  s = signe(k);
   if (!s)
   {
-    k = signe(resii(A,M)); avma=av;
-    return k? gun: gzero;
+    t = signe(resii(A,N)); avma = av;
+    return t? gun: gzero;
   }
-  if (s < 0) { A = mpinvmod(A,M); N = absi(N); }
+  if (s < 0) y = mpinvmod(A,N);
   else
   {
-    A = modii(A,M);
-    if (!signe(A)) { avma=av; return gzero; }
+    y = modii(A,N);
+    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 */
-      avma=av; return gun;
-    case 2: /* y = 2, use shifti not mulii */
-    mul = (GEN (*)(GEN,GEN))shifti; A = (GEN)1L;
+    case 1: avma = av; return gun;
+    case 2: base_is_2 = 1; break;
   }
 
   /* TODO: Move this out of here and use for general modular computations */
-  if ((k = vali(M)) && k == expi(M))
-  { /* M is a power of 2 */
-    M = (GEN)k;
-    res = (GEN(*)(GEN,GEN))resmod2n;
+  lN = lgefint(N);
+  use_montgomery = mod2(N) && lN < MONTGOMERY_LIMIT;
+  if (use_montgomery)
+  {
+    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))
-  { /* compute x % M using multiplication by 1./M */
-    M = init_remainder(M);
-    res = (GEN(*)(GEN,GEN))resiimul;
+  else if (lN > RESIIMUL_LIMIT
+       && (lgefint(k) > 3 || (((double)k[2])*expi(A)) > 2 + expi(N)))
+  {
+    D.mulred = base_is_2? &_muli2invred: &_muliiinvred;
+    D.res = &resiimul;
+    D.N = init_remainder(N);
   }
-
-  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;;)
+  else
   {
-    for (; k; man<<=1,k--)
-    {
-      y = res(sqri(y), M);
-      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;
+    D.mulred = base_is_2? &_muli2red: &_muliired;
+    D.res = &_resii;
+    D.N = N;
   }
-  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);
 }
 
 /*********************************************************************/
@@ -1457,84 +1705,115 @@ powmodulo(GEN A, GEN N, GEN M)
 /**                                                                 **/
 /*********************************************************************/
 GEN
-gnextprime(GEN n)
-{
-  return garith_proto(nextprime,n,0);
-}
+gnextprime(GEN n) { return garith_proto(nextprime,n,0); }
 
 GEN
-gprecprime(GEN n)
-{
-  return garith_proto(precprime,n,0);
-}
+gprecprime(GEN n) { return garith_proto(precprime,n,0); }
 
 GEN
 gisprime(GEN x, long flag)
 {
   switch (flag)
   {
-  case 0:
-    return arith_proto(isprime,x,1);
-  case 1:
-    return garith_proto2gs(plisprime,x,0);
-  case 2:
-    return garith_proto2gs(plisprime,x,1);
+    case 0: return arith_proto(isprime,x,1);
+    case 1: return garith_proto2gs(plisprime,x,1);
+    case 2: return arith_proto(isprimeAPRCL,x,1);
   }
   err(flagerr,"gisprime");
   return 0;
 }
 
 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
-gispsp(GEN x)
+/* assume x a BSW pseudoprime */
+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
-ispsp(GEN x)
+isprime(GEN x)
 {
-  return millerrabin(x,1);
+  return BSW_psp(x) && BSW_isprime(x);
 }
 
 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                    **/
 /**                                                                 **/
 /*********************************************************************/
-
-/* 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
-gisfundamental(GEN x)
-{
-  return arith_proto(isfundamental,x,1);
-}
+gisfundamental(GEN x) { return arith_proto(isfundamental,x,1); }
 
 long
 isfundamental(GEN x)
 {
-  long av,r;
+  long r;
   GEN p1;
 
   if (gcmp0(x)) return 0;
   r=mod4(x);
   if (!r)
   {
-    av=avma; p1=shifti(x,-2);
+    const gpmem_t av = avma;
+    p1=shifti(x,-2);
     r=mod4(p1); if (!r) return 0;
     if (signe(x)<0) r=4-r;
     r = (r==1)? 0: issquarefree(p1);
@@ -1547,17 +1826,18 @@ isfundamental(GEN x)
 GEN
 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;
 
   if (tx != t_INT && ! is_frac_t(tx)) err(arither1);
   f=factor(x); p1=(GEN)f[1]; p2=(GEN)f[2];
-  s=gun;
+  s = gun;
   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;
-  if (r>1) { tetpil=avma; s=shifti(s,2); }
-  return gerepile(av,tetpil,s);
+  if (r>1) s = shifti(s,2);
+  return gerepileuptoint(av, s);
 }
 
 /*********************************************************************/
@@ -1565,12 +1845,11 @@ quaddisc(GEN x)
 /**                              FACTORIAL                          **/
 /**                                                                 **/
 /*********************************************************************/
-GEN muluu(ulong x, ulong y);
-
 GEN
 mpfact(long n)
 {
-  long lx,k,l, av = avma;
+  const gpmem_t av = avma;
+  long lx,k,l;
   GEN x;
 
   if (n<2)
@@ -1598,10 +1877,10 @@ mpfact(long n)
 GEN
 mpfactr(long n, long prec)
 {
-  long av = avma, lim,k;
-  GEN f = cgetr(prec);
+  gpmem_t av = avma, lim;
+  long k;
+  GEN f = realun(prec);
 
-  affsr(1,f);
   if (n<2)
   {
     if (n<0) err(facter);
@@ -1629,36 +1908,37 @@ mpfactr(long n, long prec)
 void
 lucas(long n, GEN *ln, GEN *ln1)
 {
-  long taille,av;
+  gpmem_t av;
+  long taille;
   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);
-  *ln=cgeti(taille); *ln1=cgeti(taille);
-  av=avma; lucas(n/2,&z,&t);
+  taille = 3 + (long)(pariC3 * (1+labs(n)));
+  *ln = cgeti(taille);
+  *ln1= cgeti(taille);
+  av = avma; lucas(n/2, &z, &t);
   switch(n % 4)
   {
     case -3:
-      addsiz(2,sqri(z),*ln1);
-      subiiz(addsi(1,mulii(z,t)),*ln1,*ln); break;
-    case -2:
-      addsiz(2,sqri(z),*ln); addsiz(1,mulii(z,t),*ln1); break;
+      addsiz(2,sqri(z), *ln1);
+      subiiz(addsi(1,mulii(z,t)),*ln1, *ln); break;
     case -1:
-      subisz(sqri(z),2,*ln1);
-      subiiz(subis(mulii(z,t),1),*ln1,*ln); 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  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);
+      subisz(sqri(z),2, *ln1);
+      subiiz(subis(mulii(z,t),1),*ln1, *ln); 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 -2:
+    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
 fibo(long n)
 {
-  long av = avma;
+  gpmem_t av = avma;
   GEN ln,ln1;
 
   lucas(n-1,&ln,&ln1);
@@ -1670,149 +1950,141 @@ fibo(long n)
 /*                      CONTINUED FRACTIONS                        */
 /*                                                                 */
 /*******************************************************************/
-static GEN sfcont2(GEN b, GEN x, long k);
-
-GEN
-gcf(GEN x)
+static GEN
+icopy_lg(GEN x, long l)
 {
-  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
-gcf2(GEN b, GEN x)
+/* if y != NULL, stop as soon as partial quotients differ from y */
+static GEN
+Qsfcont(GEN x, GEN y, long k)
 {
-  return contfrac0(x,b,0);
-}
+  GEN  z, p1, p2, p3;
+  long i, l, ly;
 
-GEN
-gboundcf(GEN x, long k)
-{
-  return sfcont(x,x,k);
-}
+  ly = lgefint(x[2]);
+  l = 3 + (long) ((ly-2) / pariC3);
+  if (k > 0 && ++k > 0 && l > k) l = k; /* beware overflow */
+  if ((ulong)l > LGBITS) l = LGBITS;
 
-GEN
-contfrac0(GEN x, GEN b, long flag)
-{
-  long lb,tb,i;
-  GEN y;
-
-  if (!b || gcmp0(b)) return sfcont(x,x,flag);
-  tb = typ(b);
-  if (tb == t_INT) return sfcont(x,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,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;
+  p1 = (GEN)x[1];
+  p2 = (GEN)x[2];
+  z = cgetg(l,t_VEC);
+  l--;
+  if (y) {
+    gpmem_t av = avma;
+    if (l >= lg(y)) l = lg(y)-1;
+    for (i = 1; i <= l; i++)
+    {
+      z[i] = y[i];
+      p3 = p2; if (!gcmp1((GEN)z[i])) p3 = mulii((GEN)z[i], p2);
+      p3 = subii(p1, p3);
+      if (signe(p3) < 0)
+      { /* partial quotient too large */
+        p3 = addii(p3, p2);
+        if (signe(p3) >= 0) i++; /* by 1 */
+        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
-sfcont(GEN x, GEN x1, long k)
+static GEN
+sfcont(GEN x, long k)
 {
-  long  av,lx=lg(x),tx=typ(x),e,i,j,l,l1,lx1,tetpil,f;
-  GEN   y,p1,p2,p3,p4,yp;
+  gpmem_t av;
+  long lx,tx=typ(x),e,i,l;
+  GEN  y,p1,p2,p3;
 
   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)
     {
-      case t_INT:
-        y=cgetg(2,t_VEC); y[1]=lcopy(x); return y;
+      case t_INT: y = cgetg(2,t_VEC); y[1] = (long)icopy(x); return y;
       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);
-	p1[1] = (long) p2; e = bit_accuracy(lx)-1-expo(x);
-	if (e<0) err(talker,"integral part not significant in scfont");
+        e = bit_accuracy(lx)-1-expo(x);
+	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:
-        l=avma; lx1=lgefint(x[2]);
-	l1 = (long) ((double) BYTES_IN_LONG/4.0*46.093443*(lx1-2)+3);
-        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);
+        av = avma;
+        return gerepileupto(av, Qsfcont(x, NULL, k));
     }
     err(typeer,"sfcont");
   }
 
   switch(tx)
   {
-    case t_POL:
-      y=cgetg(2,t_VEC); y[1]=lcopy(x); break;
+    case t_POL: return _veccopy(x);
     case t_SER:
-      av=avma; p1=gtrunc(x); tetpil=avma;
-      y=gerepile(av,tetpil,sfcont(p1,p1,k)); break;
+      av = avma; p1 = gtrunc(x);
+      return gerepileupto(av,sfcont(p1,k));
     case t_RFRAC:
     case t_RFRACN:
-      av=avma; l1=lgef(x[1]); if (lgef(x[2])>l1) l1=lgef(x[2]);
-      if (k>0 && l1>k+1) l1=k+1;
-      yp=cgetg(l1,t_VEC); p1=(GEN)x[1]; i=0; p2=(GEN)x[2];
-      while (!gcmp0(p2) && i<=l1-2)
+      av = avma;
+      l = typ(x[1]) == t_POL? lgef(x[1]): 3;
+      if (lgef(x[2]) > l) l = lgef(x[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);
-	p1=p2; p2=p3;
+	y[i] = ldivres(p1,p2,&p3);
+        if (gcmp0(p3)) { i++; break; }
+	p1 = p2; p2 = p3;
       }
-      tetpil=avma; y=cgetg(i+1,t_VEC);
-      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 */
+      setlg(y, i); return gerepilecopy(av, y);
   }
-  return y;
+  err(typeer,"sfcont");
+  return NULL; /* not reached */
 }
 
 static GEN
 sfcont2(GEN b, GEN x, long k)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long l1 = lg(b), tx = typ(x), i;
   GEN y,p1;
 
@@ -1846,10 +2118,49 @@ sfcont2(GEN b, GEN x, long k)
   return gerepilecopy(av,y);
 }
 
+
 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)
 {
-  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;
 
   if (! is_matvec_t(tx)) err(typeer,"pnqn");
@@ -1899,7 +2210,7 @@ bestappr_mod(GEN x, GEN A, GEN B)
   {
     case t_INTMOD:
     {
-      ulong av = avma;
+      gpmem_t av = avma;
       GEN a,b,d, t = cgetg(3, t_FRAC);
       if (! ratlift((GEN)x[2], (GEN)x[1], &a,&b,A,B)) return NULL;
       if (is_pm1(b)) return icopy_av(a, (GEN)av);
@@ -1928,7 +2239,7 @@ bestappr_mod(GEN x, GEN A, GEN B)
 GEN
 bestappr(GEN x, GEN k)
 {
-  ulong av = avma, tetpil;
+  gpmem_t av = avma, tetpil;
   long tx,tk=typ(k),lx,i,e;
   GEN p0,p1,p,q0,q1,q,a,y;
 
@@ -1939,7 +2250,7 @@ bestappr(GEN x, GEN k)
     k = gcvtoi(k,&e);
   }
   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)
   {
     case t_INT:
@@ -1954,8 +2265,8 @@ bestappr(GEN x, GEN k)
       }
 
     case t_REAL:
-      p1=gun; p0=gfloor(x); q1=gzero; q0=gun; a=p0;
-      while (gcmp(q0,k)<=0)
+      p1=gun; a=p0=gfloor(x); q1=gzero; q0=gun;
+      while (cmpii(q0,k)<=0)
       {
 	x = gsub(x,a);
 	if (gcmp0(x)) { p1=p0; q1=q0; break; }
@@ -1980,7 +2291,7 @@ bestappr(GEN x, GEN k)
 GEN
 bestappr0(GEN x, GEN a, GEN b)
 {
-  ulong av;
+  gpmem_t av;
   GEN t;
   if (!b) return bestappr(x,a);
   av = avma;
@@ -2027,7 +2338,8 @@ update_f(GEN f, GEN a)
 GEN
 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;
 
   if (typ(x) != t_INT) err(arither1);
@@ -2056,9 +2368,9 @@ fundunit(GEN x)
   pol = quadpoly(x); y = 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);
-  if (signe(y[3]) < 0) { tetpil=avma; y=gneg(y); }
-  return gerepile(av,tetpil,y);
+  u1=gconj(y); y=gdiv(v1,u1);
+  if (signe(y[3]) < 0) y = gneg(y);
+  return gerepileupto(av, y);
 }
 
 GEN
@@ -2070,7 +2382,8 @@ gregula(GEN x, long prec)
 GEN
 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);
 
   if (signe(x)<=0) err(arither2);
@@ -2079,7 +2392,7 @@ regula(GEN x, long prec)
   rsqd = gsqrt(x,prec);
   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);
   u = stoi(r); v = gdeux;
   for(;;)
@@ -2103,7 +2416,7 @@ regula(GEN x, long prec)
   y = mplog(divri(reg,v));
   if (rexp)
   {
-    u1 = mulsr(rexp, glog(gdeux, prec));
+    u1 = mulsr(rexp, mplog2(prec));
     setexpo(u1, expo(u1)+1);
     y = addrr(y,u1);
   }
@@ -2174,7 +2487,7 @@ end_classno(GEN h, GEN hin, GEN forms, long lform)
   q = ground(gdiv(hin,h)); /* approximate order of G/H */
   for (i=1; i < lform; i++)
   {
-    long av = avma; 
+    gpmem_t av = avma; 
     fg = powgi((GEN)forms[i], h);
     fh = powgi(fg, q);
     a = (GEN)fh[1];
@@ -2198,6 +2511,59 @@ to_form(GEN x, long f)
   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
 two_rank(GEN x)
 {
@@ -2215,7 +2581,11 @@ two_rank(GEN x)
 }
 
 #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.
  * Assumes G is not too far from being cyclic.
@@ -2224,19 +2594,23 @@ two_rank(GEN x)
 GEN
 classno(GEN x)
 {
-  long av = avma, av2,c,lforms,k,l,i,j,j1,j2,lim,com,s, forms[MAXFORM];
-  long r2;
+  gpmem_t av = avma, av2, lim;
+  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 __x;
+  GEN Hf,D,fa;
   byteptr p = diffptr;
+  GEN __x;
 
   if (typ(x) != t_INT) err(arither1);
   s=signe(x); if (s>=0) return classno2(x);
 
   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));
   p2 = gtrunc(shiftr(gsqrt(p2,DEFAULTPREC),1));
   s = 1000;
@@ -2245,11 +2619,12 @@ classno(GEN x)
     if (is_bigint(p2)) err(talker,"discriminant too big in classno");
     s = itos(p2); if (s < 1000) s = 1000;
   }
-  r2 = two_rank(x);
+  r2 = two_rank(D);
   c=lforms=0;
   while (c <= s && *p)
   {
-    c += *p++; k = krogs(x,c);
+    NEXT_PRIME_VIADIFF(c,p);
+    k = krogs(D,c);
     if (!k) continue;
 
     av2 = avma;
@@ -2270,7 +2645,7 @@ classno(GEN x)
   tabla = new_chunk(10000);
   tablb = 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);
   for (i=0; i<s; i++)
   {
@@ -2300,8 +2675,9 @@ classno(GEN x)
           h = addii(addis(h,j2), mulss(s,com));
           forms[0] = (long)f;
           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 = mulii(h,Hf);
           return gerepileuptoint(av, shifti(h, r2));
         }
       }
@@ -2316,96 +2692,57 @@ classno(GEN x)
 GEN
 classno2(GEN x)
 {
-  long av=avma,tetpil,n,i,k,s=signe(x),fl2;
-  GEN p1,p2,p3,p4,p5,p7,p8,pi4,reg,logd,d,fd;
+  gpmem_t av = avma;
+  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 (!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];
-  n=lg(p1); d=gun;
-  for (i=1; i<n; i++)
-    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)); }
-  }
+  Hf = conductor_part(x, &D, &reg, &p1);
+  if (s < 0 && cmpis(D,-12) >= 0)
+    return gerepilecopy(av, Hf); /* |D| < 12*/
 
-  pi4 = mppi(DEFAULTPREC); logd = glog(d,DEFAULTPREC);
-  p1 = gsqrt(gdiv(gmul(d,logd),gmul2n(pi4,1)),DEFAULTPREC);
-  p4 = divri(pi4,d); p7 = ginv(mpsqrt(pi4));
+  Pi = mppi(DEFAULTPREC);
+  d = absi(D); logd = glog(d,DEFAULTPREC);
+  p1 = mpsqrt(gdiv(gmul(d,logd), gmul2n(Pi,1)));
   if (s > 0)
   {
-    reg = regula(d,DEFAULTPREC);
-    p2 = gsubsg(1, gmul2n(gdiv(glog(reg,DEFAULTPREC),logd),1));
-    p3 = gsqrt(gdivsg(2,logd),DEFAULTPREC);
-    if (gcmp(p2,p3)>=0) p1 = gmul(p2,p1);
+    p2 = subsr(1, gmul2n(divrr(mplog(reg),logd),1));
+    if (gcmp(gsqr(p2), divsr(2,logd)) >= 0) p1 = gmul(p2,p1);
   }
-  else reg = NULL; /* for gcc -Wall */
-  p1 = gtrunc(p1); n=p1[2];
-  if (lgefint(p1)!=3 || n<0)
-    err(talker,"discriminant too large in classno");
+  p1 = gtrunc(p1);
+  if (is_bigint(p1)) err(talker,"discriminant too large in classno");
+  n = itos(p1);
 
+  p4 = divri(Pi,d); p7 = ginv(mpsqrt(Pi));
   p1 = gsqrt(d,DEFAULTPREC); p3 = gzero;
   if (s > 0)
   {
     for (i=1; i<=n; i++)
     {
-      k=krogs(fd,i);
-      if (k)
-      {
-	p2 = mulir(mulss(i,i),p4);
-	p5 = subsr(1,mulrr(p7,incgam3(ghalf,p2,DEFAULTPREC)));
-	p5 = addrr(divrs(mulrr(p1,p5),i),eint1(p2,DEFAULTPREC));
-	p3 = (k>0)? addrr(p3,p5): subrr(p3,p5);
-      }
+      k = krogs(D,i); if (!k) continue;
+      p2 = mulir(mulss(i,i),p4);
+      p5 = subsr(1,mulrr(p7,incgam3(ghalf,p2,DEFAULTPREC)));
+      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);
-    if (!egalii(x,fd)) /* x != p3 */
-      p8 = gdiv(p8,ground(gdiv(regula(x,DEFAULTPREC),reg)));
   }
   else
   {
-    p1 = gdiv(p1,pi4);
+    p1 = gdiv(p1,Pi);
     for (i=1; i<=n; i++)
     {
-      k=krogs(fd,i);
-      if (k)
-      {
-	p2=mulir(mulss(i,i),p4);
-	p5=subsr(1,mulrr(p7,incgam3(ghalf,p2,DEFAULTPREC)));
-	p5=addrr(p5,divrr(divrs(p1,i),mpexp(p2)));
-	p3 = (k>0)? addrr(p3,p5): subrr(p3,p5);
-      }
+      k = krogs(D,i); if (!k) continue;
+      p2 = mulir(mulss(i,i),p4);
+      p5 = subsr(1,mulrr(p7,incgam3(ghalf,p2,DEFAULTPREC)));
+      p5 = addrr(p5, divrr(divrs(p1,i),mpexp(p2)));
+      p3 = (k>0)? addrr(p3,p5): subrr(p3,p5);
     }
   }
-  p3=ground(p3); tetpil=avma;
-  return gerepile(av,tetpil,mulii(p8,p3));
+  return gerepileuptoint(av, mulii(Hf,ground(p3)));
 }
 
 GEN