diff --git a/libinterp/corefcn/bsxfun.cc b/libinterp/corefcn/bsxfun.cc
index 8e020bd..d39e49c 100644
--- a/libinterp/corefcn/bsxfun.cc
+++ b/libinterp/corefcn/bsxfun.cc
@@ -360,6 +360,10 @@ dimensionality as the other array.\n\
           || (func.is_function_handle ()
               && ! A.is_object () && ! B.is_object ()))
         {
+          // Fast path.
+          // Use pre-built code that knows the return type of func in advance,
+          // with a separate template expansion for each combinations of types
+
           // This may break if the default behavior is overridden.  But if you
           // override arithmetic operators for builtin classes, you should
           // expect mayhem anyway (constant folding etc).  Querying
@@ -393,7 +397,8 @@ dimensionality as the other array.\n\
           for (octave_idx_type i = 0; i < nd; i++)
             if (dva (i) != dvb (i) && dva (i) != 1 && dvb (i) != 1)
               {
-                error ("bsxfun: dimensions of A and B must match");
+                error ("bsxfun: nonconformant dimensions: %s and %s",
+		       dva.str ().c_str(), dvb.str ().c_str());
                 break;
               }
 
@@ -403,11 +408,11 @@ dimensionality as the other array.\n\
               dim_vector dvc;
               dvc.resize (nd);
 
+	      // dvc (i) is max (dva (i), dvb (i)), or 0 if either is 0
+	      // Use the fact that the dva(i) and dvb(i) are either equal or 1
+	      // (the above error throws if that isn't the case, so we don't reach here.)
               for (octave_idx_type i = 0; i < nd; i++)
-                dvc (i) = (dva (i) < 1 ? dva (i)
-                                       : (dvb (i) < 1 ? dvb (i)
-                                                      : (dva (i) > dvb (i)
-                                                        ? dva (i) : dvb (i))));
+		dvc (i) = (dva (i) == 1 ? dvb (i) : dva (i));
 
               if (dva == dvb || dva.numel () == 1 || dvb.numel () == 1)
                 {
@@ -416,8 +421,8 @@ dimensionality as the other array.\n\
                   inputs (1) = B;
                   retval = func.do_multi_index_op (1, inputs);
                 }
-              else if (dvc.numel () < 1)
-                {
+              else if (dvc.numel () < 1)        // no output elements.
+                {                               // call func to get type right
                   octave_value_list inputs;
                   inputs (0) = A.resize (dvc);
                   inputs (1) = B.resize (dvc);
@@ -425,10 +430,26 @@ dimensionality as the other array.\n\
                 }
               else
                 {
+                  // Slow path.
+                  // Called for custom functions, or non-optimised combinations.
+                  // Cycle through the columns, expanding the size of each as
+                  // necessary Keep track of the type of the return value so
+                  // far.  If the type has been the same on each iteration,
+		  // it is stored in a variable of a dedicated type, and new
+		  // columns are 'insert'ed, for speed.  Otherwise, the value
+		  // is copied to  C,  a generic variable, and new columns are
+		  // appended using do_cat_op (the "slow-slow" path).
+		  //
+		  // Changes double->float and real->complex are explicitly
+		  // converted to the new type of dedicated variable.  Others,
+		  // such as those involving integers, switch to using  C
+		  // and never switch back.
                   octave_idx_type ncount = 1;
                   for (octave_idx_type i = 1; i < nd; i++)
                     ncount *= dvc (i);
 
+		  // Macro to create dedicated variables of each type, and
+		  // a flag that we are using this.  At most one flag is true.
 #define BSXDEF(T) \
                   T result_ ## T; \
                   bool have_ ## T = false;
@@ -464,6 +485,7 @@ dimensionality as the other array.\n\
                       if (maybe_update_column (Bc, B, dvb, dvc, i, idxB))
                         inputs (1) = Bc;
 
+		      // list, but only tmp(0) is ever used.
                       octave_value_list tmp = func.do_multi_index_op (1,
                                                                       inputs);
 
@@ -479,8 +501,12 @@ dimensionality as the other array.\n\
                             result_ ## T .resize (dvc); \
                         }
 
+		      // First pass. Choose which dedicated variable to use.
                       if (i == 0)
                         {
+			  // For floating point types, class_name() doesn't
+			  // tell us the complete type -- just single/double.
+			  // Check for sparsity and complex explicitly.
                           if (! tmp(0).is_sparse_type ())
                             {
                               std::string result_type = tmp(0).class_name ();
@@ -517,6 +543,8 @@ dimensionality as the other array.\n\
                                       result_ComplexNDArray.resize (dvc);
                                     }
                                 }
+			      // For integer types, class_name() tells us
+			      // the complete type, so use the above macro
                               else if BSXINIT(boolNDArray, "logical", bool)
                               else if BSXINIT(int8NDArray, "int8", int8)
                               else if BSXINIT(int16NDArray, "int16", int16)
@@ -526,24 +554,25 @@ dimensionality as the other array.\n\
                               else if BSXINIT(uint16NDArray, "uint16", uint16)
                               else if BSXINIT(uint32NDArray, "uint32", uint32)
                               else if BSXINIT(uint64NDArray, "uint64", uint64)
-                              else
+                              else  // unknown type: use slow-slow path
                                 {
                                   C = tmp (0);
                                   C = C.resize (dvc);
                                 }
                             }
-			  else
+			  else	    // sparse:  use slow-slow path
 			    {
 			      C = tmp (0);
 			      C = C.resize (dvc);
 			    }
                         }
-                      else
-                        {
+                      else	// Iteration other than the first.
+                        {	// Check if the type has changed.
+				// FIXME  This promotes single->double,
+				//        instead of double-> single.
                           update_index (ra_idx, dvc, i);
 
-                          if (have_FloatNDArray
-                              || have_FloatComplexNDArray)
+                          if (have_FloatNDArray || have_FloatComplexNDArray)
                             {
                               if (! tmp(0).is_float_type ())
                                 {
@@ -571,15 +600,12 @@ dimensionality as the other array.\n\
                                       have_FloatComplexNDArray = false;
                                       have_ComplexNDArray = true;
                                     }
-                                  else
-                                    {
-                                      result_NDArray =
-                                        NDArray (result_FloatNDArray);
-                                      result_NDArray.insert
-                                        (tmp(0).array_value (), ra_idx);
-                                      have_FloatNDArray = false;
-                                      have_NDArray = true;
-                                    }
+				  else
+				    {
+                                      result_FloatComplexNDArray.insert
+                                        (tmp(0).float_complex_array_value (),
+					 ra_idx);
+				    }
                                 }
                               else if (tmp(0).is_real_type ())
                                 result_FloatNDArray.insert
@@ -617,6 +643,9 @@ dimensionality as the other array.\n\
                                 }
                             }
 
+			  // Check if type of this column matches the previous.
+			  // If not, go to slow-slow path, using variable  C.
+			  // Otherwise, insert the column into dedicated var.
 #define BSXLOOP(T, CLS, EXTRACTOR) \
                         (have_ ## T) \
                           { \
@@ -651,6 +680,8 @@ dimensionality as the other array.\n\
                   (have_ ## T) \
                     retval(0) = result_ ## T;
 
+		  // Check if the result is still in a dedicated variable.
+		  // If so, return it.  Otherwise, return slow-slow variable, C
                   if BSXEND(NDArray)
                   else if BSXEND(ComplexNDArray)
                   else if BSXEND(FloatNDArray)
