Fixed a bunch of undef math warnings with dev snapshot OpenSCAD builds.

This commit is contained in:
Garth Minette
2020-10-03 19:50:29 -07:00
parent e3ccf482fa
commit 16ee49e8b2
15 changed files with 215 additions and 161 deletions

232
math.scad
View File

@@ -56,7 +56,7 @@ function log2(x) =
// Function: hypot()
// Usage:
// l = hypot(x,y,[z]);
// l = hypot(x,y,<z>);
// Description:
// Calculate hypotenuse length of a 2D or 3D triangle.
// Arguments:
@@ -73,7 +73,7 @@ function hypot(x,y,z=0) =
// Function: factorial()
// Usage:
// x = factorial(n,[d]);
// x = factorial(n,<d>);
// Description:
// Returns the factorial of the given integer value, or n!/d! if d is given.
// Arguments:
@@ -373,7 +373,7 @@ function modang(x) =
// Function: modrange()
// Usage:
// modrange(x, y, m, [step])
// modrange(x, y, m, <step>)
// Description:
// Returns a normalized list of numbers from `x` to `y`, by `step`, modulo `m`. Wraps if `x` > `y`.
// Arguments:
@@ -401,7 +401,7 @@ function modrange(x, y, m, step=1) =
// Function: rand_int()
// Usage:
// rand_int(minval,maxval,N,[seed]);
// rand_int(minval,maxval,N,<seed>);
// Description:
// Return a list of random integers in the range of minval to maxval, inclusive.
// Arguments:
@@ -421,7 +421,7 @@ function rand_int(minval, maxval, N, seed=undef) =
// Function: gaussian_rands()
// Usage:
// gaussian_rands(mean, stddev, [N], [seed])
// gaussian_rands(mean, stddev, <N>, <seed>)
// Description:
// Returns a random number with a gaussian/normal distribution.
// Arguments:
@@ -437,7 +437,7 @@ function gaussian_rands(mean, stddev, N=1, seed=undef) =
// Function: log_rands()
// Usage:
// log_rands(minval, maxval, factor, [N], [seed]);
// log_rands(minval, maxval, factor, <N>, <seed>);
// Description:
// Returns a single random number, with a logarithmic distribution.
// Arguments:
@@ -506,6 +506,8 @@ function lcm(a,b=[]) =
// Section: Sums, Products, Aggregate Functions.
// Function: sum()
// Usage:
// x = sum(v, <dflt>);
// Description:
// Returns the sum of all entries in the given consistent list.
// If passed an array of vectors, returns the sum the vectors.
@@ -518,8 +520,7 @@ 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_list(v) && len(v) == 0 ? dflt :
is_vector(v) || is_matrix(v)? [for(i=v) 1]*v :
v==[]? dflt :
assert(is_consistent(v), "Input to sum is non-numeric or inconsistent")
_sum(v,v[0]*0);
@@ -527,6 +528,8 @@ function _sum(v,_total,_i=0) = _i>=len(v) ? _total : _sum(v,_total+v[_i], _i+1);
// Function: cumsum()
// Usage:
// sums = cumsum(v);
// Description:
// Returns a list where each item is the cumulative sum of all items up to and including the corresponding entry in the input list.
// If passed an array of vectors, returns a list of cumulative vectors sums.
@@ -573,6 +576,8 @@ function sum_of_sines(a, sines) =
// Function: deltas()
// Usage:
// delts = deltas(v);
// 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).
@@ -588,6 +593,8 @@ function deltas(v) =
// Function: product()
// Usage:
// x = product(v);
// Description:
// Returns the product of all entries in the given list.
// If passed a list of vectors of same dimension, returns a vector of products of each part.
@@ -611,6 +618,8 @@ function _product(v, i=0, _tot) =
// Function: outer_product()
// Usage:
// x = outer_product(u,v);
// Description:
// Compute the outer product of two vectors, a matrix.
// Usage:
@@ -621,6 +630,8 @@ function outer_product(u,v) =
// Function: mean()
// Usage:
// x = mean(v);
// Description:
// 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.
@@ -660,14 +671,15 @@ function convolve(p,q) =
// Section: Matrix math
// Function: linear_solve()
// Usage: linear_solve(A,b)
// Usage:
// solv = linear_solve(A,b)
// Description:
// Solves the linear system Ax=b. If A is square and non-singular the unique solution is returned. If A is overdetermined
// the least squares solution is returned. If A is underdetermined, the minimal norm solution is returned.
// If A is rank deficient or singular then linear_solve returns []. If b is a matrix that is compatible with A
// Solves the linear system Ax=b. If `A` is square and non-singular the unique solution is returned. If `A` is overdetermined
// the least squares solution is returned. If `A` is underdetermined, the minimal norm solution is returned.
// If `A` is rank deficient or singular then linear_solve returns `[]`. If `b` is a matrix that is compatible with `A`
// then the problem is solved for the matrix valued right hand side and a matrix is returned. Note that if you
// want to solve Ax=b1 and Ax=b2 that you need to form the matrix transpose([b1,b2]) for the right hand side and then
// transpose the returned value.
// want to solve Ax=b1 and Ax=b2 that you need to form the matrix `transpose([b1,b2])` for the right hand side and then
// transpose the returned value.
function linear_solve(A,b,pivot=true) =
assert(is_matrix(A), "Input should be a matrix.")
let(
@@ -690,33 +702,32 @@ function linear_solve(A,b,pivot=true) =
// Function: matrix_inverse()
// Usage:
// matrix_inverse(A)
// mat = matrix_inverse(A)
// Description:
// Compute the matrix inverse of the square matrix A. If A is singular, returns undef.
// Note that if you just want to solve a linear system of equations you should NOT
// use this function. Instead use linear_solve, or use qr_factor. The computation
// Compute the matrix inverse of the square matrix `A`. If `A` is singular, returns `undef`.
// Note that if you just want to solve a linear system of equations you should NOT use this function.
// Instead use [[`linear_solve()`|linear_solve]], or use [[`qr_factor()`|qr_factor]]. The computation
// will be faster and more accurate.
function matrix_inverse(A) =
assert(is_matrix(A,square=true),"Input to matrix_inverse() must be a square matrix")
assert(is_matrix(A) && len(A)==len(A[0]),"Input to matrix_inverse() must be a square matrix")
linear_solve(A,ident(len(A)));
// Function: null_space()
// Usage:
// null_space(A)
// x = null_space(A)
// Description:
// Returns an orthonormal basis for the null space of A, namely the vectors {x} such that Ax=0. If the null space
// is just the origin then returns an empty list.
// Returns an orthonormal basis for the null space of `A`, namely the vectors {x} such that Ax=0.
// If the null space is just the origin then returns an empty list.
function null_space(A,eps=1e-12) =
assert(is_matrix(A))
let(
Q_R=qr_factor(transpose(A),pivot=true),
R=Q_R[1],
zrow = [for(i=idx(R)) if (all_zero(R[i],eps)) i]
Q_R = qr_factor(transpose(A),pivot=true),
R = Q_R[1],
zrow = [for(i=idx(R)) if (all_zero(R[i],eps)) i]
)
len(zrow)==0
? []
: transpose(subindex(Q_R[0],zrow));
len(zrow)==0 ? [] :
transpose(subindex(Q_R[0],zrow));
// Function: qr_factor()
@@ -730,19 +741,18 @@ function null_space(A,eps=1e-12) =
function qr_factor(A, pivot=false) =
assert(is_matrix(A), "Input must be a matrix." )
let(
m = len(A),
n = len(A[0])
m = len(A),
n = len(A[0])
)
let(
qr =_qr_factor(A, Q=ident(m),P=ident(n), pivot=pivot, column=0, m = m, n=n),
Rzero =
let( R = qr[1] )
[ for(i=[0:m-1]) [
let( ri = R[i] )
for(j=[0:n-1]) i>j ? 0 : ri[j]
qr = _qr_factor(A, Q=ident(m),P=ident(n), pivot=pivot, column=0, m = m, n=n),
Rzero = let( R = qr[1]) [
for(i=[0:m-1]) [
let( ri = R[i] )
for(j=[0:n-1]) i>j ? 0 : ri[j]
]
]
) [qr[0],Rzero,qr[2]];
]
) [qr[0], Rzero, qr[2]];
function _qr_factor(A,Q,P, pivot, column, m, n) =
column >= min(m-1,n) ? [Q,A,P] :
@@ -770,7 +780,7 @@ function _swap_matrix(n,i,j) =
// Function: back_substitute()
// Usage: back_substitute(R, b, [transpose])
// Usage: back_substitute(R, b, <transpose>)
// Description:
// Solves the problem Rx=b where R is an upper triangular square matrix. The lower triangular entries of R are
// ignored. If transpose==true then instead solve transpose(R)*x=b.
@@ -799,6 +809,8 @@ function _back_substitute(R, b, x=[]) =
// Function: det2()
// Usage:
// d = det2(M);
// Description:
// Optimized function that returns the determinant for the given 2x2 square matrix.
// Arguments:
@@ -807,11 +819,13 @@ function _back_substitute(R, b, x=[]) =
// M = [ [6,-2], [1,8] ];
// det = det2(M); // Returns: 50
function det2(M) =
assert( 0*M==[[0,0],[0,0]], "Matrix should be 2x2." )
assert(is_matrix(M,2,2), "Matrix must be 2x2.")
M[0][0] * M[1][1] - M[0][1]*M[1][0];
// Function: det3()
// Usage:
// d = det3(M);
// Description:
// Optimized function that returns the determinant for the given 3x3 square matrix.
// Arguments:
@@ -820,13 +834,15 @@ function det2(M) =
// M = [ [6,4,-2], [1,-2,8], [1,5,7] ];
// det = det3(M); // Returns: -334
function det3(M) =
assert( 0*M==[[0,0,0],[0,0,0],[0,0,0]], "Matrix should be 3x3." )
assert(is_matrix(M,3,3), "Matrix must be 3x3.")
M[0][0] * (M[1][1]*M[2][2]-M[2][1]*M[1][2]) -
M[1][0] * (M[0][1]*M[2][2]-M[2][1]*M[0][2]) +
M[2][0] * (M[0][1]*M[1][2]-M[1][1]*M[0][2]);
// Function: determinant()
// Usage:
// d = determinant(M);
// Description:
// Returns the determinant for the given square matrix.
// Arguments:
@@ -835,7 +851,7 @@ 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(is_matrix(M,square=true), "Input should be a square matrix." )
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) :
@@ -856,23 +872,22 @@ function determinant(M) =
// Function: is_matrix()
// Usage:
// is_matrix(A,[m],[n],[square])
// is_matrix(A,<m>,<n>,<square>)
// Description:
// Returns true if A is a numeric matrix of height m and width n. If m or n
// are omitted or set to undef then true is returned for any positive dimension.
// If `square` is true then the matrix is required to be square.
// specify m != n and require a square matrix then the result will always be false.
// Arguments:
// A = matrix to test
// m = optional height of matrix
// n = optional width of matrix
// square = set to true to require a square matrix. Default: false
// A = The matrix to test.
// m = Is given, requires the matrix to have the given height.
// n = Is given, requires the matrix to have the given width.
// square = If true, requires the matrix to have a width equal to its height. Default: false
function is_matrix(A,m,n,square=false) =
is_list(A[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]));
is_list(A)
&& (( is_undef(m) && len(A) ) || len(A)==m)
&& is_list(A[0])
&& (( is_undef(n) && len(A[0]) ) || len(A[0])==n)
&& (!square || len(A) == len(A[0]))
&& is_consistent(A);
// Function: norm_fro()
@@ -891,7 +906,7 @@ function norm_fro(A) =
// Function: all_zero()
// Usage:
// all_zero(x);
// x = all_zero(x, <eps>);
// Description:
// Returns true if the finite number passed to it is approximately zero, to within `eps`.
// If passed a list, recursively checks if all items in the list are approximately zero.
@@ -912,7 +927,7 @@ function all_zero(x, eps=EPSILON) =
// Function: all_nonzero()
// Usage:
// all_nonzero(x);
// x = all_nonzero(x, <eps>);
// Description:
// Returns true if the finite number passed to it is not almost zero, to within `eps`.
// If passed a list, recursively checks if all items in the list are not almost zero.
@@ -1030,7 +1045,7 @@ function all_nonnegative(x) =
// Function: approx()
// Usage:
// approx(a,b,[eps])
// b = approx(a,b,<eps>)
// Description:
// Compares two numbers or vectors, and returns true if they are closer than `eps` to each other.
// Arguments:
@@ -1044,12 +1059,9 @@ function all_nonnegative(x) =
// approx(0.3333,1/3,eps=1e-3); // Returns: true
// approx(PI,3.1415926536); // Returns: true
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] == [])
: is_num(a) && is_num(b) && (abs(a-b) <= eps);
(a==b && is_bool(a) == is_bool(b)) ||
(is_num(a) && is_num(b) && abs(a-b) <= eps) ||
(is_list(a) && is_list(b) && len(a) == len(b) && [] == [for (i=idx(a)) if (!approx(a[i],b[i],eps=eps)) 1]);
function _type_num(x) =
@@ -1063,7 +1075,7 @@ function _type_num(x) =
// Function: compare_vals()
// Usage:
// compare_vals(a, b);
// b = compare_vals(a, b);
// Description:
// Compares two values. Lists are compared recursively.
// Returns <0 if a<b. Returns >0 if a>b. Returns 0 if a==b.
@@ -1081,7 +1093,7 @@ function compare_vals(a, b) =
// Function: compare_lists()
// Usage:
// compare_lists(a, b)
// b = compare_lists(a, b)
// Description:
// Compare contents of two lists using `compare_vals()`.
// Returns <0 if `a`<`b`.
@@ -1102,6 +1114,8 @@ function compare_lists(a, b) =
// Function: any()
// Usage:
// b = any(l);
// Description:
// Returns true if any item in list `l` evaluates as true.
// If `l` is a lists of lists, `any()` is applied recursively to each sublist.
@@ -1115,17 +1129,19 @@ function compare_lists(a, b) =
// any([[0,0], [1,0]]); // Returns true.
function any(l) =
assert(is_list(l), "The input is not a list." )
_any(l, i=0, succ=false);
_any(l);
function _any(l, i=0, succ=false) =
(i>=len(l) || succ)? succ :
_any( l,
i+1,
succ = is_list(l[i]) ? _any(l[i]) : !(!l[i])
);
_any(
l, i+1,
succ = is_list(l[i]) ? _any(l[i]) : !(!l[i])
);
// Function: all()
// Usage:
// b = all(l);
// Description:
// Returns true if all items in list `l` evaluate as true.
// If `l` is a lists of lists, `all()` is applied recursively to each sublist.
@@ -1138,21 +1154,21 @@ function _any(l, i=0, succ=false) =
// all([[0,0], [0,0]]); // Returns false.
// all([[0,0], [1,0]]); // Returns false.
// all([[1,1], [1,1]]); // Returns true.
function all(l, i=0, fail=false) =
function all(l) =
assert( is_list(l), "The input is not a list." )
_all(l, i=0, fail=false);
_all(l);
function _all(l, i=0, fail=false) =
(i>=len(l) || fail)? !fail :
_all( l,
i+1,
fail = is_list(l[i]) ? !_all(l[i]) : !l[i]
) ;
_all(
l, i+1,
fail = is_list(l[i]) ? !_all(l[i]) : !l[i]
) ;
// Function: count_true()
// Usage:
// count_true(l)
// n = count_true(l)
// Description:
// Returns the number of items in `l` that evaluate as true.
// If `l` is a lists of lists, this is applied recursively to each
@@ -1170,26 +1186,22 @@ function _all(l, i=0, fail=false) =
// count_true([[0,0], [1,0]]); // Returns 1.
// count_true([[1,1], [1,1]]); // Returns 4.
// count_true([[1,1], [1,1]], nmax=3); // Returns 3.
function _count_true_rec(l, nmax, _cnt=0, _i=0) =
_i>=len(l) || (is_num(nmax) && _cnt>=nmax)? _cnt :
_count_true_rec(l, nmax, _cnt=_cnt+(l[_i]?1:0), _i=_i+1);
function count_true(l, nmax) =
!is_list(l) ? !(!l) ? 1: 0 :
let( c = [for( i = 0,
n = !is_list(l[i]) ? !(!l[i]) ? 1: 0 : undef,
c = !is_undef(n)? n : count_true(l[i], nmax),
s = c;
i<len(l) && (is_undef(nmax) || s<nmax);
i = i+1,
n = !is_list(l[i]) ? !(!l[i]) ? 1: 0 : undef,
c = !is_undef(n) || (i==len(l))? n : count_true(l[i], nmax-s),
s = s+c
) s ] )
len(c)<len(l)? nmax: c[len(c)-1];
is_undef(nmax)? len([for (x=l) if(x) 1]) :
!is_list(l) ? ( l? 1: 0) :
_count_true_rec(l, nmax);
// Section: Calculus
// Function: deriv()
// Usage: deriv(data, [h], [closed])
// Usage:
// x = deriv(data, [h], [closed])
// Description:
// Computes a numerical 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.
@@ -1253,7 +1265,8 @@ function _deriv_nonuniform(data, h, closed) =
// Function: deriv2()
// Usage: deriv2(data, [h], [closed])
// Usage:
// x = deriv2(data, <h>, <closed>)
// Description:
// Computes a numerical estimate of the second derivative of the data, which may be scalar or vector valued.
// The `h` parameter gives the step size of your sampling so the derivative can be scaled correctly.
@@ -1296,7 +1309,8 @@ function deriv2(data, h=1, closed=false) =
// Function: deriv3()
// Usage: deriv3(data, [h], [closed])
// Usage:
// x = deriv3(data, <h>, <closed>)
// Description:
// 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.
@@ -1306,6 +1320,10 @@ function deriv2(data, h=1, closed=false) =
// 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.
// Arguments:
// data = the list of the elements to compute the derivative of.
// h = the constant parametric sampling of the data.
// closed = boolean to indicate if the data set should be wrapped around from the end to the start.
function 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.")
@@ -1335,17 +1353,27 @@ function deriv3(data, h=1, closed=false) =
// Section: Complex Numbers
// Function: C_times()
// Usage: C_times(z1,z2)
// Usage:
// c = C_times(z1,z2)
// Description:
// Multiplies two complex numbers represented by 2D vectors.
// Returns a complex number as a 2D vector [REAL, IMAGINARY].
// Arguments:
// z1 = First complex number, given as a 2D vector [REAL, IMAGINARY]
// z2 = Second complex number, given as a 2D vector [REAL, IMAGINARY]
function C_times(z1,z2) =
assert( is_matrix([z1,z2],2,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)
// Usage:
// x = C_div(z1,z2)
// Description:
// Divides two complex numbers represented by 2D vectors.
// Returns a complex number as a 2D vector [REAL, IMAGINARY].
// Arguments:
// z1 = First complex number, given as a 2D vector [REAL, IMAGINARY]
// z2 = Second complex number, given as a 2D vector [REAL, IMAGINARY]
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." )
@@ -1358,16 +1386,14 @@ function C_div(z1,z2) =
// Function: quadratic_roots()
// Usage:
// roots = quadratic_roots(a,b,c,[real])
// roots = quadratic_roots(a,b,c,<real>)
// Description:
// Computes roots of the quadratic equation a*x^2+b*x+c==0, where the
// coefficients are real numbers. If real is true then returns only the
// real roots. Otherwise returns a pair of complex values. This method
// may be more reliable than the general root finder at distinguishing
// real roots from complex roots.
// https://people.csail.mit.edu/bkph/articles/Quadratics.pdf
// Algorithm from: https://people.csail.mit.edu/bkph/articles/Quadratics.pdf
function quadratic_roots(a,b,c,real=false) =
real ? [for(root = quadratic_roots(a,b,c,real=false)) if (root.y==0) root.x]
:
@@ -1393,7 +1419,7 @@ function quadratic_roots(a,b,c,real=false) =
// Function: polynomial()
// Usage:
// polynomial(p, z)
// x = 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]
@@ -1409,8 +1435,8 @@ function polynomial(p,z,k,total) =
// Function: poly_mult()
// Usage:
// polymult(p,q)
// polymult([p1,p2,p3,...])
// x = polymult(p,q)
// x = polymult([p1,p2,p3,...])
// 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.
@@ -1481,7 +1507,7 @@ function poly_add(p,q) =
// Function: poly_roots()
// Usage:
// poly_roots(p,[tol])
// poly_roots(p,<tol>)
// Description:
// Returns all complex roots of the specified real polynomial p.
// The polynomial is specified as p=[a_n, a_{n-1},...,a_1,a_0]
@@ -1551,7 +1577,7 @@ function _poly_roots(p, pderiv, s, z, tol, i=0) =
// Function: real_roots()
// Usage:
// real_roots(p, [eps], [tol])
// real_roots(p, <eps>, <tol>)
// Description:
// Returns the real roots of the specified real polynomial p.
// The polynomial is specified as p=[a_n, a_{n-1},...,a_1,a_0]