1 /*
2 SuperCollider real time audio synthesis system
3 Copyright (c) 2002 James McCartney. All rights reserved.
4 http://www.audiosynth.com
5
6 This program is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with this program; if not, write to the Free Software
18 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
19 */
20
21 #pragma once
22
23 #include "SC_Types.h"
24 #include "SC_Constants.h"
25
26 #include <cmath>
27 #include <limits>
28
29 #ifdef __SSE__
30 # include <xmmintrin.h>
31 #endif
32
33 #ifdef __SSE4_1__
34 # include <smmintrin.h>
35 #endif
36
37
38 ///////////////////////////////////////////////////////////////////////////////////////
39
sc_isnan(float x)40 inline bool sc_isnan(float x) { return std::isnan(x); }
41
sc_isnan(double x)42 inline bool sc_isnan(double x) { return std::isnan(x); }
43
sc_isfinite(float x)44 inline bool sc_isfinite(float x) { return std::isfinite(x); }
45
sc_isfinite(double x)46 inline bool sc_isfinite(double x) { return std::isfinite(x); }
47
48 ///////////////////////////////////////////////////////////////////////////////////////
49
50 // versions provided for float32 and float64
51 // did not supply template because do not want to instantiate for integers.
52 // all constants explicitly cast to prevent PowerPC frsp instruction generation.
53
54 ///////////////////////////////////////////////////////////////////////////////////////
55
56 // this is a function for preventing pathological math operations in ugens.
57 // can be used at the end of a block to fix any recirculating filter values.
zapgremlins(float32 x)58 inline float32 zapgremlins(float32 x) {
59 float32 absx = std::abs(x);
60 // very small numbers fail the first test, eliminating denormalized numbers
61 // (zero also fails the first test, but that is OK since it returns zero.)
62 // very large numbers fail the second test, eliminating infinities
63 // Not-a-Numbers fail both tests and are eliminated.
64 return (absx > (float32)1e-15 && absx < (float32)1e15) ? x : (float32)0.;
65 }
66
sc_log2(float32 x)67 inline float32 sc_log2(float32 x) { return std::log2(x); }
68
sc_log10(float32 x)69 inline float32 sc_log10(float32 x) { return std::log10(std::abs(x)); }
70
sc_midicps(float32 note)71 inline float32 sc_midicps(float32 note) {
72 return (float32)440. * std::pow((float32)2., (note - (float32)69.) * (float32)0.083333333333);
73 }
74
sc_cpsmidi(float32 freq)75 inline float32 sc_cpsmidi(float32 freq) {
76 return sc_log2(freq * (float32)0.0022727272727) * (float32)12. + (float32)69.;
77 }
78
sc_midiratio(float32 midi)79 inline float32 sc_midiratio(float32 midi) { return std::pow((float32)2., midi * (float32)0.083333333333); }
80
sc_ratiomidi(float32 ratio)81 inline float32 sc_ratiomidi(float32 ratio) { return (float32)12. * sc_log2(ratio); }
82
sc_octcps(float32 note)83 inline float32 sc_octcps(float32 note) { return (float32)440. * std::pow((float32)2., note - (float32)4.75); }
84
sc_cpsoct(float32 freq)85 inline float32 sc_cpsoct(float32 freq) { return sc_log2(freq * (float32)0.0022727272727) + (float32)4.75; }
86
sc_ampdb(float32 amp)87 inline float32 sc_ampdb(float32 amp) { return std::log10(amp) * (float32)20.; }
88
sc_dbamp(float32 db)89 inline float32 sc_dbamp(float32 db) { return std::pow((float32)10., db * (float32).05); }
90
sc_squared(float32 x)91 inline float32 sc_squared(float32 x) { return x * x; }
92
sc_cubed(float32 x)93 inline float32 sc_cubed(float32 x) { return x * x * x; }
94
sc_sqrt(float32 x)95 inline float32 sc_sqrt(float32 x) { return x < (float32)0. ? -sqrt(-x) : sqrt(x); }
96
97
sc_hanwindow(float32 x)98 inline float32 sc_hanwindow(float32 x) {
99 if (x < (float32)0. || x > (float32)1.)
100 return (float32)0.;
101 return (float32)0.5 - (float32)0.5 * static_cast<float32>(cos(x * (float32)twopi));
102 }
103
sc_welwindow(float32 x)104 inline float32 sc_welwindow(float32 x) {
105 if (x < (float32)0. || x > (float32)1.)
106 return (float32)0.;
107 return static_cast<float32>(sin(x * pi));
108 }
109
sc_triwindow(float32 x)110 inline float32 sc_triwindow(float32 x) {
111 if (x < (float32)0. || x > (float32)1.)
112 return (float32)0.;
113 if (x < (float32)0.5)
114 return (float32)2. * x;
115 else
116 return (float32)-2. * x + (float32)2.;
117 }
118
sc_bitriwindow(float32 x)119 inline float32 sc_bitriwindow(float32 x) {
120 float32 ax = (float32)1. - std::abs(x);
121 if (ax <= (float32)0.)
122 return (float32)0.;
123 return ax;
124 }
125
sc_rectwindow(float32 x)126 inline float32 sc_rectwindow(float32 x) {
127 if (x < (float32)0. || x > (float32)1.)
128 return (float32)0.;
129 return (float32)1.;
130 }
131
sc_scurve(float32 x)132 inline float32 sc_scurve(float32 x) {
133 if (x <= (float32)0.)
134 return (float32)0.;
135 if (x >= (float32)1.)
136 return (float32)1.;
137 return x * x * ((float32)3. - (float32)2. * x);
138 }
139
sc_scurve0(float32 x)140 inline float32 sc_scurve0(float32 x) {
141 // assumes that x is in range
142 return x * x * ((float32)3. - (float32)2. * x);
143 }
144
sc_ramp(float32 x)145 inline float32 sc_ramp(float32 x) {
146 if (x <= (float32)0.)
147 return (float32)0.;
148 if (x >= (float32)1.)
149 return (float32)1.;
150 return x;
151 }
152
sc_sign(float32 x)153 inline float32 sc_sign(float32 x) {
154 return x < (float32)0. ? (float32)-1. : (x > (float32)0. ? (float32)1.f : (float32)0.f);
155 }
156
sc_distort(float32 x)157 inline float32 sc_distort(float32 x) { return x / ((float32)1. + std::abs(x)); }
158
sc_distortneg(float32 x)159 inline float32 sc_distortneg(float32 x) {
160 if (x < (float32)0.)
161 return x / ((float32)1. - x);
162 else
163 return x;
164 }
165
sc_softclip(float32 x)166 inline float32 sc_softclip(float32 x) {
167 float32 absx = std::abs(x);
168 if (absx <= (float32)0.5)
169 return x;
170 else
171 return (absx - (float32)0.25) / x;
172 }
173
174 // Taylor expansion out to x**9/9! factored into multiply-adds
175 // from Phil Burk.
taylorsin(float32 x)176 inline float32 taylorsin(float32 x) {
177 // valid range from -pi/2 to +3pi/2
178 x = static_cast<float32>((float32)pi2 - std::abs(pi2 - x));
179 float32 x2 = x * x;
180 return static_cast<float32>(
181 x * (x2 * (x2 * (x2 * (x2 * (1.0 / 362880.0) - (1.0 / 5040.0)) + (1.0 / 120.0)) - (1.0 / 6.0)) + 1.0));
182 }
183
sc_trunc(float32 x)184 inline float32 sc_trunc(float32 x) { return std::trunc(x); }
185
186
sc_ceil(float32 x)187 inline float32 sc_ceil(float32 x) {
188 #ifdef __SSE4_1__
189 __m128 a = _mm_set_ss(x);
190 __m128 b = _mm_round_ss(a, a, _MM_FROUND_TO_POS_INF);
191 return _mm_cvtss_f32(b);
192 #else
193 return std::ceil(x);
194 #endif
195 }
196
sc_floor(float32 x)197 inline float32 sc_floor(float32 x) {
198 #ifdef __SSE4_1__
199 __m128 a = _mm_set_ss(x);
200 __m128 b = _mm_round_ss(a, a, _MM_FROUND_TO_NEG_INF);
201 return _mm_cvtss_f32(b);
202 #else
203 return std::floor(x);
204 #endif
205 }
206
sc_reciprocal(float32 x)207 inline float32 sc_reciprocal(float32 x) {
208 #ifdef __SSE__
209 // adapted from AP-803 Newton-Raphson Method with Streaming SIMD Extensions
210 // 23 bit accuracy (out of 24bit)
211 const __m128 arg = _mm_set_ss(x);
212 const __m128 approx = _mm_rcp_ss(arg);
213 const __m128 muls = _mm_mul_ss(_mm_mul_ss(arg, approx), approx);
214 const __m128 doubleApprox = _mm_add_ss(approx, approx);
215 const __m128 result = _mm_sub_ss(doubleApprox, muls);
216 return _mm_cvtss_f32(result);
217 #else
218 return 1.f / x;
219 #endif
220 }
221
sc_frac(float32 x)222 inline float32 sc_frac(float32 x) { return x - sc_floor(x); }
223
224
225 ////////////////////////////////
226
sc_bitNot(float32 x)227 inline float32 sc_bitNot(float32 x) { return (float32) ~(int)x; }
228
sc_lg3interp(float32 x1,float32 a,float32 b,float32 c,float32 d)229 inline float32 sc_lg3interp(float32 x1, float32 a, float32 b, float32 c, float32 d) {
230 // cubic lagrange interpolator
231 float32 x0 = x1 + 1.f;
232 float32 x2 = x1 - 1.f;
233 float32 x3 = x1 - 2.f;
234
235 float32 x03 = x0 * x3 * 0.5f;
236 float32 x12 = x1 * x2 * 0.16666666666666667f;
237
238 return x12 * (d * x0 - a * x3) + x03 * (b * x2 - c * x1);
239 }
240
sc_CalcFeedback(float32 delaytime,float32 decaytime)241 inline float32 sc_CalcFeedback(float32 delaytime, float32 decaytime) {
242 if (delaytime == 0.f || decaytime == 0.f)
243 return 0.f;
244
245 float32 absret = static_cast<float32>(std::exp(log001 * delaytime / std::abs(decaytime)));
246 float32 ret = std::copysign(absret, decaytime);
247 return ret;
248 }
249
sc_wrap1(float32 x)250 inline float32 sc_wrap1(float32 x) {
251 if (x >= (float32)1.)
252 return x + (float32)-2.;
253 if (x < (float32)-1.)
254 return x + (float32)2.;
255 return x;
256 }
257
sc_fold1(float32 x)258 inline float32 sc_fold1(float32 x) {
259 if (x >= (float32)1.)
260 return (float32)2. - x;
261 if (x < (float32)-1.)
262 return (float32)-2. - x;
263 return x;
264 }
265
266
267 ///////////////////////////////////////////////////////////////////////////////////////
268
zapgremlins(float64 x)269 inline float64 zapgremlins(float64 x) {
270 float64 absx = std::abs(x);
271 // very small numbers fail the first test, eliminating denormalized numbers
272 // (zero also fails the first test, but that is OK since it returns zero.)
273 // very large numbers fail the second test, eliminating infinities
274 // Not-a-Numbers fail both tests and are eliminated.
275 return (absx > (float64)1e-15 && absx < (float64)1e15) ? x : (float64)0.;
276 }
277
sc_log2(float64 x)278 inline float64 sc_log2(float64 x) { return std::log2(std::abs(x)); }
279
sc_log10(float64 x)280 inline float64 sc_log10(float64 x) { return std::log10(std::abs(x)); }
281
sc_midicps(float64 note)282 inline float64 sc_midicps(float64 note) {
283 return (float64)440. * std::pow((float64)2., (note - (float64)69.) * (float64)0.08333333333333333333333333);
284 }
285
sc_cpsmidi(float64 freq)286 inline float64 sc_cpsmidi(float64 freq) {
287 return sc_log2(freq * (float64)0.002272727272727272727272727) * (float64)12. + (float64)69.;
288 }
289
sc_midiratio(float64 midi)290 inline float64 sc_midiratio(float64 midi) { return std::pow((float64)2., midi * (float64)0.083333333333); }
291
sc_ratiomidi(float64 ratio)292 inline float64 sc_ratiomidi(float64 ratio) { return (float64)12. * sc_log2(ratio); }
293
sc_octcps(float64 note)294 inline float64 sc_octcps(float64 note) { return (float64)440. * std::pow((float64)2., note - (float64)4.75); }
295
sc_cpsoct(float64 freq)296 inline float64 sc_cpsoct(float64 freq) { return sc_log2(freq * (float64)0.0022727272727) + (float64)4.75; }
297
sc_ampdb(float64 amp)298 inline float64 sc_ampdb(float64 amp) { return std::log10(amp) * (float64)20.; }
299
sc_dbamp(float64 db)300 inline float64 sc_dbamp(float64 db) { return std::pow((float64)10., db * (float64).05); }
301
sc_squared(float64 x)302 inline float64 sc_squared(float64 x) { return x * x; }
303
sc_cubed(float64 x)304 inline float64 sc_cubed(float64 x) { return x * x * x; }
305
sc_sqrt(float64 x)306 inline float64 sc_sqrt(float64 x) { return x < (float64)0. ? -sqrt(-x) : sqrt(x); }
307
sc_hanwindow(float64 x)308 inline float64 sc_hanwindow(float64 x) {
309 if (x < (float64)0. || x > (float64)1.)
310 return (float64)0.;
311 return (float64)0.5 - (float64)0.5 * cos(x * twopi);
312 }
313
sc_welwindow(float64 x)314 inline float64 sc_welwindow(float64 x) {
315 if (x < (float64)0. || x > (float64)1.)
316 return (float64)0.;
317 return sin(x * pi);
318 }
319
sc_triwindow(float64 x)320 inline float64 sc_triwindow(float64 x) {
321 if (x < (float64)0. || x > (float64)1.)
322 return (float64)0.;
323 if (x < (float64)0.5)
324 return (float64)2. * x;
325 else
326 return (float64)-2. * x + (float64)2.;
327 }
328
sc_bitriwindow(float64 x)329 inline float64 sc_bitriwindow(float64 x) {
330 float64 ax = std::abs(x);
331 if (ax > (float64)1.)
332 return (float64)0.;
333 return (float64)1. - ax;
334 }
335
sc_rectwindow(float64 x)336 inline float64 sc_rectwindow(float64 x) {
337 if (x < (float64)0. || x > (float64)1.)
338 return (float64)0.;
339 return (float64)1.;
340 }
341
sc_scurve(float64 x)342 inline float64 sc_scurve(float64 x) {
343 if (x <= (float64)0.)
344 return (float64)0.;
345 if (x >= (float64)1.)
346 return (float64)1.;
347 return x * x * ((float64)3. - (float64)2. * x);
348 }
349
sc_scurve0(float64 x)350 inline float64 sc_scurve0(float64 x) {
351 // assumes that x is in range
352 return x * x * ((float64)3. - (float64)2. * x);
353 }
354
sc_ramp(float64 x)355 inline float64 sc_ramp(float64 x) {
356 if (x <= (float64)0.)
357 return (float64)0.;
358 if (x >= (float64)1.)
359 return (float64)1.;
360 return x;
361 }
362
sc_sign(float64 x)363 inline float64 sc_sign(float64 x) {
364 return x < (float64)0. ? (float64)-1. : (x > (float64)0. ? (float64)1.f : (float64)0.f);
365 }
366
sc_distort(float64 x)367 inline float64 sc_distort(float64 x) { return x / ((float64)1. + std::abs(x)); }
368
sc_distortneg(float64 x)369 inline float64 sc_distortneg(float64 x) {
370 if (x < (float64)0.)
371 return x / ((float64)1. - x);
372 else
373 return x;
374 }
375
376
sc_softclip(float64 x)377 inline float64 sc_softclip(float64 x) {
378 float64 absx = std::abs(x);
379 if (absx <= (float64)0.5)
380 return x;
381 else
382 return (absx - (float64)0.25) / x;
383 }
384
385 // Taylor expansion out to x**9/9! factored into multiply-adds
386 // from Phil Burk.
taylorsin(float64 x)387 inline float64 taylorsin(float64 x) {
388 x = pi2 - std::abs(pi2 - x);
389 float64 x2 = x * x;
390 return x * (x2 * (x2 * (x2 * (x2 * (1.0 / 362880.0) - (1.0 / 5040.0)) + (1.0 / 120.0)) - (1.0 / 6.0)) + 1.0);
391 }
392
sc_trunc(float64 x)393 inline float64 sc_trunc(float64 x) { return std::trunc(x); }
394
sc_ceil(float64 x)395 inline float64 sc_ceil(float64 x) {
396 #ifdef __SSE4_1__
397 __m128d a = _mm_set_sd(x);
398 const int cntrl = _MM_FROUND_TO_POS_INF;
399 __m128d b = _mm_round_sd(a, a, cntrl);
400 return _mm_cvtsd_f64(b);
401 #else
402 return std::ceil(x);
403 #endif
404 }
405
sc_floor(float64 x)406 inline float64 sc_floor(float64 x) {
407 #ifdef __SSE4_1__
408 __m128d a = _mm_set_sd(x);
409 const int cntrl = _MM_FROUND_TO_NEG_INF;
410 __m128d b = _mm_round_sd(a, a, cntrl);
411 return _mm_cvtsd_f64(b);
412 #else
413 return std::floor(x);
414 #endif
415 }
416
sc_reciprocal(float64 x)417 inline float64 sc_reciprocal(float64 x) { return 1. / x; }
418
419
sc_frac(float64 x)420 inline float64 sc_frac(float64 x) { return x - sc_floor(x); }
421
sc_wrap1(float64 x)422 inline float64 sc_wrap1(float64 x) {
423 if (x >= (float64)1.)
424 return x + (float64)-2.;
425 if (x < (float64)-1.)
426 return x + (float64)2.;
427 return x;
428 }
429
sc_fold1(float64 x)430 inline float64 sc_fold1(float64 x) {
431 if (x >= (float64)1.)
432 return (float64)2. - x;
433 if (x < (float64)-1.)
434 return (float64)-2. - x;
435 return x;
436 }
437
sc_grayCode(int32 x)438 inline int32 sc_grayCode(int32 x) { return x ^ (x >> 1); }
439