Merge pull request from RonaldoCMP/master

Extensive changes in arrays.scad, vectors.scad, common.scad and their regression test codes
This commit is contained in:
Revar Desmera 2020-07-29 21:44:39 -07:00 committed by GitHub
commit 5e98cbf679
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
9 changed files with 1133 additions and 686 deletions

File diff suppressed because it is too large Load diff

View file

@ -15,7 +15,8 @@
// Usage:
// typ = typeof(x);
// Description:
// Returns a string representing the type of the value. One of "undef", "boolean", "number", "nan", "string", "list", or "range"
// Returns a string representing the type of the value. One of "undef", "boolean", "number", "nan", "string", "list", "range" or "invalid".
// Some malformed "ranges", like '[0:NAN:INF]' and '[0:"a":INF]', may be classified as "undef" or "invalid".
function typeof(x) =
is_undef(x)? "undef" :
is_bool(x)? "boolean" :
@ -23,7 +24,9 @@ function typeof(x) =
is_nan(x)? "nan" :
is_string(x)? "string" :
is_list(x)? "list" :
"range";
is_range(x) ? "range" :
"invalid";
// Function: is_type()
@ -70,8 +73,8 @@ function is_str(x) = is_string(x);
// is_int(n)
// Description:
// Returns true if the given value is an integer (it is a number and it rounds to itself).
function is_int(n) = is_num(n) && n == round(n);
function is_integer(n) = is_num(n) && n == round(n);
function is_int(n) = is_finite(n) && n == round(n);
function is_integer(n) = is_finite(n) && n == round(n);
// Function: is_nan()
@ -93,7 +96,17 @@ function is_finite(v) = is_num(0*v);
// Function: is_range()
// Description:
// Returns true if its argument is a range
function is_range(x) = is_num(x[0]) && !is_list(x);
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()
@ -106,13 +119,15 @@ function is_range(x) = is_num(x[0]) && !is_list(x);
// is_list_of([3,4,5], 0); // Returns true
// is_list_of([3,4,undef], 0); // Returns false
// is_list_of([[3,4],[4,5]], [1,1]); // Returns true
// is_list_of([[3,"a"],[4,true]], [1,undef]); // Returns true
// is_list_of([[3,4], 6, [4,5]], [1,1]); // Returns false
// is_list_of([[1,[3,4]], [4,[5,6]]], [1,[2,3]]); // Returne true
// is_list_of([[1,[3,INF]], [4,[5,6]]], [1,[2,3]]); // Returne false
// is_list_of([[1,[3,4]], [4,[5,6]]], [1,[2,3]]); // Returns true
// is_list_of([[1,[3,INF]], [4,[5,6]]], [1,[2,3]]); // Returns false
// is_list_of([], [1,[2,3]]); // Returns true
function is_list_of(list,pattern) =
let(pattern = 0*pattern)
is_list(list) &&
[]==[for(entry=list) if (entry*0 != pattern) entry];
[]==[for(entry=0*list) if (entry != pattern) entry];
// Function: is_consistent()
@ -128,7 +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()
@ -139,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.
@ -313,7 +336,8 @@ function scalar_vec3(v, dflt=undef) =
// r = Radius of circle to get the number of segments for.
function segs(r) =
$fn>0? ($fn>3? $fn : 3) :
ceil(max(5, min(360/$fa, abs(r)*2*PI/$fs)));
let( r = is_finite(r)? r: 0 )
ceil(max(5, min(360/$fa, abs(r)*2*PI/$fs))) ;
@ -322,7 +346,7 @@ function segs(r) =
function _valstr(x) =
is_list(x)? str("[",str_join([for (xx=x) _valstr(xx)],","),"]") :
is_num(x)? fmt_float(x,12) : x;
is_finite(x)? fmt_float(x,12) : x;
// Module: assert_approx()

465
math.scad
View file

@ -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,7 +240,10 @@ 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)] :
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;
@ -214,7 +272,10 @@ 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)] :
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;
@ -243,7 +304,10 @@ 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)] :
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;
@ -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,11 +494,10 @@ 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)")
!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,10 +517,9 @@ 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);
function _sum(v,_total,_i=0) = _i>=len(v) ? _total : _sum(v,_total+v[_i], _i+1);
@ -495,10 +570,11 @@ 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(
assert( is_finite(a) && is_matrix(sines,undef,3), "Invalid input.")
sum([ for (s = sines)
let(
ss=point3d(s),
v=ss.x*sin(a*ss.y+ss.z)
v=ss[0]*sin(a*ss[1]+ss[2])
) v
]);
@ -506,26 +582,39 @@ function sum_of_sines(a, sines) =
// 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,7 +724,11 @@ 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()
@ -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,7 +832,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(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) :
@ -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
@ -777,8 +892,10 @@ function is_matrix(A,m,n,square=false) =
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 a<b. Returns >0 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
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];
)
cmps==[]? (len(a)-len(b)) : cmps[0];
// Function: any()
@ -843,14 +961,13 @@ 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()
// Description:
// Returns true if all items in list `l` evaluate as true.
@ -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<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];
// Section: Calculus
@ -921,21 +1052,30 @@ function count_true(l, nmax=undef, i=0, cnt=0) =
// between data[i+1] and data[i], and the data values will be linearly resampled at each corner
// to produce a uniform spacing for the derivative estimate. At the endpoints a single point method
// is used: f'(t) = (f(t+h)-f(t))/h.
// Arguments:
// data = the list of the elements to compute the derivative of.
// h = the parametric sampling of the data.
// closed = boolean to indicate if the data set should be wrapped around from the end to the start.
function deriv(data, h=1, closed=false) =
assert( is_consistent(data) , "Input list is not consistent or not numerical.")
assert( len(data)>=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] :
]
: 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]:
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
@ -947,15 +1087,13 @@ function _dnu_calc(f1,fc,f2,h1,h2) =
f1 = h2<h1 ? lerp(fc,f1,h2/h1) : f1 ,
f2 = h1<h2 ? lerp(fc,f2,h1/h2) : f2
)
(f2-f1) / 2 / min([h1,h2]);
(f2-f1) / 2 / min(h1,h2);
function _deriv_nonuniform(data, h, closed) =
assert(len(h) == len(data)-(closed?0:1),str("Vector valued h must be length ",len(data)-(closed?0:1)))
let(
L = len(data)
)
closed? [for(i=[0:1:L-1])
let( L = len(data) )
closed
? [for(i=[0:1:L-1])
_dnu_calc(data[(L+i-1)%L], data[i], data[(i+1)%L], select(h,i-1), h[i]) ]
: [
(data[1]-data[0])/h[0],
@ -967,15 +1105,23 @@ function _deriv_nonuniform(data, h, closed) =
// Function: deriv2()
// Usage: deriv2(data, [h], [closed])
// Description:
// Computes a numerical esimate of the second derivative of the data, which may be scalar or vector valued.
// 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.
// 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]. For internal points this function uses the approximation
// f''(t) = (f(t-h)-2*f(t)+f(t+h))/h^2. For the endpoints (when closed=false) the algorithm
// when sufficient points are available the method is either the four point expression
// 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) = (f(t-h)-2*f(t)+f(t+h))/h^2. For the endpoints (when closed=false),
// when sufficient points are available, the method is either the four point expression
// f''(t) = (2*f(t) - 5*f(t+h) + 4*f(t+2*h) - f(t+3*h))/h^2 or
// 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
// if five points are available.
// 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 deriv2(data, h=1, closed=false) =
assert( is_consistent(data) , "Input list is not consistent or not numerical.")
assert( len(data)>=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,16 +1185,22 @@ 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
@ -1056,14 +1211,35 @@ 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.
// 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:
@ -1074,7 +1250,9 @@ function polynomial(p, z, k, zk, total) =
// 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")
assert( is_list(p)
&& []==[for(pi=p) if( !is_vector(pi) && pi!=[]) 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)))
:
@ -1087,6 +1265,20 @@ function poly_mult(p,q) =
])
]);
function poly_mult(p,q) =
    is_undef(q) ?
       len(p)==2 ? poly_mult(p[0],p[1])
                 : poly_mult(p[0], poly_mult(select(p,1,-1)))
    :
    assert( is_vector(p) && is_vector(q),"Invalid arguments to poly_mult")
_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:
@ -1096,9 +1288,13 @@ function poly_mult(p,q) =
// 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)] :
function poly_div(n,d,q) =
is_undef(q)
? assert( is_vector(n) && is_vector(d) , "Invalid polynomials." )
let( d = _poly_trim(d) )
assert( d!=[0] , "Denominator cannot be a zero polynomial." )
poly_div(n,d,q=[])
: len(n)<len(d) ? [q,_poly_trim(n)] :
let(
t = n[0] / d[0],
newq = concat(q,[t]),
@ -1114,8 +1310,8 @@ function poly_div(n,d,q=[]) =
// 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);
let( nz = [for(i=[0:1:len(p)-1]) if ( !approx(p[i],0,eps)) i])
len(nz)==0 ? [0] : select(p,nz[0],-1);
// Function: poly_add()
@ -1124,6 +1320,7 @@ function _poly_trim(p,eps=0) =
// Description:
// Computes the sum of two polynomials.
function poly_add(p,q) =
assert( is_vector(p) && is_vector(q), "Invalid input polynomial(s)." )
let( plen = len(p),
qlen = len(q),
long = plen>qlen ? p : q,
@ -1150,11 +1347,10 @@ 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
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]]]
@ -1182,6 +1378,7 @@ function poly_roots(p,tol=1e-14,error_bound=false) =
)
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) =
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))<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

View file

@ -2,10 +2,10 @@ include<../std.scad>
include<../polyhedra.scad>
$fn=96;
if (true) {
$fn=96;
// Display of all solids with insphere, midsphere and circumsphere
for(i=[0:len(_polyhedra_)-1]) {

View file

@ -1,37 +1,14 @@
include <../std.scad>
// List/Array Ops
module test_repeat() {
assert(repeat(1, 4) == [1,1,1,1]);
assert(repeat(8, [2,3]) == [[8,8,8], [8,8,8]]);
assert(repeat(0, [2,2,3]) == [[[0,0,0],[0,0,0]], [[0,0,0],[0,0,0]]]);
assert(repeat([1,2,3],3) == [[1,2,3], [1,2,3], [1,2,3]]);
// Section: List Query Operations
module test_is_simple_list() {
assert(is_simple_list([1,2,3,4]));
assert(is_simple_list([]));
assert(!is_simple_list([1,2,[3,4]]));
}
test_repeat();
module test_in_list() {
assert(in_list("bar", ["foo", "bar", "baz"]));
assert(!in_list("bee", ["foo", "bar", "baz"]));
assert(in_list("bar", [[2,"foo"], [4,"bar"], [3,"baz"]], idx=1));
assert(!in_list(undef, [3,4,5]));
assert(in_list(undef,[3,4,undef,5]));
assert(!in_list(3,[]));
assert(!in_list(3,[4,5,[3]]));
}
test_in_list();
module test_slice() {
assert(slice([3,4,5,6,7,8,9], 3, 5) == [6,7]);
assert(slice([3,4,5,6,7,8,9], 2, -1) == [5,6,7,8,9]);
assert(slice([3,4,5,6,7,8,9], 1, 1) == []);
assert(slice([3,4,5,6,7,8,9], 6, -1) == [9]);
assert(slice([3,4,5,6,7,8,9], 2, -2) == [5,6,7,8]);
}
test_slice();
test_is_simple_list();
module test_select() {
@ -49,6 +26,74 @@ module test_select() {
test_select();
module test_slice() {
assert(slice([3,4,5,6,7,8,9], 3, 5) == [6,7]);
assert(slice([3,4,5,6,7,8,9], 2, -1) == [5,6,7,8,9]);
assert(slice([3,4,5,6,7,8,9], 1, 1) == []);
assert(slice([3,4,5,6,7,8,9], 6, -1) == [9]);
assert(slice([3,4,5,6,7,8,9], 2, -2) == [5,6,7,8]);
assert(slice([], 2, -2) == []);
}
test_slice();
module test_in_list() {
assert(in_list("bar", ["foo", "bar", "baz"]));
assert(!in_list("bee", ["foo", "bar", "baz"]));
assert(in_list("bar", [[2,"foo"], [4,"bar"], [3,"baz"]], idx=1));
assert(!in_list("bee", ["foo", "bar", ["bee"]]));
assert(in_list(NAN, [NAN])==false);
assert(!in_list(undef, [3,4,5]));
assert(in_list(undef,[3,4,undef,5]));
assert(!in_list(3,[]));
assert(!in_list(3,[4,5,[3]]));
}
test_in_list();
module test_min_index() {
assert(min_index([5,3,9,6,2,7,8,2,1])==8);
assert(min_index([5,3,9,6,2,7,8,2,7],all=true)==[4,7]);
// assert(min_index([],all=true)==[]);
}
test_min_index();
module test_max_index() {
assert(max_index([5,3,9,6,2,7,8,9,1])==2);
assert(max_index([5,3,9,6,2,7,8,9,7],all=true)==[2,7]);
// assert(max_index([],all=true)==[]);
}
test_max_index();
module test_list_increasing() {
assert(list_increasing([1,2,3,4]) == true);
assert(list_increasing([1,3,2,4]) == false);
assert(list_increasing([4,3,2,1]) == false);
}
test_list_increasing();
module test_list_decreasing() {
assert(list_decreasing([1,2,3,4]) == false);
assert(list_decreasing([4,2,3,1]) == false);
assert(list_decreasing([4,3,2,1]) == true);
}
test_list_decreasing();
// Section: Basic List Generation
module test_repeat() {
assert(repeat(1, 4) == [1,1,1,1]);
assert(repeat(8, [2,3]) == [[8,8,8], [8,8,8]]);
assert(repeat(0, [2,2,3]) == [[[0,0,0],[0,0,0]], [[0,0,0],[0,0,0]]]);
assert(repeat([1,2,3],3) == [[1,2,3], [1,2,3], [1,2,3]]);
assert(repeat(4, [2,-1]) == [[], []]);
}
test_repeat();
module test_list_range() {
assert(list_range(4) == [0,1,2,3]);
assert(list_range(n=4, step=2) == [0,2,4,6]);
@ -66,6 +111,8 @@ test_list_range();
module test_reverse() {
assert(reverse([3,4,5,6]) == [6,5,4,3]);
assert(reverse("abcd") == ["d","c","b","a"]);
assert(reverse([]) == []);
}
test_reverse();
@ -90,6 +137,8 @@ module test_deduplicate() {
assert(deduplicate(closed=true, [8,3,4,4,4,8,2,3,3,8,8]) == [8,3,4,8,2,3]);
assert(deduplicate("Hello") == ["H","e","l","o"]);
assert(deduplicate([[3,4],[7,1.99],[7,2],[1,4]],eps=0.1) == [[3,4],[7,2],[1,4]]);
assert(deduplicate([], closed=true) == []);
assert(deduplicate([[1,[1,[undef]]],[1,[1,[undef]]],[1,[2]],[1,[2,[0]]]])==[[1, [1,[undef]]],[1,[2]],[1,[2,[0]]]]);
}
test_deduplicate();
@ -148,22 +197,6 @@ module test_list_bset() {
test_list_bset();
module test_list_increasing() {
assert(list_increasing([1,2,3,4]) == true);
assert(list_increasing([1,3,2,4]) == false);
assert(list_increasing([4,3,2,1]) == false);
}
test_list_increasing();
module test_list_decreasing() {
assert(list_decreasing([1,2,3,4]) == false);
assert(list_decreasing([4,2,3,1]) == false);
assert(list_decreasing([4,3,2,1]) == true);
}
test_list_decreasing();
module test_list_shortest() {
assert(list_shortest(["foobar", "bazquxx", "abcd"]) == 4);
}
@ -315,6 +348,13 @@ test_set_intersection();
// Arrays
module test_add_scalar() {
assert(add_scalar([1,2,3],3) == [4,5,6]);
assert(add_scalar([[1,2,3],[3,4,5]],3) == [[4,5,6],[6,7,8]]);
}
test_add_scalar();
module test_subindex() {
v = [[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]];
assert(subindex(v,2) == [3, 7, 11, 15]);
@ -402,10 +442,18 @@ test_array_group();
module test_flatten() {
assert(flatten([[1,2,3], [4,5,[6,7,8]]]) == [1,2,3,4,5,[6,7,8]]);
assert(flatten([]) == []);
}
test_flatten();
module test_full_flatten() {
assert(full_flatten([[1,2,3], [4,5,[6,[7],8]]]) == [1,2,3,4,5,6,7,8]);
assert(full_flatten([]) == []);
}
test_full_flatten();
module test_array_dim() {
assert(array_dim([[[1,2,3],[4,5,6]],[[7,8,9],[10,11,12]]]) == [2,2,3]);
assert(array_dim([[[1,2,3],[4,5,6]],[[7,8,9],[10,11,12]]], 0) == 2);

View file

@ -18,6 +18,10 @@ module test_typeof() {
assert(typeof([0:1:5]) == "range");
assert(typeof([-3:2:5]) == "range");
assert(typeof([10:-2:-10]) == "range");
assert(typeof([0:NAN:INF]) == "invalid");
assert(typeof([0:"a":INF]) == "undef");
assert(typeof([0:[]:INF]) == "undef");
assert(typeof([true:1:INF]) == "undef");
}
test_typeof();
@ -102,6 +106,8 @@ module test_is_int() {
assert(!is_int(-99.1));
assert(!is_int(99.1));
assert(!is_int(undef));
assert(!is_int(INF));
assert(!is_int(NAN));
assert(!is_int(false));
assert(!is_int(true));
assert(!is_int("foo"));
@ -124,6 +130,8 @@ module test_is_integer() {
assert(!is_integer(-99.1));
assert(!is_integer(99.1));
assert(!is_integer(undef));
assert(!is_integer(INF));
assert(!is_integer(NAN));
assert(!is_integer(false));
assert(!is_integer(true));
assert(!is_integer("foo"));
@ -161,16 +169,33 @@ module test_is_range() {
assert(!is_range(5));
assert(!is_range(INF));
assert(!is_range(-INF));
assert(!is_nan(NAN));
assert(!is_range(""));
assert(!is_range("foo"));
assert(!is_range([]));
assert(!is_range([3,4,5]));
assert(!is_range([INF:4:5]));
assert(!is_range([3:NAN:5]));
assert(!is_range([3:4:"a"]));
assert(is_range([3:1:5]));
}
test_is_nan();
test_is_range();
module test_valid_range() {
assert(valid_range([0:0]));
assert(valid_range([0:1:0]));
assert(valid_range([0:1:10]));
assert(valid_range([0.1:1.1:2.1]));
assert(valid_range([0:-1:0]));
assert(valid_range([10:-1:0]));
assert(valid_range([2.1:-1.1:0.1]));
assert(!valid_range([10:1:0]));
assert(!valid_range([2.1:1.1:0.1]));
assert(!valid_range([0:-1:10]));
assert(!valid_range([0.1:-1.1:2.1]));
}
test_valid_range();
module test_is_list_of() {
assert(is_list_of([3,4,5], 0));
assert(!is_list_of([3,4,undef], 0));
@ -181,10 +206,14 @@ module test_is_list_of() {
}
test_is_list_of();
module test_is_consistent() {
assert(is_consistent([]));
assert(is_consistent([[],[]]));
assert(is_consistent([3,4,5]));
assert(is_consistent([[3,4],[4,5],[6,7]]));
assert(is_consistent([[[3],4],[[4],5]]));
assert(!is_consistent(5));
assert(!is_consistent(undef));
assert(!is_consistent([[3,4,5],[3,4]]));
assert(is_consistent([[3,[3,4,[5]]], [5,[2,9,[9]]]]));
assert(!is_consistent([[3,[3,4,[5]]], [5,[2,9,9]]]));
@ -331,11 +360,25 @@ module test_scalar_vec3() {
assert(scalar_vec3([3]) == [3,0,0]);
assert(scalar_vec3([3,4]) == [3,4,0]);
assert(scalar_vec3([3,4],dflt=1) == [3,4,1]);
assert(scalar_vec3([3,"a"],dflt=1) == [3,"a",1]);
assert(scalar_vec3([3,[2]],dflt=1) == [3,[2],1]);
assert(scalar_vec3([3],dflt=1) == [3,1,1]);
assert(scalar_vec3([3,4,5]) == [3,4,5]);
assert(scalar_vec3([3,4,5,6]) == [3,4,5]);
assert(scalar_vec3([3,4,5,6]) == [3,4,5]);
}
test_scalar_vec3();
module test_segs() {
assert_equal(segs(50,$fn=8), 8);
assert_equal(segs(50,$fa=2,$fs=2), 158);
assert(segs(1)==5);
assert(segs(11)==30);
// assert(segs(1/0)==5);
// assert(segs(0/0)==5);
// assert(segs(undef)==5);
}
test_segs();
// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap

View file

@ -110,6 +110,8 @@ module test_approx() {
assert_equal(approx(1/3, 0.3333333333), true);
assert_equal(approx(-1/3, -0.3333333333), true);
assert_equal(approx(10*[cos(30),sin(30)], 10*[sqrt(3)/2, 1/2]), true);
assert_equal(approx([1,[1,undef]], [1+1e-12,[1,true]]), false);
assert_equal(approx([1,[1,undef]], [1+1e-12,[1,undef]]), true);
}
test_approx();
@ -389,7 +391,6 @@ module test_mean() {
}
test_mean();
module test_median() {
assert_equal(median([2,3,7]), 4.5);
assert_equal(median([[1,2,3], [3,4,5], [8,9,10]]), [4.5,5.5,6.5]);
@ -397,6 +398,16 @@ module test_median() {
test_median();
module test_convolve() {
assert_equal(convolve([],[1,2,1]), []);
assert_equal(convolve([1,1],[]), []);
assert_equal(convolve([1,1],[1,2,1]), [1,3,3,1]);
assert_equal(convolve([1,2,3],[1,2,1]), [1,4,8,8,3]);
}
test_convolve();
module test_matrix_inverse() {
assert_approx(matrix_inverse(rot([20,30,40])), [[0.663413948169,0.556670399226,-0.5,0],[-0.47302145844,0.829769465589,0.296198132726,0],[0.579769465589,0.0400087565481,0.813797681349,0],[0,0,0,1]]);
}
@ -583,6 +594,24 @@ module test_factorial() {
}
test_factorial();
module test_binomial() {
assert_equal(binomial(1), [1,1]);
assert_equal(binomial(2), [1,2,1]);
assert_equal(binomial(3), [1,3,3,1]);
assert_equal(binomial(5), [1,5,10,10,5,1]);
}
test_binomial();
module test_binomial_coefficient() {
assert_equal(binomial_coefficient(2,1), 2);
assert_equal(binomial_coefficient(3,2), 3);
assert_equal(binomial_coefficient(4,2), 6);
assert_equal(binomial_coefficient(10,7), 120);
assert_equal(binomial_coefficient(10,7), binomial(10)[7]);
assert_equal(binomial_coefficient(15,4), binomial(15)[4]);
}
test_binomial_coefficient();
module test_gcd() {
assert_equal(gcd(15,25), 5);
@ -682,6 +711,7 @@ 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([1,2],[4,5,6]), [[4,5,6],[8,10,12]]);
assert_equal(outer_product([9],[7]), [[63]]);
}
test_outer_product();
@ -782,8 +812,10 @@ test_deriv3();
module test_polynomial(){
assert_equal(polynomial([],12),0);
assert_equal(polynomial([],[12,4]),[0,0]);
assert_equal(polynomial([0],12),0);
assert_equal(polynomial([0],[12,4]),[0,0]);
// 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);
@ -879,16 +911,20 @@ 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([3,2,1],[0]),[0]);
// 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]]),[]);
assert_equal(poly_mult([[1,2],[0],[5,6]]), [0]);
// assert_equal(poly_mult([[1,2],[],[5,6]]), []);
assert_equal(poly_mult([[3,4,5],[0,0,0]]),[0]);
// 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(poly_mult([4,3,3,2],[2,1,3]), [2,1,3]),[[4,3,3,2],[0]]);
// 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]]);
@ -899,7 +935,8 @@ 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]),[]);
assert_equal(poly_add([1,2,3],-[1,2,3]),[0]);
// assert_equal(poly_add([1,2,3],-[1,2,3]),[]);
}
test_poly_add();

View file

@ -9,6 +9,7 @@ module test_is_vector() {
assert(is_vector(1) == false);
assert(is_vector("foo") == false);
assert(is_vector(true) == false);
assert(is_vector([0,0,0],zero=true) == true);
assert(is_vector([0,0,0],zero=false) == false);
assert(is_vector([0,1,0],zero=true) == false);
@ -17,13 +18,6 @@ module test_is_vector() {
test_is_vector();
module test_add_scalar() {
assert(add_scalar([1,2,3],3) == [4,5,6]);
assert(add_scalar([[1,2,3],[3,4,5]],3) == [[4,5,6],[6,7,8]]);
}
test_add_scalar();
module test_vfloor() {
assert_equal(vfloor([2.0, 3.14, 18.9, 7]), [2,3,18,7]);
assert_equal(vfloor([-2.0, -3.14, -18.9, -7]), [-2,-4,-19,-7]);

View file

@ -19,8 +19,6 @@
// Arguments:
// v = The value to test to see if it is a vector.
// length = If given, make sure the vector is `length` items long.
// zero = If false, require that the length of the vector is not approximately zero. If true, require the length of the vector to be approx zero-length. Default: `undef` (don't check vector length.)
// eps = The minimum vector length that is considered non-zero. Default: `EPSILON` (`1e-9`)
// Example:
// is_vector(4); // Returns false
// is_vector([4,true,false]); // Returns false
@ -30,31 +28,16 @@
// is_vector([3,4,5],3); // Returns true
// is_vector([3,4,5],4); // Returns true
// is_vector([]); // Returns false
// is_vector([0,0,0],zero=true); // Returns true
// is_vector([0,4,0],3,zero=false); // Returns true
// is_vector([0,0,0],zero=false); // Returns false
// is_vector([0,1,0],zero=true); // Returns false
// is_vector([0,0,1],zero=false); // Returns true
// is_vector([0,0,1e-12],zero=false); // Returns false
// is_vector([],zero=false); // Returns false
function is_vector(v,length,zero,eps=EPSILON) =
is_list(v) && is_num(0*(v*v))
&& (is_undef(length) || len(v)==length)
&& (is_undef(zero) || ((norm(v) >= eps) == !zero));
// Function: add_scalar()
// Usage:
// add_scalar(v,s);
// Description:
// Given a vector and a scalar, returns the vector with the scalar added to each item in it.
// If given a list of vectors, recursively adds the scalar to the each vector.
// Arguments:
// v = The initial list of values.
// s = A scalar value to add to every item in the vector.
// Example:
// add_scalar([1,2,3],3); // Returns: [4,5,6]
// add_scalar([[1,2,3],[3,4,5]],3); // Returns: [[4,5,6],[6,7,8]]
function add_scalar(v,s) = [for (x=v) is_list(x)? add_scalar(x,s) : x+s];
// Function: vang()
// Usage:
// theta = vang([X,Y]);
@ -63,6 +46,7 @@ function add_scalar(v,s) = [for (x=v) is_list(x)? add_scalar(x,s) : x+s];
// Given a 2D vector, returns the angle in degrees counter-clockwise from X+ on the XY plane.
// Given a 3D vector, returns [THETA,PHI] where THETA is the number of degrees counter-clockwise from X+ on the XY plane, and PHI is the number of degrees up from the X+ axis along the XZ plane.
function vang(v) =
assert( is_vector(v,2) || is_vector(v,3) , "Invalid vector")
len(v)==2? atan2(v.y,v.x) :
let(res=xyz_to_spherical(v)) [res[1], 90-res[2]];
@ -70,13 +54,19 @@ function vang(v) =
// Function: vmul()
// Description:
// Element-wise vector multiplication. Multiplies each element of vector `v1` by
// the corresponding element of vector `v2`. Returns a vector of the products.
// the corresponding element of vector `v2`. The vectors should have the same dimension.
// Returns a vector of the products.
// Arguments:
// v1 = The first vector.
// v2 = The second vector.
// Example:
// vmul([3,4,5], [8,7,6]); // Returns [24, 28, 30]
function vmul(v1, v2) = [for (i = [0:1:len(v1)-1]) v1[i]*v2[i]];
function vmul(v1, v2) =
// this thighter check can be done yet because it would break other codes in the library
// assert( is_vector(v1) && is_vector(v2,len(v1)), "Incompatible vectors")
assert( is_vector(v1) && is_vector(v2), "Invalid vector(s)")
[for (i = [0:1:len(v1)-1]) v1[i]*v2[i]];
// Function: vdiv()
@ -88,7 +78,9 @@ function vmul(v1, v2) = [for (i = [0:1:len(v1)-1]) v1[i]*v2[i]];
// v2 = The second vector.
// Example:
// vdiv([24,28,30], [8,7,6]); // Returns [3, 4, 5]
function vdiv(v1, v2) = [for (i = [0:1:len(v1)-1]) v1[i]/v2[i]];
function vdiv(v1, v2) =
assert( is_vector(v1) && is_vector(v2,len(v1)), "Incompatible vectors")
[for (i = [0:1:len(v1)-1]) v1[i]/v2[i]];
// Function: vabs()
@ -97,19 +89,25 @@ function vdiv(v1, v2) = [for (i = [0:1:len(v1)-1]) v1[i]/v2[i]];
// v = The vector to get the absolute values of.
// Example:
// vabs([-1,3,-9]); // Returns: [1,3,9]
function vabs(v) = [for (x=v) abs(x)];
function vabs(v) =
assert( is_vector(v), "Invalid vector" )
[for (x=v) abs(x)];
// Function: vfloor()
// Description:
// Returns the given vector after performing a `floor()` on all items.
function vfloor(v) = [for (x=v) floor(x)];
function vfloor(v) =
assert( is_vector(v), "Invalid vector" )
[for (x=v) floor(x)];
// Function: vceil()
// Description:
// Returns the given vector after performing a `ceil()` on all items.
function vceil(v) = [for (x=v) ceil(x)];
function vceil(v) =
assert( is_vector(v), "Invalid vector" )
[for (x=v) ceil(x)];
// Function: unit()
@ -137,6 +135,7 @@ function unit(v, error=[[["ASSERT"]]]) =
// Function: vector_angle()
// Usage:
// vector_angle(v1,v2);
// vector_angle([v1,v2]);
// vector_angle(PT1,PT2,PT3);
// vector_angle([PT1,PT2,PT3]);
// Description:
@ -156,20 +155,21 @@ function unit(v, error=[[["ASSERT"]]]) =
// vector_angle([10,0,10], [0,0,0], [-10,10,0]); // Returns: 120
// vector_angle([[10,0,10], [0,0,0], [-10,10,0]]); // Returns: 120
function vector_angle(v1,v2,v3) =
let(
vecs = !is_undef(v3)? [v1-v2,v3-v2] :
!is_undef(v2)? [v1,v2] :
len(v1) == 3? [v1[0]-v1[1],v1[2]-v1[1]] :
len(v1) == 2? v1 :
assert(false, "Bad arguments to vector_angle()"),
is_valid = is_vector(vecs[0]) && is_vector(vecs[1]) && vecs[0]*0 == vecs[1]*0
assert( ( is_undef(v3) && ( is_undef(v2) || same_shape(v1,v2) ) )
|| is_consistent([v1,v2,v3]) ,
"Bad arguments.")
assert( is_vector(v1) || is_consistent(v1), "Bad arguments.")
let( vecs = ! is_undef(v3) ? [v1-v2,v3-v2] :
! is_undef(v2) ? [v1,v2] :
len(v1) == 3 ? [v1[0]-v1[1], v1[2]-v1[1]]
: v1
)
assert(is_valid, "Bad arguments to vector_angle()")
assert(is_vector(vecs[0],2) || is_vector(vecs[0],3), "Bad arguments.")
let(
norm0 = norm(vecs[0]),
norm1 = norm(vecs[1])
)
assert(norm0>0 && norm1>0,"Zero length vector given to vector_angle()")
assert(norm0>0 && norm1>0, "Zero length vector.")
// NOTE: constrain() corrects crazy FP rounding errors that exceed acos()'s domain.
acos(constrain((vecs[0]*vecs[1])/(norm0*norm1), -1, 1));
@ -177,13 +177,14 @@ function vector_angle(v1,v2,v3) =
// Function: vector_axis()
// Usage:
// vector_axis(v1,v2);
// vector_axis([v1,v2]);
// vector_axis(PT1,PT2,PT3);
// vector_axis([PT1,PT2,PT3]);
// Description:
// If given a single list of two vectors, like `vector_axis([V1,V2])`, returns the vector perpendicular the two vectors V1 and V2.
// If given a single list of three points, like `vector_axis([A,B,C])`, returns the vector perpendicular the line segments AB and BC.
// If given two vectors, like `vector_axis(V1,V1)`, returns the vector perpendicular the two vectors V1 and V2.
// If given three points, like `vector_axis(A,B,C)`, returns the vector perpendicular the line segments AB and BC.
// If given a single list of three points, like `vector_axis([A,B,C])`, returns the vector perpendicular to the plane through a, B and C.
// If given two vectors, like `vector_axis(V1,V2)`, returns the vector perpendicular to the two vectors V1 and V2.
// If given three points, like `vector_axis(A,B,C)`, returns the vector perpendicular to the plane through a, B and C.
// Arguments:
// v1 = First vector or point.
// v2 = Second vector or point.
@ -199,28 +200,23 @@ function vector_axis(v1,v2=undef,v3=undef) =
is_vector(v3)
? assert(is_consistent([v3,v2,v1]), "Bad arguments.")
vector_axis(v1-v2, v3-v2)
:
assert( is_undef(v3), "Bad arguments.")
: assert( is_undef(v3), "Bad arguments.")
is_undef(v2)
? assert( is_list(v1), "Bad arguments.")
len(v1) == 2
? vector_axis(v1[0],v1[1])
: vector_axis(v1[0],v1[1],v1[2])
:
assert(
is_vector(v1,zero=false) &&
is_vector(v2,zero=false) &&
is_consistent([v1,v2]),
"Bad arguments."
)
: assert( is_vector(v1,zero=false) && is_vector(v2,zero=false) && is_consistent([v1,v2])
, "Bad arguments.")
let(
eps = 1e-6,
w1 = point3d(v1/norm(v1)),
w2 = point3d(v2/norm(v2)),
w3 = (norm(w1-w2) > eps && norm(w1+w2) > eps) ? w2
: (norm(vabs(w2)-UP) > eps) ? UP
: (norm(vabs(w2)-UP) > eps)? UP
: RIGHT
) unit(cross(w1,w3));
// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap