grey/BOSL2
1
0
Fork 0
mirror of https://github.com/BelfrySCAD/BOSL2.git synced 2024-12-29 16:29:40 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions
BOSL2/math.scad
Adrian Mariano 262ad3e20b fix comment in math.scad
2023-07-25 16:27:47 -04:00

1660 lines
64 KiB
OpenSCAD
Raw Blame History

//////////////////////////////////////////////////////////////////////
// LibFile: math.scad
// Assorted math functions, including linear interpolation, list operations (sums, mean, products),
// convolution, quantization, log2, hyperbolic trig functions, random numbers, derivatives,
// polynomials, and root finding.
// Includes:
// include <BOSL2/std.scad>
// FileGroup: Math
// FileSummary: Math on lists, special functions, quantization, random numbers, calculus, root finding
//
// FileFootnotes: STD=Included in std.scad
//////////////////////////////////////////////////////////////////////
// Section: Math Constants
// Constant: PHI
// Synopsis: The golden ratio φ (phi). Approximately 1.6180339887
// Topics: Constants, Math
// See Also: EPSILON, INF, NAN
// Description: The golden ratio φ (phi). Approximately 1.6180339887
PHI = (1+sqrt(5))/2;
// Constant: EPSILON
// Synopsis: A tiny value to compare floating point values. `1e-9`
// Topics: Constants, Math
// See Also: PHI, EPSILON, INF, NAN
// Description: A really small value useful in comparing floating point numbers. ie: abs(a-b)<EPSILON `1e-9`
EPSILON = 1e-9;
// Constant: INF
// Synopsis: The floating point value for Infinite.
// Topics: Constants, Math
// See Also: PHI, EPSILON, INF, NAN
// Description: The value `inf`, useful for comparisons.
INF = 1/0;
// Constant: NAN
// Synopsis: The floating point value for Not a Number.
// Topics: Constants, Math
// See Also: PHI, EPSILON, INF, NAN
// Description: The value `nan`, useful for comparisons.
NAN = acos(2);
// Section: Interpolation and Counting
// Function: count()
// Synopsis: Creates a list of incrementing numbers.
// Topics: Math, Indexing
// See Also: idx()
// Usage:
// list = count(n, [s], [step], [reverse]);
// Description:
// Creates a list of `n` numbers, starting at `s`, incrementing by `step` each time.
// You can also pass a list for n and then the length of the input list is used.
// Arguments:
// n = The length of the list of numbers to create, or a list to match the length of
// s = The starting value of the list of numbers.
// step = The amount to increment successive numbers in the list.
// reverse = Reverse the list. Default: false.
// Example:
// nl1 = count(5); // Returns: [0,1,2,3,4]
// nl2 = count(5,3); // Returns: [3,4,5,6,7]
// nl3 = count(4,3,2); // Returns: [3,5,7,9]
// nl4 = count(5,reverse=true); // Returns: [4,3,2,1,0]
// nl5 = count(5,3,reverse=true); // Returns: [7,6,5,4,3]
function count(n,s=0,step=1,reverse=false) = let(n=is_list(n) ? len(n) : n)
reverse? [for (i=[n-1:-1:0]) s+i*step]
: [for (i=[0:1:n-1]) s+i*step];
// Function: lerp()
// Synopsis: Linearly interpolates between two values.
// Topics: Interpolation, Math
// See Also: v_lookup(), lerpn()
// Usage:
// x = lerp(a, b, u);
// l = lerp(a, b, LIST);
// Description:
// Interpolate between two values or vectors.
// If `u` is given as a number, returns the single interpolated value.
// 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 an extrapolation
// along the slope formed by `a` and `b`.
// Arguments:
// a = First value or vector.
// b = Second value or vector.
// u = The proportion from `a` to `b` to calculate. Standard range is 0.0 to 1.0, inclusive. If given as a list or range of values, returns a list of results.
// Example:
// x = lerp(0,20,0.3); // Returns: 6
// x = lerp(0,20,0.8); // Returns: 16
// x = lerp(0,20,-0.1); // Returns: -2
// x = lerp(0,20,1.1); // Returns: 22
// p = lerp([0,0],[20,10],0.25); // Returns [5,2.5]
// l = lerp(0,20,[0.4,0.6]); // Returns: [8,12]
// l = lerp(0,20,[0.25:0.25:0.75]); // Returns: [5,10,15]
// Example(2D):
// p1 = [-50,-20]; p2 = [50,30];
// stroke([p1,p2]);
// pts = lerp(p1, p2, [0:1/8:1]);
// // Points colored in ROYGBIV order.
// 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_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 valid range.")
[for (v = u) (1-v)*a + v*b ];
// Function: lerpn()
// Synopsis: Returns exactly `n` values, linearly interpolated between `a` and `b`.
// Topics: Interpolation, Math
// See Also: v_lookup(), lerp()
// Usage:
// x = lerpn(a, b, n);
// x = lerpn(a, b, n, [endpoint]);
// Description:
// Returns exactly `n` values, linearly interpolated between `a` and `b`.
// If `endpoint` is true, then the last value will exactly equal `b`.
// If `endpoint` is false, then the last value will be `a+(b-a)*(1-1/n)`.
// Arguments:
// a = First value or vector.
// b = Second value or vector.
// n = The number of values to return.
// endpoint = If true, the last value will be exactly `b`. If false, the last value will be one step less.
// Example:
// l = lerpn(-4,4,9); // Returns: [-4,-3,-2,-1,0,1,2,3,4]
// l = lerpn(-4,4,8,false); // Returns: [-4,-3,-2,-1,0,1,2,3]
// l = lerpn(0,1,6); // Returns: [0, 0.2, 0.4, 0.6, 0.8, 1]
// l = lerpn(0,1,5,false); // Returns: [0, 0.2, 0.4, 0.6, 0.8]
function lerpn(a,b,n,endpoint=true) =
assert(same_shape(a,b), "Bad or inconsistent inputs to lerpn")
assert(is_int(n))
assert(is_bool(endpoint))
let( d = n - (endpoint? 1 : 0) )
[for (i=[0:1:n-1]) let(u=i/d) (1-u)*a + u*b];
// Section: Miscellaneous Functions
// Function: sqr()
// Synopsis: Returns the square of the given value.
// Topics: Math
// See Also: hypot(), log2()
// Usage:
// x2 = sqr(x);
// Description:
// If given a number, returns the square of that number,
// If given a vector, returns the sum-of-squares/dot product of the vector elements.
// If given a matrix, returns the matrix multiplication of the matrix with itself.
// Example:
// sqr(3); // Returns: 9
// sqr(-4); // Returns: 16
// sqr([2,3,4]); // Returns: 29
// sqr([[1,2],[3,4]]); // Returns [[7,10],[15,22]]
function sqr(x) =
assert(is_finite(x) || is_vector(x) || is_matrix(x), "Input is not a number nor a list of numbers.")
x*x;
// Function: log2()
// Synopsis: Returns the log base 2 of the given value.
// Topics: Math
// See Also: hypot(), sqr()
// Usage:
// val = log2(x);
// Description:
// Returns the logarithm base 2 of the value given.
// Example:
// log2(0.125); // Returns: -3
// log2(16); // Returns: 4
// log2(256); // Returns: 8
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()
// Synopsis: Returns the hypotenuse length of a 2D or 3D triangle.
// Topics: Math
// See Also: hypot(), sqr(), log2()
// Usage:
// l = hypot(x, y, [z]);
// Description:
// Calculate hypotenuse length of a 2D or 3D triangle.
// Arguments:
// x = Length on the X axis.
// y = Length on the Y axis.
// z = Length on the Z axis. Optional.
// Example:
// l = hypot(3,4); // Returns: 5
// l = hypot(3,4,5); // Returns: ~7.0710678119
function hypot(x,y,z=0) =
assert( is_vector([x,y,z]), "Improper number(s).")
norm([x,y,z]);
// Function: factorial()
// Synopsis: Returns the factorial of the given integer.
// Topics: Math
// See Also: hypot(), sqr(), log2()
// Usage:
// x = factorial(n, [d]);
// Description:
// Returns the factorial of the given integer value, or n!/d! if d is given.
// Arguments:
// n = The integer number to get the factorial of. (n!)
// d = If given, the returned value will be (n! / d!)
// Example:
// x = factorial(4); // Returns: 24
// y = factorial(6); // Returns: 720
// z = factorial(9); // Returns: 362880
function factorial(n,d=0) =
assert(is_int(n) && is_int(d) && n>=0 && d>=0, "Factorial is defined only for non negative integers")
assert(d<=n, "d cannot be larger than n")
product([1,for (i=[n:-1:d+1]) i]);
// Function: binomial()
// Synopsis: Returns the binomial coefficients of the integer `n`.
// Topics: Math
// See Also: hypot(), sqr(), log2(), factorial()
// 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()
// Synopsis: Returns the `k`-th binomial coefficient of the integer `n`.
// Topics: Math
// See Also: hypot(), sqr(), log2(), factorial()
// 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: gcd()
// Synopsis: Returns the Greatest Common Divisor/Factor of two integers.
// Topics: Math
// See Also: hypot(), sqr(), log2(), factorial(), binomial(), gcd(), lcm()
// Usage:
// x = gcd(a,b)
// Description:
// Computes the Greatest Common Divisor/Factor of `a` and `b`.
function gcd(a,b) =
assert(is_int(a) && is_int(b),"Arguments to gcd must be integers")
b==0 ? abs(a) : gcd(b,a % b);
// Computes lcm for two integers
function _lcm(a,b) =
assert(is_int(a) && is_int(b), "Invalid non-integer parameters to lcm")
assert(a!=0 && b!=0, "Arguments to lcm should not be zero")
abs(a*b) / gcd(a,b);
// Computes lcm for a list of values
function _lcmlist(a) =
len(a)==1 ? a[0] :
_lcmlist(concat(lcm(a[0],a[1]),list_tail(a,2)));
// Function: lcm()
// Synopsis: Returns the Least Common Multiple of two or more integers.
// Topics: Math
// See Also: hypot(), sqr(), log2(), factorial(), binomial(), gcd(), lcm()
// Usage:
// div = lcm(a, b);
// divs = lcm(list);
// Description:
// Computes the Least Common Multiple of the two arguments or a list of arguments. Inputs should
// 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)")
_lcmlist(arglist);
// Function rational_approx()
// Usage:
// pq = rational_approx(x, maxq);
// Description:
// Finds the best rational approximation p/q to the number x so that q<=maxq. Returns
// the result as `[p,q]`. If the input is zero, then returns `[0,1]`.
// Example:
// pq1 = rational_approx(PI,10); // Returns: [22,7]
// pq2 = rational_approx(PI,10000); // Returns: [355, 113]
// pq3 = rational_approx(221/323,500); // Returns: [13,19]
// pq4 = rational_approx(0,50); // Returns: [0,1]
function rational_approx(x, maxq, cfrac=[], p, q) =
let(
next = floor(x),
fracpart = x-next,
cfrac = [each cfrac, next],
pq = _cfrac_to_pq(cfrac)
)
approx(fracpart,0) ? pq
: pq[1]>maxq ? [p,q]
: rational_approx(1/fracpart,maxq,cfrac, pq[0], pq[1]);
// Converts a continued fraction given as a list with leading integer term
// into a fraction in the form p / q, returning [p,q].
function _cfrac_to_pq(cfrac,p=0,q=1,ind) =
is_undef(ind) ? _cfrac_to_pq(cfrac,p,q,len(cfrac)-1)
: ind==0 ? [p+q*cfrac[0], q]
: _cfrac_to_pq(cfrac, q, cfrac[ind]*q+p, ind-1);
// Section: Hyperbolic Trigonometry
// Function: sinh()
// Synopsis: Returns the hyperbolic sine of the given value.
// Topics: Math, Trigonometry
// See Also: sinh(), cosh(), tanh(), asinh(), acosh(), atanh()
// Usage:
// a = sinh(x);
// 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()
// Synopsis: Returns the hyperbolic cosine of the given value.
// Topics: Math, Trigonometry
// See Also: sinh(), cosh(), tanh(), asinh(), acosh(), atanh()
// Usage:
// a = cosh(x);
// 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()
// Synopsis: Returns the hyperbolic tangent of the given value.
// Topics: Math, Trigonometry
// See Also: sinh(), cosh(), tanh(), asinh(), acosh(), atanh()
// Usage:
// a = tanh(x);
// 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()
// Synopsis: Returns the hyperbolic arc-sine of the given value.
// Topics: Math, Trigonometry
// See Also: sinh(), cosh(), tanh(), asinh(), acosh(), atanh()
// Usage:
// a = asinh(x);
// 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()
// Synopsis: Returns the hyperbolic arc-cosine of the given value.
// Topics: Math, Trigonometry
// See Also: sinh(), cosh(), tanh(), asinh(), acosh(), atanh()
// Usage:
// a = acosh(x);
// 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()
// Synopsis: Returns the hyperbolic arc-tangent of the given value.
// Topics: Math, Trigonometry
// See Also: sinh(), cosh(), tanh(), asinh(), acosh(), atanh()
// Usage:
// a = atanh(x);
// 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()
// Synopsis: Returns `x` quantized to the nearest integer multiple of `y`.
// Topics: Math, Quantization
// See Also: quant(), quantdn(), quantup()
// Usage:
// num = quant(x, y);
// Description:
// Quantize a value `x` to an integer multiple of `y`, rounding to the nearest multiple.
// The value of `y` does NOT have to be an integer. If `x` is a list, then every item
// in that list will be recursively quantized.
// Arguments:
// x = The value or list to quantize.
// y = Positive quantum to quantize to
// Example:
// a = quant(12,4); // Returns: 12
// b = quant(13,4); // Returns: 12
// c = quant(13.1,4); // Returns: 12
// d = quant(14,4); // Returns: 16
// e = quant(14.1,4); // Returns: 16
// f = quant(15,4); // Returns: 16
// g = quant(16,4); // Returns: 16
// h = quant(9,3); // Returns: 9
// i = quant(10,3); // Returns: 9
// j = quant(10.4,3); // Returns: 9
// k = quant(10.5,3); // Returns: 12
// l = quant(11,3); // Returns: 12
// m = quant(12,3); // Returns: 12
// n = quant(11,2.5); // Returns: 10
// o = quant(12,2.5); // Returns: 12.5
// p = quant([12,13,13.1,14,14.1,15,16],4); // Returns: [12,12,12,16,16,16,16]
// q = quant([9,10,10.4,10.5,11,12],3); // Returns: [9,9,9,12,12,12]
// r = quant([[9,10,10.4],[10.5,11,12]],3); // Returns: [[9,9,9],[12,12,12]]
function quant(x,y) =
assert( is_finite(y) && y>0, "The quantum `y` must be a positive value.")
is_num(x) ? round(x/y)*y
: _roundall(x/y)*y;
function _roundall(data) =
[for(x=data) is_list(x) ? _roundall(x) : round(x)];
// Function: quantdn()
// Synopsis: Returns `x` quantized down to an integer multiple of `y`.
// Topics: Math, Quantization
// See Also: quant(), quantdn(), quantup()
// Usage:
// num = quantdn(x, y);
// Description:
// Quantize a value `x` to an integer multiple of `y`, rounding down to the previous multiple.
// The value of `y` does NOT have to be an integer. If `x` is a list, then every item in that
// list will be recursively quantized down.
// Arguments:
// x = The value or list to quantize.
// y = Postive quantum to quantize to.
// Example:
// a = quantdn(12,4); // Returns: 12
// b = quantdn(13,4); // Returns: 12
// c = quantdn(13.1,4); // Returns: 12
// d = quantdn(14,4); // Returns: 12
// e = quantdn(14.1,4); // Returns: 12
// f = quantdn(15,4); // Returns: 12
// g = quantdn(16,4); // Returns: 16
// h = quantdn(9,3); // Returns: 9
// i = quantdn(10,3); // Returns: 9
// j = quantdn(10.4,3); // Returns: 9
// k = quantdn(10.5,3); // Returns: 9
// l = quantdn(11,3); // Returns: 9
// m = quantdn(12,3); // Returns: 12
// n = quantdn(11,2.5); // Returns: 10
// o = quantdn(12,2.5); // Returns: 10
// p = quantdn([12,13,13.1,14,14.1,15,16],4); // Returns: [12,12,12,12,12,12,16]
// q = quantdn([9,10,10.4,10.5,11,12],3); // Returns: [9,9,9,9,9,12]
// r = quantdn([[9,10,10.4],[10.5,11,12]],3); // Returns: [[9,9,9],[9,9,12]]
function quantdn(x,y) =
assert( is_finite(y) && y>0, "The quantum `y` must be a positive value.")
is_num(x) ? floor(x/y)*y
: _floorall(x/y)*y;
function _floorall(data) =
[for(x=data) is_list(x) ? _floorall(x) : floor(x)];
// Function: quantup()
// Synopsis: Returns `x` quantized uo to an integer multiple of `y`.
// Topics: Math, Quantization
// See Also: quant(), quantdn(), quantup()
// Usage:
// num = quantup(x, y);
// Description:
// Quantize a value `x` to an integer multiple of `y`, rounding up to the next multiple.
// The value of `y` does NOT have to be an integer. If `x` is a list, then every item in
// that list will be recursively quantized up.
// Arguments:
// x = The value or list to quantize.
// y = Positive quantum to quantize to.
// Example:
// a = quantup(12,4); // Returns: 12
// b = quantup(13,4); // Returns: 16
// c = quantup(13.1,4); // Returns: 16
// d = quantup(14,4); // Returns: 16
// e = quantup(14.1,4); // Returns: 16
// f = quantup(15,4); // Returns: 16
// g = quantup(16,4); // Returns: 16
// h = quantup(9,3); // Returns: 9
// i = quantup(10,3); // Returns: 12
// j = quantup(10.4,3); // Returns: 12
// k = quantup(10.5,3); // Returns: 12
// l = quantup(11,3); // Returns: 12
// m = quantup(12,3); // Returns: 12
// n = quantdn(11,2.5); // Returns: 12.5
// o = quantdn(12,2.5); // Returns: 12.5
// p = quantup([12,13,13.1,14,14.1,15,16],4); // Returns: [12,16,16,16,16,16,16]
// q = quantup([9,10,10.4,10.5,11,12],3); // Returns: [9,12,12,12,12,12]
// r = quantup([[9,10,10.4],[10.5,11,12]],3); // Returns: [[9,12,12],[12,12,12]]
function quantup(x,y) =
assert( is_finite(y) && y>0, "The quantum `y` must be a positive value.")
is_num(x) ? ceil(x/y)*y
: _ceilall(x/y)*y;
function _ceilall(data) =
[for(x=data) is_list(x) ? _ceilall(x) : ceil(x)];
// Section: Constraints and Modulos
// Function: constrain()
// Synopsis: Returns a value constrained between `minval` and `maxval`, inclusive.
// Topics: Math
// See Also: posmod(), modang()
// Usage:
// val = constrain(v, minval, maxval);
// Description:
// Constrains value to a range of values between minval and maxval, inclusive.
// Arguments:
// v = value to constrain.
// minval = minimum value to return, if out of range.
// maxval = maximum value to return, if out of range.
// Example:
// a = constrain(-5, -1, 1); // Returns: -1
// b = constrain(5, -1, 1); // Returns: 1
// c = constrain(0.3, -1, 1); // Returns: 0.3
// d = constrain(9.1, 0, 9); // Returns: 9
// e = constrain(-0.1, 0, 9); // Returns: 0
function constrain(v, minval, maxval) =
assert( is_finite(v+minval+maxval), "Input must be finite number(s).")
min(maxval, max(minval, v));
// Function: posmod()
// Synopsis: Returns the positive modulo of a value.
// Topics: Math
// See Also: constrain(), posmod(), modang()
// Usage:
// mod = posmod(x, m)
// Description:
// Returns the positive modulo `m` of `x`. Value returned will be in the range 0 ... `m`-1.
// Arguments:
// x = The value to constrain.
// m = Modulo value.
// Example:
// a = posmod(-700,360); // Returns: 340
// b = posmod(-270,360); // Returns: 90
// c = posmod(-120,360); // Returns: 240
// d = posmod(120,360); // Returns: 120
// e = posmod(270,360); // Returns: 270
// f = posmod(700,360); // Returns: 340
// g = posmod(3,2.5); // Returns: 0.5
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()
// Synopsis: Returns an angle normalized to between -180º and 180º.
// Topics: Math
// See Also: constrain(), posmod(), modang()
// Usage:
// ang = modang(x);
// Description:
// Takes an angle in degrees and normalizes it to an equivalent angle value between -180 and 180.
// Example:
// a1 = modang(-700,360); // Returns: 20
// a2 = modang(-270,360); // Returns: 90
// a3 = modang(-120,360); // Returns: -120
// a4 = modang(120,360); // Returns: 120
// a5 = modang(270,360); // Returns: -90
// a6 = 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;
// Section: Operations on Lists (Sums, Mean, Products)
// Function: sum()
// Synopsis: Returns the sum of a list of values.
// Topics: Math
// See Also: mean(), median(), product(), cumsum()
// Usage:
// x = sum(v, [dflt]);
// Description:
// 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.
// dflt = The default value to return if `v` is an empty list. Default: 0
// Example:
// 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) =
v==[]? dflt :
assert(is_consistent(v), "Input to sum is non-numeric or inconsistent")
is_finite(v[0]) || is_vector(v[0]) ? [for(i=v) 1]*v :
_sum(v,v[0]*0);
function _sum(v,_total,_i=0) = _i>=len(v) ? _total : _sum(v,_total+v[_i], _i+1);
// Function: mean()
// Synopsis: Returns the mean value of a list of values.
// Topics: Math, Statistics
// See Also: sum(), mean(), median(), product()
// Usage:
// x = mean(v);
// Description:
// 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) =
assert(is_list(v) && len(v)>0, "Invalid list.")
sum(v)/len(v);
// Function: median()
// Synopsis: Returns the median value of a list of values.
// Topics: Math, Statistics
// See Also: sum(), mean(), median(), product()
// Usage:
// middle = median(v)
// Description:
// Returns the median of the given vector.
function median(v) =
assert(is_vector(v), "Input to median must be a vector")
len(v)%2 ? max( list_smallest(v, ceil(len(v)/2)) ) :
let( lowest = list_smallest(v, len(v)/2 + 1),
max = max(lowest),
imax = search(max,lowest,1),
max2 = max([for(i=idx(lowest)) if(i!=imax[0]) lowest[i] ])
)
(max+max2)/2;
// Function: deltas()
// Synopsis: Returns the deltas between a list of values.
// Topics: Math, Statistics
// See Also: sum(), mean(), median(), product()
// Usage:
// delts = deltas(v,[wrap]);
// Description:
// Returns a list with the deltas of adjacent entries in the given list, optionally wrapping back to the front.
// 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.
// wrap = If true, wrap back to the start from the end. ie: return the difference between the last and first items as the last delta. Default: false
// 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, wrap=false) =
assert( is_consistent(v) && len(v)>1 , "Inconsistent list or with length<=1.")
[for (p=pair(v,wrap)) p[1]-p[0]] ;
// Function: cumsum()
// Synopsis: Returns the running cumulative sum of a list of values.
// Topics: Math, Statistics
// See Also: sum(), mean(), median(), product()
// Usage:
// sums = cumsum(v);
// 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.
// Arguments:
// v = The list to get the sum of.
// 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) =
assert(is_consistent(v), "The input is not consistent." )
len(v)<=1 ? v :
_cumsum(v,_i=1,_acc=[v[0]]);
function _cumsum(v,_i=0,_acc=[]) =
_i>=len(v) ? _acc :
_cumsum( v, _i+1, [ each _acc, _acc[len(_acc)-1] + v[_i] ] );
// Function: product()
// Synopsis: Returns the multiplicative product of a list of values.
// Topics: Math, Statistics
// See Also: sum(), mean(), median(), product(), cumsum()
// Usage:
// x = product(v);
// Description:
// Returns the product of all entries in the given list.
// 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 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) =
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])? v_mul(_tot,v[i]) : _tot*v[i] ) );
// Function: cumprod()
// Synopsis: Returns the running cumulative product of a list of values.
// Topics: Math, Statistics
// See Also: sum(), mean(), median(), product(), cumsum()
// Description:
// Returns a list where each item is the cumulative product of all items up to and including the corresponding entry in the input list.
// If passed an array of vectors, returns a list of elementwise vector products. If passed a list of square matrices returns matrix
// products multiplying on the left, so a list `[A,B,C]` will produce the output `[A,BA,CBA]`.
// Arguments:
// list = The list to get the product of.
// Example:
// cumprod([1,3,5]); // returns [1,3,15]
// cumprod([2,2,2]); // returns [2,4,8]
// cumprod([[1,2,3], [3,4,5], [5,6,7]])); // returns [[1, 2, 3], [3, 8, 15], [15, 48, 105]]
function cumprod(list) =
is_vector(list) ? _cumprod(list) :
assert(is_consistent(list), "Input must be a consistent list of scalars, vectors or square matrices")
is_matrix(list[0]) ? assert(len(list[0])==len(list[0][0]), "Matrices must be square") _cumprod(list)
: _cumprod_vec(list);
function _cumprod(v,_i=0,_acc=[]) =
_i==len(v) ? _acc :
_cumprod(
v, _i+1,
concat(
_acc,
[_i==0 ? v[_i] : v[_i]*_acc[len(_acc)-1]]
)
);
function _cumprod_vec(v,_i=0,_acc=[]) =
_i==len(v) ? _acc :
_cumprod_vec(
v, _i+1,
concat(
_acc,
[_i==0 ? v[_i] : v_mul(_acc[len(_acc)-1],v[_i])]
)
);
// Function: convolve()
// Synopsis: Returns the convolution of `p` and `q`.
// Topics: Math, Statistics
// See Also: sum(), mean(), median(), product(), cumsum()
// Usage:
// x = convolve(p,q);
// Description:
// Given two vectors, or one vector and a path or
// two paths of the same dimension, finds the convolution of them.
// If both parameter are vectors, returns the vector convolution.
// If one parameter is a vector and the other a path,
// convolves using products by scalars and returns a path.
// If both parameters are paths, convolve using scalar products
// and returns a vector.
// The returned vector or path has length len(p)+len(q)-1.
// Arguments:
// p = The first vector or path.
// q = The second vector or path.
// 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]
// c = convolve([[1,1],[2,2],[3,1]],[1,2,1])); // Returns: [[1,1],[4,4],[8,6],[8,4],[3,1]]
// d = convolve([[1,1],[2,2],[3,1]],[[1,2],[2,1]])); // Returns: [3,9,11,7]
function convolve(p,q) =
p==[] || q==[] ? [] :
assert( (is_vector(p) || is_matrix(p))
&& ( is_vector(q) || (is_matrix(q) && ( !is_vector(p[0]) || (len(p[0])==len(q[0])) ) ) ) ,
"The inputs should be vectors or paths all of the same dimension.")
let( n = len(p),
m = len(q))
[for(i=[0:n+m-2], k1 = max(0,i-n+1), k2 = min(i,m-1) )
sum([for(j=[k1:k2]) p[i-j]*q[j] ])
];
// Function: sum_of_sines()
// Synopsis: Returns the sum of one or more sine waves at a given angle.
// Topics: Math, Statistics
// See Also: sum(), mean(), median(), product(), cumsum()
// Usage:
// sum_of_sines(a,sines)
// Description:
// Given a list of sine waves, returns the sum of the sines at the given angle.
// Each sine wave is given as an `[AMPLITUDE, FREQUENCY, PHASE_ANGLE]` triplet.
// - `AMPLITUDE` is the height of the sine wave above (and below) `0`.
// - `FREQUENCY` is the number of times the sine wave repeats in 360º.
// - `PHASE_ANGLE` is the offset in degrees of the sine wave.
// Arguments:
// a = Angle to get the value for.
// sines = List of [amplitude, frequency, phase_angle] items, where the frequency is the number of times the cycle repeats around the circle.
// Example:
// v = sum_of_sines(30, [[10,3,0], [5,5.5,60]]);
function sum_of_sines(a, sines) =
assert( is_finite(a) && is_matrix(sines,undef,3), "Invalid input.")
sum([ for (s = sines)
let(
ss=point3d(s),
v=ss[0]*sin(a*ss[1]+ss[2])
) v
]);
// Section: Random Number Generation
// Function: rand_int()
// Synopsis: Returns a random integer.
// Topics: Random
// See Also: rand_int(), random_points(), gaussian_rands(), random_polygon(), spherical_random_points(), exponential_rands()
// Usage:
// rand_int(minval, maxval, n, [seed]);
// Description:
// Return a list of random integers in the range of minval to maxval, inclusive.
// Arguments:
// 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(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)];
// Function: random_points()
// Synopsis: Returns a list of random points.
// Topics: Random, Points
// See Also: rand_int(), random_points(), random_polygon(), spherical_random_points()
// Usage:
// points = random_points(n, dim, [scale], [seed]);
// Description:
// Generate `n` uniform random points of dimension `dim` with data ranging from -scale to +scale.
// The `scale` may be a number, in which case the random data lies in a cube,
// or a vector with dimension `dim`, in which case each dimension has its own scale.
// Arguments:
// n = number of points to generate. Default: 1
// dim = dimension of the points. Default: 2
// scale = the scale of the point coordinates. Default: 1
// seed = an optional seed for the random generation.
function random_points(n, dim, scale=1, seed) =
assert( is_int(n) && n>=0, "The number of points should be a non-negative integer.")
assert( is_int(dim) && dim>=1, "The point dimensions should be an integer greater than 1.")
assert( is_finite(scale) || is_vector(scale,dim), "The scale should be a number or a vector with length equal to d.")
let(
rnds = is_undef(seed)
? rands(-1,1,n*dim)
: rands(-1,1,n*dim, seed) )
is_num(scale)
? scale*[for(i=[0:1:n-1]) [for(j=[0:dim-1]) rnds[i*dim+j] ] ]
: [for(i=[0:1:n-1]) [for(j=[0:dim-1]) scale[j]*rnds[i*dim+j] ] ];
// Function: gaussian_rands()
// Synopsis: Returns a list of random numbers with a gaussian distribution.
// Topics: Random, Statistics
// See Also: rand_int(), random_points(), gaussian_rands(), random_polygon(), spherical_random_points(), exponential_rands()
// Usage:
// arr = gaussian_rands([n],[mean], [cov], [seed]);
// Description:
// Returns a random number or vector with a Gaussian/normal distribution.
// Arguments:
// n = the number of points to return. Default: 1
// mean = The average of the random value (a number or vector). Default: 0
// cov = covariance matrix of the random numbers, or variance in the 1D case. Default: 1
// seed = If given, sets the random number seed.
function gaussian_rands(n=1, mean=0, cov=1, seed=undef) =
assert(is_num(mean) || is_vector(mean))
let(
dim = is_num(mean) ? 1 : len(mean)
)
assert((dim==1 && is_num(cov)) || is_matrix(cov,dim,dim),"mean and covariance matrix not compatible")
assert(is_undef(seed) || is_finite(seed))
let(
nums = is_undef(seed)? rands(0,1,dim*n*2) : rands(0,1,dim*n*2,seed),
rdata = [for (i = count(dim*n,0,2)) sqrt(-2*ln(nums[i]))*cos(360*nums[i+1])]
)
dim==1 ? add_scalar(sqrt(cov)*rdata,mean) :
assert(is_matrix_symmetric(cov),"Supplied covariance matrix is not symmetric")
let(
L = cholesky(cov)
)
assert(is_def(L), "Supplied covariance matrix is not positive definite")
move(mean,list_to_matrix(rdata,dim)*transpose(L));
// Function: exponential_rands()
// Synopsis: Returns a list of random numbers with an exponential distribution.
// Topics: Random, Statistics
// See Also: rand_int(), random_points(), gaussian_rands(), random_polygon(), spherical_random_points()
// Usage:
// arr = exponential_rands([n], [lambda], [seed])
// Description:
// Returns random numbers with an exponential distribution with parameter lambda, and hence mean 1/lambda.
// Arguments:
// n = number of points to return. Default: 1
// lambda = distribution parameter. The mean will be 1/lambda. Default: 1
function exponential_rands(n=1, lambda=1, seed) =
assert( is_int(n) && n>=1, "The number of points should be an integer greater than zero.")
assert( is_num(lambda) && lambda>0, "The lambda parameter must be a positive number.")
let(
unif = is_def(seed) ? rands(0,1,n,seed=seed) : rands(0,1,n)
)
-(1/lambda) * [for(x=unif) x==1 ? 708.3964185322641 : ln(1-x)]; // Use ln(min_float) when x is 1
// Function: spherical_random_points()
// Synopsis: Returns a list of random points on the surface of a sphere.
// Topics: Random, Points
// See Also: rand_int(), random_points(), gaussian_rands(), random_polygon(), spherical_random_points()
// Usage:
// points = spherical_random_points([n], [radius], [seed]);
// Description:
// Generate `n` 3D uniformly distributed random points lying on a sphere centered at the origin with radius equal to `radius`.
// Arguments:
// n = number of points to generate. Default: 1
// radius = the sphere radius. Default: 1
// seed = an optional seed for the random generation.
// See https://mathworld.wolfram.com/SpherePointPicking.html
function spherical_random_points(n=1, radius=1, seed) =
assert( is_int(n) && n>=1, "The number of points should be an integer greater than zero.")
assert( is_num(radius) && radius>0, "The radius should be a non-negative number.")
let( theta = is_undef(seed)
? rands(0,360,n)
: rands(0,360,n, seed),
cosphi = rands(-1,1,n))
[for(i=[0:1:n-1]) let(
sin_phi=sqrt(1-cosphi[i]*cosphi[i])
)
radius*[sin_phi*cos(theta[i]),sin_phi*sin(theta[i]), cosphi[i]]];
// Function: random_polygon()
// Synopsis: Returns the CCW path of a simple random polygon.
// Topics: Random, Polygon
// See Also: random_points(), spherical_random_points()
// Usage:
// points = random_polygon([n], [size], [seed]);
// Description:
// Generate the `n` vertices of a random counter-clockwise simple 2d polygon
// inside a circle centered at the origin with radius `size`.
// Arguments:
// n = number of vertices of the polygon. Default: 3
// size = the radius of a circle centered at the origin containing the polygon. Default: 1
// seed = an optional seed for the random generation.
function random_polygon(n=3,size=1, seed) =
assert( is_int(n) && n>2, "Improper number of polygon vertices.")
assert( is_num(size) && size>0, "Improper size.")
let(
seed = is_undef(seed) ? rands(0,1,1)[0] : seed,
cumm = cumsum(rands(0.1,10,n+1,seed)),
angs = 360*cumm/cumm[n-1],
rads = rands(.01,size,n,seed+cumm[0])
)
[for(i=count(n)) rads[i]*[cos(angs[i]), sin(angs[i])] ];
// Section: Calculus
// Function: deriv()
// Synopsis: Returns the first derivative estimate of a list of data.
// Topics: Math, Calculus
// See Also: deriv(), deriv2(), deriv3()
// Usage:
// x = deriv(data, [h], [closed])
// Description:
// Computes a numerical 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 symetric derivative approximation
// for internal points, f'(t) = (f(t+h)-f(t-h))/2h. For the endpoints (when closed=false) the algorithm
// uses a two point method if sufficient points are available: f'(t) = (3*(f(t+h)-f(t)) - (f(t+2*h)-f(t+h)))/2h.
// .
// If `h` is a vector then it is assumed to be nonuniform, with h[i] giving the sampling distance
// 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
? [
for(i=[0:1:L-1])
(data[(i+1)%L]-data[(L+i-1)%L])/2/h
]
: let(
first = L<3 ? data[1]-data[0] :
3*(data[1]-data[0]) - (data[2]-data[1]),
last = L<3 ? data[L-1]-data[L-2]:
(data[L-3]-data[L-2])-3*(data[L-2]-data[L-1])
)
[
first/2/h,
for(i=[1:1:L-2]) (data[i+1]-data[i-1])/2/h,
last/2/h
];
function _dnu_calc(f1,fc,f2,h1,h2) =
let(
f1 = h2<h1 ? lerp(fc,f1,h2/h1) : f1 ,
f2 = h1<h2 ? lerp(fc,f2,h1/h2) : f2
)
(f2-f1) / 2 / min(h1,h2);
function _deriv_nonuniform(data, h, closed) =
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],
for(i=[1:1:L-2]) _dnu_calc(data[i-1],data[i],data[i+1], h[i-1],h[i]),
(data[L-1]-data[L-2])/h[L-2]
];
// Function: deriv2()
// Synopsis: Returns the second derivative estimate of a list of data.
// Topics: Math, Calculus
// See Also: deriv(), deriv2(), deriv3()
// Usage:
// x = deriv2(data, [h], [closed])
// Description:
// 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),
// 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( is_finite(h), "The sampling `h` must be a number." )
let( L = len(data) )
assert( L>=3, "Input list has less than 3 elements.")
closed
? [
for(i=[0:1:L-1])
(data[(i+1)%L]-2*data[i]+data[(L+i-1)%L])/h/h
]
:
let(
first =
L==3? data[0] - 2*data[1] + data[2] :
L==4? 2*data[0] - 5*data[1] + 4*data[2] - data[3] :
(35*data[0] - 104*data[1] + 114*data[2] - 56*data[3] + 11*data[4])/12,
last =
L==3? data[L-1] - 2*data[L-2] + data[L-3] :
L==4? -2*data[L-1] + 5*data[L-2] - 4*data[L-3] + data[L-4] :
(35*data[L-1] - 104*data[L-2] + 114*data[L-3] - 56*data[L-4] + 11*data[L-5])/12
) [
first/h/h,
for(i=[1:1:L-2]) (data[i+1]-2*data[i]+data[i-1])/h/h,
last/h/h
];
// Function: deriv3()
// Synopsis: Returns the third derivative estimate of a list of data.
// Topics: Math, Calculus
// See Also: deriv(), deriv2(), deriv3()
// Usage:
// x = deriv3(data, [h], [closed])
// Description:
// 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 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.
// 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 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
)
closed? [
for(i=[0:1:L-1])
(-data[(L+i-2)%L]+2*data[(L+i-1)%L]-2*data[(i+1)%L]+data[(i+2)%L])/2/h3
] :
let(
first=(-5*data[0]+18*data[1]-24*data[2]+14*data[3]-3*data[4])/2,
second=(-3*data[0]+10*data[1]-12*data[2]+6*data[3]-data[4])/2,
last=(5*data[L-1]-18*data[L-2]+24*data[L-3]-14*data[L-4]+3*data[L-5])/2,
prelast=(3*data[L-1]-10*data[L-2]+12*data[L-3]-6*data[L-4]+data[L-5])/2
) [
first/h3,
second/h3,
for(i=[2:1:L-3]) (-data[i-2]+2*data[i-1]-2*data[i+1]+data[i+2])/2/h3,
prelast/h3,
last/h3
];
// Section: Complex Numbers
// Function: complex()
// Synopsis: Replaces scalars in a list or matrix with complex number 2-vectors.
// Topics: Math, Complex Numbers
// See Also: c_mul(), c_div(), c_conj(), c_real(), c_imag(), c_ident(), c_norm()
// Usage:
// z = complex(list)
// Description:
// Converts a real valued number, vector or matrix into its complex analog
// by replacing all entries with a 2-vector that has zero imaginary part.
function complex(list) =
is_num(list) ? [list,0] :
[for(entry=list) is_num(entry) ? [entry,0] : complex(entry)];
// Function: c_mul()
// Synopsis: Multiplies two complex numbers.
// Topics: Math, Complex Numbers
// See Also: complex(), c_mul(), c_div(), c_conj(), c_real(), c_imag(), c_ident(), c_norm()
// Usage:
// c = c_mul(z1,z2)
// Description:
// Multiplies two complex numbers, vectors or matrices, where complex numbers
// or entries are represented as vectors: [REAL, IMAGINARY]. Note that all
// entries in both arguments must be complex.
// Arguments:
// z1 = First complex number, vector or matrix
// z2 = Second complex number, vector or matrix
function c_mul(z1,z2) =
is_matrix([z1,z2],2,2) ? _c_mul(z1,z2) :
_combine_complex(_c_mul(_split_complex(z1), _split_complex(z2)));
function _split_complex(data) =
is_vector(data,2) ? data
: is_num(data[0][0]) ? [data*[1,0], data*[0,1]]
: [
[for(vec=data) vec * [1,0]],
[for(vec=data) vec * [0,1]]
];
function _combine_complex(data) =
is_vector(data,2) ? data
: is_num(data[0][0]) ? [for(i=[0:len(data[0])-1]) [data[0][i],data[1][i]]]
: [for(i=[0:1:len(data[0])-1])
[for(j=[0:1:len(data[0][0])-1])
[data[0][i][j], data[1][i][j]]]];
function _c_mul(z1,z2) =
[ z1.x*z2.x - z1.y*z2.y, z1.x*z2.y + z1.y*z2.x ];
// Function: c_div()
// Synopsis: Divides two complex numbers.
// Topics: Math, Complex Numbers
// See Also: complex(), c_mul(), c_div(), c_conj(), c_real(), c_imag(), c_ident(), c_norm()
// Usage:
// x = c_div(z1,z2)
// Description:
// Divides two complex numbers represented by 2D vectors.
// Returns a complex number as a 2D vector [REAL, IMAGINARY].
// Arguments:
// z1 = First complex number, given as a 2D vector [REAL, IMAGINARY]
// z2 = Second complex number, given as a 2D vector [REAL, IMAGINARY]
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];
// Function: c_conj()
// Synopsis: Returns the complex conjugate of the input.
// Topics: Math, Complex Numbers
// See Also: complex(), c_mul(), c_div(), c_conj(), c_real(), c_imag(), c_ident(), c_norm()
// Usage:
// w = c_conj(z)
// Description:
// Computes the complex conjugate of the input, which can be a complex number,
// complex vector or complex matrix.
function c_conj(z) =
is_vector(z,2) ? [z.x,-z.y] :
[for(entry=z) c_conj(entry)];
// Function: c_real()
// Synopsis: Returns the real part of a complex number, vector or matrix..
// Topics: Math, Complex Numbers
// See Also: complex(), c_mul(), c_div(), c_conj(), c_real(), c_imag(), c_ident(), c_norm()
// Usage:
// x = c_real(z)
// Description:
// Returns real part of a complex number, vector or matrix.
function c_real(z) =
is_vector(z,2) ? z.x
: is_num(z[0][0]) ? z*[1,0]
: [for(vec=z) vec * [1,0]];
// Function: c_imag()
// Synopsis: Returns the imaginary part of a complex number, vector or matrix..
// Topics: Math, Complex Numbers
// See Also: complex(), c_mul(), c_div(), c_conj(), c_real(), c_imag(), c_ident(), c_norm()
// Usage:
// x = c_imag(z)
// Description:
// Returns imaginary part of a complex number, vector or matrix.
function c_imag(z) =
is_vector(z,2) ? z.y
: is_num(z[0][0]) ? z*[0,1]
: [for(vec=z) vec * [0,1]];
// Function: c_ident()
// Synopsis: Returns an n by n complex identity matrix.
// Topics: Math, Complex Numbers
// See Also: complex(), c_mul(), c_div(), c_conj(), c_real(), c_imag(), c_ident(), c_norm()
// Usage:
// I = c_ident(n)
// Description:
// Produce an n by n complex identity matrix
function c_ident(n) = [for (i = [0:1:n-1]) [for (j = [0:1:n-1]) (i==j)?[1,0]:[0,0]]];
// Function: c_norm()
// Synopsis: Returns the norm of a complex number or vector.
// Topics: Math, Complex Numbers
// See Also: complex(), c_mul(), c_div(), c_conj(), c_real(), c_imag(), c_ident(), c_norm()
// Usage:
// n = c_norm(z)
// Description:
// Compute the norm of a complex number or vector.
function c_norm(z) = norm_fro(z);
// Section: Polynomials
// Function: quadratic_roots()
// Synopsis: Computes roots for the quadratic equation.
// Topics: Math, Geometry, Complex Numbers
// See Also: quadratic_roots(), polynomial(), poly_mult(), poly_div(), poly_add()
// Usage:
// roots = quadratic_roots(a, b, c, [real])
// Description:
// Computes roots of the quadratic equation a*x^2+b*x+c==0, where the
// coefficients are real numbers. If real is true then returns only the
// real roots. Otherwise returns a pair of complex values. This method
// may be more reliable than the general root finder at distinguishing
// real roots from complex roots.
// Algorithm from: https://people.csail.mit.edu/bkph/articles/Quadratics.pdf
function quadratic_roots(a,b,c,real=false) =
real ? [for(root = quadratic_roots(a,b,c,real=false)) if (root.y==0) root.x]
:
is_undef(b) && is_undef(c) && is_vector(a,3) ? quadratic_roots(a[0],a[1],a[2]) :
assert(is_num(a) && is_num(b) && is_num(c))
assert(a!=0 || b!=0 || c!=0, "Quadratic must have a nonzero coefficient")
a==0 && b==0 ? [] : // No solutions
a==0 ? [[-c/b,0]] :
let(
descrim = b*b-4*a*c,
sqrt_des = sqrt(abs(descrim))
)
descrim < 0 ? // Complex case
[[-b, sqrt_des],
[-b, -sqrt_des]]/2/a :
b<0 ? // b positive
[[2*c/(-b+sqrt_des),0],
[(-b+sqrt_des)/a/2,0]]
: // b negative
[[(-b-sqrt_des)/2/a, 0],
[2*c/(-b-sqrt_des),0]];
// Function: polynomial()
// Synopsis: Calculates a polynomial equation at a given value.
// Topics: Math, Complex Numbers
// See Also: quadratic_roots(), polynomial(), poly_mult(), poly_div(), poly_add(), poly_roots()
// Usage:
// x = polynomial(p, z)
// Description:
// 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.
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_mul(total,z)+[p[k],0]);
// Function: poly_mult()
// Synopsis: Returns the polynomial result of multiplying two polynomial equations.
// Topics: Math
// See Also: quadratic_roots(), polynomial(), poly_mult(), poly_div(), poly_add(), poly_roots()
// Usage:
// x = polymult(p,q)
// x = polymult([p1,p2,p3,...])
// Description:
// Given a list of polynomials represented as real algebraic coefficient lists, with the highest degree coefficient first,
// computes the coefficient list of the product polynomial.
function poly_mult(p,q) =
is_undef(q) ?
len(p)==2
? poly_mult(p[0],p[1])
: poly_mult(p[0], poly_mult(list_tail(p)))
:
assert( is_vector(p) && is_vector(q),"Invalid arguments to poly_mult")
p*p==0 || q*q==0
? [0]
: _poly_trim(convolve(p,q));
// Function: poly_div()
// Synopsis: Returns the polynomial quotient and remainder results of dividing two polynomial equations.
// Topics: Math
// See Also: quadratic_roots(), polynomial(), poly_mult(), poly_div(), poly_add(), poly_roots()
// Usage:
// [quotient,remainder] = poly_div(n,d)
// Description:
// Computes division of the numerator polynomial by the denominator polynomial and returns
// a list of two polynomials, [quotient, remainder]. If the division has no remainder then
// the zero polynomial [0] is returned for the remainder. Similarly if the quotient is zero
// the returned quotient will be [0].
function poly_div(n,d) =
assert( is_vector(n) && is_vector(d) , "Invalid polynomials." )
let( d = _poly_trim(d),
n = _poly_trim(n) )
assert( d!=[0] , "Denominator cannot be a zero polynomial." )
n==[0]
? [[0],[0]]
: _poly_div(n,d,q=[]);
function _poly_div(n,d,q) =
len(n)<len(d) ? [q,_poly_trim(n)] :
let(
t = n[0] / d[0],
newq = concat(q,[t]),
newn = [for(i=[1:1:len(n)-1]) i<len(d) ? n[i] - t*d[i] : n[i]]
)
_poly_div(newn,d,newq);
/// Internal Function: _poly_trim()
/// Usage:
/// _poly_trim(p, [eps])
/// Description:
/// Removes leading zero terms of a polynomial. By default zeros must be exact,
/// or give epsilon for approximate zeros. Returns [0] for a zero polynomial.
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 ? [0] : list_tail(p,nz[0]);
// Function: poly_add()
// Synopsis: Returns the polynomial sum of adding two polynomial equations.
// Topics: Math
// See Also: quadratic_roots(), polynomial(), poly_mult(), poly_div(), poly_add(), poly_roots()
// Usage:
// sum = poly_add(p,q)
// 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,
short = plen>qlen ? q : p
)
_poly_trim(long + concat(repeat(0,len(long)-len(short)),short));
// Function: poly_roots()
// Synopsis: Returns all complex number roots of the given real polynomial.
// Topics: Math, Complex Numbers
// See Also: quadratic_roots(), polynomial(), poly_mult(), poly_div(), poly_add(), poly_roots()
// Usage:
// roots = poly_roots(p, [tol]);
// Description:
// Returns all complex roots of the specified real polynomial p.
// The polynomial is specified as p=[a_n, a_{n-1},...,a_1,a_0]
// where a_n is the z^n coefficient. The tol parameter gives
// the stopping tolerance for the iteration. The polynomial
// must have at least one non-zero coefficient. Convergence is poor
// if the polynomial has any repeated roots other than zero.
// Arguments:
// p = polynomial coefficients with higest power coefficient first
// tol = tolerance for iteration. Default: 1e-14
// Uses the Aberth method https://en.wikipedia.org/wiki/Aberth_method
//
// 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( 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(list_head(p),tol=tol, error_bound=error_bound))
(error_bound ? [ [[0,0], each solutions[0]], [0, each solutions[1]]]
: [[0,0], each solutions]) :
len(p)==1 ? (error_bound ? [[],[]] : []) : // Nonzero constant case has no solutions
len(p)==2 ? let( solution = [[-p[1]/p[0],0]]) // Linear case needs special handling
(error_bound ? [solution,[0]] : solution)
:
let(
n = len(p)-1, // polynomial degree
pderiv = [for(i=[0:n-1]) p[i]*(n-i)],
s = [for(i=[0:1:n]) abs(p[i])*(4*(n-i)+1)], // Error bound polynomial from Bini
// Using method from: http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/0915-24.pdf
beta = -p[1]/p[0]/n,
r = 1+pow(abs(polynomial(p,beta)/p[0]),1/n),
init = [for(i=[0:1:n-1]) // Initial guess for roots
let(angle = 360*i/n+270/n/PI)
[beta,0]+r*[cos(angle),sin(angle)]
],
roots = _poly_roots(p,pderiv,s,init,tol=tol),
error = error_bound ? [for(xi=roots) n * (norm(polynomial(p,xi))+tol*polynomial(s,norm(xi))) /
abs(norm(polynomial(pderiv,xi))-tol*polynomial(s,norm(xi)))] : 0
)
error_bound ? [roots, error] : roots;
// Internal function
// p = polynomial
// pderiv = derivative polynomial of p
// z = current guess for the roots
// tol = root tolerance
// i=iteration counter
function _poly_roots(p, pderiv, s, z, tol, i=0) =
assert(i<45, str("Polyroot exceeded iteration limit. Current solution:", z))
let(
n = len(z),
svals = [for(zk=z) tol*polynomial(s,norm(zk))],
p_of_z = [for(zk=z) polynomial(p,zk)],
done = [for(k=[0:n-1]) norm(p_of_z[k])<=svals[k]],
newton = [for(k=[0:n-1]) c_div(p_of_z[k], polynomial(pderiv,z[k]))],
zdiff = [for(k=[0:n-1]) sum([for(j=[0:n-1]) if (j!=k) c_div([1,0], z[k]-z[j])])],
w = [for(k=[0:n-1]) done[k] ? [0,0] : c_div( newton[k],
[1,0] - c_mul(newton[k], zdiff[k]))]
)
all(done) ? z : _poly_roots(p,pderiv,s,z-w,tol,i+1);
// Function: real_roots()
// Synopsis: Returns all real roots of the given real polynomial.
// Topics: Math, Complex Numbers
// See Also: quadratic_roots(), polynomial(), poly_mult(), poly_div(), poly_add(), poly_roots()
// Usage:
// roots = real_roots(p, [eps], [tol])
// Description:
// Returns the real roots of the specified real polynomial p.
// The polynomial is specified as p=[a_n, a_{n-1},...,a_1,a_0]
// where a_n is the x^n coefficient. This function works by
// computing the complex roots and returning those roots where
// the imaginary part is closed to zero. By default it uses a computed
// error bound from the polynomial solver to decide whether imaginary
// parts are zero. You can specify eps, in which case the test is
// z.y/(1+norm(z)) < eps. Because
// of poor convergence and higher error for repeated roots, such roots may
// be missed by the algorithm because their imaginary part is large.
// Arguments:
// p = polynomial to solve as coefficient list, highest power term first
// eps = used to determine whether imaginary parts of roots are zero
// 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]
: [for(i=idx(roots)) if (abs(roots[i].y)<=err[i]) roots[i].x];
// Section: Operations on Functions
// Function: root_find()
// Synopsis: Finds a root of the given continuous function.
// Topics: Math
// See Also: quadratic_roots(), polynomial(), poly_mult(), poly_div(), poly_add(), poly_roots()
// Usage:
// x = root_find(f, x0, x1, [tol])
// Description:
// Find a root of the continuous function f where the sign of f(x0) is different
// from the sign of f(x1). The function f is a function literal accepting one
// argument. You must have a version of OpenSCAD that supports function literals
// (2021.01 or newer). The tolerance (tol) specifies the accuracy of the solution:
// abs(f(x)) < tol * yrange, where yrange is the range of observed function values.
// This function can only find roots that cross the x axis: it cannot find the
// the root of x^2.
// Arguments:
// f = function literal for a scalar-valued single variable function
// x0 = endpoint of interval to search for root
// x1 = second endpoint of interval to search for root
// tol = tolerance for solution. Default: 1e-15
// The algorithm is based on Brent's method and is a combination of
// bisection and inverse quadratic approximation, where bisection occurs
// at every step, with refinement using inverse quadratic approximation
// only when that approximation gives a good result. The detail
// of how to decide when to use the quadratic came from an article
// by Crenshaw on "The World's Best Root Finder".
// https://www.embedded.com/worlds-best-root-finder/
function root_find(f,x0,x1,tol=1e-15) =
let(
y0 = f(x0),
y1 = f(x1),
yrange = y0<y1 ? [y0,y1] : [y1,y0]
)
// Check endpoints
y0==0 || _rfcheck(x0, y0,yrange,tol) ? x0 :
y1==0 || _rfcheck(x1, y1,yrange,tol) ? x1 :
assert(y0*y1<0, "Sign of function must be different at the interval endpoints")
_rootfind(f,[x0,x1],[y0,y1],yrange,tol);
function _rfcheck(x,y,range,tol) =
assert(is_finite(y), str("Function not finite at ",x))
abs(y) < tol*(range[1]-range[0]);
// xpts and ypts are arrays whose first two entries contain the
// interval bracketing the root. Extra entries are ignored.
// yrange is the total observed range of y values (used for the
// tolerance test).
function _rootfind(f, xpts, ypts, yrange, tol, i=0) =
assert(i<100, "root_find did not converge to a solution")
let(
xmid = (xpts[0]+xpts[1])/2,
ymid = f(xmid),
yrange = [min(ymid, yrange[0]), max(ymid, yrange[1])]
)
_rfcheck(xmid, ymid, yrange, tol) ? xmid :
let(
// Force root to be between x0 and midpoint
y = ymid * ypts[0] < 0 ? [ypts[0], ymid, ypts[1]]
: [ypts[1], ymid, ypts[0]],
x = ymid * ypts[0] < 0 ? [xpts[0], xmid, xpts[1]]
: [xpts[1], xmid, xpts[0]],
v = y[2]*(y[2]-y[0]) - 2*y[1]*(y[1]-y[0])
)
v <= 0 ? _rootfind(f,x,y,yrange,tol,i+1) // Root is between first two points, extra 3rd point doesn't hurt
:
let( // Do quadratic approximation
B = (x[1]-x[0]) / (y[1]-y[0]),
C = y*[-1,2,-1] / (y[2]-y[1]) / (y[2]-y[0]),
newx = x[0] - B * y[0] *(1-C*y[1]),
newy = f(newx),
new_yrange = [min(yrange[0],newy), max(yrange[1], newy)],
// select interval that contains the root by checking sign
yinterval = newy*y[0] < 0 ? [y[0],newy] : [newy,y[1]],
xinterval = newy*y[0] < 0 ? [x[0],newx] : [newx,x[1]]
)
_rfcheck(newx, newy, new_yrange, tol)
? newx
: _rootfind(f, xinterval, yinterval, new_yrange, tol, i+1);
// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap
Reference in a new issue View git blame Copy permalink