diff --git a/libsrc/conic.cc b/libsrc/conic.cc
index ba0f39b..949970f 100644
--- a/libsrc/conic.cc
+++ b/libsrc/conic.cc
@@ -585,10 +585,9 @@ void testmodsqrt()
   cout << "Enter a modulus m: ";
   cin >> m; mm=m;
   vector<bigint> plist=pdivs(mm);
-  int* flag = new int[m];
-  for(i=0; i<m; i++) flag[i]=0;
+  vector<int> flag(m,0);
   for(i=0; i<=m/2; i++) flag[(i*i)%m]=1;
-  for(i=0; i<m; i++) 
+  for(i=0; i<m; i++)
     {
       a=i; //      cout<<"a = "<<a;
 
@@ -600,12 +599,12 @@ void testmodsqrt()
 	  if((x*x-a)%mm==0) cout<<" --checks.";
 	  else cout << "--WRONG!";
 	}
-      else 
+      else
 	{
 	  cout << "\tNo solution";
 	}
 #endif
-      if(res!=flag[i]) 
+      if(res!=flag[i])
 	{
 	  cout << "WRONG ANSWER for a="<<a<<endl;
 	  ok=0;
diff --git a/libsrc/cperiods.cc b/libsrc/cperiods.cc
index a1e4791..96486d8 100644
--- a/libsrc/cperiods.cc
+++ b/libsrc/cperiods.cc
@@ -104,7 +104,7 @@ void eiperiods(bigcomplex e1, bigcomplex e2, bigcomplex e3,
 
 //#define DEBUG_CUBIC
 
-bigcomplex* solve_nonsingular_cubic(const bigint& c1, const bigint& c2, const bigint& c3) 
+vector<bigcomplex> solve_nonsingular_cubic(const bigint& c1, const bigint& c2, const bigint& c3) 
 //Returns an array of 3 complex roots.
 {
 #ifdef DEBUG_CUBIC
@@ -120,7 +120,7 @@ bigcomplex* solve_nonsingular_cubic(const bigint& c1, const bigint& c2, const bi
   bigint p3=  3*c2 - c1*c1;
   bigint q = c1*(2*sqr(c1)-9*c2)+27*c3;
   bigfloat rq=I2bigfloat(q), rp3=I2bigfloat(p3);
-  bigcomplex *roots = new bigcomplex[3];
+  vector<bigcomplex> roots(3);
   long i;
  
 #ifdef DEBUG_CUBIC
@@ -233,7 +233,7 @@ void getei(const Curvedata& E, bigcomplex& e1, bigcomplex& e2, bigcomplex& e3)
 #ifdef DEBUG
   cout<<"Solving monic cubic with coeffs "<<b2<<","<<(8*b4)<<","<<16*b6<<endl;
 #endif
-  bigcomplex* ei = solve_nonsingular_cubic(b2,8*b4,16*b6);
+  vector<bigcomplex> ei = solve_nonsingular_cubic(b2,8*b4,16*b6);
 #ifdef DEBUG
   cout<<"ei = "<<ei[0]<<","<<ei[1]<<","<<ei[2]<<endl;
 #endif
@@ -243,7 +243,6 @@ void getei(const Curvedata& E, bigcomplex& e1, bigcomplex& e2, bigcomplex& e3)
   cout<<"After rescaling,\n";
   cout<<"ei = "<<e1<<","<<e2<<","<<e3<<endl;
 #endif
-  delete [] ei;
 }
 
 Cperiods::Cperiods(const Curvedata& E)
diff --git a/libsrc/eclib/cperiods.h b/libsrc/eclib/cperiods.h
index 33d5c53..527c95b 100644
--- a/libsrc/eclib/cperiods.h
+++ b/libsrc/eclib/cperiods.h
@@ -146,8 +146,8 @@ class Cperiods {
   vector<bigcomplex> ellztopoint(const bigcomplex& z, const bigcomplex& a1, const bigcomplex& a2, const bigcomplex& a3);
 };  // end of class Cperiods def'n
 
-bigcomplex* solve_nonsingular_cubic(const bigint& c1, const bigint& c2, const bigint& c3); //Returns an array
-// Gets the 3 2-division points given the coefficients 
+vector<bigcomplex> solve_nonsingular_cubic(const bigint& c1, const bigint& c2, const bigint& c3);
+// Gets the 3 2-division points given the coefficients
 void getei(const Curvedata& E, bigcomplex& e1, bigcomplex& e2, bigcomplex& e3);
 // Reorders 3 complex nos so real parts are decreasing
 void reorder1(bigcomplex& a, bigcomplex& b, bigcomplex& c);
diff --git a/libsrc/eclib/htconst.h b/libsrc/eclib/htconst.h
index 0bddb26..edb135b 100644
--- a/libsrc/eclib/htconst.h
+++ b/libsrc/eclib/htconst.h
@@ -39,7 +39,6 @@ inline double height_constant(const Curvedata& CD)
 }
 double egr_height_constant(const Curvedata& CD);
 
-bigfloat lower_height_bound_alt(const Curvedata& CD);
 bigfloat lower_height_bound_search(const Curvedata& CD, const bigfloat& reg);
 bigfloat Gamma_n(long n); // Gamma(n) = (n-1)!
 bigfloat Gamma_n_plus_half(long n); // Gamma(n+1/2) = (2n)!sqrt(pi) / (4^n*n!)
diff --git a/libsrc/eclib/mglobsol.h b/libsrc/eclib/mglobsol.h
index 890f0cf..d73c01c 100644
--- a/libsrc/eclib/mglobsol.h
+++ b/libsrc/eclib/mglobsol.h
@@ -41,22 +41,19 @@ class quartic_sieve : public point_processor {
   int verbose, easy, use_stoll;
   long ulim;
   int num_aux;
-  long* auxs;
-  int** xgood_mod_aux;
-  int** squares;
-  long* umod;
-  long nwprimes; long* wprimes;  long* uprimes;
+  vector<long> auxs, umod, wprimes, uprimes;
+  vector<vector<int>> xgood_mod_aux, squares;
+  long nwprimes;
   long npoints, maxnpoints;
   int process(const bigint& x, const bigint& y, const bigint& z)
-  {pu=x; pv=y; pw=z; npoints++; 
-  //cout<<"[x,y,z]=["<<x<<","<<y<<","<<z<<"]\n"; 
+  {pu=x; pv=y; pw=z; npoints++;
+  //cout<<"[x,y,z]=["<<x<<","<<y<<","<<z<<"]\n";
   return (npoints>=maxnpoints);
   }
  // (x,y,z) as returned by ms's sieve; the point is (x/z,y/z^2)
 public:
   quartic_sieve(void) {;}
   quartic_sieve(quartic * gg, int moduli_option=2, int verb=0); 
-  ~quartic_sieve();
   long search(double h_lim, long maxnpts=1, int posxonly=0);
   long stoll_search(double h_lim, int posxonly=0);
   long search_range(int lower, bigfloat lower_bound, 
diff --git a/libsrc/eclib/mlocsol.h b/libsrc/eclib/mlocsol.h
index 99b06d4..1f211e3 100644
--- a/libsrc/eclib/mlocsol.h
+++ b/libsrc/eclib/mlocsol.h
@@ -52,7 +52,7 @@ int locallysoluble(const bigint& a, const bigint& c, const bigint& e,
 
 /* Samir Siksek's Local Solubility Test for odd p */
 
-int local_sol(const bigint& p,bigint *c, int verbose=0);
+int local_sol(const bigint& p, vector<bigint> c, int verbose=0);
 
  // Checks for solublility in Qp 
 int new_qpsoluble(const quartic& g, const bigint& p, int verbose=0);
diff --git a/libsrc/eclib/quadratic.h b/libsrc/eclib/quadratic.h
index 14c5492..9cb70dd 100644
--- a/libsrc/eclib/quadratic.h
+++ b/libsrc/eclib/quadratic.h
@@ -31,26 +31,16 @@
 
 class quadratic {
   friend class unimod;
-private: 
-  bigint * coeffs;  // will always have length 3
-  // init just allocates memory
-  void init();
+private:
+  vector<bigint> coeffs;
 public:
-  void set(long a, long b, long c) 
-    {coeffs[0]=a; coeffs[1]=b; coeffs[2]=c;}
-  void set(const  bigint& a, const bigint& b, const bigint& c) 
-    {coeffs[0]=a; coeffs[1]=b; coeffs[2]=c;}
-  void set(const  bigint* q)
-    {coeffs[0]=q[0]; coeffs[1]=q[1]; coeffs[2]=q[2];}
-  void set(const  quadratic& q)
-    {coeffs[0]=q.coeffs[0]; coeffs[1]=q.coeffs[1]; coeffs[2]=q.coeffs[2];}
-  quadratic() {init(); set(0,0,0);}
-  ~quadratic();
-  quadratic(const  bigint& a, const bigint& b, const bigint& c) {init(); set(a,b,c);}
-  quadratic(long a, long b, long c) {init(); set(a,b,c);}
-  quadratic(const  bigint* q) {init(); set(q);}
-  quadratic(const  quadratic& q) {init(); set(q);}
-  void operator=(const quadratic& q) {set(q);}
+  void set(const  bigint& a, const bigint& b, const bigint& c) {coeffs = {a, b, c};}
+  quadratic() { bigint zero(0); coeffs={zero,zero,zero};}
+  quadratic(const  bigint& a, const bigint& b, const bigint& c) {coeffs = {a, b, c};}
+  quadratic(long a, long b, long c)  {coeffs = {BIGINT(a), BIGINT(b), BIGINT(c)};}
+  quadratic(const  bigint* q) {coeffs = {q[0], q[1], q[2]};}
+  quadratic(const  quadratic& q) {coeffs = q.coeffs;}
+  void operator=(const quadratic& q)  {coeffs = q.coeffs;}
   bigint eval(const bigint& x, const bigint& z) const
     {return coeffs[0]*sqr(x) + coeffs[1]*x*z + coeffs[2]*sqr(z);}
   bigint operator()(const bigint& x, const bigint& z) const
@@ -59,10 +49,9 @@ class quadratic {
     {return coeffs[0]*sqr(x) + coeffs[1]*x + coeffs[2];}
   bigint operator()(const bigint& x) const
     {return coeffs[0]*sqr(x) + coeffs[1]*x + coeffs[2];}
-  bigint coeff(int i); 
-  bigint operator[](int i) const;
-  void set_coeff(int i, const bigint& a)
-    {if((i>=0)&&(i<=2)) coeffs[i]=a;}
+  bigint coeff(int i) const {return coeffs[i];}
+  bigint operator[](int i) const  {return coeffs[i];}
+  void set_coeff(int i, const bigint& a) {if((i>=0)&&(i<=2)) coeffs[i]=a;}
   friend bigint resultant(const quadratic& q1, const quadratic& q2);
   bigint disc() const  {return sqr(coeffs[1])-4*coeffs[0]*coeffs[2];}
   void output(ostream& os=cout) const
diff --git a/libsrc/htconst.cc b/libsrc/htconst.cc
index 2a5ea63..4e7e442 100644
--- a/libsrc/htconst.cc
+++ b/libsrc/htconst.cc
@@ -68,14 +68,14 @@ double silverman_bound(const Curvedata& CD)
   double realjay;
   double hjay = hj(CD,realjay);
 
-// NB the constant 1.922 = 2*0.961 below is from Bremner's correction 
+// NB the constant 1.922 = 2*0.961 below is from Bremner's correction
 // to Silverman's paper; Silverman has 0.973 givin 2*0.973 = 1.946.
 
   double mu = 1.922  + hjay/12
                      + log(fabs(I2double(delta)))/6
                      + logplus(realjay)/6
 	             + logplus(I2double(b2)/12);
-  
+
   if(b2!=0) mu += log2;
 
   return mu;
@@ -168,9 +168,8 @@ bigfloat old_calc_dv_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat&
 
 double cps_real(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6, const bigfloat& b8)
 {
-  bigfloat htc=to_bigfloat(0);
-  bigfloat dv=to_bigfloat(0);
-  bigfloat dvd=to_bigfloat(0);
+  bigfloat zero=to_bigfloat(0);
+  bigfloat htc=zero, dv=zero, dvd=zero;
   long original_prec, prec;
 #ifdef MPFP
   prec = original_prec = bit_precision();
@@ -209,17 +208,6 @@ double cps_real(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6, cons
     }
 #endif
 
-#ifdef TEST_CPS
-  bigfloat del = -b2*b2*b8-8*b4*b4*b4-27*b6*b6+9*b2*b4*b6;
-  bigfloat dv2=old_calc_dv_inf(b2,b4,b6,b8,del);
-  bigfloat dvd2=old_calc_dvd_inf(b2,b4,b6,b8,del); 
-  bigfloat diff = dv-dv2;
-  if(!is_zero(diff)) 
-    cout<<"old_calc_dv gives "<<dv2<<"\nwhile new gives "<<dv<<endl<<"diff = "<<diff<<endl;
-  diff = dvd-dvd2;
-  if(!is_zero(diff)) 
-    cout<<"old_calc_dvd gives "<<dvd2<<"\nwhile new gives "<<dvd<<endl<<"diff = "<<diff<<endl;
-#endif
 #ifdef DEBUG_CPS
   cout << "dv=" << dv << endl;
   cout << "dvd=" << dvd << endl;
@@ -280,22 +268,19 @@ vector<bigfloat> roots11( const vector<bigfloat>& coeff );
 
 inline vector<bigfloat> set_coeff(const bigfloat& c0, const bigfloat& c1, const bigfloat& c2, const bigfloat& c3, const bigfloat& c4)
 {
-  vector<bigfloat> coeff(5);
-  coeff[0]=c0;  coeff[1]=c1;  coeff[2]=c2;  coeff[3]=c3;  coeff[4]=c4;
+  vector<bigfloat> coeff = {c0,c1,c2,c3,c4};
   return coeff;
 }
 
 inline vector<bigfloat> set_coeff(const bigfloat& c0, const bigfloat& c1, const bigfloat& c2, const bigfloat& c3)
 {
-  vector<bigfloat> coeff(4);
-  coeff[0]=c0;  coeff[1]=c1;  coeff[2]=c2;  coeff[3]=c3;
+  vector<bigfloat> coeff = {c0,c1,c2,c3};
   return coeff;
 }
 
 inline vector<bigfloat> set_coeff(const bigfloat& c0, const bigfloat& c1, const bigfloat& c2)
 {
-  vector<bigfloat> coeff(3);
-  coeff[0]=c0;  coeff[1]=c1;  coeff[2]=c2;
+  vector<bigfloat> coeff = {c0,c1,c2};
   return coeff;
 }
 
@@ -310,9 +295,10 @@ vector<bigfloat> reals_in_11 ( vector<bigcomplex>& v);
 //  Procedure to calculate dv' for the infinite prime
 
 bigfloat calc_dvd_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6, const bigfloat& b8)
-{ 
+{
+  bigfloat zero=to_bigfloat(0), one=to_bigfloat(1), two=to_bigfloat(2), three=to_bigfloat(3), four=to_bigfloat(4);
   bigfloat rx;
-  bigfloat dvd=to_bigfloat(0),x,x2,Fx,Gx;
+  bigfloat dvd=zero,x,x2,Fx,Gx;
 
 #ifdef DEBUG_CPS
   cout<<"\nIn calc_dvd_inf"<<endl;
@@ -323,9 +309,9 @@ bigfloat calc_dvd_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6
   std::set<bigfloat> crit_pts; // use a set, so we don't get any duplicates
   std::set<bigfloat> F_roots;
 
-  crit_pts.insert(to_bigfloat(1));
-  crit_pts.insert(to_bigfloat(-1));
-  
+  crit_pts.insert(one);
+  crit_pts.insert(-one);
+
 #ifdef DEBUG_CPS
   cout<<"crit_pts = "<<crit_pts<<endl;
 #endif
@@ -334,8 +320,7 @@ bigfloat calc_dvd_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6
   // Keep the roots of F separate (see comment below)
 
   // Roots of F
-  //  rts=roots11(set_reverse_coeff(to_bigfloat(0),to_bigfloat(1),b2/4.0,b4/2.0,b6/4.0));
-  rts=realroots11(set_coeff(b6,2*b4,b2,to_bigfloat(4),to_bigfloat(0)));
+  rts=realroots11(set_coeff(b6,2*b4,b2, four,zero));
   F_roots.insert(rts.begin(),rts.end());
 #ifdef DEBUG_CPS
   cout<<"relevant roots of F: "<<rts<<endl;
@@ -343,8 +328,7 @@ bigfloat calc_dvd_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6
 #endif
 
   // Roots of G
-  //  rts=roots11(set_reverse_coeff(to_bigfloat(1),to_bigfloat(0),-b4,-2*b6,-b8));
-  rts=realroots11(set_coeff(-b8,-2*b6,-b4,to_bigfloat(0),to_bigfloat(1)));
+  rts=realroots11(set_coeff(-b8,-two*b6,-b4,zero,one));
   crit_pts.insert(rts.begin(),rts.end());
 #ifdef DEBUG_CPS
   cout<<"relevant roots of G: "<<rts<<endl;
@@ -352,8 +336,7 @@ bigfloat calc_dvd_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6
 #endif
 
   // Roots of F+G
-  //  rts=roots11(set_reverse_coeff(to_bigfloat(1),to_bigfloat(4),b2-b4,2.0*(b4-b6),b6-b8));
-  rts=realroots11(set_coeff(b6-b8,2.0*(b4-b6),b2-b4,to_bigfloat(4),to_bigfloat(1)));
+  rts=realroots11(set_coeff(b6-b8,two*(b4-b6),b2-b4, four,one));
   crit_pts.insert(rts.begin(),rts.end());
 #ifdef DEBUG_CPS
   cout<<"relevant roots of F+G: "<<rts<<endl;
@@ -361,8 +344,7 @@ bigfloat calc_dvd_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6
 #endif
 
   // Roots of G-F // Change from NPS
-  //  rts=roots11(set_reverse_coeff(to_bigfloat(1),to_bigfloat(-4),-(b2+b4),-2.0*(b4+b6),-(b6+b8)));
-  rts=realroots11(set_coeff(-(b6+b8),-2.0*(b4+b6),-(b2+b4),to_bigfloat(-4),to_bigfloat(1)));
+  rts=realroots11(set_coeff(-(b6+b8),-two*(b4+b6),-(b2+b4), -four,one));
   crit_pts.insert(rts.begin(),rts.end());
 #ifdef DEBUG_CPS
   cout<<"relevant roots of F-G: "<<rts<<endl;
@@ -370,8 +352,7 @@ bigfloat calc_dvd_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6
 #endif
 
   // Roots of G'
-  //  rts=roots11(set_coeff(to_bigfloat(0),-2*b8,-3*b6,-b4,to_bigfloat(0)));
-  rts=realroots11(set_coeff(-2*b8,-3*b6,-b4,to_bigfloat(0)));
+  rts=realroots11(set_coeff(-two*b8,-three*b6,-b4, zero));
   crit_pts.insert(rts.begin(),rts.end());
 #ifdef DEBUG_CPS
   cout<<"relevant roots of G': "<<rts<<endl;
@@ -379,8 +360,7 @@ bigfloat calc_dvd_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6
 #endif
 
   // Roots of F'
-  //  rts=roots11(set_coeff(to_bigfloat(0),2*b6,3*b4,b2,to_bigfloat(2)));
-  rts=realroots11(set_coeff(2*b6,3*b4,b2,to_bigfloat(2)));
+  rts=realroots11(set_coeff(two*b6,three*b4,b2,two));
   crit_pts.insert(rts.begin(),rts.end());
 #ifdef DEBUG_CPS
   cout<<"relevant roots of F': "<<rts<<endl;
@@ -398,7 +378,7 @@ bigfloat calc_dvd_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6
   for (const auto& x : F_roots)
     {
       x2=x*x;
-      Gx=abs(1-b4*x2-2.0*b6*x*x2-b8*x2*x2);
+      Gx=abs(1-b4*x2-two*b6*x*x2-b8*x2*x2);
 #ifdef DEBUG_CPS
       cout<<"x="<<x<<", Gx="<<Gx<<endl;
 #endif
@@ -411,12 +391,12 @@ bigfloat calc_dvd_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6
   for (const auto& x : crit_pts)
     {
       x2=x*x;
-      Fx=(4.0*x+b2*x2+2.0*b4*x*x2+b6*x2*x2);
+      Fx=(four*x+b2*x2+two*b4*x*x2+b6*x2*x2);
 #ifdef DEBUG_CPS
       cout<<"x="<<x<<", Fx="<<Fx<<endl;
 #endif
       if(Fx<0) continue;
-      Gx=abs(1-b4*x2-2.0*b6*x*x2-b8*x2*x2);
+      Gx=abs(one-b4*x2-two*b6*x*x2-b8*x2*x2);
 #ifdef DEBUG_CPS
       cout<<"x="<<x<<", Gx="<<Gx<<endl;
 #endif
@@ -427,7 +407,7 @@ bigfloat calc_dvd_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6
       cout<<"dvd so far ="<<dvd<<endl;
 #endif
     }
-  if(first) return to_bigfloat(-1);
+  if(first) return to_bigfloat(-one);
 
   return dvd;
 
@@ -437,8 +417,9 @@ bigfloat calc_dvd_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6
 
 bigfloat calc_dv_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6, const bigfloat& b8)
 {
+  bigfloat zero=to_bigfloat(0), one=to_bigfloat(1), two=to_bigfloat(2), six=to_bigfloat(6), four=to_bigfloat(4);
   bigfloat rx;
-  bigfloat dv=to_bigfloat(0),x,x2,fx,gx;
+  bigfloat dv=zero,x,x2,fx,gx;
 
 #ifdef DEBUG_CPS
   cout<<"\nIn calc_dv_inf"<<endl;
@@ -449,8 +430,8 @@ bigfloat calc_dv_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6,
   std::set<bigfloat> crit_pts; // use a set, so we don't get any duplicates
   std::set<bigfloat> f_roots;
 
-  crit_pts.insert(to_bigfloat(1));
-  crit_pts.insert(to_bigfloat(-1));
+  crit_pts.insert(to_bigfloat(one));
+  crit_pts.insert(to_bigfloat(-one));
   
 #ifdef DEBUG_CPS
   cout<<"crit_pts = "<<crit_pts<<endl;
@@ -460,8 +441,7 @@ bigfloat calc_dv_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6,
   // Keep the roots of f separate (see comment below)
 
   // Roots of f
-  //  rts=roots11(set_coeff(to_bigfloat(0),to_bigfloat(1),b2/4.0,b4/2.0,b6/4.0));
-  rts=realroots11(set_coeff(to_bigfloat(4),b2,2*b4,b6));
+  rts=realroots11(set_coeff(four,b2,two*b4,b6));
   f_roots.insert(rts.begin(),rts.end());
 #ifdef DEBUG_CPS
   cout<<"relevant roots of f: "<<rts<<endl;
@@ -469,8 +449,7 @@ bigfloat calc_dv_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6,
 #endif
 
   // Roots of g
-  //  rts=roots11(set_coeff(to_bigfloat(1),to_bigfloat(0),-b4,-2*b6,-b8));
-  rts=realroots11(set_coeff(to_bigfloat(1),to_bigfloat(0),-b4,-2*b6,-b8));
+  rts=realroots11(set_coeff(one, zero,-b4,-two*b6,-b8));
   crit_pts.insert(rts.begin(),rts.end());
 #ifdef DEBUG_CPS
   cout<<"relevant roots of g: "<<rts<<endl;
@@ -478,8 +457,7 @@ bigfloat calc_dv_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6,
 #endif
 
   // Roots of f+g
-  //  rts=roots11(set_coeff(to_bigfloat(1),to_bigfloat(4),b2-b4,2.0*(b4-b6),b6-b8));
-  rts=realroots11(set_coeff(to_bigfloat(1),to_bigfloat(4),b2-b4,2.0*(b4-b6),b6-b8));
+  rts=realroots11(set_coeff(one, four, b2-b4,two*(b4-b6),b6-b8));
   crit_pts.insert(rts.begin(),rts.end());
 #ifdef DEBUG_CPS
   cout<<"relevant roots of f+g: "<<rts<<endl;
@@ -487,8 +465,7 @@ bigfloat calc_dv_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6,
 #endif
 
   // Roots of g-f // Change from NPS
-  //  rts=roots11(set_coeff(to_bigfloat(1),to_bigfloat(-4),-(b2+b4),-2.0*(b4+b6),-(b6+b8)));
-  rts=realroots11(set_coeff(to_bigfloat(1),to_bigfloat(-4),-(b2+b4),-2.0*(b4+b6),-(b6+b8)));
+  rts=realroots11(set_coeff(one,-four,-(b2+b4),-two*(b4+b6),-(b6+b8)));
   crit_pts.insert(rts.begin(),rts.end());
 #ifdef DEBUG_CPS
   cout<<"relevant roots of f-g: "<<rts<<endl;
@@ -496,8 +473,7 @@ bigfloat calc_dv_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6,
 #endif
 
   // Roots of g'
-  //  rts=roots11(set_coeff(to_bigfloat(0),to_bigfloat(1),to_bigfloat(0),-b4/2.0,-b6/2.0));
-  rts=realroots11(set_coeff(to_bigfloat(2),to_bigfloat(0),-b4,-b6));
+  rts=realroots11(set_coeff(two,zero,-b4,-b6));
   crit_pts.insert(rts.begin(),rts.end());
 #ifdef DEBUG_CPS
   cout<<"relevant roots of g': "<<rts<<endl;
@@ -505,8 +481,7 @@ bigfloat calc_dv_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6,
 #endif
 
   // Roots of f'
-  //  rts=roots11(set_coeff(to_bigfloat(0),to_bigfloat(0),to_bigfloat(1),b2/6,b4/6));
-  rts=realroots11(set_coeff(to_bigfloat(6),b2,b4));
+  rts=realroots11(set_coeff(six,b2,b4));
   crit_pts.insert(rts.begin(),rts.end());
 #ifdef DEBUG_CPS
   cout<<"relevant roots of f': "<<rts<<endl;
@@ -557,532 +532,8 @@ bigfloat calc_dv_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6,
   return dv;
 }
 
-
-// PAST THIS POINT the Smart method is implemented
-
-int interval_test(const bigfloat& x, const vector<bigfloat> rts, int debug=0);
-
-vector<bigfloat> reals_in_interval ( vector<bigcomplex>& v, const vector<bigfloat> rts);
-
-// Inserts from C in to S the elements which are real and between -1
-// and +1
-void include_real_11(std::set<bigfloat>& S, const vector<bigcomplex>& C);
-
-
-// test if rr lies in one of the first numint intervals out of [i1l.i1r],[i2l.i2r]
-
-int is_in_int(const bigfloat rr,const bigfloat i1l,const bigfloat i1r,
-              const bigfloat i2l,const bigfloat i2r,const long numint)
-{
-if (numint>0)
-  { if (rr<=i1r && rr>=i1l) { return 1; }
-    if (numint==2 && rr<=i2r && rr>=i2l) { return 1; }
-  }
-return 0;
-}
-
-// test if rr lies in one of the numint intervals  [intervals[i][0],intervals[i][1]]
-
-int is_in_int2(const bigfloat rr,bigfloat** intervals,const long numint)
-{
-  long i;
-  for (i=0; i<numint; i++)
-    { if ((rr>=intervals[i][0]) && (rr<=intervals[i][1]))
-	{ return 1; }
-    }
-  return 0;
-}
-
-bigfloat old_calc_dv_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6, const bigfloat& b8, const bigfloat& del)
-{
-  bigfloat rr=to_bigfloat(0),rx,rn;
-
-#ifdef DEBUG_CPS
-  cout<<"\nIn old_calc_dv_inf"<<endl;
-  cout<<"b2="<<b2<<"\nb4="<<b4<<"\nb6="<<b6<<"\nb8="<<b8<<endl;
-#endif
-  vector<bigcomplex> rt = solvecubic(b2/4,b4/2,b6/4);
-
-  bigfloat i1l,i1u,i2l=to_bigfloat(0),i2u=to_bigfloat(0);  // Bounds on Intervals
-  long i,numint;
-  if (del<0)
-    { for (i=0; i<3; i++)
-	{ if (is_approx_zero(rt[i].imag()))
-	     { rr=rt[i].real(); i=3; }
-	}
-      numint=0; 
-      if (rr<=1.0)
-	{ numint=1;  
-          i1l=rr;  i1u=1.0;
-	  if (rr<-1.0) { i1l=-1.0; }
-	}
-    }
-  else
-    { rn=rt[0].real();
-      rx=rt[1].real();
-	// Want rn<rr<rx
-      if (rn>rx) { rr=rx; rx=rn; rn=rr; }
-      rr=rt[2].real();
-      if (rr<rn)
-        { rr=rn; rn=rt[2].real(); }
-      else if (rx<rr)
-        { rr=rx; rx=rt[2].real(); }
-
-      numint=2;
-      i1l=rn; i1u=rr; i2l=rx; i2u=1.0;
-	// Deal With Right Most Interval
-      if (i2l>i2u) { numint=1; }
-      if (i2l<-1.0) 
-	{ numint=1;
-	  i1l=-1.0; i1u=i2u;
-        }
-	// Now Deal With Left Most Interval
-      if (i1u>1.0)  { i1u=1.0; }
-      if (i1l<-1.0) { i1l=-1.0; }
-      if (i1l>i1u) 
-	{ numint=numint-1;
-	  i1l=-1.0; i1u=i2u;
-	}
-    }
-
-  if(numint==0) return to_bigfloat(-1);  // code for "infinity"
-
-  bigfloat* te = new bigfloat[36];
-  long tec=2*numint;
-  if (numint>0)
-    { te[0]=i1l; te[1]=i1u;
-      if (numint==2)
-        { te[2]=i2l; te[3]=i2u; }
-    }
-
-  if (numint!=0)
-    { // Roots of g
-      rt=solvequartic(to_bigfloat(0),-b4,-2.0*b6,-b8);
-      for (i=0; i<4; i++)
-        { if (is_approx_zero(rt[i].imag()))
-	    { rr=rt[i].real();
-              if (is_in_int(rr,i1l,i1u,i2l,i2u,numint))
-		{ te[tec]=rr; tec=tec+1; }
-            }
-        }
-
-      // Roots of f (again!)
-      rt=solvecubic(b2/4.0,b4/2.0,b6/4.0);
-      for (i=0; i<3; i++)
-        { if (is_approx_zero(rt[i].imag()))
-            { rr=rt[i].real();
-              if (is_in_int(rr,i1l,i1u,i2l,i2u,numint))
-                { te[tec]=rr; tec=tec+1; }
-            }
-        }
-
-      // Roots of f+g
-      rt=solvequartic(to_bigfloat(4),b2-b4,2.0*(b4-b6),b6-b8);
-      for (i=0; i<4; i++)
-        { if (is_approx_zero(rt[i].imag()))
-            { rr=rt[i].real();
-              if (is_in_int(rr,i1l,i1u,i2l,i2u,numint))
-                { te[tec]=rr; tec=tec+1; }
-            }
-        }
-
-      // Roots of g-f // Change from NPS
-      rt=solvequartic(to_bigfloat(-4),-(b2+b4),-2.0*(b4+b6),-(b6+b8));
-      for (i=0; i<4; i++)
-        { if (is_approx_zero(rt[i].imag()))
-            { rr=rt[i].real();
-              if (is_in_int(rr,i1l,i1u,i2l,i2u,numint))
-                { te[tec]=rr; tec=tec+1; }
-            }
-        }
-
-      // Roots of f'
-      if (-24.0*b4+b2*b2>=0.0)
-        { rn=sqrt(-24.0*b4+b2*b2);
-	  rr=(-b2+rn)/12.0;
-	  if (is_in_int(rr,i1l,i1u,i2l,i2u,numint))
-                { te[tec]=rr; tec=tec+1; }
-	  rr=(-b2-rn)/12.0;
-	  if (is_in_int(rr,i1l,i1u,i2l,i2u,numint))
-                { te[tec]=rr; tec=tec+1; }
-	}
-
-      // Roots of g'
-      rt=solvecubic(to_bigfloat(0),-b4/2.0,-b6/2.0);
-      for (i=0; i<3; i++)
-        { if (is_approx_zero(rt[i].imag()))
-            { rr=rt[i].real();
-              if (is_in_int(rr,i1l,i1u,i2l,i2u,numint))
-                { te[tec]=rr; tec=tec+1; }
-            }
-        }
-
-    }
-
-  bigfloat dv=to_bigfloat(0);
-  // Evaluate |f|, |g| at each of the tec(>0) x-values in array te:
-  for (i=0; i<tec; i++)
-    { rr=te[i]*te[i];
-      rn=4.0*rr*te[i]+b2*rr+2.0*b4*te[i]+b6;
-#ifdef DEBUG_CPS
-      cout<<"x="<<te[i]<<", fx="<<rn<<endl;
-#endif
-      rn=abs(rn);
-      rx=rr*rr-b4*rr-2.0*b6*te[i]-b8;
-      rx=abs(rx);
-#ifdef DEBUG_CPS
-      cout<<"x="<<te[i]<<", gx="<<rx<<endl;
-#endif
-      if (rn>rx) { rx=rn; }
-      if (i==0)  { dv=rx; }
-      else if (dv>rx) { dv=rx; }
-#ifdef DEBUG_CPS
-      cout<<"dv so far ="<<dv<<endl;
-#endif
-    }
-  delete[] te;
-
-  return dv;
-}
-
-bigfloat old_calc_dvd_inf(const bigfloat& b2, const bigfloat& b4, const bigfloat& b6, const bigfloat& b8, const bigfloat& del)
-{ 
-  bigfloat rr=to_bigfloat(0),rn;
-
-#ifdef DEBUG_CPS
-  cout<<"\nIn old_calc_dvd_inf"<<endl;
-  cout<<"b2="<<b2<<"\nb4="<<b4<<"\nb6="<<b6<<"\nb8="<<b8<<endl;
-#endif
-  vector<bigcomplex> rt = solvecubic(b2/4,b4/2,b6/4);
-  long i,j,numrt;
-  bigfloat rrt[4];
-  rrt[0]=to_bigfloat(0);
-  if (del<0)
-    { for (i=0; i<3; i++)
-        { if (is_approx_zero(rt[i].imag()))
-             { rr=rt[i].real(); i=3; }
-        }
-      if (is_approx_zero(rr)) { numrt=1; }
-      else                     { rrt[1]=1.0/rr; numrt=2; }
-    }
-  else
-    { numrt=1;
-      for (i=0; i<3; i++)
-	{ if (!is_approx_zero(rt[i].real()))
-	    { rrt[numrt]=1.0/rt[i].real(); 
-	      numrt++;
-            }
-	}
-    }
- 
-  // Sort the real roots 
-  long fl=0;
-  while (fl==0)
-    { fl=1;
-      for (i=0; i<numrt-1; i++)
-        { if (rrt[i]>rrt[i+1])
-	    { rr=rrt[i+1];
-	      rrt[i+1]=rrt[i];
-	      rrt[i]=rr;
-	      fl=0;
-	    }
-        }
-    }
-
-  // Make Intervals
-  bigfloat **intervals=new bigfloat*[3];
-  for (i=0; i<3; i++)
-    { intervals[i]=new bigfloat[2]; }
-
-  long numint=0;
-  if (b6<0.0)
-    { for (i=0; i<numrt; i=i+2)
-	{ intervals[numint][1]=rrt[numrt-i-1];
-	  intervals[numint][0]=rrt[numrt-i-2];
-	  numint=numint+1;
-	}
-    }
-  else
-    { intervals[numint][1]=rrt[numrt-1]+2;
-      intervals[numint][0]=rrt[numrt-1];
-      numint=numint+1;
-      for (i=1; i<numrt; i=i+2)
-        { intervals[numint][1]=rrt[numrt-i-1];
-	  if (numrt-i-2>0) 
-	    { intervals[numint][0]=rrt[numrt-i-2]; }
-	  else
-	    { intervals[numint][0]=-1.0; }
-          numint=numint+1;
-        }
-    }
-  for (i=0; i<numint; i++)
-    { if (intervals[i][1]>1.0)  { intervals[i][1]=1.0; }
-      if (intervals[i][0]<-1.0) { intervals[i][0]=-1.0; }
-    }
-  i=0;
-  while (i<numint)
-    { fl=0;
-      while (fl==0 && i<numint)
-	{ fl=1;
-          if (intervals[i][0]>intervals[i][1])
-             { fl=0; 
-	       for (j=i; j<numint-1; j++)
-		  { intervals[j][0]=intervals[j+1][0];
-		    intervals[j][1]=intervals[j+1][1];
-		  }
-	       numint=numint-1;
-	     }
-	}
-      i=i+1;
-    }
-
-  bigfloat *te=new bigfloat[36];
-  long tec=2*numint;
-  for (i=0; i<numint; i++)
-     { te[2*i]=intervals[i][0];
-       te[2*i+1]=intervals[i][1];
-     }
-
-  if (numint!=0)
-    { // Roots of g
-      rt=solvequartic(to_bigfloat(0),-b4,-2.0*b6,-b8);
-      for (i=0; i<4; i++)
-	{ if (!is_approx_zero(rt[i]))
-	    { rt[i]=to_bigfloat(1)/rt[i];
-              if (is_approx_zero(rt[i].imag()))
-                 { rr=rt[i].real();
-                   if (is_in_int2(rr,intervals,numint))
-                     { te[tec]=rr; tec=tec+1; }
-		 }
-            }
-        }
-
-      // Roots of f
-      rr=0.0;
-      if (is_in_int2(rr,intervals,numint))
-                     { te[tec]=rr; tec=tec+1; }
-      rt=solvecubic(b2/4.0,b4/2.0,b6/4.0);
-      for (i=0; i<3; i++)
-        { if (!is_approx_zero(rt[i]))
-            { rt[i]=to_bigfloat(1)/rt[i];
-              if (is_approx_zero(rt[i].imag()))
-                 { rr=rt[i].real();
-                   if (is_in_int2(rr,intervals,numint))
-                     { te[tec]=rr; tec=tec+1; }
-                 }
-            }
-        }
-
-      // Roots of f+g
-      rt=solvequartic(to_bigfloat(4),b2-b4,2.0*(b4-b6),b6-b8);
-      for (i=0; i<4; i++)
-        { if (!is_approx_zero(rt[i]))
-            { rt[i]=to_bigfloat(1)/rt[i];
-              if (is_approx_zero(rt[i].imag()))
-                 { rr=rt[i].real();
-                   if (is_in_int2(rr,intervals,numint))
-                     { te[tec]=rr; tec=tec+1; }
-		 }
-            }
-        }
-
-      // Roots of f-g
-      rt=solvequartic(to_bigfloat(-4),-b2-b4,2.0*(-b4-b6),-b6-b8);
-      for (i=0; i<4; i++)
-        { if (!is_approx_zero(rt[i]))
-            { rt[i]=to_bigfloat(1)/rt[i];
-	      if (is_approx_zero(rt[i].imag()))
-                 { rr=rt[i].real();
-                   if (is_in_int2(rr,intervals,numint))
-                     { te[tec]=rr; tec=tec+1; }
-		 }
-            }
-        }
-
-      // Roots of f'
-      rt=solvecubic(b2/2.0,3.0*b4/2.0,b6);
-      for (i=0; i<3; i++)
-        { if (!is_approx_zero(rt[i]))
-            { if (is_approx_zero(rt[i].imag()))
-                { rr=rt[i].real();
-                  if (is_in_int2(rr,intervals,numint))
-                    { te[tec]=rr; tec=tec+1; }
-		}
-            }
-        }
-
-      // Roots of g'
-      rr=0.0;
-      if (is_in_int2(rr,intervals,numint))
-          { te[tec]=rr; tec=tec+1; }
-      bigfloat dr=9.0*b6*b6-8.0*b8*b4;
-      if (dr>=0.0)
-	{ dr=sqrt(dr);
-	  rr=(3.0*b6+dr)/(-4.0*b8);
-	  if (is_in_int2(rr,intervals,numint))
-                { te[tec]=rr; tec=tec+1; }
-	  rr=(3.0*b6-dr)/(-4.0*b8);
-          if (is_in_int2(rr,intervals,numint))
-                { te[tec]=rr; tec=tec+1; }
-	}
-    }
-
-  bigfloat dvd=to_bigfloat(0);
-  bigfloat f,g;
-  for (i=0; i<tec; i++)
-    { rr=te[i];
-#ifdef DEBUG_CPS
-      cout<<"x="<<te[i]<<", fx="<<rr<<endl;
-#endif
-      f=(((b6*rr+2.0*b4)*rr+b2)*rr+4.0)*rr;
-      g=(((-b8*rr-2.0*b6)*rr-b4)*rr)*rr+1.0;
-      f=abs(f); g=abs(g);
-      rn=f;
-#ifdef DEBUG_CPS
-      cout<<"x="<<te[i]<<", gx="<<rn<<endl;
-#endif
-      if (g>f) { rn=g; }
-      if (i==0)  { dvd=rn; }
-      else if (dvd>rn) { dvd=rn; }
-#ifdef DEBUG_CPS
-      cout<<"dvd so far ="<<dvd<<endl;
-#endif
-    }
-  delete[] te;
-  for (i=0; i<3; i++)
-    { delete[] intervals[i]; }
-  delete[] intervals;
-
-  return dvd;
-
-}
-
-// coeff has length 5 but may start with leading zeros
-// returns real roots in [-1,1]
-vector<bigfloat> roots11( const vector<bigfloat>& coeff )
-{
-#ifdef DEBUG_CPS
-  cout<<"In roots11() with coeff =  "<<coeff<<endl;
-#endif
-  bigfloat a=coeff[0], a1;
-#ifdef DEBUG_CPS
-  cout<<"coeff[0] = a = "<<a<<endl;
-#endif
-  if(a!=0) 
-    {
-#ifdef DEBUG_CPS
-      cout<<"inverting nonzero a ="<<a<<endl;
-#endif
-      a1=1/a;
-      vector<bigcomplex> cr = solvequartic(a1*coeff[1],a1*coeff[2],a1*coeff[3],a1*coeff[4]);
-      return reals_in_11(cr);
-    }
-  a=coeff[1];
-#ifdef DEBUG_CPS
-  cout<<"coeff[1] = a = "<<a<<endl;
-#endif
-  if(a!=0) 
-    {
-#ifdef DEBUG_CPS
-      cout<<"inverting nonzero a ="<<a<<endl;
-#endif
-      a1=1/a;
-      vector<bigcomplex> cr = solvecubic(a1*coeff[2],a1*coeff[3],a1*coeff[4]);
-      return reals_in_11(cr);
-    }
-  vector<bigfloat> ans; // zero length
-  a=coeff[2];
-#ifdef DEBUG_CPS
-  cout<<"coeff[2] = a = "<<a<<endl;
-#endif
-  if(a==0) {cerr<<"Error in roots: degree<2"<<endl; return ans;}
-  bigfloat b=coeff[3], c=coeff[4], x;
-  bigfloat d=b*b-4*a*c;
-  if(d>=0) 
-    {
-      d=sqrt(d);  a1=1/(2*a);
-      x=a1*(-b+d);
-      if((x<=1)&&(x>=-1)) ans.push_back(x);
-      x=a1*(-b-d);
-      if((x<=1)&&(x>=-1)) ans.push_back(x);
-    }
-  return ans;
-}
-
-vector<bigfloat> reals_in ( vector<bigcomplex>& v)
-{
-  vector<bigfloat> vr;
-  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;
-  for (const auto& vi : v)
-    if(is_real(vi))
-      {
-        bigfloat x = vi.real();
-        if((x<=1)&&(x>=-1))
-          vr.push_back(x);
-      }
-  return vr;
-}
-
-int interval_test(const bigfloat& x, const vector<bigfloat> rts, int debug)
-{
-// rts will have size 1 or 3, and be ordered
-  if(debug) cout<<"Interval test("<<x<<"), rts="<<rts<<endl;
-  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)
-    ans = (x>=rts[0]);
-  else
-    ans =  ((x>=rts[0]) && (x<=rts[1])) || (x>=rts[2]);
-  if(debug) cout<<"\t returns "<<ans<<"\n";
-  return ans;
-}
-
-vector<bigfloat> reals_in_interval ( vector<bigcomplex>& v, const vector<bigfloat> rts)
-{
-  vector<bigfloat> vr;
-  bigfloat x;
-  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;
-}
-
-// Inserts from C in to S the elements which are real and between -1
-// and +1
-void include_real_11(std::set<bigfloat>& S, const vector<bigcomplex>& C)
-{
-  bigfloat x;
-  for (const auto& Ci : C)
-    {
-      if(is_real(Ci))
-	{
-	  x = Ci.real();
-	  if((x<=1)&&(x>=-1)) S.insert(x);
-	}
-    }
-}
-
-inline int in_fund_region(const bigcomplex& z)
-{
-  return (imag(z)>0)&&(abs(z)>0.999)&&(abs(2*real(z))<1.001);
-}
-
 //#define HTB_DEBUG
 
-
 // returns lower bound for height of non-torsion points, following
 // Cremona & Siksek in ANTS7.  If egr==1, return a lower bound for the
 // height of nontorsion pints with everywhere good reduction.
@@ -1101,61 +552,6 @@ bigfloat lower_height_bound(const Curvedata& CD, int egr)
   return lambda;
 }
 
-// This gives a lower bound on non-torsion points; experimental, not
-// very good (and not always positive?!!
-
-bigfloat lower_height_bound_alt(const Curvedata& CD)
-{
-  static bigfloat log2 = log(to_bigfloat(2));
-  bigcomplex w1, w2;
-  Cperiods pers(CD);
-  bigcomplex tau = gettau(pers);
-#ifdef HTB_DEBUG
-  cout<<"\nAfter normalizing, periods are: "<<pers<<endl;
-  if(!in_fund_region(tau)) 
-    cout<<"tau is not in the fundamental region!"<<endl;
-#endif  
-  bigfloat aq = abs(q(tau));
-  pers.getwi(w1,w2);  //NB w1 is the minimal period here, not necessarily the real period!
-  int realcase = abs(imag(w1))<0.0001;
-#ifdef HTB_DEBUG
-  cout<<"q  = "<<q(tau)<<endl;
-  cout<<"aq = "<<aq<<endl;
-  cout<<"w1 = "<< w1 << endl;
-  cout<<"w2 = "<< w2 << endl;
-  cout<<"real case?  ";if(realcase)cout<<"yes";else cout<<"no"; cout<<endl;
-#endif  
-
-  // Either (when w1 is real):
-  // (1/6)log|Delta| - 2log2 - (1/6)log|q| - 4(|q|)/(1-|q|)
-
-  // Or (when w1 is not real):
-  // (1/6)log|Delta| - 2log2 + (1/12)log|q| - 2(1+|q|^2)/(1-|q|)
-
-  bigfloat term1 = log(abs(I2bigfloat(getdiscr(CD))))/6;
-  bigfloat term2 = - 2*log2;
-  bigfloat term3, term4;
-
-  if(realcase) 
-    {
-      term3 = -log(aq)/6;
-      term4 = -4*(aq)/(1-aq);
-    }
-  else
-    {
-      term3 = log(aq)/12;
-      term4 = -2*(1+aq*aq)/(1-aq);
-    }
-  bigfloat bound = term1 + term2 + term3 + term4;
-#ifdef HTB_DEBUG
-  cout<<"term 1 = "<<term1<<endl;
-  cout<<"term 2 = "<<term2<<endl;
-  cout<<"term 3 = "<<term3<<endl;
-  cout<<"term 4 = "<<term4<<endl;
-#endif
-  return bound;
-}
-
 // This gives a lower bound on non-torsion points, by searching, given
 // the regulator of a known subgroup
 
@@ -1639,7 +1035,6 @@ int CurveHeightConst::test_target(const bigfloat& target, long k)
   return canonicalHeightInterval01(target,k).size()==0;
 }
 
-
 vector<Interval01> CurveHeightConst::canonicalHeightInterval01(const bigfloat& target, long k)
 {
   vector<Interval01> solution;
@@ -1760,124 +1155,3 @@ bigfloat CurveHeightConst::psi(const bigfloat& x)
   //  cout<<"z = "<<z<<endl;
   return real(z/get_real_period()); // in [0.5,1]
 }
-
-#if(0) // old code using x-intervals instead of [0,1]-intervals
-
-vector<Interval> CurveHeightConst::canonicalHeightInterval(const bigfloat& target, long k)
-{
-  vector<Interval> solution;
-  solution.push_back(Interval(e3)); // i.e. [e3,+infty]
-#ifdef debugLB
-  cout<<"testing target bound "<<target<<" with k = "<<k<<endl;
-  cout<<"Starting interval: "<<solution<<endl;
-#endif
-  for(int n=1; n<=k; n++)
-    {
-     bigfloat B=Bnmu(n,target);      // = B_n(target) in the paper
-#ifdef debugLB
-     cout<<"n = "<<n<<endl;
-     cout<<"B_"<<n<<"("<<target<<") = "<<B<<endl;
-#endif
-     if(B > 1.0e100) break;
-     vector<Interval> t=solveLEQ(n,B);
-     solution=intersect(solution,t);
-#ifdef debugLB
-     cout<<"intervals from solveLEQ: "<<t<<endl;
-     cout<<"intervals now: "<<solution<<endl;
-#endif
-     if(solution.size()==0) return solution;
-     t=solveGEQ(n,-B);
-     solution=intersect(solution,t);
-#ifdef debugLB
-     cout<<"intervals from solveGEQ: "<<t<<endl;
-     cout<<"intervals now: "<<solution<<endl;
-#endif
-     if(solution.size()==0) return solution;
-   }
- return  solution;
-}
-
-// Returns a list of intervals containing those x(P) for which x(nP) <= B
-
-vector<Interval> CurveHeightConst::solveLEQ(long n, const bigfloat& B)
-{
-  vector<Interval> ans;
-  if(B < e3) return ans;
-  vector<bigfloat> xlist=solveEllNPower(n,B);
-  if(n%2) // n is odd
-    {
-      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]));
-    }
-  return ans;
-}
-
-// Returns a list of intervals containing those x(P) for which x(nP) >= B
-
-vector<Interval> CurveHeightConst::solveGEQ(long n, const bigfloat& B)
-{
-  vector<Interval> ans;
-  if(B <= e3)
-    {
-      ans.push_back(Interval(e3)); // i.e.[e3,infty]
-      return ans;
-    }
-  vector<bigfloat> xlist=solveEllNPower(n,B);
-  int i;
-  if(n%2) // n is odd
-    {
-      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[n-1]));
-    }
-  return ans;
-}
-
-// solves x(nP)=x1 on $E_0$; returns vector of n reals in [0,1)
-// (elliptic logs of solutions)
-
-vector<bigfloat> CurveHeightConst::solveEllNPower01(long n, const bigfloat& x1)
-{
-  vector<bigfloat> ans;
-  if(x1<e3)  return ans; // empty
-  bigfloat z=psi(x1);
-#ifdef debugLB
-  cout<<"psi("<<x1<<") = "<<z<<endl;
-#endif
-  bigfloat oneovern = to_bigfloat(1)/to_bigfloat(n);
-  z*=oneovern;
-  for(int i=0; i<n; i++, z+=oneovern) ans.push_back(z);
-#ifdef debugLB
-  cout<<"solveEllNPower01("<<n<<","<<x1<<") returns "<<ans<<endl;
-#endif
-  return ans;
-}
-
-// solves x(nP)=x1 on $E_0$; returns vector of n reals
-
-vector<bigfloat> CurveHeightConst::solveEllNPower(long n, const bigfloat& x1)
-{
-  vector<bigfloat> ans=solveEllNPower01(n,x1);
-  bigfloat om=real(get_real_period());
-  for(int i=0; i<n; i++) ans[i]=real(ztopoint(om*ans[i])[0]);
-  sort(ans.begin(), ans.end());
-#ifdef debugLB
-  cout<<"solveEllNPower("<<n<<","<<x1<<") returns "<<ans<<endl;
-#endif
-  return ans;
-}
-
-#endif
diff --git a/libsrc/mglobsol.cc b/libsrc/mglobsol.cc
index 5b74859..8b81a02 100644
--- a/libsrc/mglobsol.cc
+++ b/libsrc/mglobsol.cc
@@ -100,7 +100,7 @@ quartic_sieve::quartic_sieve(quartic * gg, int moduli_option, int verb)
   if(!easy) 
     {
       nwprimes=25;
-      wprimes = new long[nwprimes];
+      wprimes.reserve(nwprimes);
       long nwp=0;
       long a8=mod(a,8); long b8=mod(2*b,8); long c8=mod(4*c,8);
       long e8=mod(e,8); long d8=mod(2*d,8);
@@ -112,24 +112,21 @@ quartic_sieve::quartic_sieve(quartic * gg, int moduli_option, int verb)
 	  for(x=1; (x<8) && !two_is_ok; x+=2)
 	    {
 	      tt = mod((t + b8*x*z0),8);
-	      if((tt==0)||(tt==1)||(tt==4)) two_is_ok=1; 
+	      if((tt==0)||(tt==1)||(tt==4)) two_is_ok=1;
 	    }
 	}
-      if(!two_is_ok) wprimes[nwp++]=2;
+      if(!two_is_ok) wprimes.push_back(2);
       primevar pr; pr++; // to start at 3
-      for(;nwp<nwprimes; pr++)
+      while(wprimes.size()<nwprimes)
 	{
-	  p=pr;
-	  if(legendre(a,p)==-1) wprimes[nwp++]=p;
-	}
-      if(verbose) 
-	{
-	  cout<<"w-primes: ";
-	  for(nwp=0; nwp<nwprimes; nwp++) cout<<wprimes[nwp]<<" ";
-	  cout<<endl;
+	  p=pr; pr++;
+	  if(legendre(a,p)==-1) wprimes.push_back(p);
 	}
+      if(verbose)
+        cout<<"w-primes: "<<wprimes<<endl;
+
 // repeat for u-primes:
-      uprimes = new long[nwprimes];
+      uprimes.reserve(nwprimes);
       nwp=0;
       two_is_ok = 0;
       for(z0=0; (z0<4) && !two_is_ok; z0++)
@@ -138,22 +135,18 @@ quartic_sieve::quartic_sieve(quartic * gg, int moduli_option, int verb)
 	  for(x=1; (x<8) && !two_is_ok; x+=2)
 	    {
 	      tt = (t + d8*x*z0)%8;
-	      if((tt==0)||(tt==1)||(tt==4)) two_is_ok=1; 
+	      if((tt==0)||(tt==1)||(tt==4)) two_is_ok=1;
 	    }
 	}
-      if(!two_is_ok) uprimes[nwp++]=2;
+      if(!two_is_ok) uprimes.push_back(2);
       pr.init(); pr++; // to start at 3
-      for(;nwp<nwprimes; pr++)
-	{
-	  p=pr;
-	  if(legendre(e,p)==-1) uprimes[nwp++]=p;
-	}
-      if(verbose) 
+      while(uprimes.size()<nwprimes)
 	{
-	  cout<<"u-primes: ";
-	  for(nwp=0; nwp<nwprimes; nwp++) cout<<uprimes[nwp]<<" ";
-	  cout<<endl;
+	  p=pr; pr++;
+	  if(legendre(e,p)==-1) uprimes.push_back(p);
 	}
+      if(verbose)
+        cout<<"u-primes: "<<uprimes<<endl;
     }
 
 // set up list of auxiliary moduli
@@ -161,61 +154,38 @@ quartic_sieve::quartic_sieve(quartic * gg, int moduli_option, int verb)
 
   switch(moduli_option) {
   case 1:
-    num_aux = 10; 
-    auxs = new long[num_aux];
-    auxs[0]=3;
-    auxs[1]=5;
-    auxs[2]=7;
-    auxs[3]=11;
-    auxs[4]=13;
-    auxs[5]=17;
-    auxs[6]=19;
-    auxs[7]=23;
-    auxs[8]=29;
-    auxs[9]=31;
+    num_aux = 10;
+    auxs = {3, 5, 7, 11, 13, 17, 19, 23, 29, 31};
     break;
 
   case 2:// the following taken from Gebel's scheme
-    num_aux = 3; 
-    auxs = new long[num_aux];
-    auxs[0]=5184;  // = (2^6)*(3^4)   // old: 6624; //  = (2^5)*(3^2)*23
-    auxs[1]=5929;  // = (7^2)*(11^2)  // old: 8075; //  = (5^2)*17*19
-    auxs[2]=4225;  // = (5^2)*(13^2)  // old: 7007; //  = (7^2)*11*13
+    num_aux = 3;
+    auxs = {5184, 5929, 4225};
+    // 5184 = (2^6)*(3^4),  5929 = (7^2)*(11^2),  4225 = (5^2)*(13^2)
     break;
 
   case 3:
   default:
     num_aux = 9;
-    auxs = new long[num_aux];
-    auxs[0]=32;
-    auxs[1]= 9;
-    auxs[2]=25;
-    auxs[3]=49;
-    auxs[4]=11;
-    auxs[5]=13;
-    auxs[6]=17;
-    auxs[7]=19;
-    auxs[8]=23;
+    auxs = {32, 9, 25, 49, 11, 13, 17, 19, 23};
     break;
   }
-  
-  xgood_mod_aux = new int*[num_aux];
-//  x1good_mod_aux = new int*[num_aux];
-  squares = new int*[num_aux];
-  umod = new long[num_aux];
+
+  xgood_mod_aux.resize(num_aux);
+  squares.resize(num_aux);
+  umod.resize(num_aux);
 
   long i,j;
-  for (i = 0; i < num_aux; i++)
+  for ( const auto& aux : auxs)
     {
-      long aux = auxs[i];
       long half_aux = ((aux + 1) / 2);
-      squares[i] = new int[aux];
+      squares[i].resize(aux);
       for (j = 0; j < aux; j++)      squares[i][j]=0;
       for (j = 0; j < half_aux; j++) squares[i][posmod( j*j, aux )]=1;
-      xgood_mod_aux[i] = new int[aux];
+      xgood_mod_aux[i].resize(aux);
     }  // end of aux loop
-  
-  if(verbose>1) 
+
+  if(verbose>1)
     {
       cout << "Finished constructing quartic_sieve, using ";
       switch(moduli_option)
@@ -228,24 +198,6 @@ quartic_sieve::quartic_sieve(quartic * gg, int moduli_option, int verb)
     }
 }
 
-quartic_sieve::~quartic_sieve()
-{
-  if(nwprimes) { delete[] wprimes; delete[] uprimes;}
-
-  if(use_stoll)  // using Stoll search so rest not needed.
-    return; 
-  
-  delete[] auxs;
-  for(long i=0; i<num_aux; i++) 
-    {
-      delete[] xgood_mod_aux[i];
-      delete[] squares[i];
-    }
-  delete[] xgood_mod_aux;
-  delete[] squares;
-  delete[] umod;
-}
-
 //#define DEBUG_RANGES
 
 void sort(bigfloat& x1, bigfloat& x2, bigfloat& x3, bigfloat& x4)
@@ -526,7 +478,7 @@ long quartic_sieve::search_range(int lower, bigfloat lower_bound,
 	  w2 = sqr(BIGINT(w));  w3 = w*w2; w4 = w2*w2;
 	  aw = a; bw = b*w; cw = c*w2; dw = d*w3; ew = e*w4;
 
-	  for ( i=0; i < num_aux; i++)  
+	  for ( i=0; i < num_aux; i++)
 	    umod[i] = posmod(first_u-ustep, auxs[i]);
 
 // set up flag array of residues coprime to w
@@ -546,19 +498,18 @@ long quartic_sieve::search_range(int lower, bigfloat lower_bound,
 	      pcw = posmod(cw, aux);
 	      pdw = posmod(dw, aux);
 	      pew = posmod(ew, aux);
-	  
+
 	      long ddddf= posmod(24*paw , aux);
 	      long dddf = posmod(36*paw + 6*pbw , aux);
 	      long ddf  = posmod(14*paw + 6*pbw + 2*pcw , aux);
 	      long df   = posmod(paw+pbw+pcw+pdw, aux);
 	      long f    = posmod(pew , aux);
-	  
-	      int* flag = xgood_mod_aux[index];
-	      int* sqs = squares[index];
-	      long x=aux;
-	      while(x--)
+
+	      vector<int> flag = xgood_mod_aux[index];
+	      vector<int> sqs = squares[index];
+	      for (long x=0; x<aux; x++)
 		{
-		  *flag++ = sqs[f];
+		  flag[x] = sqs[f];
 		  f    +=    df; if(f    >= aux) f    -= aux;
 		  df   +=   ddf; if(df   >= aux) df   -= aux;
 		  ddf  +=  dddf; if(ddf  >= aux) ddf  -= aux;
@@ -582,7 +533,7 @@ long quartic_sieve::search_range(int lower, bigfloat lower_bound,
 		  if(u_is_ok) 
 		    {
 		      u_is_ok = xgood_mod_aux[i][umodi];
-//		      if(!u_is_ok) 
+//		      if(!u_is_ok)
 //		      cout<<"(u,w)=("<<u<<","<<w<<") failed sieve mod "<<auxs[i]<<endl;
 		    }
 		}
@@ -596,13 +547,13 @@ long quartic_sieve::search_range(int lower, bigfloat lower_bound,
 	      umodw+=ustep; while(umodw>=w) umodw-=w;
 	      u_is_ok = wflag[umodw];  // true if gcd(u,w)=1
 	    }
-	  if(!u_is_ok) continue;	  
+	  if(!u_is_ok) continue;
 
 // Check that u has no impossible prime factors:
 	  for(nwp=0; (nwp<nwprimes) && u_is_ok; nwp++)
 	    u_is_ok = ndivides(uprimes[nwp],u);
 
-	  if(!u_is_ok) continue;	  
+	  if(!u_is_ok) continue;
 
 	  if(!w_vars_set)
 	    {
diff --git a/libsrc/mlocsol.cc b/libsrc/mlocsol.cc
index e7046aa..8c038de 100644
--- a/libsrc/mlocsol.cc
+++ b/libsrc/mlocsol.cc
@@ -280,11 +280,8 @@ int new_qpsoluble(const bigint& a, const bigint& b, const bigint& c, const bigin
 
 int new_zpsol(const bigint& a,const bigint& b,const bigint& c,const bigint& d,const bigint& e, const bigint& p, int verbose)
 {
-  bigint *coeff =new bigint[5];
-  coeff[0]=a;  coeff[1]=b;  coeff[2]=c;  coeff[3]=d;  coeff[4]=e;
-  int res=local_sol(p,coeff,verbose);
-  delete[] coeff;
-  return res;
+  vector<bigint> coeff = {a, b, c, d, e};
+  return local_sol(p,coeff,verbose);
 }
 
 /* Samir's Local Solubility Test for odd p */
@@ -293,11 +290,11 @@ int new_zpsol(const bigint& a,const bigint& b,const bigint& c,const bigint& d,co
 #define Term pair_ZZ_pX_long
 #define Poly ZZ_pX
 
-Factorization fact_c(bigint *c, int verbose=0)
+Factorization fact_c(vector<bigint> c, int verbose=0)
 {
-  Poly f; ZZ_p ci;
-  for (long i=0; i<5; i++) 
-    { ci=to_ZZ_p(c[i]); SetCoeff(f,i,ci); }
+  Poly f;
+  for (long i=0; i<5; i++)
+    SetCoeff(f,i,to_ZZ_p(c[i]));
   if(verbose) cout<<"Factorizing "<<f<<" after making monic: ";
   MakeMonic(f);
   if(verbose) cout<<f<<endl;
@@ -305,7 +302,7 @@ Factorization fact_c(bigint *c, int verbose=0)
   return  CanZass(f, verbose);
 }
 
-int local_sol(const bigint& p, bigint *c, int verbose)
+int local_sol(const bigint& p, vector<bigint> c, int verbose)
 {
   ZZ_p::init(p);
   if (verbose)
@@ -328,13 +325,11 @@ int local_sol(const bigint& p, bigint *c, int verbose)
       if (verbose) cout << "f is 0 mod p: Case II" << endl;
       zeromodp=1;
       for (i=0; (i<5) && zeromodp; i++) { zeromodp=div(p2,c[i]); }
-      bigint *dd=new bigint[5];
+      vector<bigint> dd(5);
       if (zeromodp)
          { for (i=0; i<5; i++)   { dd[i]=c[i]/p2; }
 	   if (verbose) cout << "f is 0 mod p^2, recursing" << endl;
-           int res=local_sol(p,dd,verbose);
-           delete[] dd;
-           return res;
+           return local_sol(p,dd,verbose);
          }
       for (i=0; i<5; i++) { dd[i]=c[i]/p; }
       Factorization fact_f=fact_c(dd);
@@ -342,13 +337,12 @@ int local_sol(const bigint& p, bigint *c, int verbose)
       if (verbose) cout << "Factorization of f/p = "<<fact_f << endl;
       for (i=0; i<fact_f.length(); i++)
         { term=fact_f[i];
-	  F=term.a; 
+	  F=term.a;
 	  e=term.b;
           if ((deg(F)==1) && (e==1))
-	    { if (verbose) 
+	    { if (verbose)
 	      cout << "Non-Repeated Root " << -ConstTerm(F)
 		   <<", returning 1"<<endl;
-              delete[] dd;
               return 1;
 	    }
         }
@@ -356,7 +350,7 @@ int local_sol(const bigint& p, bigint *c, int verbose)
       /* Go thru each repeated root and make the
          required transformation
       */
-      bigint *d=new bigint[5];
+      vector<bigint> d(5);
       fl=0;
       for (i=0; i<fact_f.length() && (!fl); i++)
         { term=fact_f[i];
@@ -364,7 +358,7 @@ int local_sol(const bigint& p, bigint *c, int verbose)
           e=term.b;
           if ((deg(F)==1) && (e!=1))
             { bigint r=-rep(ConstTerm(F));
-              if (verbose) 
+              if (verbose)
 		cout << "Repeated Root=" << r << ", recursing"<<endl;
  // Using f(pX+r)/p^2
               d[4]=dd[4]*p2*p;
@@ -375,8 +369,6 @@ int local_sol(const bigint& p, bigint *c, int verbose)
               fl=local_sol(p,d,verbose);
             }
         }
-      delete[] d;
-      delete [] dd;
       return fl;
     }
 // CASE I
@@ -450,13 +442,13 @@ int local_sol(const bigint& p, bigint *c, int verbose)
   /* For each root which is also a root of h
      transform the equations and call again
   */
-  bigint *d=new bigint[5];
+  vector<bigint> d(5);
   fl=0;
   for (i=0; i<num_r && !fl; i++)
     { te=(((h4*r[i]+h3)*r[i]+h2)*r[i])%p;
       te=((te+h1)*r[i]+h0)%p;
       if (is_zero(te))
-        { if (verbose) 
+        { if (verbose)
         	  cout << "Using " << r[i] << " and recursing"<<endl;
           d[4]=c[4]*p2;
           d[3]=p*(c[3]+4*c[4]*r[i]);
@@ -466,8 +458,6 @@ int local_sol(const bigint& p, bigint *c, int verbose)
           fl=local_sol(p,d,verbose);
         }
     }
-
-  delete[] d;
   return fl;
 }
 
diff --git a/libsrc/quadratic.cc b/libsrc/quadratic.cc
index 0b5166c..ef3fbe2 100644
--- a/libsrc/quadratic.cc
+++ b/libsrc/quadratic.cc
@@ -23,19 +23,6 @@
  
 #include <eclib/quadratic.h>
 
-quadratic::~quadratic() {delete [] coeffs;}
-
-void quadratic::init()
-{
-  coeffs = new bigint[3];
-}
-
-bigint quadratic::coeff(int i) 
-    {if((i>=0)&&(i<=2)) return coeffs[i]; else {bigint ans; return ans;}}
-
-bigint quadratic::operator[](int i) const
-    {if((i>=0)&&(i<=2)) return coeffs[i]; else {bigint ans; return ans;}}
-
 void quadratic::transform(const unimod& m)
 {
   bigint newq0 = eval(m.m11, m.m21);
@@ -93,7 +80,10 @@ void quadratic::reduce(unimod& m)
 }
 
 bigint resultant(const quadratic& q1, const quadratic& q2)
-{return sqr(q2.coeffs[0]*q1.coeffs[2]) + sqr(q1.coeffs[0]*q2.coeffs[2]) +
-   q2.coeffs[2]*q2.coeffs[0]*sqr(q1.coeffs[1]) - q2.coeffs[1]*q1.coeffs[1]*(q2.coeffs[0]*q1.coeffs[2] + q1.coeffs[0]*q2.coeffs[2]) +
-   (sqr(q2.coeffs[1])-2*q2.coeffs[0]*q2.coeffs[2])*q1.coeffs[0]*q1.coeffs[2];
+{
+  const bigint&
+    a1=q1.coeffs[0], b1=q1.coeffs[1], c1=q1.coeffs[2],
+    a2=q2.coeffs[0], b2=q2.coeffs[1], c2=q2.coeffs[2];
+  return
+    sqr(c2*a1) -a1*c2*b2*b1 + (-2*c2*a2 + sqr(b2))*c1*a1 + c2*a2*sqr(b1) - b2*a2*c1*b1 + sqr(a2*c1);
 }