diff --git a/include/igl/copyleft/quadprog.cpp b/include/igl/copyleft/quadprog.cpp
index 2d0ffa74a3..7ee417e818 100644
--- a/include/igl/copyleft/quadprog.cpp
+++ b/include/igl/copyleft/quadprog.cpp
@@ -1,4 +1,3 @@
-//#define TRACE_SOLVER
 #include "quadprog.h"
 #include "../matlab_format.h"
 #include <vector>
diff --git a/include/igl/decimate.cpp b/include/igl/decimate.cpp
index f68887a167..c5c6cf7f6b 100644
--- a/include/igl/decimate.cpp
+++ b/include/igl/decimate.cpp
@@ -186,11 +186,7 @@ IGL_INLINE bool igl::decimate(
       C.row(e) = p;
       costs(e) = cost;
     },
-#ifndef NDEBUG
-    SIZE_MAX
-#else
     10000
-#endif
     );
     for(int e = 0;e<E.rows();e++)
     {
diff --git a/include/igl/dual_contouring.cpp b/include/igl/dual_contouring.cpp
index f62eaeb68f..4a7839cd14 100644
--- a/include/igl/dual_contouring.cpp
+++ b/include/igl/dual_contouring.cpp
@@ -33,12 +33,12 @@ namespace igl
   {
     // Types
     public:
-      typedef Eigen::Matrix<Scalar,1,3> RowVector3S;
-      typedef Eigen::Matrix<Scalar,1,4> RowVector4S;
-      typedef Eigen::Matrix<Scalar,4,4> Matrix4S;
-      typedef Eigen::Matrix<Scalar,3,3> Matrix3S;
-      typedef Eigen::Matrix<Scalar,3,1> Vector3S;
-      typedef std::tuple<int, int, int> KeyTriplet;
+      using RowVector3S = Eigen::Matrix<Scalar,1,3>;
+      using RowVector4S = Eigen::Matrix<Scalar,1,4>;
+      using Matrix4S = Eigen::Matrix<Scalar,4,4>;
+      using Matrix3S = Eigen::Matrix<Scalar,3,3>;
+      using Vector3S = Eigen::Matrix<Scalar,3,1>;
+      using KeyTriplet = std::tuple<int,int,int>;
     // helper function
     public:
     // Working variables
diff --git a/include/igl/min_quad_with_fixed.cpp b/include/igl/min_quad_with_fixed.cpp
index 89d11ca372..715c7b75de 100644
--- a/include/igl/min_quad_with_fixed.cpp
+++ b/include/igl/min_quad_with_fixed.cpp
@@ -593,13 +593,6 @@ IGL_INLINE Eigen::Matrix<Scalar,n,1> igl::min_quad_with_fixed(
   const Eigen::Matrix<Scalar,m,n> & A,
   const Eigen::Matrix<Scalar,m,1> & b)
 {
-  //printf("igl::min_quad_with_fixed<S,n,m,Hpd>\n");
-  //std::cout<<igl::matlab_format(H,"H")<<std::endl;
-  //std::cout<<igl::matlab_format(f,"f")<<std::endl;
-  //std::cout<<igl::matlab_format(k,"k")<<std::endl;
-  //std::cout<<igl::matlab_format(bc,"bc")<<std::endl;
-  //std::cout<<igl::matlab_format(A,"A")<<std::endl;
-  //std::cout<<igl::matlab_format(b,"b")<<std::endl;
   const auto dyn_n = n == Eigen::Dynamic ? H.rows() : n;
   const auto dyn_m = m == Eigen::Dynamic ? A.rows() : m;
   constexpr const int nn = n == Eigen::Dynamic ? Eigen::Dynamic : n+m;
@@ -654,7 +647,6 @@ IGL_INLINE Eigen::Matrix<Scalar,n,1> igl::min_quad_with_fixed(
   const Eigen::Matrix<double,nn,1> bcbc = make_bcbc();
   const Eigen::Matrix<Scalar,nn,1> xx = 
     min_quad_with_fixed<Scalar,nn,false>(HH,ff,kk,bcbc);
-  //std::cout<<igl::matlab_format(xx.head(dyn_n).eval(),"x")<<std::endl;
   return xx.head(dyn_n);
 }
 
@@ -665,11 +657,6 @@ IGL_INLINE Eigen::Matrix<Scalar,n,1> igl::min_quad_with_fixed(
   const Eigen::Array<bool,n,1> & k,
   const Eigen::Matrix<Scalar,n,1> & bc)
 {
-  //printf("igl::min_quad_with_fixed<S,n,Hpd>\n");
-  //std::cout<<igl::matlab_format(H,"H")<<std::endl;
-  //std::cout<<igl::matlab_format(f,"f")<<std::endl;
-  //std::cout<<igl::matlab_format(k,"k")<<std::endl;
-  //std::cout<<igl::matlab_format(bc,"bc")<<std::endl;
   assert(H.isApprox(H.transpose(),1e-7)); 
   assert(H.rows() == H.cols());
   assert(H.rows() == f.size());
@@ -679,7 +666,6 @@ IGL_INLINE Eigen::Matrix<Scalar,n,1> igl::min_quad_with_fixed(
   // Everything fixed
   if(kcount == (Eigen::Dynamic?H.rows():n))
   {
-    //std::cout<<igl::matlab_format(bc,"x")<<std::endl;
     return bc;
   }
   // Nothing fixed
@@ -690,7 +676,6 @@ IGL_INLINE Eigen::Matrix<Scalar,n,1> igl::min_quad_with_fixed(
     typedef typename 
       std::conditional<Hpd,Eigen::LLT<MatrixSn>,Eigen::CompleteOrthogonalDecomposition<MatrixSn>>::type
       Solver;
-    //std::cout<<igl::matlab_format( Solver(H).solve(-f).eval() ,"x")<<std::endl;
     return Solver(H).solve(-f);
   }
   // All-but-one fixed
@@ -707,7 +692,6 @@ IGL_INLINE Eigen::Matrix<Scalar,n,1> igl::min_quad_with_fixed(
     x(u) = -f(u);
     for(int i=0;i<k.size();i++){ if(i!=u){ x(u)-=bc(i)*H(i,u); } } 
     x(u) /= H(u,u);
-    //std::cout<<igl::matlab_format(x,"x")<<std::endl;
     return x;
   }
   // Is there a smart template way to do this?
@@ -807,7 +791,6 @@ IGL_INLINE Eigen::Matrix<Scalar,n,1> igl::min_quad_with_fixed(
   const Eigen::Array<bool,n,1> & k,
   const Eigen::Matrix<Scalar,n,1> & bc)
 {
-  //printf("igl::min_quad_with_fixed<S,n,kcount,Hpd>\n");
   // 0 and n should be handle outside this function
   static_assert(kcount==Eigen::Dynamic || kcount>0                  ,"");
   static_assert(kcount==Eigen::Dynamic || kcount<n                  ,"");
@@ -893,7 +876,6 @@ IGL_INLINE Eigen::Matrix<Scalar,n,1> igl::min_quad_with_fixed(
       }
     }
   }
-  //std::cout<<igl::matlab_format(x,"x")<<std::endl;
   return x;
 }
 
diff --git a/include/igl/per_corner_normals.cpp b/include/igl/per_corner_normals.cpp
index abe9387ce3..0eac79f41f 100644
--- a/include/igl/per_corner_normals.cpp
+++ b/include/igl/per_corner_normals.cpp
@@ -10,6 +10,7 @@
 #include "vertex_triangle_adjacency.h"
 #include "per_face_normals.h"
 #include "PI.h"
+#include "doublearea.h"
 
 template <typename DerivedV, typename DerivedF, typename DerivedCN>
 IGL_INLINE void igl::per_corner_normals(
@@ -102,7 +103,6 @@ IGL_INLINE void igl::per_corner_normals(
   }
 }
 
-#include "doublearea.h"
 
 template <
   typename DerivedV, 
diff --git a/include/igl/polygon_corners.h b/include/igl/polygon_corners.h
index afdcdd6203..6766d1f6db 100644
--- a/include/igl/polygon_corners.h
+++ b/include/igl/polygon_corners.h
@@ -36,8 +36,8 @@ namespace igl
   // and sizes
   //
   // Inputs:
-  //   Q  #Q by k list of k-gon indices (Q(i,j) = -1 means that face i is a
-  //     j-gon)
+  //   Q  #Q by k list of polygon indices (ith row is a k-gon, unless Q(i,j) =
+  //     -1 then it's a j-gon)
   template <
     typename DerivedQ, 
     typename DerivedI,
diff --git a/include/igl/quadprog.cpp b/include/igl/quadprog.cpp
index ac19f75218..717b04661b 100644
--- a/include/igl/quadprog.cpp
+++ b/include/igl/quadprog.cpp
@@ -12,6 +12,23 @@ IGL_INLINE Eigen::Matrix<Scalar,n,1> igl::quadprog(
   const Eigen::Matrix<Scalar,n,1> & lb,
   const Eigen::Matrix<Scalar,n,1> & ub)
 {
+  // Alec 16/2/2021:
+  // igl::quadprog implements a very simple primal active set method. The new
+  // igl::min_quad_with_fixed is very fast for small dense problems so the
+  // iterations of igl::quadprog become very fast. Even if it ends up doing many
+  // more iterations than igl::copyleft::quadprog it would be much faster (in
+  // reality it doesn't do that many more iterations). It's a healthy 10-100x
+  // faster than igl::copyleft::quadprog for specific cases of QPs.
+  //
+  // Unfortunately, that set is limited. igl::quadprog is really good at tiny
+  // box-constrained QPs with a positive definite objective (like the kind that show
+  // up in dual contouring). igl::copyleft::quadprog handles more general problems
+  // (and also starts to beat igl::quadprog when the number of variables gets over
+  // ~20). I tried extending igl::quadprog so that we could use it for
+  // igl::copyleft::progressive_hulls and drop igl::copyleft::quadprog but it was
+  // trickier than I thought. Something like qpmad or the non GPL version of
+  // quadrog++ would be good future PR.
+  //
   typedef Eigen::Matrix<Scalar,n,1> VectorSn;
   typedef Eigen::Array<bool,n,1>    Arraybn;
   assert( (lb.array() < ub.array() ).all() );
@@ -27,7 +44,6 @@ IGL_INLINE Eigen::Matrix<Scalar,n,1> igl::quadprog(
     // Optimize for this common case.
     // Windows needs template arguments spelled out
     x = min_quad_with_fixed<Scalar,n,m>(H,f,k,bc,A,b);
-    //std::cout<<igl::matlab_format(x,"x")<<std::endl;
     // constraint violations 
     VectorSn vl = lb-x;
     VectorSn vu = x-ub;
@@ -55,10 +71,8 @@ IGL_INLINE Eigen::Matrix<Scalar,n,1> igl::quadprog(
       {
         const Scalar sign = bc(i)==ub(i)?1:-1;
         const Scalar lambda_i = sign * (H.row(i)*x+f(i));
-        //printf("  considering k(%d) (λ = %g)\n",i,lambda_i);
         if(lambda_i > worst_lambda)
         {
-          //printf("  removing k(%d) (λ = %g)\n",i,lambda_i);
           best_remove = i;
           worst_lambda = lambda_i;
         }
@@ -69,27 +83,20 @@ IGL_INLINE Eigen::Matrix<Scalar,n,1> igl::quadprog(
     {
       const auto i = best_add;
       assert(!k(i));
-      //add_lower ? printf("  adding lb(%d)\n",i) : printf("  adding lb(%d)\n",i);
       bc(i) = add_lower ? lb(i) : ub(i);
       k(i) = true;
     }else if(best_remove >= 0)
     {
       const auto i = best_remove;
       assert(k(i));
-      //printf("  removing k(%d) (λ = %g)\n",i,worst_lambda);
       k(i) = false;
     }else /*if(best_add < 0 && best_remove < 0)*/
     {
-      std::cout<<igl::matlab_format(x,"x")<<std::endl;
       return x;
     }
   }
   // Should never happen.
   assert(false && "quadprog failed after too many iterations");
-  //std::cout<<igl::eigen_format(H,"H")<<std::endl;
-  //std::cout<<igl::eigen_format(f,"f")<<std::endl;
-  //std::cout<<igl::eigen_format(lb,"lb")<<std::endl;
-  //std::cout<<igl::eigen_format(ub,"ub")<<std::endl;
   return VectorSn::Zero(dyn_n);
 }