From a75a368dd510598be0dd261d328074812a4b1740 Mon Sep 17 00:00:00 2001
From: Adrian Mariano <avm4@cornell.edu>
Date: Tue, 26 Oct 2021 16:45:14 -0400
Subject: [PATCH] Move linear algebra to linalg.scad columns->column because
 the multiindex case is handled by submatrix and also it never occurs in the
 code.

---
 arrays.scad            | 417 ++++--------------------
 attachments.scad       |  16 +-
 beziers.scad           |  20 +-
 linalg.scad            | 713 +++++++++++++++++++++++++++++++++++++++++
 math.scad              | 344 +-------------------
 paths.scad             |   4 +-
 regions.scad           |   2 +-
 rounding.scad          |  18 +-
 shapes2d.scad          |   6 +-
 shapes3d.scad          |   2 +-
 skin.scad              |   2 +-
 std.scad               |   1 +
 tests/test_arrays.scad | 100 +-----
 tests/test_math.scad   | 155 +--------
 vectors.scad           |   4 +-
 15 files changed, 812 insertions(+), 992 deletions(-)
 create mode 100644 linalg.scad

diff --git a/arrays.scad b/arrays.scad
index 5c0e7e7..47e3dde 100644
--- a/arrays.scad
+++ b/arrays.scad
@@ -5,7 +5,6 @@
 //   include <BOSL2/std.scad>
 //////////////////////////////////////////////////////////////////////
 
-
 // Terminology:
 //   **List** = An ordered collection of zero or more items.  ie: `["a", "b", "c"]`
 //   **Vector** = A list of numbers. ie: `[4, 5, 6]`
@@ -251,7 +250,7 @@ function __find_approx(val, list, eps, i=0) =
 //   list = The list to get the portion of.
 //   start = Either the index of the first item or an index range or a list of indices.
 //   end = The index of the last item when `start` is a number. When `start` is a list or a range, `end` should not be given.
-// See Also: slice(), columns(), last()
+// See Also: slice(), column(), last()
 // Example:
 //   l = [3,4,5,6,7,8,9];
 //   a = select(l, 5, 6);   // Returns [8,9]
@@ -292,7 +291,7 @@ function select(list, start, end) =
 //   list = The list to get the slice of.
 //   s = The index of the first item to return.
 //   e = The index of the last item to return.
-// See Also: select(), columns(), last()
+// See Also: select(), column(), last()
 // Example:
 //   a = slice([3,4,5,6,7,8,9], 3, 5);   // Returns [6,7,8]
 //   b = slice([3,4,5,6,7,8,9], 2, -1);  // Returns [5,6,7,8,9]
@@ -317,7 +316,7 @@ function slice(list,s=0,e=-1) =
 // Usage:
 //   item = last(list);
 // Topics: List Handling
-// See Also: select(), slice(), columns()
+// See Also: select(), slice(), column()
 // Description:
 //   Returns the last element of a list, or undef if empty.
 // Arguments:
@@ -1471,63 +1470,6 @@ function permutations(l,n=2) =
       : [for (i=idx(l), p=permutations([for (j=idx(l)) if (i!=j) l[j]], n=n-1)) concat([l[i]], p)];
 
 
-// Function: zip()
-// Usage:
-//   pairs = zip(a,b);
-//   triples = zip(a,b,c);
-//   quads = zip([LIST1,LIST2,LIST3,LIST4]);
-// Topics: List Handling, Iteration
-// See Also: zip_long()
-// Description:
-//   Zips together two or more lists into a single list.  For example, if you have two
-//   lists [3,4,5], and [8,7,6], and zip them together, you get [ [3,8],[4,7],[5,6] ].
-//   The list returned will be as long as the shortest list passed to zip().
-// Arguments:
-//   a = The first list, or a list of lists if b and c are not given.
-//   b = The second list, if given.
-//   c = The third list, if given.
-// Example:
-//   a = [9,8,7,6]; b = [1,2,3];
-//   for (p=zip(a,b)) echo(p);
-//   // ECHO: [9,1]
-//   // ECHO: [8,2]
-//   // ECHO: [7,3]
-function zip(a,b,c) =
-    b!=undef? zip([a,b,if (c!=undef) c]) :
-    let(n = min_length(a))
-    [for (i=[0:1:n-1]) [for (x=a) x[i]]];
-
-
-// Function: zip_long()
-// Usage:
-//   pairs = zip_long(a,b);
-//   triples = zip_long(a,b,c);
-//   quads = zip_long([LIST1,LIST2,LIST3,LIST4]);
-// Topics: List Handling, Iteration
-// See Also: zip()
-// Description:
-//   Zips together two or more lists into a single list.  For example, if you have two
-//   lists [3,4,5], and [8,7,6], and zip them together, you get [ [3,8],[4,7],[5,6] ].
-//   The list returned will be as long as the longest list passed to zip_long(), with
-//   shorter lists padded by the value in `fill`.
-// Arguments:
-//   a = The first list, or a list of lists if b and c are not given.
-//   b = The second list, if given.
-//   c = The third list, if given.
-//   fill = The value to pad shorter lists with.  Default: undef
-// Example:
-//   a = [9,8,7,6]; b = [1,2,3];
-//   for (p=zip_long(a,b,fill=88)) echo(p);
-//   // ECHO: [9,1]
-//   // ECHO: [8,2]
-//   // ECHO: [7,3]
-//   // ECHO: [6,88]]
-function zip_long(a,b,c,fill) =
-    b!=undef? zip_long([a,b,if (c!=undef) c],fill=fill) :
-    let(n = max_length(a))
-    [for (i=[0:1:n-1]) [for (x=a) i<len(x)? x[i] : fill]];
-
-
 
 // Section: Set Manipulation
 
@@ -1616,209 +1558,7 @@ function set_intersection(a, b) =
 
 
 
-// Section: Array Manipulation
-
-// Function: columns()
-// Usage:
-//   list = columns(M, idx);
-// Topics: Array Handling, List Handling
-// See Also: select(), slice()
-// Description:
-//   Extracts the entries listed in idx from each entry in M.  For a matrix this means
-//   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.
-// Example:
-//   M = [[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]];
-//   a = columns(M,2);      // Returns [3, 7, 11, 15]
-//   b = columns(M,[2]);    // Returns [[3], [7], [11], [15]]
-//   c = columns(M,[2,1]);  // Returns [[3, 2], [7, 6], [11, 10], [15, 14]]
-//   d = columns(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] ];
-//   e = columns(N,[0,1]);  // Returns [ [1,2], [3,undef], [4,5], [6,7] ]
-function columns(M, 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:
-//   mat = submatrix(M, idx1, idx2);
-// Topics: Matrices, Array Handling
-// See Also: columns(), block_matrix(), submatrix_set()
-// 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 columns 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: hstack()
-// Usage: 
-//   A = hstack(M1, M2)
-//   A = hstack(M1, M2, M3)
-//   A = hstack([M1, M2, M3, ...])
-// Topics: Matrices, Array Handling
-// See Also: columns(), submatrix(), block_matrix()
-// Description:
-//   Constructs a matrix by horizontally "stacking" together compatible matrices or vectors.  Vectors are treated as columsn in the stack.
-//   This command is the inverse of `columns`.  Note: strings given in vectors are broken apart into lists of characters.  Strings given
-//   in matrices are preserved as strings.  If you need to combine vectors of strings use array_group as shown below to convert the
-//   vector into a column matrix.  Also note that vertical stacking can be done directly with concat.  
-// Arguments:
-//   M1 = If given with other arguments, the first matrix (or vector) to stack.  If given alone, a list of matrices/vectors to stack. 
-//   M2 = Second matrix/vector to stack
-//   M3 = Third matrix/vector to stack.
-// Example:
-//   M = ident(3);
-//   v1 = [2,3,4];
-//   v2 = [5,6,7];
-//   v3 = [8,9,10];
-//   a = hstack(v1,v2);     // Returns [[2, 5], [3, 6], [4, 7]]
-//   b = hstack(v1,v2,v3);  // Returns [[2, 5,  8],
-//                          //          [3, 6,  9],
-//                          //          [4, 7, 10]]
-//   c = hstack([M,v1,M]);  // Returns [[1, 0, 0, 2, 1, 0, 0],
-//                          //          [0, 1, 0, 3, 0, 1, 0],
-//                          //          [0, 0, 1, 4, 0, 0, 1]]
-//   d = hstack(columns(M,0), columns(M,[1 2]));  // Returns M
-//   strvec = ["one","two"];
-//   strmat = [["three","four"], ["five","six"]];
-//   e = hstack(strvec,strvec); // Returns [["o", "n", "e", "o", "n", "e"],
-//                              //          ["t", "w", "o", "t", "w", "o"]]
-//   f = hstack(array_group(strvec,1), array_group(strvec,1));
-//                              // Returns [["one", "one"],
-//                              //          ["two", "two"]]
-//   g = hstack(strmat,strmat); //  Returns: [["three", "four", "three", "four"],
-//                              //            [ "five",  "six",  "five",  "six"]]
-function hstack(M1, M2, M3) =
-    (M3!=undef)? hstack([M1,M2,M3]) : 
-    (M2!=undef)? hstack([M1,M2]) :
-    assert(all([for(v=M1) is_list(v)]), "One of the inputs to hstack is not a list")
-    let(
-        minlen = min_length(M1),
-        maxlen = max_length(M1)
-    )
-    assert(minlen==maxlen, "Input vectors to hstack must have the same length")
-    [for(row=[0:1:minlen-1])
-        [for(matrix=M1)
-           each matrix[row]
-        ]
-    ];
-
-
-// Function: block_matrix()
-// Usage:
-//    bmat = block_matrix([[M11, M12,...],[M21, M22,...], ... ]);
-// Topics: Matrices, Array Handling
-// See Also: columns(), submatrix()
-// Description:
-//    Create a block matrix by supplying a matrix of matrices, which will
-//    be combined into one unified matrix.  Every matrix in one row
-//    must have the same height, and the combined width of the matrices
-//    in each row must be equal. Strings will stay strings. 
-// Example:
-//  A = [[1,2],
-//       [3,4]];
-//  B = ident(2);
-//  C = block_matrix([[A,B],[B,A],[A,B]]);
-//      // Returns:
-//      //        [[1, 2, 1, 0],
-//      //         [3, 4, 0, 1],
-//      //         [1, 0, 1, 2],
-//      //         [0, 1, 3, 4],
-//      //         [1, 2, 1, 0],
-//      //         [3, 4, 0, 1]]);
-//  D = block_matrix([[A,B], ident(4)]);
-//      // Returns:
-//      //        [[1, 2, 1, 0],
-//      //         [3, 4, 0, 1],
-//      //         [1, 0, 0, 0],
-//      //         [0, 1, 0, 0],
-//      //         [0, 0, 1, 0],
-//      //         [0, 0, 0, 1]]);
-//  E = [["one", "two"], [3,4]];
-//  F = block_matrix([[E,E]]);
-//      // Returns:
-//      //        [["one", "two", "one", "two"],
-//      //         [    3,     4,     3,     4]]
-function block_matrix(M) =
-    let(
-        bigM = [for(bigrow = M) each hstack(bigrow)],
-        len0 = len(bigM[0]),
-        badrows = [for(row=bigM) if (len(row)!=len0) 1]
-    )
-    assert(badrows==[], "Inconsistent or invalid input")
-    bigM;
-
-// Function: diagonal_matrix()
-// Usage:
-//   mat = diagonal_matrix(diag, [offdiag]);
-// Topics: Matrices, Array Handling
-// See Also: columns(), submatrix()
-// Description:
-//   Creates a square matrix with the items in the list `diag` on
-//   its diagonal.  The off diagonal entries are set to offdiag,
-//   which is zero by default. 
-// Arguments:
-//   diag = A list of items to put in the diagnal cells of the matrix.
-//   offdiag = Value to put in non-diagonal matrix cells.
-function diagonal_matrix(diag, offdiag=0) =
-  assert(is_list(diag) && len(diag)>0)
-  [for(i=[0:1:len(diag)-1]) [for(j=[0:len(diag)-1]) i==j?diag[i] : offdiag]];
-
-
-// Function: submatrix_set()
-// Usage:
-//   mat = submatrix_set(M, A, [m], [n]);
-// Topics: Matrices, Array Handling
-// See Also: columns(), submatrix()
-// Description:
-//   Sets a submatrix of M equal to the matrix A.  By default the top left corner of M is set to A, but
-//   you can specify offset coordinates m and n.  If A (as adjusted by m and n) extends beyond the bounds
-//   of M then the extra entries are ignored.  You can pass in A=[[]], a null matrix, and M will be
-//   returned unchanged.  Note that the input M need not be rectangular in shape.  
-// Arguments:
-//   M = Original matrix.
-//   A = Sub-matrix of parts to set.
-//   m = Row number of upper-left corner to place A at.
-//   n = Column number of upper-left corner to place A at.
-function submatrix_set(M,A,m=0,n=0) =
-    assert(is_list(M))
-    assert(is_list(A))
-    assert(is_int(m))
-    assert(is_int(n))
-    let( badrows = [for(i=idx(A)) if (!is_list(A[i])) i])
-    assert(badrows==[], str("Input submatrix malformed rows: ",badrows))
-    [for(i=[0:1:len(M)-1])
-        assert(is_list(M[i]), str("Row ",i," of input matrix is not a list"))
-        [for(j=[0:1:len(M[i])-1]) 
-            i>=m && i <len(A)+m && j>=n && j<len(A[0])+n ? A[i-m][j-n] : M[i][j]]];
+// Section: Changing list structure
 
 
 // Function: array_group()
@@ -1828,7 +1568,7 @@ function submatrix_set(M,A,m=0,n=0) =
 //   Takes a flat array of values, and groups items in sets of `cnt` length.
 //   The opposite of this is `flatten()`.
 // Topics: Matrices, Array Handling
-// See Also: columns(), submatrix(), hstack(), flatten(), full_flatten()
+// See Also: column(), submatrix(), hstack(), flatten(), full_flatten()
 // Arguments:
 //   v = The list of items to group.
 //   cnt = The number of items to put in each grouping.  Default:2
@@ -1847,7 +1587,7 @@ function array_group(v, cnt=2, dflt=0) =
 // Usage:
 //   list = flatten(l);
 // Topics: Matrices, Array Handling
-// See Also: columns(), submatrix(), hstack(), full_flatten()
+// See Also: column(), submatrix(), hstack(), full_flatten()
 // Description:
 //   Takes a list of lists and flattens it by one level.
 // Arguments:
@@ -1863,7 +1603,7 @@ function flatten(l) =
 // Usage:
 //   list = full_flatten(l);
 // Topics: Matrices, Array Handling
-// See Also: columns(), submatrix(), hstack(), flatten()
+// See Also: column(), submatrix(), hstack(), 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.
@@ -1925,110 +1665,63 @@ function array_dim(v, depth=undef) =
         :  let( dimlist = _array_dim_recurse(v))
            (depth > len(dimlist))? 0 : dimlist[depth-1] ;
            
-           
 
-
-// Function: transpose()
+// Function: zip()
 // Usage:
-//    arr = transpose(arr, [reverse]);
-// Topics: Matrices, Array Handling
-// See Also: submatrix(), block_matrix(), hstack(), flatten()
+//   pairs = zip(a,b);
+//   triples = zip(a,b,c);
+//   quads = zip([LIST1,LIST2,LIST3,LIST4]);
+// Topics: List Handling, Iteration
+// See Also: zip_long()
 // Description:
-//    Returns the transpose of the given input array.  The input should be a list of lists that are
-//    all the same length.  If you give a vector then transpose returns it unchanged.  
-//    When reverse=true, the transpose is done across to the secondary diagonal.  (See example below.)
-//    By default, reverse=false.
-// Example:
-//   arr = [
-//       ["a", "b", "c"],
-//       ["d", "e", "f"],
-//       ["g", "h", "i"]
-//   ];
-//   t = transpose(arr);
-//   // Returns:
-//   // [
-//   //     ["a", "d", "g"],
-//   //     ["b", "e", "h"],
-//   //     ["c", "f", "i"],
-//   // ]
-// Example:
-//   arr = [
-//       ["a", "b", "c"],
-//       ["d", "e", "f"]
-//   ];
-//   t = transpose(arr);
-//   // Returns:
-//   // [
-//   //     ["a", "d"],
-//   //     ["b", "e"],
-//   //     ["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 on a list of numbers returns the list unchanged
-//   transpose([3,4,5]);  // Returns: [3,4,5]
-function transpose(arr, reverse=false) =
-    assert( is_list(arr) && len(arr)>0, "Input to transpose must be a nonempty list.")
-    is_list(arr[0])
-    ?   let( len0 = len(arr[0]) )
-        assert([for(a=arr) if(!is_list(a) || len(a)!=len0) 1 ]==[], "Input to transpose has inconsistent row lengths." )
-        reverse
-        ? [for (i=[0:1:len0-1]) 
-              [ for (j=[0:1:len(arr)-1]) arr[len(arr)-1-j][len0-1-i] ] ] 
-        : [for (i=[0:1:len0-1]) 
-              [ for (j=[0:1:len(arr)-1]) arr[j][i] ] ] 
-    :  assert( is_vector(arr), "Input to transpose must be a vector or list of lists.")
-           arr;
-
-
-// Section: Matrices
-
-// Function: is_matrix_symmetric()
-// Usage:
-//   b = is_matrix_symmetric(A, [eps])
-// Description:
-//   Returns true if the input matrix is symmetric, meaning it equals its transpose.
-//   Matrix should have numerical entries.
+//   Zips together two or more lists into a single list.  For example, if you have two
+//   lists [3,4,5], and [8,7,6], and zip them together, you get [ [3,8],[4,7],[5,6] ].
+//   The list returned will be as long as the shortest list passed to zip().
 // Arguments:
-//   A = matrix to test
-//   eps = epsilon for comparing equality.  Default: 1e-12
-function is_matrix_symmetric(A,eps=1e-12) =
-    approx(A,transpose(A), eps);
+//   a = The first list, or a list of lists if b and c are not given.
+//   b = The second list, if given.
+//   c = The third list, if given.
+// Example:
+//   a = [9,8,7,6]; b = [1,2,3];
+//   for (p=zip(a,b)) echo(p);
+//   // ECHO: [9,1]
+//   // ECHO: [8,2]
+//   // ECHO: [7,3]
+function zip(a,b,c) =
+    b!=undef? zip([a,b,if (c!=undef) c]) :
+    let(n = min_length(a))
+    [for (i=[0:1:n-1]) [for (x=a) x[i]]];
 
 
-// Function&Module: echo_matrix()
+// Function: zip_long()
 // Usage:
-//    echo_matrix(M, [description=], [sig=], [eps=]);
-//    dummy = echo_matrix(M, [description=], [sig=], [eps=]),
+//   pairs = zip_long(a,b);
+//   triples = zip_long(a,b,c);
+//   quads = zip_long([LIST1,LIST2,LIST3,LIST4]);
+// Topics: List Handling, Iteration
+// See Also: zip()
 // Description:
-//    Display a numerical matrix in a readable columnar format with `sig` significant
-//    digits.  Values smaller than eps display as zero.  If you give a description
-//    it is displayed at the top.  
-function echo_matrix(M,description,sig=4,eps=1e-9) =
-  let(
-      horiz_line = chr(8213),
-      matstr = matrix_strings(M,sig=sig,eps=eps),
-      separator = str_join(repeat(horiz_line,10)),
-      dummy=echo(str(separator,"  ",is_def(description) ? description : ""))
-            [for(row=matstr) echo(row)]
-  )
-  echo(separator);
+//   Zips together two or more lists into a single list.  For example, if you have two
+//   lists [3,4,5], and [8,7,6], and zip them together, you get [ [3,8],[4,7],[5,6] ].
+//   The list returned will be as long as the longest list passed to zip_long(), with
+//   shorter lists padded by the value in `fill`.
+// Arguments:
+//   a = The first list, or a list of lists if b and c are not given.
+//   b = The second list, if given.
+//   c = The third list, if given.
+//   fill = The value to pad shorter lists with.  Default: undef
+// Example:
+//   a = [9,8,7,6]; b = [1,2,3];
+//   for (p=zip_long(a,b,fill=88)) echo(p);
+//   // ECHO: [9,1]
+//   // ECHO: [8,2]
+//   // ECHO: [7,3]
+//   // ECHO: [6,88]]
+function zip_long(a,b,c,fill) =
+    b!=undef? zip_long([a,b,if (c!=undef) c],fill=fill) :
+    let(n = max_length(a))
+    [for (i=[0:1:n-1]) [for (x=a) i<len(x)? x[i] : fill]];
 
-module echo_matrix(M,description,sig=4,eps=1e-9)
-{
-  dummy = echo_matrix(M,description,sig,eps);
-}
 
 
 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
diff --git a/attachments.scad b/attachments.scad
index 19770e2..41785b9 100644
--- a/attachments.scad
+++ b/attachments.scad
@@ -1622,8 +1622,8 @@ function _find_anchor(anchor, geom) =
                 for (face = faces)
                     let(
                         verts = select(rpts, face),
-                        ys = columns(verts,1),
-                        zs = columns(verts,2)
+                        ys = column(verts,1),
+                        zs = column(verts,2)
                     )
                     if (max(ys) >= -eps && max(zs) >= -eps &&
                         min(ys) <=  eps &&  min(zs) <=  eps)
@@ -1640,7 +1640,7 @@ function _find_anchor(anchor, geom) =
         )
         assert(len(hits)>0, "Anchor vector does not intersect with the shape.  Attachment failed.")
         let(
-            furthest = max_index(columns(hits,0)),
+            furthest = max_index(column(hits,0)),
             dist = hits[furthest][0],
             pos = hits[furthest][2],
             hitnorms = [for (hit = hits) if (approx(hit[0],dist,eps=eps)) hit[1]],
@@ -1665,7 +1665,7 @@ function _find_anchor(anchor, geom) =
         ) vnf==EMPTY_VNF? [anchor, [0,0,0], unit(anchor), 0] :
         let(
             rpts = apply(rot(from=anchor, to=RIGHT) * move(point3d(-cp)), vnf[0]),
-            maxx = max(columns(rpts,0)),
+            maxx = max(column(rpts,0)),
             idxs = [for (i = idx(rpts)) if (approx(rpts[i].x, maxx)) i],
             avep = sum(select(rpts,idxs))/len(idxs),
             mpt = approx(point2d(anchor),[0,0])? [maxx,0,0] : avep,
@@ -1716,7 +1716,7 @@ function _find_anchor(anchor, geom) =
                 )
                 if(!is_undef(isect) && !approx(isect,t[0])) [norm(isect), isect, n2]
             ],
-            maxidx = max_index(columns(isects,0)),
+            maxidx = max_index(column(isects,0)),
             isect = isects[maxidx],
             pos = point2d(cp) + isect[1],
             vec = unit(isect[2],[0,1])
@@ -1728,7 +1728,7 @@ function _find_anchor(anchor, geom) =
             anchor = point2d(anchor),
             m = rot(from=anchor, to=RIGHT) * move(-[cp.x, cp.y, 0]),
             rpts = apply(m, flatten(rgn)),
-            maxx = max(columns(rpts,0)),
+            maxx = max(column(rpts,0)),
             idxs = [for (i = idx(rpts)) if (approx(rpts[i].x, maxx)) i],
             miny = min([for (i=idxs) rpts[i].y]),
             maxy = max([for (i=idxs) rpts[i].y]),
@@ -1757,7 +1757,7 @@ function _find_anchor(anchor, geom) =
                 if(!is_undef(isect) && !approx(isect,t[0]))
                 [norm(isect), isect, n2]
             ],
-            maxidx = max_index(columns(isects,0)),
+            maxidx = max_index(column(isects,0)),
             isect = isects[maxidx],
             pos = point3d(cp) + point3d(isect[1]) + unit([0,0,anchor.z],CENTER)*l/2,
             xyvec = unit(isect[2],[0,1]),
@@ -1775,7 +1775,7 @@ function _find_anchor(anchor, geom) =
                 rot(from=xyanch, to=RIGHT, planar=true)
             ) * move(-[cp.x, cp.y]),
             rpts = apply(m, flatten(rgn)),
-            maxx = max(columns(rpts,0)),
+            maxx = max(column(rpts,0)),
             idxs = [for (i = idx(rpts)) if (approx(rpts[i].x, maxx)) i],
             ys = [for (i=idxs) rpts[i].y],
             midy = (min(ys)+max(ys))/2,
diff --git a/beziers.scad b/beziers.scad
index f4b034f..1ad0cd8 100644
--- a/beziers.scad
+++ b/beziers.scad
@@ -1043,9 +1043,9 @@ function bezier_patch_points(patch, u, v) =
     assert(is_num(u) || !is_undef(u[0]))
     assert(is_num(v) || !is_undef(v[0]))
     let(
-        vbezes = [for (i = idx(patch[0])) bezier_points(columns(patch,i), is_num(u)? [u] : u)]
+        vbezes = [for (i = idx(patch[0])) bezier_points(column(patch,i), is_num(u)? [u] : u)]
     )
-    [for (i = idx(vbezes[0])) bezier_points(columns(vbezes,i), is_num(v)? [v] : v)];
+    [for (i = idx(vbezes[0])) bezier_points(column(vbezes,i), is_num(v)? [v] : v)];
 
 
 // Function: bezier_triangle_point()
@@ -1358,7 +1358,7 @@ function bezier_patch_degenerate(patch, splinesteps=16, reverse=false, return_ed
     all(row_degen) && all(col_degen) ?  // fully degenerate case
         [EMPTY_VNF, repeat([patch[0][0]],4)] :
     all(row_degen) ?                         // degenerate to a line (top to bottom)
-        let(pts = bezier_points(columns(patch,0), samplepts))
+        let(pts = bezier_points(column(patch,0), samplepts))
         [EMPTY_VNF, [pts,pts,[pts[0]],[last(pts)]]] :
     all(col_degen) ?                         // degenerate to a line (left to right)
         let(pts = bezier_points(patch[0], samplepts))
@@ -1367,7 +1367,7 @@ function bezier_patch_degenerate(patch, splinesteps=16, reverse=false, return_ed
        let(pts = bezier_patch_points(patch, samplepts, samplepts))
        [
         vnf_vertex_array(pts, reverse=!reverse),
-        [columns(pts,0), columns(pts,len(pts)-1), pts[0], last(pts)]
+        [column(pts,0), column(pts,len(pts)-1), pts[0], last(pts)]
        ] :
     top_degen && bot_degen ?
        let(
@@ -1376,17 +1376,17 @@ function bezier_patch_degenerate(patch, splinesteps=16, reverse=false, return_ed
                         if (splinesteps%2==0) splinesteps+1,
                         each reverse(list([3:2:splinesteps]))
                        ],
-            bpatch = [for(i=[0:1:len(patch[0])-1]) bezier_points(columns(patch,i), samplepts)],
+            bpatch = [for(i=[0:1:len(patch[0])-1]) bezier_points(column(patch,i), samplepts)],
             pts = [
                   [bpatch[0][0]],
-                  for(j=[0:splinesteps-2]) bezier_points(columns(bpatch,j+1), lerpn(0,1,rowcount[j])),
+                  for(j=[0:splinesteps-2]) bezier_points(column(bpatch,j+1), lerpn(0,1,rowcount[j])),
                   [last(bpatch[0])]
                   ],
             vnf = vnf_tri_array(pts, reverse=!reverse)
          ) [
             vnf,
             [
-             columns(pts,0),
+             column(pts,0),
              [for(row=pts) last(row)],
              pts[0],
              last(pts),
@@ -1405,16 +1405,16 @@ function bezier_patch_degenerate(patch, splinesteps=16, reverse=false, return_ed
            full_degen = len(patch)>=4 && all(select(row_degen,1,ceil(len(patch)/2-1))),
            rowmax = full_degen ? count(splinesteps+1) :
                                  [for(j=[0:splinesteps]) j<=splinesteps/2 ? 2*j : splinesteps],
-           bpatch = [for(i=[0:1:len(patch[0])-1]) bezier_points(columns(patch,i), samplepts)],
+           bpatch = [for(i=[0:1:len(patch[0])-1]) bezier_points(column(patch,i), samplepts)],
            pts = [
                   [bpatch[0][0]],
-                  for(j=[1:splinesteps]) bezier_points(columns(bpatch,j), lerpn(0,1,rowmax[j]+1))
+                  for(j=[1:splinesteps]) bezier_points(column(bpatch,j), lerpn(0,1,rowmax[j]+1))
                  ],
            vnf = vnf_tri_array(pts, reverse=!reverse)
         ) [
             vnf,
             [
-             columns(pts,0),
+             column(pts,0),
              [for(row=pts) last(row)],
              pts[0],
              last(pts),
diff --git a/linalg.scad b/linalg.scad
new file mode 100644
index 0000000..235e003
--- /dev/null
+++ b/linalg.scad
@@ -0,0 +1,713 @@
+//////////////////////////////////////////////////////////////////////
+// LibFile: linalg.scad
+//   This file provides linear algebra, with support for matrix construction,
+//   solutions to linear systems of equations, QR and Cholesky factorizations, and
+//   matrix inverse.  
+// Includes:
+//   include <BOSL2/std.scad>
+//////////////////////////////////////////////////////////////////////
+
+// Section: Matrices
+//   The matrix, a rectangular array of numbers which represents a linear transformation,
+//   is the fundamental object in linear algebra.  In OpenSCAD a matrix is a list of lists of numbers
+//   with a rectangular structure.  Because OpenSCAD treats all data the same, most of the functions that
+//   index matrices or construct them will work on matrices (lists of lists) whose elements are not numbers but may be
+//   arbitrary data: strings, booleans, or even other lists.  It may even be acceptable in some cases if the structure is non-rectangular.
+//   Of course, linear algebra computations and solutions require true matrices with rectangular structure, where all the entries are
+//   finite numbers.
+//   .
+//   Matrices in OpenSCAD are lists of row vectors.  However, a potential source of confusion is that OpenSCAD
+//   treats vectors as either column vectors or row vectors as demanded by
+//   context.  Thus both `v*M` and `M*v` are valid if `M` is square and `v` has the right length.  If you want to multiply
+//   `M` on the left by `v` and `w` you can do this with `[v,w]*M` but if you want to multiply on the right side with `v` and `w` as
+//   column vectors, you now need to use {{transpose()}} because OpenSCAD doesn't adjust matrices
+//   contextually:  `A=M*transpose([v,w])`.  The solutions are now columns of A and you must extract
+//   them with {{column()}} or take the transpose of `A`.  
+
+
+// Section: Matrix testing and display
+
+// Function: is_matrix()
+// Usage:
+//   test = is_matrix(A, [m], [n], [square])
+// Description:
+//   Returns true if A is a numeric matrix of height m and width n with finite entries.  If m or n
+//   are omitted or set to undef then true is returned for any positive dimension.
+// Arguments:
+//   A = The matrix to test.
+//   m = If given, requires the matrix to have this height.
+//   n = Is given, requires the matrix to have this width.
+//   square = If true, matrix must have height equal to width. Default: false
+function is_matrix(A,m,n,square=false) =
+   is_list(A)
+   && (( is_undef(m) && len(A) ) || len(A)==m)
+   && (!square || len(A) == len(A[0]))
+   && is_vector(A[0],n)
+   && is_consistent(A);
+
+
+// Function: is_matrix_symmetric()
+// Usage:
+//   b = is_matrix_symmetric(A, [eps])
+// Description:
+//   Returns true if the input matrix is symmetric, meaning it approximately equals its transpose.  
+//   The matrix can have arbitrary entries.  
+// Arguments:
+//   A = matrix to test
+//   eps = epsilon for comparing equality.  Default: 1e-12
+function is_matrix_symmetric(A,eps=1e-12) =
+    approx(A,transpose(A), eps);
+
+
+// Function&Module: echo_matrix()
+// Usage:
+//    echo_matrix(M, [description=], [sig=], [eps=]);
+//    dummy = echo_matrix(M, [description=], [sig=], [eps=]),
+// Description:
+//    Display a numerical matrix in a readable columnar format with `sig` significant
+//    digits.  Values smaller than eps display as zero.  If you give a description
+//    it is displayed at the top.  
+function echo_matrix(M,description,sig=4,eps=1e-9) =
+  let(
+      horiz_line = chr(8213),
+      matstr = matrix_strings(M,sig=sig,eps=eps),
+      separator = str_join(repeat(horiz_line,10)),
+      dummy=echo(str(separator,"  ",is_def(description) ? description : ""))
+            [for(row=matstr) echo(row)]
+  )
+  echo(separator);
+
+module echo_matrix(M,description,sig=4,eps=1e-9)
+{
+  dummy = echo_matrix(M,description,sig,eps);
+}
+
+
+// Section: Matrix indexing
+
+// Function: column()
+// Usage:
+//   list = column(M, i);
+// Topics: Array Handling, List Handling
+// See Also: select(), slice()
+// Description:
+//   Extracts entry i from each list in M, or equivalently column i from the matrix M, and returns it as a vector.  
+//   This function will return `undef` at all entry positions indexed by i not found in M.
+// Arguments:
+//   M = The given list of lists.
+//   idx = The index, list of indices, or range of indices to fetch.
+// Example:
+//   M = [[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]];
+//   a = column(M,2);      // Returns [3, 7, 11, 15]
+//   b = column(M,0);      // Returns [1, 5, 9, 13]
+//   N = [ [1,2], [3], [4,5], [6,7,8] ];
+//   c = column(N,1);      // Returns [1,undef,5,7]
+//   data = [[1,[3,4]], [3, [9,3]], [4, [3,1]]];   // Matrix with non-numeric entries
+//   d = column(data,0);   // Returns [1,3,4]
+//   e = column(data,1);   // Returns [[3,4],[9,3],[3,1]]
+function column(M, i) =
+    assert( is_list(M), "The input is not a list." )
+    assert( is_int(i) && i>=0, "Invalid index")
+    [for(row=M) row[i]];
+
+
+
+// Function: submatrix()
+// Usage:
+//   mat = submatrix(M, idx1, idx2);
+// Topics: Matrices, Array Handling
+// See Also: column(), block_matrix(), submatrix_set()
+// 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 columns 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] ] ];
+
+
+// Section: Matrix construction and modification
+
+// Function: ident()
+// Usage:
+//   mat = ident(n);
+// Topics: Affine, Matrices
+// Description:
+//   Create an `n` by `n` square identity matrix.
+// Arguments:
+//   n = The size of the identity matrix square, `n` by `n`.
+// Example:
+//   mat = ident(3);
+//   // Returns:
+//   //   [
+//   //     [1, 0, 0],
+//   //     [0, 1, 0],
+//   //     [0, 0, 1]
+//   //   ]
+// Example:
+//   mat = ident(4);
+//   // Returns:
+//   //   [
+//   //     [1, 0, 0, 0],
+//   //     [0, 1, 0, 0],
+//   //     [0, 0, 1, 0],
+//   //     [0, 0, 0, 1]
+//   //   ]
+function ident(n) = [
+    for (i = [0:1:n-1]) [
+        for (j = [0:1:n-1]) (i==j)? 1 : 0
+    ]
+];
+
+
+// Function: diagonal_matrix()
+// Usage:
+//   mat = diagonal_matrix(diag, [offdiag]);
+// Topics: Matrices, Array Handling
+// See Also: column(), submatrix()
+// Description:
+//   Creates a square matrix with the items in the list `diag` on
+//   its diagonal.  The off diagonal entries are set to offdiag,
+//   which is zero by default. 
+// Arguments:
+//   diag = A list of items to put in the diagnal cells of the matrix.
+//   offdiag = Value to put in non-diagonal matrix cells.
+function diagonal_matrix(diag, offdiag=0) =
+  assert(is_list(diag) && len(diag)>0)
+  [for(i=[0:1:len(diag)-1]) [for(j=[0:len(diag)-1]) i==j?diag[i] : offdiag]];
+
+
+// Function: transpose()
+// Usage:
+//    M = transpose(M, [reverse]);
+// Topics: Matrices, Array Handling
+// See Also: submatrix(), block_matrix(), hstack(), flatten()
+// Description:
+//    Returns the transpose of the given input matrix.  The input can be a matrix with arbitrary entries or
+//    a numerical vector.  If you give a vector then transpose returns it unchanged.  
+//    When reverse=true, the transpose is done across to the secondary diagonal.  (See example below.)
+//    By default, reverse=false.
+// Example:
+//   M = [
+//       [1, 2, 3],
+//       [4, 5, 6],
+//       [7, 8, 9]
+//   ];
+//   t = transpose(M);
+//   // Returns:
+//   // [
+//   //     [1, 4, 7], 
+//   //     [2, 5, 8], 
+//   //     [3, 6, 9]
+//   // ]
+// Example:
+//   M = [
+//       [1, 2, 3], 
+//       [4, 5, 6]
+//   ];
+//   t = transpose(M);
+//   // Returns:
+//   // [
+//   //     [1, 4],
+//   //     [2, 5],
+//   //     [3, 6],
+//   // ]
+// Example:
+//   M = [
+//       [1, 2, 3], 
+//       [4, 5, 6], 
+//       [7, 8, 9]
+//   ];
+//   t = transpose(M, reverse=true);
+//   // Returns:
+//   // [
+//   //  [9, 6, 3],
+//   //  [8, 5, 2],
+//   //  [7, 4, 1]
+//   // ]
+// Example: Transpose on a list of numbers returns the list unchanged
+//   transpose([3,4,5]);  // Returns: [3,4,5]
+// Example: Transpose on non-numeric input
+//   arr = [
+//       [  "a",  "b", "c"],
+//       [  "d",  "e", "f"],
+//       [[1,2],[3,4],[5,6]]
+//   ];
+//   t = transpose(arr);
+//   // Returns:
+//   // [
+//   //     ["a", "d", [1,2]],
+//   //     ["b", "e", [3,4]],
+//   //     ["c", "f", [5,6]],
+//   // ]
+
+function transpose(M, reverse=false) =
+    assert( is_list(M) && len(M)>0, "Input to transpose must be a nonempty list.")
+    is_list(M[0])
+    ?   let( len0 = len(M[0]) )
+        assert([for(a=M) if(!is_list(a) || len(a)!=len0) 1 ]==[], "Input to transpose has inconsistent row lengths." )
+        reverse
+        ? [for (i=[0:1:len0-1]) 
+              [ for (j=[0:1:len(M)-1]) M[len(M)-1-j][len0-1-i] ] ] 
+        : [for (i=[0:1:len0-1]) 
+              [ for (j=[0:1:len(M)-1]) M[j][i] ] ] 
+    :  assert( is_vector(M), "Input to transpose must be a vector or list of lists.")
+           M;
+
+
+// Function: outer_product()
+// Usage:
+//   x = outer_product(u,v);
+// Description:
+//   Compute the outer product of two vectors, a matrix.  
+// Usage:
+//   M = outer_product(u,v);
+function outer_product(u,v) =
+  assert(is_vector(u) && is_vector(v), "The inputs must be vectors.")
+  [for(ui=u) ui*v];
+
+// Function: submatrix_set()
+// Usage:
+//   mat = submatrix_set(M, A, [m], [n]);
+// Topics: Matrices, Array Handling
+// See Also: column(), submatrix()
+// Description:
+//   Sets a submatrix of M equal to the matrix A.  By default the top left corner of M is set to A, but
+//   you can specify offset coordinates m and n.  If A (as adjusted by m and n) extends beyond the bounds
+//   of M then the extra entries are ignored.  You can pass in A=[[]], a null matrix, and M will be
+//   returned unchanged.  This function works on arbitrary lists of lists and the input M need not be rectangular in shape.  
+// Arguments:
+//   M = Original matrix.
+//   A = Submatrix of new values to write into M
+//   m = Row number of upper-left corner to place A at.  Default: 0
+//   n = Column number of upper-left corner to place A at.  Default: 0 
+function submatrix_set(M,A,m=0,n=0) =
+    assert(is_list(M))
+    assert(is_list(A))
+    assert(is_int(m))
+    assert(is_int(n))
+    let( badrows = [for(i=idx(A)) if (!is_list(A[i])) i])
+    assert(badrows==[], str("Input submatrix malformed rows: ",badrows))
+    [for(i=[0:1:len(M)-1])
+        assert(is_list(M[i]), str("Row ",i," of input matrix is not a list"))
+        [for(j=[0:1:len(M[i])-1]) 
+            i>=m && i <len(A)+m && j>=n && j<len(A[0])+n ? A[i-m][j-n] : M[i][j]]];
+
+
+// Function: hstack()
+// Usage: 
+//   A = hstack(M1, M2)
+//   A = hstack(M1, M2, M3)
+//   A = hstack([M1, M2, M3, ...])
+// Topics: Matrices, Array Handling
+// See Also: column(), submatrix(), block_matrix()
+// Description:
+//   Constructs a matrix by horizontally "stacking" together compatible matrices or vectors.  Vectors are treated as columsn in the stack.
+//   This command is the inverse of `column`.  Note: strings given in vectors are broken apart into lists of characters.  Strings given
+//   in matrices are preserved as strings.  If you need to combine vectors of strings use array_group as shown below to convert the
+//   vector into a column matrix.  Also note that vertical stacking can be done directly with concat.  
+// Arguments:
+//   M1 = If given with other arguments, the first matrix (or vector) to stack.  If given alone, a list of matrices/vectors to stack. 
+//   M2 = Second matrix/vector to stack
+//   M3 = Third matrix/vector to stack.
+// Example:
+//   M = ident(3);
+//   v1 = [2,3,4];
+//   v2 = [5,6,7];
+//   v3 = [8,9,10];
+//   a = hstack(v1,v2);     // Returns [[2, 5], [3, 6], [4, 7]]
+//   b = hstack(v1,v2,v3);  // Returns [[2, 5,  8],
+//                          //          [3, 6,  9],
+//                          //          [4, 7, 10]]
+//   c = hstack([M,v1,M]);  // Returns [[1, 0, 0, 2, 1, 0, 0],
+//                          //          [0, 1, 0, 3, 0, 1, 0],
+//                          //          [0, 0, 1, 4, 0, 0, 1]]
+//   d = hstack(column(M,0), submatrix(M,idx(M),[1 2]));  // Returns M
+//   strvec = ["one","two"];
+//   strmat = [["three","four"], ["five","six"]];
+//   e = hstack(strvec,strvec); // Returns [["o", "n", "e", "o", "n", "e"],
+//                              //          ["t", "w", "o", "t", "w", "o"]]
+//   f = hstack(array_group(strvec,1), array_group(strvec,1));
+//                              // Returns [["one", "one"],
+//                              //          ["two", "two"]]
+//   g = hstack(strmat,strmat); //  Returns: [["three", "four", "three", "four"],
+//                              //            [ "five",  "six",  "five",  "six"]]
+function hstack(M1, M2, M3) =
+    (M3!=undef)? hstack([M1,M2,M3]) : 
+    (M2!=undef)? hstack([M1,M2]) :
+    assert(all([for(v=M1) is_list(v)]), "One of the inputs to hstack is not a list")
+    let(
+        minlen = min_length(M1),
+        maxlen = max_length(M1)
+    )
+    assert(minlen==maxlen, "Input vectors to hstack must have the same length")
+    [for(row=[0:1:minlen-1])
+        [for(matrix=M1)
+           each matrix[row]
+        ]
+    ];
+
+
+// Function: block_matrix()
+// Usage:
+//    bmat = block_matrix([[M11, M12,...],[M21, M22,...], ... ]);
+// Topics: Matrices, Array Handling
+// See Also: column(), submatrix()
+// Description:
+//    Create a block matrix by supplying a matrix of matrices, which will
+//    be combined into one unified matrix.  Every matrix in one row
+//    must have the same height, and the combined width of the matrices
+//    in each row must be equal. Strings will stay strings. 
+// Example:
+//  A = [[1,2],
+//       [3,4]];
+//  B = ident(2);
+//  C = block_matrix([[A,B],[B,A],[A,B]]);
+//      // Returns:
+//      //        [[1, 2, 1, 0],
+//      //         [3, 4, 0, 1],
+//      //         [1, 0, 1, 2],
+//      //         [0, 1, 3, 4],
+//      //         [1, 2, 1, 0],
+//      //         [3, 4, 0, 1]]);
+//  D = block_matrix([[A,B], ident(4)]);
+//      // Returns:
+//      //        [[1, 2, 1, 0],
+//      //         [3, 4, 0, 1],
+//      //         [1, 0, 0, 0],
+//      //         [0, 1, 0, 0],
+//      //         [0, 0, 1, 0],
+//      //         [0, 0, 0, 1]]);
+//  E = [["one", "two"], [3,4]];
+//  F = block_matrix([[E,E]]);
+//      // Returns:
+//      //        [["one", "two", "one", "two"],
+//      //         [    3,     4,     3,     4]]
+function block_matrix(M) =
+    let(
+        bigM = [for(bigrow = M) each hstack(bigrow)],
+        len0 = len(bigM[0]),
+        badrows = [for(row=bigM) if (len(row)!=len0) 1]
+    )
+    assert(badrows==[], "Inconsistent or invalid input")
+    bigM;
+
+
+// Section: Solving Linear Equations and Matrix Factorizations
+
+// Function: linear_solve()
+// Usage:
+//   solv = linear_solve(A,b)
+// Description:
+//   Solves the linear system Ax=b.  If `A` is square and non-singular the unique solution is returned.  If `A` is overdetermined
+//   the least squares solution is returned. If `A` is underdetermined, the minimal norm solution is returned.
+//   If `A` is rank deficient or singular then linear_solve returns `[]`.  If `b` is a matrix that is compatible with `A`
+//   then the problem is solved for the matrix valued right hand side and a matrix is returned.  Note that if you 
+//   want to solve Ax=b1 and Ax=b2 that you need to form the matrix `transpose([b1,b2])` for the right hand side and then
+//   transpose the returned value.
+function linear_solve(A,b,pivot=true) =
+    assert(is_matrix(A), "Input should be a matrix.")
+    let(
+        m = len(A),
+        n = len(A[0])
+    )
+    assert(is_vector(b,m) || is_matrix(b,m),"Invalid right hand side or incompatible with the matrix")
+    let (
+        qr = m<n? qr_factor(transpose(A),pivot) : qr_factor(A,pivot),
+        maxdim = max(n,m),
+        mindim = min(n,m),
+        Q = submatrix(qr[0],[0:maxdim-1], [0:mindim-1]),
+        R = submatrix(qr[1],[0:mindim-1], [0:mindim-1]),
+        P = qr[2],
+        zeros = [for(i=[0:mindim-1]) if (approx(R[i][i],0)) i]
+    )
+    zeros != [] ? [] :
+    m<n ? Q*back_substitute(R,transpose(P)*b,transpose=true) // Too messy to avoid input checks here
+        : P*_back_substitute(R, transpose(Q)*b);             // Calling internal version skips input checks
+
+// Function: matrix_inverse()
+// Usage:
+//    mat = matrix_inverse(A)
+// Description:
+//    Compute the matrix inverse of the square matrix `A`.  If `A` is singular, returns `undef`.
+//    Note that if you just want to solve a linear system of equations you should NOT use this function.
+//    Instead use {{linear_solve()}}, or use {{qr_factor()}}.  The computation
+//    will be faster and more accurate.  
+function matrix_inverse(A) =
+    assert(is_matrix(A) && len(A)==len(A[0]),"Input to matrix_inverse() must be a square matrix")
+    linear_solve(A,ident(len(A)));
+
+
+// Function: rot_inverse()
+// Usage:
+//   B = rot_inverse(A)
+// Description:
+//   Inverts a 2d (3x3) or 3d (4x4) rotation matrix.  The matrix can be a rotation around any center,
+//   so it may include a translation.  
+function rot_inverse(T) =
+    assert(is_matrix(T,square=true),"Matrix must be square")
+    let( n = len(T))
+    assert(n==3 || n==4, "Matrix must be 3x3 or 4x4")
+    let(
+        rotpart =  [for(i=[0:n-2]) [for(j=[0:n-2]) T[j][i]]],
+        transpart = [for(row=[0:n-2]) T[row][n-1]]
+    )
+    assert(approx(determinant(T),1),"Matrix is not a rotation")
+    concat(hstack(rotpart, -rotpart*transpart),[[for(i=[2:n]) 0, 1]]);
+
+
+
+
+// Function: null_space()
+// Usage:
+//   x = null_space(A)
+// Description:
+//   Returns an orthonormal basis for the null space of `A`, namely the vectors {x} such that Ax=0.
+//   If the null space is just the origin then returns an empty list. 
+function null_space(A,eps=1e-12) =
+    assert(is_matrix(A))
+    let(
+        Q_R = qr_factor(transpose(A),pivot=true),
+        R = Q_R[1],
+        zrows = [for(i=idx(R)) if (all_zero(R[i],eps)) i]
+    )
+    len(zrows)==0 ? [] :
+    select(transpose(Q_R[0]), zrows);
+
+// Function: qr_factor()
+// Usage:
+//   qr = qr_factor(A,[pivot]);
+// Description:
+//   Calculates the QR factorization of the input matrix A and returns it as the list [Q,R,P].  This factorization can be
+//   used to solve linear systems of equations.  The factorization is `A = Q*R*transpose(P)`.  If pivot is false (the default)
+//   then P is the identity matrix and A = Q*R.  If pivot is true then column pivoting results in an R matrix where the diagonal
+//   is non-decreasing.  The use of pivoting is supposed to increase accuracy for poorly conditioned problems, and is necessary
+//   for rank estimation or computation of the null space, but it may be slower.  
+function qr_factor(A, pivot=false) =
+    assert(is_matrix(A), "Input must be a matrix." )
+    let(
+        m = len(A),
+        n = len(A[0])
+    )
+    let(
+        qr = _qr_factor(A, Q=ident(m),P=ident(n), pivot=pivot, col=0, m = m, n = n),
+        Rzero = let( R = qr[1]) [
+            for(i=[0:m-1]) [
+                let( ri = R[i] )
+                for(j=[0:n-1]) i>j ? 0 : ri[j]
+            ]
+        ]
+    ) [qr[0], Rzero, qr[2]];
+
+function _qr_factor(A,Q,P, pivot, col, m, n) =
+    col >= min(m-1,n) ? [Q,A,P] :
+    let(
+        swap = !pivot ? 1
+             : _swap_matrix(n,col,col+max_index([for(i=[col:n-1]) sqr([for(j=[col:m-1]) A[j][i]])])),
+        A = pivot ? A*swap : A,
+        x = [for(i=[col:1:m-1]) A[i][col]],
+        alpha = (x[0]<=0 ? 1 : -1) * norm(x),
+        u = x - concat([alpha],repeat(0,m-1)),
+        v = alpha==0 ? u : u / norm(u),
+        Qc = ident(len(x)) - 2*outer_product(v,v),
+        Qf = [for(i=[0:m-1]) [for(j=[0:m-1]) i<col || j<col ? (i==j ? 1 : 0) : Qc[i-col][j-col]]]
+    )
+    _qr_factor(Qf*A, Q*Qf, P*swap, pivot, col+1, m, n);
+
+// Produces an n x n matrix that swaps column i and j (when multiplied on the right)
+function _swap_matrix(n,i,j) =
+  assert(i<n && j<n && i>=0 && j>=0, "Swap indices out of bounds")
+  [for(y=[0:n-1]) [for (x=[0:n-1])
+     x==i ? (y==j ? 1 : 0)
+   : x==j ? (y==i ? 1 : 0)
+   : x==y ? 1 : 0]];
+
+
+
+// Function: back_substitute()
+// Usage:
+//   x = back_substitute(R, b, [transpose]);
+// Description:
+//   Solves the problem Rx=b where R is an upper triangular square matrix.  The lower triangular entries of R are
+//   ignored.  If transpose==true then instead solve transpose(R)*x=b.
+//   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 [].  
+function back_substitute(R, b, transpose = false) =
+    assert(is_matrix(R, square=true))
+    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(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(ind = n - len(x) - 1)
+      R[ind][ind] == 0 ? []
+    : let(
+          newvalue = len(x)==0
+            ? b[ind]/R[ind][ind]
+            : (b[ind]-list_tail(R[ind],ind+1) * x)/R[ind][ind]
+      )
+      _back_substitute(R, b, concat([newvalue],x));
+
+
+
+// Function: cholesky()
+// Usage:
+//   L = cholesky(A);
+// Description:
+//   Compute the cholesky factor, L, of the symmetric positive definite matrix A.
+//   The matrix L is lower triangular and `L * transpose(L) = A`.  If the A is
+//   not symmetric then an error is displayed.  If the matrix is symmetric but
+//   not positive definite then undef is returned.  
+function cholesky(A) =
+  assert(is_matrix(A,square=true),"A must be a square matrix")
+  assert(is_matrix_symmetric(A),"Cholesky factorization requires a symmetric matrix")
+  _cholesky(A,ident(len(A)), len(A));
+
+function _cholesky(A,L,n) = 
+    A[0][0]<0 ? undef :     // Matrix not positive definite
+    len(A) == 1 ? submatrix_set(L,[[sqrt(A[0][0])]], n-1,n-1):
+    let(
+        i = n+1-len(A)
+    )
+    let(
+        sqrtAii = sqrt(A[0][0]),
+        Lnext = [for(j=[0:n-1])
+                  [for(k=[0:n-1])
+                      j<i-1 || k<i-1 ?  (j==k ? 1 : 0)
+                     : j==i-1 && k==i-1 ? sqrtAii
+                     : j==i-1 ? 0
+                     : k==i-1 ? A[j-(i-1)][0]/sqrtAii
+                     : j==k ? 1 : 0]],
+        Anext = submatrix(A,[1:n-1], [1:n-1]) - outer_product(list_tail(A[0]), list_tail(A[0]))/A[0][0]
+    )
+    _cholesky(Anext,L*Lnext,n);
+
+
+// Section: Matrix Properties: Determinants, Norm, Trace
+
+// Function: det2()
+// Usage:
+//   d = det2(M);
+// Description:
+//   Rturns the determinant for the given 2x2 matrix.
+// Arguments:
+//   M = The 2x2 matrix to get the determinant of.
+// Example:
+//   M = [ [6,-2], [1,8] ];
+//   det = det2(M);  // Returns: 50
+function det2(M) = 
+    assert(is_def(M) && M*0==[[0,0],[0,0]], "Expected square matrix (2x2)")
+    cross(M[0],M[1]);
+
+
+// Function: det3()
+// Usage:
+//   d = det3(M);
+// Description:
+//   Returns the determinant for the given 3x3 matrix.
+// Arguments:
+//   M = The 3x3 square matrix to get the determinant of.
+// Example:
+//   M = [ [6,4,-2], [1,-2,8], [1,5,7] ];
+//   det = det3(M);  // Returns: -334
+function det3(M) =
+    assert(is_def(M) && M*0==[[0,0,0],[0,0,0],[0,0,0]], "Expected square matrix (3x3).")
+    M[0][0] * (M[1][1]*M[2][2]-M[2][1]*M[1][2]) -
+    M[1][0] * (M[0][1]*M[2][2]-M[2][1]*M[0][2]) +
+    M[2][0] * (M[0][1]*M[1][2]-M[1][1]*M[0][2]);
+
+// Function: det4()
+// Usage:
+//   d = det4(M);
+// Description:
+//   Returns the determinant for the given 4x4 matrix.
+// Arguments:
+//   M = The 4x4 square matrix to get the determinant of.
+// Example:
+//   M = [ [6,4,-2,1], [1,-2,8,-3], [1,5,7,4], [2,3,4,7] ];
+//   det = det4(M);  // Returns: -1773
+function det4(M) =
+    assert(is_def(M) && M*0==[[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]], "Expected square matrix (4x4).")
+    M[0][0]*M[1][1]*M[2][2]*M[3][3] + M[0][0]*M[1][2]*M[2][3]*M[3][1] + M[0][0]*M[1][3]*M[2][1]*M[3][2]
+    + M[0][1]*M[1][0]*M[2][3]*M[3][2] + M[0][1]*M[1][2]*M[2][0]*M[3][3] + M[0][1]*M[1][3]*M[2][2]*M[3][0]
+    + M[0][2]*M[1][0]*M[2][1]*M[3][3] + M[0][2]*M[1][1]*M[2][3]*M[3][0] + M[0][2]*M[1][3]*M[2][0]*M[3][1]
+    + M[0][3]*M[1][0]*M[2][2]*M[3][1] + M[0][3]*M[1][1]*M[2][0]*M[3][2] + M[0][3]*M[1][2]*M[2][1]*M[3][0]
+    - M[0][0]*M[1][1]*M[2][3]*M[3][2] - M[0][0]*M[1][2]*M[2][1]*M[3][3] - M[0][0]*M[1][3]*M[2][2]*M[3][1]
+    - M[0][1]*M[1][0]*M[2][2]*M[3][3] - M[0][1]*M[1][2]*M[2][3]*M[3][0] - M[0][1]*M[1][3]*M[2][0]*M[3][2]
+    - M[0][2]*M[1][0]*M[2][3]*M[3][1] - M[0][2]*M[1][1]*M[2][0]*M[3][3] - M[0][2]*M[1][3]*M[2][1]*M[3][0]
+    - M[0][3]*M[1][0]*M[2][1]*M[3][2] - M[0][3]*M[1][1]*M[2][2]*M[3][0] - M[0][3]*M[1][2]*M[2][0]*M[3][1];
+
+// Function: determinant()
+// Usage:
+//   d = determinant(M);
+// Description:
+//   Returns the determinant for the given square matrix.
+// Arguments:
+//   M = The NxN square matrix to get the determinant of.
+// Example:
+//   M = [ [6,4,-2,9], [1,-2,8,3], [1,5,7,6], [4,2,5,1] ];
+//   det = determinant(M);  // Returns: 2267
+function determinant(M) =
+    assert(is_list(M), "Input must be a square matrix." )  
+    len(M)==1? M[0][0] :
+    len(M)==2? det2(M) :
+    len(M)==3? det3(M) :
+    len(M)==4? det4(M) :
+    assert(is_matrix(M, square=true), "Input must be a square matrix." )    
+    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]
+                        ]
+                    ]
+                )
+        ]
+    );
+
+
+// Function: norm_fro()
+// Usage:
+//    norm_fro(A)
+// Description:
+//    Computes frobenius norm of input matrix.  The frobenius norm is the square root of the sum of the
+//    squares of all of the entries of the matrix.  On vectors it is the same as the usual 2-norm.
+//    This is an easily computed norm that is convenient for comparing two matrices.  
+function norm_fro(A) =
+    assert(is_matrix(A) || is_vector(A))
+    norm(flatten(A));
+
+
+// Function: matrix_trace()
+// Usage:
+//   matrix_trace(M)
+// Description:
+//   Computes the trace of a square matrix, the sum of the entries on the diagonal.  
+function matrix_trace(M) =
+   assert(is_matrix(M,square=true), "Input to trace must be a square matrix")
+   [for(i=[0:1:len(M)-1])1] * [for(i=[0:1:len(M)-1]) M[i][i]];
diff --git a/math.scad b/math.scad
index a80aa82..c846756 100644
--- a/math.scad
+++ b/math.scad
@@ -5,7 +5,6 @@
 //   include <BOSL2/std.scad>
 //////////////////////////////////////////////////////////////////////
 
-
 // Section: Math Constants
 
 // Constant: PHI
@@ -669,7 +668,7 @@ function lcm(a,b=[]) =
 function sum(v, dflt=0) =
     v==[]? dflt :
     assert(is_consistent(v), "Input to sum is non-numeric or inconsistent")
-    is_vector(v) || is_matrix(v) ? [for(i=v) 1]*v :
+    is_finite(v[0]) || is_vector(v[0]) ? [for(i=v) 1]*v :
     _sum(v,v[0]*0);
 
 function _sum(v,_total,_i=0) = _i>=len(v) ? _total : _sum(v,_total+v[_i], _i+1);
@@ -799,18 +798,6 @@ function _cumprod_vec(v,_i=0,_acc=[]) =
     );
 
 
-// Function: outer_product()
-// Usage:
-//   x = outer_product(u,v);
-// Description:
-//   Compute the outer product of two vectors, a matrix.  
-// Usage:
-//   M = outer_product(u,v);
-function outer_product(u,v) =
-  assert(is_vector(u) && is_vector(v), "The inputs must be vectors.")
-  [for(ui=u) ui*v];
-
-
 // Function: mean()
 // Usage:
 //   x = mean(v);
@@ -877,335 +864,6 @@ function convolve(p,q) =
 
 
 
-// Section: Matrix math
-
-// Function: ident()
-// Usage:
-//   mat = ident(n);
-// Topics: Affine, Matrices
-// Description:
-//   Create an `n` by `n` square identity matrix.
-// Arguments:
-//   n = The size of the identity matrix square, `n` by `n`.
-// Example:
-//   mat = ident(3);
-//   // Returns:
-//   //   [
-//   //     [1, 0, 0],
-//   //     [0, 1, 0],
-//   //     [0, 0, 1]
-//   //   ]
-// Example:
-//   mat = ident(4);
-//   // Returns:
-//   //   [
-//   //     [1, 0, 0, 0],
-//   //     [0, 1, 0, 0],
-//   //     [0, 0, 1, 0],
-//   //     [0, 0, 0, 1]
-//   //   ]
-function ident(n) = [
-    for (i = [0:1:n-1]) [
-        for (j = [0:1:n-1]) (i==j)? 1 : 0
-    ]
-];
-
-
-// Function: linear_solve()
-// Usage:
-//   solv = linear_solve(A,b)
-// Description:
-//   Solves the linear system Ax=b.  If `A` is square and non-singular the unique solution is returned.  If `A` is overdetermined
-//   the least squares solution is returned. If `A` is underdetermined, the minimal norm solution is returned.
-//   If `A` is rank deficient or singular then linear_solve returns `[]`.  If `b` is a matrix that is compatible with `A`
-//   then the problem is solved for the matrix valued right hand side and a matrix is returned.  Note that if you 
-//   want to solve Ax=b1 and Ax=b2 that you need to form the matrix `transpose([b1,b2])` for the right hand side and then
-//   transpose the returned value.
-function linear_solve(A,b,pivot=true) =
-    assert(is_matrix(A), "Input should be a matrix.")
-    let(
-        m = len(A),
-        n = len(A[0])
-    )
-    assert(is_vector(b,m) || is_matrix(b,m),"Invalid right hand side or incompatible with the matrix")
-    let (
-        qr = m<n? qr_factor(transpose(A),pivot) : qr_factor(A,pivot),
-        maxdim = max(n,m),
-        mindim = min(n,m),
-        Q = submatrix(qr[0],[0:maxdim-1], [0:mindim-1]),
-        R = submatrix(qr[1],[0:mindim-1], [0:mindim-1]),
-        P = qr[2],
-        zeros = [for(i=[0:mindim-1]) if (approx(R[i][i],0)) i]
-    )
-    zeros != [] ? [] :
-    m<n ? Q*back_substitute(R,transpose(P)*b,transpose=true) // Too messy to avoid input checks here
-        : P*_back_substitute(R, transpose(Q)*b);             // Calling internal version skips input checks
-
-// Function: matrix_inverse()
-// Usage:
-//    mat = matrix_inverse(A)
-// Description:
-//    Compute the matrix inverse of the square matrix `A`.  If `A` is singular, returns `undef`.
-//    Note that if you just want to solve a linear system of equations you should NOT use this function.
-//    Instead use {{linear_solve()}}, or use {{qr_factor()}}.  The computation
-//    will be faster and more accurate.  
-function matrix_inverse(A) =
-    assert(is_matrix(A) && len(A)==len(A[0]),"Input to matrix_inverse() must be a square matrix")
-    linear_solve(A,ident(len(A)));
-
-
-// Function: rot_inverse()
-// Usage:
-//   B = rot_inverse(A)
-// Description:
-//   Inverts a 2d (3x3) or 3d (4x4) rotation matrix.  The matrix can be a rotation around any center,
-//   so it may include a translation.  
-function rot_inverse(T) =
-    assert(is_matrix(T,square=true),"Matrix must be square")
-    let( n = len(T))
-    assert(n==3 || n==4, "Matrix must be 3x3 or 4x4")
-    let(
-        rotpart =  [for(i=[0:n-2]) [for(j=[0:n-2]) T[j][i]]],
-        transpart = [for(row=[0:n-2]) T[row][n-1]]
-    )
-    assert(approx(determinant(T),1),"Matrix is not a rotation")
-    concat(hstack(rotpart, -rotpart*transpart),[[for(i=[2:n]) 0, 1]]);
-
-
-
-
-// Function: null_space()
-// Usage:
-//   x = null_space(A)
-// Description:
-//   Returns an orthonormal basis for the null space of `A`, namely the vectors {x} such that Ax=0.
-//   If the null space is just the origin then returns an empty list. 
-function null_space(A,eps=1e-12) =
-    assert(is_matrix(A))
-    let(
-        Q_R = qr_factor(transpose(A),pivot=true),
-        R = Q_R[1],
-        zrow = [for(i=idx(R)) if (all_zero(R[i],eps)) i]
-    )
-    len(zrow)==0 ? [] :
-    transpose(columns(Q_R[0],zrow));
-
-
-// Function: qr_factor()
-// Usage:
-//   qr = qr_factor(A,[pivot]);
-// Description:
-//   Calculates the QR factorization of the input matrix A and returns it as the list [Q,R,P].  This factorization can be
-//   used to solve linear systems of equations.  The factorization is `A = Q*R*transpose(P)`.  If pivot is false (the default)
-//   then P is the identity matrix and A = Q*R.  If pivot is true then column pivoting results in an R matrix where the diagonal
-//   is non-decreasing.  The use of pivoting is supposed to increase accuracy for poorly conditioned problems, and is necessary
-//   for rank estimation or computation of the null space, but it may be slower.  
-function qr_factor(A, pivot=false) =
-    assert(is_matrix(A), "Input must be a matrix." )
-    let(
-        m = len(A),
-        n = len(A[0])
-    )
-    let(
-        qr = _qr_factor(A, Q=ident(m),P=ident(n), pivot=pivot, column=0, m = m, n=n),
-        Rzero = let( R = qr[1]) [
-            for(i=[0:m-1]) [
-                let( ri = R[i] )
-                for(j=[0:n-1]) i>j ? 0 : ri[j]
-            ]
-        ]
-    ) [qr[0], Rzero, qr[2]];
-
-function _qr_factor(A,Q,P, pivot, column, m, n) =
-    column >= min(m-1,n) ? [Q,A,P] :
-    let(
-        swap = !pivot ? 1
-             : _swap_matrix(n,column,column+max_index([for(i=[column:n-1]) sqr([for(j=[column:m-1]) A[j][i]])])),
-        A = pivot ? A*swap : A,
-        x = [for(i=[column:1:m-1]) A[i][column]],
-        alpha = (x[0]<=0 ? 1 : -1) * norm(x),
-        u = x - concat([alpha],repeat(0,m-1)),
-        v = alpha==0 ? u : u / norm(u),
-        Qc = ident(len(x)) - 2*outer_product(v,v),
-        Qf = [for(i=[0:m-1]) [for(j=[0:m-1]) i<column || j<column ? (i==j ? 1 : 0) : Qc[i-column][j-column]]]
-    )
-    _qr_factor(Qf*A, Q*Qf, P*swap, pivot, column+1, m, n);
-
-// Produces an n x n matrix that swaps column i and j (when multiplied on the right)
-function _swap_matrix(n,i,j) =
-  assert(i<n && j<n && i>=0 && j>=0, "Swap indices out of bounds")
-  [for(y=[0:n-1]) [for (x=[0:n-1])
-     x==i ? (y==j ? 1 : 0)
-   : x==j ? (y==i ? 1 : 0)
-   : x==y ? 1 : 0]];
-
-
-
-// Function: back_substitute()
-// Usage:
-//   x = back_substitute(R, b, [transpose]);
-// Description:
-//   Solves the problem Rx=b where R is an upper triangular square matrix.  The lower triangular entries of R are
-//   ignored.  If transpose==true then instead solve transpose(R)*x=b.
-//   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 [].  
-function back_substitute(R, b, transpose = false) =
-    assert(is_matrix(R, square=true))
-    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(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(ind = n - len(x) - 1)
-      R[ind][ind] == 0 ? []
-    : let(
-          newvalue = len(x)==0
-            ? b[ind]/R[ind][ind]
-            : (b[ind]-list_tail(R[ind],ind+1) * x)/R[ind][ind]
-      )
-      _back_substitute(R, b, concat([newvalue],x));
-
-
-
-// Function: cholesky()
-// Usage:
-//   L = cholesky(A);
-// Description:
-//   Compute the cholesky factor, L, of the symmetric positive definite matrix A.
-//   The matrix L is lower triangular and `L * transpose(L) = A`.  If the A is
-//   not symmetric then an error is displayed.  If the matrix is symmetric but
-//   not positive definite then undef is returned.  
-function cholesky(A) =
-  assert(is_matrix(A,square=true),"A must be a square matrix")
-  assert(is_matrix_symmetric(A),"Cholesky factorization requires a symmetric matrix")
-  _cholesky(A,ident(len(A)), len(A));
-
-function _cholesky(A,L,n) = 
-    A[0][0]<0 ? undef :     // Matrix not positive definite
-    len(A) == 1 ? submatrix_set(L,[[sqrt(A[0][0])]], n-1,n-1):
-    let(
-        i = n+1-len(A)
-    )
-    let(
-        sqrtAii = sqrt(A[0][0]),
-        Lnext = [for(j=[0:n-1])
-                  [for(k=[0:n-1])
-                      j<i-1 || k<i-1 ?  (j==k ? 1 : 0)
-                     : j==i-1 && k==i-1 ? sqrtAii
-                     : j==i-1 ? 0
-                     : k==i-1 ? A[j-(i-1)][0]/sqrtAii
-                     : j==k ? 1 : 0]],
-        Anext = submatrix(A,[1:n-1], [1:n-1]) - outer_product(list_tail(A[0]), list_tail(A[0]))/A[0][0]
-    )
-    _cholesky(Anext,L*Lnext,n);
-
-
-// Function: det2()
-// Usage:
-//   d = det2(M);
-// Description:
-//   Optimized function that returns the determinant for the given 2x2 square matrix.
-// Arguments:
-//   M = The 2x2 square matrix to get the determinant of.
-// Example:
-//   M = [ [6,-2], [1,8] ];
-//   det = det2(M);  // Returns: 50
-function det2(M) = 
-    assert(is_matrix(M,2,2), "Matrix must be 2x2.")
-    M[0][0] * M[1][1] - M[0][1]*M[1][0];
-
-
-// Function: det3()
-// Usage:
-//   d = det3(M);
-// Description:
-//   Optimized function that returns the determinant for the given 3x3 square matrix.
-// Arguments:
-//   M = The 3x3 square matrix to get the determinant of.
-// Example:
-//   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 must be 3x3.")
-    M[0][0] * (M[1][1]*M[2][2]-M[2][1]*M[1][2]) -
-    M[1][0] * (M[0][1]*M[2][2]-M[2][1]*M[0][2]) +
-    M[2][0] * (M[0][1]*M[1][2]-M[1][1]*M[0][2]);
-
-
-// Function: determinant()
-// Usage:
-//   d = determinant(M);
-// Description:
-//   Returns the determinant for the given square matrix.
-// Arguments:
-//   M = The NxN square matrix to get the determinant of.
-// Example:
-//   M = [ [6,4,-2,9], [1,-2,8,3], [1,5,7,6], [4,2,5,1] ];
-//   det = determinant(M);  // Returns: 2267
-function determinant(M) =
-    assert(is_matrix(M, square=true), "Input should be a square matrix." )
-    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]
-                        ]
-                    ]
-                )
-        ]
-    );
-
-
-// Function: is_matrix()
-// Usage:
-//   test = is_matrix(A, [m], [n], [square])
-// Description:
-//   Returns true if A is a numeric matrix of height m and width n.  If m or n
-//   are omitted or set to undef then true is returned for any positive dimension.
-// Arguments:
-//   A = The matrix to test.
-//   m = Is given, requires the matrix to have the given height.
-//   n = Is given, requires the matrix to have the given width.
-//   square = If true, requires the matrix to have a width equal to its height. Default: false
-function is_matrix(A,m,n,square=false) =
-   is_list(A)
-   && (( is_undef(m) && len(A) ) || len(A)==m)
-   && (!square || len(A) == len(A[0]))
-   && is_vector(A[0],n)
-   && is_consistent(A);
-
-
-// Function: norm_fro()
-// Usage:
-//    norm_fro(A)
-// Description:
-//    Computes frobenius norm of input matrix.  The frobenius norm is the square root of the sum of the
-//    squares of all of the entries of the matrix.  On vectors it is the same as the usual 2-norm.
-//    This is an easily computed norm that is convenient for comparing two matrices.  
-function norm_fro(A) =
-    assert(is_matrix(A) || is_vector(A))
-    norm(flatten(A));
-
-
-// Function: matrix_trace()
-// Usage:
-//   matrix_trace(M)
-// Description:
-//   Computes the trace of a square matrix, the sum of the entries on the diagonal.  
-function matrix_trace(M) =
-   assert(is_matrix(M,square=true), "Input to trace must be a square matrix")
-   [for(i=[0:1:len(M)-1])1] * [for(i=[0:1:len(M)-1]) M[i][i]];
 
 
 // Section: Comparisons and Logic
diff --git a/paths.scad b/paths.scad
index 89d526d..7da908e 100644
--- a/paths.scad
+++ b/paths.scad
@@ -427,7 +427,7 @@ function resample_path(path, N, spacing, closed=false) =
        distlist = lerpn(0,length,N,false), 
        cuts = _path_cut_points(path, distlist, closed=closed)
    )
-   [ each columns(cuts,0),
+   [ each column(cuts,0),
      if (!closed) last(path)     // Then add last point here
    ];
 
@@ -1174,7 +1174,7 @@ function _assemble_path_fragments(fragments, eps=EPSILON, _finished=[]) =
    len(fragments)==0? _finished :
     let(
         minxidx = min_index([
-            for (frag=fragments) min(columns(frag,0))
+            for (frag=fragments) min(column(frag,0))
         ]),
         result_l = _assemble_a_path_from_fragments(
             fragments=fragments,
diff --git a/regions.scad b/regions.scad
index 2084efe..d9f96c2 100644
--- a/regions.scad
+++ b/regions.scad
@@ -267,7 +267,7 @@ function _region_region_intersections(region1, region2, closed1=true,closed2=tru
          cornerpts = [for(i=[0:1])
                          [for(k=vector_search(points[i],eps,points[i]))
                              each if (len(k)>1) select(ptind[i],k)]],
-         risect = [for(i=[0:1]) concat(columns(intersections,i), cornerpts[i])],
+         risect = [for(i=[0:1]) concat(column(intersections,i), cornerpts[i])],
          counts = [count(len(region1)), count(len(region2))],
          pathind = [for(i=[0:1]) search(counts[i], risect[i], 0)]
        )
diff --git a/rounding.scad b/rounding.scad
index 69e2bf1..45e77e2 100644
--- a/rounding.scad
+++ b/rounding.scad
@@ -1270,7 +1270,7 @@ module convex_offset_extrude(
         // The entry r[i] is [radius,z] for a given layer
         r = move([0,bottom_height],p=concat(
                           reverse(offsets_bot), [[0,0], [0,middle]], move([0,middle], p=offsets_top)));
-        delta = [for(val=deltas(columns(r,0))) sign(val)];
+        delta = [for(val=deltas(column(r,0))) sign(val)];
         below=[-thickness,0];
         above=[0,thickness];
            // layers is a list of pairs of the relative positions for each layer, e.g. [0,thickness]
@@ -1937,8 +1937,8 @@ function rounded_prism(bottom, top, joint_bot=0, joint_top=0, joint_sides=0, k_b
      verify_vert =
        [for(i=[0:N-1],j=[0:4])
          let(
-               vline = concat(select(columns(top_patch[i],j),2,4),
-                              select(columns(bot_patch[i],j),2,4))
+               vline = concat(select(column(top_patch[i],j),2,4),
+                              select(column(bot_patch[i],j),2,4))
              )
          if (!is_collinear(vline)) [i,j]],
      //verify horiz edges
@@ -1955,8 +1955,8 @@ function rounded_prism(bottom, top, joint_bot=0, joint_top=0, joint_sides=0, k_b
           "Roundovers interfere with each other on bottom face: either input is self intersecting or top joint length is too large")
     assert(debug || (verify_vert==[] && verify_horiz==[]), "Curvature continuity failed")
     let( 
-        vnf = vnf_merge([ each columns(top_samples,0),
-                          each columns(bot_samples,0),
+        vnf = vnf_merge([ each column(top_samples,0),
+                          each column(bot_samples,0),
                           for(pts=edge_points) vnf_vertex_array(pts),
                           debug ? vnf_from_polygons(faces) 
                                 : vnf_triangulate(vnf_from_polygons(faces))
@@ -2114,7 +2114,7 @@ function _circle_mask(r) =
 //           ]),
 //           radius = [0,0,each repeat(slotradius,4),0,0], closed=false
 //       )
-//   ) apply(left(max(columns(slot,0))/2)*fwd(min(columns(slot,1))), slot);
+//   ) apply(left(max(column(slot,0))/2)*fwd(min(column(slot,1))), slot);
 //   stroke(slot(15,29,7));
 // Example: A cylindrical container with rounded edges and a rounded finger slot.
 //   function slot(slotwidth, slotheight, slotradius) = let(
@@ -2138,7 +2138,7 @@ function _circle_mask(r) =
 //           ]),
 //           radius = [0,0,each repeat(slotradius,4),0,0], closed=false
 //       )
-//   ) apply(left(max(columns(slot,0))/2)*fwd(min(columns(slot,1))), slot);
+//   ) apply(left(max(column(slot,0))/2)*fwd(min(column(slot,1))), slot);
 //   diam = 80;
 //   wall = 4;
 //   height = 40;
@@ -2162,12 +2162,12 @@ module bent_cutout_mask(r, thickness, path, radius, convexity=10)
   path = clockwise_polygon(path);
   curvepoints = arc(d=thickness, angle = [-180,0]);
   profiles = [for(pt=curvepoints) _cyl_hole(r+pt.x,apply(xscale((r+pt.x)/r), offset(path,delta=thickness/2+pt.y,check_valid=false,closed=true)))];
-  pathx = columns(path,0);
+  pathx = column(path,0);
   minangle = (min(pathx)-thickness/2)*360/(2*PI*r);
   maxangle = (max(pathx)+thickness/2)*360/(2*PI*r);
   mindist = (r+thickness/2)/cos((maxangle-minangle)/2);
   assert(maxangle-minangle<180,"Cutout angle span is too large.  Must be smaller than 180.");
-  zmean = mean(columns(path,1));
+  zmean = mean(column(path,1));
   innerzero = repeat([0,0,zmean], len(path));
   outerpt = repeat( [1.5*mindist*cos((maxangle+minangle)/2),1.5*mindist*sin((maxangle+minangle)/2),zmean], len(path));
   vnf_polyhedron(vnf_vertex_array([innerzero, each profiles, outerpt],col_wrap=true),convexity=convexity);
diff --git a/shapes2d.scad b/shapes2d.scad
index f1565d4..6936ba3 100644
--- a/shapes2d.scad
+++ b/shapes2d.scad
@@ -363,7 +363,7 @@ function regular_ngon(n=6, r, d, or, od, ir, id, side, rounding=0, realign=false
                     )
                     each arc(N=steps, cp=p, r=rounding, start=a+180/n, angle=-360/n)
                 ],
-                maxx_idx = max_index(columns(path2,0)),
+                maxx_idx = max_index(column(path2,0)),
                 path3 = polygon_shift(path2,maxx_idx)
             ) path3
         ),
@@ -1005,7 +1005,7 @@ function teardrop2d(r, ang=45, cap_h, d, anchor=CENTER, spin=0) =
                 [-cap_w,cap_h]
             ], closed=true
         ),
-        maxx_idx = max_index(columns(path,0)),
+        maxx_idx = max_index(column(path,0)),
         path2 = polygon_shift(path,maxx_idx)
     ) reorient(anchor,spin, two_d=true, path=path2, p=path2);
 
@@ -1060,7 +1060,7 @@ function glued_circles(r, spread=10, tangent=30, d, anchor=CENTER, spin=0) =
             [for (i=[0:1:lobesegs]) let(a=sa1+i*lobestep+180) r  * [cos(a),sin(a)] + cp1],
             tangent==0? [] : [for (i=[0:1:arcsegs])  let(a=ea2-i*arcstep)      r2 * [cos(a),sin(a)] + cp2]
         ),
-        maxx_idx = max_index(columns(path,0)),
+        maxx_idx = max_index(column(path,0)),
         path2 = reverse_polygon(polygon_shift(path,maxx_idx))
     ) reorient(anchor,spin, two_d=true, path=path2, extent=true, p=path2);
 
diff --git a/shapes3d.scad b/shapes3d.scad
index dcdd8cc..8f63b38 100644
--- a/shapes3d.scad
+++ b/shapes3d.scad
@@ -2266,7 +2266,7 @@ module path_text(path, text, font, size, thickness, lettersize, offset=0, revers
   usernorm = is_def(normal);
   usetop = is_def(top);
   
-  normpts = is_undef(normal) ? (reverse?1:-1)*columns(pts,3) : _cut_interp(pts,path, normal);
+  normpts = is_undef(normal) ? (reverse?1:-1)*column(pts,3) : _cut_interp(pts,path, normal);
   toppts = is_undef(top) ? undef : _cut_interp(pts,path,top);
   for(i=idx(text))
       let( tangent = pts[i][2] )
diff --git a/skin.scad b/skin.scad
index e019846..37f01cc 100644
--- a/skin.scad
+++ b/skin.scad
@@ -985,7 +985,7 @@ function path_sweep2d(shape, path, closed=false, caps, quality=1, style="min_edg
                      [for(pt = profile)
                         let( 
                             ofs = offset(path, delta=-flip*pt.x, return_faces=true,closed=closed, quality=quality),
-                            map = columns(_ofs_vmap(ofs,closed=closed),1)
+                            map = column(_ofs_vmap(ofs,closed=closed),1)
                         ) 
                         select(path3d(ofs[0],pt.y),map)
                       ]
diff --git a/std.scad b/std.scad
index 8c9066d..bde5526 100644
--- a/std.scad
+++ b/std.scad
@@ -22,6 +22,7 @@ include <paths.scad>
 include <edges.scad>
 include <arrays.scad>
 include <math.scad>
+include <linalg.scad>
 include <trigonometry.scad>
 include <vectors.scad>
 include <quaternions.scad>
diff --git a/tests/test_arrays.scad b/tests/test_arrays.scad
index a269f07..20cdc73 100644
--- a/tests/test_arrays.scad
+++ b/tests/test_arrays.scad
@@ -451,36 +451,6 @@ module test_add_scalar() {
 test_add_scalar();
 
 
-module test_columns() {
-    v = [[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]];
-    assert(columns(v,2) == [3, 7, 11, 15]);
-    assert(columns(v,[2]) == [[3], [7], [11], [15]]);
-    assert(columns(v,[2,1]) == [[3, 2], [7, 6], [11, 10], [15, 14]]);
-    assert(columns(v,[1:3]) == [[2, 3, 4], [6, 7, 8], [10, 11, 12], [14, 15, 16]]);
-}
-test_columns();
-
-
-// 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]);
@@ -535,67 +505,8 @@ module test_zip() {
     v3 = [8,9,10,11];
     assert(zip(v1,v3) == [[1,8],[2,9],[3,10],[4,11]]);
     assert(zip([v1,v3]) == [[1,8],[2,9],[3,10],[4,11]]);
-    assert(zip([v1,v2],fit="short") == [[1,5],[2,6],[3,7]]);
-    assert(zip([v1,v2],fit="long") == [[1,5],[2,6],[3,7],[4,undef]]);
-    assert(zip([v1,v2],fit="long", fill=0) == [[1,5],[2,6],[3,7],[4,0]]);
-    assert(zip([v1,v2,v3],fit="long") == [[1,5,8],[2,6,9],[3,7,10],[4,undef,11]]);
 }
-//test_zip();
-
-module test_hstack() {
-    M = ident(3);
-    v1 = [2,3,4];
-    v2 = [5,6,7];
-    v3 = [8,9,10];
-    a = hstack(v1,v2);   
-    b = hstack(v1,v2,v3);
-    c = hstack([M,v1,M]);
-    d = hstack(columns(M,0), columns(M,[1, 2]));
-    assert_equal(a,[[2, 5], [3, 6], [4, 7]]);
-    assert_equal(b,[[2, 5, 8], [3, 6, 9], [4, 7, 10]]);
-    assert_equal(c,[[1, 0, 0, 2, 1, 0, 0], [0, 1, 0, 3, 0, 1, 0], [0, 0, 1, 4, 0, 0, 1]]);
-    assert_equal(d,M);
-    strmat = [["three","four"], ["five","six"]];
-    assert_equal(hstack(strmat,strmat), [["three", "four", "three", "four"], ["five", "six", "five", "six"]]);
-    strvec = ["one","two"];
-    assert_equal(hstack(strvec,strmat),[["o", "n", "e", "three", "four"], ["t", "w", "o", "five", "six"]]);
-}
-test_hstack();
-
-
-module test_block_matrix() {
-    A = [[1,2],[3,4]];
-    B = ident(2);
-    assert_equal(block_matrix([[A,B],[B,A],[A,B]]), [[1,2,1,0],[3,4,0,1],[1,0,1,2],[0,1,3,4],[1,2,1,0],[3,4,0,1]]);
-    assert_equal(block_matrix([[A,B],ident(4)]), [[1,2,1,0],[3,4,0,1],[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]]);
-    text = [["aa","bb"],["cc","dd"]];
-    assert_equal(block_matrix([[text,B]]), [["aa","bb",1,0],["cc","dd",0,1]]);
-}
-test_block_matrix();
-
-
-module test_diagonal_matrix() {
-    assert_equal(diagonal_matrix([1,2,3]), [[1,0,0],[0,2,0],[0,0,3]]);
-    assert_equal(diagonal_matrix([1,"c",2]), [[1,0,0],[0,"c",0],[0,0,2]]);
-    assert_equal(diagonal_matrix([1,"c",2],"X"), [[1,"X","X"],["X","c","X"],["X","X",2]]);
-    assert_equal(diagonal_matrix([[1,1],[2,2],[3,3]], [0,0]), [[ [1,1],[0,0],[0,0]], [[0,0],[2,2],[0,0]], [[0,0],[0,0],[3,3]]]);
-}
-test_diagonal_matrix();
-
-module test_submatrix_set() {
-    test = [[1,2,3,4,5],[6,7,8,9,10],[11,12,13,14,15], [16,17,18,19,20]];
-    ragged = [[1,2,3,4,5],[6,7,8,9,10],[11,12], [16,17]];
-    assert_equal(submatrix_set(test,[[9,8],[7,6]]), [[9,8,3,4,5],[7,6,8,9,10],[11,12,13,14,15], [16,17,18,19,20]]);
-    assert_equal(submatrix_set(test,[[9,7],[8,6]],1),[[1,2,3,4,5],[9,7,8,9,10],[8,6,13,14,15], [16,17,18,19,20]]);
-    assert_equal(submatrix_set(test,[[9,8],[7,6]],n=1), [[1,9,8,4,5],[6,7,6,9,10],[11,12,13,14,15], [16,17,18,19,20]]);
-    assert_equal(submatrix_set(test,[[9,8],[7,6]],1,2), [[1,2,3,4,5],[6,7,9,8,10],[11,12,7,6,15], [16,17,18,19,20]]);
-    assert_equal(submatrix_set(test,[[9,8],[7,6]],-1,-1), [[6,2,3,4,5],[6,7,8,9,10],[11,12,13,14,15], [16,17,18,19,20]]);
-    assert_equal(submatrix_set(test,[[9,8],[7,6]],n=4), [[1,2,3,4,9],[6,7,8,9,7],[11,12,13,14,15], [16,17,18,19,20]]);
-    assert_equal(submatrix_set(test,[[9,8],[7,6]],7,7), [[1,2,3,4,5],[6,7,8,9,10],[11,12,13,14,15], [16,17,18,19,20]]);
-    assert_equal(submatrix_set(ragged, [["a","b"],["c","d"]], 1, 1), [[1,2,3,4,5],[6,"a","b",9,10],[11,"c"], [16,17]]);
-    assert_equal(submatrix_set(test, [[]]), test);
-}
-test_submatrix_set();
+test_zip();
 
 
 module test_array_group() {
@@ -645,13 +556,4 @@ module test_array_dim() {
 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();
-
-
 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
diff --git a/tests/test_math.scad b/tests/test_math.scad
index 2cbabf0..c5820f5 100644
--- a/tests/test_math.scad
+++ b/tests/test_math.scad
@@ -595,13 +595,13 @@ 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]);
+    assert_equal(median([2,3,7]), 3);
+    assert_equal(median([2,4,5,8]), 4.5);
 }
 test_median();
-*/
+
 
 
 module test_convolve() {
@@ -618,40 +618,6 @@ 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]]);
-}
-test_matrix_inverse();
-
-
-module test_det2() {
-    assert_equal(det2([[6,-2], [1,8]]), 50);
-    assert_equal(det2([[4,7], [3,2]]), -13);
-    assert_equal(det2([[4,3], [3,4]]), 7);
-}
-test_det2();
-
-
-module test_det3() {
-    M = [ [6,4,-2], [1,-2,8], [1,5,7] ];
-    assert_equal(det3(M), -334);
-}
-test_det3();
-
-
-module test_determinant() {
-    M = [ [6,4,-2,9], [1,-2,8,3], [1,5,7,6], [4,2,5,1] ];
-    assert_equal(determinant(M), 2267);
-}
-test_determinant();
-
-
-module test_matrix_trace() {
-    M = [ [6,4,-2,9], [1,-2,8,3], [1,5,7,6], [4,2,5,1] ];
-    assert_equal(matrix_trace(M), 6-2+7+1);
-}
-test_matrix_trace();
-
 
 // Logic
 
@@ -975,40 +941,6 @@ test_back_substitute();
 
 
 
-module test_norm_fro(){
-  assert_approx(norm_fro([[2,3,4],[4,5,6]]), 10.29563014098700);
-
-} test_norm_fro();  
-
-
-module test_linear_solve(){
-  M = [[-2,-5,-1,3],
-       [3,7,6,2],
-       [6,5,-1,-6],
-       [-7,1,2,3]];
-  assert_approx(linear_solve(M, [-3,43,-11,13]), [1,2,3,4]);
-  assert_approx(linear_solve(M, [[-5,8],[18,-61],[4,7],[-1,-12]]), [[1,-2],[1,-3],[1,-4],[1,-5]]);
-  assert_approx(linear_solve([[2]],[4]), [2]);
-  assert_approx(linear_solve([[2]],[[4,8]]), [[2, 4]]);
-  assert_approx(linear_solve(select(M,0,2), [2,4,4]), [   2.254871220604705e+00,
-                                                         -8.378819388897780e-01,
-                                                          2.330507118860985e-01,
-                                                          8.511278195488737e-01]);
-  assert_approx(linear_solve(columns(M,[0:2]), [2,4,4,4]),
-                 [-2.457142857142859e-01,
-                   5.200000000000000e-01,
-                   7.428571428571396e-02]);
-  assert_approx(linear_solve([[1,2,3,4]], [2]), [0.066666666666666, 0.13333333333, 0.2, 0.266666666666]);
-  assert_approx(linear_solve([[1],[2],[3],[4]], [4,3,2,1]), [2/3]);
-  rd = [[-2,-5,-1,3],
-        [3,7,6,2],
-        [3,7,6,2],
-        [-7,1,2,3]];
-  assert_equal(linear_solve(rd,[1,2,3,4]),[]);
-  assert_equal(linear_solve(select(rd,0,2), [2,4,4]), []);
-  assert_equal(linear_solve(transpose(select(rd,0,2)), [2,4,3,4]), []);
-}
-test_linear_solve();
 
 
 module test_cumprod(){
@@ -1186,85 +1118,6 @@ module test_quadratic_roots(){
 test_quadratic_roots();
 
 
-module test_null_space(){
-    assert_equal(null_space([[3,2,1],[3,6,3],[3,9,-3]]),[]);
-
-    function nullcheck(A,dim) =
-      let(v=null_space(A))
-        len(v)==dim && all_zero(A*transpose(v),eps=1e-12);
-    
-   A = [[-1, 2, -5, 2],[-3,-1,3,-3],[5,0,5,0],[3,-4,11,-4]];
-   assert(nullcheck(A,1));
-
-   B = [
-        [  4,    1,    8,    6,   -2,    3],
-        [ 10,    5,   10,   10,    0,    5],
-        [  8,    1,    8,    8,   -6,    1],
-        [ -8,   -8,    6,   -1,   -8,   -1],
-        [  2,    2,    0,    1,    2,    1],
-        [  2,   -3,   10,    6,   -8,    1],
-       ];
-   assert(nullcheck(B,3));
-}
-test_null_space();
-
-module test_qr_factor() {
-  // Check that R is upper triangular
-  function is_ut(R) =
-     let(bad = [for(i=[1:1:len(R)-1], j=[0:min(i-1, len(R[0])-1)]) if (!approx(R[i][j],0)) 1])
-     bad == [];
-
-  // Test the R is upper trianglar, Q is orthogonal and qr=M
-  function qrok(qr,M) =
-     is_ut(qr[1]) && approx(qr[0]*transpose(qr[0]), ident(len(qr[0]))) && approx(qr[0]*qr[1],M) && qr[2]==ident(len(qr[2]));
-
-  // Test the R is upper trianglar, Q is orthogonal, R diagonal non-increasing and qrp=M
-  function qrokpiv(qr,M) =
-       is_ut(qr[1])
-    && approx(qr[0]*transpose(qr[0]), ident(len(qr[0])))
-    && approx(qr[0]*qr[1]*transpose(qr[2]),M)
-    && is_decreasing([for(i=[0:1:min(len(qr[1]),len(qr[1][0]))-1]) abs(qr[1][i][i])]);
-
-  
-  M = [[1,2,9,4,5],
-       [6,7,8,19,10],
-       [11,12,13,14,15],
-       [1,17,18,19,20],
-       [21,22,10,24,25]];
-  
-  assert(qrok(qr_factor(M),M));
-  assert(qrok(qr_factor(select(M,0,3)),select(M,0,3)));
-  assert(qrok(qr_factor(transpose(select(M,0,3))),transpose(select(M,0,3))));
-
-  A = [[1,2,9,4,5],
-       [6,7,8,19,10],
-       [0,0,0,0,0],
-       [1,17,18,19,20],
-       [21,22,10,24,25]];
-  assert(qrok(qr_factor(A),A));
-
-  B = [[1,2,0,4,5],
-       [6,7,0,19,10],
-       [0,0,0,0,0],
-       [1,17,0,19,20],
-       [21,22,0,24,25]];
-
-  assert(qrok(qr_factor(B),B));
-  assert(qrok(qr_factor([[7]]), [[7]]));
-  assert(qrok(qr_factor([[1,2,3]]), [[1,2,3]]));
-  assert(qrok(qr_factor([[1],[2],[3]]), [[1],[2],[3]]));
-
-
-  assert(qrokpiv(qr_factor(M,pivot=true),M));
-  assert(qrokpiv(qr_factor(select(M,0,3),pivot=true),select(M,0,3)));
-  assert(qrokpiv(qr_factor(transpose(select(M,0,3)),pivot=true),transpose(select(M,0,3))));
-  assert(qrokpiv(qr_factor(B,pivot=true),B));
-  assert(qrokpiv(qr_factor([[7]],pivot=true), [[7]]));
-  assert(qrokpiv(qr_factor([[1,2,3]],pivot=true), [[1,2,3]]));
-  assert(qrokpiv(qr_factor([[1],[2],[3]],pivot=true), [[1],[2],[3]]));
-}
-test_qr_factor();
-
 
 module test_poly_mult(){
   assert_equal(poly_mult([3,2,1],[4,5,6,7]),[12,23,32,38,20,7]);
diff --git a/vectors.scad b/vectors.scad
index c57ebdd..2f2d9e2 100644
--- a/vectors.scad
+++ b/vectors.scad
@@ -543,8 +543,8 @@ function vector_nearest(query, k, target) =
             "More results are requested than the number of points.")
     tgpts
     ?   let( tree = _bt_tree(target, count(len(target))) )
-        columns(_bt_nearest( query, k, target,  tree),0)
-    :   columns(_bt_nearest( query, k, target[0], target[1]),0);
+        column(_bt_nearest( query, k, target,  tree),0)
+    :   column(_bt_nearest( query, k, target[0], target[1]),0);
 
 
 //Ball tree nearest