1 /*
2 cmath.c:
3
4 Copyright (C) 1994 Paris Smaragdis, John ffitch
5
6 This file is part of Csound.
7
8 The Csound Library is free software; you can redistribute it
9 and/or modify it under the terms of the GNU Lesser General Public
10 License as published by the Free Software Foundation; either
11 version 2.1 of the License, or (at your option) any later version.
12
13 Csound is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU Lesser General Public License for more details.
17
18 You should have received a copy of the GNU Lesser General Public
19 License along with Csound; if not, write to the Free Software
20 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
21 02110-1301 USA
22 */
23
24 /* Math functions for Csound coded by Paris Smaragdis 1994 */
25 /* Berklee College of Music Csound development team */
26
27 #include "csoundCore.h"
28 #include "cmath.h"
29 #include <math.h>
30
ipow(CSOUND * csound,POW * p)31 int32_t ipow(CSOUND *csound, POW *p) /* Power for i-rate */
32 {
33 MYFLT in = *p->in;
34 MYFLT powerOf = *p->powerOf;
35 if (UNLIKELY(in == FL(0.0) && powerOf == FL(0.0)))
36 return csound->PerfError(csound, &(p->h), Str("NaN in pow\n"));
37 else if (p->norm!=NULL && *p->norm != FL(0.0))
38 *p->sr = POWER(in, powerOf) / *p->norm;
39 else
40 *p->sr = POWER(in, powerOf);
41 return OK;
42 }
43
apow(CSOUND * csound,POW * p)44 int32_t apow(CSOUND *csound, POW *p) /* Power routine for a-rate */
45 {
46 uint32_t offset = p->h.insdshead->ksmps_offset;
47 uint32_t early = p->h.insdshead->ksmps_no_end;
48 uint32_t n, nsmps = CS_KSMPS;
49 MYFLT *in = p->in, *out = p->sr;
50 MYFLT powerOf = *p->powerOf;
51 MYFLT norm = (p->norm!=NULL ? *p->norm : FL(1.0));
52 if (norm==FL(0.0)) norm = FL(1.0);
53 memset(out, '\0', offset*sizeof(MYFLT));
54 if (UNLIKELY(early)) {
55 nsmps -= early;
56 memset(&out[nsmps], '\0', early*sizeof(MYFLT));
57 }
58 if (UNLIKELY(powerOf == FL(0.0))) {
59 MYFLT yy = FL(1.0) / norm;
60 for (n = offset; n < nsmps; n++) {
61 MYFLT xx = in[n];
62 if (UNLIKELY(xx == FL(0.0))) {
63 return csound->PerfError(csound, &(p->h),Str("NaN in pow\n"));
64 }
65 else
66 out[n] = yy;
67 }
68 }
69 else {
70 for (n = offset; n < nsmps; n++)
71 out[n] = POWER(in[n], powerOf) / norm;
72 }
73 return OK;
74 }
75
seedrand(CSOUND * csound,PRAND * p)76 int32_t seedrand(CSOUND *csound, PRAND *p)
77 {
78 uint32_t seedVal = (uint32_t)0;
79 int32_t xx = (int32_t)((double)*p->out + 0.5);
80
81 if (xx > FL(0.0))
82 seedVal = (uint32_t)xx;
83 else if (xx==0) {
84 seedVal = (uint32_t)csound->GetRandomSeedFromTime();
85 csound->Warning(csound, Str("Seeding from current time %u\n"),
86 (uint32_t)seedVal);
87 }
88 else
89 csound->Warning(csound, Str("Seeding with %u\n"), (uint32_t)seedVal);
90 csound->SeedRandMT(&(csound->randState_), NULL, seedVal);
91 csound->holdrand = (int32_t)(seedVal & (uint32_t) 0x7FFFFFFF);
92 while (seedVal >= (uint32_t)0x7FFFFFFE)
93 seedVal -= (uint32_t)0x7FFFFFFE;
94 csound->randSeed1 = ((int32_t)seedVal + 1);
95
96 return OK;
97 }
98
getseed(CSOUND * csound,GETSEED * p)99 int32_t getseed(CSOUND *csound, GETSEED *p)
100 {
101 *p->ans = csound->randSeed1;
102 return OK;
103 }
104
105 /* * * * * * RANDOM NUMBER GENERATORS * * * * * */
106
107 #define UInt32toFlt(x) ((double)(x) * (1.0 / 4294967295.03125))
108
109 #define unirand(c) ((MYFLT) UInt32toFlt(csoundRandMT(&((c)->randState_))))
110
unifrand(CSOUND * csound,MYFLT range)111 static inline MYFLT unifrand(CSOUND *csound, MYFLT range)
112 {
113 return (range * unirand(csound));
114 }
115
116 /* linear distribution routine */
117
linrand(CSOUND * csound,MYFLT range)118 static inline MYFLT linrand(CSOUND *csound, MYFLT range)
119 {
120 uint32_t r1, r2;
121
122 r1 = csoundRandMT(&(csound->randState_));
123 r2 = csoundRandMT(&(csound->randState_));
124
125 return ((MYFLT)UInt32toFlt(r1 < r2 ? r1 : r2) * range);
126 }
127
128 /* triangle distribution routine */
129
trirand(CSOUND * csound,MYFLT range)130 static inline MYFLT trirand(CSOUND *csound, MYFLT range)
131 {
132 uint64_t r1;
133
134 r1 = (uint64_t)csoundRandMT(&(csound->randState_));
135 r1 += (uint64_t)csoundRandMT(&(csound->randState_));
136
137 return ((MYFLT) ((double)((int64_t)r1 - (int64_t)0xFFFFFFFFU)
138 * (1.0 / 4294967295.03125)) * range);
139 }
140
141 /* exponential distribution routine */
142
exprand(CSOUND * csound,MYFLT lambda)143 static MYFLT exprand(CSOUND *csound, MYFLT lambda)
144 {
145 uint32_t r1;
146
147 if (UNLIKELY(lambda < FL(0.0))) return (FL(0.0)); /* for safety */
148
149 do {
150 r1 = csoundRandMT(&(csound->randState_));
151 } while (!r1);
152
153 return -((MYFLT)log(UInt32toFlt(r1)) * lambda);
154 }
155
156 /* bilateral exponential distribution routine */
157
biexprand(CSOUND * csound,MYFLT range)158 static MYFLT biexprand(CSOUND *csound, MYFLT range)
159 {
160 int32_t r1;
161
162 if (UNLIKELY(range < FL(0.0))) return (FL(0.0)); /* For safety */
163
164 while ((r1 = (int32_t)csoundRandMT(&(csound->randState_)))==0);
165
166 if (r1 < (int32_t)0) {
167 return -(LOG(-(r1) * (FL(1.0) / FL(2147483648.0))) * range);
168 }
169 return (LOG(r1 * (FL(1.0) / FL(2147483648.0))) * range);
170 }
171
172 /* gaussian distribution routine */
173
gaussrand(CSOUND * csound,MYFLT s)174 static MYFLT gaussrand(CSOUND *csound, MYFLT s)
175 {
176 int64_t r1 = -((int64_t)0xFFFFFFFFU * 6);
177 int32_t n = 12;
178 double x;
179
180 do {
181 r1 += (int64_t)csoundRandMT(&(csound->randState_));
182 } while (--n);
183 x = (double)r1;
184 return (MYFLT)(x * ((double)s * (1.0 / (3.83 * 4294967295.03125))));
185 }
186
187 /* cauchy distribution routine */
188
cauchrand(CSOUND * csound,MYFLT a)189 static MYFLT cauchrand(CSOUND *csound, MYFLT a)
190 {
191 uint32_t r1;
192 MYFLT x;
193
194 do {
195 r1 = csoundRandMT(&(csound->randState_)); /* Limit range artificially */
196 } while (r1 > (uint32_t)2143188560U && r1 < (uint32_t)2151778735U);
197 x = TAN((MYFLT)r1 * (PI_F / FL(4294967295.0))) * (FL(1.0) / FL(318.3));
198 return (x * a);
199 }
200
201 /* positive cauchy distribution routine */
202
pcauchrand(CSOUND * csound,MYFLT a)203 static MYFLT pcauchrand(CSOUND *csound, MYFLT a)
204 {
205 uint32_t r1;
206 MYFLT x;
207
208 do {
209 r1 = csoundRandMT(&(csound->randState_));
210 } while (r1 > (uint32_t)4286377121U); /* Limit range artificially */
211 x = TAN((MYFLT)r1 * HALFPI_F / FL(4294967295.0))
212 * (FL(1.0) / FL(318.3));
213 return (x * a);
214 }
215
216 /* beta distribution routine */
217
betarand(CSOUND * csound,MYFLT range,MYFLT a,MYFLT b)218 static MYFLT betarand(CSOUND *csound, MYFLT range, MYFLT a, MYFLT b)
219 {
220 double r1, r2;
221 double aa, bb;
222 if (UNLIKELY(a <= FL(0.0) || b <= FL(0.0)))
223 return FL(0.0);
224
225 aa = (double)a; bb = (double)b;
226 do {
227 uint32_t tmp;
228 do {
229 tmp = csoundRandMT(&(csound->randState_));
230 } while (!tmp);
231 r1 = pow(UInt32toFlt(tmp), 1.0 / aa);
232 do {
233 tmp = csoundRandMT(&(csound->randState_));
234 } while (!tmp);
235 r2 = r1 + pow(UInt32toFlt(tmp), 1.0 / bb);
236 } while (r2 > 1.0);
237
238 return (((MYFLT)r1 / (MYFLT)r2) * range);
239 }
240
241 /* weibull distribution routine */
242
weibrand(CSOUND * csound,MYFLT s,MYFLT t)243 static MYFLT weibrand(CSOUND *csound, MYFLT s, MYFLT t)
244 {
245 uint32_t r1;
246 double r2;
247
248 if (UNLIKELY(t <= FL(0.0))) return FL(0.0);
249
250 do {
251 r1 = csoundRandMT(&(csound->randState_));
252 } while (!r1 || r1 == (uint32_t)0xFFFFFFFFU);
253
254 r2 = 1.0 - ((double)r1 * (1.0 / 4294967295.0));
255
256 return (s * (MYFLT)pow(-(log(r2)), (1.0 / (double)t)));
257 }
258
259 /* Poisson distribution routine */
260
poissrand(CSOUND * csound,MYFLT lambda)261 static MYFLT poissrand(CSOUND *csound, MYFLT lambda)
262 {
263 MYFLT r1, r2, r3;
264
265 if (UNLIKELY(lambda < FL(0.0))) return FL(0.0);
266
267 r1 = unirand(csound);
268 r2 = EXP(-lambda);
269 r3 = FL(0.0);
270
271 while (r1 >= r2) {
272 r3++;
273 r1 *= unirand(csound);
274 }
275
276 return (r3);
277 }
278
279 /* ------------------------------------------------------------------------ */
280
auniform(CSOUND * csound,PRAND * p)281 int32_t auniform(CSOUND *csound, PRAND *p) /* Uniform distribution */
282 {
283 MYFLT *out = p->out;
284 uint32_t offset = p->h.insdshead->ksmps_offset;
285 uint32_t early = p->h.insdshead->ksmps_no_end;
286 uint32_t n, nsmps = CS_KSMPS;
287 double scale = (double)*p->arg1 * (1.0 / 4294967295.03125);
288
289 if (UNLIKELY(offset)) memset(out, '\0', offset*sizeof(MYFLT));
290 if (UNLIKELY(early)) {
291 nsmps -= early;
292 memset(&out[nsmps], '\0', early*sizeof(MYFLT));
293 }
294 for (n=offset; n<nsmps; n++) {
295 out[n] = (MYFLT)((double)csoundRandMT(&(csound->randState_)) * scale);
296 }
297 return OK;
298 }
299
ikuniform(CSOUND * csound,PRAND * p)300 int32_t ikuniform(CSOUND *csound, PRAND *p)
301 {
302 *p->out = unifrand(csound, *p->arg1);
303 return OK;
304 }
305
alinear(CSOUND * csound,PRAND * p)306 int32_t alinear(CSOUND *csound, PRAND *p) /* Linear random functions */
307 {
308 uint32_t offset = p->h.insdshead->ksmps_offset;
309 uint32_t early = p->h.insdshead->ksmps_no_end;
310 uint32_t n, nsmps = CS_KSMPS;
311 MYFLT *out = p->out;
312 MYFLT arg1 = *p->arg1;
313
314 memset(out, '\0', offset*sizeof(MYFLT));
315 if (UNLIKELY(early)) {
316 nsmps -= early;
317 memset(&out[nsmps], '\0', early*sizeof(MYFLT));
318 }
319 for (n = offset; n < nsmps; n++)
320 out[n] = linrand(csound, arg1);
321 return OK;
322
323 }
324
iklinear(CSOUND * csound,PRAND * p)325 int32_t iklinear(CSOUND *csound, PRAND *p)
326 {
327 *p->out = linrand(csound, *p->arg1);
328 return OK;
329 }
330
atrian(CSOUND * csound,PRAND * p)331 int32_t atrian(CSOUND *csound, PRAND *p) /* Triangle random functions */
332 {
333 uint32_t offset = p->h.insdshead->ksmps_offset;
334 uint32_t early = p->h.insdshead->ksmps_no_end;
335 uint32_t n, nsmps = CS_KSMPS;
336 MYFLT *out = p->out;
337 MYFLT arg1 = *p->arg1;
338
339 if (UNLIKELY(offset)) memset(out, '\0', offset*sizeof(MYFLT));
340 if (UNLIKELY(early)) {
341 nsmps -= early;
342 memset(&out[nsmps], '\0', early*sizeof(MYFLT));
343 }
344 for (n = offset; n < nsmps; n++)
345 out[n] = trirand(csound, arg1);
346 return OK;
347 }
348
iktrian(CSOUND * csound,PRAND * p)349 int32_t iktrian(CSOUND *csound, PRAND *p)
350 {
351 *p->out = trirand(csound, *p->arg1);
352 return OK;
353 }
354
exprndiset(CSOUND * csound,PRANDI * p)355 int32_t exprndiset(CSOUND *csound, PRANDI *p)
356 {
357 p->num1 = exprand(csound, *p->arg1);
358 p->num2 = exprand(csound, *p->arg1);
359 p->dfdmax = (p->num2 - p->num1) / FMAXLEN;
360 p->phs = 0;
361 p->ampcod = IS_ASIG_ARG(p->xamp) ? 1 : 0; /* (not used by krandi) */
362 p->cpscod = IS_ASIG_ARG(p->xcps) ? 1 : 0;
363 return OK;
364 }
365
kexprndi(CSOUND * csound,PRANDI * p)366 int kexprndi(CSOUND *csound, PRANDI *p)
367 { /* rslt = (num1 + diff*phs) * amp */
368 /* IV - Jul 11 2002 */
369 *p->ar = (p->num1 + (MYFLT)p->phs * p->dfdmax) * *p->xamp;
370 p->phs += (int32_t)(*p->xcps * CS_KICVT); /* phs += inc */
371 if (UNLIKELY(p->phs >= MAXLEN)) { /* when phs overflows, */
372 p->phs &= PHMASK; /* mod the phs */
373 p->num1 = p->num2; /* & new num vals */
374 p->num2 = exprand(csound, *p->arg1);
375 p->dfdmax = (p->num2 - p->num1) / FMAXLEN;
376 }
377 return OK;
378 }
379
iexprndi(CSOUND * csound,PRANDI * p)380 int32_t iexprndi(CSOUND *csound, PRANDI *p)
381 {
382 exprndiset(csound, p);
383 return kexprndi(csound, p);
384 }
385
aexprndi(CSOUND * csound,PRANDI * p)386 int32_t aexprndi(CSOUND *csound, PRANDI *p)
387 {
388 int32_t phs = p->phs, inc;
389 uint32_t offset = p->h.insdshead->ksmps_offset;
390 uint32_t early = p->h.insdshead->ksmps_no_end;
391 uint32_t n, nsmps = CS_KSMPS;
392 MYFLT *ar, *ampp, *cpsp;
393
394 cpsp = p->xcps;
395 ampp = p->xamp;
396 ar = p->ar;
397 inc = (int32_t)(cpsp[0] * csound->sicvt);
398 if (UNLIKELY(offset)) memset(ar, '\0', offset*sizeof(MYFLT));
399 if (UNLIKELY(early)) {
400 nsmps -= early;
401 memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
402 }
403 for (n=offset;n<nsmps;n++) {
404 /* IV - Jul 11 2002 */
405 if (p->ampcod)
406 ar[n] = (p->num1 + (MYFLT)phs * p->dfdmax) * ampp[n];
407 else
408 ar[n] = (p->num1 + (MYFLT)phs * p->dfdmax) * ampp[0];
409 phs += inc; /* phs += inc */
410 if (p->cpscod)
411 inc = (int32_t)(cpsp[n] * csound->sicvt); /* (nxt inc) */
412 if (UNLIKELY(phs >= MAXLEN)) { /* when phs o'flows */
413 phs &= PHMASK;
414 p->num1 = p->num2;
415 p->num2 = exprand(csound, *p->arg1);
416 p->dfdmax = (p->num2 - p->num1) / FMAXLEN;
417 }
418 }
419 p->phs = phs;
420 return OK;
421 }
422
aexp(CSOUND * csound,PRAND * p)423 int32_t aexp(CSOUND *csound, PRAND *p) /* Exponential random functions */
424 {
425 uint32_t offset = p->h.insdshead->ksmps_offset;
426 uint32_t early = p->h.insdshead->ksmps_no_end;
427 uint32_t n, nsmps = CS_KSMPS;
428 MYFLT *out = p->out;
429 MYFLT arg1 = *p->arg1;
430
431 if (UNLIKELY(offset)) memset(out, '\0', offset*sizeof(MYFLT));
432 if (UNLIKELY(early)) {
433 nsmps -= early;
434 memset(&out[nsmps], '\0', early*sizeof(MYFLT));
435 }
436 for (n = offset; n < nsmps; n++)
437 out[n] = exprand(csound, arg1);
438 return OK;
439 }
440
ikexp(CSOUND * csound,PRAND * p)441 int32_t ikexp(CSOUND *csound, PRAND *p)
442 {
443 *p->out = exprand(csound, *p->arg1);
444 return OK;
445 }
446
abiexp(CSOUND * csound,PRAND * p)447 int32_t abiexp(CSOUND *csound, PRAND *p) /* Bilateral exponential rand */
448 { /* functions */
449 uint32_t offset = p->h.insdshead->ksmps_offset;
450 uint32_t early = p->h.insdshead->ksmps_no_end;
451 uint32_t n, nsmps = CS_KSMPS;
452 MYFLT *out = p->out;
453 MYFLT arg1 = *p->arg1;
454
455 if (UNLIKELY(offset)) memset(out, '\0', offset*sizeof(MYFLT));
456 if (UNLIKELY(early)) {
457 nsmps -= early;
458 memset(&out[nsmps], '\0', early*sizeof(MYFLT));
459 }
460 for (n = offset; n < nsmps; n++)
461 out[n] = biexprand(csound, arg1);
462 return OK;
463 }
464
ikbiexp(CSOUND * csound,PRAND * p)465 int32_t ikbiexp(CSOUND *csound, PRAND *p)
466 {
467 *p->out = biexprand(csound, *p->arg1);
468 return OK;
469 }
470
agaus(CSOUND * csound,PRAND * p)471 int32_t agaus(CSOUND *csound, PRAND *p) /* Gaussian random functions */
472 {
473 uint32_t offset = p->h.insdshead->ksmps_offset;
474 uint32_t early = p->h.insdshead->ksmps_no_end;
475 uint32_t n, nsmps = CS_KSMPS;
476 MYFLT *out = p->out;
477 MYFLT arg1 = *p->arg1;
478
479 if (UNLIKELY(offset)) memset(out, '\0', offset*sizeof(MYFLT));
480 if (UNLIKELY(early)) {
481 nsmps -= early;
482 memset(&out[nsmps], '\0', early*sizeof(MYFLT));
483 }
484 for (n = offset; n < nsmps; n++)
485 out[n] = gaussrand(csound, arg1);
486 return OK;
487 }
488
ikgaus(CSOUND * csound,PRAND * p)489 int32_t ikgaus(CSOUND *csound, PRAND *p)
490 {
491 *p->out = gaussrand(csound, *p->arg1);
492 return OK;
493 }
494
acauchy(CSOUND * csound,PRAND * p)495 int32_t acauchy(CSOUND *csound, PRAND *p) /* Cauchy random functions */
496 {
497 uint32_t offset = p->h.insdshead->ksmps_offset;
498 uint32_t early = p->h.insdshead->ksmps_no_end;
499 uint32_t n, nsmps = CS_KSMPS;
500 MYFLT *out = p->out;
501 MYFLT arg1 = *p->arg1;
502
503 if (UNLIKELY(offset)) memset(out, '\0', offset*sizeof(MYFLT));
504 if (UNLIKELY(early)) {
505 nsmps -= early;
506 memset(&out[nsmps], '\0', early*sizeof(MYFLT));
507 }
508 for (n = offset; n < nsmps; n++)
509 out[n] = cauchrand(csound, arg1);
510 return OK;
511 }
512
gaussiset(CSOUND * csound,PRANDI * p)513 int32_t gaussiset(CSOUND *csound, PRANDI *p)
514 {
515 p->num1 = gaussrand(csound, *p->arg1);
516 p->num2 = gaussrand(csound, *p->arg1);
517 p->dfdmax = (p->num2 - p->num1) / FMAXLEN;
518 p->phs = 0;
519 p->ampcod = IS_ASIG_ARG(p->xamp) ? 1 : 0; /* (not used by krandi) */
520 p->cpscod = IS_ASIG_ARG(p->xcps) ? 1 : 0;
521 return OK;
522 }
523
kgaussi(CSOUND * csound,PRANDI * p)524 int32_t kgaussi(CSOUND *csound, PRANDI *p)
525 { /* rslt = (num1 + diff*phs) * amp */
526 /* IV - Jul 11 2002 */
527 *p->ar = (p->num1 + (MYFLT)p->phs * p->dfdmax) * *p->xamp;
528 p->phs += (int32_t)(*p->xcps * CS_KICVT); /* phs += inc */
529 if (UNLIKELY(p->phs >= MAXLEN)) { /* when phs overflows, */
530 p->phs &= PHMASK; /* mod the phs */
531 p->num1 = p->num2; /* & new num vals */
532 p->num2 = gaussrand(csound, *p->arg1);
533 p->dfdmax = (p->num2 - p->num1) / FMAXLEN;
534 }
535 return OK;
536 }
537
igaussi(CSOUND * csound,PRANDI * p)538 int32_t igaussi(CSOUND *csound, PRANDI *p)
539 {
540 gaussiset(csound, p);
541 return kgaussi(csound, p);
542 }
543
agaussi(CSOUND * csound,PRANDI * p)544 int32_t agaussi(CSOUND *csound, PRANDI *p)
545 {
546 int32_t phs = p->phs, inc;
547 uint32_t offset = p->h.insdshead->ksmps_offset;
548 uint32_t early = p->h.insdshead->ksmps_no_end;
549 uint32_t n, nsmps = CS_KSMPS;
550 MYFLT *ar, *ampp, *cpsp;
551
552 cpsp = p->xcps;
553 ampp = p->xamp;
554 ar = p->ar;
555 inc = (int32_t)(*cpsp * csound->sicvt);
556 if (UNLIKELY(offset)) memset(ar, '\0', offset*sizeof(MYFLT));
557 if (UNLIKELY(early)) {
558 nsmps -= early;
559 memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
560 }
561 for (n=offset;n<nsmps;n++) {
562 /* IV - Jul 11 2002 */
563 if (p->ampcod)
564 ar[n] = (p->num1 + (MYFLT)phs * p->dfdmax) * ampp[n];
565 else
566 ar[n] = (p->num1 + (MYFLT)phs * p->dfdmax) * ampp[0];
567 phs += inc; /* phs += inc */
568 if (p->cpscod)
569 inc = (int32_t)(cpsp[n] * csound->sicvt); /* (nxt inc) */
570 if (UNLIKELY(phs >= MAXLEN)) { /* when phs o'flows */
571 phs &= PHMASK;
572 p->num1 = p->num2;
573 p->num2 = gaussrand(csound, *p->arg1);
574 p->dfdmax = (p->num2 - p->num1) / FMAXLEN;
575 }
576 }
577 p->phs = phs;
578 return OK;
579 }
580
ikcauchy(CSOUND * csound,PRAND * p)581 int32_t ikcauchy(CSOUND *csound, PRAND *p)
582 {
583 *p->out = cauchrand(csound, *p->arg1);
584 return OK;
585 }
586
apcauchy(CSOUND * csound,PRAND * p)587 int32_t apcauchy(CSOUND *csound, PRAND *p) /* +ve Cauchy random functions */
588 {
589 uint32_t offset = p->h.insdshead->ksmps_offset;
590 uint32_t early = p->h.insdshead->ksmps_no_end;
591 uint32_t n, nsmps = CS_KSMPS;
592 MYFLT *out = p->out;
593 MYFLT arg1 = *p->arg1;
594
595 if (UNLIKELY(offset)) memset(out, '\0', offset*sizeof(MYFLT));
596 if (UNLIKELY(early)) {
597 nsmps -= early;
598 memset(&out[nsmps], '\0', early*sizeof(MYFLT));
599 }
600 for (n = offset; n < nsmps; n++)
601 out[n] = pcauchrand(csound, arg1);
602 return OK;
603 }
604
ikpcauchy(CSOUND * csound,PRAND * p)605 int32_t ikpcauchy(CSOUND *csound, PRAND *p)
606 {
607 *p->out = pcauchrand(csound, *p->arg1);
608 return OK;
609 }
610
cauchyiset(CSOUND * csound,PRANDI * p)611 int32_t cauchyiset(CSOUND *csound, PRANDI *p)
612 {
613 p->num1 = cauchrand(csound, *p->arg1);
614 p->num2 = cauchrand(csound, *p->arg1);
615 p->dfdmax = (p->num2 - p->num1) / FMAXLEN;
616 p->phs = 0;
617 p->ampcod = IS_ASIG_ARG(p->xamp) ? 1 : 0; /* (not used by krandi) */
618 p->cpscod = IS_ASIG_ARG(p->xcps) ? 1 : 0;
619 return OK;
620 }
621
kcauchyi(CSOUND * csound,PRANDI * p)622 int32_t kcauchyi(CSOUND *csound, PRANDI *p)
623 { /* rslt = (num1 + diff*phs) * amp */
624 /* IV - Jul 11 2002 */
625 *p->ar = (p->num1 + (MYFLT)p->phs * p->dfdmax) * *p->xamp;
626 p->phs += (int32_t)(*p->xcps * CS_KICVT); /* phs += inc */
627 if (UNLIKELY(p->phs >= MAXLEN)) { /* when phs overflows, */
628 p->phs &= PHMASK; /* mod the phs */
629 p->num1 = p->num2; /* & new num vals */
630 p->num2 = cauchrand(csound, *p->arg1);
631 p->dfdmax = (p->num2 - p->num1) / FMAXLEN;
632 }
633 return OK;
634 }
635
icauchyi(CSOUND * csound,PRANDI * p)636 int32_t icauchyi(CSOUND *csound, PRANDI *p)
637 {
638 cauchyiset(csound, p);
639 return kcauchyi(csound, p);
640 }
641
acauchyi(CSOUND * csound,PRANDI * p)642 int32_t acauchyi(CSOUND *csound, PRANDI *p)
643 {
644 int32_t phs = p->phs, inc;
645 uint32_t offset = p->h.insdshead->ksmps_offset;
646 uint32_t early = p->h.insdshead->ksmps_no_end;
647 uint32_t n, nsmps = CS_KSMPS;
648 MYFLT *ar, *ampp, *cpsp;
649
650 cpsp = p->xcps;
651 ampp = p->xamp;
652 ar = p->ar;
653 inc = (int32_t)(*cpsp * csound->sicvt);
654 if (UNLIKELY(offset)) memset(ar, '\0', offset*sizeof(MYFLT));
655 if (UNLIKELY(early)) {
656 nsmps -= early;
657 memset(&ar[nsmps], '\0', early*sizeof(MYFLT));
658 }
659 for (n=offset;n<nsmps;n++) {
660 /* IV - Jul 11 2002 */
661 if (p->ampcod)
662 ar[n] = (p->num1 + (MYFLT)phs * p->dfdmax) * ampp[n];
663 else
664 ar[n] = (p->num1 + (MYFLT)phs * p->dfdmax) * ampp[0];
665 phs += inc; /* phs += inc */
666 if (p->cpscod)
667 inc = (int32_t)(cpsp[n] * csound->sicvt); /* (nxt inc) */
668 if (UNLIKELY(phs >= MAXLEN)) { /* when phs o'flows */
669 phs &= PHMASK;
670 p->num1 = p->num2;
671 p->num2 = cauchrand(csound, *p->arg1);
672 p->dfdmax = (p->num2 - p->num1) / FMAXLEN;
673 }
674 }
675 p->phs = phs;
676 return OK;
677 }
678
abeta(CSOUND * csound,PRAND * p)679 int32_t abeta(CSOUND *csound, PRAND *p) /* Beta random functions */
680 {
681 uint32_t offset = p->h.insdshead->ksmps_offset;
682 uint32_t early = p->h.insdshead->ksmps_no_end;
683 uint32_t n, nsmps = CS_KSMPS;
684 MYFLT *out = p->out;
685 MYFLT arg1 = *p->arg1;
686 MYFLT arg2 = *p->arg2;
687 MYFLT arg3 = *p->arg3;
688
689 if (UNLIKELY(offset)) memset(out, '\0', offset*sizeof(MYFLT));
690 if (UNLIKELY(early)) {
691 nsmps -= early;
692 memset(&out[nsmps], '\0', early*sizeof(MYFLT));
693 }
694 for (n = 0; n < nsmps; n++)
695 out[n] = betarand(csound, arg1, arg2, arg3);
696 return OK;
697 }
698
ikbeta(CSOUND * csound,PRAND * p)699 int32_t ikbeta(CSOUND *csound, PRAND *p)
700 {
701 *p->out = betarand(csound, *p->arg1, *p->arg2, *p->arg3);
702 return OK;
703 }
704
aweib(CSOUND * csound,PRAND * p)705 int32_t aweib(CSOUND *csound, PRAND *p) /* Weibull randon functions */
706 {
707 uint32_t offset = p->h.insdshead->ksmps_offset;
708 uint32_t early = p->h.insdshead->ksmps_no_end;
709 uint32_t n, nsmps = CS_KSMPS;
710 MYFLT *out = p->out;
711 MYFLT arg1 = *p->arg1;
712 MYFLT arg2 = *p->arg2;
713
714 if (UNLIKELY(offset)) memset(out, '\0', offset*sizeof(MYFLT));
715 if (UNLIKELY(early)) {
716 nsmps -= early;
717 memset(&out[nsmps], '\0', early*sizeof(MYFLT));
718 }
719 for (n = offset; n < nsmps; n++)
720 out[n] = weibrand(csound, arg1, arg2);
721 return OK;
722 }
723
ikweib(CSOUND * csound,PRAND * p)724 int32_t ikweib(CSOUND *csound, PRAND *p)
725 {
726 *p->out = weibrand(csound, *p->arg1, *p->arg2);
727 return OK;
728 }
729
apoiss(CSOUND * csound,PRAND * p)730 int32_t apoiss(CSOUND *csound, PRAND *p) /* Poisson random funcions */
731 {
732 uint32_t offset = p->h.insdshead->ksmps_offset;
733 uint32_t early = p->h.insdshead->ksmps_no_end;
734 uint32_t n, nsmps = CS_KSMPS;
735 MYFLT *out = p->out;
736 MYFLT arg1 = *p->arg1;
737
738 if (UNLIKELY(offset)) memset(out, '\0', offset*sizeof(MYFLT));
739 if (UNLIKELY(early)) {
740 nsmps -= early;
741 memset(&out[nsmps], '\0', early*sizeof(MYFLT));
742 }
743 for (n = offset; n < nsmps; n++)
744 out[n] = poissrand(csound, arg1);
745 return OK;
746 }
747
ikpoiss(CSOUND * csound,PRAND * p)748 int32_t ikpoiss(CSOUND *csound, PRAND *p)
749 {
750 *p->out = poissrand(csound, *p->arg1);
751 return OK;
752 }
753
754 /* ------------------------------------------------------------------------ */
755
gen21_rand(FGDATA * ff,FUNC * ftp)756 int32_t gen21_rand(FGDATA *ff, FUNC *ftp)
757 {
758 CSOUND *csound = ff->csound;
759 int32_t i, n;
760 MYFLT *ft;
761 MYFLT scale;
762 int32_t nargs = ff->e.pcnt - 4;
763
764 ft = ftp->ftable;
765 scale = (nargs > 1 ? ff->e.p[6] : FL(1.0));
766 n = ff->flen;
767 if (ff->guardreq)
768 n++;
769 switch ((int32_t) ff->e.p[5]) {
770 case 1: /* Uniform distribution */
771 for (i = 0 ; i < n ; i++)
772 ft[i] = unifrand(csound, scale);
773 break;
774 case 2: /* Linear distribution */
775 for (i = 0 ; i < n ; i++)
776 ft[i] = linrand(csound, scale);
777 break;
778 case 3: /* Triangular about 0.5 */
779 for (i = 0 ; i < n ; i++)
780 ft[i] = trirand(csound, scale);
781 break;
782 case 4: /* Exponential */
783 for (i = 0 ; i < n ; i++)
784 ft[i] = exprand(csound, scale);
785 break;
786 case 5: /* Bilateral exponential */
787 for (i = 0 ; i < n ; i++)
788 ft[i] = biexprand(csound, scale);
789 break;
790 case 6: /* Gaussian distribution */
791 for (i = 0 ; i < n ; i++)
792 ft[i] = gaussrand(csound, scale);
793 break;
794 case 7: /* Cauchy distribution */
795 for (i = 0 ; i < n ; i++)
796 ft[i] = cauchrand(csound, scale);
797 break;
798 case 8: /* Positive Cauchy */
799 for (i = 0 ; i < n ; i++)
800 ft[i] = pcauchrand(csound, scale);
801 break;
802 case 9: /* Beta distribution */
803 if (UNLIKELY(nargs < 3)) {
804 return -1;
805 }
806 for (i = 0 ; i < n ; i++)
807 ft[i] = betarand(csound, scale, (MYFLT) ff->e.p[7], (MYFLT) ff->e.p[8]);
808 break;
809 case 10: /* Weibull Distribution */
810 if (UNLIKELY(nargs < 2)) {
811 return -1;
812 }
813 for (i = 0 ; i < n ; i++)
814 ft[i] = weibrand(csound, scale, (MYFLT) ff->e.p[7]);
815 break;
816 case 11: /* Poisson Distribution */
817 for (i = 0 ; i < n ; i++)
818 ft[i] = poissrand(csound, scale);
819 break;
820 default:
821 return -2;
822 }
823 return OK;
824 }
825
826 /* ---------------------------------------------- */
827
828 /* New Gauss generator using Box-Mueller transform
829 VL April 2020
830 */
831
gausscompute(CSOUND * csound,GAUSS * p)832 MYFLT gausscompute(CSOUND *csound, GAUSS *p) {
833 if(p->flag == 0) {
834 MYFLT u1 = unirand(csound);
835 MYFLT u2 = unirand(csound);
836 MYFLT z = SQRT(-2.*LOG(u1))*cos(2*PI*u2);
837 p->z = SQRT(-2.*LOG(u1))*sin(2*PI*u2);
838 p->flag = 1;
839 return *p->sigma*z + *p->mu;
840 } else {
841 p->flag = 0;
842 return *p->sigma*p->z + *p->mu;
843 }
844 return OK;
845 }
846
847
gauss_scalar(CSOUND * csound,GAUSS * p)848 int32_t gauss_scalar(CSOUND *csound, GAUSS *p){
849 *p->a = gausscompute(csound,p);
850 return OK;
851
852 }
853
gauss_vector(CSOUND * csound,GAUSS * p)854 int32_t gauss_vector(CSOUND *csound, GAUSS *p) {
855 uint32_t offset = p->h.insdshead->ksmps_offset;
856 uint32_t early = p->h.insdshead->ksmps_no_end;
857 uint32_t n, nsmps = CS_KSMPS;
858 MYFLT *out = p->a;
859 if (UNLIKELY(offset)) memset(out, '\0', offset*sizeof(MYFLT));
860 if (UNLIKELY(early)) {
861 nsmps -= early;
862 memset(&out[nsmps], '\0', early*sizeof(MYFLT));
863 }
864 for (n = offset; n < nsmps; n++)
865 out[n] = gausscompute(csound,p);
866 return OK;
867 }
868
869
870
871
872
873
874