diff --git a/arrays.scad b/arrays.scad
index 4f1c317..835fd65 100644
--- a/arrays.scad
+++ b/arrays.scad
@@ -18,20 +18,6 @@
 
 // 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:
@@ -73,9 +59,6 @@ function select(list, start, end=undef) =
             :   concat([for (i = [s:1:l-1]) list[i]], [for (i = [0:1:e]) list[i]]) ;
 
 
-
-
-
 // Function: slice()
 // Description:
 //   Returns a slice of a list.  The first item is index 0.
@@ -101,24 +84,23 @@ function slice(list,start,end) =
         ) [for (i=[s:1:e-1]) if (e>s) list[i]];
 
 
-
-
 // Function: in_list()
 // Description: Returns true if value `val` is in list `list`. When `val==NAN` the answer will be false for any list.
 // Arguments:
 //   val = The simple value to search for.
 //   list = The list to search.
-//   idx = If given, searches the given subindexes for matches for `val`.
+//   idx = If given, searches the given subindex 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) = 
+    assert( is_list(list) && (is_undef(idx) || is_finite(idx)),
+            "Invalid input." )
     let( s = search([val], list, num_returns_per_match=1, index_col_num=idx)[0] )
-    s==[] || s[0]==[] ? false
+    s==[] || s==[[]] ? false
     : is_undef(idx) ? val==list[s] 
     : val==list[s][idx];
-    
 
 
 // Function: min_index()
@@ -209,7 +191,6 @@ function repeat(val, n, i=0) =
     [for (j=[1:1:n[i]]) repeat(val, n, i+1)];
 
 
-
 // Function: list_range()
 // Usage:
 //   list_range(n, [s], [e])
@@ -249,7 +230,6 @@ function list_range(n=undef, s=0, e=undef, step=undef) =
     
 
 
-
 // Section: List Manipulation
 
 // Function: reverse()
@@ -315,8 +295,6 @@ function deduplicate(list, closed=false, eps=EPSILON) =
     : [for (i=[0:1:l-1]) if (i==end || !approx(list[i], list[(i+1)%l], eps)) list[i]];
 
 
-
-
 // Function: deduplicate_indexed()
 // Usage:
 //   new_idxs = deduplicate_indexed(list, indices, [closed], [eps]);
@@ -351,8 +329,6 @@ function deduplicate_indexed(list, indices, closed=false, eps=EPSILON) =
     ];
 
 
-
-
 // Function: repeat_entries()
 // Usage:
 //   newlist = repeat_entries(list, N)
@@ -390,8 +366,6 @@ function repeat_entries(list, N, exact = true) =
                : [for (val=reps_guess) round(val)]
     )
     [for(i=[0:length-1]) each repeat(list[i],reps[i])];
-    
-
 
 
 // Function: list_set()
@@ -431,7 +405,6 @@ function list_set(list=[],indices,values,dflt=0,minlen=0) =
                 dflt ,
           each repeat(dflt, minlen-max(indices))
         ];
-      
 
 
 // Function: list_insert()
@@ -465,8 +438,6 @@ function list_insert(list, indices, values, _i=0) =
         ];
 
 
-
-
 // Function: list_remove()
 // Usage:
 //   list_remove(list, indices)
@@ -494,8 +465,6 @@ function list_remove(list, indices) =
         ]; 
 
 
-
-
 // Function: list_remove_values()
 // Usage:
 //   list_remove_values(list,values,all=false) =
@@ -565,8 +534,6 @@ function list_bset(indexset, valuelist, dflt=0) =
     );
 
 
-
-
 // Section: List Length Manipulation
 
 // Function: list_shortest()
@@ -579,7 +546,6 @@ function list_shortest(array) =
     min([for (v = array) len(v)]);
 
 
-
 // Function: list_longest()
 // Description:
 //   Returns the length of the longest sublist in a list of lists.
@@ -629,7 +595,6 @@ function list_fit(array, length, fill) =
               : list_pad(array,length,fill);
 
 
-
 // Section: List Shuffling and Sorting
 
 // Function: shuffle()
@@ -684,6 +649,7 @@ function _sort_vectors2(arr) =
     ) 
     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) =
@@ -711,7 +677,6 @@ function _sort_vectors3(arr) =
     ) 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) =
@@ -742,45 +707,38 @@ function _sort_vectors4(arr) =
                                                     && y[3]>pivot[3] )))))) 
                     y ]
     ) concat( _sort_vectors4(lesser), equal, _sort_vectors4(greater) );
-    
 
 
+// when idx==undef, returns the sorted array
+// otherwise, returns the indices of the sorted array
 function _sort_general(arr, idx=undef) =
     (len(arr)<=1) ? arr :
+    is_undef(idx) 
+    ? _sort_scalar(arr)
+    : let( arrind=[for(k=[0:len(arr)-1], ark=[arr[k]]) [ k, [for (i=idx) ark[i]] ] ] )
+      _indexed_sort(arrind);
+      
+// given a list of pairs, return the first element of each pair of the list sorted by the second element of the pair
+// the sorting is done using compare_vals()
+function _indexed_sort(arrind) = 
+    arrind==[] ? [] : len(arrind)==1? [arrind[0][0]] :
+    let( pivot = arrind[floor(len(arrind)/2)][1] )
     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));
+        lesser  = [ for (entry=arrind) if (compare_vals(entry[1], pivot) <0 ) entry ],
+        equal   = [ for (entry=arrind) if (compare_vals(entry[1], pivot)==0 ) entry[0] ],
+        greater = [ for (entry=arrind) if (compare_vals(entry[1], pivot) >0 ) entry ]
+      )
+    concat(_indexed_sort(lesser), equal, _indexed_sort(greater));
+
+
+// returns true for valid index specifications idx in the interval [imin, imax) 
+// note that idx can't have any value greater or EQUAL to imax
+function _valid_idx(idx,imin,imax) =
+    is_undef(idx) 
+    || ( is_finite(idx) && idx>=imin && idx< imax )
+    || ( is_list(idx) && min(idx)>=imin && max(idx)< imax )
+    || ( valid_range(idx) && idx[0]>=imin && idx[2]< imax );
     
-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 = [
-            for (entry = arr) let(
-                val = idx==undef? entry : [for (i=idx) entry[i]],
-                cmp = compare_vals(val, pivotval)
-            ) cmp
-        ],
-        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()
 // Usage:
@@ -799,20 +757,21 @@ 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 
-    ? (
-        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);
-
+    is_def(idx) 
+    ?   assert( _valid_idx(idx,0,len(list)) , "Invalid indices.")
+        let( sarr = _sort_general(list,idx) )
+        [for(i=[0:len(sarr)-1]) list[sarr[i]] ] 
+    :   let(size = array_dim(list))
+        len(size)==1 ? _sort_scalars(list) :
+        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);
 
 
 // Function: sortidx()
@@ -838,13 +797,13 @@ function sort(list, idx=undef) =
 //   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.")
+    assert( _valid_idx(idx,0,len(list)) , "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)
+              :  0
     )
     is_undef(idx) && len(size) == 1? subindex(_sort_vectors1(aug),1) :
     is_undef(idx) && len(size) == 2 && size[1] <=4
@@ -856,25 +815,8 @@ function sortidx(list, idx=undef) =
     /*size[1]==4*/ : subindex(_sort_vectors4(aug),4)
       ) 
     :   // general case
-        subindex(_sort_general(aug, idx=list_range(s=1,n=len(aug)-1)), 0);
+        _sort_general(list,idx);
 
-function sortidx(list, idx=undef) =
-    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);
 
 // sort() does not accept strings but sortidx does; isn't inconsistent ?
 
@@ -896,7 +838,6 @@ function unique(arr) =
     ];
 
 
-
 // Function: unique_count()
 // Usage:
 //   unique_count(arr);
@@ -913,8 +854,6 @@ function unique_count(arr) =
       [ select(arr,ind), deltas( concat(ind,[len(arr)]) ) ];
 
 
-
-
 // Section: List Iteration Helpers
 
 // Function: idx()
@@ -952,10 +891,10 @@ function idx(list, step=1, end=-1,start=0) =
 //   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), "Invalid input." )
-    assert(is_undef(idx)||is_finite(idx)||is_vector(idx) ||is_range(idx), "Invalid index/indices." )
+    assert( _valid_idx(idx,0,len(l)), "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]])];
+    :   [for (i=[0:1:len(l)-1]) [ i, for (j=idx) l[i][j]] ];
 
 
 // Function: force_list()
@@ -1109,8 +1048,6 @@ function set_union(a, b, get_indices=false) =
     ) [idxs, nset];
 
 
-
-
 // Function: set_difference()
 // Usage:
 //   s = set_difference(a, b);
@@ -1130,7 +1067,6 @@ function set_difference(a, b) =
     [ for (i=idx(a)) if(found[i]==[]) a[i] ];
 
 
-
 // Function: set_intersection()
 // Usage:
 //   s = set_intersection(a, b);
@@ -1151,7 +1087,6 @@ function set_intersection(a, b) =
 
 
 
-
 // Section: Array Manipulation
 
 // Function: add_scalar()
@@ -1170,14 +1105,14 @@ 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()
 // Usage:
 //   subindex(M, idx)
 // Description:
 //   Extracts the entries listed in idx from each entry in M.  For a matrix this means
-//   selecting a specified set of columsn.  If idx is a number the return is a vector, otherwise
-//   it is a list of lists (the submatrix).  
+//   selecting a specified set of columns.  If idx is a number the return is a vector, 
+//   otherwise it is a list of lists (the submatrix).  
+//   This function will return `undef` at all entry positions indexed by idx not found in the input list M.
 // Arguments:
 //   M = The given list of lists.
 //   idx = The index, list of indices, or range of indices to fetch.
@@ -1187,11 +1122,45 @@ function add_scalar(v,s) =
 //   subindex(M,[2]);    // Returns [[3], [7], [11], [15]]
 //   subindex(M,[2,1]);  // Returns [[3, 2], [7, 6], [11, 10], [15, 14]]
 //   subindex(M,[1:3]);  // Returns [[2, 3, 4], [6, 7, 8], [10, 11, 12], [14, 15, 16]]
+//   N = [ [1,2], [3], [4,5], [6,7,8] ];
+//   subindex(N,[0,1]);  // Returns [ [1,2], [3,undef], [4,5], [6,7] ]
 function subindex(M, idx) =
-    is_num(idx)
+    assert( is_list(M), "The input is not a list." )
+    assert( !is_undef(idx) && _valid_idx(idx,0,1/0), "Invalid index input." ) 
+    is_finite(idx)
       ? [for(row=M) row[idx]]
       : [for(row=M) [for(i=idx) row[i]]];
 
+
+// Function: submatrix()
+// Usage: submatrix(M, idx1, idx2)
+// Description:
+//   The input must be a list of lists (a matrix or 2d array).  Returns a submatrix by selecting the rows listed in idx1 and columsn listed in idx2.
+// Arguments:
+//   M = Given list of lists
+//   idx1 = rows index list or range
+//   idx2 = column index list or range
+// Example:
+//   M = [[ 1, 2, 3, 4, 5],
+//        [ 6, 7, 8, 9,10],
+//        [11,12,13,14,15],
+//        [16,17,18,19,20],
+//        [21,22,23,24,25]];
+//   submatrix(M,[1:2],[3:4]);  // Returns [[9, 10], [14, 15]]
+//   submatrix(M,[1], [3,4]));  // Returns [[9,10]]
+//   submatrix(M,1, [3,4]));  // Returns [[9,10]]
+//   submatrix(M,1,3));  // Returns [[9]]
+//   submatrix(M, [3,4],1); // Returns  [[17],[22]]);
+//   submatrix(M, [1,3],[2,4]); // Returns [[8,10],[18,20]]);
+//   A = [[true,    17, "test"],
+//        [[4,2],   91, false],
+//        [6,    [3,4], undef]];
+//   submatrix(A,[0,2],[1,2]);   // Returns [[17, "test"], [[3, 4], undef]]
+
+function submatrix(M,idx1,idx2) =
+    [for(i=idx1) [for(j=idx2) M[i][j] ] ];
+
+
 // Function: zip()
 // Usage:
 //   zip(v1, v2, v3, [fit], [fill]);
@@ -1318,6 +1287,10 @@ function array_dim(v, depth=undef) =
 
 // Function: transpose()
 // Description: Returns the transposition of the given array.
+//    When reverse=true, the transposition is done in respect to the secondary diagonal, that is:
+//    .
+//    reverse(transpose(reverse(arr))) == transpose(arr, reverse=true)
+//    By default, reverse=false.
 // Example:
 //   arr = [
 //       ["a", "b", "c"],
@@ -1344,16 +1317,32 @@ function array_dim(v, depth=undef) =
 //   //     ["c", "f"],
 //   // ]
 // Example:
+//   arr = [
+//       ["a", "b", "c"],
+//       ["d", "e", "f"],
+//       ["g", "h", "i"]
+//   ];
+//   t = transpose(arr, reverse=true);
+//   // Returns:
+//   // [
+//   //  ["i", "f", "c"],
+//   //  ["h", "e", "b"],
+//   //  ["g", "d", "a"]
+//   // ]
+// Example:
 //   transpose([3,4,5]);  // Returns: [3,4,5]
-function transpose(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;
-
-
+function transpose(arr, reverse=false) =
+    assert( is_list(arr) && len(arr)>0, "The array is not a vector neither a matrix." )
+    is_list(arr[0])
+    ?   let( l0 = len(arr[0]) )
+        assert([for(a=arr) if(!is_list(a) || len(a)!=l0) 1 ]==[], "The array is not a vector neither a matrix." )
+        reverse
+        ? [for (i=[0:1:l0-1]) 
+              [ for (j=[0:1:len(arr)-1]) arr[len(arr)-1-j][l0-1-i] ] ] 
+        : [for (i=[0:1:l0-1]) 
+              [ for (j=[0:1:len(arr)-1]) arr[j][i] ] ] 
+    :  assert( is_vector(arr), "The array is not a vector neither a matrix." )
+           arr;
 
 
 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
diff --git a/common.scad b/common.scad
index d8941ee..eebd141 100644
--- a/common.scad
+++ b/common.scad
@@ -140,7 +140,7 @@ function _list_pattern(list) =
 //   is_consistent(list)
 // Description:
 //   Tests whether input is a list of entries which all have the same list structure
-//   and are filled with finite numerical data.
+//   and are filled with finite numerical data. It returns `true`for the empty list. 
 // Example:
 //   is_consistent([3,4,5]);              // Returns true
 //   is_consistent([[3,4],[4,5],[6,7]]);  // Returns true
diff --git a/math.scad b/math.scad
index e9b8579..158b1f2 100644
--- a/math.scad
+++ b/math.scad
@@ -84,7 +84,7 @@ function hypot(x,y,z=0) =
 //   y = factorial(6);  // Returns: 720
 //   z = factorial(9);  // Returns: 362880
 function factorial(n,d=0) =
-    assert(is_int(n) && is_int(d) && n>=0 && d>=0, "Factorial is not defined for negative numbers")
+    assert(is_int(n) && is_int(d) && n>=0 && d>=0, "Factorial is defined only for non negative integers")
     assert(d<=n, "d cannot be larger than n")
     product([1,for (i=[n:-1:d+1]) i]);
 
@@ -164,7 +164,7 @@ function binomial_coefficient(n,k) =
 function lerp(a,b,u) =
     assert(same_shape(a,b), "Bad or inconsistent inputs to lerp")
     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.")
+    assert(is_finite(u) || is_vector(u) || valid_range(u), "Input u to lerp must be a number, vector, or valid range.")
     [for (v = u) (1-v)*a + v*b ];
 
 
@@ -387,12 +387,13 @@ 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.")
+    assert( is_finite(x+y+step+m) && !approx(m,0), "Input must be finite numbers and the module value cannot be zero." )
     let(
         a = posmod(x, m),
         b = posmod(y, m),
-        c = step>0? (a>b? b+m : b) : (a<b? b-m : b)
-    ) [for (i=[a:step:c]) (i%m+m)%m];
+        c = step>0? (a>b? b+m : b) 
+            : (a<b? b-m : b)
+    ) [for (i=[a:step:c]) (i%m+m)%m ];
 
 
 
@@ -536,9 +537,13 @@ function _sum(v,_total,_i=0) = _i>=len(v) ? _total : _sum(v,_total+v[_i], _i+1);
 //   cumsum([2,2,2]);  // returns [2,4,6]
 //   cumsum([1,2,3]);  // returns [1,3,6]
 //   cumsum([[1,2,3], [3,4,5], [5,6,7]]);  // returns [[1,2,3], [4,6,8], [9,12,15]]
-function cumsum(v,_i=0,_acc=[]) =
+function cumsum(v) =
+    assert(is_consistent(v), "The input is not consistent." )
+    _cumsum(v,_i=0,_acc=[]);
+
+function _cumsum(v,_i=0,_acc=[]) =
     _i==len(v) ? _acc :
-    cumsum(
+    _cumsum(
         v, _i+1,
         concat(
             _acc,
@@ -598,7 +603,7 @@ function deltas(v) =
 // Description:
 //   Returns the product of all entries in the given list.
 //   If passed a list of vectors of same dimension, returns a vector of products of each part.
-//   If passed a list of square matrices, returns a the resulting product matrix.
+//   If passed a list of square matrices, returns the resulting product matrix.
 // Arguments:
 //   v = The list to get the product of.
 // Example:
@@ -606,7 +611,7 @@ function deltas(v) =
 //   product([[1,2,3], [3,4,5], [5,6,7]]);  // returns [15, 48, 105]
 function product(v) = 
     assert( is_vector(v) || is_matrix(v) || ( is_matrix(v[0],square=true) && is_consistent(v)), 
-            "Invalid input.")
+    "Invalid input.")
     _product(v, 1, v[0]);
 
 function _product(v, i=0, _tot) = 
@@ -691,9 +696,11 @@ function linear_solve(A,b) =
         zeros = [for(i=[0:mindim-1]) if (approx(R[i][i],0)) i]
     )
     zeros != [] ? [] :
-    m<n ? Q*back_substitute(R,b,transpose=true) :
-    back_substitute(R, transpose(Q)*b);
-
+    m<n 
+    // avoiding input validation in back_substitute
+    ? let( n  = len(R) )
+      Q*reverse(_back_substitute(transpose(R, reverse=true), reverse(b))) 
+    : _back_substitute(R, transpose(Q)*b);
 
 // Function: matrix_inverse()
 // Usage:
@@ -708,18 +715,6 @@ function matrix_inverse(A) =
     linear_solve(A,ident(len(A)));
 
 
-// Function: submatrix()
-// Usage: submatrix(M, ind1, ind2)
-// Description:
-//   Returns a submatrix with the specified index ranges or index sets.  
-function submatrix(M,ind1,ind2) =
-    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)
 // Description:
@@ -732,13 +727,14 @@ function qr_factor(A) =
       n = len(A[0])
     )
     let(
-        qr =_qr_factor(A, column=0, m = m, n=n, Q=ident(m)),
-        Rzero = [
-            for(i=[0:m-1]) [
-                for(j=[0:n-1])
-                i>j ? 0 : qr[1][i][j]
+        qr =_qr_factor(A, Q=ident(m), column=0, m = m, n=n),
+        Rzero = 
+          let( R = qr[1] )
+          [ for(i=[0:m-1]) [
+              let( ri = R[i] )
+              for(j=[0:n-1]) i>j ? 0 : ri[j]
             ]
-        ]
+          ]
     ) [qr[0],Rzero];
 
 function _qr_factor(A,Q, column, m, n) =
@@ -767,13 +763,12 @@ function back_substitute(R, b, transpose = false) =
     let(n=len(R))
     assert(is_vector(b,n) || is_matrix(b,n),str("R and b are not compatible in back_substitute ",n, len(b)))
     transpose
-      ? reverse(_back_substitute([for(i=[0:n-1]) [for(j=[0:n-1]) R[n-1-j][n-1-i]]],
-                                 reverse(b)))  
+      ? reverse(_back_substitute(transpose(R, reverse=true), reverse(b)))  
       : _back_substitute(R,b);
 
 function _back_substitute(R, b, x=[]) =
-      let(n=len(R))
-      len(x) == n ? x
+    let(n=len(R))
+    len(x) == n ? x
     : let(ind = n - len(x) - 1)
       R[ind][ind] == 0 ? []
     : let(
@@ -793,7 +788,7 @@ function _back_substitute(R, b, x=[]) =
 //   M = [ [6,-2], [1,8] ];
 //   det = det2(M);  // Returns: 50
 function det2(M) = 
-    assert( is_matrix(M,2,2), "Matrix should be 2x2." )
+    assert( 0*M==[[0,0],[0,0]], "Matrix should be 2x2." )
     M[0][0] * M[1][1] - M[0][1]*M[1][0];
 
 
@@ -806,7 +801,7 @@ function det2(M) =
 //   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." )
+    assert( 0*M==[[0,0,0],[0,0,0],[0,0,0]], "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]);
@@ -846,7 +841,7 @@ function determinant(M) =
 // Description:
 //   Returns true if A is a numeric matrix of height m and width n.  If m or n
 //   are omitted or set to undef then true is returned for any positive dimension.
-//   If `square` is true then the matrix is required to be square.  Note if you
+//   If `square` is true then the matrix is required to be square. 
 //   specify m != n and require a square matrix then the result will always be false.
 // Arguments:
 //   A = matrix to test
@@ -855,10 +850,10 @@ function determinant(M) =
 //   square = set to true to require a square matrix.  Default: false        
 function is_matrix(A,m,n,square=false) =
     is_list(A[0]) 
-    && ( let(v = A*A[0]) is_num(0*(v*v)) ) // a matrix of finite numbers 
-    && (is_undef(n) || len(A[0])==n )
-    && (is_undef(m) || len(A)==m )
-    && ( !square || len(A)==len(A[0]));
+    && ( 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
@@ -948,13 +943,16 @@ function compare_lists(a, b) =
 //   any([1,5,true]);       // Returns true.
 //   any([[0,0], [0,0]]);   // Returns false.
 //   any([[0,0], [1,0]]);   // Returns true.
-function any(l, i=0, succ=false) =
-    (i>=len(l) || succ)? succ :
-    any( l, 
-         i+1, 
-         succ = is_list(l[i]) ? any(l[i]) : !(!l[i])
-        );
+function any(l) =
+    assert(is_list(l), "The input is not a list." )
+    _any(l, i=0, succ=false);
 
+function _any(l, i=0, succ=false) =
+    (i>=len(l) || succ)? succ :
+    _any( l, 
+         i+1, 
+         succ = is_list(l[i]) ? _any(l[i]) : !(!l[i])
+        );
 
 
 // Function: all()
@@ -971,12 +969,15 @@ 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+1,
-         fail = is_list(l[i]) ? !all(l[i]) : !l[i]
-        ) ;
+    assert( is_list(l), "The input is not a list." )
+    _all(l, i=0, fail=false);
 
+function _all(l, i=0, fail=false) =
+    (i>=len(l) || fail)? !fail :
+    _all( l, 
+         i+1,
+         fail = is_list(l[i]) ? !_all(l[i]) : !l[i]
+        ) ;
 
 
 // Function: count_true()
@@ -999,16 +1000,6 @@ function all(l, i=0, fail=false) =
 //   count_true([[0,0], [1,0]]);   // Returns 1.
 //   count_true([[1,1], [1,1]]);   // Returns 4.
 //   count_true([[1,1], [1,1]], nmax=3);  // Returns 3.
-function count_true(l, nmax=undef, i=0, cnt=0) =
-    (i>=len(l) || (nmax!=undef && cnt>=nmax))? cnt :
-    count_true(
-        l=l, nmax=nmax, i=i+1, cnt=cnt+(
-            is_list(l[i])? count_true(l[i], nmax=nmax-cnt) :
-            (l[i]? 1 : 0)
-        )
-    );
-
-
 function count_true(l, nmax) = 
     !is_list(l) ? !(!l) ? 1: 0 :
     let( c = [for( i = 0,
@@ -1204,12 +1195,12 @@ function C_div(z1,z2) =
 //   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.
 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]);
+    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:
@@ -1219,35 +1210,15 @@ function polynomial(p,z,k,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) 
-               && []==[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]
-               ])
-         ]);
+    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")
+    p*p==0 || q*q==0
+    ? [0]
+    : _poly_trim(convolve(p,q));
 
     
 // Function: poly_div()
@@ -1258,19 +1229,23 @@ 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) =
-    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);
+function poly_div(n,d) =
+    assert( is_vector(n) && is_vector(d) , "Invalid polynomials." )
+    let( d = _poly_trim(d), 
+         n = _poly_trim(n) )
+    assert( d!=[0] , "Denominator cannot be a zero polynomial." )
+    n==[0]
+    ? [[0],[0]]
+    : _poly_div(n,d,q=[]);
+
+function _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()
@@ -1401,4 +1376,4 @@ function real_roots(p,eps=undef,tol=1e-14) =
     ? [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
+// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
\ No newline at end of file
diff --git a/paths.scad b/paths.scad
index 7d5f978..93d82fa 100644
--- a/paths.scad
+++ b/paths.scad
@@ -43,10 +43,12 @@ include <triangulation.scad>
 //   dim = list of allowed dimensions of the vectors in the path.  Default: [2,3]
 //   fast = set to true for fast check that only looks at first entry.  Default: false
 function is_path(list, dim=[2,3], fast=false) =
-    fast? is_list(list) && is_vector(list[0]) :
-    is_list(list) && is_list(list[0]) && len(list)>1 &&
-    (is_undef(dim) || in_list(len(list[0]), force_list(dim))) &&
-    is_list_of(list, repeat(0,len(list[0])));
+    fast
+    ?   is_list(list) && is_vector(list[0]) 
+    :   is_matrix(list) 
+        && len(list)>1 
+        && len(list[0])>0
+        && (is_undef(dim) || in_list(len(list[0]), force_list(dim)));
 
 
 // Function: is_closed_path()
@@ -105,32 +107,51 @@ function path_subselect(path, s1, u1, s2, u2, closed=false) =
 
 // Function: simplify_path()
 // Description:
-//   Takes a path and removes unnecessary collinear points.
+//   Takes a path and removes unnecessary subsequent collinear points.
 // Usage:
 //   simplify_path(path, [eps])
 // Arguments:
-//   path = A list of 2D path points.
+//   path = A list of path points of any dimension.
 //   eps = Largest positional variance allowed.  Default: `EPSILON` (1-e9)
 function simplify_path(path, eps=EPSILON) =
-    len(path)<=2? path : let(
-        indices = concat([0], [for (i=[1:1:len(path)-2]) if (!collinear_indexed(path, i-1, i, i+1, eps=eps)) i], [len(path)-1])
-    ) [for (i = indices) path[i]];
+    assert( is_path(path), "Invalid path." )
+    assert( is_undef(eps) || (is_finite(eps) && (eps>=0) ), "Invalid tolerance." )    
+    len(path)<=2 ? path 
+    :   let(
+            indices = [ 0,
+                        for (i=[1:1:len(path)-2]) 
+                            if (!collinear(path[i-1],path[i],path[i+1], eps=eps)) i, 
+                        len(path)-1 
+                      ]
+        ) 
+        [for (i = indices) path[i] ];
 
 
 // Function: simplify_path_indexed()
 // Description:
-//   Takes a list of points, and a path as a list of indices into `points`,
-//   and removes all path points that are unecessarily collinear.
+//   Takes a list of points, and a list of indices into `points`,
+//   and removes from the list all indices of subsequent indexed points that are unecessarily collinear.
+//   Returns the list of the remained indices.
 // Usage:
-//   simplify_path_indexed(path, eps)
+//   simplify_path_indexed(points,indices, eps)
 // Arguments:
 //   points = A list of points.
-//   path = A list of indices into `points` that forms a path.
+//   indices = A list of indices into `points` that forms a path.
 //   eps = Largest angle variance allowed.  Default: EPSILON (1-e9) degrees.
-function simplify_path_indexed(points, path, eps=EPSILON) =
-    len(path)<=2? path : let(
-        indices = concat([0], [for (i=[1:1:len(path)-2]) if (!collinear_indexed(points, path[i-1], path[i], path[i+1], eps=eps)) i], [len(path)-1])
-    ) [for (i = indices) path[i]];
+function simplify_path_indexed(points, indices, eps=EPSILON) =
+    len(indices)<=2? indices 
+    : let(
+        indices = concat( indices[0], 
+                          [for (i=[1:1:len(indices)-2]) 
+                              let( 
+                                  i1 = indices[i-1],
+                                  i2 = indices[i],
+                                  i3 = indices[i+1]
+                              )
+                              if (!collinear(points[i1],points[i2],points[i3], eps=eps)) indices[i]], 
+                          indices[len(indices)-1] )
+        ) 
+        indices;
 
 
 // Function: path_length()
diff --git a/tests/test_arrays.scad b/tests/test_arrays.scad
index 4fffb7d..8df23a8 100644
--- a/tests/test_arrays.scad
+++ b/tests/test_arrays.scad
@@ -3,14 +3,6 @@ include <../std.scad>
 
 // Section: List Query Operations
 
-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() {
     l = [3,4,5,6,7,8,9];
     assert(select(l, 5, 6) == [8,9]);
@@ -365,6 +357,27 @@ module test_subindex() {
 test_subindex();
 
 
+// Need decision about behavior for out of bounds ranges, empty ranges
+module test_submatrix(){
+  M = [[1,2,3,4,5],
+       [6,7,8,9,10],
+       [11,12,13,14,15],
+       [16,17,18,19,20],
+       [21,22,23,24,25]];
+  assert_equal(submatrix(M,[1:2], [3:4]), [[9,10],[14,15]]);
+  assert_equal(submatrix(M,[1], [3,4]), [[9,10]]);
+  assert_equal(submatrix(M,1, [3,4]), [[9,10]]);
+  assert_equal(submatrix(M, [3,4],1), [[17],[22]]);
+  assert_equal(submatrix(M, [1,3],[2,4]), [[8,10],[18,20]]);
+  assert_equal(submatrix(M, 1,3), [[9]]);
+  A = [[true,    17, "test"],
+     [[4,2],   91, false],
+     [6,    [3,4], undef]];
+  assert_equal(submatrix(A,[0,2],[1,2]),[[17, "test"], [[3, 4], undef]]);
+}
+test_submatrix();
+
+
 module test_force_list() {
     assert_equal(force_list([3,4,5]), [3,4,5]);
     assert_equal(force_list(5), [5]);
@@ -467,6 +480,7 @@ test_array_dim();
 module test_transpose() {
     assert(transpose([[1,2,3],[4,5,6],[7,8,9]]) == [[1,4,7],[2,5,8],[3,6,9]]);
     assert(transpose([[1,2,3],[4,5,6]]) == [[1,4],[2,5],[3,6]]);
+    assert(transpose([[1,2,3],[4,5,6]],reverse=true) == [[6,3], [5,2], [4,1]]);
     assert(transpose([3,4,5]) == [3,4,5]);
 }
 test_transpose();
diff --git a/tests/test_math.scad b/tests/test_math.scad
index e5b1a41..12f7d2e 100644
--- a/tests/test_math.scad
+++ b/tests/test_math.scad
@@ -853,22 +853,6 @@ module test_real_roots(){
 }
 test_real_roots();
 
-// Need decision about behavior for out of bounds ranges, empty ranges
-module test_submatrix(){
-  M = [[1,2,3,4,5],
-       [6,7,8,9,10],
-       [11,12,13,14,15],
-       [16,17,18,19,20],
-       [21,22,23,24,25]];
-  assert_equal(submatrix(M,[1:2], [3:4]), [[9,10],[14,15]]);
-  assert_equal(submatrix(M,[1], [3,4]), [[9,10]]);
-  assert_equal(submatrix(M,1, [3,4]), [[9,10]]);
-  assert_equal(submatrix(M, [3,4],1), [[17],[22]]);
-  assert_equal(submatrix(M, [1,3],[2,4]), [[8,10],[18,20]]);
-}
-test_submatrix();
-
-
 
 module test_qr_factor() {
   // Check that R is upper triangular
@@ -913,23 +897,21 @@ 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],[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([3,2,1],[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]]),[]);
+  assert_equal(poly_mult([[3,4,5],[0,0,0]]), [0]);
+  assert_equal(poly_mult([[0],[0,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],[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]]);
+  assert_equal(poly_div([0], [1,-3]), [[0],[0]]);
 }
 test_poly_div();
 
@@ -942,4 +924,4 @@ module test_poly_add(){
 }
 test_poly_add();
 
-// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
+// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
\ No newline at end of file