1//////////////////////////////////////////////////////////////////////
2// LibFile: math.scad
3// Assorted math functions, including linear interpolation, list operations (sums, mean, products),
4// convolution, quantization, log2, hyperbolic trig functions, random numbers, derivatives,
5// polynomials, and root finding.
6// Includes:
7// include <BOSL2/std.scad>
8// FileGroup: Math
9// FileSummary: Math on lists, special functions, quantization, random numbers, calculus, root finding
10//
11// FileFootnotes: STD=Included in std.scad
12//////////////////////////////////////////////////////////////////////
13
14// Section: Math Constants
15
16// Constant: PHI
17// Description: The golden ratio phi.
18PHI = (1+sqrt(5))/2;
19
20// Constant: EPSILON
21// Description: A really small value useful in comparing floating point numbers. ie: abs(a-b)<EPSILON
22EPSILON = 1e-9;
23
24// Constant: INF
25// Description: The value `inf`, useful for comparisons.
26INF = 1/0;
27
28// Constant: NAN
29// Description: The value `nan`, useful for comparisons.
30NAN = acos(2);
31
32
33
34// Section: Interpolation and Counting
35
36
37// Function: count()
38// Usage:
39// list = count(n, [s], [step], [reverse]);
40// Description:
41// Creates a list of `n` numbers, starting at `s`, incrementing by `step` each time.
42// You can also pass a list for n and then the length of the input list is used.
43// Arguments:
44// n = The length of the list of numbers to create, or a list to match the length of
45// s = The starting value of the list of numbers.
46// step = The amount to increment successive numbers in the list.
47// reverse = Reverse the list. Default: false.
48// See Also: idx()
49// Example:
50// nl1 = count(5); // Returns: [0,1,2,3,4]
51// nl2 = count(5,3); // Returns: [3,4,5,6,7]
52// nl3 = count(4,3,2); // Returns: [3,5,7,9]
53// nl4 = count(5,reverse=true); // Returns: [4,3,2,1,0]
54// nl5 = count(5,3,reverse=true); // Returns: [7,6,5,4,3]
55function count(n,s=0,step=1,reverse=false) = let(n=is_list(n) ? len(n) : n)
56 reverse? [for (i=[n-1:-1:0]) s+i*step]
57 : [for (i=[0:1:n-1]) s+i*step];
58
59
60// Function: lerp()
61// Usage:
62// x = lerp(a, b, u);
63// l = lerp(a, b, LIST);
64// Description:
65// Interpolate between two values or vectors.
66// If `u` is given as a number, returns the single interpolated value.
67// If `u` is 0.0, then the value of `a` is returned.
68// If `u` is 1.0, then the value of `b` is returned.
69// If `u` is a range, or list of numbers, returns a list of interpolated values.
70// It is valid to use a `u` value outside the range 0 to 1. The result will be an extrapolation
71// along the slope formed by `a` and `b`.
72// Arguments:
73// a = First value or vector.
74// b = Second value or vector.
75// 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.
76// Example:
77// x = lerp(0,20,0.3); // Returns: 6
78// x = lerp(0,20,0.8); // Returns: 16
79// x = lerp(0,20,-0.1); // Returns: -2
80// x = lerp(0,20,1.1); // Returns: 22
81// p = lerp([0,0],[20,10],0.25); // Returns [5,2.5]
82// l = lerp(0,20,[0.4,0.6]); // Returns: [8,12]
83// l = lerp(0,20,[0.25:0.25:0.75]); // Returns: [5,10,15]
84// Example(2D):
85// p1 = [-50,-20]; p2 = [50,30];
86// stroke([p1,p2]);
87// pts = lerp(p1, p2, [0:1/8:1]);
88// // Points colored in ROYGBIV order.
89// rainbow(pts) translate($item) circle(d=3,$fn=8);
90function lerp(a,b,u) =
91 assert(same_shape(a,b), "Bad or inconsistent inputs to lerp")
92 is_finite(u)? (1-u)*a + u*b :
93 assert(is_finite(u) || is_vector(u) || valid_range(u), "Input u to lerp must be a number, vector, or valid range.")
94 [for (v = u) (1-v)*a + v*b ];
95
96
97// Function: lerpn()
98// Usage:
99// x = lerpn(a, b, n);
100// x = lerpn(a, b, n, [endpoint]);
101// Description:
102// Returns exactly `n` values, linearly interpolated between `a` and `b`.
103// If `endpoint` is true, then the last value will exactly equal `b`.
104// If `endpoint` is false, then the last value will be `a+(b-a)*(1-1/n)`.
105// Arguments:
106// a = First value or vector.
107// b = Second value or vector.
108// n = The number of values to return.
109// endpoint = If true, the last value will be exactly `b`. If false, the last value will be one step less.
110// Example:
111// l = lerpn(-4,4,9); // Returns: [-4,-3,-2,-1,0,1,2,3,4]
112// l = lerpn(-4,4,8,false); // Returns: [-4,-3,-2,-1,0,1,2,3]
113// l = lerpn(0,1,6); // Returns: [0, 0.2, 0.4, 0.6, 0.8, 1]
114// l = lerpn(0,1,5,false); // Returns: [0, 0.2, 0.4, 0.6, 0.8]
115function lerpn(a,b,n,endpoint=true) =
116 assert(same_shape(a,b), "Bad or inconsistent inputs to lerpn")
117 assert(is_int(n))
118 assert(is_bool(endpoint))
119 let( d = n - (endpoint? 1 : 0) )
120 [for (i=[0:1:n-1]) let(u=i/d) (1-u)*a + u*b];
121
122
123
124// Section: Miscellaneous Functions
125
126// Function: sqr()
127// Usage:
128// x2 = sqr(x);
129// Description:
130// If given a number, returns the square of that number,
131// If given a vector, returns the sum-of-squares/dot product of the vector elements.
132// If given a matrix, returns the matrix multiplication of the matrix with itself.
133// Example:
134// sqr(3); // Returns: 9
135// sqr(-4); // Returns: 16
136// sqr([2,3,4]); // Returns: 29
137// sqr([[1,2],[3,4]]); // Returns [[7,10],[15,22]]
138function sqr(x) =
139 assert(is_finite(x) || is_vector(x) || is_matrix(x), "Input is not a number nor a list of numbers.")
140 x*x;
141
142
143// Function: log2()
144// Usage:
145// val = log2(x);
146// Description:
147// Returns the logarithm base 2 of the value given.
148// Example:
149// log2(0.125); // Returns: -3
150// log2(16); // Returns: 4
151// log2(256); // Returns: 8
152function log2(x) =
153 assert( is_finite(x), "Input is not a number.")
154 ln(x)/ln(2);
155
156// this may return NAN or INF; should it check x>0 ?
157
158// Function: hypot()
159// Usage:
160// l = hypot(x, y, [z]);
161// Description:
162// Calculate hypotenuse length of a 2D or 3D triangle.
163// Arguments:
164// x = Length on the X axis.
165// y = Length on the Y axis.
166// z = Length on the Z axis. Optional.
167// Example:
168// l = hypot(3,4); // Returns: 5
169// l = hypot(3,4,5); // Returns: ~7.0710678119
170function hypot(x,y,z=0) =
171 assert( is_vector([x,y,z]), "Improper number(s).")
172 norm([x,y,z]);
173
174
175// Function: factorial()
176// Usage:
177// x = factorial(n, [d]);
178// Description:
179// Returns the factorial of the given integer value, or n!/d! if d is given.
180// Arguments:
181// n = The integer number to get the factorial of. (n!)
182// d = If given, the returned value will be (n! / d!)
183// Example:
184// x = factorial(4); // Returns: 24
185// y = factorial(6); // Returns: 720
186// z = factorial(9); // Returns: 362880
187function factorial(n,d=0) =
188 assert(is_int(n) && is_int(d) && n>=0 && d>=0, "Factorial is defined only for non negative integers")
189 assert(d<=n, "d cannot be larger than n")
190 product([1,for (i=[n:-1:d+1]) i]);
191
192
193// Function: binomial()
194// Usage:
195// x = binomial(n);
196// Description:
197// Returns the binomial coefficients of the integer `n`.
198// Arguments:
199// n = The integer to get the binomial coefficients of
200// Example:
201// x = binomial(3); // Returns: [1,3,3,1]
202// y = binomial(4); // Returns: [1,4,6,4,1]
203// z = binomial(6); // Returns: [1,6,15,20,15,6,1]
204function binomial(n) =
205 assert( is_int(n) && n>0, "Input is not an integer greater than 0.")
206 [for( c = 1, i = 0;
207 i<=n;
208 c = c*(n-i)/(i+1), i = i+1
209 ) c ] ;
210
211
212// Function: binomial_coefficient()
213// Usage:
214// x = binomial_coefficient(n, k);
215// Description:
216// Returns the k-th binomial coefficient of the integer `n`.
217// Arguments:
218// n = The integer to get the binomial coefficient of
219// k = The binomial coefficient index
220// Example:
221// x = binomial_coefficient(3,2); // Returns: 3
222// y = binomial_coefficient(10,6); // Returns: 210
223function binomial_coefficient(n,k) =
224 assert( is_int(n) && is_int(k), "Some input is not a number.")
225 k < 0 || k > n ? 0 :
226 k ==0 || k ==n ? 1 :
227 let( k = min(k, n-k),
228 b = [for( c = 1, i = 0;
229 i<=k;
230 c = c*(n-i)/(i+1), i = i+1
231 ) c] )
232 b[len(b)-1];
233
234
235// Function: gcd()
236// Usage:
237// x = gcd(a,b)
238// Description:
239// Computes the Greatest Common Divisor/Factor of `a` and `b`.
240function gcd(a,b) =
241 assert(is_int(a) && is_int(b),"Arguments to gcd must be integers")
242 b==0 ? abs(a) : gcd(b,a % b);
243
244
245// Computes lcm for two integers
246function _lcm(a,b) =
247 assert(is_int(a) && is_int(b), "Invalid non-integer parameters to lcm")
248 assert(a!=0 && b!=0, "Arguments to lcm should not be zero")
249 abs(a*b) / gcd(a,b);
250
251
252// Computes lcm for a list of values
253function _lcmlist(a) =
254 len(a)==1 ? a[0] :
255 _lcmlist(concat(lcm(a[0],a[1]),list_tail(a,2)));
256
257
258// Function: lcm()
259// Usage:
260// div = lcm(a, b);
261// divs = lcm(list);
262// Description:
263// Computes the Least Common Multiple of the two arguments or a list of arguments. Inputs should
264// be non-zero integers. The output is always a positive integer. It is an error to pass zero
265// as an argument.
266function lcm(a,b=[]) =
267 !is_list(a) && !is_list(b)
268 ? _lcm(a,b)
269 : let( arglist = concat(force_list(a),force_list(b)) )
270 assert(len(arglist)>0, "Invalid call to lcm with empty list(s)")
271 _lcmlist(arglist);
272
273
274
275
276// Section: Hyperbolic Trigonometry
277
278// Function: sinh()
279// Usage:
280// a = sinh(x);
281// Description: Takes a value `x`, and returns the hyperbolic sine of it.
282function sinh(x) =
283 assert(is_finite(x), "The input must be a finite number.")
284 (exp(x)-exp(-x))/2;
285
286
287// Function: cosh()
288// Usage:
289// a = cosh(x);
290// Description: Takes a value `x`, and returns the hyperbolic cosine of it.
291function cosh(x) =
292 assert(is_finite(x), "The input must be a finite number.")
293 (exp(x)+exp(-x))/2;
294
295
296// Function: tanh()
297// Usage:
298// a = tanh(x);
299// Description: Takes a value `x`, and returns the hyperbolic tangent of it.
300function tanh(x) =
301 assert(is_finite(x), "The input must be a finite number.")
302 sinh(x)/cosh(x);
303
304
305// Function: asinh()
306// Usage:
307// a = asinh(x);
308// Description: Takes a value `x`, and returns the inverse hyperbolic sine of it.
309function asinh(x) =
310 assert(is_finite(x), "The input must be a finite number.")
311 ln(x+sqrt(x*x+1));
312
313
314// Function: acosh()
315// Usage:
316// a = acosh(x);
317// Description: Takes a value `x`, and returns the inverse hyperbolic cosine of it.
318function acosh(x) =
319 assert(is_finite(x), "The input must be a finite number.")
320 ln(x+sqrt(x*x-1));
321
322
323// Function: atanh()
324// Usage:
325// a = atanh(x);
326// Description: Takes a value `x`, and returns the inverse hyperbolic tangent of it.
327function atanh(x) =
328 assert(is_finite(x), "The input must be a finite number.")
329 ln((1+x)/(1-x))/2;
330
331
332// Section: Quantization
333
334// Function: quant()
335// Usage:
336// num = quant(x, y);
337// Description:
338// Quantize a value `x` to an integer multiple of `y`, rounding to the nearest multiple.
339// The value of `y` does NOT have to be an integer. If `x` is a list, then every item
340// in that list will be recursively quantized.
341// Arguments:
342// x = The value to quantize.
343// y = The non-zero integer quantum of the quantization.
344// Example:
345// a = quant(12,4); // Returns: 12
346// b = quant(13,4); // Returns: 12
347// c = quant(13.1,4); // Returns: 12
348// d = quant(14,4); // Returns: 16
349// e = quant(14.1,4); // Returns: 16
350// f = quant(15,4); // Returns: 16
351// g = quant(16,4); // Returns: 16
352// h = quant(9,3); // Returns: 9
353// i = quant(10,3); // Returns: 9
354// j = quant(10.4,3); // Returns: 9
355// k = quant(10.5,3); // Returns: 12
356// l = quant(11,3); // Returns: 12
357// m = quant(12,3); // Returns: 12
358// n = quant(11,2.5); // Returns: 10
359// o = quant(12,2.5); // Returns: 12.5
360// p = quant([12,13,13.1,14,14.1,15,16],4); // Returns: [12,12,12,16,16,16,16]
361// q = quant([9,10,10.4,10.5,11,12],3); // Returns: [9,9,9,12,12,12]
362// r = quant([[9,10,10.4],[10.5,11,12]],3); // Returns: [[9,9,9],[12,12,12]]
363function quant(x,y) =
364 assert( is_finite(y) && y>0, "The quantum `y` must be a non zero integer.")
365 is_list(x)
366 ? [for (v=x) quant(v,y)]
367 : assert( is_finite(x), "The input to quantize is not a number nor a list of numbers.")
368 floor(x/y+0.5)*y;
369
370
371// Function: quantdn()
372// Usage:
373// num = quantdn(x, y);
374// Description:
375// Quantize a value `x` to an integer multiple of `y`, rounding down to the previous multiple.
376// The value of `y` does NOT have to be an integer. If `x` is a list, then every item in that
377// list will be recursively quantized down.
378// Arguments:
379// x = The value to quantize.
380// y = The non-zero integer quantum of the quantization.
381// Example:
382// a = quantdn(12,4); // Returns: 12
383// b = quantdn(13,4); // Returns: 12
384// c = quantdn(13.1,4); // Returns: 12
385// d = quantdn(14,4); // Returns: 12
386// e = quantdn(14.1,4); // Returns: 12
387// f = quantdn(15,4); // Returns: 12
388// g = quantdn(16,4); // Returns: 16
389// h = quantdn(9,3); // Returns: 9
390// i = quantdn(10,3); // Returns: 9
391// j = quantdn(10.4,3); // Returns: 9
392// k = quantdn(10.5,3); // Returns: 9
393// l = quantdn(11,3); // Returns: 9
394// m = quantdn(12,3); // Returns: 12
395// n = quantdn(11,2.5); // Returns: 10
396// o = quantdn(12,2.5); // Returns: 10
397// p = quantdn([12,13,13.1,14,14.1,15,16],4); // Returns: [12,12,12,12,12,12,16]
398// q = quantdn([9,10,10.4,10.5,11,12],3); // Returns: [9,9,9,9,9,12]
399// r = quantdn([[9,10,10.4],[10.5,11,12]],3); // Returns: [[9,9,9],[9,9,12]]
400function quantdn(x,y) =
401 assert( is_finite(y) && y>0, "The quantum `y` must be a non zero integer.")
402 is_list(x)
403 ? [for (v=x) quantdn(v,y)]
404 : assert( is_finite(x), "The input to quantize must be a number or a list of numbers.")
405 floor(x/y)*y;
406
407
408// Function: quantup()
409// Usage:
410// num = quantup(x, y);
411// Description:
412// Quantize a value `x` to an integer multiple of `y`, rounding up to the next multiple.
413// The value of `y` does NOT have to be an integer. If `x` is a list, then every item in
414// that list will be recursively quantized up.
415// Arguments:
416// x = The value to quantize.
417// y = The non-zero integer quantum of the quantization.
418// Example:
419// a = quantup(12,4); // Returns: 12
420// b = quantup(13,4); // Returns: 16
421// c = quantup(13.1,4); // Returns: 16
422// d = quantup(14,4); // Returns: 16
423// e = quantup(14.1,4); // Returns: 16
424// f = quantup(15,4); // Returns: 16
425// g = quantup(16,4); // Returns: 16
426// h = quantup(9,3); // Returns: 9
427// i = quantup(10,3); // Returns: 12
428// j = quantup(10.4,3); // Returns: 12
429// k = quantup(10.5,3); // Returns: 12
430// l = quantup(11,3); // Returns: 12
431// m = quantup(12,3); // Returns: 12
432// n = quantdn(11,2.5); // Returns: 12.5
433// o = quantdn(12,2.5); // Returns: 12.5
434// p = quantup([12,13,13.1,14,14.1,15,16],4); // Returns: [12,16,16,16,16,16,16]
435// q = quantup([9,10,10.4,10.5,11,12],3); // Returns: [9,12,12,12,12,12]
436// r = quantup([[9,10,10.4],[10.5,11,12]],3); // Returns: [[9,12,12],[12,12,12]]
437function quantup(x,y) =
438 assert( is_finite(y) && y>0, "The quantum `y` must be a non zero integer.")
439 is_list(x)
440 ? [for (v=x) quantup(v,y)]
441 : assert( is_finite(x), "The input to quantize must be a number or a list of numbers.")
442 ceil(x/y)*y;
443
444
445// Section: Constraints and Modulos
446
447// Function: constrain()
448// Usage:
449// val = constrain(v, minval, maxval);
450// Description:
451// Constrains value to a range of values between minval and maxval, inclusive.
452// Arguments:
453// v = value to constrain.
454// minval = minimum value to return, if out of range.
455// maxval = maximum value to return, if out of range.
456// Example:
457// a = constrain(-5, -1, 1); // Returns: -1
458// b = constrain(5, -1, 1); // Returns: 1
459// c = constrain(0.3, -1, 1); // Returns: 0.3
460// d = constrain(9.1, 0, 9); // Returns: 9
461// e = constrain(-0.1, 0, 9); // Returns: 0
462function constrain(v, minval, maxval) =
463 assert( is_finite(v+minval+maxval), "Input must be finite number(s).")
464 min(maxval, max(minval, v));
465
466
467// Function: posmod()
468// Usage:
469// mod = posmod(x, m)
470// Description:
471// Returns the positive modulo `m` of `x`. Value returned will be in the range 0 ... `m`-1.
472// Arguments:
473// x = The value to constrain.
474// m = Modulo value.
475// Example:
476// a = posmod(-700,360); // Returns: 340
477// b = posmod(-270,360); // Returns: 90
478// c = posmod(-120,360); // Returns: 240
479// d = posmod(120,360); // Returns: 120
480// e = posmod(270,360); // Returns: 270
481// f = posmod(700,360); // Returns: 340
482// g = posmod(3,2.5); // Returns: 0.5
483function posmod(x,m) =
484 assert( is_finite(x) && is_finite(m) && !approx(m,0) , "Input must be finite numbers. The divisor cannot be zero.")
485 (x%m+m)%m;
486
487
488// Function: modang()
489// Usage:
490// ang = modang(x);
491// Description:
492// Takes an angle in degrees and normalizes it to an equivalent angle value between -180 and 180.
493// Example:
494// a1 = modang(-700,360); // Returns: 20
495// a2 = modang(-270,360); // Returns: 90
496// a3 = modang(-120,360); // Returns: -120
497// a4 = modang(120,360); // Returns: 120
498// a5 = modang(270,360); // Returns: -90
499// a6 = modang(700,360); // Returns: -20
500function modang(x) =
501 assert( is_finite(x), "Input must be a finite number.")
502 let(xx = posmod(x,360)) xx<180? xx : xx-360;
503
504
505
506// Section: Operations on Lists (Sums, Mean, Products)
507
508// Function: sum()
509// Usage:
510// x = sum(v, [dflt]);
511// Description:
512// Returns the sum of all entries in the given consistent list.
513// If passed an array of vectors, returns the sum the vectors.
514// If passed an array of matrices, returns the sum of the matrices.
515// If passed an empty list, the value of `dflt` will be returned.
516// Arguments:
517// v = The list to get the sum of.
518// dflt = The default value to return if `v` is an empty list. Default: 0
519// Example:
520// sum([1,2,3]); // returns 6.
521// sum([[1,2,3], [3,4,5], [5,6,7]]); // returns [9, 12, 15]
522function sum(v, dflt=0) =
523 v==[]? dflt :
524 assert(is_consistent(v), "Input to sum is non-numeric or inconsistent")
525 is_finite(v[0]) || is_vector(v[0]) ? [for(i=v) 1]*v :
526 _sum(v,v[0]*0);
527
528function _sum(v,_total,_i=0) = _i>=len(v) ? _total : _sum(v,_total+v[_i], _i+1);
529
530
531
532
533// Function: mean()
534// Usage:
535// x = mean(v);
536// Description:
537// Returns the arithmetic mean/average of all entries in the given array.
538// If passed a list of vectors, returns a vector of the mean of each part.
539// Arguments:
540// v = The list of values to get the mean of.
541// Example:
542// mean([2,3,4]); // returns 3.
543// mean([[1,2,3], [3,4,5], [5,6,7]]); // returns [3, 4, 5]
544function mean(v) =
545 assert(is_list(v) && len(v)>0, "Invalid list.")
546 sum(v)/len(v);
547
548
549
550// Function: median()
551// Usage:
552// middle = median(v)
553// Description:
554// Returns the median of the given vector.
555function median(v) =
556 assert(is_vector(v), "Input to median must be a vector")
557 len(v)%2 ? max( list_smallest(v, ceil(len(v)/2)) ) :
558 let( lowest = list_smallest(v, len(v)/2 + 1),
559 max = max(lowest),
560 imax = search(max,lowest,1),
561 max2 = max([for(i=idx(lowest)) if(i!=imax[0]) lowest[i] ])
562 )
563 (max+max2)/2;
564
565
566// Function: deltas()
567// Usage:
568// delts = deltas(v,[wrap]);
569// Description:
570// Returns a list with the deltas of adjacent entries in the given list, optionally wrapping back to the front.
571// The list should be a consistent list of numeric components (numbers, vectors, matrix, etc).
572// Given [a,b,c,d], returns [b-a,c-b,d-c].
573// Arguments:
574// v = The list to get the deltas of.
575// 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
576// Example:
577// deltas([2,5,9,17]); // returns [3,4,8].
578// deltas([[1,2,3], [3,6,8], [4,8,11]]); // returns [[2,4,5], [1,2,3]]
579function deltas(v, wrap=false) =
580 assert( is_consistent(v) && len(v)>1 , "Inconsistent list or with length<=1.")
581 [for (p=pair(v,wrap)) p[1]-p[0]] ;
582
583
584// Function: cumsum()
585// Usage:
586// sums = cumsum(v);
587// Description:
588// Returns a list where each item is the cumulative sum of all items up to and including the corresponding entry in the input list.
589// If passed an array of vectors, returns a list of cumulative vectors sums.
590// Arguments:
591// v = The list to get the sum of.
592// Example:
593// cumsum([1,1,1]); // returns [1,2,3]
594// cumsum([2,2,2]); // returns [2,4,6]
595// cumsum([1,2,3]); // returns [1,3,6]
596// cumsum([[1,2,3], [3,4,5], [5,6,7]]); // returns [[1,2,3], [4,6,8], [9,12,15]]
597function cumsum(v) =
598 assert(is_consistent(v), "The input is not consistent." )
599 len(v)<=1 ? v :
600 _cumsum(v,_i=1,_acc=[v[0]]);
601
602function _cumsum(v,_i=0,_acc=[]) =
603 _i>=len(v) ? _acc :
604 _cumsum( v, _i+1, [ each _acc, _acc[len(_acc)-1] + v[_i] ] );
605
606
607
608// Function: product()
609// Usage:
610// x = product(v);
611// Description:
612// Returns the product of all entries in the given list.
613// If passed a list of vectors of same dimension, returns a vector of products of each part.
614// If passed a list of square matrices, returns the resulting product matrix.
615// Arguments:
616// v = The list to get the product of.
617// Example:
618// product([2,3,4]); // returns 24.
619// product([[1,2,3], [3,4,5], [5,6,7]]); // returns [15, 48, 105]
620function product(v) =
621 assert( is_vector(v) || is_matrix(v) || ( is_matrix(v[0],square=true) && is_consistent(v)),
622 "Invalid input.")
623 _product(v, 1, v[0]);
624
625function _product(v, i=0, _tot) =
626 i>=len(v) ? _tot :
627 _product( v,
628 i+1,
629 ( is_vector(v[i])? v_mul(_tot,v[i]) : _tot*v[i] ) );
630
631
632
633// Function: cumprod()
634// Description:
635// Returns a list where each item is the cumulative product of all items up to and including the corresponding entry in the input list.
636// If passed an array of vectors, returns a list of elementwise vector products. If passed a list of square matrices returns matrix
637// products multiplying on the left, so a list `[A,B,C]` will produce the output `[A,BA,CBA]`.
638// Arguments:
639// list = The list to get the product of.
640// Example:
641// cumprod([1,3,5]); // returns [1,3,15]
642// cumprod([2,2,2]); // returns [2,4,8]
643// cumprod([[1,2,3], [3,4,5], [5,6,7]])); // returns [[1, 2, 3], [3, 8, 15], [15, 48, 105]]
644function cumprod(list) =
645 is_vector(list) ? _cumprod(list) :
646 assert(is_consistent(list), "Input must be a consistent list of scalars, vectors or square matrices")
647 is_matrix(list[0]) ? assert(len(list[0])==len(list[0][0]), "Matrices must be square") _cumprod(list)
648 : _cumprod_vec(list);
649
650function _cumprod(v,_i=0,_acc=[]) =
651 _i==len(v) ? _acc :
652 _cumprod(
653 v, _i+1,
654 concat(
655 _acc,
656 [_i==0 ? v[_i] : v[_i]*_acc[len(_acc)-1]]
657 )
658 );
659
660function _cumprod_vec(v,_i=0,_acc=[]) =
661 _i==len(v) ? _acc :
662 _cumprod_vec(
663 v, _i+1,
664 concat(
665 _acc,
666 [_i==0 ? v[_i] : v_mul(_acc[len(_acc)-1],v[_i])]
667 )
668 );
669
670
671
672// Function: convolve()
673// Usage:
674// x = convolve(p,q);
675// Description:
676// Given two vectors, or one vector and a path or
677// two paths of the same dimension, finds the convolution of them.
678// If both parameter are vectors, returns the vector convolution.
679// If one parameter is a vector and the other a path,
680// convolves using products by scalars and returns a path.
681// If both parameters are paths, convolve using scalar products
682// and returns a vector.
683// The returned vector or path has length len(p)+len(q)-1.
684// Arguments:
685// p = The first vector or path.
686// q = The second vector or path.
687// Example:
688// a = convolve([1,1],[1,2,1]); // Returns: [1,3,3,1]
689// b = convolve([1,2,3],[1,2,1])); // Returns: [1,4,8,8,3]
690// c = convolve([[1,1],[2,2],[3,1]],[1,2,1])); // Returns: [[1,1],[4,4],[8,6],[8,4],[3,1]]
691// d = convolve([[1,1],[2,2],[3,1]],[[1,2],[2,1]])); // Returns: [3,9,11,7]
692function convolve(p,q) =
693 p==[] || q==[] ? [] :
694 assert( (is_vector(p) || is_matrix(p))
695 && ( is_vector(q) || (is_matrix(q) && ( !is_vector(p[0]) || (len(p[0])==len(q[0])) ) ) ) ,
696 "The inputs should be vectors or paths all of the same dimension.")
697 let( n = len(p),
698 m = len(q))
699 [for(i=[0:n+m-2], k1 = max(0,i-n+1), k2 = min(i,m-1) )
700 sum([for(j=[k1:k2]) p[i-j]*q[j] ])
701 ];
702
703
704
705// Function: sum_of_sines()
706// Usage:
707// sum_of_sines(a,sines)
708// Description:
709// Gives the sum of a series of sines, at a given angle.
710// Arguments:
711// a = Angle to get the value for.
712// sines = List of [amplitude, frequency, offset] items, where the frequency is the number of times the cycle repeats around the circle.
713// Example:
714// v = sum_of_sines(30, [[10,3,0], [5,5.5,60]]);
715function sum_of_sines(a, sines) =
716 assert( is_finite(a) && is_matrix(sines,undef,3), "Invalid input.")
717 sum([ for (s = sines)
718 let(
719 ss=point3d(s),
720 v=ss[0]*sin(a*ss[1]+ss[2])
721 ) v
722 ]);
723
724
725
726// Section: Random Number Generation
727
728// Function: rand_int()
729// Usage:
730// rand_int(minval, maxval, n, [seed]);
731// Description:
732// Return a list of random integers in the range of minval to maxval, inclusive.
733// Arguments:
734// minval = Minimum integer value to return.
735// maxval = Maximum integer value to return.
736// N = Number of random integers to return.
737// seed = If given, sets the random number seed.
738// Example:
739// ints = rand_int(0,100,3);
740// int = rand_int(-10,10,1)[0];
741function rand_int(minval, maxval, n, seed=undef) =
742 assert( is_finite(minval+maxval+n) && (is_undef(seed) || is_finite(seed) ), "Input must be finite numbers.")
743 assert(maxval >= minval, "Max value cannot be smaller than minval")
744 let (rvect = is_def(seed) ? rands(minval,maxval+1,n,seed) : rands(minval,maxval+1,n))
745 [for(entry = rvect) floor(entry)];
746
747
748// Function: random_points()
749// Usage:
750// points = random_points(n, dim, [scale], [seed]);
751// See Also: random_polygon(), spherical_random_points()
752// Topics: Random, Points
753// Description:
754// Generate `n` uniform random points of dimension `dim` with data ranging from -scale to +scale.
755// The `scale` may be a number, in which case the random data lies in a cube,
756// or a vector with dimension `dim`, in which case each dimension has its own scale.
757// Arguments:
758// n = number of points to generate. Default: 1
759// dim = dimension of the points. Default: 2
760// scale = the scale of the point coordinates. Default: 1
761// seed = an optional seed for the random generation.
762function random_points(n, dim, scale=1, seed) =
763 assert( is_int(n) && n>=0, "The number of points should be a non-negative integer.")
764 assert( is_int(dim) && dim>=1, "The point dimensions should be an integer greater than 1.")
765 assert( is_finite(scale) || is_vector(scale,dim), "The scale should be a number or a vector with length equal to d.")
766 let(
767 rnds = is_undef(seed)
768 ? rands(-1,1,n*dim)
769 : rands(-1,1,n*dim, seed) )
770 is_num(scale)
771 ? scale*[for(i=[0:1:n-1]) [for(j=[0:dim-1]) rnds[i*dim+j] ] ]
772 : [for(i=[0:1:n-1]) [for(j=[0:dim-1]) scale[j]*rnds[i*dim+j] ] ];
773
774
775// Function: gaussian_rands()
776// Usage:
777// arr = gaussian_rands([n],[mean], [cov], [seed]);
778// Description:
779// Returns a random number or vector with a Gaussian/normal distribution.
780// Arguments:
781// n = the number of points to return. Default: 1
782// mean = The average of the random value (a number or vector). Default: 0
783// cov = covariance matrix of the random numbers, or variance in the 1D case. Default: 1
784// seed = If given, sets the random number seed.
785function gaussian_rands(n=1, mean=0, cov=1, seed=undef) =
786 assert(is_num(mean) || is_vector(mean))
787 let(
788 dim = is_num(mean) ? 1 : len(mean)
789 )
790 assert((dim==1 && is_num(cov)) || is_matrix(cov,dim,dim),"mean and covariance matrix not compatible")
791 assert(is_undef(seed) || is_finite(seed))
792 let(
793 nums = is_undef(seed)? rands(0,1,dim*n*2) : rands(0,1,dim*n*2,seed),
794 rdata = [for (i = count(dim*n,0,2)) sqrt(-2*ln(nums[i]))*cos(360*nums[i+1])]
795 )
796 dim==1 ? add_scalar(sqrt(cov)*rdata,mean) :
797 assert(is_matrix_symmetric(cov),"Supplied covariance matrix is not symmetric")
798 let(
799 L = cholesky(cov)
800 )
801 assert(is_def(L), "Supplied covariance matrix is not positive definite")
802 move(mean,list_to_matrix(rdata,dim)*transpose(L));
803
804
805// Function: exponential_rands()
806// Usage:
807// arr = exponential_rands([n], [lambda], [seed])
808// Description:
809// Returns random numbers with an exponential distribution with parameter lambda, and hence mean 1/lambda.
810// Arguments:
811// n = number of points to return. Default: 1
812// lambda = distribution parameter. The mean will be 1/lambda. Default: 1
813function exponential_rands(n=1, lambda=1, seed) =
814 assert( is_int(n) && n>=1, "The number of points should be an integer greater than zero.")
815 assert( is_num(lambda) && lambda>0, "The lambda parameter must be a positive number.")
816 let(
817 unif = is_def(seed) ? rands(0,1,n,seed=seed) : rands(0,1,n)
818 )
819 -(1/lambda) * [for(x=unif) x==1 ? 708.3964185322641 : ln(1-x)]; // Use ln(min_float) when x is 1
820
821// Function: spherical_random_points()
822// Usage:
823// points = spherical_random_points([n], [radius], [seed]);
824// See Also: random_polygon(), random_points()
825// Topics: Random, Points
826// Description:
827// Generate `n` 3D uniformly distributed random points lying on a sphere centered at the origin with radius equal to `radius`.
828// Arguments:
829// n = number of points to generate. Default: 1
830// radius = the sphere radius. Default: 1
831// seed = an optional seed for the random generation.
832
833// See https://mathworld.wolfram.com/SpherePointPicking.html
834function spherical_random_points(n=1, radius=1, seed) =
835 assert( is_int(n) && n>=1, "The number of points should be an integer greater than zero.")
836 assert( is_num(radius) && radius>0, "The radius should be a non-negative number.")
837 let( theta = is_undef(seed)
838 ? rands(0,360,n)
839 : rands(0,360,n, seed),
840 cosphi = rands(-1,1,n))
841 [for(i=[0:1:n-1]) let(
842 sin_phi=sqrt(1-cosphi[i]*cosphi[i])
843 )
844 radius*[sin_phi*cos(theta[i]),sin_phi*sin(theta[i]), cosphi[i]]];
845
846
847
848// Function: random_polygon()
849// Usage:
850// points = random_polygon([n], [size], [seed]);
851// See Also: random_points(), spherical_random_points()
852// Topics: Random, Polygon
853// Description:
854// Generate the `n` vertices of a random counter-clockwise simple 2d polygon
855// inside a circle centered at the origin with radius `size`.
856// Arguments:
857// n = number of vertices of the polygon. Default: 3
858// size = the radius of a circle centered at the origin containing the polygon. Default: 1
859// seed = an optional seed for the random generation.
860function random_polygon(n=3,size=1, seed) =
861 assert( is_int(n) && n>2, "Improper number of polygon vertices.")
862 assert( is_num(size) && size>0, "Improper size.")
863 let(
864 seed = is_undef(seed) ? rands(0,1,1)[0] : seed,
865 cumm = cumsum(rands(0.1,10,n+1,seed)),
866 angs = 360*cumm/cumm[n-1],
867 rads = rands(.01,size,n,seed+cumm[0])
868 )
869 [for(i=count(n)) rads[i]*[cos(angs[i]), sin(angs[i])] ];
870
871
872
873// Section: Calculus
874
875// Function: deriv()
876// Usage:
877// x = deriv(data, [h], [closed])
878// Description:
879// Computes a numerical derivative estimate of the data, which may be scalar or vector valued.
880// The `h` parameter gives the step size of your sampling so the derivative can be scaled correctly.
881// If the `closed` parameter is true the data is assumed to be defined on a loop with data[0] adjacent to
882// data[len(data)-1]. This function uses a symetric derivative approximation
883// for internal points, f'(t) = (f(t+h)-f(t-h))/2h. For the endpoints (when closed=false) the algorithm
884// 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.
885// .
886// If `h` is a vector then it is assumed to be nonuniform, with h[i] giving the sampling distance
887// between data[i+1] and data[i], and the data values will be linearly resampled at each corner
888// to produce a uniform spacing for the derivative estimate. At the endpoints a single point method
889// is used: f'(t) = (f(t+h)-f(t))/h.
890// Arguments:
891// data = the list of the elements to compute the derivative of.
892// h = the parametric sampling of the data.
893// closed = boolean to indicate if the data set should be wrapped around from the end to the start.
894function deriv(data, h=1, closed=false) =
895 assert( is_consistent(data) , "Input list is not consistent or not numerical.")
896 assert( len(data)>=2, "Input `data` should have at least 2 elements.")
897 assert( is_finite(h) || is_vector(h), "The sampling `h` must be a number or a list of numbers." )
898 assert( is_num(h) || len(h) == len(data)-(closed?0:1),
899 str("Vector valued `h` must have length ",len(data)-(closed?0:1)))
900 is_vector(h) ? _deriv_nonuniform(data, h, closed=closed) :
901 let( L = len(data) )
902 closed
903 ? [
904 for(i=[0:1:L-1])
905 (data[(i+1)%L]-data[(L+i-1)%L])/2/h
906 ]
907 : let(
908 first = L<3 ? data[1]-data[0] :
909 3*(data[1]-data[0]) - (data[2]-data[1]),
910 last = L<3 ? data[L-1]-data[L-2]:
911 (data[L-3]-data[L-2])-3*(data[L-2]-data[L-1])
912 )
913 [
914 first/2/h,
915 for(i=[1:1:L-2]) (data[i+1]-data[i-1])/2/h,
916 last/2/h
917 ];
918
919
920function _dnu_calc(f1,fc,f2,h1,h2) =
921 let(
922 f1 = h2<h1 ? lerp(fc,f1,h2/h1) : f1 ,
923 f2 = h1<h2 ? lerp(fc,f2,h1/h2) : f2
924 )
925 (f2-f1) / 2 / min(h1,h2);
926
927
928function _deriv_nonuniform(data, h, closed) =
929 let( L = len(data) )
930 closed
931 ? [for(i=[0:1:L-1])
932 _dnu_calc(data[(L+i-1)%L], data[i], data[(i+1)%L], select(h,i-1), h[i]) ]
933 : [
934 (data[1]-data[0])/h[0],
935 for(i=[1:1:L-2]) _dnu_calc(data[i-1],data[i],data[i+1], h[i-1],h[i]),
936 (data[L-1]-data[L-2])/h[L-2]
937 ];
938
939
940// Function: deriv2()
941// Usage:
942// x = deriv2(data, [h], [closed])
943// Description:
944// Computes a numerical estimate of the second derivative of the data, which may be scalar or vector valued.
945// The `h` parameter gives the step size of your sampling so the derivative can be scaled correctly.
946// If the `closed` parameter is true the data is assumed to be defined on a loop with data[0] adjacent to
947// data[len(data)-1]. For internal points this function uses the approximation
948// f''(t) = (f(t-h)-2*f(t)+f(t+h))/h^2. For the endpoints (when closed=false),
949// when sufficient points are available, the method is either the four point expression
950// f''(t) = (2*f(t) - 5*f(t+h) + 4*f(t+2*h) - f(t+3*h))/h^2 or
951// 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
952// if five points are available.
953// Arguments:
954// data = the list of the elements to compute the derivative of.
955// h = the constant parametric sampling of the data.
956// closed = boolean to indicate if the data set should be wrapped around from the end to the start.
957function deriv2(data, h=1, closed=false) =
958 assert( is_consistent(data) , "Input list is not consistent or not numerical.")
959 assert( is_finite(h), "The sampling `h` must be a number." )
960 let( L = len(data) )
961 assert( L>=3, "Input list has less than 3 elements.")
962 closed
963 ? [
964 for(i=[0:1:L-1])
965 (data[(i+1)%L]-2*data[i]+data[(L+i-1)%L])/h/h
966 ]
967 :
968 let(
969 first =
970 L==3? data[0] - 2*data[1] + data[2] :
971 L==4? 2*data[0] - 5*data[1] + 4*data[2] - data[3] :
972 (35*data[0] - 104*data[1] + 114*data[2] - 56*data[3] + 11*data[4])/12,
973 last =
974 L==3? data[L-1] - 2*data[L-2] + data[L-3] :
975 L==4? -2*data[L-1] + 5*data[L-2] - 4*data[L-3] + data[L-4] :
976 (35*data[L-1] - 104*data[L-2] + 114*data[L-3] - 56*data[L-4] + 11*data[L-5])/12
977 ) [
978 first/h/h,
979 for(i=[1:1:L-2]) (data[i+1]-2*data[i]+data[i-1])/h/h,
980 last/h/h
981 ];
982
983
984// Function: deriv3()
985// Usage:
986// x = deriv3(data, [h], [closed])
987// Description:
988// Computes a numerical third derivative estimate of the data, which may be scalar or vector valued.
989// The `h` parameter gives the step size of your sampling so the derivative can be scaled correctly.
990// If the `closed` parameter is true the data is assumed to be defined on a loop with data[0] adjacent to
991// data[len(data)-1]. This function uses a five point derivative estimate, so the input data must include
992// at least five points:
993// 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
994// 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
995// f'''(t) = (-3*f(t-h)+10*f(t)-12*f(t+h)+6*f(t+2*h)-f(t+3*h)) / 2h^3.
996// Arguments:
997// data = the list of the elements to compute the derivative of.
998// h = the constant parametric sampling of the data.
999// closed = boolean to indicate if the data set should be wrapped around from the end to the start.
1000function deriv3(data, h=1, closed=false) =
1001 assert( is_consistent(data) , "Input list is not consistent or not numerical.")
1002 assert( len(data)>=5, "Input list has less than 5 elements.")
1003 assert( is_finite(h), "The sampling `h` must be a number." )
1004 let(
1005 L = len(data),
1006 h3 = h*h*h
1007 )
1008 closed? [
1009 for(i=[0:1:L-1])
1010 (-data[(L+i-2)%L]+2*data[(L+i-1)%L]-2*data[(i+1)%L]+data[(i+2)%L])/2/h3
1011 ] :
1012 let(
1013 first=(-5*data[0]+18*data[1]-24*data[2]+14*data[3]-3*data[4])/2,
1014 second=(-3*data[0]+10*data[1]-12*data[2]+6*data[3]-data[4])/2,
1015 last=(5*data[L-1]-18*data[L-2]+24*data[L-3]-14*data[L-4]+3*data[L-5])/2,
1016 prelast=(3*data[L-1]-10*data[L-2]+12*data[L-3]-6*data[L-4]+data[L-5])/2
1017 ) [
1018 first/h3,
1019 second/h3,
1020 for(i=[2:1:L-3]) (-data[i-2]+2*data[i-1]-2*data[i+1]+data[i+2])/2/h3,
1021 prelast/h3,
1022 last/h3
1023 ];
1024
1025
1026// Section: Complex Numbers
1027
1028
1029// Function: complex()
1030// Usage:
1031// z = complex(list)
1032// Description:
1033// Converts a real valued number, vector or matrix into its complex analog
1034// by replacing all entries with a 2-vector that has zero imaginary part.
1035function complex(list) =
1036 is_num(list) ? [list,0] :
1037 [for(entry=list) is_num(entry) ? [entry,0] : complex(entry)];
1038
1039
1040// Function: c_mul()
1041// Usage:
1042// c = c_mul(z1,z2)
1043// Description:
1044// Multiplies two complex numbers, vectors or matrices, where complex numbers
1045// or entries are represented as vectors: [REAL, IMAGINARY]. Note that all
1046// entries in both arguments must be complex.
1047// Arguments:
1048// z1 = First complex number, vector or matrix
1049// z2 = Second complex number, vector or matrix
1050function c_mul(z1,z2) =
1051 is_matrix([z1,z2],2,2) ? _c_mul(z1,z2) :
1052 _combine_complex(_c_mul(_split_complex(z1), _split_complex(z2)));
1053
1054
1055function _split_complex(data) =
1056 is_vector(data,2) ? data
1057 : is_num(data[0][0]) ? [data*[1,0], data*[0,1]]
1058 : [
1059 [for(vec=data) vec * [1,0]],
1060 [for(vec=data) vec * [0,1]]
1061 ];
1062
1063
1064function _combine_complex(data) =
1065 is_vector(data,2) ? data
1066 : is_num(data[0][0]) ? [for(i=[0:len(data[0])-1]) [data[0][i],data[1][i]]]
1067 : [for(i=[0:1:len(data[0])-1])
1068 [for(j=[0:1:len(data[0][0])-1])
1069 [data[0][i][j], data[1][i][j]]]];
1070
1071
1072function _c_mul(z1,z2) =
1073 [ z1.x*z2.x - z1.y*z2.y, z1.x*z2.y + z1.y*z2.x ];
1074
1075
1076// Function: c_div()
1077// Usage:
1078// x = c_div(z1,z2)
1079// Description:
1080// Divides two complex numbers represented by 2D vectors.
1081// Returns a complex number as a 2D vector [REAL, IMAGINARY].
1082// Arguments:
1083// z1 = First complex number, given as a 2D vector [REAL, IMAGINARY]
1084// z2 = Second complex number, given as a 2D vector [REAL, IMAGINARY]
1085function c_div(z1,z2) =
1086 assert( is_vector(z1,2) && is_vector(z2), "Complex numbers should be represented by 2D vectors." )
1087 assert( !approx(z2,0), "The divisor `z2` cannot be zero." )
1088 let(den = z2.x*z2.x + z2.y*z2.y)
1089 [(z1.x*z2.x + z1.y*z2.y)/den, (z1.y*z2.x - z1.x*z2.y)/den];
1090
1091
1092// Function: c_conj()
1093// Usage:
1094// w = c_conj(z)
1095// Description:
1096// Computes the complex conjugate of the input, which can be a complex number,
1097// complex vector or complex matrix.
1098function c_conj(z) =
1099 is_vector(z,2) ? [z.x,-z.y] :
1100 [for(entry=z) c_conj(entry)];
1101
1102
1103// Function: c_real()
1104// Usage:
1105// x = c_real(z)
1106// Description:
1107// Returns real part of a complex number, vector or matrix.
1108function c_real(z) =
1109 is_vector(z,2) ? z.x
1110 : is_num(z[0][0]) ? z*[1,0]
1111 : [for(vec=z) vec * [1,0]];
1112
1113
1114// Function: c_imag()
1115// Usage:
1116// x = c_imag(z)
1117// Description:
1118// Returns imaginary part of a complex number, vector or matrix.
1119function c_imag(z) =
1120 is_vector(z,2) ? z.y
1121 : is_num(z[0][0]) ? z*[0,1]
1122 : [for(vec=z) vec * [0,1]];
1123
1124
1125// Function: c_ident()
1126// Usage:
1127// I = c_ident(n)
1128// Description:
1129// Produce an n by n complex identity matrix
1130function c_ident(n) = [for (i = [0:1:n-1]) [for (j = [0:1:n-1]) (i==j)?[1,0]:[0,0]]];
1131
1132
1133// Function: c_norm()
1134// Usage:
1135// n = c_norm(z)
1136// Description:
1137// Compute the norm of a complex number or vector.
1138function c_norm(z) = norm_fro(z);
1139
1140
1141// Section: Polynomials
1142
1143// Function: quadratic_roots()
1144// Usage:
1145// roots = quadratic_roots(a, b, c, [real])
1146// Description:
1147// Computes roots of the quadratic equation a*x^2+b*x+c==0, where the
1148// coefficients are real numbers. If real is true then returns only the
1149// real roots. Otherwise returns a pair of complex values. This method
1150// may be more reliable than the general root finder at distinguishing
1151// real roots from complex roots.
1152// Algorithm from: https://people.csail.mit.edu/bkph/articles/Quadratics.pdf
1153function quadratic_roots(a,b,c,real=false) =
1154 real ? [for(root = quadratic_roots(a,b,c,real=false)) if (root.y==0) root.x]
1155 :
1156 is_undef(b) && is_undef(c) && is_vector(a,3) ? quadratic_roots(a[0],a[1],a[2]) :
1157 assert(is_num(a) && is_num(b) && is_num(c))
1158 assert(a!=0 || b!=0 || c!=0, "Quadratic must have a nonzero coefficient")
1159 a==0 && b==0 ? [] : // No solutions
1160 a==0 ? [[-c/b,0]] :
1161 let(
1162 descrim = b*b-4*a*c,
1163 sqrt_des = sqrt(abs(descrim))
1164 )
1165 descrim < 0 ? // Complex case
1166 [[-b, sqrt_des],
1167 [-b, -sqrt_des]]/2/a :
1168 b<0 ? // b positive
1169 [[2*c/(-b+sqrt_des),0],
1170 [(-b+sqrt_des)/a/2,0]]
1171 : // b negative
1172 [[(-b-sqrt_des)/2/a, 0],
1173 [2*c/(-b-sqrt_des),0]];
1174
1175
1176// Function: polynomial()
1177// Usage:
1178// x = polynomial(p, z)
1179// Description:
1180// Evaluates specified real polynomial, p, at the complex or real input value, z.
1181// The polynomial is specified as p=[a_n, a_{n-1},...,a_1,a_0]
1182// where a_n is the z^n coefficient. Polynomial coefficients are real.
1183// The result is a number if `z` is a number and a complex number otherwise.
1184function polynomial(p,z,k,total) =
1185 is_undef(k)
1186 ? assert( is_vector(p) , "Input polynomial coefficients must be a vector." )
1187 assert( is_finite(z) || is_vector(z,2), "The value of `z` must be a real or a complex number." )
1188 polynomial( _poly_trim(p), z, 0, is_num(z) ? 0 : [0,0])
1189 : k==len(p) ? total
1190 : polynomial(p,z,k+1, is_num(z) ? total*z+p[k] : c_mul(total,z)+[p[k],0]);
1191
1192
1193// Function: poly_mult()
1194// Usage:
1195// x = polymult(p,q)
1196// x = polymult([p1,p2,p3,...])
1197// Description:
1198// Given a list of polynomials represented as real algebraic coefficient lists, with the highest degree coefficient first,
1199// computes the coefficient list of the product polynomial.
1200function poly_mult(p,q) =
1201 is_undef(q) ?
1202 len(p)==2
1203 ? poly_mult(p[0],p[1])
1204 : poly_mult(p[0], poly_mult(list_tail(p)))
1205 :
1206 assert( is_vector(p) && is_vector(q),"Invalid arguments to poly_mult")
1207 p*p==0 || q*q==0
1208 ? [0]
1209 : _poly_trim(convolve(p,q));
1210
1211
1212// Function: poly_div()
1213// Usage:
1214// [quotient,remainder] = poly_div(n,d)
1215// Description:
1216// Computes division of the numerator polynomial by the denominator polynomial and returns
1217// a list of two polynomials, [quotient, remainder]. If the division has no remainder then
1218// the zero polynomial [0] is returned for the remainder. Similarly if the quotient is zero
1219// the returned quotient will be [0].
1220function poly_div(n,d) =
1221 assert( is_vector(n) && is_vector(d) , "Invalid polynomials." )
1222 let( d = _poly_trim(d),
1223 n = _poly_trim(n) )
1224 assert( d!=[0] , "Denominator cannot be a zero polynomial." )
1225 n==[0]
1226 ? [[0],[0]]
1227 : _poly_div(n,d,q=[]);
1228
1229function _poly_div(n,d,q) =
1230 len(n)<len(d) ? [q,_poly_trim(n)] :
1231 let(
1232 t = n[0] / d[0],
1233 newq = concat(q,[t]),
1234 newn = [for(i=[1:1:len(n)-1]) i<len(d) ? n[i] - t*d[i] : n[i]]
1235 )
1236 _poly_div(newn,d,newq);
1237
1238
1239/// Internal Function: _poly_trim()
1240/// Usage:
1241/// _poly_trim(p, [eps])
1242/// Description:
1243/// Removes leading zero terms of a polynomial. By default zeros must be exact,
1244/// or give epsilon for approximate zeros. Returns [0] for a zero polynomial.
1245function _poly_trim(p,eps=0) =
1246 let( nz = [for(i=[0:1:len(p)-1]) if ( !approx(p[i],0,eps)) i])
1247 len(nz)==0 ? [0] : list_tail(p,nz[0]);
1248
1249
1250// Function: poly_add()
1251// Usage:
1252// sum = poly_add(p,q)
1253// Description:
1254// Computes the sum of two polynomials.
1255function poly_add(p,q) =
1256 assert( is_vector(p) && is_vector(q), "Invalid input polynomial(s)." )
1257 let( plen = len(p),
1258 qlen = len(q),
1259 long = plen>qlen ? p : q,
1260 short = plen>qlen ? q : p
1261 )
1262 _poly_trim(long + concat(repeat(0,len(long)-len(short)),short));
1263
1264
1265// Function: poly_roots()
1266// Usage:
1267// roots = poly_roots(p, [tol]);
1268// Description:
1269// Returns all complex roots of the specified real polynomial p.
1270// The polynomial is specified as p=[a_n, a_{n-1},...,a_1,a_0]
1271// where a_n is the z^n coefficient. The tol parameter gives
1272// the stopping tolerance for the iteration. The polynomial
1273// must have at least one non-zero coefficient. Convergence is poor
1274// if the polynomial has any repeated roots other than zero.
1275// Arguments:
1276// p = polynomial coefficients with higest power coefficient first
1277// tol = tolerance for iteration. Default: 1e-14
1278
1279// Uses the Aberth method https://en.wikipedia.org/wiki/Aberth_method
1280//
1281// Dario Bini. "Numerical computation of polynomial zeros by means of Aberth's Method", Numerical Algorithms, Feb 1996.
1282// https://www.researchgate.net/publication/225654837_Numerical_computation_of_polynomial_zeros_by_means_of_Aberth's_method
1283function poly_roots(p,tol=1e-14,error_bound=false) =
1284 assert( is_vector(p), "Invalid polynomial." )
1285 let( p = _poly_trim(p,eps=0) )
1286 assert( p!=[0], "Input polynomial cannot be zero." )
1287 p[len(p)-1] == 0 ? // Strip trailing zero coefficients
1288 let( solutions = poly_roots(list_head(p),tol=tol, error_bound=error_bound))
1289 (error_bound ? [ [[0,0], each solutions[0]], [0, each solutions[1]]]
1290 : [[0,0], each solutions]) :
1291 len(p)==1 ? (error_bound ? [[],[]] : []) : // Nonzero constant case has no solutions
1292 len(p)==2 ? let( solution = [[-p[1]/p[0],0]]) // Linear case needs special handling
1293 (error_bound ? [solution,[0]] : solution)
1294 :
1295 let(
1296 n = len(p)-1, // polynomial degree
1297 pderiv = [for(i=[0:n-1]) p[i]*(n-i)],
1298
1299 s = [for(i=[0:1:n]) abs(p[i])*(4*(n-i)+1)], // Error bound polynomial from Bini
1300
1301 // Using method from: http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/0915-24.pdf
1302 beta = -p[1]/p[0]/n,
1303 r = 1+pow(abs(polynomial(p,beta)/p[0]),1/n),
1304 init = [for(i=[0:1:n-1]) // Initial guess for roots
1305 let(angle = 360*i/n+270/n/PI)
1306 [beta,0]+r*[cos(angle),sin(angle)]
1307 ],
1308 roots = _poly_roots(p,pderiv,s,init,tol=tol),
1309 error = error_bound ? [for(xi=roots) n * (norm(polynomial(p,xi))+tol*polynomial(s,norm(xi))) /
1310 abs(norm(polynomial(pderiv,xi))-tol*polynomial(s,norm(xi)))] : 0
1311 )
1312 error_bound ? [roots, error] : roots;
1313
1314// Internal function
1315// p = polynomial
1316// pderiv = derivative polynomial of p
1317// z = current guess for the roots
1318// tol = root tolerance
1319// i=iteration counter
1320function _poly_roots(p, pderiv, s, z, tol, i=0) =
1321 assert(i<45, str("Polyroot exceeded iteration limit. Current solution:", z))
1322 let(
1323 n = len(z),
1324 svals = [for(zk=z) tol*polynomial(s,norm(zk))],
1325 p_of_z = [for(zk=z) polynomial(p,zk)],
1326 done = [for(k=[0:n-1]) norm(p_of_z[k])<=svals[k]],
1327 newton = [for(k=[0:n-1]) c_div(p_of_z[k], polynomial(pderiv,z[k]))],
1328 zdiff = [for(k=[0:n-1]) sum([for(j=[0:n-1]) if (j!=k) c_div([1,0], z[k]-z[j])])],
1329 w = [for(k=[0:n-1]) done[k] ? [0,0] : c_div( newton[k],
1330 [1,0] - c_mul(newton[k], zdiff[k]))]
1331 )
1332 all(done) ? z : _poly_roots(p,pderiv,s,z-w,tol,i+1);
1333
1334
1335// Function: real_roots()
1336// Usage:
1337// roots = real_roots(p, [eps], [tol])
1338// Description:
1339// Returns the real roots of the specified real polynomial p.
1340// The polynomial is specified as p=[a_n, a_{n-1},...,a_1,a_0]
1341// where a_n is the x^n coefficient. This function works by
1342// computing the complex roots and returning those roots where
1343// the imaginary part is closed to zero. By default it uses a computed
1344// error bound from the polynomial solver to decide whether imaginary
1345// parts are zero. You can specify eps, in which case the test is
1346// z.y/(1+norm(z)) < eps. Because
1347// of poor convergence and higher error for repeated roots, such roots may
1348// be missed by the algorithm because their imaginary part is large.
1349// Arguments:
1350// p = polynomial to solve as coefficient list, highest power term first
1351// eps = used to determine whether imaginary parts of roots are zero
1352// tol = tolerance for the complex polynomial root finder
1353
1354// The algorithm is based on Brent's method and is a combination of
1355// bisection and inverse quadratic approximation, where bisection occurs
1356// at every step, with refinement using inverse quadratic approximation
1357// only when that approximation gives a good result. The detail
1358// of how to decide when to use the quadratic came from an article
1359// by Crenshaw on "The World's Best Root Finder".
1360// https://www.embedded.com/worlds-best-root-finder/
1361function real_roots(p,eps=undef,tol=1e-14) =
1362 assert( is_vector(p), "Invalid polynomial." )
1363 let( p = _poly_trim(p,eps=0) )
1364 assert( p!=[0], "Input polynomial cannot be zero." )
1365 let(
1366 roots_err = poly_roots(p,error_bound=true),
1367 roots = roots_err[0],
1368 err = roots_err[1]
1369 )
1370 is_def(eps)
1371 ? [for(z=roots) if (abs(z.y)/(1+norm(z))<eps) z.x]
1372 : [for(i=idx(roots)) if (abs(roots[i].y)<=err[i]) roots[i].x];
1373
1374
1375// Section: Operations on Functions
1376
1377// Function: root_find()
1378// Usage:
1379// x = root_find(f, x0, x1, [tol])
1380// Description:
1381// Find a root of the continuous function f where the sign of f(x0) is different
1382// from the sign of f(x1). The function f is a function literal accepting one
1383// argument. You must have a version of OpenSCAD that supports function literals
1384// (2021.01 or newer). The tolerance (tol) specifies the accuracy of the solution:
1385// abs(f(x)) < tol * yrange, where yrange is the range of observed function values.
1386// This function can only find roots that cross the x axis: it cannot find the
1387// the root of x^2.
1388// Arguments:
1389// f = function literal for a scalar-valued single variable function
1390// x0 = endpoint of interval to search for root
1391// x1 = second endpoint of interval to search for root
1392// tol = tolerance for solution. Default: 1e-15
1393function root_find(f,x0,x1,tol=1e-15) =
1394 let(
1395 y0 = f(x0),
1396 y1 = f(x1),
1397 yrange = y0<y1 ? [y0,y1] : [y1,y0]
1398 )
1399 // Check endpoints
1400 y0==0 || _rfcheck(x0, y0,yrange,tol) ? x0 :
1401 y1==0 || _rfcheck(x1, y1,yrange,tol) ? x1 :
1402 assert(y0*y1<0, "Sign of function must be different at the interval endpoints")
1403 _rootfind(f,[x0,x1],[y0,y1],yrange,tol);
1404
1405function _rfcheck(x,y,range,tol) =
1406 assert(is_finite(y), str("Function not finite at ",x))
1407 abs(y) < tol*(range[1]-range[0]);
1408
1409// xpts and ypts are arrays whose first two entries contain the
1410// interval bracketing the root. Extra entries are ignored.
1411// yrange is the total observed range of y values (used for the
1412// tolerance test).
1413function _rootfind(f, xpts, ypts, yrange, tol, i=0) =
1414 assert(i<100, "root_find did not converge to a solution")
1415 let(
1416 xmid = (xpts[0]+xpts[1])/2,
1417 ymid = f(xmid),
1418 yrange = [min(ymid, yrange[0]), max(ymid, yrange[1])]
1419 )
1420 _rfcheck(xmid, ymid, yrange, tol) ? xmid :
1421 let(
1422 // Force root to be between x0 and midpoint
1423 y = ymid * ypts[0] < 0 ? [ypts[0], ymid, ypts[1]]
1424 : [ypts[1], ymid, ypts[0]],
1425 x = ymid * ypts[0] < 0 ? [xpts[0], xmid, xpts[1]]
1426 : [xpts[1], xmid, xpts[0]],
1427 v = y[2]*(y[2]-y[0]) - 2*y[1]*(y[1]-y[0])
1428 )
1429 v <= 0 ? _rootfind(f,x,y,yrange,tol,i+1) // Root is between first two points, extra 3rd point doesn't hurt
1430 :
1431 let( // Do quadratic approximation
1432 B = (x[1]-x[0]) / (y[1]-y[0]),
1433 C = y*[-1,2,-1] / (y[2]-y[1]) / (y[2]-y[0]),
1434 newx = x[0] - B * y[0] *(1-C*y[1]),
1435 newy = f(newx),
1436 new_yrange = [min(yrange[0],newy), max(yrange[1], newy)],
1437 // select interval that contains the root by checking sign
1438 yinterval = newy*y[0] < 0 ? [y[0],newy] : [newy,y[1]],
1439 xinterval = newy*y[0] < 0 ? [x[0],newx] : [newx,x[1]]
1440 )
1441 _rfcheck(newx, newy, new_yrange, tol)
1442 ? newx
1443 : _rootfind(f, xinterval, yinterval, new_yrange, tol, i+1);
1444
1445
1446
1447// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap