diff options
Diffstat (limited to 'dist_normal.c')
| -rw-r--r-- | dist_normal.c | 174 |
1 files changed, 174 insertions, 0 deletions
diff --git a/dist_normal.c b/dist_normal.c new file mode 100644 index 0000000..7847745 --- /dev/null +++ b/dist_normal.c @@ -0,0 +1,174 @@ +/* + (c) 2002:cxc@web.fm + randomix: various PRNG's + code taken from: http://remus.rutgers.edu/%7Erhoads/Code/code.html + let's check it out + */ + +#include <m_pd.h> +#include <stdlib.h> +#include <math.h> + +#define SMALLEST_RANGE .0001 + +#ifndef RAND_MAX +#define RAND_MAX 2147483647 +#endif + +static int makeseed(void) +{ + static unsigned int random1_nextseed = 1489853723; + random1_nextseed = random1_nextseed * 435898247 + 938284287; + return (random1_nextseed & 0x7fffffff); +} + +static int rand_random_fl(seed) { + int q; + double state; + + /* The following parameters are recommended settings based on research + uncomment the one you want. */ + + double a = 1389796, m = RAND_MAX; + /* static double a = 950975, m = 2147483647; */ + /* static double a = 3467255, m = 21474836472; */ + /* static double a = 657618, m = 4294967291; */ + /* static double a = 93167, m = 4294967291; */ + /* static double a = 1345659, m = 4294967291; */ + + state = seed; + state *= a; + q = state / m; + state -= q*m; + return state; +} + +/* -------------------------- dist_normal ------------------------------ */ +/* Generate a normal random variable with mean 0 and standard deviation + of 1. To adjust to some other distribution, multiply by the standard + deviation and add the mean. Box-Muller method + note: rand() is a function that returns a uniformly distributed random + number from 0 to RAND_MAX +*/ + +static t_class *dist_normal_class; + +typedef struct _dist_normal +{ + t_object x_obj; + t_float x_mn; // mean + t_float x_dv; // deviation + t_float x_u1; + t_float x_u2; + + t_float x_f; // lower limit + t_float x_g; // upper limit + unsigned int x_state; // current seed +} t_dist_normal; + +static void *dist_normal_new(t_floatarg mn, t_floatarg dv) +{ + t_dist_normal *x = (t_dist_normal *)pd_new(dist_normal_class); + x->x_mn = (mn) ? mn : 0; + x->x_dv = (dv) ? dv : 1; + x->x_u1 = 13; + x->x_u2 = 1000; + x->x_f = 0; + x->x_g = RAND_MAX; + //post("cxc/randomix.c: lolim: %f - %f, uplim: %f - %f", x->x_f, f, x->x_g, g); + x->x_state = makeseed(); + inlet_new(&x->x_obj, &x->x_obj.ob_pd, &s_float, gensym("fl1")); + // inlet_new(&x->x_obj, &x->x_obj.ob_pd, &s_float, gensym("fl2")); + //outlet_new(&x->x_obj, &s_float); + outlet_new(&x->x_obj, 0); + return (x); +} + +static double dist_normal_rand(t_dist_normal *x) +{ + int n = x->x_f; + double nval; + double m; + + m = (double)RAND_MAX * 2; + + // post("cxc/randomix.c: x_state: %d",x->x_state); + x->x_state = rand_random_fl(x->x_state); + //nval = ((double)x->x_state / m) * (double)(x->x_g - x->x_f) + (double)x->x_f; + nval = (double)x->x_state; + return (nval); +} + +static void dist_normal_doit(t_dist_normal *x) +{ + static double V2, fac; + static int phase = 0; + double S, Z, U1, U2, V1; + + if (phase) + Z = V2 * fac; + else + { + do { +/* U1 = (double)rand() / RAND_MAX; */ +/* U2 = (double)rand() / RAND_MAX; */ + U1 = (double)dist_normal_rand(x) / RAND_MAX; + U2 = (double)dist_normal_rand(x) / RAND_MAX; + + // post("cxc/randomix.c: test %f %f %f %f", x->x_u1, x->x_u2, U1, U2); + + V1 = 2 * U1 - 1; + V2 = 2 * U2 - 1; + S = V1 * V1 + V2 * V2; + } while(S >= 1); + + fac = sqrt (-2 * log(S) / S); + Z = V1 * fac; + } + + phase = 1 - phase; + + //return Z; + outlet_float(x->x_obj.ob_outlet, Z); +} + +static void dist_normal_bang(t_dist_normal *x) +{ +/* post("cxc/randomix.c: dist_normal banged"); */ +/* post("cxc/randomix.c: RAND_MAX: %d",RAND_MAX); */ +/* post("cxc/randomix.c: test: %f %f", x->x_u1, x->x_u2); */ + dist_normal_doit(x); +} + +void dist_normal_low(t_dist_normal *x, t_floatarg mn) +{ + x->x_mn = mn; +} + +void dist_normal_upp(t_dist_normal *x, t_floatarg dv) +{ + x->x_dv = dv; +} + +void dist_normal_float(t_dist_normal *x, t_floatarg r) +{ + outlet_float(x->x_obj.ob_outlet, r); +} + +void dist_normal_list(t_dist_normal *x, t_symbol *s, int argc, t_atom *argv) +{ + outlet_list(x->x_obj.ob_outlet, s, argc, argv); +} + +void dist_normal_setup(void) +{ + dist_normal_class = class_new(gensym("dist_normal"), (t_newmethod)dist_normal_new, 0, + sizeof(t_dist_normal), 0, A_DEFFLOAT, A_DEFFLOAT, 0); + class_addbang(dist_normal_class, dist_normal_bang); + class_addmethod(dist_normal_class, (t_method)dist_normal_low, gensym("fl1"), A_FLOAT, 0); + // class_addmethod(dist_normal_class, (t_method)dist_normal_upp, gensym("fl2"), A_FLOAT, 0); + class_addlist (dist_normal_class, dist_normal_list); + +/* class_addmethod(dist_normal_class, (t_method)dist_normal_seed, */ +/* gensym("seed"), A_FLOAT, 0); */ +} |