diff --git a/arrays.scad b/arrays.scad
index b49cc85..9ba0850 100644
--- a/arrays.scad
+++ b/arrays.scad
@@ -34,15 +34,19 @@
 //   is_homogeneous( [[1,["a"]], [2,[true]]], 1 ) // Returns true
 //   is_homogeneous( [[1,["a"]], [2,[true]]], 2 ) // Returns false
 //   is_homogeneous( [[1,["a"]], [true,["b"]]] )  // Returns true
-function is_homogeneous(l, depth) =
+function is_homogeneous(l, depth=10) =
     !is_list(l) || l==[] ? false :
     let( l0=l[0] )
     [] == [for(i=[1:len(l)-1]) if( ! _same_type(l[i],l0, depth+1) )  0 ];
                  
 function _same_type(a,b, depth) = 
-    (depth==0) || (a>=b) || (a==b) || (a<=b)  
-    ||  ( is_list(a) && is_list(b) && len(a)==len(b) 
-          && []==[for(i=idx(a)) if( ! _same_type(a[i],b[i],depth-1) )0] ); 
+    (depth==0) ||
+    (is_undef(a) && is_undef(b)) ||
+    (is_bool(a) && is_bool(b)) ||
+    (is_num(a) && is_num(b)) ||
+    (is_string(a) && is_string(b)) ||
+    (is_list(a) && is_list(b) && len(a)==len(b) 
+          && []==[for(i=idx(a)) if( ! _same_type(a[i],b[i],depth-1) ) 0] ); 
   
 
 // Function: select()
diff --git a/common.scad b/common.scad
index 570ff8e..5a22b61 100644
--- a/common.scad
+++ b/common.scad
@@ -90,13 +90,13 @@ function is_nan(x) = (x!=x);
 //   is_finite(x);
 // Description:
 //   Returns true if a given value `x` is a finite number.
-function is_finite(v) = is_num(0*v);
+function is_finite(x) = is_num(x) && !is_nan(0*x);
 
 
 // Function: is_range()
 // Description:
 //   Returns true if its argument is a range
-function is_range(x) = !is_list(x) && is_finite(x[0]+x[1]+x[2]) ;
+function is_range(x) = !is_list(x) && is_finite(x[0]) && is_finite(x[1]) && is_finite(x[2]) ;
 
 
 // Function: valid_range()
@@ -458,5 +458,10 @@ module shape_compare(eps=1/1024) {
 }
 
 
+function loop_start() = 0;
+function loop_done(x) = x==1;
+function looping(x) = x<2;
+function loop_next(x,b) = x>=1? 2 : (b? 0 : 1);
+
 
 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
diff --git a/geometry.scad b/geometry.scad
index cf78195..0832644 100644
--- a/geometry.scad
+++ b/geometry.scad
@@ -177,6 +177,7 @@ function line_ray_intersection(line,ray,eps=EPSILON) =
     let(
         isect = _general_line_intersection(line,ray,eps=eps)
     ) 
+    is_undef(isect[0]) ? undef :
     (isect[2]<0-eps) ? undef : isect[0];
 
 
@@ -195,7 +196,10 @@ function line_segment_intersection(line,segment,eps=EPSILON) =
     assert( _valid_line(line,  dim=2,eps=eps) &&_valid_line(segment,dim=2,eps=eps), "Invalid line or segment." )
     let(
         isect = _general_line_intersection(line,segment,eps=eps)
-    ) isect[2]<0-eps || isect[2]>1+eps ? undef : isect[0];
+    )
+    is_undef(isect[0]) ? undef :
+    isect[2]<0-eps || isect[2]>1+eps ? undef :
+    isect[0];
 
 
 // Function: ray_intersection()
@@ -214,6 +218,7 @@ function ray_intersection(r1,r2,eps=EPSILON) =
     let(
         isect = _general_line_intersection(r1,r2,eps=eps)
     ) 
+    is_undef(isect[0]) ? undef :
     isect[1]<0-eps || isect[2]<0-eps ? undef : isect[0];
 
 
@@ -233,11 +238,9 @@ function ray_segment_intersection(ray,segment,eps=EPSILON) =
     let(
         isect = _general_line_intersection(ray,segment,eps=eps)
     ) 
-    isect[1]<0-eps 
-    || isect[2]<0-eps 
-    || isect[2]>1+eps
-    ? undef 
-    : isect[0];
+    is_undef(isect[0]) ? undef :
+    isect[1]<0-eps || isect[2]<0-eps || isect[2]>1+eps ? undef :
+    isect[0];
 
 
 // Function: segment_intersection()
@@ -256,12 +259,9 @@ function segment_intersection(s1,s2,eps=EPSILON) =
     let(
         isect = _general_line_intersection(s1,s2,eps=eps)
     ) 
-    isect[1]<0-eps 
-    || isect[1]>1+eps 
-    || isect[2]<0-eps 
-    || isect[2]>1+eps 
-    ? undef 
-    : isect[0];
+    is_undef(isect[0]) ? undef :
+    isect[1]<0-eps || isect[1]>1+eps || isect[2]<0-eps || isect[2]>1+eps ? undef :
+    isect[0];
 
 
 // Function: line_closest_point()
diff --git a/math.scad b/math.scad
index 9a0a26e..759fe7e 100644
--- a/math.scad
+++ b/math.scad
@@ -56,7 +56,7 @@ function log2(x) =
 
 // Function: hypot()
 // Usage:
-//   l = hypot(x,y,[z]);
+//   l = hypot(x,y,<z>);
 // Description:
 //   Calculate hypotenuse length of a 2D or 3D triangle.
 // Arguments:
@@ -73,7 +73,7 @@ function hypot(x,y,z=0) =
 
 // Function: factorial()
 // Usage:
-//   x = factorial(n,[d]);
+//   x = factorial(n,<d>);
 // Description:
 //   Returns the factorial of the given integer value, or n!/d! if d is given.  
 // Arguments:
@@ -373,7 +373,7 @@ function modang(x) =
 
 // Function: modrange()
 // Usage:
-//   modrange(x, y, m, [step])
+//   modrange(x, y, m, <step>)
 // Description:
 //   Returns a normalized list of numbers from `x` to `y`, by `step`, modulo `m`.  Wraps if `x` > `y`.
 // Arguments:
@@ -401,7 +401,7 @@ function modrange(x, y, m, step=1) =
 
 // Function: rand_int()
 // Usage:
-//   rand_int(minval,maxval,N,[seed]);
+//   rand_int(minval,maxval,N,<seed>);
 // Description:
 //   Return a list of random integers in the range of minval to maxval, inclusive.
 // Arguments:
@@ -421,7 +421,7 @@ function rand_int(minval, maxval, N, seed=undef) =
 
 // Function: gaussian_rands()
 // Usage:
-//   gaussian_rands(mean, stddev, [N], [seed])
+//   gaussian_rands(mean, stddev, <N>, <seed>)
 // Description:
 //   Returns a random number with a gaussian/normal distribution.
 // Arguments:
@@ -437,7 +437,7 @@ function gaussian_rands(mean, stddev, N=1, seed=undef) =
 
 // Function: log_rands()
 // Usage:
-//   log_rands(minval, maxval, factor, [N], [seed]);
+//   log_rands(minval, maxval, factor, <N>, <seed>);
 // Description:
 //   Returns a single random number, with a logarithmic distribution.
 // Arguments:
@@ -506,6 +506,8 @@ function lcm(a,b=[]) =
 // Section: Sums, Products, Aggregate Functions.
 
 // Function: sum()
+// Usage:
+//   x = sum(v, <dflt>);
 // Description:
 //   Returns the sum of all entries in the given consistent list.
 //   If passed an array of vectors, returns the sum the vectors.
@@ -518,8 +520,7 @@ function lcm(a,b=[]) =
 //   sum([1,2,3]);  // returns 6.
 //   sum([[1,2,3], [3,4,5], [5,6,7]]);  // returns [9, 12, 15]
 function sum(v, dflt=0) =
-    is_list(v) && len(v) == 0 ? dflt :
-    is_vector(v) || is_matrix(v)? [for(i=v) 1]*v :
+    v==[]? dflt :
     assert(is_consistent(v), "Input to sum is non-numeric or inconsistent")
     _sum(v,v[0]*0);
 
@@ -527,6 +528,8 @@ function _sum(v,_total,_i=0) = _i>=len(v) ? _total : _sum(v,_total+v[_i], _i+1);
 
 
 // Function: cumsum()
+// Usage:
+//   sums = cumsum(v);
 // Description:
 //   Returns a list where each item is the cumulative sum of all items up to and including the corresponding entry in the input list.
 //   If passed an array of vectors, returns a list of cumulative vectors sums.
@@ -573,6 +576,8 @@ function sum_of_sines(a, sines) =
 
 
 // Function: deltas()
+// Usage:
+//   delts = deltas(v);
 // Description:
 //   Returns a list with the deltas of adjacent entries in the given list.
 //   The list should be a consistent list of numeric components (numbers, vectors, matrix, etc).
@@ -588,6 +593,8 @@ function deltas(v) =
 
 
 // Function: product()
+// Usage:
+//   x = product(v);
 // Description:
 //   Returns the product of all entries in the given list.
 //   If passed a list of vectors of same dimension, returns a vector of products of each part.
@@ -611,6 +618,8 @@ function _product(v, i=0, _tot) =
 
 
 // Function: outer_product()
+// Usage:
+//   x = outer_product(u,v);
 // Description:
 //   Compute the outer product of two vectors, a matrix.  
 // Usage:
@@ -621,6 +630,8 @@ function outer_product(u,v) =
 
 
 // Function: mean()
+// Usage:
+//   x = mean(v);
 // Description:
 //   Returns the arithmetic mean/average of all entries in the given array.
 //   If passed a list of vectors, returns a vector of the mean of each part.
@@ -660,14 +671,15 @@ function convolve(p,q) =
 // Section: Matrix math
 
 // Function: linear_solve()
-// Usage: linear_solve(A,b)
+// Usage:
+//   solv = linear_solve(A,b)
 // Description:
-//   Solves the linear system Ax=b.  If A is square and non-singular the unique solution is returned.  If A is overdetermined
-//   the least squares solution is returned. If A is underdetermined, the minimal norm solution is returned.
-//   If A is rank deficient or singular then linear_solve returns [].  If b is a matrix that is compatible with A
+//   Solves the linear system Ax=b.  If `A` is square and non-singular the unique solution is returned.  If `A` is overdetermined
+//   the least squares solution is returned. If `A` is underdetermined, the minimal norm solution is returned.
+//   If `A` is rank deficient or singular then linear_solve returns `[]`.  If `b` is a matrix that is compatible with `A`
 //   then the problem is solved for the matrix valued right hand side and a matrix is returned.  Note that if you 
-//   want to solve Ax=b1 and Ax=b2 that you need to form the matrix transpose([b1,b2]) for the right hand side and then
-//   transpose the returned value.  
+//   want to solve Ax=b1 and Ax=b2 that you need to form the matrix `transpose([b1,b2])` for the right hand side and then
+//   transpose the returned value.
 function linear_solve(A,b,pivot=true) =
     assert(is_matrix(A), "Input should be a matrix.")
     let(
@@ -690,33 +702,32 @@ function linear_solve(A,b,pivot=true) =
 
 // Function: matrix_inverse()
 // Usage:
-//    matrix_inverse(A)
+//    mat = matrix_inverse(A)
 // Description:
-//    Compute the matrix inverse of the square matrix A.  If A is singular, returns undef.
-//    Note that if you just want to solve a linear system of equations you should NOT
-//    use this function.  Instead use linear_solve, or use qr_factor.  The computation
+//    Compute the matrix inverse of the square matrix `A`.  If `A` is singular, returns `undef`.
+//    Note that if you just want to solve a linear system of equations you should NOT use this function.
+//    Instead use [[`linear_solve()`|linear_solve]], or use [[`qr_factor()`|qr_factor]].  The computation
 //    will be faster and more accurate.  
 function matrix_inverse(A) =
-    assert(is_matrix(A,square=true),"Input to matrix_inverse() must be a square matrix")
+    assert(is_matrix(A) && len(A)==len(A[0]),"Input to matrix_inverse() must be a square matrix")
     linear_solve(A,ident(len(A)));
 
 
 // Function: null_space()
 // Usage:
-//   null_space(A)
+//   x = null_space(A)
 // Description:
-//   Returns an orthonormal basis for the null space of A, namely the vectors {x} such that Ax=0.  If the null space
-//   is just the origin then returns an empty list. 
+//   Returns an orthonormal basis for the null space of `A`, namely the vectors {x} such that Ax=0.
+//   If the null space is just the origin then returns an empty list. 
 function null_space(A,eps=1e-12) =
     assert(is_matrix(A))
     let(
-      Q_R=qr_factor(transpose(A),pivot=true),
-      R=Q_R[1],
-      zrow = [for(i=idx(R)) if (all_zero(R[i],eps)) i]
+        Q_R = qr_factor(transpose(A),pivot=true),
+        R = Q_R[1],
+        zrow = [for(i=idx(R)) if (all_zero(R[i],eps)) i]
     )
-    len(zrow)==0
-      ? []
-      : transpose(subindex(Q_R[0],zrow));
+    len(zrow)==0 ? [] :
+    transpose(subindex(Q_R[0],zrow));
 
 
 // Function: qr_factor()
@@ -730,19 +741,18 @@ function null_space(A,eps=1e-12) =
 function qr_factor(A, pivot=false) =
     assert(is_matrix(A), "Input must be a matrix." )
     let(
-      m = len(A),
-      n = len(A[0])
+        m = len(A),
+        n = len(A[0])
     )
     let(
-        qr =_qr_factor(A, Q=ident(m),P=ident(n), pivot=pivot, column=0, m = m, n=n),
-        Rzero = 
-          let( R = qr[1] )
-          [ for(i=[0:m-1]) [
-              let( ri = R[i] )
-              for(j=[0:n-1]) i>j ? 0 : ri[j]
+        qr = _qr_factor(A, Q=ident(m),P=ident(n), pivot=pivot, column=0, m = m, n=n),
+        Rzero = let( R = qr[1]) [
+            for(i=[0:m-1]) [
+                let( ri = R[i] )
+                for(j=[0:n-1]) i>j ? 0 : ri[j]
             ]
-          ]
-    ) [qr[0],Rzero,qr[2]];
+        ]
+    ) [qr[0], Rzero, qr[2]];
 
 function _qr_factor(A,Q,P, pivot, column, m, n) =
     column >= min(m-1,n) ? [Q,A,P] :
@@ -770,7 +780,7 @@ function _swap_matrix(n,i,j) =
 
 
 // Function: back_substitute()
-// Usage: back_substitute(R, b, [transpose])
+// Usage: back_substitute(R, b, <transpose>)
 // Description:
 //   Solves the problem Rx=b where R is an upper triangular square matrix.  The lower triangular entries of R are
 //   ignored.  If transpose==true then instead solve transpose(R)*x=b.
@@ -799,6 +809,8 @@ function _back_substitute(R, b, x=[]) =
 
 
 // Function: det2()
+// Usage:
+//   d = det2(M);
 // Description:
 //   Optimized function that returns the determinant for the given 2x2 square matrix.
 // Arguments:
@@ -807,11 +819,13 @@ function _back_substitute(R, b, x=[]) =
 //   M = [ [6,-2], [1,8] ];
 //   det = det2(M);  // Returns: 50
 function det2(M) = 
-    assert( 0*M==[[0,0],[0,0]], "Matrix should be 2x2." )
+    assert(is_matrix(M,2,2), "Matrix must be 2x2.")
     M[0][0] * M[1][1] - M[0][1]*M[1][0];
 
 
 // Function: det3()
+// Usage:
+//   d = det3(M);
 // Description:
 //   Optimized function that returns the determinant for the given 3x3 square matrix.
 // Arguments:
@@ -820,13 +834,15 @@ function det2(M) =
 //   M = [ [6,4,-2], [1,-2,8], [1,5,7] ];
 //   det = det3(M);  // Returns: -334
 function det3(M) =
-    assert( 0*M==[[0,0,0],[0,0,0],[0,0,0]], "Matrix should be 3x3." )
+    assert(is_matrix(M,3,3), "Matrix must be 3x3.")
     M[0][0] * (M[1][1]*M[2][2]-M[2][1]*M[1][2]) -
     M[1][0] * (M[0][1]*M[2][2]-M[2][1]*M[0][2]) +
     M[2][0] * (M[0][1]*M[1][2]-M[1][1]*M[0][2]);
 
 
 // Function: determinant()
+// Usage:
+//   d = determinant(M);
 // Description:
 //   Returns the determinant for the given square matrix.
 // Arguments:
@@ -835,7 +851,7 @@ function det3(M) =
 //   M = [ [6,4,-2,9], [1,-2,8,3], [1,5,7,6], [4,2,5,1] ];
 //   det = determinant(M);  // Returns: 2267
 function determinant(M) =
-    assert(is_matrix(M,square=true), "Input should be a square matrix." )
+    assert(is_matrix(M, square=true), "Input should be a square matrix." )
     len(M)==1? M[0][0] :
     len(M)==2? det2(M) :
     len(M)==3? det3(M) :
@@ -856,23 +872,22 @@ function determinant(M) =
 
 // Function: is_matrix()
 // Usage:
-//   is_matrix(A,[m],[n],[square])
+//   is_matrix(A,<m>,<n>,<square>)
 // Description:
 //   Returns true if A is a numeric matrix of height m and width n.  If m or n
 //   are omitted or set to undef then true is returned for any positive dimension.
-//   If `square` is true then the matrix is required to be square. 
-//   specify m != n and require a square matrix then the result will always be false.
 // Arguments:
-//   A = matrix to test
-//   m = optional height of matrix
-//   n = optional width of matrix
-//   square = set to true to require a square matrix.  Default: false        
+//   A = The matrix to test.
+//   m = Is given, requires the matrix to have the given height.
+//   n = Is given, requires the matrix to have the given width.
+//   square = If true, requires the matrix to have a width equal to its height. Default: false
 function is_matrix(A,m,n,square=false) =
-    is_list(A[0]) 
-    && ( let(v = A*A[0]) is_num(0*(v*v)) ) // a matrix of finite numbers 
-    && (is_undef(n) || len(A[0])==n )
-    && (is_undef(m) || len(A)==m )
-    && ( !square || len(A)==len(A[0]));
+   is_list(A)
+   && (( is_undef(m) && len(A) ) || len(A)==m)
+   && is_list(A[0])
+   && (( is_undef(n) && len(A[0]) ) || len(A[0])==n)
+   && (!square || len(A) == len(A[0]))
+   && is_consistent(A);
 
 
 // Function: norm_fro()
@@ -891,7 +906,7 @@ function norm_fro(A) =
 
 // Function: all_zero()
 // Usage:
-//   all_zero(x);
+//   x = all_zero(x, <eps>);
 // Description:
 //   Returns true if the finite number passed to it is approximately zero, to within `eps`.
 //   If passed a list, recursively checks if all items in the list are approximately zero.
@@ -912,7 +927,7 @@ function all_zero(x, eps=EPSILON) =
 
 // Function: all_nonzero()
 // Usage:
-//   all_nonzero(x);
+//   x = all_nonzero(x, <eps>);
 // Description:
 //   Returns true if the finite number passed to it is not almost zero, to within `eps`.
 //   If passed a list, recursively checks if all items in the list are not almost zero.
@@ -1030,7 +1045,7 @@ function all_nonnegative(x) =
 
 // Function: approx()
 // Usage:
-//   approx(a,b,[eps])
+//   b = approx(a,b,<eps>)
 // Description:
 //   Compares two numbers or vectors, and returns true if they are closer than `eps` to each other.
 // Arguments:
@@ -1044,12 +1059,9 @@ function all_nonnegative(x) =
 //   approx(0.3333,1/3,eps=1e-3);  // Returns: true
 //   approx(PI,3.1415926536);     // Returns: true
 function approx(a,b,eps=EPSILON) = 
-    a==b? true :
-    a*0!=b*0? false :
-    is_list(a)
-    ? ([for (i=idx(a)) if( !approx(a[i],b[i],eps=eps)) 1] == [])
-    : is_num(a) && is_num(b) && (abs(a-b) <= eps);
-    
+    (a==b && is_bool(a) == is_bool(b)) ||
+    (is_num(a) && is_num(b) && abs(a-b) <= eps) ||
+    (is_list(a) && is_list(b) && len(a) == len(b) && [] == [for (i=idx(a)) if (!approx(a[i],b[i],eps=eps)) 1]);
 
 
 function _type_num(x) =
@@ -1063,7 +1075,7 @@ function _type_num(x) =
 
 // Function: compare_vals()
 // Usage:
-//   compare_vals(a, b);
+//   b = compare_vals(a, b);
 // Description:
 //   Compares two values.  Lists are compared recursively.
 //   Returns <0 if a<b.  Returns >0 if a>b.  Returns 0 if a==b.
@@ -1081,7 +1093,7 @@ function compare_vals(a, b) =
 
 // Function: compare_lists()
 // Usage:
-//   compare_lists(a, b)
+//   b = compare_lists(a, b)
 // Description:
 //   Compare contents of two lists using `compare_vals()`.
 //   Returns <0 if `a`<`b`.
@@ -1102,6 +1114,8 @@ function compare_lists(a, b) =
 
 
 // Function: any()
+// Usage:
+//   b = any(l);
 // Description:
 //   Returns true if any item in list `l` evaluates as true.
 //   If `l` is a lists of lists, `any()` is applied recursively to each sublist.
@@ -1115,17 +1129,19 @@ function compare_lists(a, b) =
 //   any([[0,0], [1,0]]);   // Returns true.
 function any(l) =
     assert(is_list(l), "The input is not a list." )
-    _any(l, i=0, succ=false);
+    _any(l);
 
 function _any(l, i=0, succ=false) =
     (i>=len(l) || succ)? succ :
-    _any( l, 
-         i+1, 
-         succ = is_list(l[i]) ? _any(l[i]) : !(!l[i])
-        );
+    _any(
+        l, i+1, 
+        succ = is_list(l[i]) ? _any(l[i]) : !(!l[i])
+    );
 
 
 // Function: all()
+// Usage:
+//   b = all(l);
 // Description:
 //   Returns true if all items in list `l` evaluate as true.
 //   If `l` is a lists of lists, `all()` is applied recursively to each sublist.
@@ -1138,21 +1154,21 @@ function _any(l, i=0, succ=false) =
 //   all([[0,0], [0,0]]);   // Returns false.
 //   all([[0,0], [1,0]]);   // Returns false.
 //   all([[1,1], [1,1]]);   // Returns true.
-function all(l, i=0, fail=false) =
+function all(l) =
     assert( is_list(l), "The input is not a list." )
-    _all(l, i=0, fail=false);
+    _all(l);
 
 function _all(l, i=0, fail=false) =
     (i>=len(l) || fail)? !fail :
-    _all( l, 
-         i+1,
-         fail = is_list(l[i]) ? !_all(l[i]) : !l[i]
-        ) ;
+    _all(
+        l, i+1,
+        fail = is_list(l[i]) ? !_all(l[i]) : !l[i]
+    ) ;
 
 
 // Function: count_true()
 // Usage:
-//   count_true(l)
+//   n = count_true(l)
 // Description:
 //   Returns the number of items in `l` that evaluate as true.
 //   If `l` is a lists of lists, this is applied recursively to each
@@ -1170,26 +1186,22 @@ function _all(l, i=0, fail=false) =
 //   count_true([[0,0], [1,0]]);   // Returns 1.
 //   count_true([[1,1], [1,1]]);   // Returns 4.
 //   count_true([[1,1], [1,1]], nmax=3);  // Returns 3.
+function _count_true_rec(l, nmax, _cnt=0, _i=0) =
+    _i>=len(l) || (is_num(nmax) && _cnt>=nmax)? _cnt :
+    _count_true_rec(l, nmax, _cnt=_cnt+(l[_i]?1:0), _i=_i+1);
+
 function count_true(l, nmax) = 
-    !is_list(l) ? !(!l) ? 1: 0 :
-    let( c = [for( i = 0,
-                   n = !is_list(l[i]) ? !(!l[i]) ? 1: 0 : undef,
-                   c = !is_undef(n)? n : count_true(l[i], nmax),
-                   s = c;
-                 i<len(l) && (is_undef(nmax) || s<nmax);
-                   i = i+1,
-                   n = !is_list(l[i]) ? !(!l[i]) ? 1: 0 : undef,
-                   c = !is_undef(n) || (i==len(l))? n : count_true(l[i], nmax-s),
-                   s = s+c
-                 )  s ] )
-    len(c)<len(l)? nmax: c[len(c)-1];
+    is_undef(nmax)? len([for (x=l) if(x) 1]) :
+    !is_list(l) ? ( l? 1: 0) :
+    _count_true_rec(l, nmax);
 
 
 
 // Section: Calculus
 
 // Function: deriv()
-// Usage: deriv(data, [h], [closed])
+// Usage:
+//   x = deriv(data, [h], [closed])
 // Description:
 //   Computes a numerical derivative estimate of the data, which may be scalar or vector valued.
 //   The `h` parameter gives the step size of your sampling so the derivative can be scaled correctly. 
@@ -1253,7 +1265,8 @@ function _deriv_nonuniform(data, h, closed) =
 
 
 // Function: deriv2()
-// Usage: deriv2(data, [h], [closed])
+// Usage:
+//   x = deriv2(data, <h>, <closed>)
 // Description:
 //   Computes a numerical estimate of the second derivative of the data, which may be scalar or vector valued.
 //   The `h` parameter gives the step size of your sampling so the derivative can be scaled correctly. 
@@ -1296,7 +1309,8 @@ function deriv2(data, h=1, closed=false) =
 
 
 // Function: deriv3()
-// Usage: deriv3(data, [h], [closed])
+// Usage:
+//   x = deriv3(data, <h>, <closed>)
 // Description:
 //   Computes a numerical third derivative estimate of the data, which may be scalar or vector valued.
 //   The `h` parameter gives the step size of your sampling so the derivative can be scaled correctly. 
@@ -1306,6 +1320,10 @@ function deriv2(data, h=1, closed=false) =
 //   f'''(t) = (-f(t-2*h)+2*f(t-h)-2*f(t+h)+f(t+2*h)) / 2h^3.  At the first and second points from the end
 //   the estimates are f'''(t) = (-5*f(t)+18*f(t+h)-24*f(t+2*h)+14*f(t+3*h)-3*f(t+4*h)) / 2h^3 and
 //   f'''(t) = (-3*f(t-h)+10*f(t)-12*f(t+h)+6*f(t+2*h)-f(t+3*h)) / 2h^3.
+// Arguments:
+//   data = the list of the elements to compute the derivative of.
+//   h = the constant parametric sampling of the data.
+//   closed = boolean to indicate if the data set should be wrapped around from the end to the start.
 function deriv3(data, h=1, closed=false) =
     assert( is_consistent(data) , "Input list is not consistent or not numerical.") 
     assert( len(data)>=5, "Input list has less than 5 elements.") 
@@ -1335,17 +1353,27 @@ function deriv3(data, h=1, closed=false) =
 // Section: Complex Numbers
 
 // Function: C_times()
-// Usage: C_times(z1,z2)
+// Usage:
+//   c = C_times(z1,z2)
 // Description:
 //   Multiplies two complex numbers represented by 2D vectors.  
+//   Returns a complex number as a 2D vector [REAL, IMAGINARY].
+// Arguments:
+//   z1 = First complex number, given as a 2D vector [REAL, IMAGINARY]
+//   z2 = Second complex number, given as a 2D vector [REAL, IMAGINARY]
 function C_times(z1,z2) = 
     assert( is_matrix([z1,z2],2,2), "Complex numbers should be represented by 2D vectors" )
     [ z1.x*z2.x - z1.y*z2.y, z1.x*z2.y + z1.y*z2.x ];
 
 // Function: C_div()
-// Usage: C_div(z1,z2)
+// Usage:
+//   x = C_div(z1,z2)
 // Description:
 //   Divides two complex numbers represented by 2D vectors.  
+//   Returns a complex number as a 2D vector [REAL, IMAGINARY].
+// Arguments:
+//   z1 = First complex number, given as a 2D vector [REAL, IMAGINARY]
+//   z2 = Second complex number, given as a 2D vector [REAL, IMAGINARY]
 function C_div(z1,z2) = 
     assert( is_vector(z1,2) && is_vector(z2), "Complex numbers should be represented by 2D vectors." )
     assert( !approx(z2,0), "The divisor `z2` cannot be zero." ) 
@@ -1358,16 +1386,14 @@ function C_div(z1,z2) =
 
 // Function: quadratic_roots()
 // Usage:
-//    roots = quadratic_roots(a,b,c,[real])
+//    roots = quadratic_roots(a,b,c,<real>)
 // Description:
 //    Computes roots of the quadratic equation a*x^2+b*x+c==0, where the
 //    coefficients are real numbers.  If real is true then returns only the
 //    real roots.  Otherwise returns a pair of complex values.  This method
 //    may be more reliable than the general root finder at distinguishing
 //    real roots from complex roots.  
-
-// https://people.csail.mit.edu/bkph/articles/Quadratics.pdf
-
+//    Algorithm from: https://people.csail.mit.edu/bkph/articles/Quadratics.pdf
 function quadratic_roots(a,b,c,real=false) =
   real ? [for(root = quadratic_roots(a,b,c,real=false)) if (root.y==0) root.x]
   :
@@ -1393,7 +1419,7 @@ function quadratic_roots(a,b,c,real=false) =
 
 // Function: polynomial() 
 // Usage:
-//   polynomial(p, z)
+//   x = polynomial(p, z)
 // Description:
 //   Evaluates specified real polynomial, p, at the complex or real input value, z.
 //   The polynomial is specified as p=[a_n, a_{n-1},...,a_1,a_0]
@@ -1409,8 +1435,8 @@ function polynomial(p,z,k,total) =
 
 // Function: poly_mult()
 // Usage:
-//   polymult(p,q)
-//   polymult([p1,p2,p3,...])
+//   x = polymult(p,q)
+//   x = polymult([p1,p2,p3,...])
 // Description:
 //   Given a list of polynomials represented as real coefficient lists, with the highest degree coefficient first, 
 //   computes the coefficient list of the product polynomial.  
@@ -1481,7 +1507,7 @@ function poly_add(p,q) =
 
 // Function: poly_roots()
 // Usage:
-//   poly_roots(p,[tol])
+//   poly_roots(p,<tol>)
 // Description:
 //   Returns all complex roots of the specified real polynomial p.
 //   The polynomial is specified as p=[a_n, a_{n-1},...,a_1,a_0]
@@ -1551,7 +1577,7 @@ function _poly_roots(p, pderiv, s, z, tol, i=0) =
 
 // Function: real_roots()
 // Usage:
-//   real_roots(p, [eps], [tol])
+//   real_roots(p, <eps>, <tol>)
 // Description:
 //   Returns the real roots of the specified real polynomial p.
 //   The polynomial is specified as p=[a_n, a_{n-1},...,a_1,a_0]
diff --git a/paths.scad b/paths.scad
index 3ec12b6..c94d933 100644
--- a/paths.scad
+++ b/paths.scad
@@ -1215,6 +1215,7 @@ function _sum_preserving_round(data, index=0) =
 // Arguments:
 //   path = path to subdivide
 //   N = scalar total number of points desired or with `method="segment"` can be a vector requesting `N[i]-1` points on segment i.
+//   refine = number of points to add each segment.
 //   closed = set to false if the path is open.  Default: True
 //   exact = if true return exactly the requested number of points, possibly sacrificing uniformity.  If false, return uniform point sample that may not match the number of points requested.  Default: True
 //   method = One of `"length"` or `"segment"`.  If `"length"`, adds vertices evenly along the total path length.  If `"segment"`, adds points evenly among the segments.  Default: `"length"`
@@ -1252,7 +1253,11 @@ function subdivide_path(path, N, refine, closed=true, exact=true, method="length
     assert(is_path(path))
     assert(method=="length" || method=="segment")
     assert(num_defined([N,refine]),"Must give exactly one of N and refine")
-    let(N = first_defined([N,len(path)*refine]))
+    let(
+        N = !is_undef(N)? N :
+            !is_undef(refine)? len(path) * refine :
+            undef
+    )
     assert((is_num(N) && N>0) || is_vector(N),"Parameter N to subdivide_path must be postive number or vector")
     let(
         count = len(path) - (closed?0:1), 
diff --git a/regions.scad b/regions.scad
index fb72723..d6832e2 100644
--- a/regions.scad
+++ b/regions.scad
@@ -144,7 +144,7 @@ function region_path_crossings(path, region, closed=true, eps=EPSILON) = sort([
     ) let (
         isect = _general_line_intersection(segs[si], s2, eps=eps)
     ) if (
-        !is_undef(isect) &&
+        !is_undef(isect[0]) &&
         isect[1] >= 0-eps && isect[1] < 1+eps &&
         isect[2] >= 0-eps && isect[2] < 1+eps
     )
diff --git a/shapes2d.scad b/shapes2d.scad
index 4c11efd..c337fd0 100644
--- a/shapes2d.scad
+++ b/shapes2d.scad
@@ -932,7 +932,10 @@ function oval(r, d, realign=false, circum=false, anchor=CENTER, spin=0) =
 function regular_ngon(n=6, r, d, or, od, ir, id, side, rounding=0, realign=false, anchor=CENTER, spin=0) =
     let(
         sc = 1/cos(180/n),
-        r = get_radius(r1=ir*sc, r2=or, r=r, d1=id*sc, d2=od, d=d, dflt=side/2/sin(180/n))
+        ir = is_finite(ir)? ir*sc : undef,
+        id = is_finite(id)? id*sc : undef,
+        side = is_finite(side)? side/2/sin(180/n) : undef,
+        r = get_radius(r1=ir, r2=or, r=r, d1=id, d2=od, d=d, dflt=side)
     )
     assert(!is_undef(r), "regular_ngon(): need to specify one of r, d, or, od, ir, id, side.")
     let(
@@ -970,7 +973,10 @@ function regular_ngon(n=6, r, d, or, od, ir, id, side, rounding=0, realign=false
 
 module regular_ngon(n=6, r, d, or, od, ir, id, side, rounding=0, realign=false, anchor=CENTER, spin=0) {
     sc = 1/cos(180/n);
-    r = get_radius(r1=ir*sc, r2=or, r=r, d1=id*sc, d2=od, d=d, dflt=side/2/sin(180/n));
+    ir = is_finite(ir)? ir*sc : undef;
+    id = is_finite(id)? id*sc : undef;
+    side = is_finite(side)? side/2/sin(180/n) : undef;
+    r = get_radius(r1=ir, r2=or, r=r, d1=id, d2=od, d=d, dflt=side);
     assert(!is_undef(r), "regular_ngon(): need to specify one of r, d, or, od, ir, id, side.");
     path = regular_ngon(n=n, r=r, rounding=rounding, realign=realign);
     inset = opp_ang_to_hyp(rounding, (180-360/n)/2);
diff --git a/strings.scad b/strings.scad
index 5e34fd9..26df2b7 100644
--- a/strings.scad
+++ b/strings.scad
@@ -196,7 +196,7 @@ function str_frac(str,mixed=true,improper=true,signed=true) =
     signed && str[0]=="-" ? -str_frac(substr(str,1),mixed=mixed,improper=improper,signed=false) :
     signed && str[0]=="+" ?  str_frac(substr(str,1),mixed=mixed,improper=improper,signed=false) :
     mixed ? (                      
-        str_find(str," ")>0 || is_undef(str_find(str,"/"))? (
+        !in_list(str_find(str," "), [undef,0]) || is_undef(str_find(str,"/"))? (
             let(whole = str_split(str,[" "]))
             _str_int_recurse(whole[0],10,len(whole[0])-1) + str_frac(whole[1], mixed=false, improper=improper, signed=false)
         ) : str_frac(str,mixed=false, improper=improper)
@@ -292,12 +292,12 @@ function _str_cmp_recurse(str,sindex,pattern,plen,pindex=0,) =
 // Usage:
 //   str_find(str,pattern,[last],[all],[start])
 // Description:
-//   Searches input string `str` for the string `pattern` and returns the index or indices of the matches in `str`.  
-//   By default str_find() returns the index of the first match in `str`.  If `last` is true then it returns the index of the last match.
+//   Searches input string `str` for the string `pattern` and returns the index or indices of the matches in `str`.
+//   By default `str_find()` returns the index of the first match in `str`.  If `last` is true then it returns the index of the last match.
 //   If the pattern is the empty string the first match is at zero and the last match is the last character of the `str`.
 //   If `start` is set then the search begins at index start, working either forward and backward from that position.  If you set `start`
-//   and `last` is true then the search will find the pattern if it begins at index `start`.  If no match exists, returns undef. 
-//   If you set `all` to true then all str_find() returns all of the matches in a list, or an empty list if there are no matches.  
+//   and `last` is true then the search will find the pattern if it begins at index `start`.  If no match exists, returns `undef`.
+//   If you set `all` to true then `str_find()` returns all of the matches in a list, or an empty list if there are no matches.
 // Arguments:
 //   str = String to search.
 //   pattern = string pattern to search for
diff --git a/tests/test_arrays.scad b/tests/test_arrays.scad
index 23ab265..f40b72a 100644
--- a/tests/test_arrays.scad
+++ b/tests/test_arrays.scad
@@ -8,7 +8,7 @@ module test_is_homogeneous(){
     assert(is_homogeneous([[1,["a"]], [2,[true]]])==false);
     assert(is_homogeneous([[1,["a"]], [2,[true]]],1)==true);
     assert(is_homogeneous([[1,["a"]], [2,[true]]],2)==false);
-    assert(is_homogeneous([[1,["a"]], [true,["b"]]])==true);
+    assert(is_homogeneous([[1,["a"]], [true,["b"]]])==false);
 }
 test_is_homogeneous();
 
@@ -134,12 +134,12 @@ test_list_rotate();
 
 
 module test_deduplicate() {
-    assert(deduplicate([8,3,4,4,4,8,2,3,3,8,8]) == [8,3,4,8,2,3,8]);
-    assert(deduplicate(closed=true, [8,3,4,4,4,8,2,3,3,8,8]) == [8,3,4,8,2,3]);
-    assert(deduplicate("Hello") == "Helo");
-    assert(deduplicate([[3,4],[7,1.99],[7,2],[1,4]],eps=0.1) == [[3,4],[7,2],[1,4]]);
-    assert(deduplicate([], closed=true) == []);
-    assert(deduplicate([[1,[1,[undef]]],[1,[1,[undef]]],[1,[2]],[1,[2,[0]]]])==[[1, [1,[undef]]],[1,[2]],[1,[2,[0]]]]);
+    assert_equal(deduplicate([8,3,4,4,4,8,2,3,3,8,8]), [8,3,4,8,2,3,8]);
+    assert_equal(deduplicate(closed=true, [8,3,4,4,4,8,2,3,3,8,8]), [8,3,4,8,2,3]);
+    assert_equal(deduplicate("Hello"), "Helo");
+    assert_equal(deduplicate([[3,4],[7,1.99],[7,2],[1,4]],eps=0.1), [[3,4],[7,2],[1,4]]);
+    assert_equal(deduplicate([], closed=true), []);
+    assert_equal(deduplicate([[1,[1,[undef]]],[1,[1,[undef]]],[1,[2]],[1,[2,[0]]]]), [[1, [1,[undef]]],[1,[2]],[1,[2,[0]]]]);
 }
 test_deduplicate();
 
diff --git a/tests/test_common.scad b/tests/test_common.scad
index d882d63..d51acea 100644
--- a/tests/test_common.scad
+++ b/tests/test_common.scad
@@ -208,10 +208,12 @@ module test_valid_range() {
     assert(valid_range([0:-1:0]));
     assert(valid_range([10:-1:0]));
     assert(valid_range([2.1:-1.1:0.1]));
-    assert(!valid_range([10:1:0]));
-    assert(!valid_range([2.1:1.1:0.1]));
-    assert(!valid_range([0:-1:10]));
-    assert(!valid_range([0.1:-1.1:2.1]));
+    if (version_num() < 20200600) {
+        assert(!valid_range([10:1:0]));
+        assert(!valid_range([2.1:1.1:0.1]));
+        assert(!valid_range([0:-1:10]));
+        assert(!valid_range([0.1:-1.1:2.1]));
+    }
 }
 test_valid_range();
 
diff --git a/tests/test_math.scad b/tests/test_math.scad
index 2f1faaf..9fb7c65 100644
--- a/tests/test_math.scad
+++ b/tests/test_math.scad
@@ -72,8 +72,7 @@ test_constrain();
 
 module test_is_matrix() {
     assert(is_matrix([[2,3,4],[5,6,7],[8,9,10]]));
-    assert(is_matrix([[2,3,4],[5,6,7],[8,9,10]],square=true));
-    assert(is_matrix([[2,3,4],[5,6,7],[8,9,10]],square=false));
+    assert(is_matrix([[2,3],[5,6],[8,9]],3,2));
     assert(is_matrix([[2,3],[5,6],[8,9]],m=3,n=2));
     assert(is_matrix([[2,3,4],[5,6,7]],m=2,n=3));
     assert(is_matrix([[2,3,4],[5,6,7]],2,3));
@@ -82,8 +81,6 @@ module test_is_matrix() {
     assert(is_matrix([[2,3,4],[5,6,7]],n=3));
     assert(!is_matrix([[2,3,4],[5,6,7]],m=4));
     assert(!is_matrix([[2,3,4],[5,6,7]],n=5));
-    assert(!is_matrix([[2,3,4],[5,6,7]],m=2,n=3,square=true));
-    assert(is_matrix([[2,3,4],[5,6,7],[8,9,10]],square=false));
     assert(!is_matrix([[2,3],[5,6],[8,9]],m=2,n=3));
     assert(!is_matrix([[2,3,4],[5,6,7]],m=3,n=2));
     assert(!is_matrix(undef));
@@ -688,10 +685,10 @@ module test_count_true() {
     assert_equal(count_true([1,false,undef]), 1);
     assert_equal(count_true([1,5,false]), 2);
     assert_equal(count_true([1,5,true]), 3);
-    assert_equal(count_true([[0,0], [0,0]]), 0);
-    assert_equal(count_true([[0,0], [1,0]]), 1);
-    assert_equal(count_true([[1,1], [1,1]]), 4);
-    assert_equal(count_true([[1,1], [1,1]], nmax=3), 3);
+    assert_equal(count_true([[0,0], [0,0]]), 2);
+    assert_equal(count_true([[0,0], [1,0]]), 2);
+    assert_equal(count_true([[1,1], [1,1]]), 2);
+    assert_equal(count_true([1,1,1,1,1], nmax=3), 3);
 }
 test_count_true();
 
diff --git a/tests/test_vectors.scad b/tests/test_vectors.scad
index 1aa7a91..ce8dea0 100644
--- a/tests/test_vectors.scad
+++ b/tests/test_vectors.scad
@@ -4,6 +4,8 @@ include <../std.scad>
 module test_is_vector() {
     assert(is_vector([1,2,3]) == true);
     assert(is_vector([[1,2,3]]) == false);
+    assert(is_vector([[1,2,3,4],[5,6,7,8]]) == false);
+    assert(is_vector([[1,2,3,4],[5,6]]) == false);
     assert(is_vector(["foo"]) == false);
     assert(is_vector([]) == false);
     assert(is_vector(1) == false);
diff --git a/vectors.scad b/vectors.scad
index f63fca1..d9123e1 100644
--- a/vectors.scad
+++ b/vectors.scad
@@ -38,7 +38,7 @@
 //   is_vector([1,1,1],all_nonzero=false);  // Returns true
 //   is_vector([],zero=false);              // Returns false
 function is_vector(v, length, zero, all_nonzero=false, eps=EPSILON) =
-    is_list(v) && is_num(0*(v*v))
+    is_list(v) && is_num(v[0]) && is_num(0*(v*v))
     && (is_undef(length) || len(v)==length)
     && (is_undef(zero) || ((norm(v) >= eps) == !zero))
     && (!all_nonzero || all_nonzero(v)) ;
diff --git a/version.scad b/version.scad
index 0e3bf96..447b3e7 100644
--- a/version.scad
+++ b/version.scad
@@ -8,7 +8,7 @@
 //////////////////////////////////////////////////////////////////////
 
 
-BOSL_VERSION = [2,0,436];
+BOSL_VERSION = [2,0,437];
 
 
 // Section: BOSL Library Version Functions
diff --git a/vnf.scad b/vnf.scad
index 085c281..bd5121a 100644
--- a/vnf.scad
+++ b/vnf.scad
@@ -388,17 +388,24 @@ function vnf_volume(vnf) =
 function vnf_centroid(vnf) =
     let(
         verts = vnf[0],
-        val = sum([ for(face=vnf[1], j=[1:1:len(face)-2])
-                        let(
-                            v0  = verts[face[0]],
-                            v1  = verts[face[j]],
-                            v2  = verts[face[j+1]],
-                            vol = cross(v2,v1)*v0
-                        )
-                        [ vol, (v0+v1+v2)*vol ]
-                  ])    
+        vol = sum([
+            for(face=vnf[1], j=[1:1:len(face)-2]) let(
+                v0  = verts[face[0]],
+                v1  = verts[face[j]],
+                v2  = verts[face[j+1]]
+            ) cross(v2,v1)*v0
+        ]),
+        pos = sum([
+            for(face=vnf[1], j=[1:1:len(face)-2]) let(
+                v0  = verts[face[0]],
+                v1  = verts[face[j]],
+                v2  = verts[face[j+1]],
+                vol = cross(v2,v1)*v0
+            )
+            (v0+v1+v2)*vol
+        ])
     )
-    val[1]/val[0]/4;
+    pos/vol/4;
 
 
 function _triangulate_planar_convex_polygons(polys) =