From 8a25764744cf260d775b7c81540b8449c6a1382d Mon Sep 17 00:00:00 2001 From: RonaldoCMP Date: Wed, 29 Jul 2020 22:41:02 +0100 Subject: [PATCH] Overall review, input data check, format, new function definitions New function definitions in commom.scad: 1. valid_range; 2. _list_pattern New function definitions in math.scad: 1. binomial; 2. binomial_coefficient; 3. convolve; Code review in math: 1. sum; 2. median; 3. is_matrix; 4. approx; 5. count_true; 6. doc of deriv2; 7. polynomial; 8. poly_mult; 9; poly_div; 10. _poly_trim Code change in test_common: 1. new test_valid_range; 2. test_is_consistent Code change in test_math: 1. test_approx; 2. new test_convolve; 3. new test_binomial; 4. new test_binomial_coefficient; 5. test_outer_product; 6. test_polynomial; 7. test_poly_mult; 8. test_poly_div; 9. test_poly_add; --- common.scad | 20 +- math.scad | 623 +++++++++++++++++++++++++++-------------- tests/test_common.scad | 24 +- tests/test_math.scad | 55 +++- 4 files changed, 497 insertions(+), 225 deletions(-) diff --git a/common.scad b/common.scad index 0d88e52..37ba4af 100644 --- a/common.scad +++ b/common.scad @@ -99,6 +99,16 @@ function is_finite(v) = is_num(0*v); function is_range(x) = !is_list(x) && is_finite(x[0]+x[1]+x[2]) ; +// Function: valid_range() +// Description: +// Returns true if its argument is a valid range (deprecated ranges excluded). +function valid_range(x) = + is_range(x) + && ( x[1]>0 + ? x[0]<=x[2] + : ( x[1]<0 && x[0]>=x[2] ) ); + + // Function: is_list_of() // Usage: // is_list_of(list, pattern) @@ -133,9 +143,15 @@ function is_list_of(list,pattern) = // is_consistent([[3,[3,4,[5]]], [5,[2,9,[9]]]]); // Returns true // is_consistent([[3,[3,4,[5]]], [5,[2,9,9]]]); // Returns false function is_consistent(list) = - is_list(list) && is_list_of(list, list[0]); + is_list(list) && is_list_of(list, _list_pattern(list[0])); +//Internal function +//Creates a list with the same structure of `list` with each of its elements substituted by 0. +function _list_pattern(list) = + is_list(list) + ? [for(entry=list) is_list(entry) ? _list_pattern(entry) : 0] + : 0; // Function: same_shape() @@ -146,7 +162,7 @@ function is_consistent(list) = // Example: // same_shape([3,[4,5]],[7,[3,4]]); // Returns true // same_shape([3,4,5], [7,[3,4]]); // Returns false -function same_shape(a,b) = a*0 == b*0; +function same_shape(a,b) = _list_pattern(a) == b*0; // Section: Handling `undef`s. diff --git a/math.scad b/math.scad index 7629cb1..2262cac 100644 --- a/math.scad +++ b/math.scad @@ -33,7 +33,10 @@ NAN = acos(2); // The value `nan`, useful for comparisons. // sqr([3,4]); // Returns: [9,16] // sqr([[1,2],[3,4]]); // Returns [[1,4],[9,16]] // sqr([[1,2],3]); // Returns [[1,4],9] -function sqr(x) = is_list(x) ? [for(val=x) sqr(val)] : x*x; +function sqr(x) = + is_list(x) ? [for(val=x) sqr(val)] : + is_finite(x) ? x*x : + assert(is_finite(x) || is_vector(x), "Input is not neither a number nor a list of numbers."); // Function: log2() @@ -45,8 +48,11 @@ function sqr(x) = is_list(x) ? [for(val=x) sqr(val)] : x*x; // log2(0.125); // Returns: -3 // log2(16); // Returns: 4 // log2(256); // Returns: 8 -function log2(x) = ln(x)/ln(2); +function log2(x) = + assert( is_finite(x), "Input is not a number.") + ln(x)/ln(2); +// this may return NAN or INF; should it check x>0 ? // Function: hypot() // Usage: @@ -60,7 +66,9 @@ function log2(x) = ln(x)/ln(2); // Example: // l = hypot(3,4); // Returns: 5 // l = hypot(3,4,5); // Returns: ~7.0710678119 -function hypot(x,y,z=0) = norm([x,y,z]); +function hypot(x,y,z=0) = + assert( is_vector([x,y,z]), "Improper number(s).") + norm([x,y,z]); // Function: factorial() @@ -76,11 +84,53 @@ function hypot(x,y,z=0) = norm([x,y,z]); // y = factorial(6); // Returns: 720 // z = factorial(9); // Returns: 362880 function factorial(n,d=0) = - assert(n>=0 && d>=0, "Factorial is not defined for negative numbers") + assert(is_int(n) && is_int(d) && n>=0 && d>=0, "Factorial is not defined for negative numbers") assert(d<=n, "d cannot be larger than n") product([1,for (i=[n:-1:d+1]) i]); +// Function: binomial() +// Usage: +// x = binomial(n); +// Description: +// Returns the binomial coefficients of the integer `n`. +// Arguments: +// n = The integer to get the binomial coefficients of +// Example: +// x = binomial(3); // Returns: [1,3,3,1] +// y = binomial(4); // Returns: [1,4,6,4,1] +// z = binomial(6); // Returns: [1,6,15,20,15,6,1] +function binomial(n) = + assert( is_int(n) && n>0, "Input is not an integer greater than 0.") + [for( c = 1, i = 0; + i<=n; + c = c*(n-i)/(i+1), i = i+1 + ) c ] ; + + +// Function: binomial_coefficient() +// Usage: +// x = binomial_coefficient(n,k); +// Description: +// Returns the k-th binomial coefficient of the integer `n`. +// Arguments: +// n = The integer to get the binomial coefficient of +// k = The binomial coefficient index +// Example: +// x = binomial_coefficient(3,2); // Returns: 3 +// y = binomial_coefficient(10,6); // Returns: 210 +function binomial_coefficient(n,k) = + assert( is_int(n) && is_int(k), "Some input is not a number.") + k < 0 || k > n ? 0 : + k ==0 || k ==n ? 1 : + let( k = min(k, n-k), + b = [for( c = 1, i = 0; + i<=k; + c = c*(n-i)/(i+1), i = i+1 + ) c] ) + b[len(b)-1]; + + // Function: lerp() // Usage: // x = lerp(a, b, u); @@ -91,8 +141,8 @@ function factorial(n,d=0) = // If `u` is 0.0, then the value of `a` is returned. // If `u` is 1.0, then the value of `b` is returned. // If `u` is a range, or list of numbers, returns a list of interpolated values. -// It is valid to use a `u` value outside the range 0 to 1. The result will be a predicted -// value along the slope formed by `a` and `b`, but not between those two values. +// It is valid to use a `u` value outside the range 0 to 1. The result will be an extrapolation +// along the slope formed by `a` and `b`. // Arguments: // a = First value or vector. // b = Second value or vector. @@ -113,9 +163,9 @@ function factorial(n,d=0) = // rainbow(pts) translate($item) circle(d=3,$fn=8); function lerp(a,b,u) = assert(same_shape(a,b), "Bad or inconsistent inputs to lerp") - is_num(u)? (1-u)*a + u*b : - assert(!is_undef(u)&&!is_bool(u)&&!is_string(u), "Input u to lerp must be a number, vector, or range.") - [for (v = u) lerp(a,b,v)]; + is_finite(u)? (1-u)*a + u*b : + assert(is_finite(u) || is_vector(u) || valid_range(u), "Input u to lerp must be a number, vector, or range.") + [for (v = u) (1-v)*a + v*b ]; @@ -124,40 +174,45 @@ function lerp(a,b,u) = // Function: sinh() // Description: Takes a value `x`, and returns the hyperbolic sine of it. function sinh(x) = + assert(is_finite(x), "The input must be a finite number.") (exp(x)-exp(-x))/2; // Function: cosh() // Description: Takes a value `x`, and returns the hyperbolic cosine of it. function cosh(x) = + assert(is_finite(x), "The input must be a finite number.") (exp(x)+exp(-x))/2; // Function: tanh() // Description: Takes a value `x`, and returns the hyperbolic tangent of it. function tanh(x) = + assert(is_finite(x), "The input must be a finite number.") sinh(x)/cosh(x); // Function: asinh() // Description: Takes a value `x`, and returns the inverse hyperbolic sine of it. function asinh(x) = + assert(is_finite(x), "The input must be a finite number.") ln(x+sqrt(x*x+1)); // Function: acosh() // Description: Takes a value `x`, and returns the inverse hyperbolic cosine of it. function acosh(x) = + assert(is_finite(x), "The input must be a finite number.") ln(x+sqrt(x*x-1)); // Function: atanh() // Description: Takes a value `x`, and returns the inverse hyperbolic tangent of it. function atanh(x) = + assert(is_finite(x), "The input must be a finite number.") ln((1+x)/(1-x))/2; - // Section: Quantization // Function: quant() @@ -185,8 +240,11 @@ function atanh(x) = // quant([9,10,10.4,10.5,11,12],3); // Returns: [9,9,9,12,12,12] // quant([[9,10,10.4],[10.5,11,12]],3); // Returns: [[9,9,9],[12,12,12]] function quant(x,y) = - is_list(x)? [for (v=x) quant(v,y)] : - floor(x/y+0.5)*y; + assert(is_finite(y) && !approx(y,0,eps=1e-24), "The multiple must be a non zero integer.") + is_list(x) + ? [for (v=x) quant(v,y)] + : assert( is_finite(x), "The input to quantize must be a number or a list of numbers.") + floor(x/y+0.5)*y; // Function: quantdn() @@ -214,8 +272,11 @@ function quant(x,y) = // quantdn([9,10,10.4,10.5,11,12],3); // Returns: [9,9,9,9,9,12] // quantdn([[9,10,10.4],[10.5,11,12]],3); // Returns: [[9,9,9],[9,9,12]] function quantdn(x,y) = - is_list(x)? [for (v=x) quantdn(v,y)] : - floor(x/y)*y; + assert(is_finite(y) && !approx(y,0,eps=1e-24), "The multiple must be a non zero integer.") + is_list(x) + ? [for (v=x) quantdn(v,y)] + : assert( is_finite(x), "The input to quantize must be a number or a list of numbers.") + floor(x/y)*y; // Function: quantup() @@ -243,8 +304,11 @@ function quantdn(x,y) = // quantup([9,10,10.4,10.5,11,12],3); // Returns: [9,12,12,12,12,12] // quantup([[9,10,10.4],[10.5,11,12]],3); // Returns: [[9,12,12],[12,12,12]] function quantup(x,y) = - is_list(x)? [for (v=x) quantup(v,y)] : - ceil(x/y)*y; + assert(is_finite(y) && !approx(y,0,eps=1e-24), "The multiple must be a non zero integer.") + is_list(x) + ? [for (v=x) quantup(v,y)] + : assert( is_finite(x), "The input to quantize must be a number or a list of numbers.") + ceil(x/y)*y; // Section: Constraints and Modulos @@ -264,7 +328,9 @@ function quantup(x,y) = // constrain(0.3, -1, 1); // Returns: 0.3 // constrain(9.1, 0, 9); // Returns: 9 // constrain(-0.1, 0, 9); // Returns: 0 -function constrain(v, minval, maxval) = min(maxval, max(minval, v)); +function constrain(v, minval, maxval) = + assert( is_finite(v+minval+maxval), "Input must be finite number(s).") + min(maxval, max(minval, v)); // Function: posmod() @@ -283,7 +349,9 @@ function constrain(v, minval, maxval) = min(maxval, max(minval, v)); // posmod(270,360); // Returns: 270 // posmod(700,360); // Returns: 340 // posmod(3,2.5); // Returns: 0.5 -function posmod(x,m) = (x%m+m)%m; +function posmod(x,m) = + assert( is_finite(x) && is_finite(m) && !approx(m,0) , "Input must be finite numbers. The divisor cannot be zero.") + (x%m+m)%m; // Function: modang(x) @@ -299,6 +367,7 @@ function posmod(x,m) = (x%m+m)%m; // modang(270,360); // Returns: -90 // modang(700,360); // Returns: -20 function modang(x) = + assert( is_finite(x), "Input must be a finite number.") let(xx = posmod(x,360)) xx<180? xx : xx-360; @@ -306,7 +375,7 @@ function modang(x) = // Usage: // modrange(x, y, m, [step]) // Description: -// Returns a normalized list of values from `x` to `y`, by `step`, modulo `m`. Wraps if `x` > `y`. +// Returns a normalized list of numbers from `x` to `y`, by `step`, modulo `m`. Wraps if `x` > `y`. // Arguments: // x = The start value to constrain. // y = The end value to constrain. @@ -318,6 +387,7 @@ function modang(x) = // modrange(90,270,360, step=-45); // Returns: [90,45,0,315,270] // modrange(270,90,360, step=-45); // Returns: [270,225,180,135,90] function modrange(x, y, m, step=1) = + assert( is_finite(x+y+step+m) && !approx(m,0), "Input must be finite numbers. The module value cannot be zero.") let( a = posmod(x, m), b = posmod(y, m), @@ -330,20 +400,21 @@ function modrange(x, y, m, step=1) = // Function: rand_int() // Usage: -// rand_int(min,max,N,[seed]); +// rand_int(minval,maxval,N,[seed]); // Description: -// Return a list of random integers in the range of min to max, inclusive. +// Return a list of random integers in the range of minval to maxval, inclusive. // Arguments: -// min = Minimum integer value to return. -// max = Maximum integer value to return. +// minval = Minimum integer value to return. +// maxval = Maximum integer value to return. // N = Number of random integers to return. // seed = If given, sets the random number seed. // Example: // ints = rand_int(0,100,3); // int = rand_int(-10,10,1)[0]; -function rand_int(min, max, N, seed=undef) = - assert(max >= min, "Max value cannot be smaller than min") - let (rvect = is_def(seed) ? rands(min,max+1,N,seed) : rands(min,max+1,N)) +function rand_int(minval, maxval, N, seed=undef) = + assert( is_finite(minval+maxval+N) && (is_undef(seed) || is_finite(seed) ), "Input must be finite numbers.") + assert(maxval >= minval, "Max value cannot be smaller than minval") + let (rvect = is_def(seed) ? rands(minval,maxval+1,N,seed) : rands(minval,maxval+1,N)) [for(entry = rvect) floor(entry)]; @@ -358,6 +429,7 @@ function rand_int(min, max, N, seed=undef) = // N = Number of random numbers to return. Default: 1 // seed = If given, sets the random number seed. function gaussian_rands(mean, stddev, N=1, seed=undef) = + assert( is_finite(mean+stddev+N) && (is_undef(seed) || is_finite(seed) ), "Input must be finite numbers.") let(nums = is_undef(seed)? rands(0,1,N*2) : rands(0,1,N*2,seed)) [for (i = list_range(N)) mean + stddev*sqrt(-2*ln(nums[i*2]))*cos(360*nums[i*2+1])]; @@ -374,6 +446,10 @@ function gaussian_rands(mean, stddev, N=1, seed=undef) = // N = Number of random numbers to return. Default: 1 // seed = If given, sets the random number seed. function log_rands(minval, maxval, factor, N=1, seed=undef) = + assert( is_finite(minval+maxval+N) + && (is_undef(seed) || is_finite(seed) ) + && factor>0, + "Input must be finite numbers. `factor` should be greater than zero.") assert(maxval >= minval, "maxval cannot be smaller than minval") let( minv = 1-1/pow(factor,minval), @@ -395,18 +471,18 @@ function gcd(a,b) = b==0 ? abs(a) : gcd(b,a % b); -// Computes lcm for two scalars +// Computes lcm for two integers function _lcm(a,b) = - assert(is_int(a), "Invalid non-integer parameters to lcm") - assert(is_int(b), "Invalid non-integer parameters to lcm") - assert(a!=0 && b!=0, "Arguments to lcm must be nonzero") + assert(is_int(a) && is_int(b), "Invalid non-integer parameters to lcm") + assert(a!=0 && b!=0, "Arguments to lcm must be non zero") abs(a*b) / gcd(a,b); // Computes lcm for a list of values function _lcmlist(a) = - len(a)==1 ? a[0] : - _lcmlist(concat(slice(a,0,len(a)-2),[lcm(a[len(a)-2],a[len(a)-1])])); + len(a)==1 + ? a[0] + : _lcmlist(concat(slice(a,0,len(a)-2),[lcm(a[len(a)-2],a[len(a)-1])])); // Function: lcm() @@ -418,12 +494,11 @@ function _lcmlist(a) = // be non-zero integers. The output is always a positive integer. It is an error to pass zero // as an argument. function lcm(a,b=[]) = - !is_list(a) && !is_list(b) ? _lcm(a,b) : - let( - arglist = concat(force_list(a),force_list(b)) - ) - assert(len(arglist)>0,"invalid call to lcm with empty list(s)") - _lcmlist(arglist); + !is_list(a) && !is_list(b) + ? _lcm(a,b) + : let( arglist = concat(force_list(a),force_list(b)) ) + assert(len(arglist)>0, "Invalid call to lcm with empty list(s)") + _lcmlist(arglist); @@ -431,8 +506,9 @@ function lcm(a,b=[]) = // Function: sum() // Description: -// Returns the sum of all entries in the given list. -// If passed an array of vectors, returns a vector of sums of each part. +// Returns the sum of all entries in the given consistent list. +// If passed an array of vectors, returns the sum the vectors. +// If passed an array of matrices, returns the sum of the matrices. // If passed an empty list, the value of `dflt` will be returned. // Arguments: // v = The list to get the sum of. @@ -441,11 +517,10 @@ function lcm(a,b=[]) = // sum([1,2,3]); // returns 6. // sum([[1,2,3], [3,4,5], [5,6,7]]); // returns [9, 12, 15] function sum(v, dflt=0) = - is_vector(v) ? [for(i=v) 1]*v : + is_list(v) && len(v) == 0 ? dflt : + is_vector(v) || is_matrix(v)? [for(i=v) 1]*v : assert(is_consistent(v), "Input to sum is non-numeric or inconsistent") - is_vector(v[0]) ? [for(i=v) 1]*v : - len(v) == 0 ? dflt : - _sum(v,v[0]*0); + _sum(v,v[0]*0); function _sum(v,_total,_i=0) = _i>=len(v) ? _total : _sum(v,_total+v[_i], _i+1); @@ -495,37 +570,51 @@ function sum_of_squares(v) = sum(vmul(v,v)); // Examples: // v = sum_of_sines(30, [[10,3,0], [5,5.5,60]]); function sum_of_sines(a, sines) = - sum([ - for (s = sines) let( - ss=point3d(s), - v=ss.x*sin(a*ss.y+ss.z) - ) v - ]); + assert( is_finite(a) && is_matrix(sines,undef,3), "Invalid input.") + sum([ for (s = sines) + let( + ss=point3d(s), + v=ss[0]*sin(a*ss[1]+ss[2]) + ) v + ]); // Function: deltas() // Description: // Returns a list with the deltas of adjacent entries in the given list. +// The list should be a consistent list of numeric components (numbers, vectors, matrix, etc). // Given [a,b,c,d], returns [b-a,c-b,d-c]. // Arguments: // v = The list to get the deltas of. // Example: // deltas([2,5,9,17]); // returns [3,4,8]. // deltas([[1,2,3], [3,6,8], [4,8,11]]); // returns [[2,4,5], [1,2,3]] -function deltas(v) = [for (p=pair(v)) p.y-p.x]; +function deltas(v) = + assert( is_consistent(v) && len(v)>1 , "Inconsistent list or with length<=1.") + [for (p=pair(v)) p[1]-p[0]] ; // Function: product() // Description: // Returns the product of all entries in the given list. -// If passed an array of vectors, returns a vector of products of each part. -// If passed an array of matrices, returns a the resulting product matrix. +// If passed a list of vectors of same dimension, returns a vector of products of each part. +// If passed a list of square matrices, returns a the resulting product matrix. // Arguments: // v = The list to get the product of. // Example: // product([2,3,4]); // returns 24. // product([[1,2,3], [3,4,5], [5,6,7]]); // returns [15, 48, 105] -function product(v, i=0, tot=undef) = i>=len(v)? tot : product(v, i+1, ((tot==undef)? v[i] : is_vector(v[i])? vmul(tot,v[i]) : tot*v[i])); +function product(v) = + assert( is_vector(v) || is_matrix(v) || ( is_matrix(v[0],square=true) && is_consistent(v)), + "Invalid input.") + _product(v, 1, v[0]); + +function _product(v, i=0, _tot) = + i>=len(v) ? _tot : + _product( v, + i+1, + ( is_vector(v[i])? vmul(_tot,v[i]) : _tot*v[i] ) ); + // Function: outer_product() @@ -534,21 +623,22 @@ function product(v, i=0, tot=undef) = i>=len(v)? tot : product(v, i+1, ((tot==un // Usage: // M = outer_product(u,v); function outer_product(u,v) = - assert(is_vector(u) && is_vector(v)) - assert(len(u)==len(v)) - [for(i=[0:len(u)-1]) [for(j=[0:len(u)-1]) u[i]*v[j]]]; + assert(is_vector(u) && is_vector(v), "The inputs must be vectors.") + [for(ui=u) ui*v]; // Function: mean() // Description: -// Returns the arithmatic mean/average of all entries in the given array. +// Returns the arithmetic mean/average of all entries in the given array. // If passed a list of vectors, returns a vector of the mean of each part. // Arguments: // v = The list of values to get the mean of. // Example: // mean([2,3,4]); // returns 3. // mean([[1,2,3], [3,4,5], [5,6,7]]); // returns [3, 4, 5] -function mean(v) = sum(v)/len(v); +function mean(v) = + assert(is_list(v) && len(v)>0, "Invalid list.") + sum(v)/len(v); // Function: median() @@ -556,18 +646,33 @@ function mean(v) = sum(v)/len(v); // x = median(v); // Description: // Given a list of numbers or vectors, finds the median value or midpoint. -// If passed a list of vectors, returns the vector of the median of each part. +// If passed a list of vectors, returns the vector of the median of each component. function median(v) = - assert(is_list(v)) - assert(len(v)>0) - is_vector(v[0])? ( - assert(is_consistent(v)) - [ - for (i=idx(v[0])) - let(vals = subindex(v,i)) - (min(vals)+max(vals))/2 - ] - ) : (min(v)+max(v))/2; + is_vector(v) ? (min(v)+max(v))/2 : + is_matrix(v) ? [for(ti=transpose(v)) (min(ti)+max(ti))/2 ] + : assert(false , "Invalid input."); + +// Function: convolve() +// Usage: +// x = convolve(p,q); +// Description: +// Given two vectors, finds the convolution of them. +// The length of the returned vector is len(p)+len(q)-1 . +// Arguments: +// p = The first vector. +// q = The second vector. +// Example: +// a = convolve([1,1],[1,2,1]); // Returns: [1,3,3,1] +// b = convolve([1,2,3],[1,2,1])); // Returns: [1,4,8,8,3] +function convolve(p,q) = + p==[] || q==[] ? [] : + assert( is_vector(p) && is_vector(q), "The inputs should be vectors.") + let( n = len(p), + m = len(q)) + [for(i=[0:n+m-2], k1 = max(0,i-n+1), k2 = min(i,m-1) ) + [for(j=[k1:k2]) p[i-j] ] * [for(j=[k1:k2]) q[j] ] + ]; + // Section: Matrix math @@ -582,7 +687,7 @@ function median(v) = // want to solve Ax=b1 and Ax=b2 that you need to form the matrix transpose([b1,b2]) for the right hand side and then // transpose the returned value. function linear_solve(A,b) = - assert(is_matrix(A)) + assert(is_matrix(A), "Input should be a matrix.") let( m = len(A), n = len(A[0]) @@ -619,8 +724,12 @@ function matrix_inverse(A) = // Description: // Returns a submatrix with the specified index ranges or index sets. function submatrix(M,ind1,ind2) = - [for(i=ind1) [for(j=ind2) M[i][j] ] ]; - + assert( is_matrix(M), "Input must be a matrix." ) + [for(i=ind1) + [for(j=ind2) + assert( ! is_undef(M[i][j]), "Invalid indexing." ) + M[i][j] ] ]; + // Function: qr_factor() // Usage: qr = qr_factor(A) @@ -628,7 +737,7 @@ function submatrix(M,ind1,ind2) = // Calculates the QR factorization of the input matrix A and returns it as the list [Q,R]. This factorization can be // used to solve linear systems of equations. function qr_factor(A) = - assert(is_matrix(A)) + assert(is_matrix(A), "Input must be a matrix." ) let( m = len(A), n = len(A[0]) @@ -659,8 +768,8 @@ function _qr_factor(A,Q, column, m, n) = // Function: back_substitute() // Usage: back_substitute(R, b, [transpose]) // Description: -// Solves the problem Rx=b where R is an upper triangular square matrix. No check is made that the lower triangular entries -// are actually zero. If transpose==true then instead solve transpose(R)*x=b. +// Solves the problem Rx=b where R is an upper triangular square matrix. The lower triangular entries of R are +// ignored. If transpose==true then instead solve transpose(R)*x=b. // You can supply a compatible matrix b and it will produce the solution for every column of b. Note that if you want to // solve Rx=b1 and Rx=b2 you must set b to transpose([b1,b2]) and then take the transpose of the result. If the matrix // is singular (e.g. has a zero on the diagonal) then it returns []. @@ -694,7 +803,9 @@ function back_substitute(R, b, x=[],transpose = false) = // Example: // M = [ [6,-2], [1,8] ]; // det = det2(M); // Returns: 50 -function det2(M) = M[0][0] * M[1][1] - M[0][1]*M[1][0]; +function det2(M) = + assert( is_matrix(M,2,2), "Matrix should be 2x2." ) + M[0][0] * M[1][1] - M[0][1]*M[1][0]; // Function: det3() @@ -706,6 +817,7 @@ function det2(M) = M[0][0] * M[1][1] - M[0][1]*M[1][0]; // M = [ [6,4,-2], [1,-2,8], [1,5,7] ]; // det = det3(M); // Returns: -334 function det3(M) = + assert( is_matrix(M,3,3), "Matrix should be 3x3." ) M[0][0] * (M[1][1]*M[2][2]-M[2][1]*M[1][2]) - M[1][0] * (M[0][1]*M[2][2]-M[2][1]*M[0][2]) + M[2][0] * (M[0][1]*M[1][2]-M[1][1]*M[0][2]); @@ -720,21 +832,21 @@ function det3(M) = // M = [ [6,4,-2,9], [1,-2,8,3], [1,5,7,6], [4,2,5,1] ]; // det = determinant(M); // Returns: 2267 function determinant(M) = - assert(len(M)==len(M[0])) + assert(is_matrix(M,square=true), "Input should be a square matrix." ) len(M)==1? M[0][0] : len(M)==2? det2(M) : len(M)==3? det3(M) : sum( [for (col=[0:1:len(M)-1]) ((col%2==0)? 1 : -1) * - M[col][0] * - determinant( - [for (r=[1:1:len(M)-1]) - [for (c=[0:1:len(M)-1]) - if (c!=col) M[c][r] + M[col][0] * + determinant( + [for (r=[1:1:len(M)-1]) + [for (c=[0:1:len(M)-1]) + if (c!=col) M[c][r] + ] ] - ] - ) + ) ] ); @@ -753,8 +865,11 @@ function determinant(M) = // n = optional width of matrix // square = set to true to require a square matrix. Default: false function is_matrix(A,m,n,square=false) = - is_vector(A[0],n) && is_vector(A*(0*A[0]),m) && - (!square || len(A)==len(A[0])); + is_list(A[0]) +    && ( let(v = A*A[0]) is_num(0*(v*v)) ) // a matrix of finite numbers +    && (is_undef(n) || len(A[0])==n ) +    && (is_undef(m) || len(A)==m ) +    && ( !square || len(A)==len(A[0])); // Section: Comparisons and Logic @@ -774,11 +889,13 @@ function is_matrix(A,m,n,square=false) = // approx(0.3333,1/3); // Returns: false // approx(0.3333,1/3,eps=1e-3); // Returns: true // approx(PI,3.1415926536); // Returns: true -function approx(a,b,eps=EPSILON) = +function approx(a,b,eps=EPSILON) = a==b? true : a*0!=b*0? false : - is_list(a)? ([for (i=idx(a)) if(!approx(a[i],b[i],eps=eps)) 1] == []) : - (abs(a-b) <= eps); + is_list(a) + ? ([for (i=idx(a)) if( !approx(a[i],b[i],eps=eps)) 1] == []) + : is_num(a) && is_num(b) && (abs(a-b) <= eps); + function _type_num(x) = @@ -796,7 +913,7 @@ function _type_num(x) = // Description: // Compares two values. Lists are compared recursively. // Returns <0 if a0 if a>b. Returns 0 if a==b. -// If types are not the same, then undef < bool < num < str < list < range. +// If types are not the same, then undef < bool < nan < num < str < list < range. // Arguments: // a = First value to compare. // b = Second value to compare. @@ -820,13 +937,14 @@ function compare_vals(a, b) = // a = First list to compare. // b = Second list to compare. function compare_lists(a, b) = - a==b? 0 : let( - cmps = [ - for(i=[0:1:min(len(a),len(b))-1]) let( - cmp = compare_vals(a[i],b[i]) - ) if(cmp!=0) cmp - ] - ) cmps==[]? (len(a)-len(b)) : cmps[0]; + a==b? 0 + : let( + cmps = [ for(i=[0:1:min(len(a),len(b))-1]) + let( cmp = compare_vals(a[i],b[i]) ) + if(cmp!=0) cmp + ] + ) + cmps==[]? (len(a)-len(b)) : cmps[0]; // Function: any() @@ -843,12 +961,11 @@ function compare_lists(a, b) = // any([[0,0], [1,0]]); // Returns true. function any(l, i=0, succ=false) = (i>=len(l) || succ)? succ : - any( - l, i=i+1, succ=( - is_list(l[i])? any(l[i]) : - !(!l[i]) - ) - ); + any( l, + i+1, + succ = is_list(l[i]) ? any(l[i]) : !(!l[i]) + ); + // Function: all() @@ -865,13 +982,12 @@ function any(l, i=0, succ=false) = // all([[0,0], [1,0]]); // Returns false. // all([[1,1], [1,1]]); // Returns true. function all(l, i=0, fail=false) = - (i>=len(l) || fail)? (!fail) : - all( - l, i=i+1, fail=( - is_list(l[i])? !all(l[i]) : - !l[i] - ) - ); + (i>=len(l) || fail)? !fail : + all( l, + i+1, + fail = is_list(l[i]) ? !all(l[i]) : !l[i] + ) ; + // Function: count_true() @@ -904,6 +1020,21 @@ function count_true(l, nmax=undef, i=0, cnt=0) = ); +function count_true(l, nmax) = + !is_list(l) ? !(!l) ? 1: 0 : + let( c = [for( i = 0, + n = !is_list(l[i]) ? !(!l[i]) ? 1: 0 : undef, + c = !is_undef(n)? n : count_true(l[i], nmax), + s = c; + i=2, "Input `data` should have at least 2 elements.") + assert( is_finite(h) || is_vector(h), "The sampling `h` must be a number or a list of numbers." ) + assert( is_num(h) || len(h) == len(data)-(closed?0:1), + str("Vector valued `h` must have length ",len(data)-(closed?0:1))) is_vector(h) ? _deriv_nonuniform(data, h, closed=closed) : let( L = len(data) ) - closed? [ + closed + ? [ for(i=[0:1:L-1]) (data[(i+1)%L]-data[(L+i-1)%L])/2/h - ] : - let( - first = - L<3? data[1]-data[0] : - 3*(data[1]-data[0]) - (data[2]-data[1]), - last = - L<3? data[L-1]-data[L-2]: - (data[L-3]-data[L-2])-3*(data[L-2]-data[L-1]) - ) [ + ] + : let( + first = L<3 ? data[1]-data[0] : + 3*(data[1]-data[0]) - (data[2]-data[1]), + last = L<3 ? data[L-1]-data[L-2]: + (data[L-3]-data[L-2])-3*(data[L-2]-data[L-1]) + ) + [ first/2/h, for(i=[1:1:L-2]) (data[i+1]-data[i-1])/2/h, last/2/h - ]; + ]; function _dnu_calc(f1,fc,f2,h1,h2) = let( f1 = h2

=3, "Input list has less than 3 elements.") + assert( is_finite(h), "The sampling `h` must be a number." ) let( L = len(data) ) closed? [ for(i=[0:1:L-1]) @@ -1003,16 +1149,19 @@ function deriv2(data, h=1, closed=false) = // Computes a numerical third derivative estimate of the data, which may be scalar or vector valued. // The `h` parameter gives the step size of your sampling so the derivative can be scaled correctly. // If the `closed` parameter is true the data is assumed to be defined on a loop with data[0] adjacent to -// data[len(data)-1]. This function uses a five point derivative estimate, so the input must include five points: +// data[len(data)-1]. This function uses a five point derivative estimate, so the input data must include +// at least five points: // f'''(t) = (-f(t-2*h)+2*f(t-h)-2*f(t+h)+f(t+2*h)) / 2h^3. At the first and second points from the end // the estimates are f'''(t) = (-5*f(t)+18*f(t+h)-24*f(t+2*h)+14*f(t+3*h)-3*f(t+4*h)) / 2h^3 and // f'''(t) = (-3*f(t-h)+10*f(t)-12*f(t+h)+6*f(t+2*h)-f(t+3*h)) / 2h^3. function deriv3(data, h=1, closed=false) = + assert( is_consistent(data) , "Input list is not consistent or not numerical.") + assert( len(data)>=5, "Input list has less than 5 elements.") + assert( is_finite(h), "The sampling `h` must be a number." ) let( L = len(data), h3 = h*h*h ) - assert(L>=5, "Need five points for 3rd derivative estimate") closed? [ for(i=[0:1:L-1]) (-data[(L+i-2)%L]+2*data[(L+i-1)%L]-2*data[(i+1)%L]+data[(i+2)%L])/2/h3 @@ -1036,58 +1185,101 @@ function deriv3(data, h=1, closed=false) = // Function: C_times() // Usage: C_times(z1,z2) // Description: -// Multiplies two complex numbers. -function C_times(z1,z2) = [z1.x*z2.x-z1.y*z2.y,z1.x*z2.y+z1.y*z2.x]; +// Multiplies two complex numbers represented by 2D vectors. +function C_times(z1,z2) = + assert( is_vector(z1+z2,2), "Complex numbers should be represented by 2D vectors." ) + [ z1.x*z2.x - z1.y*z2.y, z1.x*z2.y + z1.y*z2.x ]; // Function: C_div() // Usage: C_div(z1,z2) // Description: -// Divides z1 by z2. -function C_div(z1,z2) = let(den = z2.x*z2.x + z2.y*z2.y) - [(z1.x*z2.x + z1.y*z2.y)/den, (z1.y*z2.x-z1.x*z2.y)/den]; +// Divides two complex numbers represented by 2D vectors. +function C_div(z1,z2) = + assert( is_vector(z1,2) && is_vector(z2), "Complex numbers should be represented by 2D vectors." ) + assert( !approx(z2,0), "The divisor `z2` cannot be zero." ) + let(den = z2.x*z2.x + z2.y*z2.y) + [(z1.x*z2.x + z1.y*z2.y)/den, (z1.y*z2.x - z1.x*z2.y)/den]; +// For the sake of consistence with Q_mul and vmul, C_times should be called C_mul // Section: Polynomials -// Function: polynomial() +// Function: polynomial() // Usage: // polynomial(p, z) // Description: // Evaluates specified real polynomial, p, at the complex or real input value, z. // The polynomial is specified as p=[a_n, a_{n-1},...,a_1,a_0] // where a_n is the z^n coefficient. Polynomial coefficients are real. +// The result is a number if `z` is a number and a complex number otherwise. // Note: this should probably be recoded to use division by [1,-z], which is more accurate // and avoids overflow with large coefficients, but requires poly_div to support complex coefficients. -function polynomial(p, z, k, zk, total) = - is_undef(k) ? polynomial(p, z, len(p)-1, is_num(z)? 1 : [1,0], is_num(z) ? 0 : [0,0]) : - k==-1 ? total : - polynomial(p, z, k-1, is_num(z) ? zk*z : C_times(zk,z), total+zk*p[k]); +function polynomial(p, z, _k, _zk, _total) = + is_undef(_k) + ? assert( is_vector(p), "Input polynomial coefficients must be a vector." ) + let(p = _poly_trim(p)) + assert( is_finite(z) || is_vector(z,2), "The value of `z` must be a real or a complex number." ) + polynomial( p, + z, + len(p)-1, + is_num(z)? 1 : [1,0], + is_num(z) ? 0 : [0,0]) + : _k==0 + ? _total + +_zk*p[0] + : polynomial( p, + z, + _k-1, + is_num(z) ? _zk*z : C_times(_zk,z), + _total+_zk*p[_k]); +function polynomial(p,z,k,total) = +     is_undef(k) +   ?    assert( is_vector(p) , "Input polynomial coefficients must be a vector." ) +        assert( is_finite(z) || is_vector(z,2), "The value of `z` must be a real or a complex number." ) +        polynomial( _poly_trim(p), z, 0, is_num(z) ? 0 : [0,0]) +   : k==len(p) ? total +   : polynomial(p,z,k+1, is_num(z) ? total*z+p[k] : C_times(total,z)+[p[k],0]); // Function: poly_mult() // Usage // polymult(p,q) // polymult([p1,p2,p3,...]) -// Descriptoin: +// Description: // Given a list of polynomials represented as real coefficient lists, with the highest degree coefficient first, // computes the coefficient list of the product polynomial. 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=0 && j=0 && j0 && d[0]!=0 , "Denominator is zero or has leading zero coefficient") - len(n)qlen ? p : q, - short = plen>qlen ? q : p - ) - _poly_trim(long + concat(repeat(0,len(long)-len(short)),short)); + assert( is_vector(p) && is_vector(q), "Invalid input polynomial(s)." ) + let( plen = len(p), + qlen = len(q), + long = plen>qlen ? p : q, + short = plen>qlen ? q : p + ) + _poly_trim(long + concat(repeat(0,len(long)-len(short)),short)); // Function: poly_roots() @@ -1150,38 +1347,38 @@ function poly_add(p,q) = // // Dario Bini. "Numerical computation of polynomial zeros by means of Aberth's Method", Numerical Algorithms, Feb 1996. // https://www.researchgate.net/publication/225654837_Numerical_computation_of_polynomial_zeros_by_means_of_Aberth's_method - function poly_roots(p,tol=1e-14,error_bound=false) = - assert(p!=[], "Input polynomial must have a nonzero coefficient") - assert(is_vector(p), "Input must be a vector") - p[0] == 0 ? poly_roots(slice(p,1,-1),tol=tol,error_bound=error_bound) : // Strip leading zero coefficients - p[len(p)-1] == 0 ? // Strip trailing zero coefficients - let( solutions = poly_roots(select(p,0,-2),tol=tol, error_bound=error_bound)) - (error_bound ? [ [[0,0], each solutions[0]], [0, each solutions[1]]] - : [[0,0], each solutions]) : - len(p)==1 ? (error_bound ? [[],[]] : []) : // Nonzero constant case has no solutions - len(p)==2 ? let( solution = [[-p[1]/p[0],0]]) // Linear case needs special handling - (error_bound ? [solution,[0]] : solution) - : - let( - n = len(p)-1, // polynomial degree - pderiv = [for(i=[0:n-1]) p[i]*(n-i)], - - s = [for(i=[0:1:n]) abs(p[i])*(4*(n-i)+1)], // Error bound polynomial from Bini + assert( is_vector(p), "Invalid polynomial." ) + let( p = _poly_trim(p,eps=0) ) + assert( p!=[0], "Input polynomial cannot be zero." ) + p[len(p)-1] == 0 ? // Strip trailing zero coefficients + let( solutions = poly_roots(select(p,0,-2),tol=tol, error_bound=error_bound)) + (error_bound ? [ [[0,0], each solutions[0]], [0, each solutions[1]]] + : [[0,0], each solutions]) : + len(p)==1 ? (error_bound ? [[],[]] : []) : // Nonzero constant case has no solutions + len(p)==2 ? let( solution = [[-p[1]/p[0],0]]) // Linear case needs special handling + (error_bound ? [solution,[0]] : solution) + : + let( + n = len(p)-1, // polynomial degree + pderiv = [for(i=[0:n-1]) p[i]*(n-i)], + + s = [for(i=[0:1:n]) abs(p[i])*(4*(n-i)+1)], // Error bound polynomial from Bini - // Using method from: http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/0915-24.pdf - beta = -p[1]/p[0]/n, - r = 1+pow(abs(polynomial(p,beta)/p[0]),1/n), - init = [for(i=[0:1:n-1]) // Initial guess for roots - let(angle = 360*i/n+270/n/PI) - [beta,0]+r*[cos(angle),sin(angle)] - ], - roots = _poly_roots(p,pderiv,s,init,tol=tol), - error = error_bound ? [for(xi=roots) n * (norm(polynomial(p,xi))+tol*polynomial(s,norm(xi))) / - abs(norm(polynomial(pderiv,xi))-tol*polynomial(s,norm(xi)))] : 0 - ) - error_bound ? [roots, error] : roots; + // Using method from: http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/0915-24.pdf + beta = -p[1]/p[0]/n, + r = 1+pow(abs(polynomial(p,beta)/p[0]),1/n), + init = [for(i=[0:1:n-1]) // Initial guess for roots + let(angle = 360*i/n+270/n/PI) + [beta,0]+r*[cos(angle),sin(angle)] + ], + roots = _poly_roots(p,pderiv,s,init,tol=tol), + error = error_bound ? [for(xi=roots) n * (norm(polynomial(p,xi))+tol*polynomial(s,norm(xi))) / + abs(norm(polynomial(pderiv,xi))-tol*polynomial(s,norm(xi)))] : 0 + ) + error_bound ? [roots, error] : roots; +// Internal function // p = polynomial // pderiv = derivative polynomial of p // z = current guess for the roots @@ -1222,12 +1419,16 @@ function _poly_roots(p, pderiv, s, z, tol, i=0) = // tol = tolerance for the complex polynomial root finder function real_roots(p,eps=undef,tol=1e-14) = - let( + assert( is_vector(p), "Invalid polynomial." ) + let( p = _poly_trim(p,eps=0) ) + assert( p!=[0], "Input polynomial cannot be zero." ) + let( roots_err = poly_roots(p,error_bound=true), roots = roots_err[0], err = roots_err[1] - ) - is_def(eps) ? [for(z=roots) if (abs(z.y)/(1+norm(z))