mirror of
https://github.com/revarbat/BOSL2.git
synced 2025-01-16 13:50:23 +01:00
Cleaned up code formatting in math.scad.
This commit is contained in:
Revar Desmera
parent
f250f30044
commit
28114b49b5
225
math.scad
225
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
|
||||
|
||||
@ -8,7 +8,7 @@
|
||||
//////////////////////////////////////////////////////////////////////
|
||||
|
||||
|
||||
BOSL_VERSION = [2,0,140];
|
||||
BOSL_VERSION = [2,0,141];
|
||||
|
||||
|
||||
// Section: BOSL Library Version Functions
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user