diff --git a/arrays.scad b/arrays.scad
index e21ccfb..5068d3b 100644
--- a/arrays.scad
+++ b/arrays.scad
@@ -18,6 +18,20 @@
 
 // Section: List Query Operations
 
+// Function: is_simple_list()
+// Description:
+//   Returns true just when all elements of `list` are simple values.
+// Usage:
+//   is_simple_list(list)
+// Arguments:
+//   list = The list to check.
+// Example:
+//   a = is_simple_list([3,4,5,6,7,8,9]);  Returns: true
+//   b = is_simple_list([3,4,5,[6],7,8]);  Returns: false
+function is_simple_list(list) =
+		is_list(list)
+		&& []==[for(e=list) if(is_list(e)) 0];
+
 
 // Function: select()
 // Description:
@@ -44,20 +58,22 @@
 //   select(l, [1:3]);  // Returns [4,5,6]
 //   select(l, [1,3]);  // Returns [4,6]
 function select(list, start, end=undef) =
+    assert( is_list(list) || is_string(list), "Invalid list.")
     let(l=len(list))
-    end==undef? (
-        is_num(start)?
-            let(s=(start%l+l)%l) list[s] :
-            assert(is_list(start) || is_range(start), "Invalid start parameter")
-            [for (i=start) list[(i%l+l)%l]]
-    ) : (
-        assert(is_num(start), "Invalid start parameter.")
-        assert(is_num(end), "Invalid end parameter.")
-        let(s=(start%l+l)%l, e=(end%l+l)%l)
-        (s<=e)?
-            [for (i = [s:1:e]) list[i]] :
-            concat([for (i = [s:1:l-1]) list[i]], [for (i = [0:1:e]) list[i]])
-    );
+    l==0 ? []
+    :   end==undef? 
+            is_num(start)?
+                list[ (start%l+l)%l ]
+            :   assert( is_list(start) || is_range(start), "Invalid start parameter")
+                [for (i=start) list[ (i%l+l)%l ] ]
+        :   assert(is_num(start), "Invalid start parameter.")
+            assert(is_num(end), "Invalid end parameter.")
+            let( s = (start%l+l)%l, e = (end%l+l)%l )
+            (s <= e)? [for (i = [s:1:e]) list[i]]
+            :   concat([for (i = [s:1:l-1]) list[i]], [for (i = [0:1:e]) list[i]]) ;
+
+
+
 
 
 // Function: slice()
@@ -65,8 +81,8 @@ function select(list, start, end=undef) =
 //   Returns a slice of a list.  The first item is index 0.
 //   Negative indexes are counted back from the end.  The last item is -1.
 // Arguments:
-//   arr = The array/list to get the slice of.
-//   st = The index of the first item to return.
+//   list = The array/list to get the slice of.
+//   start = The index of the first item to return.
 //   end = The index after the last item to return, unless negative, in which case the last item to return.
 // Example:
 //   slice([3,4,5,6,7,8,9], 3, 5);   // Returns [6,7]
@@ -74,28 +90,35 @@ function select(list, start, end=undef) =
 //   slice([3,4,5,6,7,8,9], 1, 1);   // Returns []
 //   slice([3,4,5,6,7,8,9], 6, -1);  // Returns [9]
 //   slice([3,4,5,6,7,8,9], 2, -2);  // Returns [5,6,7,8]
-function slice(arr,st,end) = let(
-        l=len(arr),
-        s=st<0?(l+st):st,
-        e=end<0?(l+end+1):end
-    ) [for (i=[s:1:e-1]) if (e>s) arr[i]];
+function slice(list,start,end) =
+    assert( is_list(list), "Invalid list" )
+    assert( is_finite(start) && is_finite(end), "Invalid number(s)" )
+    let( l = len(list) )
+    l==0 ? []
+    :   let(
+            s = start<0? (l+start): start,
+            e = end<0? (l+end+1): end
+        ) [for (i=[s:1:e-1]) if (e>s) list[i]];
+
+
 
 
 // Function: in_list()
-// Description: Returns true if value `x` is in list `l`.
+// Description: Returns true if value `val` is in list `list`. When `val==NAN` the answer will be false for any list.
 // Arguments:
-//   x = The value to search for.
-//   l = The list to search.
-//   idx = If given, searches the given subindexes for matches for `x`.
+//   val = The simple value to search for.
+//   list = The list to search.
+//   idx = If given, searches the given subindexes for matches for `val`.
 // Example:
 //   in_list("bar", ["foo", "bar", "baz"]);  // Returns true.
 //   in_list("bee", ["foo", "bar", "baz"]);  // Returns false.
 //   in_list("bar", [[2,"foo"], [4,"bar"], [3,"baz"]], idx=1);  // Returns true.
 function in_list(val,list,idx=undef) = 
     let( s = search([val], list, num_returns_per_match=1, index_col_num=idx)[0] )
-    s==[] ? false
+    s==[] || s[0]==[] ? false
     : is_undef(idx) ? val==list[s] 
     : val==list[s][idx];
+    
 
 
 // Function: min_index()
@@ -108,9 +131,10 @@ function in_list(val,list,idx=undef) =
 //   vals = vector of values
 //   all = set to true to return indices of all occurences of the minimum.  Default: false
 // Example:
-//   min_index([5,3,9,6,2,7,8,2,1]); // Returns: 4
-//   min_index([5,3,9,6,2,7,8,2,1],all=true); // Returns: [4,7]
+//   min_index([5,3,9,6,2,7,8,2,1]); // Returns: 8
+//   min_index([5,3,9,6,2,7,8,2,7],all=true); // Returns: [4,7]
 function min_index(vals, all=false) =
+    assert( is_vector(vals) && len(vals)>0 , "Invalid or empty list of numbers.")
     all ? search(min(vals),vals,0) : search(min(vals), vals)[0];
 
 
@@ -127,6 +151,7 @@ function min_index(vals, all=false) =
 //   max_index([5,3,9,6,2,7,8,9,1]); // Returns: 2
 //   max_index([5,3,9,6,2,7,8,9,1],all=true); // Returns: [2,7]
 function max_index(vals, all=false) =
+    assert( is_vector(vals) && len(vals)>0 , "Invalid or empty list of numbers.")
     all ? search(max(vals),vals,0) : search(max(vals), vals)[0];
 
 
@@ -179,10 +204,12 @@ function list_decreasing(list) =
 //   repeat([1,2,3],3);   // Returns [[1,2,3], [1,2,3], [1,2,3]]
 function repeat(val, n, i=0) =
     is_num(n)? [for(j=[1:1:n]) val] :
+    assert( is_list(n), "Invalid count number.")
     (i>=len(n))? val :
     [for (j=[1:1:n[i]]) repeat(val, n, i+1)];
 
 
+
 // Function: list_range()
 // Usage:
 //   list_range(n, [s], [e])
@@ -192,7 +219,7 @@ function repeat(val, n, i=0) =
 // Description:
 //   Returns a list, counting up from starting value `s`, by `step` increments,
 //   until either `n` values are in the list, or it reaches the end value `e`.
-//   If both `n` and `e` are given, returns `n` values evenly spread fron `s`
+//   If both `n` and `e` are given, returns `n` values evenly spread from `s`
 //   to `e`, and `step` is ignored.
 // Arguments:
 //   n = Desired number of values in returned list, if given.
@@ -205,25 +232,26 @@ function repeat(val, n, i=0) =
 //   list_range(n=4, s=3, step=3);   // Returns [3,6,9,12]
 //   list_range(n=5, s=0, e=10);     // Returns [0, 2.5, 5, 7.5, 10]
 //   list_range(e=3);                // Returns [0,1,2,3]
-//   list_range(e=6, step=2);        // Returns [0,2,4,6]
+//   list_range(e=7, step=2);        // Returns [0,2,4,6]
 //   list_range(s=3, e=5);           // Returns [3,4,5]
 //   list_range(s=3, e=8, step=2);   // Returns [3,5,7]
-//   list_range(s=4, e=8, step=2);   // Returns [4,6,8]
+//   list_range(s=4, e=8.3, step=2); // Returns [4,6,8]
 //   list_range(n=4, s=[3,4], step=[2,3]);  // Returns [[3,4], [5,7], [7,10], [9,13]]
 function list_range(n=undef, s=0, e=undef, step=undef) =
-    (n!=undef && e!=undef)? (
-        assert(is_undef(n) || is_undef(e) || is_undef(step), "At most 2 of n, e, and step can be given.")
-        [for (i=[0:1:n-1]) s+(e-s)*i/(n-1)]
-    ) : let(step = default(step,1))
-    (n!=undef)? [for (i=[0:1:n-1]) let(v=s+step*i) v] :
-    (e!=undef)? [for (v=[s:step:e]) v] :
-    assert(e!=undef||n!=undef, "Must supply one of `n` or `e`.");
+    assert( is_undef(n) || is_finite(n), "Parameter `n` must be a number.")
+    assert( is_undef(n) || is_undef(e) || is_undef(step), "At most 2 of n, e, and step can be given.")
+    let( step = (n!=undef && e!=undef)? (e-s)/(n-1) : default(step,1) )
+    is_undef(e) ? 
+        assert( is_consistent([s, step]), "Incompatible data.")
+        [for (i=[0:1:n-1]) s+step*i ]
+    :   assert( is_vector([s,step,e]), "Start `s`, step `step` and end `e` must be numbers.")
+        [for (v=[s:step:e]) v] ;
+    
 
 
 
 // Section: List Manipulation
 
-
 // Function: reverse()
 // Description: Reverses a list/array.
 // Arguments:
@@ -255,34 +283,38 @@ function reverse(list) =
 //   l8 = list_rotate([1,2,3,4,5],5);  // Returns: [1,2,3,4,5]
 //   l9 = list_rotate([1,2,3,4,5],6);  // Returns: [2,3,4,5,1]
 function list_rotate(list,n=1) =
-    assert(is_list(list)||is_string(list))
-    assert(is_num(n))
+    assert(is_list(list)||is_string(list), "Invalid list or string.")
+    assert(is_finite(n), "Invalid number")
     select(list,n,n+len(list)-1);
 
 
 // Function: deduplicate()
 // Usage:
-//   deduplicate(list);
+//   deduplicate(list,[close],[eps]);
 // Description:
 //   Removes consecutive duplicate items in a list.
+//   When `eps` is zero, the comparison between consecutive items is exact.
+//   Otherwise, when all list items and subitems are numbers, the comparison is within the tolerance `eps`.
 //   This is different from `unique()` in that the list is *not* sorted.
 // Arguments:
 //   list = The list to deduplicate.
 //   closed = If true, drops trailing items if they match the first list item.
-//   eps = The maximum difference to allow between numbers or vectors.
+//   eps = The maximum tolerance between items.
 // Examples:
 //   deduplicate([8,3,4,4,4,8,2,3,3,8,8]);  // Returns: [8,3,4,8,2,3,8]
 //   deduplicate(closed=true, [8,3,4,4,4,8,2,3,3,8,8]);  // Returns: [8,3,4,8,2,3]
 //   deduplicate("Hello");  // Returns: ["H","e","l","o"]
 //   deduplicate([[3,4],[7,2],[7,1.99],[1,4]],eps=0.1);  // Returns: [[3,4],[7,2],[1,4]]
+//   deduplicate([[7,undef],[7,undef],[1,4],[1,4+1e-12]],eps=0);    // Returns: [[7,undef],[1,4],[1,4+1e-12]]
 function deduplicate(list, closed=false, eps=EPSILON) =
     assert(is_list(list)||is_string(list))
-    let(
-        l = len(list),
-        end = l-(closed?0:1)
-    ) (is_num(list[0]) || is_vector(list[0]))?
-        [for (i=[0:1:l-1]) if (i==end || !approx(list[i], list[(i+1)%l], eps)) list[i]] :
-        [for (i=[0:1:l-1]) if (i==end || list[i] != list[(i+1)%l]) list[i]];
+    let( l = len(list),
+         end = l-(closed?0:1) )
+    is_string(list) || (eps==0)
+    ? [for (i=[0:1:l-1]) if (i==end || list[i] != list[(i+1)%l]) list[i]]
+    : [for (i=[0:1:l-1]) if (i==end || !approx(list[i], list[(i+1)%l], eps)) list[i]];
+
+
 
 
 // Function: deduplicate_indexed()
@@ -299,26 +331,28 @@ function deduplicate(list, closed=false, eps=EPSILON) =
 // Examples:
 //   deduplicate_indexed([8,6,4,6,3], [1,4,3,1,2,2,0,1]);  // Returns: [1,4,3,2,0,1]
 //   deduplicate_indexed([8,6,4,6,3], [1,4,3,1,2,2,0,1], closed=true);  // Returns: [1,4,3,2,0]
+//   deduplicate_indexed([[7,undef],[7,undef],[1,4],[1,4],[1,4+1e-12]],eps=0);    // Returns: [0,2,4]
 function deduplicate_indexed(list, indices, closed=false, eps=EPSILON) =
-    assert(is_list(list)||is_string(list))
-    assert(indices==[] || is_vector(indices))
+    assert(is_list(list)||is_string(list), "Improper list or string.")
     indices==[]? [] :
-    let(
-        l = len(indices),
-        end = l-(closed?0:1)
-    ) [
-        for (i = [0:1:l-1]) let(
-            a = list[indices[i]],
-            b = list[indices[(i+1)%l]],
-            eq = (a == b)? true :
-                (a*0 != b*0)? false :
-                is_num(a)? approx(a, b, eps=eps) :
-                is_vector(a)? approx(a, b, eps=eps) :
-                false
-        ) if (i==end || !eq) indices[i]
+    assert(is_vector(indices), "Indices must be a list of numbers.")
+    let( l = len(indices),
+         end = l-(closed?0:1) ) 
+    [ for (i = [0:1:l-1]) 
+        let(
+           a = list[indices[i]],
+           b = list[indices[(i+1)%l]],
+           eq = (a == b)? true :
+                (a*0 != b*0) || (eps==0)? false :
+                is_num(a) || is_vector(a) ? approx(a, b, eps=eps) 
+                : false
+        ) 
+        if (i==end || !eq) indices[i]
     ];
 
 
+
+
 // Function: repeat_entries()
 // Usage:
 //   newlist = repeat_entries(list, N)
@@ -345,17 +379,19 @@ function deduplicate_indexed(list, indices, closed=false, eps=EPSILON) =
 //   echo(repeat_entries(list, 6, exact=false));  // Ouputs [0,0,1,1,2,2,3,3]
 //   echo(repeat_entries(list, [1,1,2,1], exact=false));  // Ouputs [0,1,2,2,3]
 function repeat_entries(list, N, exact = true) =
-    assert(is_list(list))
-    assert((is_num(N) && N>0) || is_vector(N),"Parameter N to repeat_entries must be postive number or vector")
+    assert(is_list(list) && len(list)>0, "The list cannot be void.")
+    assert((is_finite(N) && N>0) || is_vector(N,len(list)),
+            "Parameter N must be a number greater than zero or vector with the same length of `list`")
     let(
         length = len(list),
-        reps_guess = is_list(N)?
-            assert(len(N)==len(list), "Vector parameter N to repeat_entries has the wrong length")
-            N : repeat(N/length,length),
-        reps = exact? _sum_preserving_round(reps_guess) :
-            [for (val=reps_guess) round(val)]
+        reps_guess = is_list(N)? N : repeat(N/length,length),
+        reps = exact ?
+                 _sum_preserving_round(reps_guess) 
+               : [for (val=reps_guess) round(val)]
     )
     [for(i=[0:length-1]) each repeat(list[i],reps[i])];
+    
+
 
 
 // Function: list_set()
@@ -363,12 +399,11 @@ function repeat_entries(list, N, exact = true) =
 //   list_set(list, indices, values, [dflt], [minlen])
 // Description:
 //   Takes the input list and returns a new list such that `list[indices[i]] = values[i]` for all of
-//   the (index,value) pairs supplied.  If you supply `indices` that are beyond the length of the list
-//   then the list is extended and filled in with the `dflt` value.  If you set `minlen` then the list is
-//   lengthed, if necessary, by padding with `dflt` to that length.  The `indices` list can be in any
-//   order but run time will be (much) faster for long lists if it is already sorted.  Reptitions are
-//   not allowed.  If `indices` is given as a non-list scalar, then that index of the given `list` will
-//   be set to the value of `values`.
+//   the (index,value) pairs supplied and unchanged for other indices.  If you supply `indices` that are 
+//   beyond the length of the list then the list is extended and filled in with the `dflt` value.  
+//   If you set `minlen` then the list is lengthed, if necessary, by padding with `dflt` to that length.  
+//   Repetitions in `indices` are not allowed. The lists `indices` and `values` must have the same length.  
+//   If `indices` is given as a scalar, then that index of the given `list` will be set to the scalar value of `values`.
 // Arguments:
 //   list = List to set items in.  Default: []
 //   indices = List of indices into `list` to set.
@@ -378,92 +413,82 @@ function repeat_entries(list, N, exact = true) =
 // Examples:
 //   list_set([2,3,4,5], 2, 21);  // Returns: [2,3,21,5]
 //   list_set([2,3,4,5], [1,3], [81,47]);  // Returns: [2,81,4,47]
-function list_set(list=[],indices,values,dflt=0,minlen=0) =
+function list_set(list=[],indices,values,dflt=0,minlen=0) = 
     assert(is_list(list)||is_string(list))
     !is_list(indices)? (
-        (is_num(indices) && indices<len(list))? [for (i=idx(list)) i==indices? values : list[i]] :
-        list_set(list,[indices],[values],dflt)
-    ) :
-    assert(len(indices)==len(values),"Index list and value list must have the same length")
-    let(
-        sortind = list_increasing(indices) ? list_range(len(indices)) : sortidx(indices),
-        lastind = len(indices)==0 ? -1 : indices[select(sortind,-1)]
-    )
-    concat(
-        [for(j=[0:1:indices[sortind[0]]-1]) j>=len(list) ? dflt : list[j]],
-        [values[sortind[0]]], 
-        [for(i=[1:1:len(sortind)-1]) each
-            assert(indices[sortind[i]]!=indices[sortind[i-1]],"Repeated index")
-            concat(
-                [for(j=[1+indices[sortind[i-1]]:1:indices[sortind[i]]-1]) j>=len(list) ? dflt : list[j]],
-                [values[sortind[i]]]
-            )
-        ],
-        slice(list,1+lastind, len(list)),
-        repeat(dflt, minlen-lastind-1)
-    );
+        (is_finite(indices) && indices<len(list))? 
+           [for (i=idx(list)) i==indices? values : list[i]]
+        :  list_set(list,[indices],[values],dflt) )
+    :   assert(is_vector(indices) && is_list(values) && len(values)==len(indices) ,
+               "Index list and value list must have the same length")
+        let( midx = max(len(list)-1, max(indices)) )
+        [ for(i=[0:midx] ) 
+            let( j = search(i,indices,0),
+                 k = j[0] )
+            assert( len(j)<2, "Repeated indices are not acceptable." )
+            k!=undef ? values[k] :
+            i<len(list) ? list[i]:
+                dflt ,
+          each repeat(dflt, minlen-max(indices))
+        ];
+      
 
 
 // Function: list_insert()
 // Usage:
-//   list_insert(list, pos, elements);
+//   list_insert(list, indices, values);
 // Description:
-//   Insert `elements` into `list` before position `pos`.
+//   Insert `values` into `list` before position `indices`.
 // Example:
 //   list_insert([3,6,9,12],1,5);  // Returns [3,5,6,9,12]
 //   list_insert([3,6,9,12],[1,3],[5,11]);  // Returns [3,5,6,9,11,12]
-function list_insert(list, pos, elements, _i=0) =
+function list_insert(list, indices, values, _i=0) = 
     assert(is_list(list)||is_string(list))
-    is_list(pos)? (
-        assert(len(pos)==len(elements))
-        let(
-            idxs = sortidx(pos),
-            lastidx = pos[idxs[len(idxs)-1]]
-        )
-        concat(
-            [
-                for(i=idx(idxs)) each concat(
-                    assert(pos[idxs[i]]<=len(list), "Indices in pos must be <= len(list)")
-                    [for (j=[(i==0?0:pos[idxs[i-1]]):1:pos[idxs[i]]-1]) list[j]],
-                    [elements[idxs[i]]]
-                )
-            ],
-            [for (j=[lastidx:1:len(list)-1]) list[j]]
-        )
-    ) : (
-        assert(pos<=len(list), "Indices in pos must be <= len(list)")
-        concat(
-            slice(list,0,pos),
-            elements,
-            (pos<len(list)? slice(list,pos,-1) : [])
-        )
-    );
+    ! is_list(indices)? 
+        assert( is_finite(indices) && is_finite(values), "Invalid indices/values." ) 
+        assert( indices<=len(list), "Indices must be <= len(list) ." )
+        [for (i=idx(list)) each ( i==indices?  [ values, list[i] ] : [ list[i] ] ) ]
+    :   assert( is_vector(indices) && is_list(values) && len(values)==len(indices) ,
+               "Index list and value list must have the same length")
+        assert( max(indices)<=len(list), "Indices must be <= len(list) ." )
+        let( maxidx = max(indices),
+             minidx  = min(indices) )
+        [ for(i=[0:1:minidx-1] ) list[i],
+          for(i=[minidx: min(maxidx, len(list)-1)] ) 
+            let( j = search(i,indices,0),
+                 k = j[0],
+                 x = assert( len(j)<2, "Repeated indices are not acceptable." )
+              )
+            each ( k != undef  ? [ values[k], list[i] ] : [ list[i] ] ),
+          for(i=[min(maxidx, len(list)-1)+1:1:len(list)-1] ) list[i],
+          if(maxidx==len(list)) values[max_index(indices)]
+        ];
+
+
 
 
 // Function: list_remove()
 // Usage:
-//   list_remove(list, elements)
+//   list_remove(list, indices)
 // Description:
-//   Remove all items from `list` whose indexes are in `elements`.
+//   Remove all items from `list` whose indexes are in `indices`.
 // Arguments:
 //   list = The list to remove items from.
-//   elements = The list of indexes of items to remove.
+//   indices = The list of indexes of items to remove.
 // Example:
 //   list_insert([3,6,9,12],1);      // Returns: [3,9,12]
 //   list_insert([3,6,9,12],[1,3]);  // Returns: [3,9]
-function list_remove(list, elements) =
-    assert(is_list(list)||is_string(list))
-    !is_list(elements) ? list_remove(list,[elements]) :
-    len(elements)==0 ? list :
-    let(
-        sortind = list_increasing(elements) ? list_range(len(elements)) : sortidx(elements),
-        lastind = elements[select(sortind,-1)]
-    )
-    assert(lastind<len(list),"Element index beyond list end")
-    concat(slice(list, 0, elements[sortind[0]]),
-        [for(i=[1:1:len(sortind)-1]) each slice(list,1+elements[sortind[i-1]], elements[sortind[i]])],
-        slice(list,1+lastind, len(list))
-    );
+function list_remove(list, indices) =
+    assert(is_list(list)||is_string(list), "Invalid list/string." )
+    is_finite(indices) 
+    ?   [  for(i=[0:1:min(indices, len(list)-1)-1]) list[i],
+           for(i=[min(indices, len(list)-1)+1:1:len(list)-1]) list[i]  ]
+    :   assert( is_vector(indices), "Invalid list `indices`." )
+        len(indices)==0 ? list :
+        [ for(i=[0:len(list)-1])
+            if ( []==search(i,indices,1) ) list[i] ]; 
+
+
 
 
 // Function: list_remove_values()
@@ -503,8 +528,8 @@ function list_remove_values(list,values=[],all=false) =
 // Example:
 //   bselect([3,4,5,6,7], [false,true,true,false,true]);  // Returns: [4,5,7]
 function bselect(array,index) =
-    assert(is_list(array)||is_string(array))
-    assert(is_list(index))
+    assert(is_list(array)||is_string(array), "Improper array." )
+    assert(is_list(index) && len(index)>=len(array) , "Improper index list." )
     [for(i=[0:len(array)-1]) if (index[i]) array[i]];
 
 
@@ -514,8 +539,8 @@ function bselect(array,index) =
 // Description:
 //   Opposite of `bselect()`.  Returns a list the same length as `indexlist`, where each item will
 //   either be 0 if the corresponding item in `indexset` is false, or the next sequential value
-//   from `valuelist` if true.  The number of `true` values in `indexset` must be equal to the length
-//   of `valuelist`.
+//   from `valuelist` if the item is true.  The number of `true` values in `indexset` must be equal 
+//   to the length of `valuelist`.
 // Arguments:
 //   indexset = A list of boolean values.
 //   valuelist = The list of values to set into the returned list.
@@ -524,73 +549,79 @@ function bselect(array,index) =
 //   list_bset([false,true,false,true,false], [3,4]);  // Returns: [0,3,0,4,0]
 //   list_bset([false,true,false,true,false], [3,4],dflt=1);  // Returns: [1,3,1,4,1]
 function list_bset(indexset, valuelist, dflt=0) =
-    assert(is_list(indexset))
-    assert(is_list(valuelist))
-    let(
-        trueind = search([true], indexset,0)[0]
-    ) concat(
+    assert(is_list(indexset), "The index set is not a list." )
+    assert(is_list(valuelist), "The `valuelist` is not a list." )
+    let( trueind = search([true], indexset,0)[0] )
+    assert( !(len(trueind)>len(valuelist)), str("List `valuelist` too short; its length should be ",len(trueind)) )
+    assert( !(len(trueind)<len(valuelist)), str("List `valuelist` too long; its length should be ",len(trueind)) )
+    concat(
         list_set([],trueind, valuelist, dflt=dflt),    // Fill in all of the values
         repeat(dflt,len(indexset)-max(trueind)-1)  // Add trailing values so length matches indexset
     );
 
 
 
+
 // Section: List Length Manipulation
 
 // Function: list_shortest()
 // Description:
 //   Returns the length of the shortest sublist in a list of lists.
 // Arguments:
-//   vecs = A list of lists.
-function list_shortest(vecs) =
-    assert(is_list(vecs)||is_string(list))
-    min([for (v = vecs) len(v)]);
+//   array = A list of lists.
+function list_shortest(array) =
+    assert(is_list(array)||is_string(list), "Invalid input." )
+    min([for (v = array) len(v)]);
+
 
 
 // Function: list_longest()
 // Description:
 //   Returns the length of the longest sublist in a list of lists.
 // Arguments:
-//   vecs = A list of lists.
-function list_longest(vecs) =
-    assert(is_list(vecs)||is_string(list))
-    max([for (v = vecs) len(v)]);
+//   array = A list of lists.
+function list_longest(array) =
+    assert(is_list(array)||is_string(list), "Invalid input." )
+    max([for (v = array) len(v)]);
 
 
 // Function: list_pad()
 // Description:
-//   If the list `v` is shorter than `minlen` length, pad it to length with the value given in `fill`.
+//   If the list `array` is shorter than `minlen` length, pad it to length with the value given in `fill`.
 // Arguments:
-//   v = A list.
+//   array = A list.
 //   minlen = The minimum length to pad the list to.
 //   fill = The value to pad the list with.
-function list_pad(v, minlen, fill=undef) =
-    assert(is_list(v)||is_string(list))
-    concat(v,repeat(fill,minlen-len(v)));
+function list_pad(array, minlen, fill=undef) =
+    assert(is_list(array)||is_string(list), "Invalid input." )
+    concat(array,repeat(fill,minlen-len(array)));
 
 
 // Function: list_trim()
 // Description:
-//   If the list `v` is longer than `maxlen` length, truncates it to be `maxlen` items long.
+//   If the list `array` is longer than `maxlen` length, truncates it to be `maxlen` items long.
 // Arguments:
-//   v = A list.
+//   array = A list.
 //   minlen = The minimum length to pad the list to.
-function list_trim(v, maxlen) =
-    assert(is_list(v)||is_string(list))
-    [for (i=[0:1:min(len(v),maxlen)-1]) v[i]];
+function list_trim(array, maxlen) =
+    assert(is_list(array)||is_string(list), "Invalid input." )
+    [for (i=[0:1:min(len(array),maxlen)-1]) array[i]];
 
 
 // Function: list_fit()
 // Description:
-//   If the list `v` is longer than `length` items long, truncates it to be exactly `length` items long.
-//   If the list `v` is shorter than `length` items long, pad it to length with the value given in `fill`.
+//   If the list `array` is longer than `length` items long, truncates it to be exactly `length` items long.
+//   If the list `array` is shorter than `length` items long, pad it to length with the value given in `fill`.
 // Arguments:
-//   v = A list.
+//   array = A list.
 //   minlen = The minimum length to pad the list to.
 //   fill = The value to pad the list with.
-function list_fit(v, length, fill) =
-    assert(is_list(v)||is_string(list))
-    let(l=len(v)) (l==length)? v : (l>length)? list_trim(v,length) : list_pad(v,length,fill);
+function list_fit(array, length, fill) =
+    assert(is_list(array)||is_string(list), "Invalid input." )
+    let(l=len(array)) 
+    l==length ? array : 
+    l> length ? list_trim(array,length) 
+              : list_pad(array,length,fill);
 
 
 
@@ -600,132 +631,131 @@ function list_fit(v, length, fill) =
 // Description:
 //   Shuffles the input list into random order.
 function shuffle(list) =
-    assert(is_list(list)||is_string(list))
+    assert(is_list(list)||is_string(list), "Invalid input." )
     len(list)<=1 ? list :
     let (
         rval = rands(0,1,len(list)),
         left  = [for (i=[0:len(list)-1]) if (rval[i]< 0.5) list[i]],
         right = [for (i=[0:len(list)-1]) if (rval[i]>=0.5) list[i]]
-    ) concat(shuffle(left), shuffle(right));
+    ) 
+    concat(shuffle(left), shuffle(right));
 
 
 // Sort a vector of scalar values
 function _sort_scalars(arr) =
-    len(arr)<=1 ? arr : let(
+    len(arr)<=1 ? arr : 
+    let(
         pivot   = arr[floor(len(arr)/2)],
         lesser  = [ for (y = arr) if (y  < pivot) y ],
         equal   = [ for (y = arr) if (y == pivot) y ],
         greater = [ for (y = arr) if (y  > pivot) y ]
-    ) concat( _sort_scalars(lesser), equal, _sort_scalars(greater) );
+    ) 
+    concat( _sort_scalars(lesser), equal, _sort_scalars(greater) );
 
 
 // Sort a vector of vectors based on the first entry only of each vector
 function _sort_vectors1(arr) =
     len(arr)<=1 ? arr :
-    !(len(arr)>0) ? [] : let(
+    !(len(arr)>0) ? [] : 
+    let(
         pivot   = arr[floor(len(arr)/2)],
         lesser  = [ for (y = arr) if (y[0]  < pivot[0]) y ],
         equal   = [ for (y = arr) if (y[0] == pivot[0]) y ],
         greater = [ for (y = arr) if (y[0]  > pivot[0]) y ]
-    ) concat( _sort_vectors1(lesser), equal, _sort_vectors1(greater) );
+    ) 
+    concat( _sort_vectors1(lesser), equal, _sort_vectors1(greater) );
 
 
 // Sort a vector of vectors based on the first two entries of each vector
 // Lexicographic order, remaining entries of vector ignored
 function _sort_vectors2(arr) =
     len(arr)<=1 ? arr :
-    !(len(arr)>0) ? [] : let(
+    !(len(arr)>0) ? [] : 
+    let(
         pivot   = arr[floor(len(arr)/2)],
         lesser  = [ for (y = arr) if (y[0] < pivot[0] || (y[0]==pivot[0] && y[1]<pivot[1])) y ],
         equal   = [ for (y = arr) if (y[0] == pivot[0] && y[1]==pivot[1]) y ],
-        greater  = [ for (y = arr) if (y[0] > pivot[0] || (y[0]==pivot[0] && y[1]>pivot[1])) y ]
-    ) concat( _sort_vectors2(lesser), equal, _sort_vectors2(greater) );
+        greater = [ for (y = arr) if (y[0] > pivot[0] || (y[0]==pivot[0] && y[1]>pivot[1])) y ]
+    ) 
+    concat( _sort_vectors2(lesser), equal, _sort_vectors2(greater) );
 
 // Sort a vector of vectors based on the first three entries of each vector
 // Lexicographic order, remaining entries of vector ignored
 function _sort_vectors3(arr) =
     len(arr)<=1 ? arr : let(
         pivot   = arr[floor(len(arr)/2)],
-        lesser  = [
-            for (y = arr) if (
-                y[0] < pivot[0] || (
-                    y[0]==pivot[0] && (
-                        y[1]<pivot[1] || (
-                            y[1]==pivot[1] &&
-                            y[2]<pivot[2]
-                        )
-                    )
-                )
-            ) y
-        ],
-        equal = [
-            for (y = arr) if (
-                y[0] == pivot[0] && y[1]== pivot[1] && y[2]==pivot[2]
-            ) y
-        ],
-        greater = [
-            for (y = arr) if (
-                y[0] > pivot[0] || (
-                    y[0]==pivot[0] && (
-                        y[1]>pivot[1] || (
-                            y[1]==pivot[1] &&
-                            y[2]>pivot[2]
-                        )
-                    )
-                )
-            ) y
-        ]
+        lesser  = [ for (y = arr) 
+                      if ( y[0] < pivot[0] 
+                           || ( y[0]==pivot[0] 
+                                && ( y[1]<pivot[1] 
+                                     || ( y[1]==pivot[1] 
+                                          && y[2]<pivot[2] ))))
+                    y ],
+        equal = [ for (y = arr) 
+                    if ( y[0] == pivot[0] 
+                         && y[1]== pivot[1] 
+                         && y[2]==pivot[2] )
+                  y ],
+        greater = [ for (y = arr) 
+                      if ( y[0] > pivot[0] 
+                           || ( y[0]==pivot[0] 
+                                && ( y[1] > pivot[1] 
+                                     || ( y[1]==pivot[1] 
+                                          && y[2] > pivot[2] )))) 
+                    y ]
     ) concat( _sort_vectors3(lesser), equal, _sort_vectors3(greater) );
 
 
+
 // Sort a vector of vectors based on the first four entries of each vector
 // Lexicographic order, remaining entries of vector ignored
 function _sort_vectors4(arr) =
     len(arr)<=1 ? arr : let(
         pivot = arr[floor(len(arr)/2)],
-        lesser = [
-            for (y = arr) if (
-                y[0] < pivot[0] || (
-                    y[0]==pivot[0] && (
-                        y[1]<pivot[1] || (
-                            y[1]==pivot[1] && (
-                                y[2]<pivot[2] || (
-                                    y[2]==pivot[2] &&
-                                    y[3]<pivot[3]
-                                )
-                            )
-                        )
-                    )
-                )
-            ) y
-        ],
-        equal = [
-            for (y = arr) if (
-                y[0] == pivot[0] &&
-                y[1] == pivot[1] &&
-                y[2] == pivot[2] &&
-                y[3] == pivot[3]
-            ) y
-        ],
-        greater = [
-            for (y = arr) if (
-                y[0] > pivot[0] || (
-                    y[0]==pivot[0] && (
-                        y[1]>pivot[1] || (
-                            y[1]==pivot[1] && (
-                                y[2]>pivot[2] || (
-                                    y[2]==pivot[2] &&
-                                    y[3]>pivot[3]
-                                )
-                            )
-                        )
-                    )
-                )
-            ) y
-        ]
+        lesser = [  for (y = arr) 
+                      if ( y[0] < pivot[0] 
+                           || ( y[0]==pivot[0] 
+                                && ( y[1]<pivot[1] 
+                                     || ( y[1]==pivot[1] 
+                                         && ( y[2]<pivot[2] 
+                                              || ( y[2]==pivot[2] 
+                                                   && y[3]<pivot[3] ))))))
+                      y ],
+        equal = [ for (y = arr) 
+                    if (  y[0] == pivot[0] 
+                          && y[1] == pivot[1] 
+                          && y[2] == pivot[2] 
+                          && y[3] == pivot[3] ) 
+                  y  ],
+        greater = [ for (y = arr) 
+                      if ( y[0] > pivot[0] 
+                           || ( y[0]==pivot[0] 
+                                && ( y[1]>pivot[1] 
+                                     || ( y[1]==pivot[1] 
+                                          && ( y[2]>pivot[2] 
+                                               || ( y[2]==pivot[2] 
+                                                    && y[3]>pivot[3] )))))) 
+                    y ]
     ) concat( _sort_vectors4(lesser), equal, _sort_vectors4(greater) );
+    
 
 
+function _sort_general(arr, idx=undef) =
+    (len(arr)<=1) ? arr :
+    let(
+        pivot = arr[floor(len(arr)/2)],
+        pivotval = idx==undef? pivot : [for (i=idx) pivot[i]],
+        compare = 
+            is_undef(idx) ? [for(entry=arr) compare_vals(entry, pivotval) ] :
+            [ for (entry = arr) 
+                  let( val = [for (i=idx) entry[i] ] )
+                  compare_vals(val, pivotval) ] ,
+        lesser  = [ for (i = [0:1:len(arr)-1]) if (compare[i] < 0) arr[i] ],
+        equal   = [ for (i = [0:1:len(arr)-1]) if (compare[i] ==0) arr[i] ],
+        greater = [ for (i = [0:1:len(arr)-1]) if (compare[i] > 0) arr[i] ]
+    )
+    concat(_sort_general(lesser,idx), equal, _sort_general(greater,idx));
+    
 function _sort_general(arr, idx=undef) =
     (len(arr)<=1) ? arr :
     let(
@@ -744,6 +774,9 @@ function _sort_general(arr, idx=undef) =
     concat(_sort_general(lesser,idx), equal, _sort_general(greater,idx));
 
 
+
+
+
 // Function: sort()
 // Usage:
 //   sort(list, [idx])
@@ -761,16 +794,20 @@ function _sort_general(arr, idx=undef) =
 //   sorted = sort(l);  // Returns [2,3,8,9,12,16,23,34,37,45,89]
 function sort(list, idx=undef) =
     !is_list(list) || len(list)<=1 ? list :
+    assert( is_undef(idx) || is_finite(idx) || is_vector(idx) || is_range(idx) , "Invalid indices.")
     is_def(idx) ? _sort_general(list,idx) :
     let(size = array_dim(list))
     len(size)==1 ? _sort_scalars(list) :
-    len(size)==2 && size[1] <=4 ? (
+    len(size)==2 && size[1] <=4 
+    ? (
         size[1]==0 ? list :
         size[1]==1 ? _sort_vectors1(list) :
         size[1]==2 ? _sort_vectors2(list) :
-        size[1]==3 ? _sort_vectors3(list) :
-        /*size[1]==4*/ _sort_vectors4(list)
-    ) : _sort_general(list);
+        size[1]==3 ? _sort_vectors3(list)
+   /*size[1]==4*/  : _sort_vectors4(list)
+      ) 
+    : _sort_general(list);
+
 
 
 // Function: sortidx()
@@ -794,6 +831,28 @@ function sort(list, idx=undef) =
 //   idxs1 = sortidx(lst, idx=1); // Returns: [3,0,2,1]
 //   idxs2 = sortidx(lst, idx=0); // Returns: [1,2,0,3]
 //   idxs3 = sortidx(lst, idx=[1,3]); // Returns: [3,0,2,1]
+function sortidx(list, idx=undef) =
+    assert( is_list(list) || is_string(list) , "Invalid input to sort." )
+    assert( is_undef(idx) || is_finite(idx) || is_vector(idx) , "Invalid indices.")
+    list==[] ? [] : 
+    let(
+        size = array_dim(list),
+        aug = is_undef(idx) && (len(size) == 1 || (len(size) == 2 && size[1]<=4))
+              ?  zip(list, list_range(len(list)))
+              :  enumerate(list,idx=idx)
+    )
+    is_undef(idx) && len(size) == 1? subindex(_sort_vectors1(aug),1) :
+    is_undef(idx) && len(size) == 2 && size[1] <=4
+    ? (
+        size[1]==0 ? list_range(len(arr)) :
+        size[1]==1 ? subindex(_sort_vectors1(aug),1) :
+        size[1]==2 ? subindex(_sort_vectors2(aug),2) :
+        size[1]==3 ? subindex(_sort_vectors3(aug),3)
+    /*size[1]==4*/ : subindex(_sort_vectors4(aug),4)
+      ) 
+    :   // general case
+        subindex(_sort_general(aug, idx=list_range(s=1,n=len(aug)-1)), 0);
+
 function sortidx(list, idx=undef) =
     list==[] ? [] : let(
         size = array_dim(list),
@@ -812,6 +871,8 @@ function sortidx(list, idx=undef) =
     // general case
     subindex(_sort_general(aug, idx=list_range(s=1,n=len(aug)-1)), 0);
 
+// sort() does not accept strings but sortidx does; isn't inconsistent ?
+
 
 // Function: unique()
 // Usage:
@@ -821,16 +882,16 @@ function sortidx(list, idx=undef) =
 // Arguments:
 //   arr = The list to uniquify.
 function unique(arr) =
-    assert(is_list(arr)||is_string(arr))
-    len(arr)<=1? arr : let(
-        sorted = sort(arr)
-    ) [
-        for (i=[0:1:len(sorted)-1])
+    assert(is_list(arr)||is_string(arr), "Invalid input." )
+    len(arr)<=1? arr : 
+    let( sorted = sort(arr))
+    [   for (i=[0:1:len(sorted)-1])
             if (i==0 || (sorted[i] != sorted[i-1]))
                 sorted[i]
     ];
 
 
+
 // Function: unique_count()
 // Usage:
 //   unique_count(arr);
@@ -840,12 +901,14 @@ function unique(arr) =
 // Arguments:
 //   arr = The list to analyze. 
 function unique_count(arr) =
-      assert(is_list(arr) || is_string(arr))
+      assert(is_list(arr) || is_string(arr), "Invalid input." )
       arr == [] ? [[],[]] : 
       let( arr=sort(arr) )
-      let(ind = [0,for(i=[1:1:len(arr)-1]) if (arr[i]!=arr[i-1]) i])
-      [select(arr,ind),
-       deltas(concat(ind,[len(arr)]))];
+      let( ind = [0, for(i=[1:1:len(arr)-1]) if (arr[i]!=arr[i-1]) i] )
+      [ select(arr,ind), deltas( concat(ind,[len(arr)]) ) ];
+
+
+
 
 // Section: List Iteration Helpers
 
@@ -864,7 +927,7 @@ function unique_count(arr) =
 //   colors = ["red", "green", "blue"];
 //   for (i=idx(colors)) right(20*i) color(colors[i]) circle(d=10);
 function idx(list, step=1, end=-1,start=0) =
-    assert(is_list(list)||is_string(list))
+    assert(is_list(list)||is_string(list), "Invalid input." )
     [start : step : len(list)+end];
 
 
@@ -883,10 +946,11 @@ function idx(list, step=1, end=-1,start=0) =
 //   colors = ["red", "green", "blue"];
 //   for (p=enumerate(colors)) right(20*p[0]) color(p[1]) circle(d=10);
 function enumerate(l,idx=undef) =
-    assert(is_list(l)||is_string(list))
-    (idx==undef)?
-        [for (i=[0:1:len(l)-1]) [i,l[i]]] :
-        [for (i=[0:1:len(l)-1]) concat([i], [for (j=idx) l[i][j]])];
+    assert(is_list(l)||is_string(list), "Invalid input." )
+    assert(is_undef(idx)||is_finite(idx)||is_vector(idx) ||is_range(idx), "Invalid index/indices." )
+    (idx==undef)
+    ?   [for (i=[0:1:len(l)-1]) [i,l[i]]]
+    :   [for (i=[0:1:len(l)-1]) concat([i], [for (j=idx) l[i][j]])];
 
 
 // Function: force_list()
@@ -909,8 +973,7 @@ function enumerate(l,idx=undef) =
 //   w = force_list(4, n=3, fill=1);  // Returns: [4,1,1]
 function force_list(value, n=1, fill) =
     is_list(value) ? value :
-    is_undef(fill)? [for (i=[1:1:n]) value] :
-    [value, for (i=[2:1:n]) fill];
+    is_undef(fill)? [for (i=[1:1:n]) value] : [value, for (i=[2:1:n]) fill];
 
 
 // Function: pair()
@@ -927,7 +990,7 @@ function force_list(value, n=1, fill) =
 //   l = ["A","B","C","D"];
 //   echo([for (p=pair(l)) str(p.y,p.x)]);  // Outputs: ["BA", "CB", "DC"]
 function pair(v) =
-    assert(is_list(v)||is_string(v))
+    assert(is_list(v)||is_string(v), "Invalid input." )
     [for (i=[0:1:len(v)-2]) [v[i],v[i+1]]];
 
 
@@ -935,8 +998,8 @@ function pair(v) =
 // Usage:
 //   pair_wrap(v)
 // Description:
-//   Takes a list, and returns a list of adjacent pairss from it, wrapping around from the end to the start of the list.
-// Example(2D): Note that the last point and first point DO get paired together.
+//   Takes a list, and returns a list of adjacent pairs from it, wrapping around from the end to the start of the list.
+// Example(2D): 
 //   for (p = pair_wrap(circle(d=20, $fn=12)))
 //       move(p[0])
 //           rot(from=BACK, to=p[1]-p[0])
@@ -945,7 +1008,7 @@ function pair(v) =
 //   l = ["A","B","C","D"];
 //   echo([for (p=pair_wrap(l)) str(p.y,p.x)]);  // Outputs: ["BA", "CB", "DC", "AD"]
 function pair_wrap(v) =
-    assert(is_list(v)||is_string(v))
+    assert(is_list(v)||is_string(v), "Invalid input." )
     [for (i=[0:1:len(v)-1]) [v[i],v[(i+1)%len(v)]]];
 
 
@@ -958,7 +1021,7 @@ function pair_wrap(v) =
 //   l = ["A","B","C","D","E"];
 //   echo([for (p=triplet(l)) str(p.z,p.y,p.x)]);  // Outputs: ["CBA", "DCB", "EDC"]
 function triplet(v) =
-    assert(is_list(v)||is_string(v))
+    assert(is_list(v)||is_string(v), "Invalid input." )
     [for (i=[0:1:len(v)-3]) [v[i],v[i+1],v[i+2]]];
 
 
@@ -971,7 +1034,7 @@ function triplet(v) =
 //   l = ["A","B","C","D"];
 //   echo([for (p=triplet_wrap(l)) str(p.z,p.y,p.x)]);  // Outputs: ["CBA", "DCB", "ADC", "BAD"]
 function triplet_wrap(v) =
-    assert(is_list(v)||is_string(v))
+    assert(is_list(v)||is_string(v), "Invalid input." )
     [for (i=[0:1:len(v)-1]) [v[i],v[(i+1)%len(v)],v[(i+2)%len(v)]]];
 
 
@@ -991,10 +1054,11 @@ function triplet_wrap(v) =
 // Example(2D):
 //   for (p=permute(regular_ngon(n=7,d=100))) stroke(p);
 function permute(l,n=2,_s=0) =
-    assert(is_list(l))
-    assert(len(l)-_s >= n)
-    n==1? [for (i=[_s:1:len(l)-1]) [l[i]]] :
-    [for (i=[_s:1:len(l)-n], p=permute(l,n=n-1,_s=i+1)) concat([l[i]], p)];
+    assert(is_list(l), "Invalid list." )
+    assert( is_finite(n) && n>=1 && n<=len(l), "Invalid number `n`." )
+    n==1
+    ?   [for (i=[_s:1:len(l)-1]) [l[i]]] 
+    :   [for (i=[_s:1:len(l)-n], p=permute(l,n=n-1,_s=i+1)) concat([l[i]], p)];
 
 
 
@@ -1020,27 +1084,28 @@ function permute(l,n=2,_s=0) =
 //   set_v = set_union(set_a, set_b, get_indices=true);
 //   // set_v now equals [[5,0,1,2,6], [2,3,5,7,11,1,8]]
 function set_union(a, b, get_indices=false) =
+    assert( is_list(a) && is_list(b), "Invalid sets." )
     let(
         found1 = search(b, a),
         found2 = search(b, b),
-        c = [
-            for (i=idx(b))
-            if (found1[i] == [] && found2[i] == i)
-            b[i]
-        ],
+        c = [ for (i=idx(b))
+                if (found1[i] == [] && found2[i] == i)
+                    b[i] 
+            ],
         nset = concat(a, c)
-    ) !get_indices? nset :
+    ) 
+    ! get_indices ? nset :
     let(
         la = len(a),
         found3 = search(b, c),
-        idxs = [
-            for (i=idx(b))
-            (found1[i] != [])? found1[i] :
-            la + found3[i]
-        ]
+        idxs =  [ for (i=idx(b))
+                    (found1[i] != [])? found1[i] : la + found3[i]
+                ]
     ) [idxs, nset];
 
 
+
+
 // Function: set_difference()
 // Usage:
 //   s = set_difference(a, b);
@@ -1055,9 +1120,10 @@ function set_union(a, b, get_indices=false) =
 //   set_d = set_difference(set_a, set_b);
 //   // set_d now equals [7,11]
 function set_difference(a, b) =
-    let(
-        found = search(a, b, num_returns_per_match=1)
-    ) [ for (i=idx(a)) if(found[i]==[]) a[i] ];
+    assert( is_list(a) && is_list(b), "Invalid sets." )
+    let( found = search(a, b, num_returns_per_match=1) )
+    [ for (i=idx(a)) if(found[i]==[]) a[i] ];
+
 
 
 // Function: set_intersection()
@@ -1074,14 +1140,32 @@ function set_difference(a, b) =
 //   set_i = set_intersection(set_a, set_b);
 //   // set_i now equals [2,3,5]
 function set_intersection(a, b) =
-    let(
-        found = search(a, b, num_returns_per_match=1)
-    ) [ for (i=idx(a)) if(found[i]!=[]) a[i] ];
+    assert( is_list(a) && is_list(b), "Invalid sets." )
+    let( found = search(a, b, num_returns_per_match=1) )
+    [ for (i=idx(a)) if(found[i]!=[]) a[i] ];
+
 
 
 
 // Section: Array Manipulation
 
+// Function: add_scalar()
+// Usage:  
+//   add_scalar(v,s);
+// Description:
+//   Given an array and a scalar, returns the array with the scalar added to each item in it.
+//   If given a list of arrays, recursively adds the scalar to the each array.
+// Arguments:
+//   v = The initial array.
+//   s = A scalar value to add to every item in the array.
+// Example:
+//   add_scalar([1,2,3],3);            // Returns: [4,5,6]
+//   add_scalar([[1,2,3],[3,4,5]],3);  // Returns: [[4,5,6],[6,7,8]]
+function add_scalar(v,s) = 
+    is_finite(s) ? [for (x=v) is_list(x)? add_scalar(x,s) : is_finite(x) ? x+s: x] : v;
+
+
+
 // Function: subindex()
 // Description:
 //   For each array item, return the indexed subitem.
@@ -1096,10 +1180,11 @@ function set_intersection(a, b) =
 //   subindex(v,2);      // Returns [3, 7, 11, 15]
 //   subindex(v,[2,1]);  // Returns [[3, 2], [7, 6], [11, 10], [15, 14]]
 //   subindex(v,[1:3]);  // Returns [[2, 3, 4], [6, 7, 8], [10, 11, 12], [14, 15, 16]]
-function subindex(v, idx) = [
-    for(val=v) let(value=[for(i=idx) val[i]])
+function subindex(v, idx) = 
+    [ for(val=v) 
+        let( value=[for(i=idx) val[i]] )
         len(value)==1 ? value[0] : value
-];
+    ];
 
 
 // Function: zip()
@@ -1131,15 +1216,16 @@ function subindex(v, idx) = [
 function zip(vecs, v2, v3, fit=false, fill=undef) =
     (v3!=undef)? zip([vecs,v2,v3], fit=fit, fill=fill) :
     (v2!=undef)? zip([vecs,v2], fit=fit, fill=fill) :
-    assert(in_list(fit, [false, "short", "long"]))
+    assert(in_list(fit, [false, "short", "long"]), "Invalid fit value." )
     assert(all([for(v=vecs) is_list(v)]), "One of the inputs to zip is not a list")
     let(
         minlen = list_shortest(vecs),
-        maxlen = list_longest(vecs),
-        dummy = (fit==false)? assert(minlen==maxlen, "Input vectors to zip must have the same length") : 0
-    ) (fit == "long")?
-        [for(i=[0:1:maxlen-1]) [for(v=vecs) for(x=(i<len(v)? v[i] : (fill==undef)? [fill] : fill)) x] ] :
-        [for(i=[0:1:minlen-1]) [for(v=vecs) for(x=v[i]) x] ];
+        maxlen = list_longest(vecs)
+    )
+    assert(fit!=false || minlen==maxlen, "Input vectors to zip must have the same length")
+    (fit == "long")
+    ?   [for(i=[0:1:maxlen-1]) [for(v=vecs) for(x=(i<len(v)? v[i] : (fill==undef)? [fill] : fill)) x] ] 
+    :   [for(i=[0:1:minlen-1]) [for(v=vecs) for(x=v[i]) x] ];
 
 
 // Function: array_group()
@@ -1167,17 +1253,29 @@ function array_group(v, cnt=2, dflt=0) = [for (i = [0:cnt:len(v)-1]) [for (j = [
 function flatten(l) = [for (a = l) each a];
 
 
+// Function: full_flatten()
+// Description: 
+//   Collects in a list all elements recursively found in any level of the given list.
+//   The output list is ordered in depth first order.
+// Arguments:
+//   l = List to flatten.
+// Example:
+//   full_flatten([[1,2,3], [4,5,[6,7,8]]]) returns [1,2,3,4,5,6,7,8]
+function full_flatten(l) = [for(a=l) if(is_list(a)) (each full_flatten(a)) else a ];
+
+
 // Internal.  Not exposed.
 function _array_dim_recurse(v) =
-    !is_list(v[0])?  (
-        sum( [for(entry=v) is_list(entry) ? 1 : 0]) == 0 ? [] : [undef]
-    ) : let(
-        firstlen = len(v[0]),
-        first = sum( [for(entry = v) len(entry) == firstlen  ? 0 : 1]   ) == 0 ? firstlen : undef,
-        leveldown = flatten(v)
-    ) is_list(leveldown[0])? (
-        concat([first],_array_dim_recurse(leveldown))
-    ) : [first];
+    !is_list(v[0])
+    ?   sum( [for(entry=v) is_list(entry) ? 1 : 0] ) == 0 ? [] : [undef]
+    :   let(
+          firstlen = len(v[0]),
+          first = sum( [for(entry = v) len(entry) == firstlen  ? 0 : 1]   ) == 0 ? firstlen : undef,
+          leveldown = flatten(v)
+        ) 
+        is_list(leveldown[0])
+        ?  concat([first],_array_dim_recurse(leveldown))
+        : [first];
 
 
 // Function: array_dim()
@@ -1201,15 +1299,16 @@ function _array_dim_recurse(v) =
 //   array_dim([[[1,2,3],[4,5,6]],[[7,8,9],[10,11,12]]], 2);  // Returns 3
 //   array_dim([[[1,2,3],[4,5,6]],[[7,8,9]]]);                // Returns [2,undef,3]
 function array_dim(v, depth=undef) =
-    (depth == undef)? (
-        concat([len(v)], _array_dim_recurse(v))
-    ) : (depth == 0)? (
-        len(v)
-    ) : (
-        let(dimlist = _array_dim_recurse(v))
-        (depth > len(dimlist))? 0 : dimlist[depth-1]
-    );
+    assert( is_undef(depth) || ( is_finite(depth) && depth>=0 ), "Invalid depth.")
+    ! is_list(v) ? 0 :
+    (depth == undef)
+    ?   concat([len(v)], _array_dim_recurse(v))
+    :   (depth == 0)
+        ?  len(v)
+        :  let( dimlist = _array_dim_recurse(v))
+           (depth > len(dimlist))? 0 : dimlist[depth-1] ;
 
+// This function may return undef!
 
 
 // Function: transpose()
@@ -1242,7 +1341,12 @@ function array_dim(v, depth=undef) =
 // Example:
 //   transpose([3,4,5]);  // Returns: [3,4,5]
 function transpose(arr) =
-    is_list(arr[0])? [for (i=[0:1:len(arr[0])-1]) [for (j=[0:1:len(arr)-1]) arr[j][i]]] : arr;
+    let( a0 = arr[0] )
+    is_list(a0)
+    ?   assert([for(a=arr) if(len(a)!=len(a0)) 1]==[], "The array is not a matrix." )
+        [for (i=[0:1:len(a0)-1]) 
+          [ for (j=[0:1:len(arr)-1]) arr[j][i] ] ] 
+    :  arr;
 
 
 
diff --git a/common.scad b/common.scad
index 648f3c7..c72c77b 100644
--- a/common.scad
+++ b/common.scad
@@ -15,7 +15,8 @@
 // Usage:
 //   typ = typeof(x);
 // Description:
-//   Returns a string representing the type of the value.  One of "undef", "boolean", "number", "nan", "string", "list", or "range"
+//   Returns a string representing the type of the value.  One of "undef", "boolean", "number", "nan", "string", "list", "range" or "invalid".
+//   Some malformed "ranges", like '[0:NAN:INF]' and '[0:"a":INF]', may be classified as "undef" or "invalid".
 function typeof(x) =
     is_undef(x)? "undef" :
     is_bool(x)? "boolean" :
@@ -23,7 +24,9 @@ function typeof(x) =
     is_nan(x)? "nan" :
     is_string(x)? "string" :
     is_list(x)? "list" :
-    "range";
+    is_range(x) ? "range" :
+    "invalid";
+
 
 
 // Function: is_type()
@@ -70,8 +73,8 @@ function is_str(x) = is_string(x);
 //   is_int(n)
 // Description:
 //   Returns true if the given value is an integer (it is a number and it rounds to itself).  
-function is_int(n) = is_num(n) && n == round(n);
-function is_integer(n) = is_num(n) && n == round(n);
+function is_int(n) = is_finite(n) && n == round(n);
+function is_integer(n) = is_finite(n) && n == round(n);
 
 
 // Function: is_nan()
@@ -93,7 +96,17 @@ function is_finite(v) = is_num(0*v);
 // Function: is_range()
 // Description:
 //   Returns true if its argument is a range
-function is_range(x) = is_num(x[0]) && !is_list(x);
+function is_range(x) = !is_list(x) && is_finite(x[0]+x[1]+x[2]) ;
+
+
+// Function: valid_range()
+// Description:
+//   Returns true if its argument is a valid range (deprecated ranges excluded).
+function valid_range(x) = 
+    is_range(x) 
+    && ( x[1]>0 
+         ? x[0]<=x[2]
+         : ( x[1]<0 && x[0]>=x[2] ) );
 
 
 // Function: is_list_of()
@@ -106,13 +119,15 @@ function is_range(x) = is_num(x[0]) && !is_list(x);
 //   is_list_of([3,4,5], 0);            // Returns true
 //   is_list_of([3,4,undef], 0);        // Returns false
 //   is_list_of([[3,4],[4,5]], [1,1]);  // Returns true
+//   is_list_of([[3,"a"],[4,true]], [1,undef]);  // Returns true
 //   is_list_of([[3,4], 6, [4,5]], [1,1]);  // Returns false
-//   is_list_of([[1,[3,4]], [4,[5,6]]], [1,[2,3]]);    // Returne true
-//   is_list_of([[1,[3,INF]], [4,[5,6]]], [1,[2,3]]);  // Returne false
+//   is_list_of([[1,[3,4]], [4,[5,6]]], [1,[2,3]]);    // Returns true
+//   is_list_of([[1,[3,INF]], [4,[5,6]]], [1,[2,3]]);  // Returns false
+//   is_list_of([], [1,[2,3]]);                        // Returns true
 function is_list_of(list,pattern) =
     let(pattern = 0*pattern)
     is_list(list) &&
-    []==[for(entry=list) if (entry*0 != pattern) entry];
+    []==[for(entry=0*list) if (entry != pattern) entry];
 
 
 // Function: is_consistent()
@@ -128,7 +143,15 @@ function is_list_of(list,pattern) =
 //   is_consistent([[3,[3,4,[5]]], [5,[2,9,[9]]]]); // Returns true
 //   is_consistent([[3,[3,4,[5]]], [5,[2,9,9]]]);   // Returns false
 function is_consistent(list) =
-  is_list(list) && is_list_of(list, list[0]);
+  is_list(list) && is_list_of(list, _list_pattern(list[0]));
+
+
+//Internal function
+//Creates a list with the same structure of `list` with each of its elements substituted by 0.
+function _list_pattern(list) =
+  is_list(list) 
+  ? [for(entry=list) is_list(entry) ? _list_pattern(entry) : 0]
+  : 0;
 
 
 // Function: same_shape()
@@ -139,7 +162,7 @@ function is_consistent(list) =
 // Example:
 //   same_shape([3,[4,5]],[7,[3,4]]);   // Returns true
 //   same_shape([3,4,5], [7,[3,4]]);    // Returns false
-function same_shape(a,b) = a*0 == b*0;
+function same_shape(a,b) = _list_pattern(a) == b*0;
 
 
 // Section: Handling `undef`s.
@@ -311,9 +334,10 @@ function scalar_vec3(v, dflt=undef) =
 //   Calculate the standard number of sides OpenSCAD would give a circle based on `$fn`, `$fa`, and `$fs`.
 // Arguments:
 //   r = Radius of circle to get the number of segments for.
-function segs(r) =
+function segs(r) = 
     $fn>0? ($fn>3? $fn : 3) :
-    ceil(max(5, min(360/$fa, abs(r)*2*PI/$fs)));
+    let( r = is_finite(r)? r: 0 ) 
+    ceil(max(5, min(360/$fa, abs(r)*2*PI/$fs))) ;
 
 
 
@@ -322,7 +346,7 @@ function segs(r) =
 
 function _valstr(x) =
     is_list(x)? str("[",str_join([for (xx=x) _valstr(xx)],","),"]") :
-    is_num(x)? fmt_float(x,12) : x;
+    is_finite(x)? fmt_float(x,12) : x;
 
 
 // Module: assert_approx()
diff --git a/math.scad b/math.scad
index 82b08bb..3ff62d7 100644
--- a/math.scad
+++ b/math.scad
@@ -33,7 +33,10 @@ NAN = acos(2);  // The value `nan`, useful for comparisons.
 //   sqr([3,4]); // Returns: [9,16]
 //   sqr([[1,2],[3,4]]);  // Returns [[1,4],[9,16]]
 //   sqr([[1,2],3]);      // Returns [[1,4],9]
-function sqr(x) = is_list(x) ? [for(val=x) sqr(val)] : x*x;
+function sqr(x) = 
+    is_list(x) ? [for(val=x) sqr(val)] : 
+    is_finite(x) ? x*x :
+    assert(is_finite(x) || is_vector(x), "Input is not neither a number nor a list of numbers.");
 
 
 // Function: log2()
@@ -45,8 +48,11 @@ function sqr(x) = is_list(x) ? [for(val=x) sqr(val)] : x*x;
 //   log2(0.125);  // Returns: -3
 //   log2(16);     // Returns: 4
 //   log2(256);    // Returns: 8
-function log2(x) = ln(x)/ln(2);
+function log2(x) = 
+    assert( is_finite(x), "Input is not a number.")
+    ln(x)/ln(2);
 
+// this may return NAN or INF; should it check x>0 ?
 
 // Function: hypot()
 // Usage:
@@ -60,7 +66,9 @@ function log2(x) = ln(x)/ln(2);
 // Example:
 //   l = hypot(3,4);  // Returns: 5
 //   l = hypot(3,4,5);  // Returns: ~7.0710678119
-function hypot(x,y,z=0) = norm([x,y,z]);
+function hypot(x,y,z=0) = 
+    assert( is_vector([x,y,z]), "Improper number(s).")
+    norm([x,y,z]);
 
 
 // Function: factorial()
@@ -76,11 +84,53 @@ function hypot(x,y,z=0) = norm([x,y,z]);
 //   y = factorial(6);  // Returns: 720
 //   z = factorial(9);  // Returns: 362880
 function factorial(n,d=0) =
-    assert(n>=0 && d>=0, "Factorial is not defined for negative numbers")
+    assert(is_int(n) && is_int(d) && n>=0 && d>=0, "Factorial is not defined for negative numbers")
     assert(d<=n, "d cannot be larger than n")
     product([1,for (i=[n:-1:d+1]) i]);
 
 
+// Function: binomial()
+// Usage:
+//   x = binomial(n);
+// Description:
+//   Returns the binomial coefficients of the integer `n`.  
+// Arguments:
+//   n = The integer to get the binomial coefficients of
+// Example:
+//   x = binomial(3);  // Returns: [1,3,3,1]
+//   y = binomial(4);  // Returns: [1,4,6,4,1]
+//   z = binomial(6);  // Returns: [1,6,15,20,15,6,1]
+function binomial(n) =
+    assert( is_int(n) && n>0, "Input is not an integer greater than 0.")
+    [for( c = 1, i = 0; 
+        i<=n; 
+         c = c*(n-i)/(i+1), i = i+1
+        ) c ] ;
+
+
+// Function: binomial_coefficient()
+// Usage:
+//   x = binomial_coefficient(n,k);
+// Description:
+//   Returns the k-th binomial coefficient of the integer `n`.  
+// Arguments:
+//   n = The integer to get the binomial coefficient of
+//   k = The binomial coefficient index
+// Example:
+//   x = binomial_coefficient(3,2);  // Returns: 3
+//   y = binomial_coefficient(10,6); // Returns: 210
+function binomial_coefficient(n,k) =
+    assert( is_int(n) && is_int(k), "Some input is not a number.")
+    k < 0 || k > n ? 0 :
+    k ==0 || k ==n ? 1 :
+    let( k = min(k, n-k),
+         b = [for( c = 1, i = 0; 
+                   i<=k; 
+                   c = c*(n-i)/(i+1), i = i+1
+                 ) c] )
+    b[len(b)-1];
+
+
 // Function: lerp()
 // Usage:
 //   x = lerp(a, b, u);
@@ -91,8 +141,8 @@ function factorial(n,d=0) =
 //   If `u` is 0.0, then the value of `a` is returned.
 //   If `u` is 1.0, then the value of `b` is returned.
 //   If `u` is a range, or list of numbers, returns a list of interpolated values.
-//   It is valid to use a `u` value outside the range 0 to 1.  The result will be a predicted
-//   value along the slope formed by `a` and `b`, but not between those two values.
+//   It is valid to use a `u` value outside the range 0 to 1.  The result will be an extrapolation
+//   along the slope formed by `a` and `b`.
 // Arguments:
 //   a = First value or vector.
 //   b = Second value or vector.
@@ -113,9 +163,9 @@ function factorial(n,d=0) =
 //   rainbow(pts) translate($item) circle(d=3,$fn=8);
 function lerp(a,b,u) =
     assert(same_shape(a,b), "Bad or inconsistent inputs to lerp")
-    is_num(u)? (1-u)*a + u*b :
-    assert(!is_undef(u)&&!is_bool(u)&&!is_string(u), "Input u to lerp must be a number, vector, or range.")
-    [for (v = u) lerp(a,b,v)];
+    is_finite(u)? (1-u)*a + u*b :
+    assert(is_finite(u) || is_vector(u) || valid_range(u), "Input u to lerp must be a number, vector, or range.")
+    [for (v = u) (1-v)*a + v*b ];
 
 
 
@@ -124,40 +174,45 @@ function lerp(a,b,u) =
 // Function: sinh()
 // Description: Takes a value `x`, and returns the hyperbolic sine of it.
 function sinh(x) =
+    assert(is_finite(x), "The input must be a finite number.")
     (exp(x)-exp(-x))/2;
 
 
 // Function: cosh()
 // Description: Takes a value `x`, and returns the hyperbolic cosine of it.
 function cosh(x) =
+    assert(is_finite(x), "The input must be a finite number.")
     (exp(x)+exp(-x))/2;
 
 
 // Function: tanh()
 // Description: Takes a value `x`, and returns the hyperbolic tangent of it.
 function tanh(x) =
+    assert(is_finite(x), "The input must be a finite number.")
     sinh(x)/cosh(x);
 
 
 // Function: asinh()
 // Description: Takes a value `x`, and returns the inverse hyperbolic sine of it.
 function asinh(x) =
+    assert(is_finite(x), "The input must be a finite number.")
     ln(x+sqrt(x*x+1));
 
 
 // Function: acosh()
 // Description: Takes a value `x`, and returns the inverse hyperbolic cosine of it.
 function acosh(x) =
+    assert(is_finite(x), "The input must be a finite number.")
     ln(x+sqrt(x*x-1));
 
 
 // Function: atanh()
 // Description: Takes a value `x`, and returns the inverse hyperbolic tangent of it.
 function atanh(x) =
+    assert(is_finite(x), "The input must be a finite number.")
     ln((1+x)/(1-x))/2;
 
 
-
 // Section: Quantization
 
 // Function: quant()
@@ -185,8 +240,11 @@ function atanh(x) =
 //   quant([9,10,10.4,10.5,11,12],3);      // Returns: [9,9,9,12,12,12]
 //   quant([[9,10,10.4],[10.5,11,12]],3);  // Returns: [[9,9,9],[12,12,12]]
 function quant(x,y) =
-    is_list(x)? [for (v=x) quant(v,y)] :
-    floor(x/y+0.5)*y;
+    assert(is_finite(y) && !approx(y,0,eps=1e-24), "The multiple must be a non zero integer.")
+    is_list(x)
+    ?   [for (v=x) quant(v,y)]
+    :   assert( is_finite(x), "The input to quantize must be a number or a list of numbers.")
+        floor(x/y+0.5)*y;
 
 
 // Function: quantdn()
@@ -214,8 +272,11 @@ function quant(x,y) =
 //   quantdn([9,10,10.4,10.5,11,12],3);      // Returns: [9,9,9,9,9,12]
 //   quantdn([[9,10,10.4],[10.5,11,12]],3);  // Returns: [[9,9,9],[9,9,12]]
 function quantdn(x,y) =
-    is_list(x)? [for (v=x) quantdn(v,y)] :
-    floor(x/y)*y;
+    assert(is_finite(y) && !approx(y,0,eps=1e-24), "The multiple must be a non zero integer.")
+    is_list(x)
+    ?    [for (v=x) quantdn(v,y)]
+    :    assert( is_finite(x), "The input to quantize must be a number or a list of numbers.")
+        floor(x/y)*y;
 
 
 // Function: quantup()
@@ -243,8 +304,11 @@ function quantdn(x,y) =
 //   quantup([9,10,10.4,10.5,11,12],3);      // Returns: [9,12,12,12,12,12]
 //   quantup([[9,10,10.4],[10.5,11,12]],3);  // Returns: [[9,12,12],[12,12,12]]
 function quantup(x,y) =
-    is_list(x)? [for (v=x) quantup(v,y)] :
-    ceil(x/y)*y;
+    assert(is_finite(y) && !approx(y,0,eps=1e-24), "The multiple must be a non zero integer.")
+    is_list(x)
+    ?    [for (v=x) quantup(v,y)]
+    :    assert( is_finite(x), "The input to quantize must be a number or a list of numbers.")
+        ceil(x/y)*y;
 
 
 // Section: Constraints and Modulos
@@ -264,7 +328,9 @@ function quantup(x,y) =
 //   constrain(0.3, -1, 1);  // Returns: 0.3
 //   constrain(9.1, 0, 9);   // Returns: 9
 //   constrain(-0.1, 0, 9);  // Returns: 0
-function constrain(v, minval, maxval) = min(maxval, max(minval, v));
+function constrain(v, minval, maxval) = 
+    assert( is_finite(v+minval+maxval), "Input must be finite number(s).")
+    min(maxval, max(minval, v));
 
 
 // Function: posmod()
@@ -283,7 +349,9 @@ function constrain(v, minval, maxval) = min(maxval, max(minval, v));
 //   posmod(270,360);   // Returns: 270
 //   posmod(700,360);   // Returns: 340
 //   posmod(3,2.5);     // Returns: 0.5
-function posmod(x,m) = (x%m+m)%m;
+function posmod(x,m) = 
+    assert( is_finite(x) && is_finite(m) && !approx(m,0) , "Input must be finite numbers. The divisor cannot be zero.")
+    (x%m+m)%m;
 
 
 // Function: modang(x)
@@ -299,6 +367,7 @@ function posmod(x,m) = (x%m+m)%m;
 //   modang(270,360);   // Returns: -90
 //   modang(700,360);   // Returns: -20
 function modang(x) =
+    assert( is_finite(x), "Input must be a finite number.")
     let(xx = posmod(x,360)) xx<180? xx : xx-360;
 
 
@@ -306,7 +375,7 @@ function modang(x) =
 // Usage:
 //   modrange(x, y, m, [step])
 // Description:
-//   Returns a normalized list of values from `x` to `y`, by `step`, modulo `m`.  Wraps if `x` > `y`.
+//   Returns a normalized list of numbers from `x` to `y`, by `step`, modulo `m`.  Wraps if `x` > `y`.
 // Arguments:
 //   x = The start value to constrain.
 //   y = The end value to constrain.
@@ -318,6 +387,7 @@ function modang(x) =
 //   modrange(90,270,360, step=-45);  // Returns: [90,45,0,315,270]
 //   modrange(270,90,360, step=-45);  // Returns: [270,225,180,135,90]
 function modrange(x, y, m, step=1) =
+    assert( is_finite(x+y+step+m) && !approx(m,0), "Input must be finite numbers. The module value cannot be zero.")
     let(
         a = posmod(x, m),
         b = posmod(y, m),
@@ -330,20 +400,21 @@ function modrange(x, y, m, step=1) =
 
 // Function: rand_int()
 // Usage:
-//   rand_int(min,max,N,[seed]);
+//   rand_int(minval,maxval,N,[seed]);
 // Description:
-//   Return a list of random integers in the range of min to max, inclusive.
+//   Return a list of random integers in the range of minval to maxval, inclusive.
 // Arguments:
-//   min = Minimum integer value to return.
-//   max = Maximum integer value to return.
+//   minval = Minimum integer value to return.
+//   maxval = Maximum integer value to return.
 //   N = Number of random integers to return.
 //   seed = If given, sets the random number seed.
 // Example:
 //   ints = rand_int(0,100,3);
 //   int = rand_int(-10,10,1)[0];
-function rand_int(min, max, N, seed=undef) =
-    assert(max >= min, "Max value cannot be smaller than min")
-    let (rvect = is_def(seed) ? rands(min,max+1,N,seed) : rands(min,max+1,N))
+function rand_int(minval, maxval, N, seed=undef) =
+    assert( is_finite(minval+maxval+N) && (is_undef(seed) || is_finite(seed) ), "Input must be finite numbers.")
+    assert(maxval >= minval, "Max value cannot be smaller than minval")
+    let (rvect = is_def(seed) ? rands(minval,maxval+1,N,seed) : rands(minval,maxval+1,N))
     [for(entry = rvect) floor(entry)];
 
 
@@ -358,6 +429,7 @@ function rand_int(min, max, N, seed=undef) =
 //   N = Number of random numbers to return.  Default: 1
 //   seed = If given, sets the random number seed.
 function gaussian_rands(mean, stddev, N=1, seed=undef) =
+    assert( is_finite(mean+stddev+N) && (is_undef(seed) || is_finite(seed) ), "Input must be finite numbers.")
     let(nums = is_undef(seed)? rands(0,1,N*2) : rands(0,1,N*2,seed))
     [for (i = list_range(N)) mean + stddev*sqrt(-2*ln(nums[i*2]))*cos(360*nums[i*2+1])];
 
@@ -374,6 +446,10 @@ function gaussian_rands(mean, stddev, N=1, seed=undef) =
 //   N = Number of random numbers to return.  Default: 1
 //   seed = If given, sets the random number seed.
 function log_rands(minval, maxval, factor, N=1, seed=undef) =
+    assert( is_finite(minval+maxval+N) 
+            && (is_undef(seed) || is_finite(seed) )
+            && factor>0, 
+            "Input must be finite numbers. `factor` should be greater than zero.")
     assert(maxval >= minval, "maxval cannot be smaller than minval")
     let(
         minv = 1-1/pow(factor,minval),
@@ -395,18 +471,18 @@ function gcd(a,b) =
     b==0 ? abs(a) : gcd(b,a % b);
 
 
-// Computes lcm for two scalars
+// Computes lcm for two integers
 function _lcm(a,b) =
-    assert(is_int(a), "Invalid non-integer parameters to lcm")
-    assert(is_int(b), "Invalid non-integer parameters to lcm")
-    assert(a!=0 && b!=0, "Arguments to lcm must be nonzero")
+    assert(is_int(a) && is_int(b), "Invalid non-integer parameters to lcm")
+    assert(a!=0 && b!=0, "Arguments to lcm must be non zero")
     abs(a*b) / gcd(a,b);
 
 
 // Computes lcm for a list of values
 function _lcmlist(a) =
-    len(a)==1 ? a[0] :
-    _lcmlist(concat(slice(a,0,len(a)-2),[lcm(a[len(a)-2],a[len(a)-1])]));
+    len(a)==1 
+    ?   a[0] 
+    :   _lcmlist(concat(slice(a,0,len(a)-2),[lcm(a[len(a)-2],a[len(a)-1])]));
 
 
 // Function: lcm()
@@ -418,12 +494,11 @@ function _lcmlist(a) =
 //   be non-zero integers.  The output is always a positive integer.  It is an error to pass zero
 //   as an argument.  
 function lcm(a,b=[]) =
-    !is_list(a) && !is_list(b) ? _lcm(a,b) : 
-    let(
-        arglist = concat(force_list(a),force_list(b))
-    )
-    assert(len(arglist)>0,"invalid call to lcm with empty list(s)")
-    _lcmlist(arglist);
+    !is_list(a) && !is_list(b) 
+    ?   _lcm(a,b) 
+    :   let( arglist = concat(force_list(a),force_list(b)) )
+        assert(len(arglist)>0, "Invalid call to lcm with empty list(s)")
+        _lcmlist(arglist);
 
 
 
@@ -431,8 +506,9 @@ function lcm(a,b=[]) =
 
 // Function: sum()
 // Description:
-//   Returns the sum of all entries in the given list.
-//   If passed an array of vectors, returns a vector of sums of each part.
+//   Returns the sum of all entries in the given consistent list.
+//   If passed an array of vectors, returns the sum the vectors.
+//   If passed an array of matrices, returns the sum of the matrices.
 //   If passed an empty list, the value of `dflt` will be returned.
 // Arguments:
 //   v = The list to get the sum of.
@@ -441,11 +517,10 @@ 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_vector(v) ? [for(i=v) 1]*v :
+    is_list(v) && len(v) == 0 ? dflt :
+    is_vector(v) || is_matrix(v)? [for(i=v) 1]*v :
     assert(is_consistent(v), "Input to sum is non-numeric or inconsistent")
-    is_vector(v[0]) ? [for(i=v) 1]*v :
-    len(v) == 0 ? dflt :
-                  _sum(v,v[0]*0);
+    _sum(v,v[0]*0);
 
 function _sum(v,_total,_i=0) = _i>=len(v) ? _total : _sum(v,_total+v[_i], _i+1);
 
@@ -495,37 +570,51 @@ function sum_of_squares(v) = sum(vmul(v,v));
 // Examples:
 //   v = sum_of_sines(30, [[10,3,0], [5,5.5,60]]);
 function sum_of_sines(a, sines) =
-    sum([
-        for (s = sines) let(
-            ss=point3d(s),
-            v=ss.x*sin(a*ss.y+ss.z)
-        ) v
-    ]);
+    assert( is_finite(a) && is_matrix(sines,undef,3), "Invalid input.")
+    sum([ for (s = sines) 
+            let(
+              ss=point3d(s),
+              v=ss[0]*sin(a*ss[1]+ss[2])
+            ) v
+        ]);
 
 
 // Function: deltas()
 // 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).
 //   Given [a,b,c,d], returns [b-a,c-b,d-c].
 // Arguments:
 //   v = The list to get the deltas of.
 // Example:
 //   deltas([2,5,9,17]);  // returns [3,4,8].
 //   deltas([[1,2,3], [3,6,8], [4,8,11]]);  // returns [[2,4,5], [1,2,3]]
-function deltas(v) = [for (p=pair(v)) p.y-p.x];
+function deltas(v) = 
+    assert( is_consistent(v) && len(v)>1 , "Inconsistent list or with length<=1.")
+    [for (p=pair(v)) p[1]-p[0]] ;
 
 
 // Function: product()
 // Description:
 //   Returns the product of all entries in the given list.
-//   If passed an array of vectors, returns a vector of products of each part.
-//   If passed an array of matrices, returns a the resulting product matrix.
+//   If passed a list of vectors of same dimension, returns a vector of products of each part.
+//   If passed a list of square matrices, returns a the resulting product matrix.
 // Arguments:
 //   v = The list to get the product of.
 // Example:
 //   product([2,3,4]);  // returns 24.
 //   product([[1,2,3], [3,4,5], [5,6,7]]);  // returns [15, 48, 105]
-function product(v, i=0, tot=undef) = i>=len(v)? tot : product(v, i+1, ((tot==undef)? v[i] : is_vector(v[i])? vmul(tot,v[i]) : tot*v[i]));
+function product(v) = 
+    assert( is_vector(v) || is_matrix(v) || ( is_matrix(v[0],square=true) && is_consistent(v)), 
+            "Invalid input.")
+    _product(v, 1, v[0]);
+
+function _product(v, i=0, _tot) = 
+    i>=len(v) ? _tot :
+    _product( v, 
+              i+1, 
+              ( is_vector(v[i])? vmul(_tot,v[i]) : _tot*v[i] ) );
+               
 
 
 // Function: outer_product()
@@ -534,21 +623,22 @@ function product(v, i=0, tot=undef) = i>=len(v)? tot : product(v, i+1, ((tot==un
 // Usage:
 //   M = outer_product(u,v);
 function outer_product(u,v) =
-  assert(is_vector(u) && is_vector(v))
-  assert(len(u)==len(v))
-  [for(i=[0:len(u)-1]) [for(j=[0:len(u)-1]) u[i]*v[j]]];
+  assert(is_vector(u) && is_vector(v), "The inputs must be vectors.")
+  [for(ui=u) ui*v];
 
 
 // Function: mean()
 // Description:
-//   Returns the arithmatic mean/average of all entries in the given array.
+//   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.
 // Arguments:
 //   v = The list of values to get the mean of.
 // Example:
 //   mean([2,3,4]);  // returns 3.
 //   mean([[1,2,3], [3,4,5], [5,6,7]]);  // returns [3, 4, 5]
-function mean(v) = sum(v)/len(v);
+function mean(v) = 
+    assert(is_list(v) && len(v)>0, "Invalid list.")
+    sum(v)/len(v);
 
 
 // Function: median()
@@ -556,18 +646,33 @@ function mean(v) = sum(v)/len(v);
 //   x = median(v);
 // Description:
 //   Given a list of numbers or vectors, finds the median value or midpoint.
-//   If passed a list of vectors, returns the vector of the median of each part.
+//   If passed a list of vectors, returns the vector of the median of each component.
 function median(v) =
-    assert(is_list(v))
-    assert(len(v)>0)
-    is_vector(v[0])? (
-        assert(is_consistent(v))
-        [
-            for (i=idx(v[0]))
-            let(vals = subindex(v,i))
-            (min(vals)+max(vals))/2
-        ]
-    ) : (min(v)+max(v))/2;
+    is_vector(v) ? (min(v)+max(v))/2 :
+    is_matrix(v) ? [for(ti=transpose(v))  (min(ti)+max(ti))/2 ]
+    :   assert(false , "Invalid input.");
+
+// Function: convolve()
+// Usage:
+//   x = convolve(p,q);
+// Description:
+//   Given two vectors, finds the convolution of them.
+//   The length of the returned vector is len(p)+len(q)-1 .
+// Arguments:
+//   p = The first vector.
+//   q = The second vector.
+// Example:
+//   a = convolve([1,1],[1,2,1]); // Returns: [1,3,3,1]
+//   b = convolve([1,2,3],[1,2,1])); // Returns: [1,4,8,8,3]
+function convolve(p,q) =
+    p==[] || q==[] ? [] :
+    assert( is_vector(p) && is_vector(q), "The inputs should be vectors.")
+    let( n = len(p),
+         m = len(q))
+    [for(i=[0:n+m-2], k1 = max(0,i-n+1), k2 = min(i,m-1) )
+       [for(j=[k1:k2]) p[i-j] ] * [for(j=[k1:k2]) q[j] ] 
+    ];
+
 
 
 // Section: Matrix math
@@ -582,7 +687,7 @@ function median(v) =
 //   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) =
-    assert(is_matrix(A))
+    assert(is_matrix(A), "Input should be a matrix.")
     let(
         m = len(A),
         n = len(A[0])
@@ -619,8 +724,12 @@ function matrix_inverse(A) =
 // Description:
 //   Returns a submatrix with the specified index ranges or index sets.  
 function submatrix(M,ind1,ind2) =
-    [for(i=ind1) [for(j=ind2) M[i][j] ] ];
-
+    assert( is_matrix(M), "Input must be a matrix." )
+    [for(i=ind1) 
+        [for(j=ind2) 
+            assert( ! is_undef(M[i][j]), "Invalid indexing." )
+            M[i][j] ] ];
+    
 
 // Function: qr_factor()
 // Usage: qr = qr_factor(A)
@@ -628,7 +737,7 @@ function submatrix(M,ind1,ind2) =
 //   Calculates the QR factorization of the input matrix A and returns it as the list [Q,R].  This factorization can be
 //   used to solve linear systems of equations.  
 function qr_factor(A) =
-    assert(is_matrix(A))
+    assert(is_matrix(A), "Input must be a matrix." )
     let(
       m = len(A),
       n = len(A[0])
@@ -659,8 +768,8 @@ function _qr_factor(A,Q, column, m, n) =
 // Function: back_substitute()
 // Usage: back_substitute(R, b, [transpose])
 // Description:
-//   Solves the problem Rx=b where R is an upper triangular square matrix.  No check is made that the lower triangular entries
-//   are actually zero.  If transpose==true then instead solve transpose(R)*x=b.
+//   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.
 //   You can supply a compatible matrix b and it will produce the solution for every column of b.  Note that if you want to
 //   solve Rx=b1 and Rx=b2 you must set b to transpose([b1,b2]) and then take the transpose of the result.  If the matrix
 //   is singular (e.g. has a zero on the diagonal) then it returns [].  
@@ -694,7 +803,9 @@ function back_substitute(R, b, x=[],transpose = false) =
 // Example:
 //   M = [ [6,-2], [1,8] ];
 //   det = det2(M);  // Returns: 50
-function det2(M) = M[0][0] * M[1][1] - M[0][1]*M[1][0];
+function det2(M) = 
+    assert( is_matrix(M,2,2), "Matrix should be 2x2." )
+    M[0][0] * M[1][1] - M[0][1]*M[1][0];
 
 
 // Function: det3()
@@ -706,6 +817,7 @@ function det2(M) = M[0][0] * M[1][1] - M[0][1]*M[1][0];
 //   M = [ [6,4,-2], [1,-2,8], [1,5,7] ];
 //   det = det3(M);  // Returns: -334
 function det3(M) =
+    assert( is_matrix(M,3,3), "Matrix should 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]);
@@ -720,21 +832,21 @@ 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(len(M)==len(M[0]))
+    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) :
     sum(
         [for (col=[0:1:len(M)-1])
             ((col%2==0)? 1 : -1) *
-            M[col][0] *
-            determinant(
-                [for (r=[1:1:len(M)-1])
-                    [for (c=[0:1:len(M)-1])
-                        if (c!=col) M[c][r]
+                M[col][0] *
+                determinant(
+                    [for (r=[1:1:len(M)-1])
+                        [for (c=[0:1:len(M)-1])
+                            if (c!=col) M[c][r]
+                        ]
                     ]
-                ]
-            )
+                )
         ]
     );
 
@@ -753,8 +865,11 @@ function determinant(M) =
 //   n = optional width of matrix
 //   square = set to true to require a square matrix.  Default: false        
 function is_matrix(A,m,n,square=false) =
-    is_vector(A[0],n) && is_vector(A*(0*A[0]),m) &&
-    (!square || len(A)==len(A[0]));
+    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]));
 
 
 // Section: Comparisons and Logic
@@ -774,11 +889,13 @@ function is_matrix(A,m,n,square=false) =
 //   approx(0.3333,1/3);          // Returns: false
 //   approx(0.3333,1/3,eps=1e-3);  // Returns: true
 //   approx(PI,3.1415926536);     // Returns: true
-function approx(a,b,eps=EPSILON) =
+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] == []) :
-    (abs(a-b) <= eps);
+    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);
+    
 
 
 function _type_num(x) =
@@ -796,7 +913,7 @@ function _type_num(x) =
 // Description:
 //   Compares two values.  Lists are compared recursively.
 //   Returns <0 if a<b.  Returns >0 if a>b.  Returns 0 if a==b.
-//   If types are not the same, then undef < bool < num < str < list < range.
+//   If types are not the same, then undef < bool < nan < num < str < list < range.
 // Arguments:
 //   a = First value to compare.
 //   b = Second value to compare.
@@ -820,13 +937,14 @@ function compare_vals(a, b) =
 //   a = First list to compare.
 //   b = Second list to compare.
 function compare_lists(a, b) =
-    a==b? 0 : let(
-        cmps = [
-            for(i=[0:1:min(len(a),len(b))-1]) let(
-                cmp = compare_vals(a[i],b[i])
-            ) if(cmp!=0) cmp
-        ]
-    ) cmps==[]? (len(a)-len(b)) : cmps[0];
+    a==b? 0 
+    :   let(
+          cmps = [ for(i=[0:1:min(len(a),len(b))-1]) 
+                      let( cmp = compare_vals(a[i],b[i]) )
+                      if(cmp!=0) cmp 
+                 ]
+           ) 
+        cmps==[]? (len(a)-len(b)) : cmps[0];
 
 
 // Function: any()
@@ -843,12 +961,11 @@ function compare_lists(a, b) =
 //   any([[0,0], [1,0]]);   // Returns true.
 function any(l, i=0, succ=false) =
     (i>=len(l) || succ)? succ :
-    any(
-        l, i=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()
@@ -865,13 +982,12 @@ function any(l, i=0, succ=false) =
 //   all([[0,0], [1,0]]);   // Returns false.
 //   all([[1,1], [1,1]]);   // Returns true.
 function all(l, i=0, fail=false) =
-    (i>=len(l) || fail)? (!fail) :
-    all(
-        l, i=i+1, fail=(
-            is_list(l[i])? !all(l[i]) :
-            !l[i]
-        )
-    );
+    (i>=len(l) || fail)? !fail :
+    all( l, 
+         i+1,
+         fail = is_list(l[i]) ? !all(l[i]) : !l[i]
+        ) ;
+
 
 
 // Function: count_true()
@@ -904,6 +1020,21 @@ function count_true(l, nmax=undef, i=0, cnt=0) =
     );
 
 
+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];
+
+
 
 // Section: Calculus
 
@@ -921,42 +1052,49 @@ function count_true(l, nmax=undef, i=0, cnt=0) =
 //   between data[i+1] and data[i], and the data values will be linearly resampled at each corner
 //   to produce a uniform spacing for the derivative estimate.  At the endpoints a single point method
 //   is used: f'(t) = (f(t+h)-f(t))/h.  
+// Arguments:
+//   data = the list of the elements to compute the derivative of.
+//   h = the parametric sampling of the data.
+//   closed = boolean to indicate if the data set should be wrapped around from the end to the start.
 function deriv(data, h=1, closed=false) =
+    assert( is_consistent(data) , "Input list is not consistent or not numerical.") 
+    assert( len(data)>=2, "Input `data` should have at least 2 elements.") 
+    assert( is_finite(h) || is_vector(h), "The sampling `h` must be a number or a list of numbers." )
+    assert( is_num(h) || len(h) == len(data)-(closed?0:1),
+            str("Vector valued `h` must have length ",len(data)-(closed?0:1)))
     is_vector(h) ? _deriv_nonuniform(data, h, closed=closed) :
     let( L = len(data) )
-    closed? [
+    closed
+    ? [
         for(i=[0:1:L-1])
         (data[(i+1)%L]-data[(L+i-1)%L])/2/h
-    ] :
-    let(
-        first =
-            L<3? data[1]-data[0] : 
-            3*(data[1]-data[0]) - (data[2]-data[1]),
-        last =
-            L<3? data[L-1]-data[L-2]:
-            (data[L-3]-data[L-2])-3*(data[L-2]-data[L-1])
-    ) [
+      ]
+    : let(
+        first = L<3 ? data[1]-data[0] : 
+                3*(data[1]-data[0]) - (data[2]-data[1]),
+        last = L<3 ? data[L-1]-data[L-2]:
+               (data[L-3]-data[L-2])-3*(data[L-2]-data[L-1])
+         ) 
+      [
         first/2/h,
         for(i=[1:1:L-2]) (data[i+1]-data[i-1])/2/h,
         last/2/h
-    ];
+      ];
 
 
 function _dnu_calc(f1,fc,f2,h1,h2) =
     let(
         f1 = h2<h1 ? lerp(fc,f1,h2/h1) : f1 , 
         f2 = h1<h2 ? lerp(fc,f2,h1/h2) : f2
-    )
-    (f2-f1) / 2 / min([h1,h2]);
+       )
+    (f2-f1) / 2 / min(h1,h2);
 
 
 function _deriv_nonuniform(data, h, closed) =
-    assert(len(h) == len(data)-(closed?0:1),str("Vector valued h must be length ",len(data)-(closed?0:1)))
-    let(
-      L = len(data)
-    )
-    closed? [for(i=[0:1:L-1])
-    	        _dnu_calc(data[(L+i-1)%L], data[i], data[(i+1)%L], select(h,i-1), h[i]) ]
+    let( L = len(data) )
+    closed
+    ? [for(i=[0:1:L-1])
+          _dnu_calc(data[(L+i-1)%L], data[i], data[(i+1)%L], select(h,i-1), h[i]) ]
     : [
         (data[1]-data[0])/h[0],
         for(i=[1:1:L-2]) _dnu_calc(data[i-1],data[i],data[i+1], h[i-1],h[i]),
@@ -967,15 +1105,23 @@ function _deriv_nonuniform(data, h, closed) =
 // Function: deriv2()
 // Usage: deriv2(data, [h], [closed])
 // Description:
-//   Computes a numerical esimate of the second derivative of the data, which may be scalar or vector valued.
+//   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. 
 //   If the `closed` parameter is true the data is assumed to be defined on a loop with data[0] adjacent to
 //   data[len(data)-1].  For internal points this function uses the approximation 
-//   f''(t) = (f(t-h)-2*f(t)+f(t+h))/h^2.  For the endpoints (when closed=false) the algorithm
-//   when sufficient points are available the method is either the four point expression
-//   f''(t) = (2*f(t) - 5*f(t+h) + 4*f(t+2*h) - f(t+3*h))/h^2 or if five points are available
+//   f''(t) = (f(t-h)-2*f(t)+f(t+h))/h^2.  For the endpoints (when closed=false),
+//   when sufficient points are available, the method is either the four point expression
+//   f''(t) = (2*f(t) - 5*f(t+h) + 4*f(t+2*h) - f(t+3*h))/h^2 or 
 //   f''(t) = (35*f(t) - 104*f(t+h) + 114*f(t+2*h) - 56*f(t+3*h) + 11*f(t+4*h)) / 12h^2
+//   if five points are available.
+// 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 deriv2(data, h=1, closed=false) =
+    assert( is_consistent(data) , "Input list is not consistent or not numerical.") 
+    assert( len(data)>=3, "Input list has less than 3 elements.") 
+    assert( is_finite(h), "The sampling `h` must be a number." )
     let( L = len(data) )
     closed? [
         for(i=[0:1:L-1])
@@ -1003,16 +1149,19 @@ function deriv2(data, h=1, closed=false) =
 //   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. 
 //   If the `closed` parameter is true the data is assumed to be defined on a loop with data[0] adjacent to
-//   data[len(data)-1].  This function uses a five point derivative estimate, so the input must include five points:
+//   data[len(data)-1].  This function uses a five point derivative estimate, so the input data must include 
+//   at least five points:
 //   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.
 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.") 
+    assert( is_finite(h), "The sampling `h` must be a number." )
     let(
         L = len(data),
         h3 = h*h*h
     )
-    assert(L>=5, "Need five points for 3rd derivative estimate")
     closed? [
         for(i=[0:1:L-1])
         (-data[(L+i-2)%L]+2*data[(L+i-1)%L]-2*data[(i+1)%L]+data[(i+2)%L])/2/h3
@@ -1036,34 +1185,61 @@ function deriv3(data, h=1, closed=false) =
 // Function: C_times()
 // Usage: C_times(z1,z2)
 // Description:
-//   Multiplies two complex numbers.  
-function C_times(z1,z2) = [z1.x*z2.x-z1.y*z2.y,z1.x*z2.y+z1.y*z2.x];
+//   Multiplies two complex numbers represented by 2D vectors.  
+function C_times(z1,z2) = 
+    assert( is_vector(z1+z2,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)
 // Description:
-//   Divides z1 by z2.  
-function C_div(z1,z2) = let(den = z2.x*z2.x + z2.y*z2.y)
-   [(z1.x*z2.x + z1.y*z2.y)/den, (z1.y*z2.x-z1.x*z2.y)/den];
+//   Divides two complex numbers represented by 2D vectors.  
+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." ) 
+    let(den = z2.x*z2.x + z2.y*z2.y)
+    [(z1.x*z2.x + z1.y*z2.y)/den, (z1.y*z2.x - z1.x*z2.y)/den];
 
+// For the sake of consistence with Q_mul and vmul, C_times should be called C_mul
 
 // Section: Polynomials
 
-// Function: polynomial()
+// Function: polynomial() 
 // Usage:
 //   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]
 //   where a_n is the z^n coefficient.  Polynomial coefficients are real.
+//   The result is a number if `z` is a number and a complex number otherwise.
 
 // Note: this should probably be recoded to use division by [1,-z], which is more accurate
 // and avoids overflow with large coefficients, but requires poly_div to support complex coefficients.  
-function polynomial(p, z, k, zk, total) =
-   is_undef(k) ? polynomial(p, z, len(p)-1, is_num(z)? 1 : [1,0], is_num(z) ? 0 : [0,0]) :
-   k==-1 ? total :
-   polynomial(p, z, k-1, is_num(z) ? zk*z : C_times(zk,z), total+zk*p[k]);
+function polynomial(p, z, _k, _zk, _total) =
+    is_undef(_k)  
+    ?   assert( is_vector(p), "Input polynomial coefficients must be a vector." )
+        let(p = _poly_trim(p))
+        assert( is_finite(z) || is_vector(z,2), "The value of `z` must be a real or a complex number." )
+        polynomial( p, 
+                    z, 
+                    len(p)-1, 
+                    is_num(z)? 1 : [1,0], 
+                    is_num(z) ? 0 : [0,0]) 
+    :   _k==0 
+        ? _total + +_zk*p[0]
+        : polynomial( p, 
+                      z, 
+                      _k-1, 
+                      is_num(z) ? _zk*z : C_times(_zk,z), 
+                      _total+_zk*p[_k]);
 
+function polynomial(p,z,k,total) =
+     is_undef(k)
+   ?    assert( is_vector(p) , "Input polynomial coefficients must be a vector." )
+        assert( is_finite(z) || is_vector(z,2), "The value of `z` must be a real or a complex number." )
+        polynomial( _poly_trim(p), z, 0, is_num(z) ? 0 : [0,0])
+   : k==len(p) ? total
+   : polynomial(p,z,k+1, is_num(z) ? total*z+p[k] : C_times(total,z)+[p[k],0]);
 
 // Function: poly_mult()
 // Usage:
@@ -1073,21 +1249,37 @@ function polynomial(p, z, k, zk, total) =
 //   Given a list of polynomials represented as real coefficient lists, with the highest degree coefficient first, 
 //   computes the coefficient list of the product polynomial.  
 function poly_mult(p,q) = 
-  is_undef(q) ?
-     assert(is_list(p) && (is_vector(p[0]) || p[0]==[]), "Invalid arguments to poly_mult")
-     len(p)==2 ? poly_mult(p[0],p[1]) 
-               : poly_mult(p[0], poly_mult(select(p,1,-1)))
-  :
-  _poly_trim(
-  [
-  for(n = [len(p)+len(q)-2:-1:0])
-      sum( [for(i=[0:1:len(p)-1])
-           let(j = len(p)+len(q)- 2 - n - i)
-           if (j>=0 && j<len(q)) p[i]*q[j]
-               ])
-   ]);        
-
+    is_undef(q) ?
+       assert( is_list(p) 
+               && []==[for(pi=p) if( !is_vector(pi) && pi!=[]) 0], 
+               "Invalid arguments to poly_mult")
+       len(p)==2 ? poly_mult(p[0],p[1]) 
+                 : poly_mult(p[0], poly_mult(select(p,1,-1)))
+    :
+    _poly_trim(
+    [
+    for(n = [len(p)+len(q)-2:-1:0])
+        sum( [for(i=[0:1:len(p)-1])
+             let(j = len(p)+len(q)- 2 - n - i)
+             if (j>=0 && j<len(q)) p[i]*q[j]
+                 ])
+     ]);        
+     
+function poly_mult(p,q) = 
+    is_undef(q) ?
+       len(p)==2 ? poly_mult(p[0],p[1]) 
+                 : poly_mult(p[0], poly_mult(select(p,1,-1)))
+    :
+    assert( is_vector(p) && is_vector(q),"Invalid arguments to poly_mult")
+    _poly_trim( [
+                  for(n = [len(p)+len(q)-2:-1:0])
+                      sum( [for(i=[0:1:len(p)-1])
+                           let(j = len(p)+len(q)- 2 - n - i)
+                           if (j>=0 && j<len(q)) p[i]*q[j]
+                               ])
+                   ]);
 
+    
 // Function: poly_div()
 // Usage:
 //    [quotient,remainder] = poly_div(n,d)
@@ -1096,15 +1288,19 @@ function poly_mult(p,q) =
 //    a list of two polynomials, [quotient, remainder].  If the division has no remainder then
 //    the zero polynomial [] is returned for the remainder.  Similarly if the quotient is zero
 //    the returned quotient will be [].  
-function poly_div(n,d,q=[]) =
-    assert(len(d)>0 && d[0]!=0 , "Denominator is zero or has leading zero coefficient")
-    len(n)<len(d) ? [q,_poly_trim(n)] : 
-    let(
-      t = n[0] / d[0],
-      newq = concat(q,[t]),
-      newn =  [for(i=[1:1:len(n)-1]) i<len(d) ? n[i] - t*d[i] : n[i]]
-    )  
-    poly_div(newn,d,newq);
+function poly_div(n,d,q) =
+    is_undef(q) 
+    ?   assert( is_vector(n) && is_vector(d) , "Invalid polynomials." )
+        let( d = _poly_trim(d) )
+        assert( d!=[0] , "Denominator cannot be a zero polynomial." )
+        poly_div(n,d,q=[])
+    :   len(n)<len(d) ? [q,_poly_trim(n)] : 
+        let(
+          t = n[0] / d[0], 
+          newq = concat(q,[t]),
+          newn = [for(i=[1:1:len(n)-1]) i<len(d) ? n[i] - t*d[i] : n[i]]
+        )  
+        poly_div(newn,d,newq);
 
 
 // Internal Function: _poly_trim()
@@ -1114,8 +1310,8 @@ function poly_div(n,d,q=[]) =
 //    Removes leading zero terms of a polynomial.  By default zeros must be exact,
 //    or give epsilon for approximate zeros.  
 function _poly_trim(p,eps=0) =
-  let(  nz = [for(i=[0:1:len(p)-1]) if (!approx(p[i],0,eps)) i])
-  len(nz)==0 ? [] : select(p,nz[0],-1);
+    let( nz = [for(i=[0:1:len(p)-1]) if ( !approx(p[i],0,eps)) i])
+    len(nz)==0 ? [0] : select(p,nz[0],-1);
 
 
 // Function: poly_add()
@@ -1124,12 +1320,13 @@ function _poly_trim(p,eps=0) =
 // Description:
 //    Computes the sum of two polynomials.  
 function poly_add(p,q) = 
-  let(  plen = len(p),
-        qlen = len(q),
-        long = plen>qlen ? p : q,
-        short = plen>qlen ? q : p
-     )
-   _poly_trim(long + concat(repeat(0,len(long)-len(short)),short));
+    assert( is_vector(p) && is_vector(q), "Invalid input polynomial(s)." )
+    let(  plen = len(p),
+          qlen = len(q),
+          long = plen>qlen ? p : q,
+          short = plen>qlen ? q : p
+       )
+     _poly_trim(long + concat(repeat(0,len(long)-len(short)),short));
 
 
 // Function: poly_roots()
@@ -1150,38 +1347,38 @@ function poly_add(p,q) =
 //
 // Dario Bini. "Numerical computation of polynomial zeros by means of Aberth's Method", Numerical Algorithms, Feb 1996.
 // https://www.researchgate.net/publication/225654837_Numerical_computation_of_polynomial_zeros_by_means_of_Aberth's_method
-
 function poly_roots(p,tol=1e-14,error_bound=false) =
-  assert(p!=[], "Input polynomial must have a nonzero coefficient")
-  assert(is_vector(p), "Input must be a vector")
-  p[0] == 0 ? poly_roots(slice(p,1,-1),tol=tol,error_bound=error_bound) :    // Strip leading zero coefficients
-  p[len(p)-1] == 0 ?                                       // Strip trailing zero coefficients
-      let( solutions = poly_roots(select(p,0,-2),tol=tol, error_bound=error_bound))
-      (error_bound ? [ [[0,0], each solutions[0]], [0, each solutions[1]]]
-                  : [[0,0], each solutions]) :
-  len(p)==1 ? (error_bound ? [[],[]] : []) :               // Nonzero constant case has no solutions
-  len(p)==2 ? let( solution = [[-p[1]/p[0],0]])            // Linear case needs special handling
-              (error_bound ? [solution,[0]] : solution)
-  : 
-  let(
-      n = len(p)-1,   // polynomial degree
-      pderiv = [for(i=[0:n-1]) p[i]*(n-i)],
-         
-      s = [for(i=[0:1:n]) abs(p[i])*(4*(n-i)+1)],  // Error bound polynomial from Bini
+    assert( is_vector(p), "Invalid polynomial." )
+    let( p = _poly_trim(p,eps=0) )
+    assert( p!=[0], "Input polynomial cannot be zero." )
+    p[len(p)-1] == 0 ?                                       // Strip trailing zero coefficients
+        let( solutions = poly_roots(select(p,0,-2),tol=tol, error_bound=error_bound))
+        (error_bound ? [ [[0,0], each solutions[0]], [0, each solutions[1]]]
+                    : [[0,0], each solutions]) :
+    len(p)==1 ? (error_bound ? [[],[]] : []) :               // Nonzero constant case has no solutions
+    len(p)==2 ? let( solution = [[-p[1]/p[0],0]])            // Linear case needs special handling
+                (error_bound ? [solution,[0]] : solution)
+    : 
+    let(
+        n = len(p)-1,   // polynomial degree
+        pderiv = [for(i=[0:n-1]) p[i]*(n-i)],
+           
+        s = [for(i=[0:1:n]) abs(p[i])*(4*(n-i)+1)],  // Error bound polynomial from Bini
 
-      // Using method from: http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/0915-24.pdf
-      beta = -p[1]/p[0]/n,
-      r = 1+pow(abs(polynomial(p,beta)/p[0]),1/n),
-      init = [for(i=[0:1:n-1])                // Initial guess for roots       
-               let(angle = 360*i/n+270/n/PI)
-               [beta,0]+r*[cos(angle),sin(angle)]
-             ],
-      roots = _poly_roots(p,pderiv,s,init,tol=tol),
-      error = error_bound ? [for(xi=roots) n * (norm(polynomial(p,xi))+tol*polynomial(s,norm(xi))) /
-                                abs(norm(polynomial(pderiv,xi))-tol*polynomial(s,norm(xi)))] : 0
-    )
-    error_bound ? [roots, error] : roots;
+        // Using method from: http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/0915-24.pdf
+        beta = -p[1]/p[0]/n,
+        r = 1+pow(abs(polynomial(p,beta)/p[0]),1/n),
+        init = [for(i=[0:1:n-1])                // Initial guess for roots       
+                 let(angle = 360*i/n+270/n/PI)
+                 [beta,0]+r*[cos(angle),sin(angle)]
+               ],
+        roots = _poly_roots(p,pderiv,s,init,tol=tol),
+        error = error_bound ? [for(xi=roots) n * (norm(polynomial(p,xi))+tol*polynomial(s,norm(xi))) /
+                                  abs(norm(polynomial(pderiv,xi))-tol*polynomial(s,norm(xi)))] : 0
+      )
+      error_bound ? [roots, error] : roots;
 
+// Internal function
 // p = polynomial
 // pderiv = derivative polynomial of p
 // z = current guess for the roots
@@ -1222,12 +1419,16 @@ function _poly_roots(p, pderiv, s, z, tol, i=0) =
 //   tol = tolerance for the complex polynomial root finder
 
 function real_roots(p,eps=undef,tol=1e-14) =
-   let( 
+    assert( is_vector(p), "Invalid polynomial." )
+    let( p = _poly_trim(p,eps=0) )
+    assert( p!=[0], "Input polynomial cannot be zero." )
+    let( 
        roots_err = poly_roots(p,error_bound=true),
        roots = roots_err[0],
        err = roots_err[1]
-   )
-   is_def(eps) ? [for(z=roots) if (abs(z.y)/(1+norm(z))<eps) z.x]
-               : [for(i=idx(roots)) if (abs(roots[i].y)<=err[i]) roots[i].x];
+    )
+    is_def(eps) 
+    ? [for(z=roots) if (abs(z.y)/(1+norm(z))<eps) z.x]
+    : [for(i=idx(roots)) if (abs(roots[i].y)<=err[i]) roots[i].x];
 
 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
diff --git a/tests/polyhedra.scad b/tests/polyhedra.scad
index 3a77fbc..8c66f2e 100644
--- a/tests/polyhedra.scad
+++ b/tests/polyhedra.scad
@@ -2,10 +2,10 @@ include<../std.scad>
 include<../polyhedra.scad>
 
 
-$fn=96;
-
 if (true) {
 
+   $fn=96;
+
   // Display of all solids with insphere, midsphere and circumsphere
     
   for(i=[0:len(_polyhedra_)-1]) {
diff --git a/tests/test_arrays.scad b/tests/test_arrays.scad
index 756e433..f621cd0 100644
--- a/tests/test_arrays.scad
+++ b/tests/test_arrays.scad
@@ -1,37 +1,14 @@
 include <../std.scad>
 
-// List/Array Ops
 
-module test_repeat() {
-    assert(repeat(1, 4) == [1,1,1,1]);
-    assert(repeat(8, [2,3]) == [[8,8,8], [8,8,8]]);
-    assert(repeat(0, [2,2,3]) == [[[0,0,0],[0,0,0]], [[0,0,0],[0,0,0]]]);
-    assert(repeat([1,2,3],3) == [[1,2,3], [1,2,3], [1,2,3]]);
-}
-test_repeat();
+// Section: List Query Operations
 
-
-module test_in_list() {
-    assert(in_list("bar", ["foo", "bar", "baz"]));
-    assert(!in_list("bee", ["foo", "bar", "baz"]));
-    assert(in_list("bar", [[2,"foo"], [4,"bar"], [3,"baz"]], idx=1));
-    assert(!in_list(undef, [3,4,5]));
-    assert(in_list(undef,[3,4,undef,5]));
-    assert(!in_list(3,[]));
-    assert(!in_list(3,[4,5,[3]]));
-    
-}
-test_in_list();
-
-
-module test_slice() {
-    assert(slice([3,4,5,6,7,8,9], 3, 5) == [6,7]);
-    assert(slice([3,4,5,6,7,8,9], 2, -1) == [5,6,7,8,9]);
-    assert(slice([3,4,5,6,7,8,9], 1, 1) == []);
-    assert(slice([3,4,5,6,7,8,9], 6, -1) == [9]);
-    assert(slice([3,4,5,6,7,8,9], 2, -2) == [5,6,7,8]);
-}
-test_slice();
+module test_is_simple_list() {
+		assert(is_simple_list([1,2,3,4]));
+		assert(is_simple_list([]));
+		assert(!is_simple_list([1,2,[3,4]]));
+} 
+test_is_simple_list();
 
 
 module test_select() {
@@ -49,6 +26,74 @@ module test_select() {
 test_select();
 
 
+module test_slice() {
+    assert(slice([3,4,5,6,7,8,9], 3, 5) == [6,7]);
+    assert(slice([3,4,5,6,7,8,9], 2, -1) == [5,6,7,8,9]);
+    assert(slice([3,4,5,6,7,8,9], 1, 1) == []);
+    assert(slice([3,4,5,6,7,8,9], 6, -1) == [9]);
+    assert(slice([3,4,5,6,7,8,9], 2, -2) == [5,6,7,8]);
+    assert(slice([], 2, -2) == []);
+}
+test_slice();
+
+
+module test_in_list() {
+    assert(in_list("bar", ["foo", "bar", "baz"]));
+    assert(!in_list("bee", ["foo", "bar", "baz"]));
+    assert(in_list("bar", [[2,"foo"], [4,"bar"], [3,"baz"]], idx=1));
+    assert(!in_list("bee", ["foo", "bar", ["bee"]]));
+    assert(in_list(NAN, [NAN])==false);
+    assert(!in_list(undef, [3,4,5]));
+    assert(in_list(undef,[3,4,undef,5]));
+    assert(!in_list(3,[]));
+    assert(!in_list(3,[4,5,[3]]));
+}
+test_in_list();
+
+
+module test_min_index() {
+    assert(min_index([5,3,9,6,2,7,8,2,1])==8);
+    assert(min_index([5,3,9,6,2,7,8,2,7],all=true)==[4,7]);
+//    assert(min_index([],all=true)==[]);
+}
+test_min_index();
+
+
+module test_max_index() {
+    assert(max_index([5,3,9,6,2,7,8,9,1])==2);
+    assert(max_index([5,3,9,6,2,7,8,9,7],all=true)==[2,7]);
+//    assert(max_index([],all=true)==[]);
+}
+test_max_index();
+
+
+module test_list_increasing() {
+    assert(list_increasing([1,2,3,4]) == true);
+    assert(list_increasing([1,3,2,4]) == false);
+    assert(list_increasing([4,3,2,1]) == false);
+}
+test_list_increasing();
+
+
+module test_list_decreasing() {
+    assert(list_decreasing([1,2,3,4]) == false);
+    assert(list_decreasing([4,2,3,1]) == false);
+    assert(list_decreasing([4,3,2,1]) == true);
+}
+test_list_decreasing();
+
+// Section: Basic List Generation
+
+module test_repeat() {
+    assert(repeat(1, 4) == [1,1,1,1]);
+    assert(repeat(8, [2,3]) == [[8,8,8], [8,8,8]]);
+    assert(repeat(0, [2,2,3]) == [[[0,0,0],[0,0,0]], [[0,0,0],[0,0,0]]]);
+    assert(repeat([1,2,3],3) == [[1,2,3], [1,2,3], [1,2,3]]);
+    assert(repeat(4, [2,-1]) == [[], []]);
+}
+test_repeat();
+
+
 module test_list_range() {
     assert(list_range(4) == [0,1,2,3]);
     assert(list_range(n=4, step=2) == [0,2,4,6]);
@@ -66,6 +111,8 @@ test_list_range();
 
 module test_reverse() {
     assert(reverse([3,4,5,6]) == [6,5,4,3]);
+    assert(reverse("abcd") == ["d","c","b","a"]);
+    assert(reverse([]) == []);
 }
 test_reverse();
 
@@ -90,6 +137,8 @@ module test_deduplicate() {
     assert(deduplicate(closed=true, [8,3,4,4,4,8,2,3,3,8,8]) == [8,3,4,8,2,3]);
     assert(deduplicate("Hello") == ["H","e","l","o"]);
     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]]]]);
 }
 test_deduplicate();
 
@@ -148,22 +197,6 @@ module test_list_bset() {
 test_list_bset();
 
 
-module test_list_increasing() {
-    assert(list_increasing([1,2,3,4]) == true);
-    assert(list_increasing([1,3,2,4]) == false);
-    assert(list_increasing([4,3,2,1]) == false);
-}
-test_list_increasing();
-
-
-module test_list_decreasing() {
-    assert(list_decreasing([1,2,3,4]) == false);
-    assert(list_decreasing([4,2,3,1]) == false);
-    assert(list_decreasing([4,3,2,1]) == true);
-}
-test_list_decreasing();
-
-
 module test_list_shortest() {
     assert(list_shortest(["foobar", "bazquxx", "abcd"]) == 4);
 }
@@ -315,6 +348,13 @@ test_set_intersection();
 // Arrays
 
 
+module test_add_scalar() {
+    assert(add_scalar([1,2,3],3) == [4,5,6]);
+    assert(add_scalar([[1,2,3],[3,4,5]],3) == [[4,5,6],[6,7,8]]);
+}
+test_add_scalar();
+
+
 module test_subindex() {
     v = [[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]];
     assert(subindex(v,2) == [3, 7, 11, 15]);
@@ -402,10 +442,18 @@ test_array_group();
 
 module test_flatten() {
     assert(flatten([[1,2,3], [4,5,[6,7,8]]]) == [1,2,3,4,5,[6,7,8]]);
+    assert(flatten([]) == []);
 }
 test_flatten();
 
 
+module test_full_flatten() {
+    assert(full_flatten([[1,2,3], [4,5,[6,[7],8]]]) == [1,2,3,4,5,6,7,8]);
+    assert(full_flatten([]) == []);
+}
+test_full_flatten();
+
+
 module test_array_dim() {
     assert(array_dim([[[1,2,3],[4,5,6]],[[7,8,9],[10,11,12]]]) == [2,2,3]);
     assert(array_dim([[[1,2,3],[4,5,6]],[[7,8,9],[10,11,12]]], 0) == 2);
diff --git a/tests/test_common.scad b/tests/test_common.scad
index fa6bdbd..871eac2 100644
--- a/tests/test_common.scad
+++ b/tests/test_common.scad
@@ -18,6 +18,10 @@ module test_typeof() {
     assert(typeof([0:1:5]) == "range");
     assert(typeof([-3:2:5]) == "range");
     assert(typeof([10:-2:-10]) == "range");
+    assert(typeof([0:NAN:INF]) == "invalid");
+    assert(typeof([0:"a":INF]) == "undef"); 
+    assert(typeof([0:[]:INF]) == "undef"); 
+    assert(typeof([true:1:INF]) == "undef"); 
 }
 test_typeof();
 
@@ -102,6 +106,8 @@ module test_is_int() {
     assert(!is_int(-99.1));
     assert(!is_int(99.1));
     assert(!is_int(undef));
+    assert(!is_int(INF));
+    assert(!is_int(NAN));
     assert(!is_int(false));
     assert(!is_int(true));
     assert(!is_int("foo"));
@@ -124,6 +130,8 @@ module test_is_integer() {
     assert(!is_integer(-99.1));
     assert(!is_integer(99.1));
     assert(!is_integer(undef));
+    assert(!is_integer(INF));
+    assert(!is_integer(NAN));
     assert(!is_integer(false));
     assert(!is_integer(true));
     assert(!is_integer("foo"));
@@ -161,16 +169,33 @@ module test_is_range() {
     assert(!is_range(5));
     assert(!is_range(INF));
     assert(!is_range(-INF));
-    assert(!is_nan(NAN));
     assert(!is_range(""));
     assert(!is_range("foo"));
     assert(!is_range([]));
     assert(!is_range([3,4,5]));
+    assert(!is_range([INF:4:5]));
+    assert(!is_range([3:NAN:5]));
+    assert(!is_range([3:4:"a"]));
     assert(is_range([3:1:5]));
 }
-test_is_nan();
+test_is_range();
 
 
+module test_valid_range() {
+    assert(valid_range([0:0]));
+    assert(valid_range([0:1:0]));
+    assert(valid_range([0:1:10]));
+    assert(valid_range([0.1:1.1:2.1]));
+    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]));
+}
+test_valid_range();
+
 module test_is_list_of() {
     assert(is_list_of([3,4,5], 0));
     assert(!is_list_of([3,4,undef], 0));
@@ -181,10 +206,14 @@ module test_is_list_of() {
 }
 test_is_list_of();
 
-
 module test_is_consistent() {
+    assert(is_consistent([]));
+    assert(is_consistent([[],[]]));
     assert(is_consistent([3,4,5]));
     assert(is_consistent([[3,4],[4,5],[6,7]]));
+    assert(is_consistent([[[3],4],[[4],5]]));
+    assert(!is_consistent(5));
+    assert(!is_consistent(undef));
     assert(!is_consistent([[3,4,5],[3,4]]));
     assert(is_consistent([[3,[3,4,[5]]], [5,[2,9,[9]]]]));
     assert(!is_consistent([[3,[3,4,[5]]], [5,[2,9,9]]]));
@@ -331,11 +360,25 @@ module test_scalar_vec3() {
     assert(scalar_vec3([3]) == [3,0,0]);
     assert(scalar_vec3([3,4]) == [3,4,0]);
     assert(scalar_vec3([3,4],dflt=1) == [3,4,1]);
+    assert(scalar_vec3([3,"a"],dflt=1) == [3,"a",1]);
+    assert(scalar_vec3([3,[2]],dflt=1) == [3,[2],1]);
     assert(scalar_vec3([3],dflt=1) == [3,1,1]);
     assert(scalar_vec3([3,4,5]) == [3,4,5]);
     assert(scalar_vec3([3,4,5,6]) == [3,4,5]);
+    assert(scalar_vec3([3,4,5,6]) == [3,4,5]);
 }
 test_scalar_vec3();
 
 
+module test_segs() {
+    assert_equal(segs(50,$fn=8), 8);
+    assert_equal(segs(50,$fa=2,$fs=2), 158);
+    assert(segs(1)==5);
+    assert(segs(11)==30);
+  //  assert(segs(1/0)==5);
+  //  assert(segs(0/0)==5);
+  //  assert(segs(undef)==5);
+}
+test_segs();
+
 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
diff --git a/tests/test_math.scad b/tests/test_math.scad
index 3233e6f..f1b36b9 100644
--- a/tests/test_math.scad
+++ b/tests/test_math.scad
@@ -110,6 +110,8 @@ module test_approx() {
     assert_equal(approx(1/3, 0.3333333333), true);
     assert_equal(approx(-1/3, -0.3333333333), true);
     assert_equal(approx(10*[cos(30),sin(30)], 10*[sqrt(3)/2, 1/2]), true);
+    assert_equal(approx([1,[1,undef]], [1+1e-12,[1,true]]), false);
+    assert_equal(approx([1,[1,undef]], [1+1e-12,[1,undef]]), true);
 }
 test_approx();
 
@@ -389,7 +391,6 @@ module test_mean() {
 }
 test_mean();
 
-
 module test_median() {
     assert_equal(median([2,3,7]), 4.5);
     assert_equal(median([[1,2,3], [3,4,5], [8,9,10]]), [4.5,5.5,6.5]);
@@ -397,6 +398,16 @@ module test_median() {
 test_median();
 
 
+module test_convolve() {
+    assert_equal(convolve([],[1,2,1]), []);
+    assert_equal(convolve([1,1],[]), []);
+    assert_equal(convolve([1,1],[1,2,1]), [1,3,3,1]);
+    assert_equal(convolve([1,2,3],[1,2,1]), [1,4,8,8,3]);
+}
+test_convolve();
+
+
+
 module test_matrix_inverse() {
     assert_approx(matrix_inverse(rot([20,30,40])), [[0.663413948169,0.556670399226,-0.5,0],[-0.47302145844,0.829769465589,0.296198132726,0],[0.579769465589,0.0400087565481,0.813797681349,0],[0,0,0,1]]);
 }
@@ -583,6 +594,24 @@ module test_factorial() {
 }
 test_factorial();
 
+module test_binomial() {
+    assert_equal(binomial(1), [1,1]);
+    assert_equal(binomial(2), [1,2,1]);
+    assert_equal(binomial(3), [1,3,3,1]);
+    assert_equal(binomial(5), [1,5,10,10,5,1]);
+}
+test_binomial();
+
+module test_binomial_coefficient() {
+    assert_equal(binomial_coefficient(2,1), 2);
+    assert_equal(binomial_coefficient(3,2), 3);
+    assert_equal(binomial_coefficient(4,2), 6);
+    assert_equal(binomial_coefficient(10,7), 120);
+    assert_equal(binomial_coefficient(10,7), binomial(10)[7]);
+    assert_equal(binomial_coefficient(15,4), binomial(15)[4]);
+}
+test_binomial_coefficient();
+
 
 module test_gcd() {
     assert_equal(gcd(15,25), 5);
@@ -682,6 +711,7 @@ test_linear_solve();
 
 module test_outer_product(){
   assert_equal(outer_product([1,2,3],[4,5,6]), [[4,5,6],[8,10,12],[12,15,18]]);
+  assert_equal(outer_product([1,2],[4,5,6]), [[4,5,6],[8,10,12]]);
   assert_equal(outer_product([9],[7]), [[63]]);
 }
 test_outer_product();
@@ -782,8 +812,10 @@ test_deriv3();
 
 
 module test_polynomial(){
-  assert_equal(polynomial([],12),0);
-  assert_equal(polynomial([],[12,4]),[0,0]);
+  assert_equal(polynomial([0],12),0);
+  assert_equal(polynomial([0],[12,4]),[0,0]);
+//  assert_equal(polynomial([],12),0);
+//  assert_equal(polynomial([],[12,4]),[0,0]);
   assert_equal(polynomial([1,2,3,4],3),58);
   assert_equal(polynomial([1,2,3,4],[3,-1]),[47,-41]);
   assert_equal(polynomial([0,0,2],4),2);
@@ -879,16 +911,20 @@ test_qr_factor();
 
 module test_poly_mult(){
   assert_equal(poly_mult([3,2,1],[4,5,6,7]),[12,23,32,38,20,7]);
-  assert_equal(poly_mult([3,2,1],[]),[]);
+  assert_equal(poly_mult([3,2,1],[0]),[0]);
+//  assert_equal(poly_mult([3,2,1],[]),[]);
   assert_equal(poly_mult([[1,2],[3,4],[5,6]]), [15,68,100,48]);
-  assert_equal(poly_mult([[1,2],[],[5,6]]), []);
-  assert_equal(poly_mult([[3,4,5],[0,0,0]]),[]);
+  assert_equal(poly_mult([[1,2],[0],[5,6]]), [0]);
+//  assert_equal(poly_mult([[1,2],[],[5,6]]), []);
+  assert_equal(poly_mult([[3,4,5],[0,0,0]]),[0]);
+//  assert_equal(poly_mult([[3,4,5],[0,0,0]]),[]);
 }
 test_poly_mult();
 
- 
+
 module test_poly_div(){
-  assert_equal(poly_div(poly_mult([4,3,3,2],[2,1,3]), [2,1,3]),[[4,3,3,2],[]]);
+  assert_equal(poly_div(poly_mult([4,3,3,2],[2,1,3]), [2,1,3]),[[4,3,3,2],[0]]);
+//  assert_equal(poly_div(poly_mult([4,3,3,2],[2,1,3]), [2,1,3]),[[4,3,3,2],[]]);
   assert_equal(poly_div([1,2,3,4],[1,2,3,4,5]), [[], [1,2,3,4]]);
   assert_equal(poly_div(poly_add(poly_mult([1,2,3,4],[2,0,2]), [1,1,2]), [1,2,3,4]), [[2,0,2],[1,1,2]]);
   assert_equal(poly_div([1,2,3,4], [1,-3]), [[1,5,18],[58]]);
@@ -899,7 +935,8 @@ test_poly_div();
 module test_poly_add(){
   assert_equal(poly_add([2,3,4],[3,4,5,6]),[3,6,8,10]);
   assert_equal(poly_add([1,2,3,4],[-1,-2,3,4]), [6,8]);
-  assert_equal(poly_add([1,2,3],-[1,2,3]),[]);
+  assert_equal(poly_add([1,2,3],-[1,2,3]),[0]);
+//  assert_equal(poly_add([1,2,3],-[1,2,3]),[]);
 }
 test_poly_add();
 
diff --git a/tests/test_vectors.scad b/tests/test_vectors.scad
index 42f0a04..743112f 100644
--- a/tests/test_vectors.scad
+++ b/tests/test_vectors.scad
@@ -9,6 +9,7 @@ module test_is_vector() {
     assert(is_vector(1) == false);
     assert(is_vector("foo") == false);
     assert(is_vector(true) == false);
+
     assert(is_vector([0,0,0],zero=true) == true);
     assert(is_vector([0,0,0],zero=false) == false);
     assert(is_vector([0,1,0],zero=true) == false);
@@ -17,13 +18,6 @@ module test_is_vector() {
 test_is_vector();
 
 
-module test_add_scalar() {
-    assert(add_scalar([1,2,3],3) == [4,5,6]);
-    assert(add_scalar([[1,2,3],[3,4,5]],3) == [[4,5,6],[6,7,8]]);
-}
-test_add_scalar();
-
-
 module test_vfloor() {
     assert_equal(vfloor([2.0, 3.14, 18.9, 7]), [2,3,18,7]);
     assert_equal(vfloor([-2.0, -3.14, -18.9, -7]), [-2,-4,-19,-7]);
diff --git a/vectors.scad b/vectors.scad
index bb28ead..51f4735 100644
--- a/vectors.scad
+++ b/vectors.scad
@@ -19,42 +19,25 @@
 // Arguments:
 //   v = The value to test to see if it is a vector.
 //   length = If given, make sure the vector is `length` items long.
-//   zero = If false, require that the length of the vector is not approximately zero.  If true, require the length of the vector to be approx zero-length.  Default: `undef` (don't check vector length.)
-//   eps = The minimum vector length that is considered non-zero.  Default: `EPSILON` (`1e-9`)
 // Example:
-//   is_vector(4);              // Returns false
-//   is_vector([4,true,false]); // Returns false
-//   is_vector([3,4,INF,5]);    // Returns false
-//   is_vector([3,4,5,6]);      // Returns true
-//   is_vector([3,4,undef,5]);  // Returns false
-//   is_vector([3,4,5],3);      // Returns true
-//   is_vector([3,4,5],4);      // Returns true
-//   is_vector([]);             // Returns false
-//   is_vector([0,0,0],zero=true);   // Returns true
-//   is_vector([0,0,0],zero=false);  // Returns false
-//   is_vector([0,1,0],zero=true);   // Returns false
-//   is_vector([0,0,1],zero=false);  // Returns true
+//   is_vector(4);                          // Returns false
+//   is_vector([4,true,false]);             // Returns false
+//   is_vector([3,4,INF,5]);                // Returns false
+//   is_vector([3,4,5,6]);                  // Returns true
+//   is_vector([3,4,undef,5]);              // Returns false
+//   is_vector([3,4,5],3);                  // Returns true
+//   is_vector([3,4,5],4);                  // Returns true
+//   is_vector([]);                         // Returns false
+//   is_vector([0,4,0],3,zero=false);     // Returns true
+//   is_vector([0,0,0],zero=false);       // Returns false
+//   is_vector([0,0,1e-12],zero=false);   // Returns false
+//   is_vector([],zero=false);            // Returns false
 function is_vector(v,length,zero,eps=EPSILON) =
     is_list(v) && is_num(0*(v*v))
     && (is_undef(length) || len(v)==length)
     && (is_undef(zero) || ((norm(v) >= eps) == !zero));
 
 
-// Function: add_scalar()
-// Usage:
-//   add_scalar(v,s);
-// Description:
-//   Given a vector and a scalar, returns the vector with the scalar added to each item in it.
-//   If given a list of vectors, recursively adds the scalar to the each vector.
-// Arguments:
-//   v = The initial list of values.
-//   s = A scalar value to add to every item in the vector.
-// Example:
-//   add_scalar([1,2,3],3);            // Returns: [4,5,6]
-//   add_scalar([[1,2,3],[3,4,5]],3);  // Returns: [[4,5,6],[6,7,8]]
-function add_scalar(v,s) = [for (x=v) is_list(x)? add_scalar(x,s) : x+s];
-
-
 // Function: vang()
 // Usage:
 //   theta = vang([X,Y]);
@@ -63,6 +46,7 @@ function add_scalar(v,s) = [for (x=v) is_list(x)? add_scalar(x,s) : x+s];
 //   Given a 2D vector, returns the angle in degrees counter-clockwise from X+ on the XY plane.
 //   Given a 3D vector, returns [THETA,PHI] where THETA is the number of degrees counter-clockwise from X+ on the XY plane, and PHI is the number of degrees up from the X+ axis along the XZ plane.
 function vang(v) =
+    assert( is_vector(v,2) || is_vector(v,3) , "Invalid vector")
     len(v)==2? atan2(v.y,v.x) :
     let(res=xyz_to_spherical(v)) [res[1], 90-res[2]];
 
@@ -70,13 +54,19 @@ function vang(v) =
 // Function: vmul()
 // Description:
 //   Element-wise vector multiplication.  Multiplies each element of vector `v1` by
-//   the corresponding element of vector `v2`.  Returns a vector of the products.
+//   the corresponding element of vector `v2`. The vectors should have the same dimension.
+//   Returns a vector of the products.
 // Arguments:
 //   v1 = The first vector.
 //   v2 = The second vector.
 // Example:
 //   vmul([3,4,5], [8,7,6]);  // Returns [24, 28, 30]
-function vmul(v1, v2) = [for (i = [0:1:len(v1)-1]) v1[i]*v2[i]];
+function vmul(v1, v2) = 
+// this thighter check can be done yet because it would break other codes in the library
+//    assert( is_vector(v1) && is_vector(v2,len(v1)), "Incompatible vectors")
+    assert( is_vector(v1) && is_vector(v2), "Invalid vector(s)")
+    [for (i = [0:1:len(v1)-1]) v1[i]*v2[i]];
+    
 
 
 // Function: vdiv()
@@ -88,7 +78,9 @@ function vmul(v1, v2) = [for (i = [0:1:len(v1)-1]) v1[i]*v2[i]];
 //   v2 = The second vector.
 // Example:
 //   vdiv([24,28,30], [8,7,6]);  // Returns [3, 4, 5]
-function vdiv(v1, v2) = [for (i = [0:1:len(v1)-1]) v1[i]/v2[i]];
+function vdiv(v1, v2) = 
+    assert( is_vector(v1) && is_vector(v2,len(v1)), "Incompatible vectors")
+    [for (i = [0:1:len(v1)-1]) v1[i]/v2[i]];
 
 
 // Function: vabs()
@@ -97,19 +89,25 @@ function vdiv(v1, v2) = [for (i = [0:1:len(v1)-1]) v1[i]/v2[i]];
 //   v = The vector to get the absolute values of.
 // Example:
 //   vabs([-1,3,-9]);  // Returns: [1,3,9]
-function vabs(v) = [for (x=v) abs(x)];
+function vabs(v) =
+    assert( is_vector(v), "Invalid vector" ) 
+    [for (x=v) abs(x)];
 
 
 // Function: vfloor()
 // Description:
 //   Returns the given vector after performing a `floor()` on all items.
-function vfloor(v) = [for (x=v) floor(x)];
+function vfloor(v) =
+    assert( is_vector(v), "Invalid vector" ) 
+    [for (x=v) floor(x)];
 
 
 // Function: vceil()
 // Description:
 //   Returns the given vector after performing a `ceil()` on all items.
-function vceil(v) = [for (x=v) ceil(x)];
+function vceil(v) =
+    assert( is_vector(v), "Invalid vector" ) 
+    [for (x=v) ceil(x)];
 
 
 // Function: unit()
@@ -137,6 +135,7 @@ function unit(v, error=[[["ASSERT"]]]) =
 // Function: vector_angle()
 // Usage:
 //   vector_angle(v1,v2);
+//   vector_angle([v1,v2]);
 //   vector_angle(PT1,PT2,PT3);
 //   vector_angle([PT1,PT2,PT3]);
 // Description:
@@ -156,34 +155,36 @@ function unit(v, error=[[["ASSERT"]]]) =
 //   vector_angle([10,0,10], [0,0,0], [-10,10,0]);  // Returns: 120
 //   vector_angle([[10,0,10], [0,0,0], [-10,10,0]]);  // Returns: 120
 function vector_angle(v1,v2,v3) =
-    let(
-        vecs = !is_undef(v3)? [v1-v2,v3-v2] :
-            !is_undef(v2)? [v1,v2] :
-            len(v1) == 3? [v1[0]-v1[1],v1[2]-v1[1]] :
-            len(v1) == 2? v1 :
-            assert(false, "Bad arguments to vector_angle()"),
-        is_valid = is_vector(vecs[0]) && is_vector(vecs[1]) && vecs[0]*0 == vecs[1]*0
+    assert( ( is_undef(v3) && ( is_undef(v2) || same_shape(v1,v2) ) )
+            || is_consistent([v1,v2,v3]) ,
+            "Bad arguments.")
+    assert( is_vector(v1) || is_consistent(v1), "Bad arguments.") 
+    let( vecs = ! is_undef(v3) ? [v1-v2,v3-v2] :
+                ! is_undef(v2) ? [v1,v2] :
+                len(v1) == 3   ? [v1[0]-v1[1], v1[2]-v1[1]] 
+                               : v1
     )
-    assert(is_valid, "Bad arguments to vector_angle()")
+    assert(is_vector(vecs[0],2) || is_vector(vecs[0],3), "Bad arguments.")
     let(
         norm0 = norm(vecs[0]),
         norm1 = norm(vecs[1])
     )
-    assert(norm0>0 && norm1>0,"Zero length vector given to vector_angle()")
+    assert(norm0>0 && norm1>0, "Zero length vector.")
     // NOTE: constrain() corrects crazy FP rounding errors that exceed acos()'s domain.
     acos(constrain((vecs[0]*vecs[1])/(norm0*norm1), -1, 1));
-
+    
 
 // Function: vector_axis()
 // Usage:
 //   vector_axis(v1,v2);
+//   vector_axis([v1,v2]);
 //   vector_axis(PT1,PT2,PT3);
 //   vector_axis([PT1,PT2,PT3]);
 // Description:
 //   If given a single list of two vectors, like `vector_axis([V1,V2])`, returns the vector perpendicular the two vectors V1 and V2.
-//   If given a single list of three points, like `vector_axis([A,B,C])`, returns the vector perpendicular the line segments AB and BC.
-//   If given two vectors, like `vector_axis(V1,V1)`, returns the vector perpendicular the two vectors V1 and V2.
-//   If given three points, like `vector_axis(A,B,C)`, returns the vector perpendicular the line segments AB and BC.
+//   If given a single list of three points, like `vector_axis([A,B,C])`, returns the vector perpendicular to the plane through a, B and C.
+//   If given two vectors, like `vector_axis(V1,V2)`, returns the vector perpendicular to the two vectors V1 and V2.
+//   If given three points, like `vector_axis(A,B,C)`, returns the vector perpendicular to the plane through a, B and C.
 // Arguments:
 //   v1 = First vector or point.
 //   v2 = Second vector or point.
@@ -199,28 +200,23 @@ function vector_axis(v1,v2=undef,v3=undef) =
     is_vector(v3)
     ?   assert(is_consistent([v3,v2,v1]), "Bad arguments.")
         vector_axis(v1-v2, v3-v2)
-    :
-    assert( is_undef(v3), "Bad arguments.")
-    is_undef(v2)
-    ?   assert( is_list(v1), "Bad arguments.")
-        len(v1) == 2
-            ? vector_axis(v1[0],v1[1])
-            : vector_axis(v1[0],v1[1],v1[2])
-    :
-    assert(
-        is_vector(v1,zero=false) &&
-            is_vector(v2,zero=false) &&
-            is_consistent([v1,v2]),
-        "Bad arguments."
-    )
-    let(
-        eps = 1e-6,
-        w1 = point3d(v1/norm(v1)),
-        w2 = point3d(v2/norm(v2)),
-        w3 = (norm(w1-w2) > eps && norm(w1+w2) > eps) ? w2
-           : (norm(vabs(w2)-UP) > eps) ? UP
-           : RIGHT
-    ) unit(cross(w1,w3));
+    :   assert( is_undef(v3), "Bad arguments.")
+        is_undef(v2)
+        ?   assert( is_list(v1), "Bad arguments.")
+            len(v1) == 2 
+            ?   vector_axis(v1[0],v1[1]) 
+            :   vector_axis(v1[0],v1[1],v1[2])
+        :   assert( is_vector(v1,zero=false) && is_vector(v2,zero=false) && is_consistent([v1,v2])
+                    , "Bad arguments.")  
+            let(
+              eps = 1e-6,
+              w1 = point3d(v1/norm(v1)),
+              w2 = point3d(v2/norm(v2)),
+              w3 = (norm(w1-w2) > eps && norm(w1+w2) > eps) ? w2 
+                   : (norm(vabs(w2)-UP) > eps)? UP 
+                   : RIGHT
+            ) unit(cross(w1,w3));
+
 
 
 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap