LRgaus.c

Go to the documentation of this file.

/*
 * Copyright    2019    R.K. Owen, Ph.D.
 * License      see lgpl.md (Gnu Lesser General Public License)
 */
#ifdef __cplusplus
extern "C" {
#endif

#include <math.h>
#include "libran.h"

/* double */
double LRd_gausbm_RAN(LR_obj *o) {
        double one = 1.0, ntwo = -2.0, twopi = 2.0 * M_PI;
        double c, z1, z2;

        if (isnan(o->x.d)) {
                /* generate new pair - 1-U avoids possible overflow */
                c = sqrt(ntwo * log(one - o->ud(o)));
                z2 = twopi * o->ud(o);
                z1 = c * sin(z2);
                o->x.d = z2 = c * cos(z2);
                return o->m.d + o->s.d * z1;
        } else {
                /* return saved variate */
                c = o->x.d;
                o->x.d = NAN;
                return o->m.d + o->s.d * c;
        }
}

/* double */
double LRd_gausmar_RAN(LR_obj *o) {
        double one = 1.0, two = 2.0;
        double s, z1, z2;

        if (isnan(o->x.d)) {
                do {
                        z1 = two*o->ud(o) - one;
                        z2 = two*o->ud(o) - one;
                        s = z1*z1 + z2*z2;
                } while (s > one);
                s = sqrt(-two*log(s)/s);
                o->x.d = z2 * s;
                return o->m.d + o->s.d * z1 * s;
        } else {
                /* return saved variate */
                s = o->x.d;
                o->x.d = NAN;
                return o->m.d + o->s.d * s;
        }
}

double LRd_gaus_PDF(LR_obj *o, double x) {
        double  one = 1.0, half = .5,
                is = one / o->s.d,
                cc = half * M_2_SQRTPI * M_SQRT1_2 * is,
                xm = (x - o->m.d) * is,
                xx = xm * xm;

        return cc * exp(-half * xx);
}

double LRd_gaus_CDF(LR_obj *o, double x) {
        double  one = 1.0,
                half = 0.5,
                is = M_SQRT1_2 / o->s.d,
                xm = (x - o->m.d) * is;

        return half * (one + erf(xm));
}

/* float */
float LRf_gausbm_RAN(LR_obj *o) {
        float one = 1.0, ntwo = -2.0, twopi = 2.0 * M_PI;
        float c, z1, z2;

        if (isnan(o->x.f)) {
                /* generate new pair - 1-U avoids possible overflow */
                c = sqrt(ntwo * log(one - o->uf(o)));
                z2 = twopi * o->uf(o);
                z1 = c * sinf(z2);
                o->x.f = z2 = c * cosf(z2);
                return o->m.f + o->s.f * z1;
        } else {
                /* return saved variate */
                c = o->x.f;
                o->x.f = NAN;
                return o->m.f + o->s.f * c;
        }
}

float LRf_gausmar_RAN(LR_obj *o) {
        float one = 1.0, two = 2.0;
        float s, z1, z2;

        if (isnan(o->x.f)) {
                do {
                        z1 = two*o->ud(o) - one;
                        z2 = two*o->ud(o) - one;
                        s = z1*z1 + z2*z2;
                } while (s > one);
                s = sqrtf(-two*logf(s)/s);
                o->x.f = z2 * s;
                return o->m.f + o->s.f * z1 * s;
        } else {
                /* return saved variate */
                s = o->x.f;
                o->x.f = NAN;
                return o->m.f + o->s.f * s;
        }
}

float LRf_gaus_PDF(LR_obj *o, float x) {
        float   one = 1.0, half = .5,
                is = one / o->s.f,
                cc = half * M_2_SQRTPI * M_SQRT1_2 * is,
                xm = (x - o->m.f) * is,
                xx = xm * xm;

        return cc * expf(-half * xx);
}

float LRf_gaus_CDF(LR_obj *o, float x) {
        float   one = 1.0,
                half = 0.5,
                is = M_SQRT1_2 / o->s.f,
                xm = (x - o->m.f) * is;

        return half * (one + erff(xm));
}

#ifdef __cplusplus
}
#endif