mirror of
https://github.com/BelfrySCAD/BOSL2.git
synced 2025-01-01 09:49:45 +00:00
input checks, general review, adding new functions
Defines functions: 1. binomial() 2. convolve() 3. all_numeric() 4. is_addable() Input data check in all functions. Recoding of: 1. cumsum (a non-recursive version) 2. median() 3. count_true() (more efficient version) 4. polynomial functions so that a zero polynomial is represented by [0] and not by []
This commit is contained in:
RonaldoCMP
parent
3d226f1ffa
commit
51be74b5ec
1 changed files with 449 additions and 213 deletions
500
math.scad
500
math.scad
|
|
@ -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,51 @@ 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 +139,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.
|
||||
|
|
@ -114,46 +162,86 @@ function factorial(n,d=0) =
|
|||
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.")
|
||||
assert(is_finite(u) || is_vector(u) || is_range(u), "Input u to lerp must be a number, vector, or range.")
|
||||
[for (v = u) lerp(a,b,v)];
|
||||
|
||||
|
||||
|
||||
// Function: all_numeric()
|
||||
// Usage:
|
||||
// x = all_numeric(list);
|
||||
// Description:
|
||||
// Returns true if all simple value in `list` is a finite number and the list is not empty. Uses recursion.
|
||||
// Arguments:
|
||||
// list = The list to check
|
||||
// Example:
|
||||
// x = all_numeric([1,[2,3,[4]]]); // Returns: true
|
||||
// y = all_numeric([1,[2,3,[1/0]]]); // Returns: false
|
||||
function all_numeric(list) =
|
||||
!is_list(list) ? is_num(0*list)
|
||||
: let(v = list*list[0])
|
||||
is_num(0*(v*v)) // true just for vectors and matrices
|
||||
|| []==[for(vi=list) if( !all_numeric(vi)) 0] ;
|
||||
|
||||
|
||||
// Function: is_addable()
|
||||
// Usage:
|
||||
// x = is_addable(list);
|
||||
// Description:
|
||||
// Returns true if `list` is both consistent and numerical.
|
||||
// Arguments:
|
||||
// list = The list to check
|
||||
// Example:
|
||||
// x = is_addable([[[1],2],[[0],2]])); // Returns: true
|
||||
// y = is_addable([[[1],2],[[0],[]]])); // Returns: false
|
||||
function is_addable(list) = // consistent and numerical
|
||||
is_list_of(list,list[0])
|
||||
&& ( ( let(v = list*list[0]) is_num(0*(v*v)) )
|
||||
|| []==[for(vi=list) if( !all_numeric(vi)) 0] );
|
||||
|
||||
|
||||
|
||||
// Section: Hyperbolic Trigonometry
|
||||
|
||||
// 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;
|
||||
|
||||
|
||||
|
|
@ -185,10 +273,14 @@ 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_int(y), "The multiple must be an 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;
|
||||
|
||||
|
||||
|
||||
// Function: quantdn()
|
||||
// Description:
|
||||
// Quantize a value `x` to an integer multiple of `y`, rounding down to the previous multiple.
|
||||
|
|
@ -214,7 +306,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_int(y), "The multiple must be an 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 +338,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_int(y), "The multiple must be an 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 +362,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 +383,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_int(m), "Input must be finite numbers.")
|
||||
(x%m+m)%m;
|
||||
|
||||
|
||||
// Function: modang(x)
|
||||
|
|
@ -299,6 +401,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 +409,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 +421,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) && is_int(m), "Input must be finite numbers.")
|
||||
let(
|
||||
a = posmod(x, m),
|
||||
b = posmod(y, m),
|
||||
|
|
@ -330,20 +434,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+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 +463,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+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 +480,7 @@ 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+factor+N+seed), "Input must be finite numbers.")
|
||||
assert(maxval >= minval, "maxval cannot be smaller than minval")
|
||||
let(
|
||||
minv = 1-1/pow(factor,minval),
|
||||
|
|
@ -395,18 +502,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(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 +525,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);
|
||||
|
||||
|
||||
|
|
@ -441,10 +547,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);
|
||||
|
|
@ -452,24 +557,30 @@ function _sum(v,_total,_i=0) = _i>=len(v) ? _total : _sum(v,_total+v[_i], _i+1);
|
|||
|
||||
// Function: cumsum()
|
||||
// Description:
|
||||
// Returns a list where each item is the cumulative sum of all items up to and including the corresponding entry in the input list.
|
||||
// If passed an array of vectors, returns a list of cumulative vectors sums.
|
||||
// Returns a list where each item is the cumulative sum of all items up to and including the corresponding entry
|
||||
// in the input list `v`. If passed an array of vectors, returns a list of cumulative vectors sums.
|
||||
// if a compatible value `off` is given, it is added to each item in the output list.
|
||||
// Arguments:
|
||||
// v = The list to get the sum of.
|
||||
// off = An offset to add to each output sum.
|
||||
// Example:
|
||||
// cumsum([1,1,1]); // returns [1,2,3]
|
||||
// cumsum([2,2,2]); // returns [2,4,6]
|
||||
// cumsum([1,2,3]); // returns [1,3,6]
|
||||
// cumsum([[1,2,3], [3,4,5], [5,6,7]]); // returns [[1,2,3], [4,6,8], [9,12,15]]
|
||||
function cumsum(v,_i=0,_acc=[]) =
|
||||
_i==len(v) ? _acc :
|
||||
cumsum(
|
||||
v, _i+1,
|
||||
concat(
|
||||
_acc,
|
||||
[_i==0 ? v[_i] : select(_acc,-1)+v[_i]]
|
||||
)
|
||||
);
|
||||
// cumsum([-2,-1,0,1,2],3); // returns [1,0,0,1,3]
|
||||
function cumsum(v, off) =
|
||||
assert( is_list(v), "The input is not a list.")
|
||||
v==[] ? [] :
|
||||
let( off = is_undef(off)? 0*v[0]: off )
|
||||
assert( is_consistent(v) , "Improper inputs.")
|
||||
[for( i = 0,
|
||||
S = v[0]+off;
|
||||
i<=len(v)-1;
|
||||
i = i+1,
|
||||
S = S + v[i]
|
||||
) S ];
|
||||
|
||||
|
||||
|
||||
// Function: sum_of_squares()
|
||||
|
|
@ -495,10 +606,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), "Invalid inputs.")
|
||||
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,13 +618,17 @@ 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) =
|
||||
(is_vector(v) || is_matrix(v)) && len(v)>1 ? [for (p=pair(v)) p[1]-p[0]] :
|
||||
all_numeric(v) ? [for (p=pair(v)) p[1]-p[0]] :
|
||||
assert( false, "Inconsistent or non numeric input list.");
|
||||
|
||||
|
||||
// Function: product()
|
||||
|
|
@ -525,7 +641,16 @@ function deltas(v) = [for (p=pair(v)) p.y-p.x];
|
|||
// 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_consistent(v), "Inconsistent 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 +659,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 +682,35 @@ 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) && len(v)>0 ? (min(v)+max(v))/2 :
|
||||
is_matrix(v) && len(v)>0
|
||||
? [for(ti=transpose(v)) (min(ti)+max(ti))/2 ]
|
||||
: assert(false , "Invalid input or empty list.");
|
||||
|
||||
// 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 +725,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,6 +762,7 @@ function matrix_inverse(A) =
|
|||
// Description:
|
||||
// Returns a submatrix with the specified index ranges or index sets.
|
||||
function submatrix(M,ind1,ind2) =
|
||||
assert( is_matrix(M), "Input must be a matrix." )
|
||||
[for(i=ind1) [for(j=ind2) M[i][j] ] ];
|
||||
|
||||
|
||||
|
|
@ -628,7 +772,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 +803,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 +838,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 +852,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 +867,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 +900,16 @@ 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_vector(A[0],n)
|
||||
&& is_vector(A*(0*A[0]),m)
|
||||
&& ( !square || len(A)==len(A[0]));
|
||||
|
||||
function is_matrix(A,m,n,square=false) =
|
||||
is_list(A[0])
|
||||
&& ( let(v = A*A[0]) is_num(0*(v*v)) ) // a matrix of finite numbers
|
||||
&& (is_undef(n) || len(A[0])==n )
|
||||
&& (is_undef(m) || len(A)==m )
|
||||
&& ( !square || len(A)==len(A[0]));
|
||||
|
||||
|
||||
// Section: Comparisons and Logic
|
||||
|
|
@ -777,8 +932,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 +953,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 +977,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 +1001,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,15 +1022,14 @@ 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()
|
||||
// Usage:
|
||||
// count_true(l)
|
||||
|
|
@ -904,6 +1060,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 +1092,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_addable(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 +1127,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 +1145,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_addable(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 +1189,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_addable(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 +1225,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)
|
||||
// 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,25 +1251,42 @@ 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: poly_mult()
|
||||
// Usage
|
||||
// polymult(p,q)
|
||||
// polymult([p1,p2,p3,...])
|
||||
// Descriptoin:
|
||||
// Description:
|
||||
// Given a list of polynomials represented as real coefficient lists, with the highest degree coefficient first,
|
||||
// computes the coefficient list of the product polynomial.
|
||||
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 +1299,21 @@ 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 +1323,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]),
|
||||
|
|
@ -1115,7 +1346,7 @@ function poly_div(n,d,q=[]) =
|
|||
// 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);
|
||||
len(nz)==0 ? [0] : select(p,nz[0],-1);
|
||||
|
||||
|
||||
// Function: poly_add()
|
||||
|
|
@ -1124,6 +1355,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 +1382,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 +1413,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 +1454,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
|
||||
|
|
|
|||
Loading…
Reference in a new issue