diff --git a/libsrc/arith.cc b/libsrc/arith.cc
index cab9414..3b819d7 100644
--- a/libsrc/arith.cc
+++ b/libsrc/arith.cc
@@ -117,7 +117,7 @@ long primeclass::number(long n)
     }
   if(!ok)
     {
-      cout<<"Not enough primes in primeclass.number("<<n<<") !"<<endl;
+      cerr<<"Not enough primes in primeclass.number("<<n<<") !"<<endl;
     }
   return p_val;
 }
@@ -135,7 +135,7 @@ vector<long> primeclass::getfirst (long n)  /* returns list of first n primes */
     }
   if(!ok)
     {
-      cout<<"Not enough primes in primeclass.getfirst("<<n<<") !"<<endl;
+      cerr<<"Not enough primes in primeclass.getfirst("<<n<<") !"<<endl;
     }
   return ans;
 }
@@ -194,16 +194,15 @@ vector<long> posdivs(long a, const vector<long>& plist)
 // cout << plist.size() << " primes: " <<endl; cout << plist << endl;
  long j,k,p,e,nd = 1;
  vector<long> dlist(1,1);  // cout << "Divisor 0 = 1" << endl;
- vector<long>::const_iterator pr = plist.begin();
- while(pr!=plist.end())
+ for (const auto& p : plist)
    {
-     p=*pr++; e = val(p,a);
+     e = val(p,a);
      dlist.resize((e+1)*dlist.size());
      for (j=0; j<e; j++)
        for (k=0; k<nd; k++)
 	 dlist[nd*(j+1)+k] = p*dlist[nd*j+k];
      nd*=(e+1);
-   } 
+   }
  return dlist;
 }
 
@@ -213,10 +212,9 @@ vector<long> alldivs(long a, const vector<long>& plist)
  long j,k,p,e,nd = 2;
  vector<long> dlist(1,1);  // cout << "Divisor 0 =  1" << endl;
  dlist.push_back(-1);      // cout << "Divisor 1 = -1" << endl;
- vector<long>::const_iterator pr = plist.begin();
- while(pr!=plist.end())
+ for (const auto& p : plist)
  {
-   p = *pr++; e = val(p,a);
+   e = val(p,a);
    dlist.resize((e+1)*dlist.size());
    for (j=0; j<e; j++)
      for (k=0; k<nd; k++)
@@ -230,10 +228,9 @@ vector<long> sqdivs(long a, const vector<long>& plist)
 {
  long j,k,p,e,nd = 1;
  vector<long> dlist(1,1);
- vector<long>::const_iterator pr = plist.begin();
- while(pr!=plist.end())
+ for (const auto& p : plist)
    {
-     p = *pr++; e = val(p,a)/2;
+     e = val(p,a)/2;
      dlist.resize((e+1)*dlist.size());
      for (j=0; j<e; j++)
        for (k=0; k<nd; k++)
@@ -247,10 +244,9 @@ vector<long> sqfreedivs(long a, const vector<long>& plist)
 {
  long j,k,p,e,nd = 1;
  vector<long> dlist(1,1);
- vector<long>::const_iterator pr = plist.begin();
- while(pr!=plist.end())
+ for (const auto& p : plist)
  {
-  p = *pr++; e = 1;
+  e = 1;
   dlist.resize((e+1)*dlist.size());
   for (j=0; j<e; j++)
     for (k=0; k<nd; k++)
@@ -507,40 +503,26 @@ int is_squarefree(long n)
   if(n%9==0) return 0;
   if(n%25==0) return 0;
   if(n%49==0) return 0;
-  vector<long>plist = pdivs(n);
-  for(unsigned int i=0; i<plist.size(); i++) 
-    if(val(plist[i],n)>1) return 0;
+  auto plist = pdivs(n);
+  for( const auto p : plist)
+    if (val(p,n)>1)
+      return 0;
   return 1;
 }
 
 int is_valid_conductor(long n)
 {
-  long m=n, p, e;
+  if (n<11) return 0;
+  long m=n, e;
   e=0; while(!(m&1)) {e++; m>>=1;}  if(e>8) return 0;
   e=0; while(!(m%3)) {e++; m/=3;}   if(e>5) return 0;
-  vector<long>plist = pdivs(m);
-  for(unsigned int i=0; i<plist.size(); i++) 
+  auto plist = pdivs(m);
+  for( const auto& p : plist)
     {
-      p = plist[i]; 
-      e=0; while(!(m%p)) {e++; m/=p;} if(e>2) return 0;
+      if (val(p,m)>2)
+        return 0;
     }
-  return 1;
-}
-
-
-
-
-// The following is no longer used
-
-long old_kronecker(long d, long n)
-{ long ans=0; long m=n, d4=d%4; if(d4<0)d4+=4;
-  if ((gcd(d,n)==1) && ((d4==0)||(d4==1)))
-  { ans=1;
-    while (!(m%4)) m/=4;
-    if (!(m%2)) {m/=2; ans*=((d-1)%8 ? -1 : 1);}
-    ans *= legendre(d,m);
-  }
-  return ans;
+      return 1;
 }
 
 /* END OF FILE */
diff --git a/libsrc/conic.cc b/libsrc/conic.cc
index e27203a..ba0f39b 100644
--- a/libsrc/conic.cc
+++ b/libsrc/conic.cc
@@ -60,19 +60,8 @@ int solve_conic(const bigint& a, const bigint& b, const bigint& c, const bigint&
                 bigint& x, bigint& y, bigint& z, int method)
 {
   vector<bigint> factorbase = pdivs(2*d);
-  //  cout<<"factorbase(1) = "<<factorbase<<endl;
-  if(is_zero(b))
-    {
-      factorbase=vector_union(factorbase,pdivs(a));
-      factorbase=vector_union(factorbase,pdivs(c));
-    }
-  else
-    {
-      bigint disc = b*b-4*a*c;
-      factorbase=vector_union(factorbase,pdivs(a));
-      factorbase=vector_union(factorbase,pdivs(disc));
-    }
-  //  cout<<"factorbase(2) = "<<factorbase<<endl;
+  factorbase=vector_union(factorbase,pdivs(a));
+  factorbase=vector_union(factorbase,pdivs((is_zero(b)? c : b*b-4*a*c)));
   return solve_conic(a,b,c,d,factorbase,x,y,z,method);
 }
 
@@ -112,7 +101,7 @@ int solve_conic(const bigint& a, const bigint& b, const bigint& c, const bigint&
       vector<bigint> aplist, bplist, cplist, dplist;
 
       int nondiag=!is_zero(b);
-      bb=a*d;  
+      bb=a*d;
       aa=-a*c; if(nondiag) aa=sqr(b)-4*aa;
       aplist=factorbase;
       bplist=factorbase;
@@ -123,7 +112,7 @@ int solve_conic(const bigint& a, const bigint& b, const bigint& c, const bigint&
       if(solve_conic_diag(a1,aplist,b1,bplist,x,y,z,method))
 	{
 	  conic_diag_reduce(a1,b1,x,y,z,verb);
-	  x*=(a2*b2); y*=a2; z*=b2; 
+	  x*=(a2*b2); y*=a2; z*=b2;
 	  if(nondiag) x-=b*z;
 	  y *= a;
 	  z *= a; if(nondiag) zz<<=1;
@@ -138,7 +127,7 @@ int solve_conic(const bigint& a, const bigint& b, const bigint& c, const bigint&
 	}
     }
 }
-  
+
 int solve_conic_diag_nontriv(const bigint& a, const vector<bigint>& aplist,
 			     const bigint& b, const vector<bigint>& bplist,
 			     bigint& x, bigint& y, bigint& z,
@@ -529,56 +518,22 @@ int testlocsol(const bigint& a, const vector<bigint>& alist,
 // coprime and square-free, their prime factors being in alist etc.
 {
   int as=sign(a), bs=sign(b), cs=sign(c);
-  if((as==bs)&&(bs==cs)) 
-    {
-      //      cout<<"testlocsol("<<a<<","<<b<<","<<c<<") returning 0 because of signs\n";
-      return 0;
-    }  
-  bigint p, two; two=2;
-  bigint mab=-a*b; 
-  vector<bigint>::const_iterator pr;
-  pr=clist.begin(); 
-  while(pr!=clist.end())
-    {
-      p=*pr++;
-      if(p==two) continue;
-      if(legendre(mab,p)!=1) 
-	{
-// 	  cout<<"testlocsol fails legendre(mab,p) with "
-// 	      <<"(a,b,p)=("<<a<<","<<b<<","<<p<<")\n"; 
-	  return 0;
-	}
-    }
-  bigint mbc=-b*c; 
-  pr=alist.begin();
-  while(pr!=alist.end())
-    {
-      p=*pr++;
-      if(p==two) continue;
-      if(legendre(mbc,p)!=1)
-	{
-// 	  cout<<"testlocsol fails legendre(mbc,p) with "
-// 	      <<"(b,c,p)=("<<b<<","<<c<<","<<p<<")\n"; 
-	  return 0;
-	}
-    }
-  bigint mca=-c*a;
-  pr=blist.begin();
-  while(pr!=blist.end())
-    {
-      p=*pr++;
-      if(p==two) continue;
-      if(legendre(mca,p)!=1)
-	{
-// 	  cout<<"testlocsol fails legendre(mca,p) with "
-// 	      <<"(c,a,p)=("<<c<<","<<a<<","<<p<<")\n"; 
-	  return 0;
-	}
-    }
+  if((as==bs)&&(bs==cs))
+    return 0;
+  auto test = [](const vector<bigint>& plist, const bigint& m)
+              {
+                return !std::any_of(plist.begin(), plist.end(),
+                                    [m](const auto& p)
+                                    {
+                                      return (p>2) && (legendre(m,p)!=1);
+                                    }
+                                    );
+              };
+  return test(clist, -a*b) && test(alist, -b*c) && test(blist, -a*c);
   return 1;
-} 
+}
 
-int testlocsol(const bigint& a, const vector<bigint>& alist, 
+int testlocsol(const bigint& a, const vector<bigint>& alist,
 	       const bigint& b, const vector<bigint>& blist)
 // tests if ax^2+by^2=z^2 is soluble, where a, b are
 // square-free, their prime factors being in alist and blist.
@@ -589,20 +544,17 @@ int testlocsol(const bigint& a, const vector<bigint>& alist,
   vector<bigint> a0list, b0list, clist;
 
   long sa=sign(a), sb=sign(b);
-  if((sa<0)&&(sb<0)) 
+  if((sa<0)&&(sb<0))
     {
 //       cout<<"testlocsol("<<a<<","<<b<<") returning 0 because of signs\n";
       return 0;  // nothing more to do as no real solution
     }
   if(sa<0) ::negate(a0);
   if(sb<0) ::negate(b0);
-  
-  vector<bigint>::const_iterator pr;
-  pr=alist.begin(); 
-  while(pr!=alist.end())
+
+  for (const auto& p : alist)
     {
-      p=*pr++;
-      if(div(p,b)) 
+      if(div(p,b))
 	{
 	  c*=p; clist.push_back(p);
 	}
@@ -611,26 +563,16 @@ int testlocsol(const bigint& a, const vector<bigint>& alist,
 	  a0*=p; a0list.push_back(p);
 	}
     }
-  pr=blist.begin(); 
-  while(pr!=blist.end())
+  for (const auto& p : blist)
     {
-      p=*pr++;
-      if(!div(p,c)) {b0*=p; b0list.push_back(p);}
-    }  
-  
-#if(0)
-  if((a!=-a0*c)||(b!=-b0*c)||(abs(c)!=gcd(a,b)))
-    {
-      cout<<"Error: (a,b)=("<<a<<","<<b<<") gives a0="<<a0<<", b0="<<b0<<", c="<<c<<endl;
+      if(!div(p,c))
+        {
+          b0*=p; b0list.push_back(p);
+        }
     }
-  //  cout<<"Calling testlocsol(a,b,c) with a,b,c="<<a0<<","<<b0<<","<<c<<"\n";
-  //  cout<<"alist="<<a0list<<endl;
-  //  cout<<"blist="<<b0list<<endl;
-  //  cout<<"clist="<<clist<<endl;
-#endif
 
   // Now a=a0*c, b=b0*c, c=gcd(a,b), and a0list, b0list, clist hold their primes
-  
+
   return testlocsol(a0,a0list,b0,b0list,c,clist);
 
 }
diff --git a/libsrc/curve.cc b/libsrc/curve.cc
index 425bb03..456df79 100644
--- a/libsrc/curve.cc
+++ b/libsrc/curve.cc
@@ -29,7 +29,7 @@
 
 int valid_invariants(const bigint& c4, const bigint& c6)  
 {
-  bigint disc = c4*c4*c4; disc -= c6*c6;
+  bigint disc = c4*c4*c4 - c6*c6;
   if (sign(disc)==0) return 0;
   if (ndiv(1728,disc)) return 0;  // need c4^3-c6^2=1728D, with D|=0
   long x6= mod(c6,27);
@@ -83,10 +83,8 @@ void minimise_c4c6(const bigint& c4, const bigint& c6, const bigint& discr,
   g = gcd( c62, newdiscr );  if(is_one(g)) return;
   const vector<bigint>& p_list = pdivs(g);
 //  cout<<"g = "<<g<<endl;
-  vector<bigint>::const_iterator pr = p_list.begin();
-  while ( pr!=p_list.end() )
+  for (const auto& p : p_list)
   {
-    p = *pr++;
     d = (long)floor(val(p,g)/12.0);
 //    cout<<"With p="<<p<<", initial d="<<d<<endl;
     if (p==2)
@@ -268,40 +266,3 @@ void Curve::tex_print(ostream &os) const
 
 
 // end of file: curve.cc
-
-//
-// The following functions from OM's code are never used
-//
-// /* number of solutions to y^2+ay=b mod p */
-// long quadroots(const bigint& a, const bigint& b, long p ) 
-// {
-//         if (p == 2) { 
-//            if (!odd(a))  return(1) ; 
-//            if (!odd(b))  return(2) ;
-//            return(0) ;
-//                 } 
-//         else {
-//          long d = I2long((a*a+4*b)%p);
-//          if (d==0) return 1;
-//          if (legendre(d,p)==1) return 2; 
-//          return 0;
-//               }
-// }
-
-// // for finding number of points mod 2 and 3
-// long pointsmod(long p, const Curve& E)
-// {
-//   bigint a1,a2,a3,a4,a6;
-//   E.getai(a1,a2,a3,a4,a6);
-//   if (p == 2)
-//     return( 1 + quadroots(a3, a6, 2) +
-//            quadroots(a1 + a3, 1 + a2 + a4 + a6, 2)) ;
-//   if (p ==3)
-//     return( 1 + quadroots(a3, a6, 3) +
-//            quadroots(a1 + a3, 1 + a2 + a4 + a6, 3)
-//            + quadroots(-a1 + a3, -1 + a2 - a4 + a6, 3)) ;
-//  // Now p>3
-//   long count=0;
-//   for (long x=0;x<p;x++) count+=quadroots(x*a1+a3,x*(x*(x+a2)+a4)+a6, p);
-//   return  count;
-// }
diff --git a/libsrc/curvedata.cc b/libsrc/curvedata.cc
index cbd058c..2e424e6 100644
--- a/libsrc/curvedata.cc
+++ b/libsrc/curvedata.cc
@@ -67,9 +67,8 @@ Curvedata::Curvedata(const vector<bigrational>& qai, bigint& scale)
 #ifdef DEBUG_Q_INPUT
   cout<<"Denominator primes: "<<plist<<endl;
 #endif
-  for(vector<bigint>::const_iterator pp=plist.begin(); pp!=plist.end(); pp++)
+  for( const auto& p : plist)
     {
-      bigint p=*pp;
 #ifdef DEBUG_Q_INPUT
       cout<<"p =  "<<p<<endl;
 #endif
@@ -225,11 +224,10 @@ void Curvedata::minimalize()
 	{
 	  bigint p;
 	  vector<bigint> new_bad_primes;
-	  vector<bigint>::iterator pr = the_bad_primes.begin();
-	  while(pr!=the_bad_primes.end())
+          for (const auto& p : the_bad_primes)
 	    {
-	      bigint p=*pr++;
-	      if(div(p,discr)) new_bad_primes.push_back(p);
+	      if(div(p,discr))
+                new_bad_primes.push_back(p);
 	    }
 	  the_bad_primes=new_bad_primes;
 	}
@@ -314,7 +312,6 @@ void Curvedata::output(ostream& os) const
   else os<<"#torsion not yet computed"<<endl;
 }
 
-
 Curvedata opt_x_shift(const Curvedata& C, bigint& k)
 {
   bigint b2,b4,b6,b8,four,zero; four=BIGINT(4); zero=BIGINT(0);
@@ -326,86 +323,4 @@ Curvedata opt_x_shift(const Curvedata& C, bigint& k)
   return CD;
 }
 
-#if(0)
-
-// next the invariants package
-
-bigfloat agm( bigfloat a, bigfloat b)
-{
-  bigfloat a0, b0 ;
-  do{
-    a0 = (a + b) / 2.0 ;
-    b0 = sqrt( a * b ) ;
-    a = a0 ; b = b0 ;
-  } while ( !is_approx_zero(a-b) ) ;
-  return ( a ) ;
-}
- 
-bigfloat cuberoot(bigfloat x ) /* computes x^(1.0/3.0) for all x */
-{ 
-  bigfloat third(1); third/=3;
-  if (x >= 0.0) 
-    return (pow(x, third) ) ;
-  else 
-    return( - pow(-x, third) );
-}
-
-void  sort(bigfloat roots[3] )
-{
-        bigfloat x0, x1, x2, temp ;
-        x0 = roots[0] ; x1 = roots[1] ; x2 = roots[2] ;
-        if (x0 > x1) { temp = x0; x0 = x1 ; x1 = temp ; }  /* now x0 < x1 */
-        if (x1 > x2) { temp = x1; x1 = x2 ; x2 = temp ; }  /* now x1 < x2, x0 < x2 */
-        if (x0 > x1) { temp = x0; x0 = x1 ; x1 = temp ; }  /* now x0 < x1 < x2 */
-        roots[0] = x0 ; roots[1] = x1; roots[2] = x2 ;
-}
-
-CurvedataExtra::CurvedataExtra(const Curvedata& c)
-: Curvedata(c)
-{
-  if ( isnull() ) {roots[0]=roots[1]=roots[2]=period=0; return; }
-  bigfloat fb2 = I2bigfloat(b2);
-  bigfloat fc6 = I2bigfloat(c6);
-  bigfloat fdiscr = I2bigfloat(discr);
-  if (conncomp == 1) {
-    bigfloat sqdisc  = 24.0 * sqrt(-3.0 * fdiscr) ;
-    roots[2] = ( -fb2 + cuberoot(fc6 - sqdisc) +
-                cuberoot(fc6 + sqdisc) ) / 12.0 ;
-  } else {
-    bigfloat theta = atan2( 24.0 * sqrt(3.0*fdiscr),  fc6) ;
-    bigfloat fc4 = I2bigfloat(c4);
-    bigfloat sqc4 = sqrt( fc4 ) ;
-    roots[0] = (-fb2 + 2.0 * sqc4 * cos( theta/3.0 )) / 12.0 ;
-    roots[1] = (-fb2 + 2.0 * sqc4 * cos( (theta+2.0*Pi())/3.0)) / 12.0;
-    roots[2] = (-fb2 + 2.0 * sqc4 * cos( (theta+4.0*Pi())/3.0)) / 12.0;
-    sort(roots) ;
-  }
-  // note that roots[2] (the THIRD one!) is always the largest root 
-  // now compute the period
-  bigfloat a, b ; // parameters to feed to agm()
-  if (conncomp == 2){
-    a = sqrt(roots[2] - roots[0]) ;
-    b = sqrt(roots[2] - roots[1]) ;
-  } else {
-    bigfloat fb4 = I2bigfloat(b4); // bigfloat version of b4
-    bigfloat fp = sqrt( fb4/2.0 + roots[2] * (fb2/2.0 + 3.0 * roots[2]) );
-                                  // f', where f is the cubic 
-    a = 2.0 * sqrt( fp ) ;
-    b = sqrt( 3.0 * roots[2] + fb2/4.0 + 2.0 * fp ) ;
-  }
-  period =  (2.0 * Pi() / agm(a, b) );
-}
-
-void CurvedataExtra::output(ostream& os) const
-{
-  Curvedata::output(os);
-  if(conncomp == 1) 
-     os << "The real two-division point is " << roots[2] << endl ;
-  else 
-  if(conncomp == 2) 
-     os << "The real two-division points are " 
-        << roots[0] << ", " << roots[1] << ", "<< roots[2] << endl ; 
-  os << "The real period is " << period << endl;
-}
-#endif
 
diff --git a/libsrc/curvered.cc b/libsrc/curvered.cc
index a7875de..dafbfbc 100644
--- a/libsrc/curvered.cc
+++ b/libsrc/curvered.cc
@@ -115,7 +115,7 @@ void CurveRed::operator=(const CurveRed& E)
   conncomp=E.conncomp;
   ntorsion=E.ntorsion;
   reduct_array = E.reduct_array;
-  N=E.N; 
+  N=E.N;
 }
 
 CurveRed::CurveRed(const Curvedata& E)
@@ -136,10 +136,9 @@ CurveRed::CurveRed(const Curvedata& E)
   // Because the curve is minimal, Tate's algorithm reduce-loop is not needed
 
 //cout<<"Running Tate's algorithm"<<endl;
-  vector<bigint>::const_iterator pi= the_bad_primes.begin();
-  while(pi!=the_bad_primes.end())
+
+  for (const auto& p : the_bad_primes)
     {
-    p = *pi++;
     ord_p_discr = val(p,discr);
     ord_p_j = ord_p_discr - 3*val(p,c4);
     if (ord_p_j < 0) ord_p_j = 0;
@@ -166,7 +165,7 @@ CurveRed::CurveRed(const Curvedata& E)
     r = mod(r,p);
     t = mod(t,p);
     C.transform(r,zero,t);
-    
+
     // test for Types In, II, III, IV
     if ( ndiv(p,c4) )
       {temp=-C.a2;
@@ -301,7 +300,7 @@ CurveRed::CurveRed(const Curvedata& E)
 	  reduct_array[p] = Reduction_type
 	    (ord_p_discr, ord_p_discr - 6, ord_p_j, 5, c_p);
 	  break; }
-      
+
       // change coords so that p^3|C.a3, p^5|C.a6
       if ( pdiv2 ) t = -p*p*root(a6t,2,p);
       else t = p*p*mod(-a3t*halfmodp, p);
@@ -325,11 +324,8 @@ CurveRed::CurveRed(const Curvedata& E)
     };  // end switch
   }     // end primes for-loop
   N = BIGINT(1);
-  map<bigint,Reduction_type>::const_iterator ri; 
-  for(ri = reduct_array.begin(); ri!=reduct_array.end(); ri++)
-    {
-      N *= pow((ri->first), (ri->second).ord_p_N);
-    }
+  for( const auto& ri : reduct_array)
+    N *= pow((ri.first), (ri.second).ord_p_N);
   return;
 }     // end of Tate's algorithm
 
@@ -338,7 +334,7 @@ CurveRed::CurveRed(const Curvedata& E)
 bigint local_Tamagawa_exponent(CurveRed& c, const bigint& p)
 {
   if (is_zero(p)) return BIGINT(c.conncomp);
-  map<bigint,Reduction_type>::const_iterator ri = c.reduct_array.find(p);
+  auto ri = c.reduct_array.find(p);
   if (ri == c.reduct_array.end())
     return BIGINT(1);
   Reduction_type info = ri->second;
@@ -362,9 +358,9 @@ bigint global_Tamagawa_exponent(CurveRed& c, int real_too)
   const bigint two = BIGINT(2);
   bigint ans = ((real_too && (c.conncomp==2))? two: one);
 
-  for(map<bigint,Reduction_type>::const_iterator ri = c.reduct_array.begin(); ri!=c.reduct_array.end(); ri++)
+  for( const auto&  ri : c.reduct_array)
     {
-      Reduction_type info = ri->second;
+      Reduction_type info = ri.second;
       int code = info.Kcode.code;
       int ep = (code%20==1? 2: info.c_p); // Type I*m, m even: [2,2]
       ans = lcm(ans,BIGINT(ep));
@@ -377,8 +373,8 @@ vector<long> tamagawa_primes(CurveRed& C, int real_too)
 {
   vector<bigint> T = pdivs(global_Tamagawa_exponent(C, real_too));
   vector<long> t;
-  for(vector<bigint>::const_iterator ti = T.begin(); ti!=T.end(); ++ti)
-    t.push_back(I2long(*ti));
+  for( const auto&  ti : T)
+    t.push_back(I2long(ti));
   return t;
 }
 
@@ -393,28 +389,28 @@ vector<long> tamagawa_primes(CurveRed& C, int real_too)
 
 int getord_p_discr(const CurveRed& c, const bigint& p)
 {
-  map<bigint,Reduction_type>::const_iterator ri = c.reduct_array.find(p);
+  auto ri = c.reduct_array.find(p);
   if(ri==c.reduct_array.end()) return 0;
   return (ri->second).ord_p_discr;
 }
 
 int getord_p_N(const CurveRed& c, const bigint& p)
 {
-  map<bigint,Reduction_type>::const_iterator ri = c.reduct_array.find(p);
+  auto ri = c.reduct_array.find(p);
   if(ri==c.reduct_array.end()) return 0;
   return (ri->second).ord_p_N;
 }
 
 int getord_p_j_denom(const CurveRed& c, const bigint& p)
 {
-  map<bigint,Reduction_type>::const_iterator ri = c.reduct_array.find(p);
+  auto ri = c.reduct_array.find(p);
   if(ri==c.reduct_array.end()) return 0;
   return (ri->second).ord_p_j_denom;
 }
 
 int getc_p(const CurveRed& c, const bigint& p)
 {
-  map<bigint,Reduction_type>::const_iterator ri = c.reduct_array.find(p);
+  auto ri = c.reduct_array.find(p);
   if(ri==c.reduct_array.end()) return 1;
   return (ri->second).c_p;
 }
@@ -422,11 +418,8 @@ int getc_p(const CurveRed& c, const bigint& p)
 bigint prodcp(const CurveRed& c)
 {
   bigint ans = BIGINT(1);
-  map<bigint,Reduction_type>::const_iterator ri; 
-  for(ri = c.reduct_array.begin(); ri!=c.reduct_array.end(); ri++)
-    {
-      ans *= (ri->second).c_p;
-    }
+  for( const auto& ri : c.reduct_array)
+    ans *= (ri.second).c_p;
   return ans;
 }
 
@@ -444,7 +437,7 @@ bigint global_Tamagawa_number(CurveRed& c, int real_too)
 
 Kodaira_code getKodaira_code(const CurveRed& c, const bigint& p)
 {
-  map<bigint,Reduction_type>::const_iterator ri = c.reduct_array.find(p);
+  auto ri = c.reduct_array.find(p);
   if(ri==c.reduct_array.end()) return Kodaira_code();
   return (ri->second).Kcode;
 }
@@ -474,12 +467,11 @@ void CurveRed::display(ostream& os)
   os << "Global Root Number = " << GlobalRootNumber(*this) << endl;
   os << "Reduction type at bad primes:\n";
   os <<"p\tord(d)\tord(N)\tord(j)\tKodaira\tc_p\troot_number\n";
-  map<bigint,Reduction_type>::const_iterator ri;  
-  for(ri = reduct_array.begin(); ri!=reduct_array.end(); ri++)
+  for( const auto& ri : reduct_array)
     {
-      if((ri->second).local_root_number==0) 
-	setLocalRootNumber(ri->first);     
-      os << (ri->first) << "\t" << (ri->second) << endl;
+      if((ri.second).local_root_number==0)
+	setLocalRootNumber(ri.first);
+      os << (ri.first) << "\t" << (ri.second) << endl;
     }
 }
 
@@ -503,9 +495,9 @@ void CurveRed::display(ostream& os)
 int LocalRootNumber(CurveRed& c, const bigint& p)
 {
   if(is_zero(p)) return -1;  // the infinite prime
-  map<bigint,Reduction_type>::const_iterator ri = c.reduct_array.find(p);
+  auto ri = c.reduct_array.find(p);
   if(ri==c.reduct_array.end()) return 1; // good reduction case
-  if((ri->second).local_root_number==0) 
+  if((ri->second).local_root_number==0)
     c.setLocalRootNumber(p);
   return (ri->second).local_root_number;
 }
@@ -513,12 +505,11 @@ int LocalRootNumber(CurveRed& c, const bigint& p)
 int GlobalRootNumber(CurveRed& c)
 {
   int ans=-1;
-  map<bigint,Reduction_type>::const_iterator ri; 
-  for(ri = c.reduct_array.begin(); ri!=c.reduct_array.end(); ri++)
+  for( const auto& ri : c.reduct_array)
     {
-      if((ri->second).local_root_number==0)  
-	c.setLocalRootNumber(ri->first);
-      ans *= (ri->second).local_root_number;
+      if((ri.second).local_root_number==0)
+	c.setLocalRootNumber(ri.first);
+      ans *= (ri.second).local_root_number;
     }
   return ans;
 }
@@ -537,7 +528,7 @@ int kro_3(long x); // kronecker(x,3)
 void CurveRed::setLocalRootNumber2()
 {
   static const bigint two = BIGINT(2);
-  map<bigint,Reduction_type>::iterator ri = reduct_array.find(two);
+  auto ri = reduct_array.find(two);
   if(ri==reduct_array.end()) return;
   Reduction_type& info = ri->second;
   int kod = PariKodairaCode(info.Kcode);
@@ -734,7 +725,7 @@ void CurveRed::setLocalRootNumber2()
 void CurveRed::setLocalRootNumber3()
 {
   static const bigint three = BIGINT(3);
-  map<bigint,Reduction_type>::iterator ri = reduct_array.find(three);
+  auto ri = reduct_array.find(three);
   if(ri==reduct_array.end()) return;
   Reduction_type& info = ri->second;
   int kod = PariKodairaCode(info.Kcode);
@@ -841,7 +832,7 @@ void CurveRed::setLocalRootNumber3()
 
 void CurveRed::setLocalRootNumber_not_2_or_3(const bigint& p)
 {
-  map<bigint,Reduction_type>::iterator ri = reduct_array.find(p);
+  auto ri = reduct_array.find(p);
   if(ri==reduct_array.end()) return;
   Reduction_type& info = ri->second;
   if (info.ord_p_N == 1) 
diff --git a/libsrc/curvesort.cc b/libsrc/curvesort.cc
index 38fec8a..a8302fb 100644
--- a/libsrc/curvesort.cc
+++ b/libsrc/curvesort.cc
@@ -267,9 +267,8 @@ int booknumber(int level, int form)  // permutes numbers starting from 1
 int codeletter_to_int(string code)  // i counts from 0!
 {
   int n=0;
-  string::iterator c;
-  for (c=code.begin(); c<code.end(); c++)
-    n = 26*n + ((*c)-'a');
+  for ( const auto& c : code)
+    n = 26*n + (c-'a');
   return n;
   // int b = code[0]-'a';
   // if(code[1]=='\0')  return b;
diff --git a/libsrc/eclib/arith.h b/libsrc/eclib/arith.h
index 36c9daf..f6204b4 100644
--- a/libsrc/eclib/arith.h
+++ b/libsrc/eclib/arith.h
@@ -97,7 +97,7 @@ class primevar {
    for(pr.init(n); pr.ok(); pr++) {p=pr; ... ;}  // iff pr is existing primevar
 
 ---To loop through all primes:
-   
+
    primevar pr; //or for existing primevar, pr.init();
    long p;
    while(pr.ok()) {p=pr; pr=++; ...;}
@@ -110,29 +110,25 @@ vector<long> pdivs(long);      /* list of prime divisors */
 vector<long> posdivs(long, const vector<long>& plist);  // all positive divisors
 inline vector<long> posdivs(long n)
 {
-  const vector<long>& plist = pdivs(n);
-  return posdivs(n,plist);
+  return posdivs(n, pdivs(n));
 }
 
 vector<long> alldivs(long, const vector<long>& plist);  // absolutely all divisors
 inline vector<long> alldivs(long n)
 {
-  const vector<long>& plist = pdivs(n);
-  return alldivs(n,plist);
+  return alldivs(n, pdivs(n));
 }
 
 vector<long> sqdivs(long, const vector<long>& plist);   //  divisors whose square divides
 inline vector<long> sqdivs(long n)
 {
-  const vector<long>& plist = pdivs(n);
-  return sqdivs(n,plist);
+  return sqdivs(n, pdivs(n));
 }
 
 vector<long> sqfreedivs(long, const vector<long>& plist); // square-free divisors
 inline vector<long> sqfreedivs(long n)
 {
-  const vector<long>& plist = pdivs(n);
-  return sqfreedivs(n,plist);
+  return sqfreedivs(n, pdivs(n));
 }
 
 
diff --git a/libsrc/elog.cc b/libsrc/elog.cc
index 5a131ac..c652a9f 100644
--- a/libsrc/elog.cc
+++ b/libsrc/elog.cc
@@ -308,13 +308,11 @@ vector<Point> division_points_by2(Curvedata& E,  const Point& P)
   cout<<"Possible x-coordinates:"<<xans<<endl;
 #endif
   vector<Point> ans;
-  for(vector<bigrational>::const_iterator x=xans.begin(); x!=xans.end(); x++)
+  for( const auto & x : xans)
     {
-      vector<Point> x_points = points_from_x(E,*x);
-      for(vector<Point>::const_iterator Qi=x_points.begin(); 
-          Qi!=x_points.end(); Qi++)
+      vector<Point> x_points = points_from_x(E,x);
+      for( const auto& Q : x_points)
         {
-          Point Q = *Qi;
           if(2*Q==P) // as it might = -P
             {
 #ifdef DEBUG_DIVBY2
diff --git a/libsrc/fixc6.cc b/libsrc/fixc6.cc
index 95bd08c..1706b67 100644
--- a/libsrc/fixc6.cc
+++ b/libsrc/fixc6.cc
@@ -57,11 +57,11 @@ fixc6::fixc6()
 void fixc6::operator()(long N, int i, bigint& c4, bigint& c6)
 {
   pair<long,int> key(N,i+1);
-  map<pair<long,int>,bigint>::const_iterator j = fixc6table.find(key); 
+  auto j = fixc6table.find(key);
   if(j!=fixc6table.end()) c6=j->second;
-  j = fixc4table.find(key); 
+  j = fixc4table.find(key);
   if(j!=fixc4table.end()) c4=j->second;
-  return; 
+  return;
 }
 
 map< pair<long,int>, bigint > fixc6::fixc6table;
diff --git a/libsrc/gf.cc b/libsrc/gf.cc
index ff4a530..7aa0920 100644
--- a/libsrc/gf.cc
+++ b/libsrc/gf.cc
@@ -26,17 +26,17 @@
 map<ZZ,ZZ_pContext> ZZ_pContextCache;
 
 galois_field::galois_field(void)  //dummy
- :q(to_ZZ(2))  
+ :q(to_ZZ(2))
 {
   //  cout<<"In galois_field constructor, calling default ZZ_p::init(2)"<<endl;
   ZZ_p::init(q);
 }
 
-galois_field::galois_field(const ZZ& qq) 
- :q(qq) 
+galois_field::galois_field(const ZZ& qq)
+ :q(qq)
 {
   //  cout<<"In galois_field constructor with q="<<qq<<endl;
-  map<ZZ,ZZ_pContext>::iterator t = ZZ_pContextCache.find(q);
+  auto t = ZZ_pContextCache.find(q);
   if(t==ZZ_pContextCache.end())
     {
       //      cout<<"Calling ZZ_p::init("<<q<<")"<<endl;
diff --git a/libsrc/hilbert.cc b/libsrc/hilbert.cc
index 9d74367..d419e1a 100644
--- a/libsrc/hilbert.cc
+++ b/libsrc/hilbert.cc
@@ -79,10 +79,8 @@ int global_hilbert(const bigint& a, const bigint& b, const vector<bigint>& plist
 #ifdef DEBUG_HILBERT
   cout<<"Passed local condition at infinity..."<<endl;
 #endif
-  vector<bigint>::const_iterator pr = plist.begin();
-  while(pr!=plist.end())
+  for (const auto& badp : plist)
     {
-      badp=*pr++;
 #ifdef DEBUG_HILBERT
       cout<<"Testing local condition at "<<badp<<"..."<<endl;
 #endif
diff --git a/libsrc/homspace.cc b/libsrc/homspace.cc
index c0d7b6a..c3e4469 100644
--- a/libsrc/homspace.cc
+++ b/libsrc/homspace.cc
@@ -366,19 +366,16 @@ if (verbose>1)
 
    // Compute the number of cusps
    long maxncusps =0, dd, pp, nc;
-   vector<long>::const_iterator d,p;
-   for(d=(dlist).begin();d!=(dlist).end();d++)
+   for( const auto& d : dlist)
      {
-       dd = *d;
-       nc = ::gcd(dd,modulus/dd);
-       for(p=plist.begin();p!=plist.end();p++) // computing phi(dd)
+       nc = ::gcd(dd,modulus/d);
+       for( const auto& p : plist) // computing phi(dd)
          {
-           pp = *p;
-           if ((nc%pp)==0) 
-             nc = nc*(pp-1)/pp;
+           if ((nc%p)==0)
+             nc = nc*(p-1)/p;
          }
        maxncusps += nc;
-     }  
+     }
    if (verbose) cout << "Number of cusps is "<<maxncusps<<endl;
 
    cusplist cusps(maxncusps, this);
@@ -507,7 +504,7 @@ vec homspace::extend_coords(const vec& v)
   vec ans(nsymb);
   int i,j;
   for(i=1; i<=nsymb; i++)
-    {      
+    {
       j = coordindex[i-1];
       if (j==0) ans[i] = 0;
       else if (j>0) ans[i] =  v*coord_vecs[j];
diff --git a/libsrc/htconst.cc b/libsrc/htconst.cc
index 5577ffb..2a5ea63 100644
--- a/libsrc/htconst.cc
+++ b/libsrc/htconst.cc
@@ -107,22 +107,21 @@ double cps_bound(const Curvedata& CD)
 #endif
   CurveRed CR(CD);
   vector<bigint> plist = getbad_primes((Curvedata&)CD);
-  for(unsigned int i=0; i<plist.size(); i++)
+  for( const auto& q : plist)
     {
-      bigint q=plist[i];
-      if(getc_p(CR,q)==1) 
+      if(getc_p(CR,q)==1)
 	{
 #ifdef DEBUG_CPS
       cout<<"q = "<<q<<", alpha = 0 since c_q=0"<<endl;
       cout<<"sum so far = "<<bd<<endl;
 #endif
-	  continue;      
+	  continue;
 	}
       double alpha =0;
       int m, Kc = getKodaira_code(CR,q).code;
       switch (Kc%10){
       case 0: // Im
-	m = Kc/10; 
+	m = Kc/10;
 	alpha = (m%2 ? double(m*m-1)/double(4*m) : double(m)/4);
 	break;
       case 1: // I*m
@@ -396,10 +395,8 @@ bigfloat calc_dvd_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6
   // Samir 25/07/04
 
   int first=1;
-  std::set<bigfloat>::const_iterator xi = F_roots.begin();
-  while(xi!=F_roots.end())
+  for (const auto& x : F_roots)
     {
-      x=*xi++;
       x2=x*x;
       Gx=abs(1-b4*x2-2.0*b6*x*x2-b8*x2*x2);
 #ifdef DEBUG_CPS
@@ -411,10 +408,8 @@ bigfloat calc_dvd_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6
       cout<<"dvd so far ="<<dvd<<endl;
 #endif
     }
-  xi = crit_pts.begin();
-  while(xi!=crit_pts.end())
-    { 
-      x=*xi++;
+  for (const auto& x : crit_pts)
+    {
       x2=x*x;
       Fx=(4.0*x+b2*x2+2.0*b4*x*x2+b6*x2*x2);
 #ifdef DEBUG_CPS
@@ -526,10 +521,8 @@ bigfloat calc_dv_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6,
   // Samir 25/07/04
 
   int first=1;
-  std::set<bigfloat>::const_iterator xi = f_roots.begin();
-  while(xi!=f_roots.end())
+  for( const auto& x : f_roots)
     {
-      x=*xi++;
       x2=x*x;
       gx=abs(x2*x2-b4*x2-2.0*b6*x-b8);
 #ifdef DEBUG_CPS
@@ -541,10 +534,8 @@ bigfloat calc_dv_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6,
       cout<<"dv so far ="<<dv<<endl;
 #endif
     }
-  xi = crit_pts.begin();
-  while(xi!=crit_pts.end())
-    { 
-      x=*xi++;
+  for( const auto& x : crit_pts)
+    {
       x2=x*x;
       fx=(4.0*x*x2+b2*x2+2.0*b4*x+b6);
 #ifdef DEBUG_CPS
@@ -1021,28 +1012,22 @@ vector<bigfloat> roots11( const vector<bigfloat>& coeff )
 vector<bigfloat> reals_in ( vector<bigcomplex>& v)
 {
   vector<bigfloat> vr;
-  vector<bigcomplex>::iterator vi = v.begin();
-  while(vi!=v.end()) 
-    {
-      if(is_real(*vi)) vr.push_back((*vi).real());
-      vi++;
-    }
+  for (const auto& vi : v)
+    if(is_real(vi))
+      vr.push_back(vi.real());
   return vr;
 }
 
 vector<bigfloat> reals_in_11 ( vector<bigcomplex>& v)
 {
   vector<bigfloat> vr;
-  vector<bigcomplex>::iterator vi = v.begin();
-  while(vi!=v.end()) 
-    {
-      if(is_real(*vi)) 
-	{
-	  bigfloat x = (*vi).real();
-	  if((x<=1)&&(x>=-1)) vr.push_back(x);
-	}
-      vi++;
-    }
+  for (const auto& vi : v)
+    if(is_real(vi))
+      {
+        bigfloat x = vi.real();
+        if((x<=1)&&(x>=-1))
+          vr.push_back(x);
+      }
   return vr;
 }
 
@@ -1053,10 +1038,10 @@ int interval_test(const bigfloat& x, const vector<bigfloat> rts, int debug)
   if(x>1) {if(debug) cout<<"\t returns 0\n"; return 0;}
   if(x<-1) {if(debug) cout<<"\t returns 0\n";  return 0;}
   int ans;
-  if(rts.size()==1) 
+  if(rts.size()==1)
     ans = (x>=rts[0]);
-  else 
-    ans =  ((x>=rts[0]) && (x<=rts[1])) || (x>=rts[2]); 
+  else
+    ans =  ((x>=rts[0]) && (x<=rts[1])) || (x>=rts[2]);
   if(debug) cout<<"\t returns "<<ans<<"\n";
   return ans;
 }
@@ -1064,17 +1049,14 @@ int interval_test(const bigfloat& x, const vector<bigfloat> rts, int debug)
 vector<bigfloat> reals_in_interval ( vector<bigcomplex>& v, const vector<bigfloat> rts)
 {
   vector<bigfloat> vr;
-  vector<bigcomplex>::iterator vi = v.begin();
   bigfloat x;
-  while(vi!=v.end()) 
-    {
-      if(is_real(*vi)) 
-	{
-	  x=(*vi).real();
-	  if(interval_test(x,rts,1)) vr.push_back(x);
-	}
-      vi++;
-    }
+  for( const auto& vi : v)
+    if(is_real(vi))
+      {
+        bigfloat x=vi.real();
+        if(interval_test(x,rts,1))
+          vr.push_back(x);
+      }
   return vr;
 }
 
@@ -1083,15 +1065,13 @@ vector<bigfloat> reals_in_interval ( vector<bigcomplex>& v, const vector<bigfloa
 void include_real_11(std::set<bigfloat>& S, const vector<bigcomplex>& C)
 {
   bigfloat x;
-  vector<bigcomplex>::const_iterator Ci = C.begin();
-  while(Ci!=C.end())
+  for (const auto& Ci : C)
     {
-      if(is_real(*Ci))
+      if(is_real(Ci))
 	{
-	  x = (*Ci).real();
+	  x = Ci.real();
 	  if((x<=1)&&(x>=-1)) S.insert(x);
 	}
-      Ci++;
     }
 }
 
@@ -1423,12 +1403,12 @@ void Interval::intersect(const Interval& I)
 vector<Interval> intersect(const vector<Interval>& L1, const vector<Interval>& L2)
 {
   vector<Interval> ans;
-  vector<Interval>::const_iterator I, J;
-  for(I=L1.begin(); I!=L1.end(); I++)
-    for(J=L2.begin(); J!=L2.end(); J++)
+  for( const auto& I : L1)
+    for( const auto& J : L2)
       {
-	Interval K = intersect(*I,*J);
-	if(!K.is_empty()) ans.push_back(K);
+	Interval K = intersect(I,J);
+	if(!K.is_empty())
+          ans.push_back(K);
       }
   return ans;
 }
@@ -1459,15 +1439,14 @@ void Interval01::intersect(const Interval01& I)
   empty=(lh>rh);
 }
 
-vector<Interval01> intersect(const vector<Interval01>& L1, 
+vector<Interval01> intersect(const vector<Interval01>& L1,
 			     const vector<Interval01>& L2)
 {
   vector<Interval01> ans;
-  vector<Interval01>::const_iterator I, J;
-  for(I=L1.begin(); I!=L1.end(); I++)
-    for(J=L2.begin(); J!=L2.end(); J++)
+  for( const auto& I : L1)
+    for( const auto& J : L2)
       {
-	Interval01 K = intersect(*I,*J);
+	Interval01 K = intersect(I,J);
 	if(!K.is_empty()) ans.push_back(K);
       }
   return ans;
@@ -1554,11 +1533,9 @@ CurveHeightConst::CurveHeightConst(CurveRed& CR)
 
 long CurveHeightConst::e_p(long p)
 {
-  map<long, long>::iterator pe = ann.find(p);
+  auto pe = ann.find(p);
   if (pe!=ann.end())
-    {
-      return pe->second;
-    }
+    return pe->second;
   long e = exponent(*this,p);
   ann[p] = e;
   return e;
@@ -1566,7 +1543,7 @@ long CurveHeightConst::e_p(long p)
 
 bigfloat CurveHeightConst::D(long n) // D_E(n) in the paper
 {
-  map<long, bigfloat>::iterator DEn = DE.find(n);
+  auto DEn = DE.find(n);
   if (DEn!=DE.end())
     {
 #ifdef debugLB
@@ -1832,12 +1809,12 @@ vector<Interval> CurveHeightConst::solveLEQ(long n, const bigfloat& B)
       ans.push_back(Interval(e3,xlist[0]));
       for(int i=1; i<n; i+=2)
 	ans.push_back(Interval(xlist[i],xlist[i+1]));
-    }	
+    }
   else // n is even
     {
       for(int i=0; i<n; i+=2)
-	ans.push_back(Interval(xlist[i],xlist[i+1]));      
-    }  
+	ans.push_back(Interval(xlist[i],xlist[i+1]));
+    }
   return ans;
 }
 
@@ -1846,7 +1823,7 @@ vector<Interval> CurveHeightConst::solveLEQ(long n, const bigfloat& B)
 vector<Interval> CurveHeightConst::solveGEQ(long n, const bigfloat& B)
 {
   vector<Interval> ans;
-  if(B <= e3) 
+  if(B <= e3)
     {
       ans.push_back(Interval(e3)); // i.e.[e3,infty]
       return ans;
@@ -1858,14 +1835,14 @@ vector<Interval> CurveHeightConst::solveGEQ(long n, const bigfloat& B)
       for(i=1; i<n; i+=2)
 	ans.push_back(Interval(xlist[i-1],xlist[i]));
       ans.push_back(Interval(xlist[n-1]));
-    }	
+    }
   else // n is even
     {
       ans.push_back(Interval(e3,xlist[0]));
       for(i=2; i<n; i+=2)
-	ans.push_back(Interval(xlist[i-1],xlist[i]));      
+	ans.push_back(Interval(xlist[i-1],xlist[i]));
       ans.push_back(Interval(xlist[n-1]));
-    }  
+    }
   return ans;
 }
 
diff --git a/libsrc/isogs.cc b/libsrc/isogs.cc
index 0cddc8f..776b6a5 100644
--- a/libsrc/isogs.cc
+++ b/libsrc/isogs.cc
@@ -175,9 +175,8 @@ int semistable(const CurveRed& CR)
 {
   int ans=1;
   vector<bigint> plist = getbad_primes(CR);
-  vector<bigint>::iterator pvar = plist.begin();
-  while(pvar!=plist.end())
-    if(getord_p_N(CR,*pvar++)>1) return 0;
+  for (const auto& p : plist)
+    if(getord_p_N(CR,p)>1) return 0;
   return ans;
 }
 
@@ -271,33 +270,29 @@ void IsogenyClass::process(long i)  // process i'th curve
   CurveRed thisc = curves[i];
   if (verb) cout << "Working on curve " << i+1 << ": " << (Curve)thisc << endl;
   Cperiods pers(thisc);
-  vector<long>::iterator lvar=llist.begin();
-  long il=0, ell, n;
-  while(lvar!=llist.end())
+  long il=0, n;
+  for (const auto& ell : llist)
     {
-      ell = *lvar++;
       if (verb) cout << "trying l = " << ell << "..." << flush;
 
       vector<CurveRed> lisogcurves = lisog(thisc,pers,ell,verb);
 
       if (verb) cout << lisogcurves.size() << " isogenous curves found." << endl;
-      if(i==0) 
+      if(i==0)
 	{
 // 	  cout<<"setting lworks["<<il<<"] to "<< !lisogcurves.empty() << endl;
 // 	  cout<<"where llist["<<il<<"] = "<<llist[il]<<endl;
 	  lworks[il++] = !lisogcurves.empty();
 	}
-      vector<CurveRed>::iterator Ci=lisogcurves.begin();
       n=0;
-      while(Ci!=lisogcurves.end())
-	{ 
-	  CurveRed newc = *Ci++; n++;
+      for (const auto& newc : lisogcurves)
+	{
+	  n++;
 	  if (verb) cout << "\t"<<n<<": "<<(Curve)newc<<"\t: ";
 	  int j=0, isnew = 1;
-	  vector<CurveRed>::iterator oldCi=curves.begin();
-	  while(oldCi!=curves.end())
-	    {	      
-	      if ((Curve)newc==(Curve)(*oldCi++))
+          for (const auto& oldCi : curves)
+	    {
+	      if ((Curve)newc==(Curve)(oldCi))
 		{
 		  isnew=0;
 		  matset(i,j,ell);
@@ -305,7 +300,7 @@ void IsogenyClass::process(long i)  // process i'th curve
 		}
 	      j++;
 	    }
-	  if (isnew) 
+	  if (isnew)
 	    {
 	      curves.push_back(newc);
 	      fromlist.push_back(i);
@@ -351,12 +346,14 @@ void IsogenyClass::displaycurves(ostream& os)const
 {
   os << endl << ncurves << " curve(s) in the isogeny class"<<endl<<endl;
   if(ncurves==0) return;
-  long i;
-  for (i=0; i<ncurves; i++)
-    { Curve ci = (Curve)curves[i];
-      os << (i+1) << ": " << ci;
-      if (i>0) os << "  is "<< isoglist[i]<<"-isogenous to curve "<<fromlist[i]+1;
+  long i=1;
+  for ( const auto& C : curves)
+    {
+      os << i << ": " << (Curve)C;
+      if (i>1)
+        os << "  is "<< isoglist[i-1]<<"-isogenous to curve "<<fromlist[i-1]+1;
       os<<endl;
+      i++;
     }
   os<<endl;
 }
@@ -366,10 +363,14 @@ void IsogenyClass::displaymat(ostream& os)const
   if(ncurves==0) return;
   long i,j;
   os << "Isogeny matrix:\n";
-  os << "\t"; for(j=0; j<ncurves; j++) os<<(j+1)<<"\t";  os<<"\n";
+  os << "\t";
+  for(j=0; j<ncurves; j++)
+    os<<(j+1)<<"\t";  os<<"\n";
   for(i=0; i<ncurves; i++)
     {
-      os<<(i+1)<<"\t"; for(j=0; j<ncurves; j++) os<<mat_entry(i,j)<<"\t";  os<<"\n";
+      os<<(i+1)<<"\t";
+      for(j=0; j<ncurves; j++)
+        os<<mat_entry(i,j)<<"\t";  os<<"\n";
     }
   os<<endl;
 }
@@ -391,8 +392,7 @@ void IsogenyClass::dumpdata(ostream& os, long rank)
 //cout << "After C.get_ntorsion(), C.ntorsion = " << C.ntorsion << endl;
       if(sign(C.discr)>0) os << "+1" <<sep;
       else  os << "-1" <<sep;
-      vector<bigint>::const_iterator pi;
-      pi=C.the_bad_primes.begin();
+      auto pi=C.the_bad_primes.begin();
       while(pi!=C.the_bad_primes.end()) os << C.reduct_array[*pi++].ord_p_discr << sep;
       pi=C.the_bad_primes.begin();
       while(pi!=C.the_bad_primes.end()) os << C.reduct_array[*pi++].ord_p_j_denom << sep;
diff --git a/libsrc/lambda.cc b/libsrc/lambda.cc
index a2ba6a1..2751fbb 100644
--- a/libsrc/lambda.cc
+++ b/libsrc/lambda.cc
@@ -80,10 +80,8 @@ vector<bigfloat> lambda_bad(const CurveRed& C, long& nlambda, int verbose)
   bigint discr = getdiscr(C);
   vector<bigint> plist = getbad_primes(C);
   long i, j, nl, nnl;
-  vector<bigint>::const_iterator pr=plist.begin();
-  while(pr!=plist.end())
+  for (const auto& p : plist)
     {
-      bigint p = *pr++;
       if (ndiv(p*p,discr))
 	{
 	  if(verbose)
@@ -167,10 +165,9 @@ int make_point_from_x(Curvedata* CD, const bigfloat& x, long maxdd, Point* P)
 int make_point_from_x_and_ht(Curvedata* CD, vector<bigfloat> lambdas, const bigfloat& xp, const bigfloat& ht, Point* P)
 {
   bigfloat rh = realheight(xp,CD);
-  vector<bigfloat>::const_iterator lam = lambdas.begin();
-  while(lam!=lambdas.end())
+  for (const auto& lam : lambdas)
     {
-      bigfloat logd = (ht-rh-(*lam++))/2;
+      bigfloat logd = (ht-rh-(lam))/2;
       bigfloat approxd = exp(logd);
       bigint xa, xd2, xd = Iround(approxd);
       if(xd>0) 
diff --git a/libsrc/marith.cc b/libsrc/marith.cc
index 7724391..6fe544f 100644
--- a/libsrc/marith.cc
+++ b/libsrc/marith.cc
@@ -241,17 +241,16 @@ int sqrt_mod_m(bigint& x, const bigint& a, const bigint& m, const vector<bigint>
   bigint mm, p, xp, q; int e;
   x=0;  mm=1;
   
- for(vector<bigint>::const_iterator pr = mpdivs.begin(); pr!=mpdivs.end(); pr++)
+  for( const auto& p : mpdivs)
     {
-      p=*pr;
       e = val(p,m);
       if(e==0) continue;
-      if(p==2) 
+      if(p==2)
 	{if(!sqrt_mod_2_power(xp,a,e)) return 0;}
       else
 	{if(!sqrt_mod_p_power(xp,a,p,e)) return 0;}
       q=pow(p,e);
-      if(pr==mpdivs.begin())
+      if(p==mpdivs.front())
 	x=xp;
       else
 	x=chrem(x,xp,mm,q);
@@ -270,10 +269,9 @@ int modsqrt(const bigint& a, const vector<bigint>& bplist, bigint& x)
 {
   // Assumes b square-free, primes factors in bplist
   bigint u, v, p, amodp, xmodp, m;
-  int res=1;  x=0; m=1;
-  for(vector<bigint>::const_iterator pr = bplist.begin(); res&&(pr!=bplist.end()); pr++)
+  x=0; m=1;
+  for( const auto& p : bplist)
     {
-      p=*pr;
       if(p==2)
 	{
 	  xmodp=odd(a);
@@ -362,11 +360,10 @@ vector<bigint> pdivs_use_factorbase(bigint& n, const std::set<bigint> factor_bas
 {
   vector<bigint> plist;
   if(n<2) return plist;
-  std::set<bigint>::const_iterator pri = factor_base.begin(); 
-  while((n>1)&&(pri!=factor_base.end()))
+  for (const auto& p : factor_base)
     {
-      bigint p=*pri++;
-      if(divide_out(n,p)) 	  
+      if (n==1) break;
+      if(divide_out(n,p))
 	plist.push_back(p);
     }
   return plist;
@@ -466,9 +463,8 @@ factor(const bigint& n, int proof=1)
   oss<<n;
   vector<bigint> plist =  read_vec_from_string(factor(oss.str()));
   if(proof)
-    for(vector<bigint>::const_iterator pi=plist.begin(); pi!=plist.end(); pi++)
+    for( const auto& p : plist)
       {
-	bigint p =*pi;
 	if(!is_prime(p))
 	  {
 	    cout<<"WARNING:  pari's factor() returned p="<<p
@@ -550,10 +546,9 @@ vector<bigint> posdivs(const bigint& number, const vector<bigint>& plist)
  int e, nu = 1; int nd=nu;
  vector<int> elist;
  elist.reserve(np);
- vector<bigint>::const_iterator pr = plist.begin();
- while(pr!=plist.end())
+ for (const auto& p : plist)
    {
-     e=val(*pr++,number);
+     e=val(p,number);
      elist.push_back(e);
      nd*=(1+e);
    }
@@ -562,12 +557,10 @@ vector<bigint> posdivs(const bigint& number, const vector<bigint>& plist)
  // cout<<"In posdivs (1) : dlist = "<<dlist<<endl;
  dlist.resize(nd);
  // cout<<"In posdivs (2) : dlist = "<<dlist<<endl;
- pr=plist.begin();
- vector<int>::iterator ei = elist.begin();
+ auto ei = elist.begin();
  nd=nu;
- while(pr!=plist.end())
+ for (const auto& p : plist)
    {
-     bigint p=*pr++;
      e=*ei++;
      for (int j=0; j<e; j++)
        for (int k=0; k<nd; k++)
@@ -592,10 +585,9 @@ vector<bigint> alldivs(const bigint& number, const vector<bigint>& plist)
  int e, nu = 2; int nd=nu;
  vector<int> elist;
  elist.reserve(np);
- vector<bigint>::const_iterator pr = plist.begin();
- while(pr!=plist.end())
+ for (const auto& p : plist)
    {
-     e=val(*pr++,number);
+     e=val(p,number);
      elist.push_back(e);
      nd*=(1+e);
    }
@@ -603,11 +595,9 @@ vector<bigint> alldivs(const bigint& number, const vector<bigint>& plist)
  dlist.push_back(BIGINT(-1));
  dlist.resize(nd);
  nd=nu;
- pr=plist.begin();
- vector<int>::iterator ei = elist.begin();
- while(pr!=plist.end())
+ auto ei = elist.begin();
+ for (const auto& p : plist)
    {
-     bigint p=*pr++;
      e=*ei++;
      for (int j=0; j<e; j++)
        for (int k=0; k<nd; k++)
@@ -629,21 +619,18 @@ vector<bigint> sqdivs(const bigint& number, const vector<bigint>& plist)
  int e, nu = 1; int nd=nu;
  vector<int> elist;
  elist.reserve(np);
- vector<bigint>::const_iterator pr = plist.begin();
- while(pr!=plist.end())
+ for (const auto& p : plist)
    {
-     e=val(*pr++,number)/2;
+     e=val(p,number)/2;
      elist.push_back(e);
      nd*=(1+e);
    }
  vector<bigint> dlist(1,BIGINT(1)); 
  dlist.resize(nd);
  nd=nu;
- pr=plist.begin();
- vector<int>::iterator ei = elist.begin();
- while(pr!=plist.end())
+ auto ei = elist.begin();
+ for (const auto& p : plist)
    {
-     bigint p=*pr++;
      e=*ei++;
      for (int j=0; j<e; j++)
        for (int k=0; k<nd; k++)
@@ -662,24 +649,14 @@ vector<bigint> sqfreedivs(const bigint& number)
 vector<bigint> sqfreedivs(const bigint& number, const vector<bigint>& plist)
 {
  int np = plist.size();
- int e, nu = 1; int nd=nu;
- vector<int> elist;
- elist.reserve(np);
- vector<bigint>::const_iterator pr = plist.begin();
- while(pr!=plist.end())
-   {
-     e=1; pr++;
-     elist.push_back(e);
-     nd*=(1+e);
-   }
- vector<bigint> dlist(1,BIGINT(1)); 
+ int e = 1, nu = 1, nd=pow(2,np);
+ vector<int> elist(np,1);
+ vector<bigint> dlist(1,BIGINT(1));
  dlist.resize(nd);
  nd=nu;
- pr=plist.begin();
- vector<int>::iterator ei=elist.begin();
- while(pr!=plist.end())
+ auto ei=elist.begin();
+ for (const auto& p : plist)
    {
-     bigint p=*pr++;
      e=*ei++;
      for (int j=0; j<e; j++)
        for (int k=0; k<nd; k++)
@@ -705,10 +682,8 @@ void sqfdecomp(const bigint& a, vector<bigint>& plist, bigint& a1, bigint& a2)
   long j;
   vector<bigint> aplist;
   a1=1;  a2=1;
-  vector<bigint>::const_iterator pr = plist.begin();
-  while(pr!=plist.end())
+ for (const auto& p : plist)
     {
-      bigint p = *pr++;
       long e = val(p,a);
       if(e==0) continue;
       if(e&1) {a1*=p; aplist.push_back(p);}
diff --git a/libsrc/mequiv.cc b/libsrc/mequiv.cc
index 0bc1480..e4f8758 100644
--- a/libsrc/mequiv.cc
+++ b/libsrc/mequiv.cc
@@ -24,8 +24,6 @@
 #include <eclib/compproc.h> // for is_small, is_real
 #include <eclib/mequiv.h>
 
-//#ifdef NEW_EQUIV
-
 int new_equiv(quartic* q1, quartic* q2, int info)
 {
   if(info)
@@ -85,56 +83,6 @@ int new_equiv(quartic* q1, quartic* q2, int info)
 }
 
 
-#if(0) // Old code: compute the complex roots and test for integrality
-
-  // The roots are sqrt(z1*w1)+sqrt(z2*w2)+sqrt(z3*w3) with some choice of signs,
-  // where the wi are the z-values of the second quartic.  We compute the zi
-  // from the roots of the first, then to ensure that conjugates are matched
-  // we compute the wi from the zi (via, invisibly, the phi_i).
-
-  bigcomplex* roots1 = q1->getroots();
-  bigfloat xa1=I2bigfloat(a1), xa2=3*I2bigfloat(a2);
-  bigcomplex rz1=xa1*(roots1[0]+roots1[1]-roots1[2]-roots1[3]);
-  bigcomplex rz2=xa1*(roots1[0]-roots1[1]+roots1[2]-roots1[3]);
-  bigcomplex rz3=xa1*(roots1[0]-roots1[1]-roots1[2]+roots1[3]);
-  bigcomplex t = I2bigfloat(a1*p2-a2*p1), tt=1/I2bigfloat(3*a1);
-  bigcomplex rzw1 = rz1*sqrt((xa2*rz1*rz1+t)*tt);
-  bigcomplex rzw2 = rz2*sqrt((xa2*rz2*rz2+t)*tt);
-  bigcomplex rzw3 = rz3*sqrt((xa2*rz3*rz3+t)*tt);
-
-  bigcomplex u;
-
-  for(long k=0; k<4; k++)
-    {
-      switch(k) {
-      case 0:  u=rzw1+rzw2+rzw3; break;
-      case 1:  u=rzw1+rzw2-rzw3; break;
-      case 2:  u=rzw1-rzw2+rzw3; break;
-      case 3:  u=rzw1-rzw2-rzw3; break;
-      }
-      if(is_approx_zero(imag(u)))
-	{
-	  bigint uu = Iround(real(u));
-	  bigint uu2=uu*uu;
-	  bigint fuu1 = (uu2-2*p)*uu2+s;
-	  bigint fuu2 = 8*r*uu;
-	  if(fuu1==fuu2) 
-	    {
-	      if(info) cout<<"Root u = "<<uu<<endl;
-	      return 1;
-	    }
-	  if(fuu1==-fuu2) 
-	    {
-	      if(info) cout<<"Root u = "<<-uu<<endl;
-	      return 1;
-	    }
-	}
-    }
-
-#endif
-
-//#else // not using new_equiv so no need to compile this stuff
-
 int testd(const bigint& a, const bigint& b, const bigint& c, 
 	  const bigint& d, const bigint& e, const bigint& as, 
 	  const bigint& bs, const bigint& cs, const bigint& ds, 
@@ -193,7 +141,7 @@ bigcomplex crossratio(const bigcomplex& x1,const bigcomplex& x2,const bigcomplex
 }
 
 int rootsequiv(const quartic* q1, const quartic* q2, int i, const vector<bigint>& dlist, int info)
-{  int ans=0;
+{
    bigcomplex *x = q1->getroots();
    bigcomplex *y = q2->getroots();
    bigcomplex x1=x[0], x2=x[1], x3=x[2], x4=x[3];
@@ -247,18 +195,15 @@ int rootsequiv(const quartic* q1, const quartic* q2, int i, const vector<bigint>
           delta = real(cdelta);
    bigfloat det = alpha*delta-beta*gamma;
 
-   if (info) 
+   if (info)
      {
        cout << "Real transformation has alpha, beta, gamam, delta = ";
        cout <<alpha<<" "<<beta<<" "<<gamma<<" "<<delta<<endl;
        cout << "Testing divisors of "<<q1->getdisc()<<":\n";
 //       cout << dlist << endl;
      }
-   bigint d;
-   vector<bigint>::const_iterator dvar;
-   for (dvar=dlist.begin(); dvar!=dlist.end() && (!ans); dvar++)
+   for ( const auto& d : dlist)
      {
-       d = *dvar;
        bigfloat rscale = sqrt(abs(I2bigfloat(d)/det));
        bigfloat rscaler = floor(rscale+0.5);
        if(abs(rscale-rscaler)<0.001)
@@ -275,16 +220,17 @@ int rootsequiv(const quartic* q1, const quartic* q2, int i, const vector<bigint>
 		  cout<<"rscale = "<<rscale<<endl;
 		  cout << "al,be,ga,de = "<<al<<" "<<be<<" "<<ga<<" "<<de<<endl;
 		}
-	       ans=testd(q1->geta(),q1->getb(),q1->getcc(),q1->getd(),q1->gete(),
+	       if (testd(q1->geta(),q1->getb(),q1->getcc(),q1->getd(),q1->gete(),
 			 q2->geta(),q2->getb(),q2->getcc(),q2->getd(),q2->gete(),
-			 d,al,be,ga,de,info);
+			 d,al,be,ga,de,info))
+                 return 1;
 	     }
 	 }
      }
-   return ans ;
+   return 0;
  }    // of rootsequiv()
- 
- 
+
+
 int equiv(const quartic* q1, const quartic* q2, const vector<bigint>& dlist, int info)
 {
    bigint iiq1 = q1->getI(), jjq1 = q1->getJ(), discq1 = q1->getdisc();
@@ -307,8 +253,8 @@ int equiv(const quartic* q1, const quartic* q2, const vector<bigint>& dlist, int
        if(info) {if(!ans) cout << "Not "; cout<<"equiv\n";}
        return ans;
      }
-   else 
-     { 
+   else
+     {
        if (info) {cout << "equiv failed on first test!\n";
                   cout << "First  has I="<<iiq1<<", J="<<jjq1<<",";
                   cout << " disc="<<discq1<<", type="<<typeq1<<endl;
@@ -318,6 +264,3 @@ int equiv(const quartic* q1, const quartic* q2, const vector<bigint>& dlist, int
        return 0;
      }
  }  // of equiv()
-
-
-//#endif
diff --git a/libsrc/mlocsol.cc b/libsrc/mlocsol.cc
index 3df2f22..e7046aa 100644
--- a/libsrc/mlocsol.cc
+++ b/libsrc/mlocsol.cc
@@ -125,7 +125,7 @@ int zpsol(const bigint& a,const bigint& b,const bigint& c,const bigint& d,const
        {
 	 if(p==2) if(lemma7(a,b,c,d,e,p,nu+1, x)==+1) return 1;
 	 if(lemma6(a,b,c,d,e,p,nu+1, x)==+1) return 1;
-       }    
+       }
    }
 
  for(i=0; i<p; ++i, x+=pnu)
@@ -212,8 +212,8 @@ int locallysoluble(const bigint& a, const bigint& b, const bigint& c, const bigi
       bigint D = c*c-4*a*e;
       if(global_hilbert(a,D,plist,badp)) return 0;
     }
-  for (vector<bigint>::const_iterator p = plist.begin(); (p!=plist.end())&&sol; p++)
-   { badp = *p;
+  for ( const auto& badp : plist)
+   {
      sol = new_qpsoluble(a,b,c,d,e,badp,0);
 #ifdef CHECK_LOC_SOL
      cout<<"Checking solubility with old method (p="<<p<<")...\t";
@@ -226,6 +226,7 @@ int locallysoluble(const bigint& a, const bigint& b, const bigint& c, const bigi
        }
      else cout << "...OK.\n";
 #endif
+     if (!sol) break;
    }
    return sol;
 }  /* end of locallysoluble */
diff --git a/libsrc/newforms.cc b/libsrc/newforms.cc
index 40afb3c..7a1a00a 100644
--- a/libsrc/newforms.cc
+++ b/libsrc/newforms.cc
@@ -49,10 +49,10 @@ int less_ap(long a, long b, int old=0)
 int less_apvec(const vector<long>& v, const vector<long>& w, int old=0); 
 int less_apvec(const vector<long>& v, const vector<long>& w, int old) 
 {
-  vector<long>::const_iterator vi=v.begin(), wi=w.begin();
-  while(vi!=v.end())
+  auto wi=w.begin();
+  for ( const auto& vi : v)
     {
-      int s = less_ap(*vi++,*wi++,old); 
+      int s = less_ap(vi,*wi++,old);
       if(s) return s;
     }
   return 0;
@@ -94,25 +94,24 @@ vector<long> eiglist(const newform& f, int oldorder)
   vector<long> eigs;
   primevar pr;
   long N = (f.nf)->modulus;
-  vector<long>::const_iterator aqi=f.aqlist.begin();
-  vector<long>::const_iterator api=f.aplist.begin();
-  vector<long>::iterator eigsi;
+  auto aqi=f.aqlist.begin();
+  auto api=f.aplist.begin();
   if(oldorder)
     {
       eigs.resize(f.aplist.size());
-      eigsi=eigs.begin();
-      while(aqi!=f.aqlist.end()) 
+      auto eigsi=eigs.begin();
+      while(aqi!=f.aqlist.end())
 	*eigsi++ = *aqi++;
       while(api!=f.aplist.end())
 	{
-	  if(ndivides(pr,N)) *eigsi++ = *api;     
+	  if(ndivides(pr,N)) *eigsi++ = *api;
 	  api++; pr++;
 	}
     }
   else
     {
       eigs=f.aplist; // copy; now adjust the aq:
-      eigsi=eigs.begin();
+      auto eigsi=eigs.begin();
       while((aqi!=f.aqlist.end())&&(eigsi!=eigs.end()))
 	{
 	  if(::divides(pr.value(),N)) *eigsi = (*aqi++);
@@ -254,16 +253,17 @@ void newform::fixup_eigs()
 { 
   long denom = nf->h1->h1denom();
   aqlist.resize(nf->npdivs);
-  vector<long>::iterator api=aplist.begin(), pi=nf->plist.begin();
-  vector<long>::iterator aqi=aqlist.begin();
+  auto api=aplist.begin();
+  auto pi=nf->plist.begin();
+  auto aqi=aqlist.begin();
   primevar pr;   long q, i;
   long n = nf->modulus;
   while((api!=aplist.end())&&(aqi!=aqlist.end()))
     {
       q=pr.value(); pr++;
-      if(::divides(q,n)) 
+      if(::divides(q,n))
 	{
-	  *aqi++=*api; 
+	  *aqi++=*api;
 	  *api=(::divides(q*q,n)? 0: -*api);
 	  pi++;
 	}
@@ -274,9 +274,9 @@ void newform::fixup_eigs()
       long piv;
       ssubspace espace;
       if(sign==-1)
-        espace=make1d(bminus,piv); 
+        espace=make1d(bminus,piv);
       else
-        espace=make1d(bplus,piv); 
+        espace=make1d(bplus,piv);
       piv*=denom;
       while(aqi!=aqlist.end()) // compute missing aq
 	{
@@ -289,7 +289,7 @@ void newform::fixup_eigs()
 	}
     }
   if(nf->verbose) cout<<"aqlist = "<<aqlist<<endl;
-  
+
   //Compute sfe:
 
   sfe=-1;
@@ -300,14 +300,14 @@ void newform::fixup_eigs()
 // Before recovering eigenbases, we need to put back the aq into the
 // aplist (and resort, for efficiency).
 void newform::unfix_eigs()
-{ 
-  vector<long>::iterator api=aplist.begin();
-  vector<long>::iterator aqi=aqlist.begin();
+{
+  auto api=aplist.begin();
+  auto aqi=aqlist.begin();
   primevar pr;
   long n = nf->modulus;
   while((api!=aplist.end())&&(aqi!=aqlist.end()))
     {
-      if(::divides(pr.value(),n)) *api=*aqi++; 
+      if(::divides(pr.value(),n)) *api=*aqi++;
       api++;
       pr++;
     }
@@ -315,14 +315,14 @@ void newform::unfix_eigs()
 
 // After recovering eigenbases, we need to replace the ap for bad p
 void newform::refix_eigs()
-{ 
-  vector<long>::iterator api=aplist.begin();
+{
+  auto api=aplist.begin();
   primevar pr;
   long n = nf->modulus, np = nf->npdivs, ip=0, q;
   while((api!=aplist.end())&&(ip<np))
     {
       q=pr.value();
-      if(::divides(q,n)) 
+      if(::divides(q,n))
 	{
 	  *api=(::divides(q*q,n)? 0: -*api);
 	  ip++;
@@ -336,8 +336,8 @@ void newform::find_bsd_ratio()
 {
   // get ap for p=p0:
 
-  primevar pr; 
-  vector<long>::const_iterator api=aplist.begin();
+  primevar pr;
+  auto api = aplist.begin();
   while(pr.value()!=nf->p0) {pr++; api++;}
   ap0=*api;
   np0 = 1 + (nf->p0) - ap0;
@@ -364,13 +364,13 @@ void newform::find_bsd_ratio()
       if(denomplus>1)
 	{
 	  if(::divides(denomplus,dp0))  dp0/=denomplus;
-	  else 
+	  else
 	    cout<<"newform constructor error: dp0 not divisible by denomplus!"
-		<<endl; 
+		<<endl;
 	}
     }
-  loverp = rational(dp0,np0);   
-  if(nf->verbose) 
+  loverp = rational(dp0,np0);
+  if(nf->verbose)
     {
       cout<<"pdot = "<<pdot<<endl;
       cout<<"dp0 = "<<dp0<<endl;
@@ -494,15 +494,15 @@ void newform::find_twisting_primes()
           mplus =0;
         }
     }
-  if(sign!=+1) 
+  if(sign!=+1)
     {
       lminus=0;
       mminus=0;
-    }  
+    }
   if(nf->squarelevel) return;
 
   long n = nf->modulus;
-  for (primevar lvar; lvar.ok() && 
+  for (primevar lvar; lvar.ok() &&
 	 (((sign!=-1)&&(mplus==0)) ||
 	  ((sign!=+1)&&(mminus==0))); lvar++)
     {
@@ -515,9 +515,9 @@ void newform::find_twisting_primes()
 
       if((sign!=-1)&&(mplus==0)&&(l%4==1))
 	{
-	  lplus = l;   
+	  lplus = l;
 	  //cout << "Trying lplus = " << l << "\n";
-	  map<long,vec>::const_iterator vi = nf->mvlplusvecs.find(l);
+	  auto  vi = nf->mvlplusvecs.find(l);
 	  if(vi==nf->mvlplusvecs.end())
 	    mplus = (nf->mvlplusvecs[l]=nf->h1->manintwist(l))*bplus;
 	  else
@@ -535,16 +535,16 @@ void newform::find_twisting_primes()
 	  if((denomplus>1)&&(mplus!=0))
 	    {
 	      if(::divides(denomplus,mplus))  mplus/=denomplus;
-	      else 
+	      else
 		cout<<"Warning in newform constructor: mplus not divisible by denomplus!"
-		    <<endl; 
+		    <<endl;
 	    }
 	}
       if((sign!=+1)&&(mminus==0)&&(l%4==3))
 	{
-	  lminus = l;   
+	  lminus = l;
 	  //cout << "Trying lminus = " << l << "\n";
-	  map<long,vec>::const_iterator vi = nf->mvlminusvecs.find(l);
+	  auto vi = nf->mvlminusvecs.find(l);
 	  if(vi==nf->mvlminusvecs.end())
 	    mminus = (nf->mvlminusvecs[l]=nf->h1->manintwist(l))*bminus;
 	  else
@@ -1550,7 +1550,7 @@ vector<long> aqlist(vector<long> aplist, long N)
 {
   long iq=0, p, naq = pdivs(N).size();
   //cout << "Setting aq of size "<<naq<<endl;
-  vector<long>::const_iterator api = aplist.begin();
+  auto api = aplist.begin();
   vector<long> aq(naq);
   for(primevar pr; (iq<naq)&&pr.ok(); pr++)
     {
@@ -1911,10 +1911,9 @@ vector<long> newforms::apvec(long p) //  computes a[p] for each newform
   //cout<<"Computing images of M-symbols"<<endl<<flush;
   //cout<<"jlist = "<<jlist<<endl;
 
-  for(std::set<long>::const_iterator jj=jlist.begin(); jj!=jlist.end(); jj++)
+  for( const auto& j : jlist)
     {
       imagej=vec(n1ds); // initialised to 0
-      j=*jj;
       symb s = h1->symbol(h1->freegens[j-1]);
       //cout<<"Computing image of "<<j<<"'th M-symbol "<<s<<endl;
       //cout<<" = "<<s<<"..."<<flush;
@@ -2075,18 +2074,17 @@ vector<int> newforms::showcurves(vector<int> forms, int verbose, string filename
   bigfloat rperiod;
   bigint a1,a2,a3,a4,a6, N;
   vector<int> badcurves; // will hold the indices of forms for which we fail to find a curve
-  vector<int>::const_iterator inf; // will iterate through the forms to be used
 
-  for(inf=forms.begin(); inf!=forms.end(); inf++)
+  for( const auto& inf : forms)
    {
      if(verbose)
-       cout<<"\n"<<"Form number "<<*inf+1<<"\n";
-     else cout<<(*inf+1)<<" ";
+       cout<<"\n"<<"Form number "<<inf+1<<"\n";
+     else cout<<(inf+1)<<" ";
 
      if (output_curves)
-       curve_out << modulus << " "<< codeletter(*inf) << " 1 ";
+       curve_out << modulus << " "<< codeletter(inf) << " 1 ";
 
-     Curve C = getcurve(*inf,-1,rperiod,verbose);
+     Curve C = getcurve(inf,-1,rperiod,verbose);
      Curvedata CD(C,1);  // The 1 causes minimalization
      if(verbose) cout << "\nCurve = \t";
      cout << (Curve)CD << "\t";
@@ -2098,7 +2096,7 @@ vector<int> newforms::showcurves(vector<int> forms, int verbose, string filename
      if(N!=modulus)
        {
 	 cout<<"No curve found"<<endl;
-	 badcurves.push_back(*inf);
+	 badcurves.push_back(inf);
          if (output_curves)
            curve_out<<endl;
        }
@@ -2108,7 +2106,7 @@ vector<int> newforms::showcurves(vector<int> forms, int verbose, string filename
            C.getai(a1,a2,a3,a4,a6);
            curve_out<<"["<<a1<<","<<a2<<","<<a3<<","<<a4<<","<<a6<<"]";
            int nt = CD.get_ntorsion();
-           int r = nflist[*inf].rank(); // analytic rank
+           int r = nflist[inf].rank(); // analytic rank
            curve_out<<" "<<r<<" "<<nt<<" 0"<<endl;
          }
    }
diff --git a/libsrc/oldforms.cc b/libsrc/oldforms.cc
index cd25a39..1d13f1e 100644
--- a/libsrc/oldforms.cc
+++ b/libsrc/oldforms.cc
@@ -129,10 +129,8 @@ oldforms::oldforms(long intp, const level* iN, int verbose, int plus)
    nap = ntp;
    noldclasses = totalolddim = 0;
    plusflag=plus;
-   vector<long>::const_iterator d;
-   for(d=(N->dlist).begin();d!=(N->dlist).end();d++)
+   for (const auto& M : N->dlist)
      {
-       long M=*d;
        if(M<11) continue;
        if(M==(N->modulus)) continue;
        getoldclasses(M,verbose);
@@ -170,36 +168,32 @@ void oldforms::getoldclasses(long d, int verbose)
       long m = n/d;
       long k=0, xmult, mult, j, beta;
       vector<long> betalist; // =new long[N->npdivs];
-      vector<long>::const_iterator qj=(N->plist).begin();
-      while(qj!=(N->plist).end())
+      for (const auto& qj : N->plist)
 	{
-	  beta=val(*qj++,m);
+	  beta=val(qj,m);
 	  if(beta>0) k++;
 	  betalist.push_back(beta);
 	}
       if(verbose>1) cout<<"betas: "<<betalist<<endl;
       vector<long> nextoldformap(nap);
-      vector<long>::iterator betai;
-      vector<long>::const_iterator aqj;
       primevar pr; long iform, c, ip, p, aq; int bit;
       for(iform=0; iform<nforms; iform++)
-	{ 
+	{
 	  vector<long>& aqlist=olddata.nflist[iform].aqlist;
 	  nextoldformap = olddata.nflist[iform].aplist;
-	  if(verbose>1) 
+	  if(verbose>1)
 	    {
 	      cout<<"form #"<<(iform+1)<<": "<<"aqlist="<<aqlist<<endl;
 	      cout<<"aplist before adjusting (size "<<oldnap<<") ="<<nextoldformap<<endl;
 	    }
 	  for (c=0; c<(1<<k); c++) // 2^k different oldclasses
-	    {  
+	    {
 	      if(verbose>1) cout<<"c="<<c<<endl;
-	      mult=1; j=0; 
-	      betai=betalist.begin();
-	      aqj=aqlist.begin();
-              for (qj=(N->plist).begin(); qj!=(N->plist).end(); qj++)
+	      mult=1; j=0;
+	      auto betai=betalist.begin();
+	      auto aqj=aqlist.begin();
+              for ( const auto& p : N->plist)
 		{
-		  p = *qj;
                   ip = prime_pi(p);
 		  if(verbose>1) cout<<"p="<<p<<" (ip="<<ip<<")"<<endl;
                   beta = *betai++;
diff --git a/libsrc/periods.cc b/libsrc/periods.cc
index d5819ad..0ba6ee5 100644
--- a/libsrc/periods.cc
+++ b/libsrc/periods.cc
@@ -80,7 +80,7 @@ vector<long> resort_aplist(const level* iN,
   long i, j, p, ip = iN->npdivs;
   for(i=0; i<nap; i++)
     { p = primelist[i];
-    vector<long>::const_iterator pi = find(iN->plist.begin(),iN->plist.end(),p);
+    auto pi = find(iN->plist.begin(),iN->plist.end(),p);
     if(pi==iN->plist.end()) // then p is good
       {
 	aplist.push_back( apl[ip++]);
@@ -180,9 +180,9 @@ void summer::use2357(long n, long an)
 		{
 		  am2357=am235*a7p_cache[i7];
 		  use(m2357,am2357);
-		}      	      
-	    }      
-	}      
+		}
+	    }
+	}
     }
 }
 
diff --git a/libsrc/points.cc b/libsrc/points.cc
index 5398cc0..967ff46 100644
--- a/libsrc/points.cc
+++ b/libsrc/points.cc
@@ -246,27 +246,13 @@ int Point::isvalid() const // P on its curve ?
   if((sign(X)==0)&&(sign(Y)==0)&&(sign(Z)==0)) return 0;   
 // Point at infinity is on a curve
   if((sign(X)==0)&&(sign(Z)==0)) return 1 ;
-  else{
-    // should calculate 
-    //                       y^2 +a1 x y + a3 y 
-    //  and
-    //                       x^3 + a2 x^2 + a4 x + a6
-    // where 
-    //          x = X/Z, y = Y/Z
-    // and verify equality.
-    //
-    // In homogeneous coordinates:
-    //
-    // Lhs = Y^2 Z + a1 XYZ + a3 YZ^2 = (YZ) *(Y + a1 X + a3 Z)
-    //
-    //
-    // Rhs = X^3 +a2 X^2 Z + a4 X Z^2 + a6 Z^3
-    //
-    bigint A1,A2,A3,A4,A6;    E->getai(A1,A2,A3,A4,A6);
-    const bigint& Lhs = Y*Z*(Y + A1*X + A3*Z) ;
-    const bigint& Rhs = A6*pow(Z,3) + X*(A4*Z*Z + X*(A2*Z + X)) ;
-    return Lhs == Rhs ;
-  }
+  else
+    {
+      bigint A1,A2,A3,A4,A6;    E->getai(A1,A2,A3,A4,A6);
+      const bigint& Lhs = Y*Z*(Y + A1*X + A3*Z) ;
+      const bigint& Rhs = A6*pow(Z,3) + X*(A4*Z*Z + X*(A2*Z + X)) ;
+      return Lhs == Rhs ;
+    }
 }
 
 // Find all points with a given rational x-coordinate:
@@ -387,17 +373,11 @@ vector<Point> two_torsion(Curvedata& E, int exact)
       b2=a2; b4=a4; b6=a6;
     }
   vector<bigint> xlist = Introotscubic(b2,b4,b6);
-  int n, n2tors = xlist.size();
-
-  // If there are 3 points of order 2, we order them for consistency:
-  if(n2tors==3)  sort(xlist.begin(),xlist.end());
-
   vector<Point> two_tors;
   if (!exact)
     two_tors.push_back(Point(E)) ;     // zero point
-  for(n=0; n<n2tors; n++)
+  for( const auto& ei : xlist)
     {
-      bigint ei = xlist[n];
       if(scaled)
 	two_tors.push_back(Point(E,2*ei,-a1*ei-4*a3,BIGINT(8)));
       else
@@ -452,9 +432,8 @@ vector<Point> three_torsion(Curvedata& E, int exact)
   vector<Point> three_tors;
   if (!exact)
     three_tors.push_back(Point(E)) ;     // zero point
-  for(unsigned int n=0; n<xlist.size(); n++)
+  for( auto& xi : xlist)
     {
-      xi = xlist[n];
       if(xi%3==0) // 3-torsion must be integral
 	{
 	  xi/=3;
@@ -499,15 +478,14 @@ vector<Point> m_torsion(Curvedata& E, long m, int exact)
 #endif
 
   // accumulate the points with each x-coorindate:
-  for(auto xi = xs.begin(); xi!=xs.end(); ++xi)
+  for( const auto& xi : xs)
     {
-      vector<Point> Ps = points_from_x(E, *xi);
+      vector<Point> Ps = points_from_x(E, xi);
 #ifdef DEBUG_TORSION
-      cout<<" x = "<<(*xi)<< " gives points "<<Ps<<" of order dividing "<<m<<": "<<Ps<<"\n";
+      cout<<" x = "<<xi<< " gives points "<<Ps<<" of order dividing "<<m<<": "<<Ps<<"\n";
 #endif
-      for (auto Pi=Ps.begin(); Pi!=Ps.end(); ++Pi)
+      for (auto P : Ps)
         {
-          Point P = *Pi;
           if ((!exact)||(order(P)==m))
             m_tors.push_back(P);
         }
@@ -592,18 +570,13 @@ vector<Point> torsion_points(Curvedata& E)
   // Now we have structure C2xCn for n = 4, 6, or 8
 
   // Find a point of order 2 not in the second factor and add it in
-  Point T1 = cycle[exponent/2], T2 = two_tors[0];
-  for( auto Ti=two_tors.begin()+1; Ti!=two_tors.end(); ++Ti)
-    {
-      T2 = *Ti;
-      if (T2 != T1)
-        break;
-    }
+  Point T1 = cycle[exponent/2];
+  auto Ti = std::find_if_not(two_tors.begin()+1, two_tors.end(), [T1]( const auto& T ){return T==T1;});
+  Point T2 = *Ti;
   vector<Point> all_torsion = cycle;
-  for( auto Pi=cycle.begin(); Pi!=cycle.end(); ++Pi)
-    {
-      all_torsion.push_back(T2+*Pi);
-    }
+  for( auto P : cycle)
+    all_torsion.push_back(T2+P);
+
   ::sort(all_torsion.begin(), all_torsion.end(), Point_cmp);
  return all_torsion;
 }
@@ -616,16 +589,14 @@ vector<Point> Point::division_points(int m) // list of Q s.t. n*Q=this
 #endif
   vector<Point> ans;
   vector<Point> Qs;
-  Point Q;
   if (is_torsion())
     {
       Qs = torsion_points(*E);
 #ifdef DEBUG_DIVISION_POINTS
       cout << "Testing all torsion points " << Qs << endl;
 #endif
-      for (vector<Point>::iterator Qi = Qs.begin(); Qi!=Qs.end(); Qi++)
+      for ( const auto& Q : Qs)
         {
-          Q = *Qi;
 #ifdef DEBUG_DIVISION_POINTS
           cout << "   Q = "<<Q<<", "<<m<<"*Q="<<m*Q<<endl;
 #endif
@@ -642,23 +613,25 @@ vector<Point> Point::division_points(int m) // list of Q s.t. n*Q=this
 #endif
       return ans;
     }
+
+  // Now in the non-torsion case
+
   ZPoly pol = division_points_X_pol(E->a1,E->a2,E->a3,E->a4,E->a6, m, X, Z);
 #ifdef DEBUG_DIVISION_POINTS
   cout << "division poly = " << pol << endl;
 #endif
-
   vector<bigrational> xQs = roots(pol);
 #ifdef DEBUG_DIVISION_POINTS
   cout << " with roots " << xQs << endl;
 #endif
-  for(vector<bigrational>::iterator xQi = xQs.begin(); xQi!=xQs.end(); xQi++)
+  for( const auto& xQ : xQs)
     {
-      Qs = points_from_x(*E, *xQi);
+      Qs = points_from_x(*E, xQ);
       // will have length 0 or 2 since non-torsion, and we only want
       // exctly one when there two so, must check which works
       if (Qs.size()>0)
         {
-          Q = Qs[0];
+          Point Q = Qs[0];
           if (m*Q == *this)
             {
               ans.push_back(Q);
@@ -738,11 +711,11 @@ int Point::has_good_reduction(const vector<bigint>& plist, bigint& p0, int check
         p0 = BIGINT(0);
         return 0;
       }
-  for(vector<bigint>::const_iterator pi = plist.begin(); pi!=plist.end(); ++pi)
+  for( const auto& p : plist)
     {
-      if(!has_good_reduction(*pi))
+      if(!has_good_reduction(p))
         {
-          p0 = *pi;
+          p0 = p;
           return 0;
         }
     }
diff --git a/libsrc/pointsmod.cc b/libsrc/pointsmod.cc
index 6ca5d51..a46de34 100644
--- a/libsrc/pointsmod.cc
+++ b/libsrc/pointsmod.cc
@@ -558,7 +558,7 @@ bigint my_bg_algorithm(const pointmodq& PP,
       
       // look in table to see if H2= i*P for a suitable i
 
-      map<bigint,long>::iterator HTi = HT.find(LiftGF(H2.get_x()));
+      auto HTi = HT.find(LiftGF(H2.get_x()));
       if(HTi!=HT.end())
 	{
 	  i = HTi->second;
@@ -580,12 +580,11 @@ bigint my_bg_algorithm(const pointmodq& PP,
 bigint my_order_point(const pointmodq& P, const bigint& mult)
 {
   vector<bigint> plist = pdivs(mult);
-  unsigned int i; bigint m, p, ans = BIGINT(1);
+  bigint ans = BIGINT(1);
   if(P.is_zero()) return ans;
-  for(i=0; i<plist.size(); i++)
+  for( const auto& p : plist)
     {
-      p = plist[i];
-      m = mult;  
+      bigint m = mult;
       divide_out(m,p);
       pointmodq Q = m*P;
       while(!Q.is_zero()) {Q=p*Q; ans*=p;}
diff --git a/libsrc/polys.cc b/libsrc/polys.cc
index 3d10176..e5f5bd9 100644
--- a/libsrc/polys.cc
+++ b/libsrc/polys.cc
@@ -111,12 +111,9 @@ vector<bigint> introots(const ZPoly& f)
   vector<bigint> ans;
   if (ratroots.size()==0)
     return ans;
-  for(vector<bigrational>::iterator ri=ratroots.begin(); ri!=ratroots.end(); ri++)
-    {
-      bigrational r = *ri;
-      if (den(r)==1)
-        ans.push_back(num(r));
-    }
+  for( const auto& r : ratroots)
+    if (den(r)==1)
+      ans.push_back(num(r));
   sort(ans.begin(), ans.end());
   return ans;
 }
diff --git a/libsrc/saturate.cc b/libsrc/saturate.cc
index 7e60a01..a83ca5a 100644
--- a/libsrc/saturate.cc
+++ b/libsrc/saturate.cc
@@ -136,27 +136,13 @@ void saturator::nextq()
 
       // First just check the order of E mod q, skip this q if not a multiple of p
 
-      map<bigint,bigint>::iterator Eqoi = Emodq_order.find(q);
       bigint order_mod_q;
+      auto Eqoi = Emodq_order.find(q);
       if(Eqoi==Emodq_order.end())
 	{
 	  if(verbose>2) cout<<"Computing order mod q =  "<<q<<": "<<endl;
           curvemodq Eq(*E,q);
-          if(0) // use orders of some random points as a proxy
-            {
-              bigint upper, lower; // bounds on group order
-              set_hasse_bounds(q,lower,upper);
-              pointmodq P1 = Eq.random_point();
-              bigint n1 = my_order_point(P1,lower,upper);
-              if (verbose>2)
-                cout<<"q="<<q<<"\tn1 = "<<n1<<endl;
-              order_mod_q = n1;
-            }
-          else
-            {
-              order_mod_q = Eq.group_order();
-            }
-
+          order_mod_q = Eq.group_order();
           Emodq_order[q] = order_mod_q;
           if (verbose>2)
             cout<<"Setting order mod "<<q<<" to "<<order_mod_q<<endl;
@@ -180,7 +166,7 @@ void saturator::nextq()
         }
       // next compute the structure of q if not yet known
 
-      map<bigint,curvemodqbasis>::iterator Eqi = Emodq.find(q);
+      auto Eqi = Emodq.find(q);
       if(Eqi==Emodq.end())
 	{
 	  if(verbose>2) cout<<"Initializing q =  "<<q<<": "<<endl;
@@ -491,9 +477,8 @@ int saturator::saturate(vector<long>& unsat, long& index,
     {
       if (verbose)
         cout << "Tamagawa index primes are " << tam_primes << endl;
-      for (vector<long>::iterator pi = tam_primes.begin(); pi!=tam_primes.end(); pi++)
+      for ( const auto& p : tam_primes)
         {
-          p = *pi;
           if ((p > ib) && ((sat_bd==-1) || (p <= sat_bd)))
             {
               if (verbose)
@@ -513,22 +498,14 @@ int saturator::saturate(vector<long>& unsat, long& index,
 void saturator::show_q_tally()
 {
   cout << "Summary of auxiliary primes used" <<endl;
-  int num_q_used = 0;
-  map<bigint,int>::iterator qcount;
   cout << "Number of q used: " << q_tally.size() << endl;
   cout << "Maximum   q used: " << maxq << " (used for p="<<maxp<<")"<<endl;
   if (verbose<2)
     return;
   cout << "Counts of how many times each q was used:" << endl;
-  bigint q;
-  int c;
-  for (qcount = q_tally.begin(); qcount!=q_tally.end(); qcount++)
-    {
-      q = qcount->first;
-      c = qcount->second;
-      if (c)
-        cout << q << "\t" << c <<endl;
-    }
+  for ( const auto& qcount : q_tally)
+    if (qcount.second)
+      cout << qcount.first << "\t" << qcount.second <<endl;
 }
 
 
diff --git a/libsrc/sifter.cc b/libsrc/sifter.cc
index 1ce2dee..2d7ed65 100644
--- a/libsrc/sifter.cc
+++ b/libsrc/sifter.cc
@@ -174,11 +174,10 @@ int* sifter::eps(const bigint& x, const bigint& z2)
 
 void sifter::process(const vector<Point>& Plist) 
 {
-  vector<Point>::const_iterator P=Plist.begin();
-  while(P!=Plist.end()) 
+  for ( const auto& P : Plist)
     {
-      if(verbose) cout<<"Processing point "<<*P<<endl;
-      process(*P++);
+      if(verbose) cout<<"Processing point "<<P<<endl;
+      process(P);
     }
 }
 
diff --git a/libsrc/smat.cc b/libsrc/smat.cc
index 57e1d6d..cbac61f 100644
--- a/libsrc/smat.cc
+++ b/libsrc/smat.cc
@@ -200,11 +200,10 @@ void smat::setrow ( int i, const svec& v) // i counts from 1
   }
 
   coli++;
-  for(map<int,scalar>::const_iterator vi=v.entries.begin();
-      vi!=v.entries.end(); vi++)
+  for( const auto& vi : v.entries)
     {
-      *coli++ = vi->first;
-      *vali++ = vi->second;
+      *coli++ = vi.first;
+      *vali++ = vi.second;
     }
 
   // check
@@ -580,56 +579,34 @@ vec operator*  (smat& m, const vec& v)
 svec operator* ( const svec& v, const smat& A )
 {
   svec prod(A.ncols());
-  if( v.d != A.nrows() ) 
-    { 
-      cerr << "incompatible sizes in v*A\n"; 
-      cerr << "Dimensions "<<v.d<<" and "<<dim(A)<<endl;
-    }
-  else
+  if( v.d != A.nrows() )
     {
-      map<int,scalar>::const_iterator vi;
-      for(vi=v.entries.begin(); vi!=v.entries.end(); vi++)
-        prod += (vi->second)*(A.row(vi->first));
-    }
-  return prod;
-}
-
-#if(0)
-svec mult_mod_p( const svec& v, const smat& A, const scalar& p  )
-{
-  svec prod(A.ncols());
-  if( v.d != A.nrows() ) 
-    { 
-      cerr << "incompatible sizes in v*A\n"; 
-      cerr << "Dimensions "<<v.d<<" and "<<dim(A)<<endl;
+      cerr << "incompatible sizes in v*A\n"
+           << "Dimensions "<<v.d<<" and "<<dim(A)<<endl;
     }
   else
     {
-      map<int,scalar>::const_iterator vi;
-      for(vi=v.entries.begin(); vi!=v.entries.end(); vi++)
-        prod.add_scalar_times_mod_p(A.row(vi->first), vi->second,p);
+      for( const auto& vi : v.entries)
+        prod += (vi.second)*(A.row(vi.first));
     }
   return prod;
 }
 
-#else
-
 svec mult_mod_p( const svec& v, const smat& A, const scalar& p  )
 {
   vec prod(A.ncols());
   if( v.d != A.nrows() )
     {
-      cerr << "incompatible sizes in v*A\n";
-      cerr << "Dimensions "<<v.d<<" and "<<dim(A)<<endl;
+      cerr << "incompatible sizes in v*A\n"
+           << "Dimensions "<<v.d<<" and "<<dim(A)<<endl;
     }
   else
     {
-      map<int,scalar>::const_iterator vi;
-      for(vi=v.entries.begin(); vi!=v.entries.end(); vi++)
+      for( const auto& vi : v.entries)
         {
           // prod.add_scalar_times_mod_p(A.row(vi->first), vi->second,p);
-          int i = (vi->first)-1;        // the row of A to use (from 0)
-          scalar c = vi->second;        // the coefficient tomultiply it by
+          int i = (vi.first)-1;         // the row of A to use (from 0)
+          scalar c = vi.second;         // the coefficient to multiply it by
           int d = A.col[i][0];          // #nonzero entries in this row
           int *posi = A.col[i] +1;      // pointer to array of columns
           scalar *values = A.val[i];    // pointer to array of values
@@ -639,7 +616,7 @@ svec mult_mod_p( const svec& v, const smat& A, const scalar& p  )
     }
   return svec(prod);
 }
-#endif
+
 
 svec mult_mod_p( const smat& A, const svec& v, const scalar& p  )
 {
diff --git a/libsrc/smat_elim.cc b/libsrc/smat_elim.cc
index df65b7b..b7a3575 100644
--- a/libsrc/smat_elim.cc
+++ b/libsrc/smat_elim.cc
@@ -1273,17 +1273,14 @@ void smat_elim::step5dense()
       return;
     }
   mat dmat(nrr, nrc);
-  map<int,scalar>::const_iterator vi;
-  vector<int>::const_iterator rci;
   for (i=0; i<nrr; i++)
     {
       svec v = row(remaining_rows[i]);
       j=0;
-      for(vi=v.entries.begin();
-          vi!=v.entries.end(); vi++)
+      for( const auto& vi : v.entries)
         {
-          while (remaining_cols[j]<(vi->first)) j++;
-          dmat.set(i+1,j+1,vi->second);
+          while (remaining_cols[j]<(vi.first)) j++;
+          dmat.set(i+1,j+1,vi.second);
         }
     }
 
diff --git a/libsrc/svec.cc b/libsrc/svec.cc
index a38029b..e68c3c1 100644
--- a/libsrc/svec.cc
+++ b/libsrc/svec.cc
@@ -45,56 +45,47 @@ svec::svec (int dim, scalar* a)     // conversion constructor
 vec svec::as_vec( ) const
 {
   vec v(d); // initializes to 0
-  for(map<int,scalar>::const_iterator wi=entries.begin(); 
-      wi!=entries.end(); wi++)
-    v[wi->first] = wi->second;
+  for (const auto& wi : entries)
+    v[wi.first] = wi.second;
   return v;
 }
 
 scalar svec::elem(int i)  const   // returns i'th entry
 {
-  map<int,scalar>::const_iterator vi = entries.find(i);
+  auto vi = entries.find(i);
   if(vi==entries.end()) return 0;
   return vi->second;
 }
 
-void svec::set(int i, scalar a) 
+void svec::set(int i, scalar a)
 {
   if(a) entries[i]=a;
-#if(0) // 
-  else // a==0 so kill the i'entry if it is there
-    { 
-      map<int,scalar>::iterator vi = entries.find(i);   
-      if(vi!=entries.end()) entries.erase(vi);
-    }
-#endif
 }
 
 void svec::erase(int i)
 {
-  map<int,scalar>::iterator vi = entries.find(i);   
-  if(vi==entries.end()) 
+  auto vi = entries.find(i);
+  if(vi==entries.end())
     {
       cerr<<"Error in svec::erase(): cannot delete missing entry #"<<i
-	  <<" from v = "<<(*this)<<endl; 
-    } 
+	  <<" from v = "<<(*this)<<endl;
+    }
   else entries.erase(vi);
 }
 
 std::set<int> svec::support() const
 {
   std::set<int> ans;
-  map<int,scalar>::const_iterator vi;
-  for(vi=entries.begin(); vi!=entries.end(); vi++) 
-    ans.insert(vi->first);
+  for (const auto& vi : entries)
+    ans.insert(vi.first);
   return ans;
 }
 
 void svec::add(int i, scalar a)   // adds a to i'th entry
 {
   if(!a) return;
-  map<int,scalar>::iterator vi = entries.find(i);   
-  if(vi==entries.end()) 
+  auto  vi = entries.find(i);
+  if(vi==entries.end())
     entries[i]=a;
   else
     {
@@ -107,8 +98,8 @@ void svec::add(int i, scalar a)   // adds a to i'th entry
 void svec::sub(int i, scalar a)   // subtracts a from i'th entry
 {
   if(!a) return;
-  map<int,scalar>::iterator vi = entries.find(i);   
-  if(vi==entries.end()) 
+  auto vi = entries.find(i);
+  if(vi==entries.end())
     entries[i]=-a;
   else
     {
@@ -122,8 +113,8 @@ void svec::add_mod_p(int i, scalar a, const scalar& p)
 {
   a=xmod(a,p);
   if(!a) return;
-  map<int,scalar>::iterator vi = entries.find(i);   
-  if(vi==entries.end()) 
+  auto vi = entries.find(i);
+  if(vi==entries.end())
     entries[i]=a;
   else
     {
@@ -136,8 +127,8 @@ void svec::sub_mod_p(int i, scalar a, const scalar& p)
 {
   a=xmod(-a,p);
   if(!a) return;
-  map<int,scalar>::iterator vi = entries.find(i);   
-  if(vi==entries.end()) 
+  auto vi = entries.find(i);
+  if(vi==entries.end())
     entries[i]=a;
   else
     {
@@ -154,8 +145,8 @@ svec& svec::operator+=(const svec& w)
       cerr << "Incompatible svecs in svec::operator+=()"<<endl;
       return *this;
     }
-  map<int,scalar>::const_iterator  wi=w.entries.begin();
-  map<int,scalar>::iterator vi=entries.begin();
+  auto wi=w.entries.begin();
+  auto vi=entries.begin();
   while(wi!=w.entries.end())
     {
       if(vi==entries.end())
@@ -164,17 +155,17 @@ svec& svec::operator+=(const svec& w)
 	    {
 	      entries[wi->first]=wi->second;
 	      wi++;
-	    } 	    
+	    }
 	}
       else
 	{
-	  if((vi->first)<(wi->first)) {vi++;} 
+	  if((vi->first)<(wi->first)) {vi++;}
 	  else
-	    if((wi->first)<(vi->first)) 
+	    if((wi->first)<(vi->first))
 	      {
 		entries[wi->first]=wi->second;
 		wi++;
-	      } 
+	      }
 	    else
 	      {
 		scalar sum = (vi->second) + (wi->second);
@@ -194,8 +185,8 @@ svec& svec::operator-=(const svec& w)
       cerr << "Incompatible svecs in svec::operator-=()"<<endl;
       return *this;
     }
-  map<int,scalar>::const_iterator  wi=w.entries.begin();
-  map<int,scalar>::iterator vi=entries.begin();
+  auto wi=w.entries.begin();
+  auto vi=entries.begin();
   while(wi!=w.entries.end())
     {
       if(vi==entries.end())
@@ -204,17 +195,17 @@ svec& svec::operator-=(const svec& w)
 	    {
 	      entries[wi->first]=-wi->second;
 	      wi++;
-	    } 	    
+	    }
 	}
       else
 	{
-	  if((vi->first)<(wi->first))  {vi++;} 
+	  if((vi->first)<(wi->first))  {vi++;}
 	  else
-	    if((wi->first)<(vi->first)) 
+	    if((wi->first)<(vi->first))
 	      {
 		entries[wi->first]=-wi->second;
 		wi++;
-	      } 
+	      }
 	    else
 	      {
 		scalar diff = (vi->second) - (wi->second);
@@ -235,8 +226,8 @@ svec& svec::add_scalar_times(const svec& w, scalar a)
       return *this;
     }
   if(a==0) return *this;
-  map<int,scalar>::const_iterator  wi=w.entries.begin();
-  map<int,scalar>::iterator vi=entries.begin();
+  auto wi=w.entries.begin();
+  auto vi=entries.begin();
   while(wi!=w.entries.end())
     {
       if(vi==entries.end())
@@ -245,17 +236,17 @@ svec& svec::add_scalar_times(const svec& w, scalar a)
 	    {
 	      entries[wi->first]=a*(wi->second);
 	      wi++;
-	    } 	    
+	    }
 	}
       else
 	{
-	  if((vi->first)<(wi->first)) {vi++;} 
+	  if((vi->first)<(wi->first)) {vi++;}
 	  else
-	    if((wi->first)<(vi->first)) 
+	    if((wi->first)<(vi->first))
 	      {
 		entries[wi->first]=a*(wi->second);
 		wi++;
-	      } 
+	      }
 	    else
 	      {
 		scalar sum = (vi->second) + a* (wi->second);
@@ -271,21 +262,27 @@ svec& svec::add_scalar_times(const svec& w, scalar a)
 svec& svec::operator*=(scalar scal)
 {
   //  if(scal==0) cout<<"Attempt to multiply svec by 0\n"<<endl;
-  for( map<int,scalar>::iterator vi=entries.begin(); 
-       vi != entries.end(); vi++)
-    (vi->second)*=scal;
+  for ( auto& vi : entries)
+    (vi.second)*=scal;
   return *this;
 }
 
 void svec::reduce_mod_p(const scalar& p)
 {
   scalar a;
-  map<int,scalar>::iterator vi, vi2;
-  for( vi=entries.begin(); vi != entries.end(); )
+  auto vi = entries.begin();
+  while( vi != entries.end() )
     {
       a = mod(vi->second,p);
-      if(a) {(vi->second)=a; vi++;}
-      else {vi2=vi; vi++; entries.erase(vi2->first);}
+      if(a)
+        {
+          (vi->second)=a;
+          vi++;
+        }
+      else
+        {
+          vi = entries.erase(vi);
+        }
     }
 }
 
@@ -293,9 +290,8 @@ svec& svec::mult_by_scalar_mod_p(scalar scal, const scalar& p)
 {
   //  if(xmod(scal,p)==0) cout<<"Attempt to multiply svec by 0\n"<<endl;
   if(scal!=1)
-    for( map<int,scalar>::iterator vi=entries.begin(); 
-	 vi != entries.end(); vi++)
-      (vi->second)=xmodmul(vi->second,scal,p);
+    for( auto& vi : entries)
+      (vi.second)=xmodmul(vi.second,scal,p);
   return *this;
 }
 
@@ -307,8 +303,8 @@ svec& svec::add_scalar_times_mod_p(const svec& w, scalar a, const scalar& p)
       return *this;
     }
   if(a==0) return *this;
-  map<int,scalar>::const_iterator  wi=w.entries.begin();
-  map<int,scalar>::iterator vi=entries.begin();
+  auto wi=w.entries.begin();
+  auto vi=entries.begin();
   while(wi!=w.entries.end())
     {
       if(vi==entries.end())
@@ -317,17 +313,17 @@ svec& svec::add_scalar_times_mod_p(const svec& w, scalar a, const scalar& p)
 	    {
 	      entries[wi->first]=xmodmul(a,(wi->second),p);
 	      wi++;
-	    } 	    
+	    }
 	}
       else
 	{
-	  if((vi->first)<(wi->first)) {vi++;} 
+	  if((vi->first)<(wi->first)) {vi++;}
 	  else
-	    if((wi->first)<(vi->first)) 
+	    if((wi->first)<(vi->first))
 	      {
 		entries[wi->first]=xmodmul(a,(wi->second),p);
 		wi++;
-	      } 
+	      }
 	    else
 	      {
 		scalar sum = xmod((vi->second) + xmodmul(a, (wi->second),p),p);
@@ -352,8 +348,8 @@ svec& svec::add_scalar_times_mod_p(const svec& w, scalar a, std::set<int>& ons,
       return *this;
     }
   if(a==0) return *this;
-  map<int,scalar>::const_iterator  wi=w.entries.begin();
-  map<int,scalar>::iterator vi=entries.begin();
+  auto wi=w.entries.begin();
+  auto vi=entries.begin();
   while(wi!=w.entries.end())
     {
       if(vi==entries.end())
@@ -363,18 +359,18 @@ svec& svec::add_scalar_times_mod_p(const svec& w, scalar a, std::set<int>& ons,
 	      entries[wi->first]=xmodmul(a,(wi->second),p);
 	      ons.insert(wi->first);
 	      wi++;
-	    } 	    
+	    }
 	}
       else
 	{
-	  if((vi->first)<(wi->first)) {vi++;} 
+	  if((vi->first)<(wi->first)) {vi++;}
 	  else
-	    if((wi->first)<(vi->first)) 
+	    if((wi->first)<(vi->first))
 	      {
 		entries[wi->first]=xmodmul(a,(wi->second),p);
 		ons.insert(wi->first);
 		wi++;
-	      } 
+	      }
 	    else
 	      {
 		scalar sum = xmod((vi->second) + xmodmul(a, (wi->second),p),p);
@@ -395,8 +391,8 @@ svec& svec::add_mod_p(const svec& w, const scalar& p)
       cerr << "Incompatible svecs in svec::add_scalar_times()"<<endl;
       return *this;
     }
-  map<int,scalar>::const_iterator  wi=w.entries.begin();
-  map<int,scalar>::iterator vi=entries.begin();
+  auto wi=w.entries.begin();
+  auto vi=entries.begin();
   while(wi!=w.entries.end())
     {
       if(vi==entries.end())
@@ -405,7 +401,7 @@ svec& svec::add_mod_p(const svec& w, const scalar& p)
 	    {
 	      entries[wi->first]=wi->second;
 	      wi++;
-	    } 	    
+	    }
 	}
       else
 	{
@@ -435,8 +431,8 @@ svec& svec::sub_mod_p(const svec& w, const scalar& p)
       cerr << "Incompatible svecs in svec::add_scalar_times()"<<endl;
       return *this;
     }
-  map<int,scalar>::const_iterator  wi=w.entries.begin();
-  map<int,scalar>::iterator vi=entries.begin();
+  auto wi=w.entries.begin();
+  auto vi=entries.begin();
   while(wi!=w.entries.end())
     {
       if(vi==entries.end())
@@ -445,17 +441,17 @@ svec& svec::sub_mod_p(const svec& w, const scalar& p)
 	    {
 	      entries[wi->first]=-(wi->second);
 	      wi++;
-	    } 	    
+	    }
 	}
       else
 	{
-	  if((vi->first)<(wi->first)) {vi++;} 
+	  if((vi->first)<(wi->first)) {vi++;}
 	  else
-	    if((wi->first)<(vi->first)) 
+	    if((wi->first)<(vi->first))
 	      {
 		entries[wi->first]=-(wi->second);
 		wi++;
-	      } 
+	      }
 	    else
 	      {
 		scalar sum = xmod((vi->second) - (wi->second),p);
@@ -470,13 +466,12 @@ svec& svec::sub_mod_p(const svec& w, const scalar& p)
 
 svec& svec::operator/=(scalar scal)
 {
-  if(scal==0) 
+  if(scal==0)
     {
       cerr<<"Attempt to divide svec by 0"<<endl;
     }
-  for( map<int,scalar>::iterator vi=entries.begin(); 
-       vi != entries.end(); vi++)
-    (vi->second)/=scal;
+  for( auto& vi : entries)
+    (vi.second)/=scal;
   return *this;
 }
 
@@ -485,11 +480,10 @@ svec& svec::operator/=(scalar scal)
 int eqmodp(const svec& v1, const svec& v2, const scalar& p)
 {
   if(v1.d!=v2.d) return 0;
-  map<int,scalar>::const_iterator vi;
-  for(vi=v1.entries.begin(); vi!=v1.entries.end(); vi++)
-    if(xmod((vi->second)-(v2.elem(vi->first)),p)!=0) return 0;
-  for(vi=v2.entries.begin(); vi!=v2.entries.end(); vi++)
-    if(xmod((vi->second)-(v1.elem(vi->first)),p)!=0) return 0;
+  for( const auto& vi : v1.entries)
+    if(xmod((vi.second)-(v2.elem(vi.first)),p)!=0) return 0;
+  for( const auto& vi : v2.entries)
+    if(xmod((vi.second)-(v1.elem(vi.first)),p)!=0) return 0;
   return true;
 }
 
@@ -502,9 +496,8 @@ int operator==(const svec& v1, const vec& v2)
 
 ostream& operator << (ostream& s, const svec& v)
 {
-  map<int,scalar>::const_iterator vi;
   s<<"[";
-  for(vi=v.entries.begin(); vi!=v.entries.end(); vi++)
+  for(auto vi=v.entries.begin(); vi!=v.entries.end(); vi++)
     {
       if(vi!=v.entries.begin()) s<<",";
       s << "("<<vi->first<<":"<<vi->second<<")";
@@ -517,7 +510,7 @@ scalar operator*(const svec& v, const svec& w) //dot prod
 {
   scalar ans=0;
   if((v.entries.size()==0)||(w.entries.size()==0)) return ans;
-  map<int,scalar>::const_iterator vi=v.entries.begin(), wi=w.entries.begin();
+  auto vi=v.entries.begin(), wi=w.entries.begin();
   while((vi!=v.entries.end())&&(wi!=w.entries.end()))
     {
       if((vi->first)<(wi->first)) {vi++;} else
@@ -534,18 +527,16 @@ scalar operator*(const svec& v, const vec& w) //dot prod
 {
   scalar ans=0;
   if((v.entries.size()==0)) return ans;
-  map<int,scalar>::const_iterator vi;
-  for(vi=v.entries.begin(); vi!=v.entries.end(); vi++)
-    ans+=(vi->second)* w[vi->first];
+  for( const auto& vi : v.entries)
+    ans+=(vi.second)* w[vi.first];
   return ans;
 }
 
 scalar dotmodp(const svec& v, const vec& w, scalar pr)
 {
   scalar ans=0;
-  map<int,scalar>::const_iterator vi;
-  for(vi=v.entries.begin(); vi!=v.entries.end(); vi++)
-    ans=xmod(ans+xmodmul(vi->second,w[vi->first],pr),pr);
+  for( const auto& vi : v.entries)
+    ans=xmod(ans+xmodmul(vi.second,w[vi.first],pr),pr);
   return ans;
 }
 
@@ -553,7 +544,7 @@ scalar dotmodp(const svec& v, const svec& w, scalar pr)
 {
   scalar ans=0;
   if((v.entries.size()==0)||(w.entries.size()==0)) return ans;
-  map<int,scalar>::const_iterator vi=v.entries.begin(), wi=w.entries.begin();
+  auto vi=v.entries.begin(), wi=w.entries.begin();
   while((vi!=v.entries.end())&&(wi!=w.entries.end()))
     {
       if((vi->first)<(wi->first)) {vi++;} else
@@ -569,9 +560,12 @@ scalar dotmodp(const svec& v, const svec& w, scalar pr)
 scalar content(const svec& v)
 {
   scalar ans=0;
-  map<int,scalar>::const_iterator vi;
-  for(vi=v.entries.begin(); (ans!=1)&&(vi!=v.entries.end()); vi++)
-    ans=gcd(ans,vi->second);
+  for( const auto & vi : v.entries)
+    {
+      ans=gcd(ans,vi.second);
+      if (ans==1)
+        break;
+    }
   return ans;
 }
 
diff --git a/libsrc/symb.cc b/libsrc/symb.cc
index dcd5db2..e158629 100644
--- a/libsrc/symb.cc
+++ b/libsrc/symb.cc
@@ -108,23 +108,13 @@ void symblist::add(const symb& s, long start)
  }
 }
 
-/*
-long symblist::index(const symb& s, long start) const
-{
- long i,ans;
- for (i=start,ans=-1; ((i<num)&&(ans==-1)); i++) if (list[i]==s) ans=i;
- return ans;
-}
-*/
-
 long symblist::index(const symb& s, long start) const
 {
   //  cout<<"index of "<<s;
  symb ss = s.normalize();
- long c = ss.cee(), d=ss.dee(); 
- map<pair<long,long>,long>::const_iterator 
-   j = hashtable.find(pair<long,long>(c,d)); 
- if(j==hashtable.end()) 
+ long c = ss.cee(), d=ss.dee();
+ auto j = hashtable.find(pair<long,long>(c,d));
+ if(j==hashtable.end())
    return -1;
  // cout<<" is "<<j->second<<endl;
  return j->second;
diff --git a/libsrc/timer.cc b/libsrc/timer.cc
index b4addf4..68ef5b0 100644
--- a/libsrc/timer.cc
+++ b/libsrc/timer.cc
@@ -172,13 +172,9 @@ void timer::stop( string name ) {
  * Stop all timers.
  */
 void timer::stopAll() {
-  timers::iterator it;
-
   // TODO Check if timer is active via flags
-
-  for( it = times_.begin(); it != times_.end(); it++ ) {
-    stop( it -> first );
-  }
+  for ( const auto& t : times_)
+    stop(t.first);
 }
 
 /**
@@ -226,11 +222,8 @@ void timer::show( int nline, string name, int idx1, int idx2 ) {
  * Loop through all timers, and write out.
  */
 void timer::showAll( int nline ) {
-  timers::iterator it;
- 
-  for( it = times_.begin(); it != times_.end(); it++ ) {
-    show( nline, it -> first );
-  }
+  for (const auto& t : times_)
+    show( nline, t.first );
 }
 
 /**
@@ -248,11 +241,8 @@ void timer::clear( string name ) {
  * Removes all times from all subtimers.
  */
 void timer::clearAll() {
-  timers::iterator it;
-
-  for( it = times_.begin(); it != times_.end(); it++ ) {
-    clear( it -> first );
-  }
+  for (const auto& t : times_)
+    clear( t.first );
 }
 
 /**
@@ -261,13 +251,9 @@ void timer::clearAll() {
  * List names of all available timers
  */
 void timer::list() {
-  timers::iterator it;
-
   string message;
-
-  for( it = times_.begin(); it != times_.end(); it++ ) {
-    message += it -> first + " ";
-  }
+  for (const auto& t : times_)
+    message += t.first + " ";
 
   message += "\n";
   s_ -> write( message.c_str(), message.size() );
@@ -290,14 +276,9 @@ int timer::count( string name ) {
  * Return total time of a given timer.
  */
 double timer::total( string name ) {
-  // Cache size of vector
-  vector<double>::size_type size = times_[name].size();
   double total = 0;
-  
-  for( int i = 0; i < (int)size; i++ ) {
-    total += times_[name][i];
-  }
-
+  for ( const auto& t : times_[name])
+    total += t;
   return total;
 }
 
@@ -307,10 +288,7 @@ double timer::total( string name ) {
  * Return average time of a given timer.
  */
 double timer::average( string name ) {
-  // Cache size of vector
-  vector<double>::size_type size = times_[name].size();
-
-  return total( name ) / size;
+  return total( name ) /  times_[name].size();
 }
 
 /**
diff --git a/libsrc/xsplit.cc b/libsrc/xsplit.cc
index 7df057a..fba7087 100644
--- a/libsrc/xsplit.cc
+++ b/libsrc/xsplit.cc
@@ -552,7 +552,7 @@ void form_finder::find( ff_data &data ) {
 
   // The recursive part:
   vector<long> t_eigs = h->eigrange(depth);
-  vector<long>::const_iterator apvar = t_eigs.begin();
+  auto apvar = t_eigs.begin();
 
   stringstream t_eigs_ss;
   std::copy(t_eigs.begin(),t_eigs.end(),ostream_iterator<long>(t_eigs_ss," "));
diff --git a/libsrc/xsplit_data.cc b/libsrc/xsplit_data.cc
index 3eaa46a..b8302e2 100644
--- a/libsrc/xsplit_data.cc
+++ b/libsrc/xsplit_data.cc
@@ -179,14 +179,8 @@ ff_data* ff_data::child( long eig ) {
  * Returns number of completed children for current node.
  */
 int ff_data::numCompleteChildren() {
-  int completeCount = 0;
-
-  vector< childstatus >::iterator it;
-  for( it = completedChildren_.begin(); it != completedChildren_.end(); it++ ) {
-    if( *it != NOT_COMPLETE ) completeCount++; 
-  }
-
-  return completeCount;
+  return std::count_if(completedChildren_.begin(), completedChildren_.end(),
+                       [](childstatus s) {return s!=NOT_COMPLETE;});
 }
 
 /**