From 4e2092f347007681fb105da828ba3309875d399a Mon Sep 17 00:00:00 2001
From: Adrian Mariano <avm4@cornell.edu>
Date: Sat, 11 Jul 2020 12:22:44 -0400
Subject: [PATCH 1/3] Added a bunch of tests

---
 tests/test_math.scad | 952 +++++++++++++++++++++++--------------------
 1 file changed, 508 insertions(+), 444 deletions(-)

diff --git a/tests/test_math.scad b/tests/test_math.scad
index 16f4988..dcbad7d 100644
--- a/tests/test_math.scad
+++ b/tests/test_math.scad
@@ -1,586 +1,650 @@
-include <BOSL2/std.scad>
+include <BOSL/constants.scad>
+include <BOSL/math.scad>
+
+eps = 1e-9;
 
 // Simple Calculations
 
 module test_quant() {
-    assert_equal(quant(-4,3), -3);
-    assert_equal(quant(-3,3), -3);
-    assert_equal(quant(-2,3), -3);
-    assert_equal(quant(-1,3), 0);
-    assert_equal(quant(0,3), 0);
-    assert_equal(quant(1,3), 0);
-    assert_equal(quant(2,3), 3);
-    assert_equal(quant(3,3), 3);
-    assert_equal(quant(4,3), 3);
-    assert_equal(quant(7,3), 6);
-    assert_equal(quant([12,13,13.1,14,14.1,15,16],4), [12,12,12,16,16,16,16]);
-    assert_equal(quant([9,10,10.4,10.5,11,12],3), [9,9,9,12,12,12]);
-    assert_equal(quant([[9,10,10.4],[10.5,11,12]],3), [[9,9,9],[12,12,12]]);
+	assert(quant(-4,3) == -3);
+	assert(quant(-3,3) == -3);
+	assert(quant(-2,3) == -3);
+	assert(quant(-1,3) == 0);
+	assert(quant(0,3) == 0);
+	assert(quant(1,3) == 0);
+	assert(quant(2,3) == 3);
+	assert(quant(3,3) == 3);
+	assert(quant(4,3) == 3);
+	assert(quant(7,3) == 6);
 }
 test_quant();
 
 
 module test_quantdn() {
-    assert_equal(quantdn(-4,3), -6);
-    assert_equal(quantdn(-3,3), -3);
-    assert_equal(quantdn(-2,3), -3);
-    assert_equal(quantdn(-1,3), -3);
-    assert_equal(quantdn(0,3), 0);
-    assert_equal(quantdn(1,3), 0);
-    assert_equal(quantdn(2,3), 0);
-    assert_equal(quantdn(3,3), 3);
-    assert_equal(quantdn(4,3), 3);
-    assert_equal(quantdn(7,3), 6);
-    assert_equal(quantdn([12,13,13.1,14,14.1,15,16],4), [12,12,12,12,12,12,16]);
-    assert_equal(quantdn([9,10,10.4,10.5,11,12],3), [9,9,9,9,9,12]);
-    assert_equal(quantdn([[9,10,10.4],[10.5,11,12]],3), [[9,9,9],[9,9,12]]);
+	assert(quantdn(-4,3) == -6);
+	assert(quantdn(-3,3) == -3);
+	assert(quantdn(-2,3) == -3);
+	assert(quantdn(-1,3) == -3);
+	assert(quantdn(0,3) == 0);
+	assert(quantdn(1,3) == 0);
+	assert(quantdn(2,3) == 0);
+	assert(quantdn(3,3) == 3);
+	assert(quantdn(4,3) == 3);
+	assert(quantdn(7,3) == 6);
 }
 test_quantdn();
 
 
 module test_quantup() {
-    assert_equal(quantup(-4,3), -3);
-    assert_equal(quantup(-3,3), -3);
-    assert_equal(quantup(-2,3), 0);
-    assert_equal(quantup(-1,3), 0);
-    assert_equal(quantup(0,3), 0);
-    assert_equal(quantup(1,3), 3);
-    assert_equal(quantup(2,3), 3);
-    assert_equal(quantup(3,3), 3);
-    assert_equal(quantup(4,3), 6);
-    assert_equal(quantup(7,3), 9);
-    assert_equal(quantup([12,13,13.1,14,14.1,15,16],4), [12,16,16,16,16,16,16]);
-    assert_equal(quantup([9,10,10.4,10.5,11,12],3), [9,12,12,12,12,12]);
-    assert_equal(quantup([[9,10,10.4],[10.5,11,12]],3), [[9,12,12],[12,12,12]]);
+	assert(quantup(-4,3) == -3);
+	assert(quantup(-3,3) == -3);
+	assert(quantup(-2,3) == 0);
+	assert(quantup(-1,3) == 0);
+	assert(quantup(0,3) == 0);
+	assert(quantup(1,3) == 3);
+	assert(quantup(2,3) == 3);
+	assert(quantup(3,3) == 3);
+	assert(quantup(4,3) == 6);
+	assert(quantup(7,3) == 9);
 }
 test_quantup();
 
 
 module test_constrain() {
-    assert_equal(constrain(-2,-1,1), -1);
-    assert_equal(constrain(-1.75,-1,1), -1);
-    assert_equal(constrain(-1,-1,1), -1);
-    assert_equal(constrain(-0.75,-1,1), -0.75);
-    assert_equal(constrain(0,-1,1), 0);
-    assert_equal(constrain(0.75,-1,1), 0.75);
-    assert_equal(constrain(1,-1,1), 1);
-    assert_equal(constrain(1.75,-1,1), 1);
-    assert_equal(constrain(2,-1,1), 1);
+	assert(constrain(-2,-1,1) == -1);
+	assert(constrain(-1.75,-1,1) == -1);
+	assert(constrain(-1,-1,1) == -1);
+	assert(constrain(-0.75,-1,1) == -0.75);
+	assert(constrain(0,-1,1) == 0);
+	assert(constrain(0.75,-1,1) == 0.75);
+	assert(constrain(1,-1,1) == 1);
+	assert(constrain(1.75,-1,1) == 1);
+	assert(constrain(2,-1,1) == 1);
 }
 test_constrain();
 
 
-module test_is_matrix() {
-    assert(is_matrix([[2,3,4],[5,6,7],[8,9,10]]));
-    assert(is_matrix([[2,3,4],[5,6,7],[8,9,10]],square=true));
-    assert(is_matrix([[2,3,4],[5,6,7],[8,9,10]],square=false));
-    assert(is_matrix([[2,3],[5,6],[8,9]],m=3,n=2));
-    assert(is_matrix([[2,3,4],[5,6,7]],m=2,n=3));
-    assert(!is_matrix([[2,3,4],[5,6,7]],m=2,n=3,square=true));
-    assert(is_matrix([[2,3,4],[5,6,7],[8,9,10]],square=false));
-    assert(!is_matrix([[2,3],[5,6],[8,9]],m=2,n=3));
-    assert(!is_matrix([[2,3,4],[5,6,7]],m=3,n=2));
-    assert(!is_matrix(undef));
-    assert(!is_matrix(NAN));
-    assert(!is_matrix(INF));
-    assert(!is_matrix(-5));
-    assert(!is_matrix(0));
-    assert(!is_matrix(5));
-    assert(!is_matrix(""));
-    assert(!is_matrix("foo"));
-    assert(!is_matrix([3,4,5]));
-    assert(!is_matrix([]));
-}
-test_is_matrix();
-
-
-module test_approx() {
-    assert_equal(approx(PI, 3.141592653589793236), true);
-    assert_equal(approx(PI, 3.1415926), false);
-    assert_equal(approx(PI, 3.1415926, eps=1e-6), true);
-    assert_equal(approx(-PI, -3.141592653589793236), true);
-    assert_equal(approx(-PI, -3.1415926), false);
-    assert_equal(approx(-PI, -3.1415926, eps=1e-6), true);
-    assert_equal(approx(1/3, 0.3333333333), true);
-    assert_equal(approx(-1/3, -0.3333333333), true);
-    assert_equal(approx(10*[cos(30),sin(30)], 10*[sqrt(3)/2, 1/2]), true);
-}
-test_approx();
-
-
-module test_min_index() {
-    vals = rands(-100,100,100);
-    minval = min(vals);
-    minidx = min_index(vals);
-    assert_equal(vals[minidx], minval);
-    assert_equal(min_index([3,4,5,6]), 0);
-    assert_equal(min_index([4,3,5,6]), 1);
-    assert_equal(min_index([4,5,3,6]), 2);
-    assert_equal(min_index([4,5,6,3]), 3);
-    assert_equal(min_index([6,5,4,3]), 3);
-    assert_equal(min_index([6,3,4,5]), 1);
-    assert_equal(min_index([-56,72,-874,5]), 2);
-}
-test_min_index();
-
-
-module test_max_index() {
-    vals = rands(-100,100,100);
-    maxval = max(vals);
-    maxidx = max_index(vals);
-    assert_equal(vals[maxidx], maxval);
-    assert_equal(max_index([3,4,5,6]), 3);
-    assert_equal(max_index([3,4,6,5]), 2);
-    assert_equal(max_index([3,6,4,5]), 1);
-    assert_equal(max_index([6,3,4,5]), 0);
-    assert_equal(max_index([5,6,4,3]), 1);
-    assert_equal(max_index([-56,72,-874,5]), 1);
-}
-test_max_index();
-
-
 module test_posmod() {
-    assert_equal(posmod(-5,3), 1);
-    assert_equal(posmod(-4,3), 2);
-    assert_equal(posmod(-3,3), 0);
-    assert_equal(posmod(-2,3), 1);
-    assert_equal(posmod(-1,3), 2);
-    assert_equal(posmod(0,3), 0);
-    assert_equal(posmod(1,3), 1);
-    assert_equal(posmod(2,3), 2);
-    assert_equal(posmod(3,3), 0);
+	assert(posmod(-5,3) == 1);
+	assert(posmod(-4,3) == 2);
+	assert(posmod(-3,3) == 0);
+	assert(posmod(-2,3) == 1);
+	assert(posmod(-1,3) == 2);
+	assert(posmod(0,3) == 0);
+	assert(posmod(1,3) == 1);
+	assert(posmod(2,3) == 2);
+	assert(posmod(3,3) == 0);
 }
 test_posmod();
 
 
-module test_modang() {
-    assert_equal(modang(-700), 20);
-    assert_equal(modang(-270), 90);
-    assert_equal(modang(-120), -120);
-    assert_equal(modang(120), 120);
-    assert_equal(modang(270), -90);
-    assert_equal(modang(700), -20);
-}
-test_modang();
-
-
 module test_modrange() {
-    assert_equal(modrange(-5,5,3), [1,2]);
-    assert_equal(modrange(-1,4,3), [2,0,1]);
-    assert_equal(modrange(1,8,10,step=2), [1,3,5,7]);
-    assert_equal(modrange(5,12,10,step=2), [5,7,9,1]);
+	assert(modrange(-5,5,3) == [1,2]);
+	assert(modrange(-1,4,3) == [2,0,1]);
+	assert(modrange(1,8,10,step=2) == [1,3,5,7]);
+	assert(modrange(5,12,10,step=2) == [5,7,9,1]);
 }
 test_modrange();
 
 
-module test_sqr() {
-    assert_equal(sqr(-3), 9);
-    assert_equal(sqr(0), 0);
-    assert_equal(sqr(1), 1);
-    assert_equal(sqr(2), 4);
-    assert_equal(sqr(2.5), 6.25);
-    assert_equal(sqr(3), 9);
-    assert_equal(sqr(16), 256);
-}
-test_sqr();
-
-
-module test_log2() {
-    assert_equal(log2(0.125), -3);
-    assert_equal(log2(16), 4);
-    assert_equal(log2(256), 8);
-}
-test_log2();
-
-
-module test_rand_int() {
-    nums = rand_int(-100,100,1000,seed=2134);
-    assert_equal(len(nums), 1000);
-    for (num = nums) {
-        assert(num>=-100);
-        assert(num<=100);
-        assert_equal(num, floor(num));
-    }
-}
-test_rand_int();
-
-
-module test_gaussian_rands() {
-    nums1 = gaussian_rands(0,10,1000,seed=2132);
-    nums2 = gaussian_rands(0,10,1000,seed=2130);
-    nums3 = gaussian_rands(0,10,1000,seed=2132);
-    assert_equal(len(nums1), 1000);
-    assert_equal(len(nums2), 1000);
-    assert_equal(len(nums3), 1000);
-    assert_equal(nums1, nums3);
-    assert(nums1!=nums2);
-}
-test_gaussian_rands();
-
-
-module test_log_rands() {
-    nums1 = log_rands(0,100,10,1000,seed=2189);
-    nums2 = log_rands(0,100,10,1000,seed=2310);
-    nums3 = log_rands(0,100,10,1000,seed=2189);
-    assert_equal(len(nums1), 1000);
-    assert_equal(len(nums2), 1000);
-    assert_equal(len(nums3), 1000);
-    assert_equal(nums1, nums3);
-    assert(nums1!=nums2);
-}
-test_log_rands();
-
-
 module test_segs() {
-    assert_equal(segs(50,$fn=8), 8);
-    assert_equal(segs(50,$fa=2,$fs=2), 158);
+	assert(segs(50,$fn=8) == 8);
+	assert(segs(50,$fa=2,$fs=2) == 158);
 }
 test_segs();
 
 
 module test_lerp() {
-    assert_equal(lerp(-20,20,0), -20);
-    assert_equal(lerp(-20,20,0.25), -10);
-    assert_equal(lerp(-20,20,0.5), 0);
-    assert_equal(lerp(-20,20,0.75), 10);
-    assert_equal(lerp(-20,20,1), 20);
-    assert_equal(lerp(-20,20,[0,0.25,0.5,0.75,1]), [-20,-10,0,10,20]);
-    assert_equal(lerp(-20,20,[0:0.25:1]), [-20,-10,0,10,20]);
-    assert_equal(lerp([10,10],[30,-10],0.5), [20,0]);
+	assert(lerp(-20,20,0) == -20);
+	assert(lerp(-20,20,0.25) == -10);
+	assert(lerp(-20,20,0.5) == 0);
+	assert(lerp(-20,20,0.75) == 10);
+	assert(lerp(-20,20,1) == 20);
+	assert(lerp([10,10],[30,-10],0.5) == [20,0]);
 }
 test_lerp();
 
 
 module test_hypot() {
-    assert_approx(hypot(20,30), norm([20,30]));
+	assert(hypot(20,30) == norm([20,30]));
 }
 test_hypot();
 
 
 module test_sinh() {
-    assert_approx(sinh(-2), -3.6268604078);
-    assert_approx(sinh(-1), -1.1752011936);
-    assert_approx(sinh(0), 0);
-    assert_approx(sinh(1), 1.1752011936);
-    assert_approx(sinh(2), 3.6268604078);
+	assert(abs(sinh(-2)+3.6268604078) < eps);
+	assert(abs(sinh(-1)+1.1752011936) < eps);
+	assert(abs(sinh(0)) < eps);
+	assert(abs(sinh(1)-1.1752011936) < eps);
+	assert(abs(sinh(2)-3.6268604078) < eps);
 }
 test_sinh();
 
 
 module test_cosh() {
-    assert_approx(cosh(-2), 3.7621956911);
-    assert_approx(cosh(-1), 1.5430806348);
-    assert_approx(cosh(0), 1);
-    assert_approx(cosh(1), 1.5430806348);
-    assert_approx(cosh(2), 3.7621956911);
+	assert(abs(cosh(-2)-3.7621956911) < eps);
+	assert(abs(cosh(-1)-1.5430806348) < eps);
+	assert(abs(cosh(0)-1) < eps);
+	assert(abs(cosh(1)-1.5430806348) < eps);
+	assert(abs(cosh(2)-3.7621956911) < eps);
 }
 test_cosh();
 
 
 module test_tanh() {
-    assert_approx(tanh(-2), -0.9640275801);
-    assert_approx(tanh(-1), -0.761594156);
-    assert_approx(tanh(0), 0);
-    assert_approx(tanh(1), 0.761594156);
-    assert_approx(tanh(2), 0.9640275801);
+	assert(abs(tanh(-2)+0.9640275801) < eps);
+	assert(abs(tanh(-1)+0.761594156) < eps);
+	assert(abs(tanh(0)) < eps);
+	assert(abs(tanh(1)-0.761594156) < eps);
+	assert(abs(tanh(2)-0.9640275801) < eps);
 }
 test_tanh();
 
 
 module test_asinh() {
-    assert_approx(asinh(sinh(-2)), -2);
-    assert_approx(asinh(sinh(-1)), -1);
-    assert_approx(asinh(sinh(0)), 0);
-    assert_approx(asinh(sinh(1)), 1);
-    assert_approx(asinh(sinh(2)), 2);
+	assert(abs(asinh(sinh(-2))+2) < eps);
+	assert(abs(asinh(sinh(-1))+1) < eps);
+	assert(abs(asinh(sinh(0))) < eps);
+	assert(abs(asinh(sinh(1))-1) < eps);
+	assert(abs(asinh(sinh(2))-2) < eps);
 }
 test_asinh();
 
 
 module test_acosh() {
-    assert_approx(acosh(cosh(-2)), 2);
-    assert_approx(acosh(cosh(-1)), 1);
-    assert_approx(acosh(cosh(0)), 0);
-    assert_approx(acosh(cosh(1)), 1);
-    assert_approx(acosh(cosh(2)), 2);
+	assert(abs(acosh(cosh(-2))-2) < eps);
+	assert(abs(acosh(cosh(-1))-1) < eps);
+	assert(abs(acosh(cosh(0))) < eps);
+	assert(abs(acosh(cosh(1))-1) < eps);
+	assert(abs(acosh(cosh(2))-2) < eps);
 }
 test_acosh();
 
 
 module test_atanh() {
-    assert_approx(atanh(tanh(-2)), -2);
-    assert_approx(atanh(tanh(-1)), -1);
-    assert_approx(atanh(tanh(0)), 0);
-    assert_approx(atanh(tanh(1)), 1);
-    assert_approx(atanh(tanh(2)), 2);
+	assert(abs(atanh(tanh(-2))+2) < eps);
+	assert(abs(atanh(tanh(-1))+1) < eps);
+	assert(abs(atanh(tanh(0))) < eps);
+	assert(abs(atanh(tanh(1))-1) < eps);
+	assert(abs(atanh(tanh(2))-2) < eps);
 }
 test_atanh();
 
 
 module test_sum() {
-    assert_equal(sum([]), 0);
-    assert_equal(sum([],dflt=undef), undef);
-    assert_equal(sum([1,2,3]), 6);
-    assert_equal(sum([-2,-1,0,1,2]), 0);
-    assert_equal(sum([[1,2,3], [3,4,5], [5,6,7]]), [9,12,15]);
+	assert(sum([1,2,3]) == 6);
+	assert(sum([-2,-1,0,1,2]) == 0);
+	assert(sum([[1,2,3], [3,4,5], [5,6,7]]) == [9,12,15]);
 }
 test_sum();
 
 
-module test_cumsum() {
-    assert_equal(cumsum([]), []);
-    assert_equal(cumsum([1,1,1]), [1,2,3]);
-    assert_equal(cumsum([2,2,2]), [2,4,6]);
-    assert_equal(cumsum([1,2,3]), [1,3,6]);
-    assert_equal(cumsum([-2,-1,0,1,2]), [-2,-3,-3,-2,0]);
-    assert_equal(cumsum([[1,2,3], [3,4,5], [5,6,7]]), [[1,2,3],[4,6,8],[9,12,15]]);
-}
-test_cumsum();
-
-
 module test_sum_of_squares() {
-    assert_equal(sum_of_squares([1,2,3]), 14);
-    assert_equal(sum_of_squares([1,2,4]), 21);
-    assert_equal(sum_of_squares([-3,-2,-1]), 14);
+	assert(sum_of_squares([1,2,3]) == 14);
+	assert(sum_of_squares([1,2,4]) == 21);
+	assert(sum_of_squares([-3,-2,-1]) == 14);
 }
 test_sum_of_squares();
 
 
 module test_sum_of_sines() {
-    assert_equal(sum_of_sines(0, [[3,4,0],[2,2,0]]), 0);
-    assert_equal(sum_of_sines(45, [[3,4,0],[2,2,0]]), 2);
-    assert_equal(sum_of_sines(90, [[3,4,0],[2,2,0]]), 0);
-    assert_equal(sum_of_sines(135, [[3,4,0],[2,2,0]]), -2);
-    assert_equal(sum_of_sines(180, [[3,4,0],[2,2,0]]), 0);
+	assert(sum_of_sines(0, [[3,4,0],[2,2,0]]) == 0);
+	assert(sum_of_sines(45, [[3,4,0],[2,2,0]]) == 2);
+	assert(sum_of_sines(90, [[3,4,0],[2,2,0]]) == 0);
+	assert(sum_of_sines(135, [[3,4,0],[2,2,0]]) == -2);
+	assert(sum_of_sines(180, [[3,4,0],[2,2,0]]) == 0);
 }
 test_sum_of_sines();
 
 
-module test_deltas() {
-    assert_equal(deltas([2,5,9,17]), [3,4,8]);
-    assert_equal(deltas([[1,2,3], [3,6,8], [4,8,11]]), [[2,4,5], [1,2,3]]);
-}
-test_deltas();
-
-
-module test_product() {
-    assert_equal(product([2,3,4]), 24);
-    assert_equal(product([[1,2,3], [3,4,5], [5,6,7]]), [15, 48, 105]);
-    m1 = [[2,3,4],[4,5,6],[6,7,8]];
-    m2 = [[4,1,2],[3,7,2],[8,7,4]];
-    m3 = [[3,7,8],[9,2,4],[5,8,3]];
-    assert_equal(product([m1,m2,m3]), m1*m2*m3);
-}
-test_product();
-
-
 module test_mean() {
-    assert_equal(mean([2,3,4]), 3);
-    assert_equal(mean([[1,2,3], [3,4,5], [5,6,7]]), [3,4,5]);
+	assert(mean([2,3,4]) == 3);
+	assert(mean([[1,2,3], [3,4,5], [5,6,7]]) == [3,4,5]);
 }
 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]);
-}
-test_median();
-
-
-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();
-
-
 // Logic
 
 
 module test_compare_vals() {
-    assert(compare_vals(-10,0) < 0);
-    assert(compare_vals(10,0) > 0);
-    assert(compare_vals(10,10) == 0);
+	assert(compare_vals(-10,0) == -1);
+	assert(compare_vals(10,0) == 1);
+	assert(compare_vals(10,10) == 0);
 
-    assert(compare_vals("abc","abcd") < 0);
-    assert(compare_vals("abcd","abc") > 0);
-    assert(compare_vals("abcd","abcd") == 0);
+	assert(compare_vals("abc","abcd") == -1);
+	assert(compare_vals("abcd","abc") == 1);
+	assert(compare_vals("abcd","abcd") == 0);
 
-    assert(compare_vals(false,false) == 0);
-    assert(compare_vals(true,false) > 0);
-    assert(compare_vals(false,true) < 0);
-    assert(compare_vals(true,true) == 0);
+	assert(compare_vals(false,false) == 0);
+	assert(compare_vals(true,false) == 1);
+	assert(compare_vals(false,true) == -1);
+	assert(compare_vals(true,true) == 0);
 
-    assert(compare_vals([2,3,4], [2,3,4,5]) < 0);
-    assert(compare_vals([2,3,4,5], [2,3,4,5]) == 0);
-    assert(compare_vals([2,3,4,5], [2,3,4]) > 0);
-    assert(compare_vals([2,3,4,5], [2,3,5,5]) < 0);
-    assert(compare_vals([[2,3,4,5]], [[2,3,5,5]]) < 0);
+	assert(compare_vals([2,3,4], [2,3,4,5]) == -1);
+	assert(compare_vals([2,3,4,5], [2,3,4,5]) == 0);
+	assert(compare_vals([2,3,4,5], [2,3,4]) == 1);
+	assert(compare_vals([2,3,4,5], [2,3,5,5]) == -1);
+	assert(compare_vals([[2,3,4,5]], [[2,3,5,5]]) == -1);
 
-    assert(compare_vals([[2,3,4],[3,4,5]], [[2,3,4], [3,4,5]]) == 0);
-    assert(compare_vals([[2,3,4],[3,4,5]], [[2,3,4,5], [3,4,5]]) < 0);
-    assert(compare_vals([[2,3,4],[3,4,5]], [[2,3,4], [3,4,5,6]]) < 0);
-    assert(compare_vals([[2,3,4,5],[3,4,5]], [[2,3,4], [3,4,5]]) > 0);
-    assert(compare_vals([[2,3,4],[3,4,5,6]], [[2,3,4], [3,4,5]]) > 0);
-    assert(compare_vals([[2,3,4],[3,5,5]], [[2,3,4], [3,4,5]]) > 0);
-    assert(compare_vals([[2,3,4],[3,4,5]], [[2,3,4], [3,5,5]]) < 0);
-
-    assert(compare_vals(undef, undef) == 0);
-    assert(compare_vals(undef, true) < 0);
-    assert(compare_vals(undef, 0) < 0);
-    assert(compare_vals(undef, "foo") < 0);
-    assert(compare_vals(undef, [2,3,4]) < 0);
-    assert(compare_vals(undef, [0:3]) < 0);
-
-    assert(compare_vals(true, undef) > 0);
-    assert(compare_vals(true, true) == 0);
-    assert(compare_vals(true, 0) < 0);
-    assert(compare_vals(true, "foo") < 0);
-    assert(compare_vals(true, [2,3,4]) < 0);
-    assert(compare_vals(true, [0:3]) < 0);
-
-    assert(compare_vals(0, undef) > 0);
-    assert(compare_vals(0, true) > 0);
-    assert(compare_vals(0, 0) == 0);
-    assert(compare_vals(0, "foo") < 0);
-    assert(compare_vals(0, [2,3,4]) < 0);
-    assert(compare_vals(0, [0:3]) < 0);
-
-    assert(compare_vals(1, undef) > 0);
-    assert(compare_vals(1, true) > 0);
-    assert(compare_vals(1, 1) == 0);
-    assert(compare_vals(1, "foo") < 0);
-    assert(compare_vals(1, [2,3,4]) < 0);
-    assert(compare_vals(1, [0:3]) < 0);
-
-    assert(compare_vals("foo", undef) > 0);
-    assert(compare_vals("foo", true) > 0);
-    assert(compare_vals("foo", 1) > 0);
-    assert(compare_vals("foo", "foo") == 0);
-    assert(compare_vals("foo", [2,3,4]) < 0);
-    assert(compare_vals("foo", [0:3]) < 0);
-
-    assert(compare_vals([2,3,4], undef) > 0);
-    assert(compare_vals([2,3,4], true) > 0);
-    assert(compare_vals([2,3,4], 1) > 0);
-    assert(compare_vals([2,3,4], "foo") > 0);
-    assert(compare_vals([2,3,4], [2,3,4]) == 0);
-    assert(compare_vals([2,3,4], [0:3]) < 0);
-
-    assert(compare_vals([0:3], undef) > 0);
-    assert(compare_vals([0:3], true) > 0);
-    assert(compare_vals([0:3], 1) > 0);
-    assert(compare_vals([0:3], "foo") > 0);
-    assert(compare_vals([0:3], [2,3,4]) > 0);
-    assert(compare_vals([0:3], [0:3]) == 0);
+	assert(compare_vals([[2,3,4],[3,4,5]], [[2,3,4], [3,4,5]]) == 0);
+	assert(compare_vals([[2,3,4],[3,4,5]], [[2,3,4,5], [3,4,5]]) == -1);
+	assert(compare_vals([[2,3,4],[3,4,5]], [[2,3,4], [3,4,5,6]]) == -1);
+	assert(compare_vals([[2,3,4,5],[3,4,5]], [[2,3,4], [3,4,5]]) == 1);
+	assert(compare_vals([[2,3,4],[3,4,5,6]], [[2,3,4], [3,4,5]]) == 1);
+	assert(compare_vals([[2,3,4],[3,5,5]], [[2,3,4], [3,4,5]]) == 1);
+	assert(compare_vals([[2,3,4],[3,4,5]], [[2,3,4], [3,5,5]]) == -1);
 }
 test_compare_vals();
 
 
 module test_compare_lists() {
-    assert(compare_lists([2,3,4], [2,3,4,5]) < 0);
-    assert(compare_lists([2,3,4,5], [2,3,4,5]) == 0);
-    assert(compare_lists([2,3,4,5], [2,3,4]) > 0);
-    assert(compare_lists([2,3,4,5], [2,3,5,5]) < 0);
+	assert(compare_lists([2,3,4], [2,3,4,5]) == -1);
+	assert(compare_lists([2,3,4,5], [2,3,4,5]) == 0);
+	assert(compare_lists([2,3,4,5], [2,3,4]) == 1);
+	assert(compare_lists([2,3,4,5], [2,3,5,5]) == -1);
 
-    assert(compare_lists([[2,3,4],[3,4,5]], [[2,3,4], [3,4,5]]) == 0);
-    assert(compare_lists([[2,3,4],[3,4,5]], [[2,3,4,5], [3,4,5]]) < 0);
-    assert(compare_lists([[2,3,4],[3,4,5]], [[2,3,4], [3,4,5,6]]) < 0);
-    assert(compare_lists([[2,3,4,5],[3,4,5]], [[2,3,4], [3,4,5]]) > 0);
-    assert(compare_lists([[2,3,4],[3,4,5,6]], [[2,3,4], [3,4,5]]) > 0);
-    assert(compare_lists([[2,3,4],[3,5,5]], [[2,3,4], [3,4,5]]) > 0);
-    assert(compare_lists([[2,3,4],[3,4,5]], [[2,3,4], [3,5,5]]) < 0);
+	assert(compare_lists([[2,3,4],[3,4,5]], [[2,3,4], [3,4,5]]) == 0);
+	assert(compare_lists([[2,3,4],[3,4,5]], [[2,3,4,5], [3,4,5]]) == -1);
+	assert(compare_lists([[2,3,4],[3,4,5]], [[2,3,4], [3,4,5,6]]) == -1);
+	assert(compare_lists([[2,3,4,5],[3,4,5]], [[2,3,4], [3,4,5]]) == 1);
+	assert(compare_lists([[2,3,4],[3,4,5,6]], [[2,3,4], [3,4,5]]) == 1);
+	assert(compare_lists([[2,3,4],[3,5,5]], [[2,3,4], [3,4,5]]) == 1);
+	assert(compare_lists([[2,3,4],[3,4,5]], [[2,3,4], [3,5,5]]) == -1);
 
-    assert(compare_lists("cat", "bat") > 0);
-    assert(compare_lists(["cat"], ["bat"]) > 0);
+	assert(compare_lists("cat", "bat") == 1);
+	assert(compare_lists(["cat"], ["bat"]) == 1);
 }
 test_compare_lists();
 
 
 module test_any() {
-    assert_equal(any([0,false,undef]), false);
-    assert_equal(any([1,false,undef]), true);
-    assert_equal(any([1,5,true]), true);
-    assert_equal(any([[0,0], [0,0]]), false);
-    assert_equal(any([[0,0], [1,0]]), true);
+	assert(any([0,false,undef]) == false);
+	assert(any([1,false,undef]) == true);
+	assert(any([1,5,true]) == true);
+	assert(any([[0,0], [0,0]]) == false);
+	assert(any([[0,0], [1,0]]) == true);
 }
 test_any();
 
 
 module test_all() {
-    assert_equal(all([0,false,undef]), false);
-    assert_equal(all([1,false,undef]), false);
-    assert_equal(all([1,5,true]), true);
-    assert_equal(all([[0,0], [0,0]]), false);
-    assert_equal(all([[0,0], [1,0]]), false);
-    assert_equal(all([[1,1], [1,1]]), true);
+	assert(all([0,false,undef]) == false);
+	assert(all([1,false,undef]) == false);
+	assert(all([1,5,true]) == true);
+	assert(all([[0,0], [0,0]]) == false);
+	assert(all([[0,0], [1,0]]) == false);
+	assert(all([[1,1], [1,1]]) == true);
 }
 test_all();
 
 
 module test_count_true() {
-    assert_equal(count_true([0,false,undef]), 0);
-    assert_equal(count_true([1,false,undef]), 1);
-    assert_equal(count_true([1,5,false]), 2);
-    assert_equal(count_true([1,5,true]), 3);
-    assert_equal(count_true([[0,0], [0,0]]), 0);
-    assert_equal(count_true([[0,0], [1,0]]), 1);
-    assert_equal(count_true([[1,1], [1,1]]), 4);
-    assert_equal(count_true([[1,1], [1,1]], nmax=3), 3);
+	assert(count_true([0,false,undef]) == 0);
+	assert(count_true([1,false,undef]) == 1);
+	assert(count_true([1,5,false]) == 2);
+	assert(count_true([1,5,true]) == 3);
+	assert(count_true([[0,0], [0,0]]) == 0);
+	assert(count_true([[0,0], [1,0]]) == 1);
+	assert(count_true([[1,1], [1,1]]) == 4);
+	assert(count_true([[1,1], [1,1]], nmax=3) == 3);
 }
 test_count_true();
 
 
-module test_factorial() {
-    assert_equal(factorial(1), 1);
-    assert_equal(factorial(2), 2);
-    assert_equal(factorial(3), 6);
-    assert_equal(factorial(4), 24);
-    assert_equal(factorial(5), 120);
-    assert_equal(factorial(6), 720);
-    assert_equal(factorial(7), 5040);
-    assert_equal(factorial(8), 40320);
+
+// List/Array Ops
+
+module test_cdr() {
+	assert(cdr([]) == []);
+	assert(cdr([88]) == []);
+	assert(cdr([1,2,3]) == [2,3]);
+	assert(cdr(["a","b","c"]) == ["b","c"]);
 }
-test_factorial();
+test_cdr();
 
 
-module test_gcd() {
-    assert_equal(gcd(15,25), 5);
-    assert_equal(gcd(15,27), 3);
-    assert_equal(gcd(270,405), 135);
+module test_replist() {
+	assert(replist(1, 4) == [1,1,1,1]);
+	assert(replist(8, [2,3]) == [[8,8,8], [8,8,8]]);
+	assert(replist(0, [2,2,3]) == [[[0,0,0],[0,0,0]], [[0,0,0],[0,0,0]]]);
+	assert(replist([1,2,3],3) == [[1,2,3], [1,2,3], [1,2,3]]);
 }
-test_gcd();
+test_replist();
 
 
-module test_lcm() {
-    assert_equal(lcm(15,25), 75);
-    assert_equal(lcm(15,27), 135);
-    assert_equal(lcm(270,405), 810);
+module test_in_list() {
+	assert(in_list("bar", ["foo", "bar", "baz"]));
+	assert(!in_list("bee", ["foo", "bar", "baz"]));
+	assert(in_list("bar", [[2,"foo"], [4,"bar"], [3,"baz"]], idx=1));
 }
-test_lcm();
+test_in_list();
 
 
+module test_slice() {
+	assert(slice([3,4,5,6,7,8,9], 3, 5) == [6,7]);
+	assert(slice([3,4,5,6,7,8,9], 2, -1) == [5,6,7,8,9]);
+	assert(slice([3,4,5,6,7,8,9], 1, 1) == []);
+	assert(slice([3,4,5,6,7,8,9], 6, -1) == [9]);
+	assert(slice([3,4,5,6,7,8,9], 2, -2) == [5,6,7,8]);
+}
+test_slice();
 
-// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
+
+module test_select() {
+	l = [3,4,5,6,7,8,9];
+	assert(select(l, 5, 6) == [8,9]);
+	assert(select(l, 5, 8) == [8,9,3,4]);
+	assert(select(l, 5, 2) == [8,9,3,4,5]);
+	assert(select(l, -3, -1) == [7,8,9]);
+	assert(select(l, 3, 3) == [6]);
+	assert(select(l, 4) == 7);
+	assert(select(l, -2) == 8);
+	assert(select(l, [1:3]) == [4,5,6]);
+	assert(select(l, [1,3]) == [4,6]);
+}
+test_select();
+
+
+module test_reverse() {
+	assert(reverse([3,4,5,6]) == [6,5,4,3]);
+}
+test_reverse();
+
+
+module test_array_subindex() {
+	v = [[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]];
+	assert(array_subindex(v,2) == [3, 7, 11, 15]);
+	assert(array_subindex(v,[2,1]) == [[3, 2], [7, 6], [11, 10], [15, 14]]);
+	assert(array_subindex(v,[1:3]) == [[2, 3, 4], [6, 7, 8], [10, 11, 12], [14, 15, 16]]);
+}
+test_array_subindex();
+
+
+module test_list_range() {
+	assert(list_range(4) == [0,1,2,3]);
+	assert(list_range(n=4, step=2) == [0,2,4,6]);
+	assert(list_range(n=4, s=3, step=3) == [3,6,9,12]);
+	assert(list_range(n=4, s=3, e=9, step=3) == [3,6,9]);
+	assert(list_range(e=3) == [0,1,2,3]);
+	assert(list_range(e=6, step=2) == [0,2,4,6]);
+	assert(list_range(s=3, e=5) == [3,4,5]);
+	assert(list_range(s=3, e=8, step=2) == [3,5,7]);
+	assert(list_range(s=4, e=8, step=2) == [4,6,8]);
+	assert(list_range(n=4, s=[3,4], step=[2,3]) == [[3,4], [5,7], [7,10], [9,13]]);
+}
+test_list_range();
+
+
+module test_array_shortest() {
+	assert(array_shortest(["foobar", "bazquxx", "abcd"]) == 4);
+}
+test_array_shortest();
+
+
+module test_array_longest() {
+	assert(array_longest(["foobar", "bazquxx", "abcd"]) == 7);
+}
+test_array_longest();
+
+
+module test_array_pad() {
+	assert(array_pad([4,5,6], 5, 8) == [4,5,6,8,8]);
+	assert(array_pad([4,5,6,7,8], 5, 8) == [4,5,6,7,8]);
+	assert(array_pad([4,5,6,7,8,9], 5, 8) == [4,5,6,7,8,9]);
+}
+test_array_pad();
+
+
+module test_array_trim() {
+	assert(array_trim([4,5,6], 5) == [4,5,6]);
+	assert(array_trim([4,5,6,7,8], 5) == [4,5,6,7,8]);
+	assert(array_trim([3,4,5,6,7,8,9], 5) == [3,4,5,6,7]);
+}
+test_array_trim();
+
+
+module test_array_fit() {
+	assert(array_fit([4,5,6], 5, 8) == [4,5,6,8,8]);
+	assert(array_fit([4,5,6,7,8], 5, 8) == [4,5,6,7,8]);
+	assert(array_fit([3,4,5,6,7,8,9], 5, 8) == [3,4,5,6,7]);
+}
+test_array_fit();
+
+
+module test_enumerate() {
+	assert(enumerate(["a","b","c"]) == [[0,"a"], [1,"b"], [2,"c"]]);
+	assert(enumerate([[88,"a"],[76,"b"],[21,"c"]], idx=1) == [[0,"a"], [1,"b"], [2,"c"]]);
+	assert(enumerate([["cat","a",12],["dog","b",10],["log","c",14]], idx=[1:2]) == [[0,"a",12], [1,"b",10], [2,"c",14]]);
+}
+test_enumerate();
+
+
+module test_array_zip() {
+	v1 = [1,2,3,4];
+	v2 = [5,6,7];
+	v3 = [8,9,10,11];
+	assert(array_zip(v1,v3) == [[1,8],[2,9],[3,10],[4,11]]);
+	assert(array_zip([v1,v3]) == [[1,8],[2,9],[3,10],[4,11]]);
+	assert(array_zip([v1,v2],fit="short") == [[1,5],[2,6],[3,7]]);
+	assert(array_zip([v1,v2],fit="long") == [[1,5],[2,6],[3,7],[4,undef]]);
+	assert(array_zip([v1,v2],fit="long", fill=0) == [[1,5],[2,6],[3,7],[4,0]]);
+	assert(array_zip([v1,v2,v3],fit="long") == [[1,5,8],[2,6,9],[3,7,10],[4,undef,11]]);
+}
+test_array_zip();
+
+
+module test_array_group() {
+	v = [1,2,3,4,5,6];
+	assert(array_group(v,2) == [[1,2], [3,4], [5,6]]);
+	assert(array_group(v,3) == [[1,2,3], [4,5,6]]);
+	assert(array_group(v,4,0) == [[1,2,3,4], [5,6,0,0]]);
+}
+test_array_group();
+
+
+module test_flatten() {
+	assert(flatten([[1,2,3], [4,5,[6,7,8]]]) == [1,2,3,4,5,[6,7,8]]);
+}
+test_flatten();
+
+
+module test_sort() {
+	assert(sort([7,3,9,4,3,1,8]) == [1,3,3,4,7,8,9]);
+	assert(sort(["cat", "oat", "sat", "bat", "vat", "rat", "pat", "mat", "fat", "hat", "eat"]) == ["bat", "cat", "eat", "fat", "hat", "mat", "oat", "pat", "rat", "sat", "vat"]);
+	assert(sort(enumerate([[2,3,4],[1,2,3],[2,4,3]]),idx=1)==[[1,[1,2,3]], [0,[2,3,4]], [2,[2,4,3]]]);
+}
+test_sort();
+
+
+module test_sortidx() {
+	lst1 = ["d","b","e","c"];
+	assert(sortidx(lst1) == [1,3,0,2]);
+	lst2 = [
+		["foo", 88, [0,0,1], false],
+		["bar", 90, [0,1,0], true],
+		["baz", 89, [1,0,0], false],
+		["qux", 23, [1,1,1], true]
+	];
+	assert(sortidx(lst2, idx=1) == [3,0,2,1]);
+	assert(sortidx(lst2, idx=0) == [1,2,0,3]);
+	assert(sortidx(lst2, idx=[1,3]) == [3,0,2,1]);
+	lst3 = [[-4, 0, 0], [0, 0, -4], [0, -4, 0], [-4, 0, 0], [0, -4, 0], [0, 0, 4], [0, 0, -4], [0, 4, 0], [4, 0, 0], [0, 0, 4], [0, 4, 0], [4, 0, 0]];
+	assert(sortidx(lst3)==[0,3,2,4,1,6,5,9,7,10,8,11]);
+}
+test_sortidx();
+
+
+module test_unique() {
+	assert(unique([]) == []);
+	assert(unique([8]) == [8]);
+	assert(unique([7,3,9,4,3,1,8]) == [1,3,4,7,8,9]);
+}
+test_unique();
+
+
+module test_array_dim() {
+	assert(array_dim([[[1,2,3],[4,5,6]],[[7,8,9],[10,11,12]]]) == [2,2,3]);
+	assert(array_dim([[[1,2,3],[4,5,6]],[[7,8,9],[10,11,12]]], 0) == 2);
+	assert(array_dim([[[1,2,3],[4,5,6]],[[7,8,9],[10,11,12]]], 2) == 3);
+	assert(array_dim([[[1,2,3],[4,5,6]],[[7,8,9]]]) == [2,undef,3]);
+}
+test_array_dim();
+
+
+module test_vmul() {
+	assert(vmul([3,4,5], [8,7,6]) == [24,28,30]);
+	assert(vmul([1,2,3], [4,5,6]) == [4,10,18]);
+}
+test_vmul();
+
+
+module test_vdiv() {
+	assert(vdiv([24,28,30], [8,7,6]) == [3, 4, 5]);
+}
+test_vdiv();
+
+
+module test_vabs() {
+	assert(vabs([2,4,8]) == [2,4,8]);
+	assert(vabs([-2,-4,-8]) == [2,4,8]);
+	assert(vabs([-2,4,8]) == [2,4,8]);
+	assert(vabs([2,-4,8]) == [2,4,8]);
+	assert(vabs([2,4,-8]) == [2,4,8]);
+}
+test_vabs();
+
+
+module test_normalize() {
+	assert(normalize([10,0,0]) == [1,0,0]);
+	assert(normalize([0,10,0]) == [0,1,0]);
+	assert(normalize([0,0,10]) == [0,0,1]);
+	assert(abs(norm(normalize([10,10,10]))-1) < eps);
+	assert(abs(norm(normalize([-10,-10,-10]))-1) < eps);
+	assert(abs(norm(normalize([-10,0,0]))-1) < eps);
+	assert(abs(norm(normalize([0,-10,0]))-1) < eps);
+	assert(abs(norm(normalize([0,0,-10]))-1) < eps);
+}
+test_normalize();
+
+
+module test_vector_angle() {
+	vecs = [[10,0,0], [-10,0,0], [0,10,0], [0,-10,0], [0,0,10], [0,0,-10]];
+	for (a=vecs, b=vecs) {
+		if(a==b) {
+			assert(vector_angle(a,b)==0);
+		} else if(a==-b) {
+			assert(vector_angle(a,b)==180);
+		} else {
+			assert(vector_angle(a,b)==90);
+		}
+	}
+	assert(abs(vector_angle([10,10,0],[10,0,0])-45) < eps);
+}
+test_vector_angle();
+
+
+module test_vector_axis() {
+	assert(norm(vector_axis([10,0,0],[10,10,0]) - [0,0,1]) < eps);
+	assert(norm(vector_axis([10,0,0],[0,10,0]) - [0,0,1]) < eps);
+	assert(norm(vector_axis([0,10,0],[10,0,0]) - [0,0,-1]) < eps);
+	assert(norm(vector_axis([0,0,10],[10,0,0]) - [0,1,0]) < eps);
+	assert(norm(vector_axis([10,0,0],[0,0,10]) - [0,-1,0]) < eps);
+	assert(norm(vector_axis([10,0,10],[0,-10,0]) - [sin(45),0,-sin(45)]) < eps);
+}
+test_vector_axis();
+
+
+module test_point2d() {
+	assert(point2d([1,2,3])==[1,2]);
+	assert(point2d([2,3])==[2,3]);
+	assert(point2d([1])==[1,0]);
+}
+test_point2d();
+
+
+module test_path2d() {
+	assert(path2d([[1], [1,2], [1,2,3], [1,2,3,4], [1,2,3,4,5]])==[[1,0],[1,2],[1,2],[1,2],[1,2]]);
+}
+test_path2d();
+
+
+module test_point3d() {
+	assert(point3d([1,2,3,4,5])==[1,2,3]);
+	assert(point3d([1,2,3,4])==[1,2,3]);
+	assert(point3d([1,2,3])==[1,2,3]);
+	assert(point3d([2,3])==[2,3,0]);
+	assert(point3d([1])==[1,0,0]);
+}
+test_point3d();
+
+
+module test_path3d() {
+	assert(path3d([[1], [1,2], [1,2,3], [1,2,3,4], [1,2,3,4,5]])==[[1,0,0],[1,2,0],[1,2,3],[1,2,3],[1,2,3]]);
+}
+test_path3d();
+
+
+module test_translate_points() {
+	pts = [[0,0,1], [0,1,0], [1,0,0], [0,0,-1], [0,-1,0], [-1,0,0]];
+	assert(translate_points(pts, v=[1,2,3]) == [[1,2,4], [1,3,3], [2,2,3], [1,2,2], [1,1,3], [0,2,3]]);
+	assert(translate_points(pts, v=[-1,-2,-3]) == [[-1,-2,-2], [-1,-1,-3], [0,-2,-3], [-1,-2,-4], [-1,-3,-3], [-2,-2,-3]]);
+}
+test_translate_points();
+
+
+module test_scale_points() {
+	pts = [[0,0,1], [0,1,0], [1,0,0], [0,0,-1], [0,-1,0], [-1,0,0]];
+	assert(scale_points(pts, v=[2,3,4]) == [[0,0,4], [0,3,0], [2,0,0], [0,0,-4], [0,-3,0], [-2,0,0]]);
+	assert(scale_points(pts, v=[-2,-3,-4]) == [[0,0,-4], [0,-3,0], [-2,0,0], [0,0,4], [0,3,0], [2,0,0]]);
+	assert(scale_points(pts, v=[1,1,1]) == [[0,0,1], [0,1,0], [1,0,0], [0,0,-1], [0,-1,0], [-1,0,0]]);
+	assert(scale_points(pts, v=[-1,-1,-1]) == [[0,0,-1], [0,-1,0], [-1,0,0], [0,0,1], [0,1,0], [1,0,0]]);
+}
+test_scale_points();
+
+
+module test_rotate_points2d() {
+	pts = [[0,1], [1,0], [0,-1], [-1,0]];
+	s = sin(45);
+	assert(rotate_points2d(pts,45) == [[-s,s],[s,s],[s,-s],[-s,-s]]);
+	assert(rotate_points2d(pts,90) == [[-1,0],[0,1],[1,0],[0,-1]]);
+	assert(rotate_points2d(pts,90,cp=[1,0]) == [[0,-1],[1,0],[2,-1],[1,-2]]);
+}
+test_rotate_points2d();
+
+
+module test_rotate_points3d() {
+	pts = [[0,0,1], [0,1,0], [1,0,0], [0,0,-1], [0,-1,0], [-1,0,0]];
+	assert(rotate_points3d(pts, [90,0,0]) == [[0,-1,0], [0,0,1], [1,0,0], [0,1,0], [0,0,-1], [-1,0,0]]);
+	assert(rotate_points3d(pts, [0,90,0]) == [[1,0,0], [0,1,0], [0,0,-1], [-1,0,0], [0,-1,0], [0,0,1]]);
+	assert(rotate_points3d(pts, [0,0,90]) == [[0,0,1], [-1,0,0], [0,1,0], [0,0,-1], [1,0,0], [0,-1,0]]);
+	assert(rotate_points3d(pts, [0,0,90],cp=[2,0,0]) == [[2,-2,1], [1,-2,0], [2,-1,0], [2,-2,-1], [3,-2,0], [2,-3,0]]);
+	assert(rotate_points3d(pts, 90, axis=V_UP) == [[0,0,1], [-1,0,0], [0,1,0], [0,0,-1], [1,0,0], [0,-1,0]]);
+	assert(rotate_points3d(pts, 90, axis=V_DOWN) == [[0,0,1], [1,0,0], [0,-1,0], [0,0,-1], [-1,0,0], [0,1,0]]);
+	assert(rotate_points3d(pts, 90, axis=V_RIGHT)  == [[0,-1,0], [0,0,1], [1,0,0], [0,1,0], [0,0,-1], [-1,0,0]]);
+	assert(rotate_points3d(pts, from=V_UP, to=V_BACK) == [[0,1,0], [0,0,-1], [1,0,0], [0,-1,0], [0,0,1], [-1,0,0]]);
+	assert(rotate_points3d(pts, 90, from=V_UP, to=V_BACK), [[0,1,0], [-1,0,0], [0,0,-1], [0,-1,0], [1,0,0], [0,0,1]]);
+	assert(rotate_points3d(pts, from=V_UP, to=V_UP*2) == [[0,0,1], [0,1,0], [1,0,0], [0,0,-1], [0,-1,0], [-1,0,0]]);
+	assert(rotate_points3d(pts, from=V_UP, to=V_DOWN*2) == [[0,0,-1], [0,1,0], [-1,0,0], [0,0,1], [0,-1,0], [1,0,0]]);
+}
+test_rotate_points3d();
+
+
+module test_simplify_path()
+{
+	path = [[-20,10],[-10,0],[-5,0],[0,0],[5,0],[10,0], [10,10]];
+	assert(simplify_path(path) == [[-20,10],[-10,0],[10,0], [10,10]]);
+}
+test_simplify_path();
+
+
+module test_simplify_path_indexed()
+{
+	points = [[-20,10],[-10,0],[-5,0],[0,0],[5,0],[10,0], [10,10]];
+	path = list_range(len(points));
+	assert(simplify_path_indexed(points, path) == [0,1,5,6]);
+}
+test_simplify_path_indexed();
+
+
+// vim: noexpandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap

From 3864f14333c36f2e91f99b6037b3b455def2b215 Mon Sep 17 00:00:00 2001
From: Adrian Mariano <avm4@cornell.edu>
Date: Sat, 11 Jul 2020 12:22:59 -0400
Subject: [PATCH 2/3] Added error bounds to poly_roots, added poly_div,
 poly_mult and poly_add, fixed bugs in factorial.

---
 math.scad | 147 +++++++++++++++++++++++++++++++++++++++++++++---------
 1 file changed, 124 insertions(+), 23 deletions(-)

diff --git a/math.scad b/math.scad
index 09e3b9c..7629cb1 100644
--- a/math.scad
+++ b/math.scad
@@ -67,7 +67,7 @@ function hypot(x,y,z=0) = norm([x,y,z]);
 // Usage:
 //   x = factorial(n,[d]);
 // Description:
-//   Returns the factorial of the given integer value.
+//   Returns the factorial of the given integer value, or n!/d! if d is given.  
 // Arguments:
 //   n = The integer number to get the factorial of.  (n!)
 //   d = If given, the returned value will be (n! / d!)
@@ -75,7 +75,10 @@ function hypot(x,y,z=0) = norm([x,y,z]);
 //   x = factorial(4);  // Returns: 24
 //   y = factorial(6);  // Returns: 720
 //   z = factorial(9);  // Returns: 362880
-function factorial(n,d=1) = product([for (i=[n:-1:d]) i]);
+function factorial(n,d=0) =
+    assert(n>=0 && d>=0, "Factorial is not defined for negative numbers")
+    assert(d<=n, "d cannot be larger than n")
+    product([1,for (i=[n:-1:d+1]) i]);
 
 
 // Function: lerp()
@@ -525,6 +528,17 @@ function deltas(v) = [for (p=pair(v)) p.y-p.x];
 function product(v, i=0, tot=undef) = i>=len(v)? tot : product(v, i+1, ((tot==undef)? v[i] : is_vector(v[i])? vmul(tot,v[i]) : tot*v[i]));
 
 
+// Function: outer_product()
+// 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))
+  assert(len(u)==len(v))
+  [for(i=[0:len(u)-1]) [for(j=[0:len(u)-1]) u[i]*v[j]]];
+
+
 // Function: mean()
 // Description:
 //   Returns the arithmatic mean/average of all entries in the given array.
@@ -563,7 +577,7 @@ function median(v) =
 // 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 `undef`.  If b is a matrix that is compatible with A
+//   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.  
@@ -582,7 +596,7 @@ function linear_solve(A,b) =
         R = submatrix(qr[1],[0:mindim-1], [0:mindim-1]),
         zeros = [for(i=[0:mindim-1]) if (approx(R[i][i],0)) i]
     )
-    zeros != [] ? undef :
+    zeros != [] ? [] :
     m<n ? Q*back_substitute(R,b,transpose=true) :
     back_substitute(R, transpose(Q)*b);
 
@@ -604,7 +618,8 @@ function matrix_inverse(A) =
 // Usage: submatrix(M, ind1, ind2)
 // Description:
 //   Returns a submatrix with the specified index ranges or index sets.  
-function submatrix(M,ind1,ind2) = [for(i=ind1) [for(j=ind2) M[i][j] ] ];
+function submatrix(M,ind1,ind2) =
+    [for(i=ind1) [for(j=ind2) M[i][j] ] ];
 
 
 // Function: qr_factor()
@@ -635,7 +650,7 @@ function _qr_factor(A,Q, column, m, n) =
         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*transpose([v])*[v],
+        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, column+1, m, n);
@@ -647,11 +662,12 @@ function _qr_factor(A,Q, column, m, n) =
 //   Solves the problem Rx=b where R is an upper triangular square matrix.  No check is made that the lower triangular entries
 //   are actually zero.  If transpose==true then instead solve transpose(R)*x=b.
 //   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. 
+//   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, x=[],transpose = false) =
     assert(is_matrix(R, square=true))
     let(n=len(R))
-    assert(is_vector(b,n) || is_matrix(b,n),"R and b are not compatible in back_substitute")
+    assert(is_vector(b,n) || is_matrix(b,n),str("R and b are not compatible in back_substitute ",n, len(b)))
     !is_vector(b) ? transpose([for(i=[0:len(b[0])-1]) back_substitute(R,subindex(b,i),transpose=transpose)]) :
     transpose?
         reverse(back_substitute(
@@ -660,7 +676,10 @@ function back_substitute(R, b, x=[],transpose = false) =
         )) :
     len(x) == n ? x :
     let(
-        ind = n - len(x) - 1,
+        ind = n - len(x) - 1
+    )
+    R[ind][ind] == 0 ? [] :
+    let(
         newvalue =
             len(x)==0? b[ind]/R[ind][ind] : 
             (b[ind]-select(R[ind],ind+1,-1) * x)/R[ind][ind]
@@ -733,7 +752,7 @@ function determinant(M) =
 //   m = optional height of matrix
 //   n = optional width of matrix
 //   square = set to true to require a square matrix.  Default: false        
-function is_matrix(A,n,m,square=false) =
+function is_matrix(A,m,n,square=false) =
     is_vector(A[0],n) && is_vector(A*(0*A[0]),m) &&
     (!square || len(A)==len(A[0]));
 
@@ -1037,12 +1056,82 @@ function C_div(z1,z2) = let(den = z2.x*z2.x + z2.y*z2.y)
 //   Evaluates specified real polynomial, p, at the complex or real input value, z.
 //   The polynomial is specified as p=[a_n, a_{n-1},...,a_1,a_0]
 //   where a_n is the z^n coefficient.  Polynomial coefficients are real.
+
+// Note: this should probably be recoded to use division by [1,-z], which is more accurate
+// and avoids overflow with large coefficients, but requires poly_div to support complex coefficients.  
 function polynomial(p, z, k, zk, total) =
    is_undef(k) ? polynomial(p, z, len(p)-1, is_num(z)? 1 : [1,0], is_num(z) ? 0 : [0,0]) :
    k==-1 ? total :
    polynomial(p, z, k-1, is_num(z) ? zk*z : C_times(zk,z), total+zk*p[k]);
 
 
+// Function: poly_mult()
+// Usage
+//   polymult(p,q)
+//   polymult([p1,p2,p3,...])
+// Descriptoin:
+//   Given a list of polynomials represented as real coefficient lists, with the highest degree coefficient first, 
+//   computes the coefficient list of the product polynomial.  
+function poly_mult(p,q) = 
+  is_undef(q) ?
+     assert(is_list(p) && (is_vector(p[0]) || p[0]==[]), "Invalid arguments to poly_mult")
+     len(p)==2 ? poly_mult(p[0],p[1]) 
+               : poly_mult(p[0], poly_mult(select(p,1,-1)))
+  :
+  _poly_trim(
+  [
+  for(n = [len(p)+len(q)-2:-1:0])
+      sum( [for(i=[0:1:len(p)-1])
+           let(j = len(p)+len(q)- 2 - n - i)
+           if (j>=0 && j<len(q)) p[i]*q[j]
+               ])
+   ]);        
+
+
+// Function: poly_div()
+// Usage:
+//    [quotient,remainder] = poly_div(n,d)
+// Description:
+//    Computes division of the numerator polynomial by the denominator polynomial and returns
+//    a list of two polynomials, [quotient, remainder].  If the division has no remainder then
+//    the zero polynomial [] is returned for the remainder.  Similarly if the quotient is zero
+//    the returned quotient will be [].  
+function poly_div(n,d,q=[]) =
+    assert(len(d)>0 && d[0]!=0 , "Denominator is zero or has leading zero coefficient")
+    len(n)<len(d) ? [q,_poly_trim(n)] : 
+    let(
+      t = n[0] / d[0],
+      newq = concat(q,[t]),
+      newn =  [for(i=[1:1:len(n)-1]) i<len(d) ? n[i] - t*d[i] : n[i]]
+    )  
+    poly_div(newn,d,newq);
+
+
+// Internal Function: _poly_trim()
+// Usage:
+//    _poly_trim(p,[eps])
+// Description:
+//    Removes leading zero terms of a polynomial.  By default zeros must be exact,
+//    or give epsilon for approximate zeros.  
+function _poly_trim(p,eps=0) =
+  let(  nz = [for(i=[0:1:len(p)-1]) if (!approx(p[i],0,eps)) i])
+  len(nz)==0 ? [] : select(p,nz[0],-1);
+
+
+// Function: poly_add()
+// Usage:
+//    sum = poly_add(p,q)
+// Description:
+//    Computes the sum of two polynomials.  
+function poly_add(p,q) = 
+  let(  plen = len(p),
+        qlen = len(q),
+        long = plen>qlen ? p : q,
+        short = plen>qlen ? q : p
+     )
+   _poly_trim(long + concat(repeat(0,len(long)-len(short)),short));
+
+
 // Function: poly_roots()
 // Usage:
 //   poly_roots(p,[tol])
@@ -1062,14 +1151,18 @@ function polynomial(p, z, k, zk, total) =
 // Dario Bini. "Numerical computation of polynomial zeros by means of Aberth's Method", Numerical Algorithms, Feb 1996.
 // https://www.researchgate.net/publication/225654837_Numerical_computation_of_polynomial_zeros_by_means_of_Aberth's_method
 
-function poly_roots(p,tol=1e-14) =
+function poly_roots(p,tol=1e-14,error_bound=false) =
   assert(p!=[], "Input polynomial must have a nonzero coefficient")
   assert(is_vector(p), "Input must be a vector")
-  p[0] == 0 ? poly_roots(slice(p,1,-1)) :    // Strip leading zero coefficients
-  p[len(p)-1] == 0 ? [[0,0],                 // Strip trailing zero coefficients
-                      each poly_roots(select(p,0,-2))] :
-  len(p)==1 ? [] :                  // Nonzero constant case has no solutions
-  len(p)==2 ? [[-p[1]/p[0],0]] :    // Linear case needs special handling
+  p[0] == 0 ? poly_roots(slice(p,1,-1),tol=tol,error_bound=error_bound) :    // Strip leading zero coefficients
+  p[len(p)-1] == 0 ?                                       // Strip trailing zero coefficients
+      let( solutions = poly_roots(select(p,0,-2),tol=tol, error_bound=error_bound))
+      (error_bound ? [ [[0,0], each solutions[0]], [0, each solutions[1]]]
+                  : [[0,0], each solutions]) :
+  len(p)==1 ? (error_bound ? [[],[]] : []) :               // Nonzero constant case has no solutions
+  len(p)==2 ? let( solution = [[-p[1]/p[0],0]])            // Linear case needs special handling
+              (error_bound ? [solution,[0]] : solution)
+  : 
   let(
       n = len(p)-1,   // polynomial degree
       pderiv = [for(i=[0:n-1]) p[i]*(n-i)],
@@ -1082,9 +1175,12 @@ function poly_roots(p,tol=1e-14) =
       init = [for(i=[0:1:n-1])                // Initial guess for roots       
                let(angle = 360*i/n+270/n/PI)
                [beta,0]+r*[cos(angle),sin(angle)]
-             ]
+             ],
+      roots = _poly_roots(p,pderiv,s,init,tol=tol),
+      error = error_bound ? [for(xi=roots) n * (norm(polynomial(p,xi))+tol*polynomial(s,norm(xi))) /
+                                abs(norm(polynomial(pderiv,xi))-tol*polynomial(s,norm(xi)))] : 0
     )
-   _poly_roots(p,pderiv,s,init,tol=tol);
+    error_bound ? [roots, error] : roots;
 
 // p = polynomial
 // pderiv = derivative polynomial of p
@@ -1114,7 +1210,10 @@ function _poly_roots(p, pderiv, s, z, tol, i=0) =
 //   The polynomial is specified as p=[a_n, a_{n-1},...,a_1,a_0]
 //   where a_n is the x^n coefficient.  This function works by
 //   computing the complex roots and returning those roots where
-//   the imaginary part is closed to zero, specifically: z.y/(1+norm(z)) < eps.  Because
+//   the imaginary part is closed to zero.  By default it uses a computed
+//   error bound from the polynomial solver to decide whether imaginary
+//   parts are zero.  You can specify eps, in which case the test is
+//   z.y/(1+norm(z)) < eps.  Because
 //   of poor convergence and higher error for repeated roots, such roots may
 //   be missed by the algorithm because their imaginary part is large.  
 // Arguments:
@@ -1122,11 +1221,13 @@ function _poly_roots(p, pderiv, s, z, tol, i=0) =
 //   eps = used to determine whether imaginary parts of roots are zero
 //   tol = tolerance for the complex polynomial root finder
 
-function real_roots(p,eps=EPSILON,tol=1e-14) =
+function real_roots(p,eps=undef,tol=1e-14) =
    let( 
-       roots = poly_roots(p)
+       roots_err = poly_roots(p,error_bound=true),
+       roots = roots_err[0],
+       err = roots_err[1]
    )
-   [for(z=roots) if (abs(z.y)/(1+norm(z))<eps) z.x];
-
+   is_def(eps) ? [for(z=roots) if (abs(z.y)/(1+norm(z))<eps) z.x]
+               : [for(i=idx(roots)) if (abs(roots[i].y)<=err[i]) roots[i].x];
 
 // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap

From f190b50c4dc81fb1a77a45d8b3dd6c16a2a49855 Mon Sep 17 00:00:00 2001
From: Adrian Mariano <avm4@cornell.edu>
Date: Sun, 12 Jul 2020 07:10:56 -0400
Subject: [PATCH 3/3] substituted correct file

---
 tests/test_math.scad | 1192 +++++++++++++++++++++++++-----------------
 1 file changed, 724 insertions(+), 468 deletions(-)

diff --git a/tests/test_math.scad b/tests/test_math.scad
index dcbad7d..8b34da7 100644
--- a/tests/test_math.scad
+++ b/tests/test_math.scad
@@ -1,650 +1,906 @@
-include <BOSL/constants.scad>
-include <BOSL/math.scad>
-
-eps = 1e-9;
+include <BOSL2/std.scad>
 
 // Simple Calculations
 
 module test_quant() {
-	assert(quant(-4,3) == -3);
-	assert(quant(-3,3) == -3);
-	assert(quant(-2,3) == -3);
-	assert(quant(-1,3) == 0);
-	assert(quant(0,3) == 0);
-	assert(quant(1,3) == 0);
-	assert(quant(2,3) == 3);
-	assert(quant(3,3) == 3);
-	assert(quant(4,3) == 3);
-	assert(quant(7,3) == 6);
+    assert_equal(quant(-4,3), -3);
+    assert_equal(quant(-3,3), -3);
+    assert_equal(quant(-2,3), -3);
+    assert_equal(quant(-1,3), 0);
+    assert_equal(quant(0,3), 0);
+    assert_equal(quant(1,3), 0);
+    assert_equal(quant(2,3), 3);
+    assert_equal(quant(3,3), 3);
+    assert_equal(quant(4,3), 3);
+    assert_equal(quant(7,3), 6);
+    assert_equal(quant([12,13,13.1,14,14.1,15,16],4), [12,12,12,16,16,16,16]);
+    assert_equal(quant([9,10,10.4,10.5,11,12],3), [9,9,9,12,12,12]);
+    assert_equal(quant([[9,10,10.4],[10.5,11,12]],3), [[9,9,9],[12,12,12]]);
 }
 test_quant();
 
 
 module test_quantdn() {
-	assert(quantdn(-4,3) == -6);
-	assert(quantdn(-3,3) == -3);
-	assert(quantdn(-2,3) == -3);
-	assert(quantdn(-1,3) == -3);
-	assert(quantdn(0,3) == 0);
-	assert(quantdn(1,3) == 0);
-	assert(quantdn(2,3) == 0);
-	assert(quantdn(3,3) == 3);
-	assert(quantdn(4,3) == 3);
-	assert(quantdn(7,3) == 6);
+    assert_equal(quantdn(-4,3), -6);
+    assert_equal(quantdn(-3,3), -3);
+    assert_equal(quantdn(-2,3), -3);
+    assert_equal(quantdn(-1,3), -3);
+    assert_equal(quantdn(0,3), 0);
+    assert_equal(quantdn(1,3), 0);
+    assert_equal(quantdn(2,3), 0);
+    assert_equal(quantdn(3,3), 3);
+    assert_equal(quantdn(4,3), 3);
+    assert_equal(quantdn(7,3), 6);
+    assert_equal(quantdn([12,13,13.1,14,14.1,15,16],4), [12,12,12,12,12,12,16]);
+    assert_equal(quantdn([9,10,10.4,10.5,11,12],3), [9,9,9,9,9,12]);
+    assert_equal(quantdn([[9,10,10.4],[10.5,11,12]],3), [[9,9,9],[9,9,12]]);
 }
 test_quantdn();
 
 
 module test_quantup() {
-	assert(quantup(-4,3) == -3);
-	assert(quantup(-3,3) == -3);
-	assert(quantup(-2,3) == 0);
-	assert(quantup(-1,3) == 0);
-	assert(quantup(0,3) == 0);
-	assert(quantup(1,3) == 3);
-	assert(quantup(2,3) == 3);
-	assert(quantup(3,3) == 3);
-	assert(quantup(4,3) == 6);
-	assert(quantup(7,3) == 9);
+    assert_equal(quantup(-4,3), -3);
+    assert_equal(quantup(-3,3), -3);
+    assert_equal(quantup(-2,3), 0);
+    assert_equal(quantup(-1,3), 0);
+    assert_equal(quantup(0,3), 0);
+    assert_equal(quantup(1,3), 3);
+    assert_equal(quantup(2,3), 3);
+    assert_equal(quantup(3,3), 3);
+    assert_equal(quantup(4,3), 6);
+    assert_equal(quantup(7,3), 9);
+    assert_equal(quantup([12,13,13.1,14,14.1,15,16],4), [12,16,16,16,16,16,16]);
+    assert_equal(quantup([9,10,10.4,10.5,11,12],3), [9,12,12,12,12,12]);
+    assert_equal(quantup([[9,10,10.4],[10.5,11,12]],3), [[9,12,12],[12,12,12]]);
 }
 test_quantup();
 
 
 module test_constrain() {
-	assert(constrain(-2,-1,1) == -1);
-	assert(constrain(-1.75,-1,1) == -1);
-	assert(constrain(-1,-1,1) == -1);
-	assert(constrain(-0.75,-1,1) == -0.75);
-	assert(constrain(0,-1,1) == 0);
-	assert(constrain(0.75,-1,1) == 0.75);
-	assert(constrain(1,-1,1) == 1);
-	assert(constrain(1.75,-1,1) == 1);
-	assert(constrain(2,-1,1) == 1);
+    assert_equal(constrain(-2,-1,1), -1);
+    assert_equal(constrain(-1.75,-1,1), -1);
+    assert_equal(constrain(-1,-1,1), -1);
+    assert_equal(constrain(-0.75,-1,1), -0.75);
+    assert_equal(constrain(0,-1,1), 0);
+    assert_equal(constrain(0.75,-1,1), 0.75);
+    assert_equal(constrain(1,-1,1), 1);
+    assert_equal(constrain(1.75,-1,1), 1);
+    assert_equal(constrain(2,-1,1), 1);
 }
 test_constrain();
 
 
+module test_is_matrix() {
+    assert(is_matrix([[2,3,4],[5,6,7],[8,9,10]]));
+    assert(is_matrix([[2,3,4],[5,6,7],[8,9,10]],square=true));
+    assert(is_matrix([[2,3,4],[5,6,7],[8,9,10]],square=false));
+    assert(is_matrix([[2,3],[5,6],[8,9]],m=3,n=2));
+    assert(is_matrix([[2,3,4],[5,6,7]],m=2,n=3));
+    assert(is_matrix([[2,3,4],[5,6,7]],2,3));
+    assert(is_matrix([[2,3,4],[5,6,7]],m=2));
+    assert(is_matrix([[2,3,4],[5,6,7]],2));
+    assert(is_matrix([[2,3,4],[5,6,7]],n=3));
+    assert(!is_matrix([[2,3,4],[5,6,7]],m=4));
+    assert(!is_matrix([[2,3,4],[5,6,7]],n=5));
+    assert(!is_matrix([[2,3,4],[5,6,7]],m=2,n=3,square=true));
+    assert(is_matrix([[2,3,4],[5,6,7],[8,9,10]],square=false));
+    assert(!is_matrix([[2,3],[5,6],[8,9]],m=2,n=3));
+    assert(!is_matrix([[2,3,4],[5,6,7]],m=3,n=2));
+    assert(!is_matrix(undef));
+    assert(!is_matrix(NAN));
+    assert(!is_matrix(INF));
+    assert(!is_matrix(-5));
+    assert(!is_matrix(0));
+    assert(!is_matrix(5));
+    assert(!is_matrix(""));
+    assert(!is_matrix("foo"));
+    assert(!is_matrix([3,4,5]));
+    assert(!is_matrix([]));
+}
+test_is_matrix();
+
+
+module test_approx() {
+    assert_equal(approx(PI, 3.141592653589793236), true);
+    assert_equal(approx(PI, 3.1415926), false);
+    assert_equal(approx(PI, 3.1415926, eps=1e-6), true);
+    assert_equal(approx(-PI, -3.141592653589793236), true);
+    assert_equal(approx(-PI, -3.1415926), false);
+    assert_equal(approx(-PI, -3.1415926, eps=1e-6), true);
+    assert_equal(approx(1/3, 0.3333333333), true);
+    assert_equal(approx(-1/3, -0.3333333333), true);
+    assert_equal(approx(10*[cos(30),sin(30)], 10*[sqrt(3)/2, 1/2]), true);
+}
+test_approx();
+
+
+module test_min_index() {
+    vals = rands(-100,100,100);
+    minval = min(vals);
+    minidx = min_index(vals);
+    assert_equal(vals[minidx], minval);
+    assert_equal(min_index([3,4,5,6]), 0);
+    assert_equal(min_index([4,3,5,6]), 1);
+    assert_equal(min_index([4,5,3,6]), 2);
+    assert_equal(min_index([4,5,6,3]), 3);
+    assert_equal(min_index([6,5,4,3]), 3);
+    assert_equal(min_index([6,3,4,5]), 1);
+    assert_equal(min_index([-56,72,-874,5]), 2);
+}
+test_min_index();
+
+
+module test_max_index() {
+    vals = rands(-100,100,100);
+    maxval = max(vals);
+    maxidx = max_index(vals);
+    assert_equal(vals[maxidx], maxval);
+    assert_equal(max_index([3,4,5,6]), 3);
+    assert_equal(max_index([3,4,6,5]), 2);
+    assert_equal(max_index([3,6,4,5]), 1);
+    assert_equal(max_index([6,3,4,5]), 0);
+    assert_equal(max_index([5,6,4,3]), 1);
+    assert_equal(max_index([-56,72,-874,5]), 1);
+}
+test_max_index();
+
+
 module test_posmod() {
-	assert(posmod(-5,3) == 1);
-	assert(posmod(-4,3) == 2);
-	assert(posmod(-3,3) == 0);
-	assert(posmod(-2,3) == 1);
-	assert(posmod(-1,3) == 2);
-	assert(posmod(0,3) == 0);
-	assert(posmod(1,3) == 1);
-	assert(posmod(2,3) == 2);
-	assert(posmod(3,3) == 0);
+    assert_equal(posmod(-5,3), 1);
+    assert_equal(posmod(-4,3), 2);
+    assert_equal(posmod(-3,3), 0);
+    assert_equal(posmod(-2,3), 1);
+    assert_equal(posmod(-1,3), 2);
+    assert_equal(posmod(0,3), 0);
+    assert_equal(posmod(1,3), 1);
+    assert_equal(posmod(2,3), 2);
+    assert_equal(posmod(3,3), 0);
 }
 test_posmod();
 
 
+module test_modang() {
+    assert_equal(modang(-700), 20);
+    assert_equal(modang(-270), 90);
+    assert_equal(modang(-120), -120);
+    assert_equal(modang(120), 120);
+    assert_equal(modang(270), -90);
+    assert_equal(modang(700), -20);
+}
+test_modang();
+
+
 module test_modrange() {
-	assert(modrange(-5,5,3) == [1,2]);
-	assert(modrange(-1,4,3) == [2,0,1]);
-	assert(modrange(1,8,10,step=2) == [1,3,5,7]);
-	assert(modrange(5,12,10,step=2) == [5,7,9,1]);
+    assert_equal(modrange(-5,5,3), [1,2]);
+    assert_equal(modrange(-1,4,3), [2,0,1]);
+    assert_equal(modrange(1,8,10,step=2), [1,3,5,7]);
+    assert_equal(modrange(5,12,10,step=2), [5,7,9,1]);
 }
 test_modrange();
 
 
+module test_sqr() {
+    assert_equal(sqr(-3), 9);
+    assert_equal(sqr(0), 0);
+    assert_equal(sqr(1), 1);
+    assert_equal(sqr(2), 4);
+    assert_equal(sqr(2.5), 6.25);
+    assert_equal(sqr(3), 9);
+    assert_equal(sqr(16), 256);
+    assert_equal(sqr([2,3,4]), [4,9,16]);
+    assert_equal(sqr([[2,3,4],[3,5,7]]), [[4,9,16],[9,25,49]]);
+    assert_equal(sqr([]),[]);
+}
+test_sqr();
+
+
+module test_log2() {
+    assert_equal(log2(0.125), -3);
+    assert_equal(log2(16), 4);
+    assert_equal(log2(256), 8);
+}
+test_log2();
+
+
+module test_rand_int() {
+    nums = rand_int(-100,100,1000,seed=2134);
+    assert_equal(len(nums), 1000);
+    for (num = nums) {
+        assert(num>=-100);
+        assert(num<=100);
+        assert_equal(num, floor(num));
+    }
+}
+test_rand_int();
+
+
+module test_gaussian_rands() {
+    nums1 = gaussian_rands(0,10,1000,seed=2132);
+    nums2 = gaussian_rands(0,10,1000,seed=2130);
+    nums3 = gaussian_rands(0,10,1000,seed=2132);
+    assert_equal(len(nums1), 1000);
+    assert_equal(len(nums2), 1000);
+    assert_equal(len(nums3), 1000);
+    assert_equal(nums1, nums3);
+    assert(nums1!=nums2);
+}
+test_gaussian_rands();
+
+
+module test_log_rands() {
+    nums1 = log_rands(0,100,10,1000,seed=2189);
+    nums2 = log_rands(0,100,10,1000,seed=2310);
+    nums3 = log_rands(0,100,10,1000,seed=2189);
+    assert_equal(len(nums1), 1000);
+    assert_equal(len(nums2), 1000);
+    assert_equal(len(nums3), 1000);
+    assert_equal(nums1, nums3);
+    assert(nums1!=nums2);
+}
+test_log_rands();
+
+
 module test_segs() {
-	assert(segs(50,$fn=8) == 8);
-	assert(segs(50,$fa=2,$fs=2) == 158);
+    assert_equal(segs(50,$fn=8), 8);
+    assert_equal(segs(50,$fa=2,$fs=2), 158);
 }
 test_segs();
 
 
 module test_lerp() {
-	assert(lerp(-20,20,0) == -20);
-	assert(lerp(-20,20,0.25) == -10);
-	assert(lerp(-20,20,0.5) == 0);
-	assert(lerp(-20,20,0.75) == 10);
-	assert(lerp(-20,20,1) == 20);
-	assert(lerp([10,10],[30,-10],0.5) == [20,0]);
+    assert_equal(lerp(-20,20,0), -20);
+    assert_equal(lerp(-20,20,0.25), -10);
+    assert_equal(lerp(-20,20,0.5), 0);
+    assert_equal(lerp(-20,20,0.75), 10);
+    assert_equal(lerp(-20,20,1), 20);
+    assert_equal(lerp(-20,20,[0,0.25,0.5,0.75,1]), [-20,-10,0,10,20]);
+    assert_equal(lerp(-20,20,[0:0.25:1]), [-20,-10,0,10,20]);
+    assert_equal(lerp([10,10],[30,-10],0.5), [20,0]);
 }
 test_lerp();
 
 
 module test_hypot() {
-	assert(hypot(20,30) == norm([20,30]));
+    assert_approx(hypot(20,30), norm([20,30]));
 }
 test_hypot();
 
 
 module test_sinh() {
-	assert(abs(sinh(-2)+3.6268604078) < eps);
-	assert(abs(sinh(-1)+1.1752011936) < eps);
-	assert(abs(sinh(0)) < eps);
-	assert(abs(sinh(1)-1.1752011936) < eps);
-	assert(abs(sinh(2)-3.6268604078) < eps);
+    assert_approx(sinh(-2), -3.6268604078);
+    assert_approx(sinh(-1), -1.1752011936);
+    assert_approx(sinh(0), 0);
+    assert_approx(sinh(1), 1.1752011936);
+    assert_approx(sinh(2), 3.6268604078);
 }
 test_sinh();
 
 
 module test_cosh() {
-	assert(abs(cosh(-2)-3.7621956911) < eps);
-	assert(abs(cosh(-1)-1.5430806348) < eps);
-	assert(abs(cosh(0)-1) < eps);
-	assert(abs(cosh(1)-1.5430806348) < eps);
-	assert(abs(cosh(2)-3.7621956911) < eps);
+    assert_approx(cosh(-2), 3.7621956911);
+    assert_approx(cosh(-1), 1.5430806348);
+    assert_approx(cosh(0), 1);
+    assert_approx(cosh(1), 1.5430806348);
+    assert_approx(cosh(2), 3.7621956911);
 }
 test_cosh();
 
 
 module test_tanh() {
-	assert(abs(tanh(-2)+0.9640275801) < eps);
-	assert(abs(tanh(-1)+0.761594156) < eps);
-	assert(abs(tanh(0)) < eps);
-	assert(abs(tanh(1)-0.761594156) < eps);
-	assert(abs(tanh(2)-0.9640275801) < eps);
+    assert_approx(tanh(-2), -0.9640275801);
+    assert_approx(tanh(-1), -0.761594156);
+    assert_approx(tanh(0), 0);
+    assert_approx(tanh(1), 0.761594156);
+    assert_approx(tanh(2), 0.9640275801);
 }
 test_tanh();
 
 
 module test_asinh() {
-	assert(abs(asinh(sinh(-2))+2) < eps);
-	assert(abs(asinh(sinh(-1))+1) < eps);
-	assert(abs(asinh(sinh(0))) < eps);
-	assert(abs(asinh(sinh(1))-1) < eps);
-	assert(abs(asinh(sinh(2))-2) < eps);
+    assert_approx(asinh(sinh(-2)), -2);
+    assert_approx(asinh(sinh(-1)), -1);
+    assert_approx(asinh(sinh(0)), 0);
+    assert_approx(asinh(sinh(1)), 1);
+    assert_approx(asinh(sinh(2)), 2);
 }
 test_asinh();
 
 
 module test_acosh() {
-	assert(abs(acosh(cosh(-2))-2) < eps);
-	assert(abs(acosh(cosh(-1))-1) < eps);
-	assert(abs(acosh(cosh(0))) < eps);
-	assert(abs(acosh(cosh(1))-1) < eps);
-	assert(abs(acosh(cosh(2))-2) < eps);
+    assert_approx(acosh(cosh(-2)), 2);
+    assert_approx(acosh(cosh(-1)), 1);
+    assert_approx(acosh(cosh(0)), 0);
+    assert_approx(acosh(cosh(1)), 1);
+    assert_approx(acosh(cosh(2)), 2);
 }
 test_acosh();
 
 
 module test_atanh() {
-	assert(abs(atanh(tanh(-2))+2) < eps);
-	assert(abs(atanh(tanh(-1))+1) < eps);
-	assert(abs(atanh(tanh(0))) < eps);
-	assert(abs(atanh(tanh(1))-1) < eps);
-	assert(abs(atanh(tanh(2))-2) < eps);
+    assert_approx(atanh(tanh(-2)), -2);
+    assert_approx(atanh(tanh(-1)), -1);
+    assert_approx(atanh(tanh(0)), 0);
+    assert_approx(atanh(tanh(1)), 1);
+    assert_approx(atanh(tanh(2)), 2);
 }
 test_atanh();
 
 
 module test_sum() {
-	assert(sum([1,2,3]) == 6);
-	assert(sum([-2,-1,0,1,2]) == 0);
-	assert(sum([[1,2,3], [3,4,5], [5,6,7]]) == [9,12,15]);
+    assert_equal(sum([]), 0);
+    assert_equal(sum([],dflt=undef), undef);
+    assert_equal(sum([1,2,3]), 6);
+    assert_equal(sum([-2,-1,0,1,2]), 0);
+    assert_equal(sum([[1,2,3], [3,4,5], [5,6,7]]), [9,12,15]);
 }
 test_sum();
 
 
+module test_cumsum() {
+    assert_equal(cumsum([]), []);
+    assert_equal(cumsum([1,1,1]), [1,2,3]);
+    assert_equal(cumsum([2,2,2]), [2,4,6]);
+    assert_equal(cumsum([1,2,3]), [1,3,6]);
+    assert_equal(cumsum([-2,-1,0,1,2]), [-2,-3,-3,-2,0]);
+    assert_equal(cumsum([[1,2,3], [3,4,5], [5,6,7]]), [[1,2,3],[4,6,8],[9,12,15]]);
+}
+test_cumsum();
+
+
 module test_sum_of_squares() {
-	assert(sum_of_squares([1,2,3]) == 14);
-	assert(sum_of_squares([1,2,4]) == 21);
-	assert(sum_of_squares([-3,-2,-1]) == 14);
+    assert_equal(sum_of_squares([1,2,3]), 14);
+    assert_equal(sum_of_squares([1,2,4]), 21);
+    assert_equal(sum_of_squares([-3,-2,-1]), 14);
 }
 test_sum_of_squares();
 
 
 module test_sum_of_sines() {
-	assert(sum_of_sines(0, [[3,4,0],[2,2,0]]) == 0);
-	assert(sum_of_sines(45, [[3,4,0],[2,2,0]]) == 2);
-	assert(sum_of_sines(90, [[3,4,0],[2,2,0]]) == 0);
-	assert(sum_of_sines(135, [[3,4,0],[2,2,0]]) == -2);
-	assert(sum_of_sines(180, [[3,4,0],[2,2,0]]) == 0);
+    assert_equal(sum_of_sines(0, [[3,4,0],[2,2,0]]), 0);
+    assert_equal(sum_of_sines(45, [[3,4,0],[2,2,0]]), 2);
+    assert_equal(sum_of_sines(90, [[3,4,0],[2,2,0]]), 0);
+    assert_equal(sum_of_sines(135, [[3,4,0],[2,2,0]]), -2);
+    assert_equal(sum_of_sines(180, [[3,4,0],[2,2,0]]), 0);
 }
 test_sum_of_sines();
 
 
+module test_deltas() {
+    assert_equal(deltas([2,5,9,17]), [3,4,8]);
+    assert_equal(deltas([[1,2,3], [3,6,8], [4,8,11]]), [[2,4,5], [1,2,3]]);
+}
+test_deltas();
+
+
+module test_product() {
+    assert_equal(product([2,3,4]), 24);
+    assert_equal(product([[1,2,3], [3,4,5], [5,6,7]]), [15, 48, 105]);
+    m1 = [[2,3,4],[4,5,6],[6,7,8]];
+    m2 = [[4,1,2],[3,7,2],[8,7,4]];
+    m3 = [[3,7,8],[9,2,4],[5,8,3]];
+    assert_equal(product([m1,m2,m3]), m1*m2*m3);
+}
+test_product();
+
+
 module test_mean() {
-	assert(mean([2,3,4]) == 3);
-	assert(mean([[1,2,3], [3,4,5], [5,6,7]]) == [3,4,5]);
+    assert_equal(mean([2,3,4]), 3);
+    assert_equal(mean([[1,2,3], [3,4,5], [5,6,7]]), [3,4,5]);
 }
 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]);
+}
+test_median();
+
+
+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();
+
+
 // Logic
 
 
 module test_compare_vals() {
-	assert(compare_vals(-10,0) == -1);
-	assert(compare_vals(10,0) == 1);
-	assert(compare_vals(10,10) == 0);
+    assert(compare_vals(-10,0) < 0);
+    assert(compare_vals(10,0) > 0);
+    assert(compare_vals(10,10) == 0);
 
-	assert(compare_vals("abc","abcd") == -1);
-	assert(compare_vals("abcd","abc") == 1);
-	assert(compare_vals("abcd","abcd") == 0);
+    assert(compare_vals("abc","abcd") < 0);
+    assert(compare_vals("abcd","abc") > 0);
+    assert(compare_vals("abcd","abcd") == 0);
 
-	assert(compare_vals(false,false) == 0);
-	assert(compare_vals(true,false) == 1);
-	assert(compare_vals(false,true) == -1);
-	assert(compare_vals(true,true) == 0);
+    assert(compare_vals(false,false) == 0);
+    assert(compare_vals(true,false) > 0);
+    assert(compare_vals(false,true) < 0);
+    assert(compare_vals(true,true) == 0);
 
-	assert(compare_vals([2,3,4], [2,3,4,5]) == -1);
-	assert(compare_vals([2,3,4,5], [2,3,4,5]) == 0);
-	assert(compare_vals([2,3,4,5], [2,3,4]) == 1);
-	assert(compare_vals([2,3,4,5], [2,3,5,5]) == -1);
-	assert(compare_vals([[2,3,4,5]], [[2,3,5,5]]) == -1);
+    assert(compare_vals([2,3,4], [2,3,4,5]) < 0);
+    assert(compare_vals([2,3,4,5], [2,3,4,5]) == 0);
+    assert(compare_vals([2,3,4,5], [2,3,4]) > 0);
+    assert(compare_vals([2,3,4,5], [2,3,5,5]) < 0);
+    assert(compare_vals([[2,3,4,5]], [[2,3,5,5]]) < 0);
 
-	assert(compare_vals([[2,3,4],[3,4,5]], [[2,3,4], [3,4,5]]) == 0);
-	assert(compare_vals([[2,3,4],[3,4,5]], [[2,3,4,5], [3,4,5]]) == -1);
-	assert(compare_vals([[2,3,4],[3,4,5]], [[2,3,4], [3,4,5,6]]) == -1);
-	assert(compare_vals([[2,3,4,5],[3,4,5]], [[2,3,4], [3,4,5]]) == 1);
-	assert(compare_vals([[2,3,4],[3,4,5,6]], [[2,3,4], [3,4,5]]) == 1);
-	assert(compare_vals([[2,3,4],[3,5,5]], [[2,3,4], [3,4,5]]) == 1);
-	assert(compare_vals([[2,3,4],[3,4,5]], [[2,3,4], [3,5,5]]) == -1);
+    assert(compare_vals([[2,3,4],[3,4,5]], [[2,3,4], [3,4,5]]) == 0);
+    assert(compare_vals([[2,3,4],[3,4,5]], [[2,3,4,5], [3,4,5]]) < 0);
+    assert(compare_vals([[2,3,4],[3,4,5]], [[2,3,4], [3,4,5,6]]) < 0);
+    assert(compare_vals([[2,3,4,5],[3,4,5]], [[2,3,4], [3,4,5]]) > 0);
+    assert(compare_vals([[2,3,4],[3,4,5,6]], [[2,3,4], [3,4,5]]) > 0);
+    assert(compare_vals([[2,3,4],[3,5,5]], [[2,3,4], [3,4,5]]) > 0);
+    assert(compare_vals([[2,3,4],[3,4,5]], [[2,3,4], [3,5,5]]) < 0);
+
+    assert(compare_vals(undef, undef) == 0);
+    assert(compare_vals(undef, true) < 0);
+    assert(compare_vals(undef, 0) < 0);
+    assert(compare_vals(undef, "foo") < 0);
+    assert(compare_vals(undef, [2,3,4]) < 0);
+    assert(compare_vals(undef, [0:3]) < 0);
+
+    assert(compare_vals(true, undef) > 0);
+    assert(compare_vals(true, true) == 0);
+    assert(compare_vals(true, 0) < 0);
+    assert(compare_vals(true, "foo") < 0);
+    assert(compare_vals(true, [2,3,4]) < 0);
+    assert(compare_vals(true, [0:3]) < 0);
+
+    assert(compare_vals(0, undef) > 0);
+    assert(compare_vals(0, true) > 0);
+    assert(compare_vals(0, 0) == 0);
+    assert(compare_vals(0, "foo") < 0);
+    assert(compare_vals(0, [2,3,4]) < 0);
+    assert(compare_vals(0, [0:3]) < 0);
+
+    assert(compare_vals(1, undef) > 0);
+    assert(compare_vals(1, true) > 0);
+    assert(compare_vals(1, 1) == 0);
+    assert(compare_vals(1, "foo") < 0);
+    assert(compare_vals(1, [2,3,4]) < 0);
+    assert(compare_vals(1, [0:3]) < 0);
+
+    assert(compare_vals("foo", undef) > 0);
+    assert(compare_vals("foo", true) > 0);
+    assert(compare_vals("foo", 1) > 0);
+    assert(compare_vals("foo", "foo") == 0);
+    assert(compare_vals("foo", [2,3,4]) < 0);
+    assert(compare_vals("foo", [0:3]) < 0);
+
+    assert(compare_vals([2,3,4], undef) > 0);
+    assert(compare_vals([2,3,4], true) > 0);
+    assert(compare_vals([2,3,4], 1) > 0);
+    assert(compare_vals([2,3,4], "foo") > 0);
+    assert(compare_vals([2,3,4], [2,3,4]) == 0);
+    assert(compare_vals([2,3,4], [0:3]) < 0);
+
+    assert(compare_vals([0:3], undef) > 0);
+    assert(compare_vals([0:3], true) > 0);
+    assert(compare_vals([0:3], 1) > 0);
+    assert(compare_vals([0:3], "foo") > 0);
+    assert(compare_vals([0:3], [2,3,4]) > 0);
+    assert(compare_vals([0:3], [0:3]) == 0);
 }
 test_compare_vals();
 
 
 module test_compare_lists() {
-	assert(compare_lists([2,3,4], [2,3,4,5]) == -1);
-	assert(compare_lists([2,3,4,5], [2,3,4,5]) == 0);
-	assert(compare_lists([2,3,4,5], [2,3,4]) == 1);
-	assert(compare_lists([2,3,4,5], [2,3,5,5]) == -1);
+    assert(compare_lists([2,3,4], [2,3,4,5]) < 0);
+    assert(compare_lists([2,3,4,5], [2,3,4,5]) == 0);
+    assert(compare_lists([2,3,4,5], [2,3,4]) > 0);
+    assert(compare_lists([2,3,4,5], [2,3,5,5]) < 0);
 
-	assert(compare_lists([[2,3,4],[3,4,5]], [[2,3,4], [3,4,5]]) == 0);
-	assert(compare_lists([[2,3,4],[3,4,5]], [[2,3,4,5], [3,4,5]]) == -1);
-	assert(compare_lists([[2,3,4],[3,4,5]], [[2,3,4], [3,4,5,6]]) == -1);
-	assert(compare_lists([[2,3,4,5],[3,4,5]], [[2,3,4], [3,4,5]]) == 1);
-	assert(compare_lists([[2,3,4],[3,4,5,6]], [[2,3,4], [3,4,5]]) == 1);
-	assert(compare_lists([[2,3,4],[3,5,5]], [[2,3,4], [3,4,5]]) == 1);
-	assert(compare_lists([[2,3,4],[3,4,5]], [[2,3,4], [3,5,5]]) == -1);
+    assert(compare_lists([[2,3,4],[3,4,5]], [[2,3,4], [3,4,5]]) == 0);
+    assert(compare_lists([[2,3,4],[3,4,5]], [[2,3,4,5], [3,4,5]]) < 0);
+    assert(compare_lists([[2,3,4],[3,4,5]], [[2,3,4], [3,4,5,6]]) < 0);
+    assert(compare_lists([[2,3,4,5],[3,4,5]], [[2,3,4], [3,4,5]]) > 0);
+    assert(compare_lists([[2,3,4],[3,4,5,6]], [[2,3,4], [3,4,5]]) > 0);
+    assert(compare_lists([[2,3,4],[3,5,5]], [[2,3,4], [3,4,5]]) > 0);
+    assert(compare_lists([[2,3,4],[3,4,5]], [[2,3,4], [3,5,5]]) < 0);
 
-	assert(compare_lists("cat", "bat") == 1);
-	assert(compare_lists(["cat"], ["bat"]) == 1);
+    assert(compare_lists("cat", "bat") > 0);
+    assert(compare_lists(["cat"], ["bat"]) > 0);
 }
 test_compare_lists();
 
 
 module test_any() {
-	assert(any([0,false,undef]) == false);
-	assert(any([1,false,undef]) == true);
-	assert(any([1,5,true]) == true);
-	assert(any([[0,0], [0,0]]) == false);
-	assert(any([[0,0], [1,0]]) == true);
+    assert_equal(any([0,false,undef]), false);
+    assert_equal(any([1,false,undef]), true);
+    assert_equal(any([1,5,true]), true);
+    assert_equal(any([[0,0], [0,0]]), false);
+    assert_equal(any([[0,0], [1,0]]), true);
+    assert_equal(any([[false,false],[[false,[false],[[[true]]]],false],[false,false]]), true);
+    assert_equal(any([[false,false],[[false,[false],[[[false]]]],false],[false,false]]), false);
+    assert_equal(any([]), false);
 }
 test_any();
 
 
 module test_all() {
-	assert(all([0,false,undef]) == false);
-	assert(all([1,false,undef]) == false);
-	assert(all([1,5,true]) == true);
-	assert(all([[0,0], [0,0]]) == false);
-	assert(all([[0,0], [1,0]]) == false);
-	assert(all([[1,1], [1,1]]) == true);
+    assert_equal(all([0,false,undef]), false);
+    assert_equal(all([1,false,undef]), false);
+    assert_equal(all([1,5,true]), true);
+    assert_equal(all([[0,0], [0,0]]), false);
+    assert_equal(all([[0,0], [1,0]]), false);
+    assert_equal(all([[1,1], [1,1]]), true);
+    assert_equal(all([[true,true],[[true,[true],[[[true]]]],true],[true,true]]), true);    
+    assert_equal(all([[true,true],[[true,[true],[[[false]]]],true],[true,true]]), false);
+    assert_equal(all([]), true);
 }
 test_all();
 
 
 module test_count_true() {
-	assert(count_true([0,false,undef]) == 0);
-	assert(count_true([1,false,undef]) == 1);
-	assert(count_true([1,5,false]) == 2);
-	assert(count_true([1,5,true]) == 3);
-	assert(count_true([[0,0], [0,0]]) == 0);
-	assert(count_true([[0,0], [1,0]]) == 1);
-	assert(count_true([[1,1], [1,1]]) == 4);
-	assert(count_true([[1,1], [1,1]], nmax=3) == 3);
+    assert_equal(count_true([0,false,undef]), 0);
+    assert_equal(count_true([1,false,undef]), 1);
+    assert_equal(count_true([1,5,false]), 2);
+    assert_equal(count_true([1,5,true]), 3);
+    assert_equal(count_true([[0,0], [0,0]]), 0);
+    assert_equal(count_true([[0,0], [1,0]]), 1);
+    assert_equal(count_true([[1,1], [1,1]]), 4);
+    assert_equal(count_true([[1,1], [1,1]], nmax=3), 3);
 }
 test_count_true();
 
 
-
-// List/Array Ops
-
-module test_cdr() {
-	assert(cdr([]) == []);
-	assert(cdr([88]) == []);
-	assert(cdr([1,2,3]) == [2,3]);
-	assert(cdr(["a","b","c"]) == ["b","c"]);
+module test_factorial() {
+    assert_equal(factorial(0), 1);
+    assert_equal(factorial(1), 1);
+    assert_equal(factorial(2), 2);
+    assert_equal(factorial(3), 6);
+    assert_equal(factorial(4), 24);
+    assert_equal(factorial(5), 120);
+    assert_equal(factorial(6), 720);
+    assert_equal(factorial(7), 5040);
+    assert_equal(factorial(8), 40320);
+    assert_equal(factorial(25,21), 303600);
+    assert_equal(factorial(25,25), 1);
 }
-test_cdr();
+test_factorial();
 
 
-module test_replist() {
-	assert(replist(1, 4) == [1,1,1,1]);
-	assert(replist(8, [2,3]) == [[8,8,8], [8,8,8]]);
-	assert(replist(0, [2,2,3]) == [[[0,0,0],[0,0,0]], [[0,0,0],[0,0,0]]]);
-	assert(replist([1,2,3],3) == [[1,2,3], [1,2,3], [1,2,3]]);
+module test_gcd() {
+    assert_equal(gcd(15,25), 5);
+    assert_equal(gcd(15,27), 3);
+    assert_equal(gcd(270,405), 135);
+    assert_equal(gcd(39, 101),1);
+    assert_equal(gcd(15,-25), 5);
+    assert_equal(gcd(-15,25), 5);
+    assert_equal(gcd(5,0),5);
+    assert_equal(gcd(0,5),5);
 }
-test_replist();
+test_gcd();
 
 
-module test_in_list() {
-	assert(in_list("bar", ["foo", "bar", "baz"]));
-	assert(!in_list("bee", ["foo", "bar", "baz"]));
-	assert(in_list("bar", [[2,"foo"], [4,"bar"], [3,"baz"]], idx=1));
+module test_lcm() {
+    assert_equal(lcm(15,25), 75);
+    assert_equal(lcm(15,27), 135);
+    assert_equal(lcm(270,405), 810);
+    assert_equal(lcm([3,5,15,25,35]),525);
 }
-test_in_list();
+test_lcm();
 
 
-module test_slice() {
-	assert(slice([3,4,5,6,7,8,9], 3, 5) == [6,7]);
-	assert(slice([3,4,5,6,7,8,9], 2, -1) == [5,6,7,8,9]);
-	assert(slice([3,4,5,6,7,8,9], 1, 1) == []);
-	assert(slice([3,4,5,6,7,8,9], 6, -1) == [9]);
-	assert(slice([3,4,5,6,7,8,9], 2, -2) == [5,6,7,8]);
+module test_C_times() {
+    assert_equal(C_times([4,5],[9,-4]), [56,29]);
+    assert_equal(C_times([-7,2],[24,3]), [-174, 27]);
 }
-test_slice();
+test_C_times();
 
 
-module test_select() {
-	l = [3,4,5,6,7,8,9];
-	assert(select(l, 5, 6) == [8,9]);
-	assert(select(l, 5, 8) == [8,9,3,4]);
-	assert(select(l, 5, 2) == [8,9,3,4,5]);
-	assert(select(l, -3, -1) == [7,8,9]);
-	assert(select(l, 3, 3) == [6]);
-	assert(select(l, 4) == 7);
-	assert(select(l, -2) == 8);
-	assert(select(l, [1:3]) == [4,5,6]);
-	assert(select(l, [1,3]) == [4,6]);
+module test_C_div() {
+    assert_equal(C_div([56,29],[9,-4]), [4,5]);
+    assert_equal(C_div([-174,27],[-7,2]), [24,3]);
+}    
+test_C_div();
+
+
+module test_back_substitute(){
+   R = [[12,4,3,2],
+        [0,2,-4,2],
+        [0,0,4,5],
+        [0,0,0,15]];
+   assert_approx(back_substitute(R, [1,2,3,3]), [-0.675, 1.8, 0.5, 0.2]);
+   assert_approx(back_substitute(R, [6, 3, 3.5, 37], transpose=true), [0.5, 0.5, 1, 2]);
+   assert_approx(back_substitute(R, [[38,101],[-6,-16], [31, 71], [45, 105]]), [[1, 4],[2,3],[4,9],[3,7]]);
+   assert_approx(back_substitute(R, [[12,48],[8,22],[11,36],[71,164]],transpose=true), [[1, 4],[2,3],[4,9],[3,7]]);
+   assert_approx(back_substitute([[2]], [4]), [2]);
+   sing1 =[[0,4,3,2],
+         [0,3,-4,2],
+         [0,0,4,5],
+         [0,0,0,15]];
+   sing2 =[[12,4,3,2],
+         [0,0,-4,2],
+         [0,0,4,5],
+         [0,0,0,15]];
+   sing3 = [[12,4,3,2],
+        [0,2,-4,2],
+        [0,0,4,5],
+        [0,0,0,0]];
+   assert_approx(back_substitute(sing1, [1,2,3,4]), []);
+   assert_approx(back_substitute(sing2, [1,2,3,4]), []);
+   assert_approx(back_substitute(sing3, [1,2,3,4]), []);
 }
-test_select();
+test_back_substitute();
 
 
-module test_reverse() {
-	assert(reverse([3,4,5,6]) == [6,5,4,3]);
+
+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(subindex(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_reverse();
+test_linear_solve();
 
 
-module test_array_subindex() {
-	v = [[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]];
-	assert(array_subindex(v,2) == [3, 7, 11, 15]);
-	assert(array_subindex(v,[2,1]) == [[3, 2], [7, 6], [11, 10], [15, 14]]);
-	assert(array_subindex(v,[1:3]) == [[2, 3, 4], [6, 7, 8], [10, 11, 12], [14, 15, 16]]);
+module test_outer_product(){
+  assert_equal(outer_product([1,2,3],[4,5,6]), [[4,5,6],[8,10,12],[12,15,18]]);
+  assert_equal(outer_product([9],[7]), [[63]]);
 }
-test_array_subindex();
+test_outer_product();
 
 
-module test_list_range() {
-	assert(list_range(4) == [0,1,2,3]);
-	assert(list_range(n=4, step=2) == [0,2,4,6]);
-	assert(list_range(n=4, s=3, step=3) == [3,6,9,12]);
-	assert(list_range(n=4, s=3, e=9, step=3) == [3,6,9]);
-	assert(list_range(e=3) == [0,1,2,3]);
-	assert(list_range(e=6, step=2) == [0,2,4,6]);
-	assert(list_range(s=3, e=5) == [3,4,5]);
-	assert(list_range(s=3, e=8, step=2) == [3,5,7]);
-	assert(list_range(s=4, e=8, step=2) == [4,6,8]);
-	assert(list_range(n=4, s=[3,4], step=[2,3]) == [[3,4], [5,7], [7,10], [9,13]]);
+module test_deriv(){
+  pent = [for(x=[0:70:359]) [cos(x), sin(x)]];
+  assert_approx(deriv(pent,closed=true), 
+        [[-0.321393804843,0.556670399226],
+         [-0.883022221559,0.321393804843],
+         [-0.604022773555,-0.719846310393],
+         [0.469846310393,-0.813797681349],
+         [0.925416578398,0.163175911167],
+         [0.413175911167,0.492403876506]]);
+  assert_approx(deriv(pent,closed=true,h=2), 
+     0.5*[[-0.321393804843,0.556670399226],
+         [-0.883022221559,0.321393804843],
+         [-0.604022773555,-0.719846310393],
+         [0.469846310393,-0.813797681349],
+         [0.925416578398,0.163175911167],
+         [0.413175911167,0.492403876506]]);
+  assert_approx(deriv(pent,closed=false),
+        [[-0.432937491789,1.55799143673],
+         [-0.883022221559,0.321393804843],
+         [-0.604022773555,-0.719846310393],
+         [0.469846310393,-0.813797681349],
+         [0.925416578398,0.163175911167],
+         [0.696902572292,1.45914323952]]);
+  spent = yscale(8,pent);
+  lens = path_segment_lengths(spent,closed=true);
+  assert_approx(deriv(spent, closed=true, h=lens),
+         [[-0.0381285841663,0.998065839726],
+          [-0.254979378104,0.0449763331253],
+          [-0.216850793938,-0.953089506601],
+          [0.123993253223,-0.982919228715],
+          [0.191478335034,0.0131898128456],
+          [0.0674850818111,0.996109041561]]);
+  assert_approx(deriv(spent, closed=false, h=select(lens,0,-2)),
+         [[-0.0871925973657,0.996191473044],
+          [-0.254979378104,0.0449763331253],
+          [-0.216850793938,-0.953089506601],
+          [0.123993253223,-0.982919228715],
+          [0.191478335034,0.0131898128456],
+          [0.124034734589,0.992277876714]]);
 }
-test_list_range();
+test_deriv();
 
 
-module test_array_shortest() {
-	assert(array_shortest(["foobar", "bazquxx", "abcd"]) == 4);
+module test_deriv2(){
+    oct = [for(x=[0:45:359]) [cos(x), sin(x)]];
+    assert_approx(deriv2(oct),
+           [[-0.828427124746,0.0719095841794],[-0.414213562373,-0.414213562373],[0,-0.585786437627],
+            [0.414213562373,-0.414213562373],[0.585786437627,0],[0.414213562373,0.414213562373],
+            [0,0.585786437627],[-0.636634192232,0.534938683021]]);
+    assert_approx(deriv2(oct,closed=false),
+           [[-0.828427124746,0.0719095841794],[-0.414213562373,-0.414213562373],[0,-0.585786437627],
+            [0.414213562373,-0.414213562373],[0.585786437627,0],[0.414213562373,0.414213562373],
+            [0,0.585786437627],[-0.636634192232,0.534938683021]]);
+    assert_approx(deriv2(oct,closed=true),
+           [[-0.585786437627,0],[-0.414213562373,-0.414213562373],[0,-0.585786437627],
+            [0.414213562373,-0.414213562373],[0.585786437627,0],[0.414213562373,0.414213562373],
+            [0,0.585786437627],[-0.414213562373,0.414213562373]]);
+    assert_approx(deriv2(oct,closed=false,h=2),
+         0.25*[[-0.828427124746,0.0719095841794],[-0.414213562373,-0.414213562373],[0,-0.585786437627],
+            [0.414213562373,-0.414213562373],[0.585786437627,0],[0.414213562373,0.414213562373],
+            [0,0.585786437627],[-0.636634192232,0.534938683021]]);
+    assert_approx(deriv2(oct,closed=true,h=2),
+         0.25* [[-0.585786437627,0],[-0.414213562373,-0.414213562373],[0,-0.585786437627],
+            [0.414213562373,-0.414213562373],[0.585786437627,0],[0.414213562373,0.414213562373],
+            [0,0.585786437627],[-0.414213562373,0.414213562373]]);
 }
-test_array_shortest();
+test_deriv2();
 
 
-module test_array_longest() {
-	assert(array_longest(["foobar", "bazquxx", "abcd"]) == 7);
+module test_deriv3(){
+    oct = [for(x=[0:45:359]) [cos(x), sin(x)]];
+    assert_approx(deriv3(oct),
+           [[0.414213562373,-0.686291501015],[0.414213562373,-0.343145750508],[0.414213562373,0],
+            [0.292893218813,0.292893218813],[0,0.414213562373],[-0.292893218813,0.292893218813],
+            [-0.535533905933,0.0502525316942],[-0.778174593052,-0.192388155425]]);
+    assert_approx(deriv3(oct,closed=false),
+           [[0.414213562373,-0.686291501015],[0.414213562373,-0.343145750508],[0.414213562373,0],
+            [0.292893218813,0.292893218813],[0,0.414213562373],[-0.292893218813,0.292893218813],
+            [-0.535533905933,0.0502525316942],[-0.778174593052,-0.192388155425]]);
+    assert_approx(deriv3(oct,closed=false,h=2),
+           [[0.414213562373,-0.686291501015],[0.414213562373,-0.343145750508],[0.414213562373,0],
+            [0.292893218813,0.292893218813],[0,0.414213562373],[-0.292893218813,0.292893218813],
+            [-0.535533905933,0.0502525316942],[-0.778174593052,-0.192388155425]]/8);
+    assert_approx(deriv3(oct,closed=true),
+           [[0,-0.414213562373],[0.292893218813,-0.292893218813],[0.414213562373,0],[0.292893218813,0.292893218813],
+            [0,0.414213562373],[-0.292893218813,0.292893218813],[-0.414213562373,0],[-0.292893218813,-0.292893218813]]);
+    assert_approx(deriv3(oct,closed=true,h=2),
+           [[0,-0.414213562373],[0.292893218813,-0.292893218813],[0.414213562373,0],[0.292893218813,0.292893218813],
+            [0,0.414213562373],[-0.292893218813,0.292893218813],[-0.414213562373,0],[-0.292893218813,-0.292893218813]]/8);
 }
-test_array_longest();
+test_deriv3();
+  
 
 
-module test_array_pad() {
-	assert(array_pad([4,5,6], 5, 8) == [4,5,6,8,8]);
-	assert(array_pad([4,5,6,7,8], 5, 8) == [4,5,6,7,8]);
-	assert(array_pad([4,5,6,7,8,9], 5, 8) == [4,5,6,7,8,9]);
+module test_polynomial(){
+  assert_equal(polynomial([],12),0);
+  assert_equal(polynomial([],[12,4]),[0,0]);
+  assert_equal(polynomial([1,2,3,4],3),58);
+  assert_equal(polynomial([1,2,3,4],[3,-1]),[47,-41]);
+  assert_equal(polynomial([0,0,2],4),2);
 }
-test_array_pad();
+test_polynomial();
 
 
-module test_array_trim() {
-	assert(array_trim([4,5,6], 5) == [4,5,6]);
-	assert(array_trim([4,5,6,7,8], 5) == [4,5,6,7,8]);
-	assert(array_trim([3,4,5,6,7,8,9], 5) == [3,4,5,6,7]);
+module test_poly_roots(){
+   // Fifth roots of unity
+   assert_approx(
+        poly_roots([1,0,0,0,0,-1]),
+        [[1,0],[0.309016994375,0.951056516295],[-0.809016994375,0.587785252292],
+         [-0.809016994375,-0.587785252292],[0.309016994375,-0.951056516295]]);
+   assert_approx(poly_roots(poly_mult([[1,-2,5],[12,-24,24],[-2, -12, -20],[1,-10,50]])),
+               [[1, 1], [5, 5], [1, 2], [-3, 1], [-3, -1], [1, -1], [1, -2], [5, -5]]);
+   assert_approx(poly_roots([.124,.231,.942, -.334]),
+                 [[0.3242874219074053,0],[-1.093595323856930,2.666477428660098], [-1.093595323856930,-2.666477428660098]]);
 }
-test_array_trim();
+test_poly_roots();
 
-
-module test_array_fit() {
-	assert(array_fit([4,5,6], 5, 8) == [4,5,6,8,8]);
-	assert(array_fit([4,5,6,7,8], 5, 8) == [4,5,6,7,8]);
-	assert(array_fit([3,4,5,6,7,8,9], 5, 8) == [3,4,5,6,7]);
+module test_real_roots(){
+   // Wilkinson polynomial is a nasty test:
+   assert_approx(
+       sort(real_roots(poly_mult([[1,-1],[1,-2],[1,-3],[1,-4],[1,-5],[1,-6],[1,-7],[1,-8],[1,-9],[1,-10]]))),
+       list_range(n=10,s=1));
+   assert_equal(real_roots([3]), []);
+   assert_equal(real_roots(poly_mult([[1,-2,5],[12,-24,24],[-2, -12, -20],[1,-10,50]])),[]);
+   assert_equal(real_roots(poly_mult([[1,-2,5],[12,-24,24],[-2, -12, -20],[1,-10,50],[1,0,0]])),[0,0]);
+   assert_approx(real_roots(poly_mult([[1,-2,5],[12,-24,24],[-2, -12, -20],[1,-10,50],[1,4]])),[-4]);
+   assert(approx(real_roots([1,-10,25]),[5,5],eps=5e-6));
+   assert_approx(real_roots([4,-3]), [0.75]);
+   assert_approx(real_roots([0,0,0,4,-3]), [0.75]);
 }
-test_array_fit();
+test_real_roots();
 
-
-module test_enumerate() {
-	assert(enumerate(["a","b","c"]) == [[0,"a"], [1,"b"], [2,"c"]]);
-	assert(enumerate([[88,"a"],[76,"b"],[21,"c"]], idx=1) == [[0,"a"], [1,"b"], [2,"c"]]);
-	assert(enumerate([["cat","a",12],["dog","b",10],["log","c",14]], idx=[1:2]) == [[0,"a",12], [1,"b",10], [2,"c",14]]);
+// 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_enumerate();
+test_submatrix();
 
 
-module test_array_zip() {
-	v1 = [1,2,3,4];
-	v2 = [5,6,7];
-	v3 = [8,9,10,11];
-	assert(array_zip(v1,v3) == [[1,8],[2,9],[3,10],[4,11]]);
-	assert(array_zip([v1,v3]) == [[1,8],[2,9],[3,10],[4,11]]);
-	assert(array_zip([v1,v2],fit="short") == [[1,5],[2,6],[3,7]]);
-	assert(array_zip([v1,v2],fit="long") == [[1,5],[2,6],[3,7],[4,undef]]);
-	assert(array_zip([v1,v2],fit="long", fill=0) == [[1,5],[2,6],[3,7],[4,0]]);
-	assert(array_zip([v1,v2,v3],fit="long") == [[1,5,8],[2,6,9],[3,7,10],[4,undef,11]]);
+
+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);
+  
+  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]]));
 }
-test_array_zip();
+test_qr_factor();
 
 
-module test_array_group() {
-	v = [1,2,3,4,5,6];
-	assert(array_group(v,2) == [[1,2], [3,4], [5,6]]);
-	assert(array_group(v,3) == [[1,2,3], [4,5,6]]);
-	assert(array_group(v,4,0) == [[1,2,3,4], [5,6,0,0]]);
+module test_poly_mult(){
+  assert_equal(poly_mult([3,2,1],[4,5,6,7]),[12,23,32,38,20,7]);
+  assert_equal(poly_mult([3,2,1],[]),[]);
+  assert_equal(poly_mult([[1,2],[3,4],[5,6]]), [15,68,100,48]);
+  assert_equal(poly_mult([[1,2],[],[5,6]]), []);
+  assert_equal(poly_mult([[3,4,5],[0,0,0]]),[]);
 }
-test_array_group();
+test_poly_mult();
 
-
-module test_flatten() {
-	assert(flatten([[1,2,3], [4,5,[6,7,8]]]) == [1,2,3,4,5,[6,7,8]]);
+ 
+module test_poly_div(){
+  assert_equal(poly_div(poly_mult([4,3,3,2],[2,1,3]), [2,1,3]),[[4,3,3,2],[]]);
+  assert_equal(poly_div([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]]);
 }
-test_flatten();
+test_poly_div();
 
 
-module test_sort() {
-	assert(sort([7,3,9,4,3,1,8]) == [1,3,3,4,7,8,9]);
-	assert(sort(["cat", "oat", "sat", "bat", "vat", "rat", "pat", "mat", "fat", "hat", "eat"]) == ["bat", "cat", "eat", "fat", "hat", "mat", "oat", "pat", "rat", "sat", "vat"]);
-	assert(sort(enumerate([[2,3,4],[1,2,3],[2,4,3]]),idx=1)==[[1,[1,2,3]], [0,[2,3,4]], [2,[2,4,3]]]);
+module test_poly_add(){
+  assert_equal(poly_add([2,3,4],[3,4,5,6]),[3,6,8,10]);
+  assert_equal(poly_add([1,2,3,4],[-1,-2,3,4]), [6,8]);
+  assert_equal(poly_add([1,2,3],-[1,2,3]),[]);
 }
-test_sort();
+test_poly_add();
 
-
-module test_sortidx() {
-	lst1 = ["d","b","e","c"];
-	assert(sortidx(lst1) == [1,3,0,2]);
-	lst2 = [
-		["foo", 88, [0,0,1], false],
-		["bar", 90, [0,1,0], true],
-		["baz", 89, [1,0,0], false],
-		["qux", 23, [1,1,1], true]
-	];
-	assert(sortidx(lst2, idx=1) == [3,0,2,1]);
-	assert(sortidx(lst2, idx=0) == [1,2,0,3]);
-	assert(sortidx(lst2, idx=[1,3]) == [3,0,2,1]);
-	lst3 = [[-4, 0, 0], [0, 0, -4], [0, -4, 0], [-4, 0, 0], [0, -4, 0], [0, 0, 4], [0, 0, -4], [0, 4, 0], [4, 0, 0], [0, 0, 4], [0, 4, 0], [4, 0, 0]];
-	assert(sortidx(lst3)==[0,3,2,4,1,6,5,9,7,10,8,11]);
-}
-test_sortidx();
-
-
-module test_unique() {
-	assert(unique([]) == []);
-	assert(unique([8]) == [8]);
-	assert(unique([7,3,9,4,3,1,8]) == [1,3,4,7,8,9]);
-}
-test_unique();
-
-
-module test_array_dim() {
-	assert(array_dim([[[1,2,3],[4,5,6]],[[7,8,9],[10,11,12]]]) == [2,2,3]);
-	assert(array_dim([[[1,2,3],[4,5,6]],[[7,8,9],[10,11,12]]], 0) == 2);
-	assert(array_dim([[[1,2,3],[4,5,6]],[[7,8,9],[10,11,12]]], 2) == 3);
-	assert(array_dim([[[1,2,3],[4,5,6]],[[7,8,9]]]) == [2,undef,3]);
-}
-test_array_dim();
-
-
-module test_vmul() {
-	assert(vmul([3,4,5], [8,7,6]) == [24,28,30]);
-	assert(vmul([1,2,3], [4,5,6]) == [4,10,18]);
-}
-test_vmul();
-
-
-module test_vdiv() {
-	assert(vdiv([24,28,30], [8,7,6]) == [3, 4, 5]);
-}
-test_vdiv();
-
-
-module test_vabs() {
-	assert(vabs([2,4,8]) == [2,4,8]);
-	assert(vabs([-2,-4,-8]) == [2,4,8]);
-	assert(vabs([-2,4,8]) == [2,4,8]);
-	assert(vabs([2,-4,8]) == [2,4,8]);
-	assert(vabs([2,4,-8]) == [2,4,8]);
-}
-test_vabs();
-
-
-module test_normalize() {
-	assert(normalize([10,0,0]) == [1,0,0]);
-	assert(normalize([0,10,0]) == [0,1,0]);
-	assert(normalize([0,0,10]) == [0,0,1]);
-	assert(abs(norm(normalize([10,10,10]))-1) < eps);
-	assert(abs(norm(normalize([-10,-10,-10]))-1) < eps);
-	assert(abs(norm(normalize([-10,0,0]))-1) < eps);
-	assert(abs(norm(normalize([0,-10,0]))-1) < eps);
-	assert(abs(norm(normalize([0,0,-10]))-1) < eps);
-}
-test_normalize();
-
-
-module test_vector_angle() {
-	vecs = [[10,0,0], [-10,0,0], [0,10,0], [0,-10,0], [0,0,10], [0,0,-10]];
-	for (a=vecs, b=vecs) {
-		if(a==b) {
-			assert(vector_angle(a,b)==0);
-		} else if(a==-b) {
-			assert(vector_angle(a,b)==180);
-		} else {
-			assert(vector_angle(a,b)==90);
-		}
-	}
-	assert(abs(vector_angle([10,10,0],[10,0,0])-45) < eps);
-}
-test_vector_angle();
-
-
-module test_vector_axis() {
-	assert(norm(vector_axis([10,0,0],[10,10,0]) - [0,0,1]) < eps);
-	assert(norm(vector_axis([10,0,0],[0,10,0]) - [0,0,1]) < eps);
-	assert(norm(vector_axis([0,10,0],[10,0,0]) - [0,0,-1]) < eps);
-	assert(norm(vector_axis([0,0,10],[10,0,0]) - [0,1,0]) < eps);
-	assert(norm(vector_axis([10,0,0],[0,0,10]) - [0,-1,0]) < eps);
-	assert(norm(vector_axis([10,0,10],[0,-10,0]) - [sin(45),0,-sin(45)]) < eps);
-}
-test_vector_axis();
-
-
-module test_point2d() {
-	assert(point2d([1,2,3])==[1,2]);
-	assert(point2d([2,3])==[2,3]);
-	assert(point2d([1])==[1,0]);
-}
-test_point2d();
-
-
-module test_path2d() {
-	assert(path2d([[1], [1,2], [1,2,3], [1,2,3,4], [1,2,3,4,5]])==[[1,0],[1,2],[1,2],[1,2],[1,2]]);
-}
-test_path2d();
-
-
-module test_point3d() {
-	assert(point3d([1,2,3,4,5])==[1,2,3]);
-	assert(point3d([1,2,3,4])==[1,2,3]);
-	assert(point3d([1,2,3])==[1,2,3]);
-	assert(point3d([2,3])==[2,3,0]);
-	assert(point3d([1])==[1,0,0]);
-}
-test_point3d();
-
-
-module test_path3d() {
-	assert(path3d([[1], [1,2], [1,2,3], [1,2,3,4], [1,2,3,4,5]])==[[1,0,0],[1,2,0],[1,2,3],[1,2,3],[1,2,3]]);
-}
-test_path3d();
-
-
-module test_translate_points() {
-	pts = [[0,0,1], [0,1,0], [1,0,0], [0,0,-1], [0,-1,0], [-1,0,0]];
-	assert(translate_points(pts, v=[1,2,3]) == [[1,2,4], [1,3,3], [2,2,3], [1,2,2], [1,1,3], [0,2,3]]);
-	assert(translate_points(pts, v=[-1,-2,-3]) == [[-1,-2,-2], [-1,-1,-3], [0,-2,-3], [-1,-2,-4], [-1,-3,-3], [-2,-2,-3]]);
-}
-test_translate_points();
-
-
-module test_scale_points() {
-	pts = [[0,0,1], [0,1,0], [1,0,0], [0,0,-1], [0,-1,0], [-1,0,0]];
-	assert(scale_points(pts, v=[2,3,4]) == [[0,0,4], [0,3,0], [2,0,0], [0,0,-4], [0,-3,0], [-2,0,0]]);
-	assert(scale_points(pts, v=[-2,-3,-4]) == [[0,0,-4], [0,-3,0], [-2,0,0], [0,0,4], [0,3,0], [2,0,0]]);
-	assert(scale_points(pts, v=[1,1,1]) == [[0,0,1], [0,1,0], [1,0,0], [0,0,-1], [0,-1,0], [-1,0,0]]);
-	assert(scale_points(pts, v=[-1,-1,-1]) == [[0,0,-1], [0,-1,0], [-1,0,0], [0,0,1], [0,1,0], [1,0,0]]);
-}
-test_scale_points();
-
-
-module test_rotate_points2d() {
-	pts = [[0,1], [1,0], [0,-1], [-1,0]];
-	s = sin(45);
-	assert(rotate_points2d(pts,45) == [[-s,s],[s,s],[s,-s],[-s,-s]]);
-	assert(rotate_points2d(pts,90) == [[-1,0],[0,1],[1,0],[0,-1]]);
-	assert(rotate_points2d(pts,90,cp=[1,0]) == [[0,-1],[1,0],[2,-1],[1,-2]]);
-}
-test_rotate_points2d();
-
-
-module test_rotate_points3d() {
-	pts = [[0,0,1], [0,1,0], [1,0,0], [0,0,-1], [0,-1,0], [-1,0,0]];
-	assert(rotate_points3d(pts, [90,0,0]) == [[0,-1,0], [0,0,1], [1,0,0], [0,1,0], [0,0,-1], [-1,0,0]]);
-	assert(rotate_points3d(pts, [0,90,0]) == [[1,0,0], [0,1,0], [0,0,-1], [-1,0,0], [0,-1,0], [0,0,1]]);
-	assert(rotate_points3d(pts, [0,0,90]) == [[0,0,1], [-1,0,0], [0,1,0], [0,0,-1], [1,0,0], [0,-1,0]]);
-	assert(rotate_points3d(pts, [0,0,90],cp=[2,0,0]) == [[2,-2,1], [1,-2,0], [2,-1,0], [2,-2,-1], [3,-2,0], [2,-3,0]]);
-	assert(rotate_points3d(pts, 90, axis=V_UP) == [[0,0,1], [-1,0,0], [0,1,0], [0,0,-1], [1,0,0], [0,-1,0]]);
-	assert(rotate_points3d(pts, 90, axis=V_DOWN) == [[0,0,1], [1,0,0], [0,-1,0], [0,0,-1], [-1,0,0], [0,1,0]]);
-	assert(rotate_points3d(pts, 90, axis=V_RIGHT)  == [[0,-1,0], [0,0,1], [1,0,0], [0,1,0], [0,0,-1], [-1,0,0]]);
-	assert(rotate_points3d(pts, from=V_UP, to=V_BACK) == [[0,1,0], [0,0,-1], [1,0,0], [0,-1,0], [0,0,1], [-1,0,0]]);
-	assert(rotate_points3d(pts, 90, from=V_UP, to=V_BACK), [[0,1,0], [-1,0,0], [0,0,-1], [0,-1,0], [1,0,0], [0,0,1]]);
-	assert(rotate_points3d(pts, from=V_UP, to=V_UP*2) == [[0,0,1], [0,1,0], [1,0,0], [0,0,-1], [0,-1,0], [-1,0,0]]);
-	assert(rotate_points3d(pts, from=V_UP, to=V_DOWN*2) == [[0,0,-1], [0,1,0], [-1,0,0], [0,0,1], [0,-1,0], [1,0,0]]);
-}
-test_rotate_points3d();
-
-
-module test_simplify_path()
-{
-	path = [[-20,10],[-10,0],[-5,0],[0,0],[5,0],[10,0], [10,10]];
-	assert(simplify_path(path) == [[-20,10],[-10,0],[10,0], [10,10]]);
-}
-test_simplify_path();
-
-
-module test_simplify_path_indexed()
-{
-	points = [[-20,10],[-10,0],[-5,0],[0,0],[5,0],[10,0], [10,10]];
-	path = list_range(len(points));
-	assert(simplify_path_indexed(points, path) == [0,1,5,6]);
-}
-test_simplify_path_indexed();
-
-
-// vim: noexpandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
+// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap