Beast - Music Synthesizer and Composer  0.11.1+10.g2da35
bsemath.hh
Go to the documentation of this file.
1  // Licensed GNU LGPL v2.1 or later: http://www.gnu.org/licenses/lgpl.html
2 #ifndef __BSE_MATH_H__
3 #define __BSE_MATH_H__
4 
5 #include <bse/bsedefs.hh>
6 #include <bse/bseieee754.hh> /* provides math.h */
7 
8 /* --- constants --- */
9 /* PI is defined in bseieee754.hh */
10 #define BSE_1_DIV_PI (0.3183098861837906715377675267450287240689) // 1/pi
11 #define BSE_PI_DIV_2 (1.570796326794896619231321691639751442099) // pi/2
12 #define BSE_2_DIV_PI (0.6366197723675813430755350534900574481378) // 2/pi
13 #define BSE_2_DIV_SQRT_PI (1.128379167095512573896158903121545171688) // 2/sqrt(pi)
14 #define BSE_2_DIV_LN2 (2.88539008177792681471984936200378427485329) // 2/ln(2)
15 #define BSE_PI_DIV_4 (0.7853981633974483096156608458198757210493) // pi/4
16 #define BSE_E (2.718281828459045235360287471352662497757) // e^1
17 #define BSE_LOG2E (1.442695040888963407359924681001892137427) // log_2(e^1)
18 #define BSE_LOG10E (0.4342944819032518276511289189166050822944) // log_10(e^1)
19 #define BSE_LN2 (0.6931471805599453094172321214581765680755) // ln(2)
20 #define BSE_SQRT2 (1.41421356237309504880168872420969807857) // sqrt(2)
21 #define BSE_1_DIV_SQRT2 (0.7071067811865475244008443621048490392848) // 1/sqrt(2)
22 #define BSE_LN4 (1.386294361119890618834464242916353136151) // ln(4)
23 #define BSE_LN10 (2.302585092994045684017991454684364207601) // ln(10)
24 #define BSE_LOG2_10 (3.321928094887362347870319429489390175865) // log_2(10)
25 #define BSE_LOG2POW20_10 (0.1660964047443681173935159714744695087932) // log_2(10)/20
26 #define BSE_2_POW_1_DIV_12 (1.059463094359295264561825294946341700779) // 2^(1/12)
27 #define BSE_LN_2_POW_1_DIV_12 (5.776226504666210911810267678818138067296e-2) // ln(2^(1/12))
28 #define BSE_LN_2_POW_1_DIV_1200_d (5.776226504666210911810267678818138067296e-4) // ln(2^(1/1200))
29 #define BSE_2_POW_1_DIV_72 (1.009673533228510862192521401118605073603) // 2^(1/72)
30 #define BSE_LN_2_POW_1_DIV_72 (9.62704417444368485301711279803023011216e-3) // ln(2^(1/72))
31 #define BSE_DECIBEL20_FACTOR (8.68588963806503655302257837833210164588794) // 20.0 / ln (10.0)
32 #define BSE_DECIBEL10_FACTOR (4.34294481903251827651128918916605082294397) // 10.0 / ln (10.0)
33 #define BSE_1_DIV_DECIBEL20_FACTOR (0.1151292546497022842008995727342182103801) // ln (10) / 20
34 #define BSE_COMPLEX_ONE (bse_complex (1, 0))
35 
36 /* --- structures --- */
37 typedef struct {
38  double re;
39  double im;
40 } BseComplex;
41 
42 /* --- complex numbers --- */
43 static inline BseComplex bse_complex (double re,
44  double im);
45 static inline BseComplex bse_complex_polar (double abs,
46  double arg);
47 static inline BseComplex bse_complex_add (BseComplex c1,
48  BseComplex c2);
49 static inline BseComplex bse_complex_add3 (BseComplex c1,
50  BseComplex c2,
51  BseComplex c3);
52 static inline BseComplex bse_complex_sub (BseComplex c1,
53  BseComplex c2);
54 static inline BseComplex bse_complex_sub3 (BseComplex c1,
55  BseComplex c2,
56  BseComplex c3);
57 static inline BseComplex bse_complex_scale (BseComplex c1,
58  double scale);
59 static inline BseComplex bse_complex_mul (BseComplex c1,
60  BseComplex c2);
61 static inline BseComplex bse_complex_mul3 (BseComplex c1,
62  BseComplex c2,
63  BseComplex c3);
64 static inline BseComplex bse_complex_div (BseComplex a,
65  BseComplex b);
66 static inline BseComplex bse_complex_reciprocal (BseComplex c);
67 static inline BseComplex bse_complex_sqrt (BseComplex z);
68 static inline BseComplex bse_complex_conj (BseComplex c); /* {re, -im} */
69 static inline BseComplex bse_complex_id (BseComplex c);
70 static inline BseComplex bse_complex_inv (BseComplex c); /* {-re, -im} */
71 static inline double bse_complex_abs (BseComplex c);
72 static inline double bse_complex_arg (BseComplex c);
73 static inline BseComplex bse_complex_sin (BseComplex c);
74 static inline BseComplex bse_complex_cos (BseComplex c);
75 static inline BseComplex bse_complex_tan (BseComplex c);
76 static inline BseComplex bse_complex_sinh (BseComplex c);
77 static inline BseComplex bse_complex_cosh (BseComplex c);
78 static inline BseComplex bse_complex_tanh (BseComplex c);
79 std::string bse_complex_str (BseComplex c);
80 std::string bse_complex_list (uint n_points,
81  BseComplex *points,
82  const std::string &indent);
83 void bse_complex_gnuplot (const char *file_name,
84  uint n_points,
85  BseComplex *points);
86 std::string bse_string_from_double (long double value);
87 
88 /* --- polynomials --- */
89 /* example, degree=2: 5+2x+7x^2 => a[0..degree] = { 5, 2, 7 } */
90 static inline void bse_poly_add (uint degree,
91  double *a, /* a[0..degree] */
92  double *b);
93 static inline void bse_poly_sub (uint order,
94  double *a, /* [0..degree] */
95  double *b);
96 static inline void bse_poly_mul (double *p, /* out:[0..aorder+border] */
97  uint aorder,
98  const double *a, /* in:[0..aorder] */
99  uint border,
100  const double *b); /* in:[0..border] */
101 static inline void bse_poly_scale (uint order,
102  double *a, /* [0..degree] */
103  double scale);
104 static inline void bse_poly_xscale (uint order,
105  double *a, /* [0..degree] */
106  double xscale);
107 static inline double bse_poly_eval (uint degree,
108  double *a, /* [0..degree] */
109  double x);
110 void bse_poly_complex_roots (uint poly_degree,
111  double *a, /* [0..degree] (degree+1 elements) */
112  BseComplex *roots); /* [degree] */
113 void bse_poly_from_re_roots (uint poly_degree,
114  double *a, /* [0..degree] */
115  BseComplex *roots);
116 void bse_cpoly_from_roots (uint poly_degree,
117  BseComplex *c, /* [0..degree] */
118  BseComplex *roots);
119 static inline void bse_cpoly_mul_monomial (uint degree, /* _new_ degree */
120  BseComplex *c, /* in:[0..degree-1] out:[0..degree] */
121  BseComplex root); /* c(x) *= (x^1 - root) */
122 static inline void bse_cpoly_mul_reciprocal (uint degree, /* _new_ degree */
123  BseComplex *c, /* in:[0..degree-1] out:[0..degree] */
124  BseComplex root); /* c(x) *= (1 - root * x^-1) */
125 static inline void bse_cpoly_mul (BseComplex *p, /* out:[0..aorder+border] */
126  uint aorder,
127  BseComplex *a, /* in:[0..aorder] */
128  uint border,
129  BseComplex *b); /* in:[0..border] */
130 gboolean bse_poly2_droots (gdouble roots[2],
131  gdouble a,
132  gdouble b,
133  gdouble c);
134 std::string bse_poly_str (uint degree,
135  double *a,
136  const std::string &var);
137 std::string bse_poly_str1 (uint degree,
138  double *a,
139  const std::string &var);
140 
141 /* --- transformations --- */
142 double bse_temp_freq (double kammer_freq,
143  int semitone_delta);
144 
145 /* --- miscellaneous --- */
146 double bse_bit_depth_epsilon (guint n_bits); /* 1..32 */
147 gint bse_rand_int (void); /* +-G_MAXINT */
148 gfloat bse_rand_float (void); /* -1.0..1.0 */
149 gint bse_rand_bool (void); /* random bit */
150 void bse_float_gnuplot (const char *file_name,
151  double xstart,
152  double xstep,
153  uint n_ypoints,
154  const float *ypoints);
155 
156 
157 /* --- implementations --- */
158 static inline BseComplex
159 bse_complex (double re,
160  double im)
161 {
162  BseComplex r;
163  r.re = re;
164  r.im = im;
165  return r;
166 }
167 static inline BseComplex
168 bse_complex_polar (double abs,
169  double arg)
170 {
171  return bse_complex (abs * cos (arg), abs * sin (arg));
172 }
173 static inline BseComplex
174 bse_complex_add (BseComplex c1,
175  BseComplex c2)
176 {
177  return bse_complex (c1.re + c2.re, c1.im + c2.im);
178 }
179 static inline BseComplex
180 bse_complex_add3 (BseComplex c1,
181  BseComplex c2,
182  BseComplex c3)
183 {
184  return bse_complex (c1.re + c2.re + c3.re, c1.im + c2.im + c3.im);
185 }
186 static inline BseComplex
187 bse_complex_sub (BseComplex c1,
188  BseComplex c2)
189 {
190  return bse_complex (c1.re - c2.re, c1.im - c2.im);
191 }
192 static inline BseComplex
193 bse_complex_sub3 (BseComplex c1,
194  BseComplex c2,
195  BseComplex c3)
196 {
197  return bse_complex (c1.re - c2.re - c3.re, c1.im - c2.im - c3.im);
198 }
199 static inline BseComplex
200 bse_complex_scale (BseComplex c1,
201  double scale)
202 {
203  return bse_complex (c1.re * scale, c1.im * scale);
204 }
205 static inline BseComplex
206 bse_complex_mul (BseComplex c1,
207  BseComplex c2)
208 {
209  return bse_complex (c1.re * c2.re - c1.im * c2.im, c1.re * c2.im + c1.im * c2.re);
210 }
211 static inline BseComplex
212 bse_complex_mul3 (BseComplex c1,
213  BseComplex c2,
214  BseComplex c3)
215 {
216  double aec = c1.re * c2.re * c3.re;
217  double bde = c1.im * c2.im * c3.re;
218  double adf = c1.re * c2.im * c3.im;
219  double bcf = c1.im * c2.re * c3.im;
220  double ade = c1.re * c2.im * c3.re;
221  double bce = c1.im * c2.re * c3.re;
222  double acf = c1.re * c2.re * c3.im;
223  double bdf = c1.im * c2.im * c3.im;
224 
225  return bse_complex (aec - bde - adf - bcf, ade + bce + acf - bdf);
226 }
227 static inline BseComplex
228 bse_complex_div (BseComplex a,
229  BseComplex b)
230 {
231  BseComplex c;
232  if (fabs (b.re) >= fabs (b.im))
233  {
234  double r = b.im / b.re, den = b.re + r * b.im;
235  c.re = (a.re + r * a.im) / den;
236  c.im = (a.im - r * a.re) / den;
237  }
238  else
239  {
240  double r = b.re / b.im, den = b.im + r * b.re;
241  c.re = (a.re * r + a.im) / den;
242  c.im = (a.im * r - a.re) / den;
243  }
244  return c;
245 }
246 static inline BseComplex
247 bse_complex_reciprocal (BseComplex c)
248 {
249  if (fabs (c.re) >= fabs (c.im))
250  {
251  double r = c.im / c.re, den = c.re + r * c.im;
252  c.re = 1. / den;
253  c.im = - r / den;
254  }
255  else
256  {
257  double r = c.re / c.im, den = c.im + r * c.re;
258  c.re = r / den;
259  c.im = - 1. / den;
260  }
261  return c;
262 }
263 static inline BseComplex
264 bse_complex_sqrt (BseComplex z)
265 {
266  if (z.re == 0.0 && z.im == 0.0)
267  return z;
268  else
269  {
270  BseComplex c;
271  double w, x = fabs (z.re), y = fabs (z.im);
272  if (x >= y)
273  {
274  double r = y / x;
275  w = sqrt (x) * sqrt (0.5 * (1.0 + sqrt (1.0 + r * r)));
276  }
277  else
278  {
279  double r = x / y;
280  w = sqrt (y) * sqrt (0.5 * (r + sqrt (1.0 + r * r)));
281  }
282  if (z.re >= 0.0)
283  {
284  c.re = w;
285  c.im = z.im / (2.0 * w);
286  }
287  else
288  {
289  c.im = z.im >= 0 ? w : -w;
290  c.re = z.im / (2.0 * c.im);
291  }
292  return c;
293  }
294 }
295 static inline BseComplex
296 bse_complex_conj (BseComplex c)
297 {
298  return bse_complex (c.re, -c.im);
299 }
300 static inline BseComplex
301 bse_complex_inv (BseComplex c)
302 {
303  return bse_complex (-c.re, -c.im);
304 }
305 static inline BseComplex
306 bse_complex_id (BseComplex c)
307 {
308  return c;
309 }
310 static inline double
311 bse_complex_abs (BseComplex c)
312 {
313  /* compute (a^2 + b^2)^(1/2) without destructive underflow or overflow */
314  double absa = fabs (c.re), absb = fabs (c.im);
315  return (absa > absb ?
316  absb == 0.0 ? absa :
317  absa * sqrt (1.0 + (absb / absa) * (absb / absa)) :
318  absb == 0.0 ? 0.0 :
319  absb * sqrt (1.0 + (absa / absb) * (absa / absb)));
320 }
321 static inline double
322 bse_complex_arg (BseComplex c)
323 {
324  double a = atan2 (c.im, c.re);
325  return a;
326 }
327 static inline BseComplex
328 bse_complex_sin (BseComplex c)
329 {
330  return bse_complex (sin (c.re) * cosh (c.im), cos (c.re) * sinh (c.im));
331 }
332 static inline BseComplex
333 bse_complex_cos (BseComplex c)
334 {
335  return bse_complex (cos (c.re) * cosh (c.im), - sin (c.re) * sinh (c.im));
336 }
337 static inline BseComplex
338 bse_complex_tan (BseComplex c)
339 {
340  return bse_complex_div (bse_complex (tan (c.re), tanh (c.im)),
341  bse_complex (1.0, -tan (c.re) * tanh (c.im)));
342 }
343 static inline BseComplex
344 bse_complex_sinh (BseComplex c)
345 {
346  return bse_complex (sinh (c.re) * cos (c.im), cosh (c.re) * sin (c.im));
347 }
348 static inline BseComplex
349 bse_complex_cosh (BseComplex c)
350 {
351  return bse_complex (cosh (c.re) * cos (c.im), sinh (c.re) * sin (c.im));
352 }
353 static inline BseComplex
354 bse_complex_tanh (BseComplex c)
355 {
356  return bse_complex_div (bse_complex_sinh (c),
357  bse_complex_cosh (c));
358 }
359 static inline void
360 bse_poly_add (uint degree,
361  double *a,
362  double *b)
363 {
364  uint i;
365 
366  for (i = 0; i <= degree; i++)
367  a[i] += b[i];
368 }
369 static inline void
370 bse_poly_sub (uint degree,
371  double *a,
372  double *b)
373 {
374  uint i;
375 
376  for (i = 0; i <= degree; i++)
377  a[i] -= b[i];
378 }
379 static inline void
380 bse_poly_mul (double *p, /* out:[0..aorder+border] */
381  uint aorder,
382  const double *a, /* in:[0..aorder] */
383  uint border,
384  const double *b) /* in:[0..border] */
385 {
386  uint i;
387 
388  for (i = aorder + border; i > 0; i--)
389  {
390  uint j;
391  double t = 0;
392 
393  for (j = i - MIN (border, i); j <= MIN (aorder, i); j++)
394  t += a[j] * b[i - j];
395  p[i] = t;
396  }
397  p[0] = a[0] * b[0];
398 }
399 static inline void
400 bse_cpoly_mul_monomial (uint degree,
401  BseComplex *c,
402  BseComplex root)
403 {
404  uint j;
405 
406  c[degree] = c[degree - 1];
407  for (j = degree - 1; j >= 1; j--)
408  c[j] = bse_complex_sub (c[j - 1], bse_complex_mul (c[j], root));
409  c[0] = bse_complex_mul (c[0], bse_complex_inv (root));
410 }
411 static inline void
412 bse_cpoly_mul_reciprocal (uint degree,
413  BseComplex *c,
414  BseComplex root)
415 {
416  uint j;
417 
418  c[degree] = bse_complex_mul (c[degree - 1], bse_complex_inv (root));
419  for (j = degree - 1; j >= 1; j--)
420  c[j] = bse_complex_sub (c[j], bse_complex_mul (c[j - 1], root));
421  /* c[0] = c[0]; */
422 }
423 static inline void
424 bse_cpoly_mul (BseComplex *p, /* [0..aorder+border] */
425  uint aorder,
426  BseComplex *a,
427  uint border,
428  BseComplex *b)
429 {
430  uint i;
431 
432  for (i = aorder + border; i > 0; i--)
433  {
434  BseComplex t;
435  uint j;
436 
437  t = bse_complex (0, 0);
438  for (j = i - MIN (i, border); j <= MIN (aorder, i); j++)
439  t = bse_complex_add (t, bse_complex_mul (a[j], b[i - j]));
440  p[i] = t;
441  }
442  p[0] = bse_complex_mul (a[0], b[0]);
443 }
444 static inline void
445 bse_poly_scale (uint degree,
446  double *a,
447  double scale)
448 {
449  uint i;
450 
451  for (i = 0; i <= degree; i++)
452  a[i] *= scale;
453 }
454 static inline void
455 bse_poly_xscale (uint degree,
456  double *a,
457  double xscale)
458 {
459  double scale = xscale;
460  uint i;
461 
462  for (i = 1; i <= degree; i++)
463  {
464  a[i] *= scale;
465  scale *= xscale;
466  }
467 }
468 static inline double
469 bse_poly_eval (uint degree,
470  double *a,
471  double x)
472 {
473  double sum = a[degree];
474 
475  while (degree--)
476  sum = sum * x + a[degree];
477  return sum;
478 }
479 
480 #endif /* __BSE_MATH_H__ */ /* vim: set ts=8 sw=2 sts=2: */
cosh
sqrt
Definition: bsemath.hh:37
abs
tanh
STL class.
sinh
fabs
atan2
sin
tan
cos