grey/BOSL2
1
0
Fork 0
mirror of https://github.com/BelfrySCAD/BOSL2.git synced 2024-12-29 16:29:40 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge pull request #201 from adrianVmariano/master

Browse source 
Updates
This commit is contained in:
Revar Desmera 2020-07-13 09:35:21 -07:00 committed by GitHub
commit 0ade5d5baa
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
2 changed files with 444 additions and 23 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

147
math.scad
View file

@ -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

320
tests/test_math.scad
View file

@ -76,6 +76,12 @@ module test_is_matrix() {
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));
@ -181,6 +187,9 @@ module test_sqr() {
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();
@ -525,6 +534,9 @@ module test_any() {
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();
@ -536,6 +548,9 @@ module test_all() {
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();
@ -554,6 +569,7 @@ test_count_true();
module test_factorial() {
assert_equal(factorial(0), 1);
assert_equal(factorial(1), 1);
assert_equal(factorial(2), 2);
assert_equal(factorial(3), 6);
@ -562,6 +578,8 @@ module test_factorial() {
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_factorial();
@ -570,6 +588,11 @@ 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_gcd();
@ -578,9 +601,306 @@ 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_lcm();
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_C_times();
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_back_substitute();
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_linear_solve();
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_outer_product();
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_deriv();
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_deriv2();
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_deriv3();
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_polynomial();
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_poly_roots();
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_real_roots();
// 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_submatrix();
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_qr_factor();
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_poly_mult();
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_poly_div();
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_poly_add();
// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap