diff --git a/math.scad b/math.scad index 010c307..2733f06 100644 --- a/math.scad +++ b/math.scad @@ -377,9 +377,6 @@ function log_rands(minval, maxval, factor, N=1, seed=undef) = // Section: GCD/GCF, LCM -// If argument is a list return it. Otherwise return a singleton list containing the argument. -function _force_list(x) = is_list(x) ? x : [x]; - // Function: gcd() // Usage: // gcd(a,b) @@ -415,7 +412,7 @@ function _lcmlist(a) = function lcm(a,b=[]) = !is_list(a) && !is_list(b) ? _lcm(a,b) : let( - arglist = concat(_force_list(a),_force_list(b)) + arglist = concat((is_list(a)?a:[a]), (is_list(b)?b:[b])) ) assert(len(arglist)>0,"invalid call to lcm with empty list(s)") _lcmlist(arglist); @@ -539,34 +536,29 @@ function mean(v) = sum(v)/len(v); // 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`. function linear_solve(A,b) = - let( - dim = array_dim(A), - m=dim[0], n=dim[1] - ) - assert(len(b)==m,str("Incompatible matrix and vector",dim,len(b))) - let ( - qr = m<n ? qr_factor(transpose(A)) : qr_factor(A), - maxdim = max(n,m), - mindim = min(n,m), - Q = submatrix(qr[0],[0:maxdim-1], [0:mindim-1]), - R = submatrix(qr[1],[0:mindim-1], [0:mindim-1]), - zeros = [for(i=[0:mindim-1]) if (approx(R[i][i],0)) i] - ) - zeros != [] ? undef : - m<n ? Q*back_substitute(R,b,transpose=true) : - back_substitute(R, transpose(Q)*b); + let( + dim = array_dim(A), + m=dim[0], n=dim[1] + ) + assert(len(b)==m,str("Incompatible matrix and vector",dim,len(b))) + let ( + qr = m<n ? qr_factor(transpose(A)) : qr_factor(A), + maxdim = max(n,m), + mindim = min(n,m), + Q = submatrix(qr[0],[0:maxdim-1], [0:mindim-1]), + R = submatrix(qr[1],[0:mindim-1], [0:mindim-1]), + zeros = [for(i=[0:mindim-1]) if (approx(R[i][i],0)) i] + ) + zeros != [] ? undef : + m<n ? Q*back_substitute(R,b,transpose=true) : + back_substitute(R, transpose(Q)*b); // Function: submatrix() // 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() @@ -575,30 +567,34 @@ 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) = - let( - dim = array_dim(A), - m = dim[0], - n = dim[1] - ) - assert(len(dim)==2) - let( - qr =_qr_factor(A, column=0, m = m, n=m, Q=ident(m)), - Rzero = [for(i=[0:m-1]) [for(j=[0:n-1]) i>j ? 0 : qr[1][i][j]]] - ) - [qr[0],Rzero]; + let( + dim = array_dim(A), + m = dim[0], + n = dim[1] + ) + assert(len(dim)==2) + let( + qr =_qr_factor(A, column=0, m = m, n=m, Q=ident(m)), + Rzero = [ + for(i=[0:m-1]) [ + for(j=[0:n-1]) + i>j ? 0 : qr[1][i][j] + ] + ] + ) [qr[0],Rzero]; function _qr_factor(A,Q, column, m, n) = - column >= min(m-1,n) ? [Q,A] : - let( - x = [for(i=[column:1:m-1]) A[i][column]], - alpha = (x[0]<=0 ? 1 : -1) * norm(x), - u = x - concat([alpha],replist(0,m-1)), - v = u / norm(u), - Qc = ident(len(x)) - 2*transpose([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); - + column >= min(m-1,n) ? [Q,A] : + let( + x = [for(i=[column:1:m-1]) A[i][column]], + alpha = (x[0]<=0 ? 1 : -1) * norm(x), + u = x - concat([alpha],replist(0,m-1)), + v = u / norm(u), + Qc = ident(len(x)) - 2*transpose([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); + // Function: back_substitute() @@ -607,18 +603,19 @@ 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. function back_substitute(R, b, x=[],transpose = false) = - let(n=len(b)) - transpose ? - reverse(back_substitute( - [for(i=[0:n-1]) [for(j=[0:n-1]) R[n-1-j][n-1-i]]], - reverse(b), x, false)) : - len(x) == n ? x : - let( - ind = n - len(x) - 1, - newvalue = len(x)==0 ? b[ind]/R[ind][ind] : - (b[ind]-select(R[ind],ind+1,-1) * x)/R[ind][ind] - ) - back_substitute(R, b, concat([newvalue],x)); + let(n=len(b)) + transpose? + reverse(back_substitute( + [for(i=[0:n-1]) [for(j=[0:n-1]) R[n-1-j][n-1-i]]], + reverse(b), x, false + )) : + len(x) == n ? x : + let( + ind = n - len(x) - 1, + newvalue = + len(x)==0? b[ind]/R[ind][ind] : + (b[ind]-select(R[ind],ind+1,-1) * x)/R[ind][ind] + ) back_substitute(R, b, concat([newvalue],x)); // Function: det2() @@ -814,6 +811,8 @@ function count_true(l, nmax=undef, i=0, cnt=0) = ) ); + + // Section: Calculus // Function: deriv() @@ -826,19 +825,23 @@ function count_true(l, nmax=undef, i=0, cnt=0) = // for internal points, f'(t) = (f(t+h)-f(t-h))/2h. For the endpoints (when closed=false) the algorithm // uses a two point method if sufficient points are available: f'(t) = (3*(f(t+h)-f(t)) - (f(t+2*h)-f(t+h)))/2h. function deriv(data, h=1, closed=false) = - let( L = len(data) ) - 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]) - ) - [ first/2/h, - for(i=[1:1:L-2]) (data[i+1]-data[i-1])/2/h, - last/2/h]; + let( L = len(data) ) + 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]) + ) [ + first/2/h, + for(i=[1:1:L-2]) (data[i+1]-data[i-1])/2/h, + last/2/h + ]; // Function: deriv2() @@ -853,22 +856,25 @@ function deriv(data, h=1, closed=false) = // f''(t) = (2*f(t) - 5*f(t+h) + 4*f(t+2*h) - f(t+3*h))/h^2 or if five points are available // f''(t) = (35*f(t) - 104*f(t+h) + 114*f(t+2*h) - 56*f(t+3*h) + 11*f(t+4*h)) / 12h^2 function deriv2(data, h=1, closed=false) = - let( L = len(data) ) - closed ? - [ for(i=[0:1:L-1]) (data[(i+1)%L]-2*data[i]+data[(L+i-1)%L])/h/h ] : - let( first = L<3 ? undef : - L==3 ? data[0] - 2*data[1] + data[2] : - L==4 ? 2*data[0] - 5*data[1] + 4*data[2] - data[3] : - (35*data[0] - 104*data[1] + 114*data[2] - 56*data[3] + 11*data[4])/12, - last = L<3 ? undef : - L==3 ? data[L-1] - 2*data[L-2] + data[L-3] : - L==4 ? -2*data[L-1] + 5*data[L-2] - 4*data[L-3] + data[L-4] : - (35*data[L-1] - 104*data[L-2] + 114*data[L-3] - 56*data[L-4] + 11*data[L-5])/12 - ) - [ first/h/h, - for(i=[1:1:L-2]) (data[i+1]-2*data[i]+data[i-1])/h/h, - last/h/h]; - + let( L = len(data) ) + closed? [ + for(i=[0:1:L-1]) + (data[(i+1)%L]-2*data[i]+data[(L+i-1)%L])/h/h + ] : + let( + first = L<3? undef : + L==3? data[0] - 2*data[1] + data[2] : + L==4? 2*data[0] - 5*data[1] + 4*data[2] - data[3] : + (35*data[0] - 104*data[1] + 114*data[2] - 56*data[3] + 11*data[4])/12, + last = L<3? undef : + L==3? data[L-1] - 2*data[L-2] + data[L-3] : + L==4? -2*data[L-1] + 5*data[L-2] - 4*data[L-3] + data[L-4] : + (35*data[L-1] - 104*data[L-2] + 114*data[L-3] - 56*data[L-4] + 11*data[L-5])/12 + ) [ + first/h/h, + for(i=[1:1:L-2]) (data[i+1]-2*data[i]+data[i-1])/h/h, + last/h/h + ]; // Function: deriv3() @@ -882,25 +888,28 @@ function deriv2(data, h=1, closed=false) = // 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) = - 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] : - let( - first=(-5*data[0]+18*data[1]-24*data[2]+14*data[3]-3*data[4])/2, - second=(-3*data[0]+10*data[1]-12*data[2]+6*data[3]-data[4])/2, - last=(5*data[L-1]-18*data[L-2]+24*data[L-3]-14*data[L-4]+3*data[L-5])/2, - prelast=(3*data[L-1]-10*data[L-2]+12*data[L-3]-6*data[L-4]+data[L-5])/2 - ) - [ - first/h3, - second/h3, - for(i=[2:1:L-3]) (-data[i-2]+2*data[i-1]-2*data[i+1]+data[i+2])/2/h3, - prelast/h3, - last/h3 - ]; + 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 + ] : + let( + first=(-5*data[0]+18*data[1]-24*data[2]+14*data[3]-3*data[4])/2, + second=(-3*data[0]+10*data[1]-12*data[2]+6*data[3]-data[4])/2, + last=(5*data[L-1]-18*data[L-2]+24*data[L-3]-14*data[L-4]+3*data[L-5])/2, + prelast=(3*data[L-1]-10*data[L-2]+12*data[L-3]-6*data[L-4]+data[L-5])/2 + ) [ + first/h3, + second/h3, + for(i=[2:1:L-3]) (-data[i-2]+2*data[i-1]-2*data[i+1]+data[i+2])/2/h3, + prelast/h3, + last/h3 + ]; + // vim: noexpandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap diff --git a/version.scad b/version.scad index 80fa107..889a910 100644 --- a/version.scad +++ b/version.scad @@ -8,7 +8,7 @@ ////////////////////////////////////////////////////////////////////// -BOSL_VERSION = [2,0,140]; +BOSL_VERSION = [2,0,141]; // Section: BOSL Library Version Functions