diff --git a/src/include/clProbDist/clProbDist.h b/src/include/clProbDist/clProbDist.h index af8e077..17a5f0d 100644 --- a/src/include/clProbDist/clProbDist.h +++ b/src/include/clProbDist/clProbDist.h @@ -869,47 +869,47 @@ typedef enum clprobdistStatus_ { extern "C" { #endif - /*! @brief Retrieve the last error message. - * - * The buffer containing the error message is internally allocated and must - * not be freed by the client. - * - * @return Error message or `NULL`. - */ - const char* clprobdistGetErrorString(); +/*! @brief Retrieve the last error message. +* +* The buffer containing the error message is internally allocated and must +* not be freed by the client. +* +* @return Error message or `NULL`. +*/ +const char* clprobdistGetErrorString(); - /*! @brief Generate an include option string for use with the OpenCL C compiler - * - * Generate and return "-I${CLPROBDIST_ROOT}/include", where \c ${CLPROBDIST_ROOT} is - * the value of the \c CLPROBDIST_ROOT environment variable. - * This string is meant to be passed as an option to the OpenCL C compiler for - * programs that make use of the clProbDist device-side headers. - * If the \c CLPROBDIST_ROOT environment variable is not defined, it - * defaults to `/usr` if the file - * `/usr/include/clProbDist/clProbDist.h` exists, else to the current - * directory of execution of the program. - * - * A static buffer is return and need not be released; it could change upon - * successive calls to the function. - * - * An error is returned in \c err if the preallocated buffer is too small to - * contain the include string. - * - * @param[out] err Error status variable, or `NULL`. - * - * @return An OpenCL C compiler option to indicate where to find the - * device-side clProbDist headers. - */ - const char* clprobdistGetLibraryDeviceIncludes(cl_int* err); +/*! @brief Generate an include option string for use with the OpenCL C compiler +* +* Generate and return "-I${CLPROBDIST_ROOT}/include", where \c ${CLPROBDIST_ROOT} is +* the value of the \c CLPROBDIST_ROOT environment variable. +* This string is meant to be passed as an option to the OpenCL C compiler for +* programs that make use of the clProbDist device-side headers. +* If the \c CLPROBDIST_ROOT environment variable is not defined, it + * defaults to `/usr` if the file + * `/usr/include/clProbDist/clProbDist.h` exists, else to the current +* directory of execution of the program. +* +* A static buffer is return and need not be released; it could change upon +* successive calls to the function. +* +* An error is returned in \c err if the preallocated buffer is too small to +* contain the include string. +* +* @param[out] err Error status variable, or `NULL`. +* +* @return An OpenCL C compiler option to indicate where to find the +* device-side clProbDist headers. +*/ +const char* clprobdistGetLibraryDeviceIncludes(cl_int* err); - /*! @brief Retrieve the library installation path - * - * @return Value of the CLPROBDIST_ROOT environment variable, if - * defined; else, `/usr` if the file - * `/usr/include/clProbDist/clProbDist.h` exists; or, the current - * directory (.) of execution of the program otherwise. - */ - const char* clprobdistGetLibraryRoot(); +/*! @brief Retrieve the library installation path +* + * @return Value of the CLPROBDIST_ROOT environment variable, if + * defined; else, `/usr` if the file + * `/usr/include/clProbDist/clProbDist.h` exists; or, the current +* directory (.) of execution of the program otherwise. +*/ +const char* clprobdistGetLibraryRoot(); #ifdef __cplusplus } @@ -917,4 +917,4 @@ extern "C" { #endif /* CLPROBDIST_H * * vim: syntax=c.doxygen spell spelllang=en fdm=syntax fdls=0 expandtab -*/ +*/ \ No newline at end of file diff --git a/src/include/clProbDist/clProbDist_template.h b/src/include/clProbDist/clProbDist_template.h index 8fc76c5..b4cf6bd 100644 --- a/src/include/clProbDist/clProbDist_template.h +++ b/src/include/clProbDist/clProbDist_template.h @@ -30,6 +30,10 @@ #include +#ifdef __cplusplus +extern "C" +{ +#endif /*! @file clProbDist_template.h * @brief Template of the API for specific probability distributions (not to be included as is!) @@ -66,7 +70,7 @@ * * ### Host and Device APIs * - * The functions described here are all available on the host, in all implementations, + * The functions described here are all available on the host, in all implementations, * unless specified otherwise. Only some of the functions and types are also * available on the device in addition to the host; they are tagged with * [**device**]. @@ -245,7 +249,7 @@ cl_double clprobdistStdDeviationWithObject(const clprobdistObject* dist, clprobd * * Copy the distribution object @c srcDist located in global memory into the * buffer @c destDist located in local or private memory. - * This function *does not* allocate memory for the distribution object, + * This function *does not* allocate memory for the distribution object, * it assumes that this has already been done. * * When the distribution API is configured to use distribution objects stored @@ -278,7 +282,7 @@ clprobdistStatus clprobdistCopyOverFromGlobal(_CLPROBDIST__OBJ_MEM clprobd * They take the distribution parameters as their first arguments * (see the description of \ref DIST_PARAMS in \ref distributions). * They cannot take advantage of precomputed values for the distributions. - * + * * @{ */ @@ -389,6 +393,10 @@ cl_double clprobdistStdDeviation(DIST_PARAMS, clprobdistStatus* err); /*! @} */ +#ifdef __cplusplus +} +#endif + #endif /* CLPROBDIST_H * * vim: syntax=c.doxygen spell spelllang=en fdm=syntax fdls=0 expandtab */ diff --git a/src/include/clProbDist/exponential.h b/src/include/clProbDist/exponential.h index 48c9f1a..a6cad40 100644 --- a/src/include/clProbDist/exponential.h +++ b/src/include/clProbDist/exponential.h @@ -36,6 +36,11 @@ #include "clProbDist/clProbDist.h" #include "clProbDist/continuous.h" +#ifdef __cplusplus +extern "C" +{ +#endif + /*! @brief Exponential distribution object [**device**] * * A structure that represents an exponential distribution object. @@ -59,6 +64,7 @@ typedef struct _clprobdistExponential clprobdistExponential; * @param[out] err Error status variable, or \c NULL. * @return New distribution object. */ + clprobdistExponential* clprobdistExponentialCreate(cl_double lambda, size_t* bufSize, clprobdistStatus* err); /*! @copydoc clprobdistDestroy() @@ -194,5 +200,8 @@ cl_double clprobdistExponentialStdDeviation(cl_double lambda, clprobdistStatus* /*! @} */ +#ifdef __cplusplus +} +#endif #endif diff --git a/src/include/clProbDist/gamma.h b/src/include/clProbDist/gamma.h index be26bca..c2fce36 100644 --- a/src/include/clProbDist/gamma.h +++ b/src/include/clProbDist/gamma.h @@ -44,6 +44,11 @@ #include "clProbDist/clProbDist.h" #include "clProbDist/continuous.h" +#ifdef __cplusplus +extern "C" +{ +#endif + /*! @brief Gamma distribution object [**device**] * * A structure that represents a gamma distribution object. @@ -252,4 +257,8 @@ cl_double clprobdistGammaComplCDF_1(cl_double alpha, int d, cl_double x, clprobd */ cl_double clprobdistGammaInverseCDF_1(cl_double alpha, int d, cl_double u, clprobdistStatus* err); +#ifdef __cplusplus +} +#endif + #endif diff --git a/src/include/clProbDist/lognormal.h b/src/include/clProbDist/lognormal.h index ad21db6..098fa99 100644 --- a/src/include/clProbDist/lognormal.h +++ b/src/include/clProbDist/lognormal.h @@ -37,6 +37,11 @@ #include "clProbDist/clProbDist.h" #include "clProbDist/continuous.h" +#ifdef __cplusplus +extern "C" +{ +#endif + /*! @brief Lognormal distribution object [**device**] * * A structure that represents a lognormal distribution object. @@ -221,6 +226,8 @@ cl_double clprobdistLognormalStdDeviation(cl_double mu, cl_double sigma, clprobd /*! @} */ - +#ifdef __cplusplus +} +#endif #endif diff --git a/src/include/clProbDist/normal.h b/src/include/clProbDist/normal.h index b264a6d..9a3a5d8 100644 --- a/src/include/clProbDist/normal.h +++ b/src/include/clProbDist/normal.h @@ -36,6 +36,11 @@ #include "clProbDist/clProbDist.h" #include "clProbDist/continuous.h" +#ifdef __cplusplus +extern "C" +{ +#endif + /*! @brief Normal distribution object [**device**] * * A structure that represents a normal distribution object. @@ -246,6 +251,8 @@ cl_double clprobdistNormalStdDeviation(cl_double mu, cl_double sigma, clprobdist /*! @} */ - +#ifdef __cplusplus +} +#endif #endif diff --git a/src/include/clProbDist/poisson.h b/src/include/clProbDist/poisson.h index c01a2d5..5225726 100644 --- a/src/include/clProbDist/poisson.h +++ b/src/include/clProbDist/poisson.h @@ -36,6 +36,10 @@ #include "clProbDist/clProbDist.h" #include +#ifdef __cplusplus +extern "C" +{ +#endif /*! @brief Poisson distribution object [**device**] * diff --git a/src/include/clProbDist/private/continuous.c.h b/src/include/clProbDist/private/continuous.c.h index 8d21b64..be542b2 100644 --- a/src/include/clProbDist/private/continuous.c.h +++ b/src/include/clProbDist/private/continuous.c.h @@ -61,7 +61,7 @@ clprobdistStatus clprobdistContinuousSetXsup(cl_double xb, clprobdistContinuous* /*********************************** * Abstract functions ***********************************/ -cl_double clprobdistContinuousCDF(cl_double x, clprobdistStatus* err){ +cl_double clprobdistContinuousCDF(cl_double x, clprobdistStatus* err) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): Abstract function", __func__); return -1; } @@ -75,7 +75,7 @@ clprobdistStatus clprobdistContinuousFindInterval(cl_double u, cl_double* iv, cl const cl_double XLIM = DBL_MAX / 2.0; const cl_double B0 = 8.0; - + cl_double b = B0; while (b < XLIM && u > clprobdistContinuousCDF(b, err)) b *= 2.0; @@ -86,7 +86,7 @@ clprobdistStatus clprobdistContinuousFindInterval(cl_double u, cl_double* iv, cl } cl_double a = -B0; - while (a > -XLIM && u < clprobdistContinuousCDF(a,err)) + while (a > -XLIM && u < clprobdistContinuousCDF(a, err)) a *= 2.0; if (a < -B0) { iv[1] = a / 2.0; @@ -99,10 +99,10 @@ clprobdistStatus clprobdistContinuousFindInterval(cl_double u, cl_double* iv, cl return CLPROBDIST_SUCCESS; } -cl_double clprobdistContinuousDensity(cl_double x, clprobdistStatus* err){ +cl_double clprobdistContinuousDensity(cl_double x, clprobdistStatus* err) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): Abstract function", __func__); - return -1; -} + return -1; +} cl_double clprobdistContinuousGetMean(clprobdistStatus* err) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): Abstract function", __func__); @@ -129,7 +129,7 @@ cl_double clprobdistContinuousComplCDF(cl_double x, clprobdistStatus* err) { } cl_double clprobdistContinuousInverseBrent(cl_double a, cl_double b, cl_double u, cl_double tol, clprobdistContinuous* distObj, clprobdistStatus* err) { - if (u > 1.0 || u < 0.0){ + if (u > 1.0 || u < 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): u not in [0, 1]", __func__); return -1; } @@ -146,15 +146,15 @@ cl_double clprobdistContinuousInverseBrent(cl_double a, cl_double b, cl_double u } int MAXITER = 50; // Maximum number of iterations tol += clprobdistEPSARRAY[distObj->decPrec] + DBL_EPSILON; // in case tol is too small - + cl_double ua = clprobdistContinuousCDF(a, err) - u; - if (ua > 0.0){ + if (ua > 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): u < cdf(a)", __func__); return -1; } cl_double ub = clprobdistContinuousCDF(b, err) - u; - if (ub < 0.0){ + if (ub < 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): u > cdf(b)", __func__); return -1; } @@ -179,8 +179,8 @@ cl_double clprobdistContinuousInverseBrent(cl_double a, cl_double b, cl_double u } int i; for (i = 0; i < MAXITER; ++i) { - cl_double tol1 = tol + 4.0*DBL_EPSILON*fabs(b); - cl_double xm = 0.5*(c - b); + cl_double tol1 = tol + 4.0 * DBL_EPSILON * fabs(b); + cl_double xm = 0.5 * (c - b); /*if (DEBUG) { System.out.println(PrintfFormat.d(3, i) + " " + PrintfFormat.g(18, decPrec, b) + " " + @@ -200,32 +200,29 @@ cl_double clprobdistContinuousInverseBrent(cl_double a, cl_double b, cl_double u s = ub / ua; q = 1.0 - s; p = 2.0 * xm * s; - } - else { + } else { // quadratic interpolation q = ua / uc; r = ub / uc; s = ub / ua; - p = s*(2.0*xm*q*(q - r) - (b - a)*(r - 1.0)); - q = (q - 1.0)*(r - 1.0)* (s - 1.0); + p = s * (2.0 * xm * q * (q - r) - (b - a) * (r - 1.0)); + q = (q - 1.0) * (r - 1.0) * (s - 1.0); } if (p > 0.0) q = -q; p = fabs(p); // Accept interpolation? - if ((2.0*p >= (3.0*xm*q - fabs(q*tol1))) || - (p >= fabs(0.5*t*q))) { + if ((2.0 * p >= (3.0 * xm * q - fabs(q * tol1))) || + (p >= fabs(0.5 * t * q))) { len = xm; t = len; - } - else { + } else { t = len; len = p / q; } - } - else { + } else { len = xm; t = len; } @@ -241,13 +238,12 @@ cl_double clprobdistContinuousInverseBrent(cl_double a, cl_double b, cl_double u ub = clprobdistContinuousCDF(b, err) - u; - if (ub*(uc / fabs(uc)) > 0.0) { + if (ub * (uc / fabs(uc)) > 0.0) { c = a; uc = ua; len = b - a; t = len; - } - else if (fabs(uc) < fabs(ub)) { + } else if (fabs(uc) < fabs(ub)) { a = b; b = c; c = a; ua = ub; ub = uc; uc = ua; } @@ -256,23 +252,23 @@ cl_double clprobdistContinuousInverseBrent(cl_double a, cl_double b, cl_double u /* String lineSep = System.getProperty("line.separator"); System.out.println(lineSep +"*********** inverseBrent: no convergence after " + MAXITER +" iterations");*/ } - return b; + return b; } cl_double clprobdistContinuousInverseBisection(cl_double u, clprobdistContinuous* distObj, clprobdistStatus* err) { int MAXITER = 100; // Maximum number of iterations cl_double EPSILON = clprobdistEPSARRAY[distObj->decPrec]; // Absolute precision - if (u > 1.0 || u < 0.0){ + if (u > 1.0 || u < 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): u not in [0, 1]", __func__); return -1; } - + if (distObj->decPrec > DBL_DIG) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): decPrec too large", __func__); return -1; } - + if (distObj->decPrec <= 0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): decPrec <= 0", __func__); return -1; @@ -289,7 +285,7 @@ cl_double clprobdistContinuousInverseBisection(cl_double u, clprobdistContinuous return x; } - cl_double iv[2] = { 0, 0 }; + cl_double iv[2] = { 0, 0 }; clprobdistContinuousFindInterval(u, iv, distObj, err); cl_double xa = iv[0]; cl_double xb = iv[1]; @@ -302,16 +298,13 @@ cl_double clprobdistContinuousInverseBisection(cl_double u, clprobdistContinuous while (!fini) { x = (xa + xb) / 2.0; - y = clprobdistContinuousCDF(x,err) - u; + y = clprobdistContinuousCDF(x, err) - u; if ((y == 0.0) || - (fabs((xb - xa) / (x + DBL_EPSILON)) <= EPSILON)) - { + (fabs((xb - xa) / (x + DBL_EPSILON)) <= EPSILON)) { fini = CL_TRUE; - } - else if (y*ya < 0.0){ + } else if (y * ya < 0.0) { xb = x; - } - else + } else xa = x; ++i; @@ -332,4 +325,4 @@ cl_double clprobdistContinuousInverseCDF(cl_double u, clprobdistContinuous* dist } #endif -#endif +#endif \ No newline at end of file diff --git a/src/include/clProbDist/private/exponential.c.h b/src/include/clProbDist/private/exponential.c.h index 910e617..a3f8be3 100644 --- a/src/include/clProbDist/private/exponential.c.h +++ b/src/include/clProbDist/private/exponential.c.h @@ -30,6 +30,11 @@ #define _CLPROBDIST_EXPONENTIAL_OBJ_MEM #endif +#ifdef __cplusplus +extern "C" +{ +#endif + struct _clprobdistExponential { clprobdistContinuous continuousDistObj; cl_double lambda; @@ -39,8 +44,7 @@ struct _clprobdistExponential { constant const cl_double clprobdistExponentialXBIG = 100.0; constant const cl_double clprobdistExponentialXBIGM = 1000.0; -clprobdistStatus clprobdistExponentialSetLambda(_CLPROBDIST_EXPONENTIAL_OBJ_MEM clprobdistExponential* distObj, cl_double newlambda) -{ +clprobdistStatus clprobdistExponentialSetLambda(_CLPROBDIST_EXPONENTIAL_OBJ_MEM clprobdistExponential* distObj, cl_double newlambda) { clprobdistStatus err = CLPROBDIST_SUCCESS; if (distObj->lambda <= 0) @@ -52,11 +56,10 @@ clprobdistStatus clprobdistExponentialSetLambda(_CLPROBDIST_EXPONENTIAL_OBJ_MEM return err; } -cl_double clprobdistExponentialDensity(cl_double lambda, cl_double x, clprobdistStatus* err) -{ +cl_double clprobdistExponentialDensity(cl_double lambda, cl_double x, clprobdistStatus* err) { if (err) *err = CLPROBDIST_SUCCESS; - if (lambda <= 0){ + if (lambda <= 0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): lambda <= 0", __func__); return -1; } @@ -64,11 +67,10 @@ cl_double clprobdistExponentialDensity(cl_double lambda, cl_double x, clprobdist return x < 0 ? 0 : lambda * exp(-lambda * x); } -cl_double clprobdistExponentialCDF(cl_double lambda, cl_double x, clprobdistStatus* err) -{ +cl_double clprobdistExponentialCDF(cl_double lambda, cl_double x, clprobdistStatus* err) { if (err) *err = CLPROBDIST_SUCCESS; - if (lambda <= 0){ + if (lambda <= 0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): lambda <= 0", __func__); return -1; } @@ -83,27 +85,26 @@ cl_double clprobdistExponentialCDF(cl_double lambda, cl_double x, clprobdistStat return -expm1(-y); } -cl_double clprobdistExponentialComplCDF(cl_double lambda, cl_double x, clprobdistStatus* err) -{ +cl_double clprobdistExponentialComplCDF(cl_double lambda, cl_double x, clprobdistStatus* err) { if (err) *err = CLPROBDIST_SUCCESS; - if (lambda <= 0){ + if (lambda <= 0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): lambda <= 0", __func__); return -1; } if (x <= 0.0) return 1.0; - if (lambda*x >= clprobdistExponentialXBIGM) + if (lambda * x >= clprobdistExponentialXBIGM) return 0.0; - return exp(-lambda*x); + return exp(-lambda * x); } -cl_double clprobdistExponentialInverseCDF(cl_double lambda, cl_double u, clprobdistStatus* err){ +cl_double clprobdistExponentialInverseCDF(cl_double lambda, cl_double u, clprobdistStatus* err) { if (err) *err = CLPROBDIST_SUCCESS; - if (lambda <= 0){ + if (lambda <= 0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): lambda <= 0", __func__); return -1; } @@ -122,10 +123,9 @@ cl_double clprobdistExponentialInverseCDF(cl_double lambda, cl_double u, clprobd return -log1p(-u) / lambda; // log1p is defined on OpenCL C } -cl_double clprobdistExponentialMean(cl_double lambda, clprobdistStatus* err) -{ +cl_double clprobdistExponentialMean(cl_double lambda, clprobdistStatus* err) { if (err) *err = CLPROBDIST_SUCCESS; - if (lambda <= 0.0){ + if (lambda <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): lambda <= 0", __func__); return -1; } @@ -133,10 +133,9 @@ cl_double clprobdistExponentialMean(cl_double lambda, clprobdistStatus* err) return (1 / lambda); } -cl_double clprobdistExponentialVariance(cl_double lambda, clprobdistStatus* err) -{ +cl_double clprobdistExponentialVariance(cl_double lambda, clprobdistStatus* err) { if (err) *err = CLPROBDIST_SUCCESS; - if (lambda <= 0.0){ + if (lambda <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): lambda <= 0", __func__); return -1; } @@ -144,10 +143,9 @@ cl_double clprobdistExponentialVariance(cl_double lambda, clprobdistStatus* err) return (1 / (lambda * lambda)); } -cl_double clprobdistExponentialStdDeviation(cl_double lambda, clprobdistStatus* err) -{ +cl_double clprobdistExponentialStdDeviation(cl_double lambda, clprobdistStatus* err) { *err = CLPROBDIST_SUCCESS; - if (lambda <= 0.0){ + if (lambda <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): lambda <= 0", __func__); return -1; } @@ -156,8 +154,7 @@ cl_double clprobdistExponentialStdDeviation(cl_double lambda, clprobdistStatus* } -cl_double clprobdistExponentialGetLambda(_CLPROBDIST_EXPONENTIAL_OBJ_MEM const clprobdistExponential* distObj, clprobdistStatus* err) -{ +cl_double clprobdistExponentialGetLambda(_CLPROBDIST_EXPONENTIAL_OBJ_MEM const clprobdistExponential* distObj, clprobdistStatus* err) { if (err) *err = CLPROBDIST_SUCCESS; return distObj->lambda; } @@ -180,17 +177,20 @@ cl_double clprobdistExponentialInverseCDFWithObject(_CLPROBDIST_EXPONENTIAL_OBJ_ } -cl_double clprobdistExponentialMeanWithObject(_CLPROBDIST_EXPONENTIAL_OBJ_MEM const clprobdistExponential* distObj, clprobdistStatus* err){ - return clprobdistExponentialMean(distObj->lambda, err); +cl_double clprobdistExponentialMeanWithObject(_CLPROBDIST_EXPONENTIAL_OBJ_MEM const clprobdistExponential* distObj, clprobdistStatus* err) { + return clprobdistExponentialMean(distObj->lambda, err); } -cl_double clprobdistExponentialVarianceWithObject(_CLPROBDIST_EXPONENTIAL_OBJ_MEM const clprobdistExponential* distObj, clprobdistStatus* err){ - return clprobdistExponentialVariance(distObj->lambda, err); +cl_double clprobdistExponentialVarianceWithObject(_CLPROBDIST_EXPONENTIAL_OBJ_MEM const clprobdistExponential* distObj, clprobdistStatus* err) { + return clprobdistExponentialVariance(distObj->lambda, err); } -cl_double clprobdistExponentialStdDeviationWithObject(_CLPROBDIST_EXPONENTIAL_OBJ_MEM const clprobdistExponential* distObj, clprobdistStatus* err){ - return clprobdistExponentialStdDeviation(distObj->lambda, err); +cl_double clprobdistExponentialStdDeviationWithObject(_CLPROBDIST_EXPONENTIAL_OBJ_MEM const clprobdistExponential* distObj, clprobdistStatus* err) { + return clprobdistExponentialStdDeviation(distObj->lambda, err); } +#ifdef __cplusplus +} +#endif #endif diff --git a/src/include/clProbDist/private/gamma.c.h b/src/include/clProbDist/private/gamma.c.h index 1245feb..561393f 100644 --- a/src/include/clProbDist/private/gamma.c.h +++ b/src/include/clProbDist/private/gamma.c.h @@ -30,9 +30,14 @@ #define _CLPROBDIST_GAMMA_OBJ_MEM #endif +#ifdef __cplusplus +extern "C" +{ +#endif + #define signum(a) (((a) < 0) ? -1 : (((a) > 0) ? 1 : 0)) -typedef struct _clprobdistGammaParams{ +typedef struct _clprobdistGammaParams { cl_double alpha; cl_double lambda; cl_double logFactor; @@ -61,18 +66,17 @@ constant const cl_double clprobdistGammaEPSARRAY[] = { /* * This is the function 1 + (1 - x*x + 2*x*log(x)) / ((1 - x)*(1 - x)) */ -static cl_double clprobdistGamma_mybelog(cl_double x) -{ +static cl_double clprobdistGamma_mybelog(cl_double x) { if (x < 1.0e-30) return 0.0; if (x > 1.0e30) - return 2.0*(log(x) - 1.0) / x; + return 2.0 * (log(x) - 1.0) / x; if (x == 1.0) return 1.0; cl_double t = 1.0 - x; if (x < 0.9 || x > 1.1) { - cl_double w = (t + x*log(x)) / (t*t); + cl_double w = (t + x * log(x)) / (t * t); return 2.0 * w; } @@ -88,7 +92,7 @@ static cl_double clprobdistGamma_mybelog(cl_double x) sum += term; j++; } while (fabs(term / sum) > EPS); - return 2.0*sum; + return 2.0 * sum; } //************************************ @@ -188,8 +192,7 @@ cl_double clprobdistRootFinderBrentDekker(cl_double a, cl_double b, clprobdistRo s = fb / fa; p = s * (2.0 * xm * q * (q - r) - (b - a) * (r - 1.0)); q = (q - 1.0) * (r - 1.0) * (s - 1.0); - } - else { + } else { // Linear interpolation s = fb / fa; p = 2.0 * xm * s; @@ -203,16 +206,14 @@ cl_double clprobdistRootFinderBrentDekker(cl_double a, cl_double b, clprobdistRo // Is interpolation acceptable ? if (((2.0 * p) >= (3.0 * xm * q - fabs(tol1 * q))) - || (p >= fabs(0.5 * e * q))) { + || (p >= fabs(0.5 * e * q))) { d = xm; e = d; - } - else { + } else { e = d; d = p / q; } - } - else { + } else { // Bisection necessary d = xm; e = d; @@ -226,13 +227,12 @@ cl_double clprobdistRootFinderBrentDekker(cl_double a, cl_double b, clprobdistRo b -= tol1; else b += tol1; - fb = clprobdistRootFinderFuncEvaluate(b,f, err); + fb = clprobdistRootFinderFuncEvaluate(b, f, err); if (fb * (signum(fc)) > 0.0) { c = a; fc = fa; d = e = b - a; - } - else { + } else { a = b; b = c; c = a; @@ -274,7 +274,7 @@ cl_double clprobdistRootFinderBisection(cl_double a, cl_double b, clprobdistRoot } cl_double xa = a; cl_double xb = b; - cl_double ya = clprobdistRootFinderFuncEvaluate(a,f, err); + cl_double ya = clprobdistRootFinderFuncEvaluate(a, f, err); cl_double x = 0, y = 0; int MAXITER = 1200; // Maximum number of iterations //bool DEBUG = false; @@ -286,10 +286,10 @@ cl_double clprobdistRootFinderBisection(cl_double a, cl_double b, clprobdistRoot int i = 0; while (!fini) { x = (xa + xb) / 2.0; - y = clprobdistRootFinderFuncEvaluate(x,f, err); + y = clprobdistRootFinderFuncEvaluate(x, f, err); if ((fabs(y) <= clprobdistGamma_MINVAL) || - (fabs(xb - xa) <= tol * fabs(x)) || - (fabs(xb - xa) <= clprobdistGamma_MINVAL)) { + (fabs(xb - xa) <= tol * fabs(x)) || + (fabs(xb - xa) <= clprobdistGamma_MINVAL)) { if (fabs(x) > clprobdistGamma_MINVAL) return x; else @@ -312,15 +312,14 @@ cl_double clprobdistRootFinderBisection(cl_double a, cl_double b, clprobdistRoot //Static functions -cl_double clprobdistGammaDensity(cl_double alpha, cl_double lambda, int decprec, cl_double x, clprobdistStatus* err) -{ +cl_double clprobdistGammaDensity(cl_double alpha, cl_double lambda, int decprec, cl_double x, clprobdistStatus* err) { if (err) *err = CLPROBDIST_SUCCESS; - if (alpha <= 0.0){ + if (alpha <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): alpha <= 0", __func__); return -1; } - if (lambda <= 0.0){ + if (lambda <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): lambda <= 0", __func__); return -1; } @@ -328,7 +327,7 @@ cl_double clprobdistGammaDensity(cl_double alpha, cl_double lambda, int decprec, if (x <= 0) return 0.0; - cl_double z = alpha * log(lambda*x) - lambda*x - lgamma(alpha); + cl_double z = alpha * log(lambda * x) - lambda * x - lgamma(alpha); if (z > -clprobdistGammaXBIGM) return exp(z) / x; else @@ -337,12 +336,12 @@ cl_double clprobdistGammaDensity(cl_double alpha, cl_double lambda, int decprec, cl_double clprobdistGammaStdDeviation(cl_double alpha, cl_double lambda, int decprec, clprobdistStatus* err) { if (err) *err = CLPROBDIST_SUCCESS; - if (alpha <= 0.0){ + if (alpha <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): alpha <= 0", __func__); return -1; } - if (lambda <= 0.0){ + if (lambda <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): lambda <= 0", __func__); return -1; } @@ -351,12 +350,12 @@ cl_double clprobdistGammaStdDeviation(cl_double alpha, cl_double lambda, int dec } cl_double clprobdistGammaInverseCDF_1(cl_double alpha, int d, cl_double u, clprobdistStatus* err) { if (err) *err = CLPROBDIST_SUCCESS; - if (alpha <= 0.0){ + if (alpha <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): alpha <= 0", __func__); return -1; } - if (u > 1.0 || u < 0.0){ + if (u > 1.0 || u < 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): u not in [0,1]", __func__); return -1; } @@ -367,7 +366,7 @@ cl_double clprobdistGammaInverseCDF_1(cl_double alpha, int d, cl_double u, clpro if (u >= 1.0) return DBL_MAX; - if (d <= 0){ + if (d <= 0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): d <= 0", __func__); return -1; } @@ -400,8 +399,7 @@ cl_double clprobdistGammaInverseCDF_1(cl_double alpha, int d, cl_double u, clpro return clprobdistRootFinderBisection(x, xmax, &f, EPS, err); else return clprobdistRootFinderBisection(0, x, &f, EPS, err); - } - else { + } else { if (cdf < u) return clprobdistRootFinderBrentDekker(x, xmax, &f, EPS, err); else @@ -412,12 +410,12 @@ cl_double clprobdistGammaInverseCDF_1(cl_double alpha, int d, cl_double u, clpro cl_double clprobdistGammaMean(cl_double alpha, cl_double lambda, int decprec, clprobdistStatus* err) { if (err) *err = CLPROBDIST_SUCCESS; - if (alpha <= 0.0){ + if (alpha <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): alpha <= 0", __func__); return -1; } - if (lambda <= 0.0){ + if (lambda <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): lambda <= 0", __func__); return -1; } @@ -427,12 +425,12 @@ cl_double clprobdistGammaMean(cl_double alpha, cl_double lambda, int decprec, cl } cl_double clprobdistGammaVariance(cl_double alpha, cl_double lambda, int decprec, clprobdistStatus* err) { if (err) *err = CLPROBDIST_SUCCESS; - if (alpha <= 0.0){ + if (alpha <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): alpha <= 0", __func__); return -1; } - if (lambda <= 0.0){ + if (lambda <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): lambda <= 0", __func__); return -1; } @@ -443,12 +441,12 @@ cl_double clprobdistGammaVariance(cl_double alpha, cl_double lambda, int decprec cl_double clprobdistGammaComplCDF_1(cl_double alpha, int d, cl_double x, clprobdistStatus* err) { if (err) *err = CLPROBDIST_SUCCESS; - if (alpha <= 0.0){ + if (alpha <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): alpha <= 0", __func__); return -1; } - if (d <= 0){ + if (d <= 0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): d <= 0", __func__); return -1; } @@ -461,8 +459,7 @@ cl_double clprobdistGammaComplCDF_1(cl_double alpha, int d, cl_double x, clprobd if (alpha >= 70.0) { if (x >= alpha * clprobdistGammaXBIG) return 0.0; - } - else { + } else { if (x >= clprobdistGammaXBIGM) return 0.0; } @@ -484,7 +481,7 @@ cl_double clprobdistGammaComplCDF_1(cl_double alpha, int d, cl_double x, clprobd cl_double RENORM = 1.0E100; cl_double R, dif; int i; - cl_double factor = exp(alpha*log(x) - x - lgamma(alpha)); + cl_double factor = exp(alpha * log(x) - x - lgamma(alpha)); cl_double A = 1.0 - alpha; cl_double B = A + x + 1.0; @@ -504,8 +501,8 @@ cl_double clprobdistGammaComplCDF_1(cl_double alpha, int d, cl_double x, clprobd if (V[5] != 0.0) { R = V[4] / V[5]; dif = fabs(res - R); - if (dif <= EPS*R){ - return factor*res; + if (dif <= EPS * R) { + return factor * res; } res = R; @@ -519,17 +516,17 @@ cl_double clprobdistGammaComplCDF_1(cl_double alpha, int d, cl_double x, clprobd } while (CL_TRUE); } cl_double clprobdistGammaComplCDF(cl_double alpha, cl_double lambda, int d, cl_double x, clprobdistStatus* err) { - return clprobdistGammaComplCDF_1(alpha, d, lambda*x, err); + return clprobdistGammaComplCDF_1(alpha, d, lambda * x, err); } -cl_double clprobdistGammaCDF_1(cl_double alpha, int d, cl_double x, clprobdistStatus* err){ +cl_double clprobdistGammaCDF_1(cl_double alpha, int d, cl_double x, clprobdistStatus* err) { cl_double ret = clprobdistGammaCDF_2(alpha, d, x, err); if (ret < 0.0) return 1.0 - clprobdistGammaComplCDF_1(alpha, d, x, err); return ret; } -cl_double clprobdistGammaCDF_2(cl_double alpha, int d, cl_double x, clprobdistStatus* err){ +cl_double clprobdistGammaCDF_2(cl_double alpha, int d, cl_double x, clprobdistStatus* err) { if (err) *err = CLPROBDIST_SUCCESS; if (x <= 0.0) @@ -541,8 +538,7 @@ cl_double clprobdistGammaCDF_2(cl_double alpha, int d, cl_double x, clprobdistSt if (alpha > 10.0) { if (x > alpha * 10.0) return 1.0; - } - else { + } else { if (x > clprobdistGammaXBIG) return 1.0; } @@ -557,7 +553,7 @@ cl_double clprobdistGammaCDF_2(cl_double alpha, int d, cl_double x, clprobdistSt if (x <= 1.0 || x < alpha) { cl_double factor, z, rn, term; - factor = exp(alpha*log(x) - x - lgamma(alpha)); + factor = exp(alpha * log(x) - x - lgamma(alpha)); cl_double EPS = clprobdistGammaEPSARRAY[d]; z = 1.0; term = 1.0; @@ -567,7 +563,7 @@ cl_double clprobdistGammaCDF_2(cl_double alpha, int d, cl_double x, clprobdistSt term *= x / rn; z += term; } while (term >= EPS * z); - return z*factor / alpha; + return z * factor / alpha; } // this is a workaround to avoid birecusivity on GPU's @@ -576,23 +572,20 @@ cl_double clprobdistGammaCDF_2(cl_double alpha, int d, cl_double x, clprobdistSt } // Dynamic functions -cl_double clprobdistGammaCDF(cl_double alpha, cl_double lambda, int d, cl_double x, clprobdistStatus* err){ - return clprobdistGammaCDF_1(alpha, d, lambda*x, err); +cl_double clprobdistGammaCDF(cl_double alpha, cl_double lambda, int d, cl_double x, clprobdistStatus* err) { + return clprobdistGammaCDF_1(alpha, d, lambda * x, err); } -cl_double clprobdistGammaCDFWithObject(_CLPROBDIST_GAMMA_OBJ_MEM const clprobdistGamma* distObj, cl_double x, clprobdistStatus* err) -{ - return clprobdistGammaCDF(distObj->params.alpha,distObj->params.lambda,distObj->contDistObj.decPrec,x,err); +cl_double clprobdistGammaCDFWithObject(_CLPROBDIST_GAMMA_OBJ_MEM const clprobdistGamma* distObj, cl_double x, clprobdistStatus* err) { + return clprobdistGammaCDF(distObj->params.alpha, distObj->params.lambda, distObj->contDistObj.decPrec, x, err); } -cl_double clprobdistGammaInverseCDF(cl_double alpha, cl_double lambda, int d, cl_double u, clprobdistStatus* err) -{ +cl_double clprobdistGammaInverseCDF(cl_double alpha, cl_double lambda, int d, cl_double u, clprobdistStatus* err) { return clprobdistGammaInverseCDF_1(alpha, d, u, err) / lambda; } -cl_double clprobdistGammaInverseCDFWithObject(_CLPROBDIST_GAMMA_OBJ_MEM const clprobdistGamma* distObj, cl_double u, clprobdistStatus* err) -{ +cl_double clprobdistGammaInverseCDFWithObject(_CLPROBDIST_GAMMA_OBJ_MEM const clprobdistGamma* distObj, cl_double u, clprobdistStatus* err) { return clprobdistGammaInverseCDF_1(distObj->params.alpha, distObj->contDistObj.decPrec, u, err) / distObj->params.lambda; } @@ -629,9 +622,13 @@ cl_double clprobdistGammaGetAlpha(_CLPROBDIST_GAMMA_OBJ_MEM const clprobdistGamm return distObj->params.alpha; } cl_double clprobdistGammaGetLambda(_CLPROBDIST_GAMMA_OBJ_MEM const clprobdistGamma* distObj, clprobdistStatus* err) { - + if (err) *err = CLPROBDIST_SUCCESS; return distObj->params.lambda; } +#ifdef __cplusplus +} +#endif + #endif diff --git a/src/include/clProbDist/private/lognormal.c.h b/src/include/clProbDist/private/lognormal.c.h index ed7b442..9ca3ff0 100644 --- a/src/include/clProbDist/private/lognormal.c.h +++ b/src/include/clProbDist/private/lognormal.c.h @@ -30,6 +30,11 @@ #define _CLPROBDIST_LOGNORMAL_OBJ_MEM #endif +#ifdef __cplusplus +extern "C" +{ +#endif + struct _clprobdistlLognormal { clprobdistContinuous contDistObj; cl_double mu; @@ -47,7 +52,7 @@ constant const cl_double clprobdistLognormalXBIGM = 1000.0; //*********************************************************************** cl_double clprobdistLognormalDensity(cl_double mu, cl_double sigma, cl_double x, clprobdistStatus* err) { - if (sigma <= 0.0){ + if (sigma <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): sigma <= 0", __func__); return 0; } @@ -56,11 +61,11 @@ cl_double clprobdistLognormalDensity(cl_double mu, cl_double sigma, cl_double x, cl_double diff = log(x) - mu; - return exp(-diff*diff / (2 * sigma*sigma)) / (sqrt(2 * clprobdistLognormalPI)*sigma*x); + return exp(-diff * diff / (2 * sigma * sigma)) / (sqrt(2 * clprobdistLognormalPI) * sigma * x); } cl_double clprobdistLognormalCDF(cl_double mu, cl_double sigma, cl_double x, clprobdistStatus* err) { - if (sigma <= 0.0){ + if (sigma <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): sigma <= 0", __func__); return 0; } @@ -70,7 +75,7 @@ cl_double clprobdistLognormalCDF(cl_double mu, cl_double sigma, cl_double x, clp return clprobdistStdNormalCDF((log(x) - mu) / sigma, err); } cl_double clprobdistLognormalComplCDF(cl_double mu, cl_double sigma, cl_double x, clprobdistStatus* err) { - if (sigma <= 0.0){ + if (sigma <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): sigma <= 0", __func__); return 0; } @@ -81,12 +86,12 @@ cl_double clprobdistLognormalComplCDF(cl_double mu, cl_double sigma, cl_double x cl_double clprobdistLognormalInverseCDF(cl_double mu, cl_double sigma, cl_double u, clprobdistStatus* err) { cl_double t, v; - if (sigma <= 0.0){ + if (sigma <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): sigma <= 0", __func__); return 0; } - if (u > 1.0 || u < 0.0){ + if (u > 1.0 || u < 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): u is not in [0, 1]", __func__); return -1; } @@ -102,13 +107,13 @@ cl_double clprobdistLognormalInverseCDF(cl_double mu, cl_double sigma, cl_double if ((t >= clprobdistLognormalXBIG) || (v >= DBL_MAX_EXP * clprobdistLognormalLN2)) return DBL_MAX; - if ((t <= -clprobdistLognormalXBIG) || (v <= -DBL_MAX_EXP* clprobdistLognormalLN2)) + if ((t <= -clprobdistLognormalXBIG) || (v <= -DBL_MAX_EXP * clprobdistLognormalLN2)) return 0.0; return exp(v); } cl_double clprobdistLognormalMean(cl_double mu, cl_double sigma, clprobdistStatus* err) { - if (sigma <= 0.0){ + if (sigma <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): sigma <= 0", __func__); return 0; } @@ -116,7 +121,7 @@ cl_double clprobdistLognormalMean(cl_double mu, cl_double sigma, clprobdistStatu return (exp(mu + (sigma * sigma) / 2.0)); } cl_double clprobdistLognormalVariance(cl_double mu, cl_double sigma, clprobdistStatus* err) { - if (sigma <= 0.0){ + if (sigma <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): sigma <= 0", __func__); return 0; } @@ -124,7 +129,7 @@ cl_double clprobdistLognormalVariance(cl_double mu, cl_double sigma, clprobdistS return (exp(2.0 * mu + sigma * sigma) * (exp(sigma * sigma) - 1.0)); } cl_double clprobdistLognormalStdDeviation(cl_double mu, cl_double sigma, clprobdistStatus* err) { - + return sqrt(clprobdistLognormalVariance(mu, sigma, err)); } @@ -169,4 +174,9 @@ cl_double clprobdistLognormalGetSigma(_CLPROBDIST_LOGNORMAL_OBJ_MEM const clprob if (err) *err = CLPROBDIST_SUCCESS; return dist->sigma; } + +#ifdef __cplusplus +} +#endif + #endif diff --git a/src/include/clProbDist/private/normal.c.h b/src/include/clProbDist/private/normal.c.h index b816c34..0da47f3 100644 --- a/src/include/clProbDist/private/normal.c.h +++ b/src/include/clProbDist/private/normal.c.h @@ -30,6 +30,11 @@ #define _CLPROBDIST_NORMAL_OBJ_MEM #endif +#ifdef __cplusplus +extern "C" +{ +#endif + struct _clprobdistNormal { clprobdistContinuous contDistObj; cl_double mu; @@ -79,19 +84,19 @@ constant cl_double clprobdistNormal_AbarF[] = { }; static cl_double clprobdistNormalEvalCheby(constant cl_double* a, int n, cl_double x, clprobdistStatus* err) { - if (fabs(x) > 1.0){ + if (fabs(x) > 1.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): Chebychev polynomial evaluated at x outside [-1, 1]", __func__); return 0; } - cl_double xx = 2.0*x; + cl_double xx = 2.0 * x; cl_double b0 = 0.0; cl_double b1 = 0.0; cl_double b2 = 0.0; for (int j = n; j >= 0; j--) { b2 = b1; b1 = b0; - b0 = (xx*b1 - b2) + a[j]; + b0 = (xx * b1 - b2) + a[j]; } if (err) *err = CLPROBDIST_SUCCESS; @@ -103,17 +108,16 @@ static cl_double clprobdistNormalEvalCheby(constant cl_double* a, int n, cl_doub //*********************************************************************** cl_double clprobdistStdNormalDensity(cl_double x, clprobdistStatus* err) { if (err) *err = CLPROBDIST_SUCCESS; - return exp(-0.5*x*x) / clprobdistNormal_SQRT2PI; + return exp(-0.5 * x * x) / clprobdistNormal_SQRT2PI; } cl_double clprobdistNormalDensity(cl_double mu, cl_double sigma, cl_double x, clprobdistStatus* err) { - if (sigma <= 0.0){ + if (sigma <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): sigma <= 0", __func__); return 0; - } - else{ + } else { cl_double z = (x - mu) / sigma; - return exp(-0.5*z*z) / (clprobdistNormal_SQRT2PI*sigma); + return exp(-0.5 * z * z) / (clprobdistNormal_SQRT2PI * sigma); } } @@ -185,8 +189,7 @@ cl_double clprobdistStdNormalCDF(cl_double x, clprobdistStatus* err) { x = -x; t = (x - 3.75) / (x + 3.75); r = 1.0 - 0.5 * exp(-x * x) * clprobdistNormalEvalCheby(clprobdistNormal_NORMAL2_A, clprobdistNormal_COEFFMAX, t, err); - } - else { + } else { t = (x - 3.75) / (x + 3.75); r = 0.5 * exp(-x * x) * clprobdistNormalEvalCheby(clprobdistNormal_NORMAL2_A, clprobdistNormal_COEFFMAX, t, err); } @@ -194,11 +197,10 @@ cl_double clprobdistStdNormalCDF(cl_double x, clprobdistStatus* err) { } cl_double clprobdistNormalCDF(cl_double mu, cl_double sigma, cl_double x, clprobdistStatus* err) { - if (sigma <= 0.0){ + if (sigma <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): sigma <= 0", __func__); return 0; - } - else + } else return clprobdistStdNormalCDF((x - mu) / sigma, err); } @@ -238,11 +240,10 @@ cl_double clprobdistStdNormalComplCDF(cl_double x, clprobdistStatus* err) { } cl_double clprobdistNormalComplCDF(cl_double mu, cl_double sigma, cl_double x, clprobdistStatus* err) { - if (sigma <= 0.0){ + if (sigma <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): sigma <= 0", __func__); return 0; - } - else + } else return clprobdistStdNormalComplCDF((x - mu) / sigma, err); } @@ -328,8 +329,7 @@ cl_double clprobdistStdNormalInverseCDF(cl_double u, clprobdistStatus* err) { cl_double y, z, v, w; cl_double x = u; - if (u < 0.0 || u > 1.0) - { + if (u < 0.0 || u > 1.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): u is not in [0, 1]", __func__); return -1; } @@ -343,8 +343,7 @@ cl_double clprobdistStdNormalInverseCDF(cl_double u, clprobdistStatus* err) { if (x < 0.0) { x = -x; negatif = CL_TRUE; - } - else { + } else { negatif = CL_FALSE; } @@ -357,8 +356,7 @@ cl_double clprobdistStdNormalInverseCDF(cl_double u, clprobdistStatus* err) { } z = (v / w) * x; - } - else if (x <= 0.9375) { + } else if (x <= 0.9375) { y = x * x - 0.87890625; v = w = 0.0; for (i = 7; i >= 0; i--) { @@ -367,8 +365,7 @@ cl_double clprobdistStdNormalInverseCDF(cl_double u, clprobdistStatus* err) { } z = (v / w) * x; - } - else { + } else { if (u > 0.5) y = 1.0 / sqrt(-log(1.0 - x)); else @@ -400,39 +397,35 @@ cl_double clprobdistStdNormalInverseCDF(cl_double u, clprobdistStatus* err) { } return -(z * clprobdistNormal_SQRT2); - } - else + } else return z * clprobdistNormal_SQRT2; } cl_double clprobdistNormalInverseCDF(cl_double mu, cl_double sigma, cl_double u, clprobdistStatus* err) { - if (sigma <= 0.0){ + if (sigma <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): sigma <= 0", __func__); return 0; - } - else - return mu + sigma*clprobdistStdNormalInverseCDF(u, err); + } else + return mu + sigma * clprobdistStdNormalInverseCDF(u, err); } cl_double clprobdistNormalMean(cl_double mu, cl_double sigma, clprobdistStatus* err) { - if (sigma <= 0.0){ + if (sigma <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): sigma <= 0", __func__); return 0; - } - else{ + } else { if (err) *err = CLPROBDIST_SUCCESS; return mu; } } cl_double clprobdistNormalVariance(cl_double mu, cl_double sigma, clprobdistStatus* err) { - if (sigma <= 0.0){ + if (sigma <= 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): sigma <= 0", __func__); return 0; - } - else{ + } else { if (err) *err = CLPROBDIST_SUCCESS; return (sigma * sigma); } @@ -450,7 +443,7 @@ cl_double clprobdistNormalStdDeviation(cl_double mu, cl_double sigma, clprobdist cl_double clprobdistNormalDensityWithObject(_CLPROBDIST_NORMAL_OBJ_MEM const clprobdistNormal* dist, cl_double x, clprobdistStatus* err) { if (err) *err = CLPROBDIST_SUCCESS; cl_double z = (x - dist->mu) / dist->sigma; - return exp(-0.5*z*z) / (clprobdistNormal_SQRT2PI*dist->sigma); + return exp(-0.5 * z * z) / (clprobdistNormal_SQRT2PI * dist->sigma); } cl_double clprobdistNormalCDFWithObject(_CLPROBDIST_NORMAL_OBJ_MEM const clprobdistNormal* dist, cl_double x, clprobdistStatus* err) { @@ -486,4 +479,9 @@ cl_double clprobdistNormalGetSigma(_CLPROBDIST_NORMAL_OBJ_MEM const clprobdistNo if (err) *err = CLPROBDIST_SUCCESS; return dist->sigma; } + +#ifdef __cplusplus +} +#endif + #endif diff --git a/src/include/clProbDist/private/poisson.c.h b/src/include/clProbDist/private/poisson.c.h index 4a1412e..e31578b 100644 --- a/src/include/clProbDist/private/poisson.c.h +++ b/src/include/clProbDist/private/poisson.c.h @@ -34,6 +34,11 @@ #define _CLPROBDIST_POISSON_OBJ_MEM #endif +#ifdef __cplusplus +extern "C" +{ +#endif + #ifndef __OPENCL_C_VERSION__ #include #endif @@ -63,17 +68,17 @@ constant const cl_double clprobdistPoisson_EPSILON = 1.0e-16; cl_double clprobdistPoissonCDF_1(cl_double lambda, cl_int x, clprobdistStatus* err); -static cl_double clprobdistPoisson_CDF(cl_int x, _CLPROBDIST_POISSON_OBJ_MEM const clprobdistPoisson* dist){ +static cl_double clprobdistPoisson_CDF(cl_int x, _CLPROBDIST_POISSON_OBJ_MEM const clprobdistPoisson* dist) { return dist->cdf[x + (dist->params.len - 1)]; } static cl_double clprobdistPoisson_factorial(cl_int n, clprobdistStatus* err) { if (err) *err = CLPROBDIST_SUCCESS; - if (n < 0){ + if (n < 0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): n < 0", __func__); return -1; } - + cl_double T = 1; for (cl_int j = 2; j <= n; j++) T *= j; @@ -82,28 +87,26 @@ static cl_double clprobdistPoisson_factorial(cl_int n, clprobdistStatus* err) { cl_double clprobdistPoissonProb(cl_double lambda, cl_int x, clprobdistStatus* err) { if (err) *err = CLPROBDIST_SUCCESS; - if (x < 0) + if (x < 0) return 0.0; if (lambda >= 100.0) { - if ((cl_double)x >= 10.0*lambda) + if ((cl_double)x >= 10.0 * lambda) return 0.0; - } - else if (lambda >= 3.0) { - if ((cl_double)x >= 100.0*lambda) - return 0.0; - } - else { - if ((cl_double)x >= 200.0*fmax(1.0, lambda)) + } else if (lambda >= 3.0) { + if ((cl_double)x >= 100.0 * lambda) + return 0.0; + } else { + if ((cl_double)x >= 200.0 * fmax(1.0, lambda)) return 0.0; } cl_double lambdaLIM = 20.0; cl_double Res; if (lambda < lambdaLIM && x <= 100) - Res = exp(-lambda)*pow(lambda, x) / clprobdistPoisson_factorial(x, err); + Res = exp(-lambda) * pow(lambda, x) / clprobdistPoisson_factorial(x, err); else { - cl_double y = x*log(lambda) - lgamma(x + 1.0) - lambda; + cl_double y = x * log(lambda) - lgamma(x + 1.0) - lambda; Res = exp(y); } return Res; @@ -114,21 +117,20 @@ cl_double clprobdistPoissonCDF(cl_double lambda, cl_int x, clprobdistStatus* err cl_double clprobdistPoissonComplCDF(cl_double lambda, cl_int x, clprobdistStatus* err) { if (err) *err = CLPROBDIST_SUCCESS; //check Params - if (lambda < 0){ + if (lambda < 0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): lambda < 0", __func__); return -1; } - if (x <= 0) + if (x <= 0) return 1.0; if (lambda >= 100.0) { - if ((cl_double)x >= 10.0*lambda) + if ((cl_double)x >= 10.0 * lambda) return 0.0; - } - else { + } else { if ((cl_double)x >= 100 + 100.0 * fmax(1.0, lambda)) - return 0.0; + return 0.0; } /* If lambda > lambdaLIM, we use the Chi2 distribution according to the @@ -145,7 +147,7 @@ cl_double clprobdistPoissonComplCDF(cl_double lambda, cl_int x, clprobdistStatus if (x <= lambda) return 1.0 - clprobdistPoissonCDF_1(lambda, x - 1, err); - + // Naive computation: sum all prob. from i = x to i = oo cl_double term, sum; @@ -153,7 +155,7 @@ cl_double clprobdistPoissonComplCDF(cl_double lambda, cl_int x, clprobdistStatus // Sum at least IMAX prob. terms from i = s to i = oo sum = term = clprobdistPoissonProb(lambda, x, err); - + cl_int i = x + 1; while (term > clprobdistPoisson_EPSILON || i <= x + IMAX) { term *= lambda / i; @@ -177,24 +179,23 @@ cl_double clprobdistPoissonCDF_1(cl_double lambda, cl_int x, clprobdistStatus* e * naive computation for dist->params.lambdalim > 200.0, slower for dist->params.lambdalim < 200.0 */ if (err) *err = CLPROBDIST_SUCCESS; - if (lambda < 0.0){ + if (lambda < 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): lambda < 0", __func__); return -1; } if (lambda == 0.0) return 1.0; - + if (x < 0) return 0.0; - + if (lambda >= 100.0) { - if ((cl_double)x >= 10.0*lambda) - return 1.0; - } - else { - if ((cl_double)x >= 100.0*fmax(1.0, lambda)) + if ((cl_double)x >= 10.0 * lambda) + return 1.0; + } else { + if ((cl_double)x >= 100.0 * fmax(1.0, lambda)) return 1.0; } @@ -221,17 +222,17 @@ cl_double clprobdistPoissonCDF_1(cl_double lambda, cl_int x, clprobdistStatus* e term *= lambda / j; sum += term; } - return sum*exp(-lambda); + return sum * exp(-lambda); } cl_int clprobdistPoissonInverseCDF(cl_double lambda, cl_double u, clprobdistStatus* err) { //Check Params - if (u < 0.0 || u > 1.0){ + if (u < 0.0 || u > 1.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): u is not in range [0,1]", __func__); return -1; } - if (lambda < 0.0){ + if (lambda < 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): lambda < 0", __func__); return -1; } @@ -239,15 +240,14 @@ cl_int clprobdistPoissonInverseCDF(cl_double lambda, cl_double u, clprobdistStat if (u >= 1.0) return INT_MAX; - if (u <= clprobdistPoissonProb(lambda, 0, err)) + if (u <= clprobdistPoissonProb(lambda, 0, err)) return 0; //Calculate InverseCDF cl_int i; cl_double lambdaLIM = 700.0; - if (lambda < lambdaLIM) - { + if (lambda < lambdaLIM) { cl_double sumprev = -1.0; cl_double term = exp(-lambda); cl_double sum = term; @@ -261,8 +261,7 @@ cl_int clprobdistPoissonInverseCDF(cl_double lambda, cl_double u, clprobdistStat return i; - } - else { + } else { // When lambda is very large, the probabilities are empirically // negligible when i is far from lambda. We start with a binary search // over [0,lambda] for the smallest value of i for which the probability @@ -296,7 +295,7 @@ cl_int clprobdistPoissonInverseCDF(cl_double lambda, cl_double u, clprobdistStat // probabilities for i <= mid. // Probabilities below clprobdistPoisson_EPSILON*u are deemed // negligible. - while (term >= clprobdistPoisson_EPSILON*u && i > 0) { + while (term >= clprobdistPoisson_EPSILON * u && i > 0) { term *= i / lambda; sum += term; i--; @@ -313,8 +312,7 @@ cl_int clprobdistPoissonInverseCDF(cl_double lambda, cl_double u, clprobdistStat prev = sum; sum += term; } - } - else { + } else { // The computed CDF is too big so we substract from it. sum -= term; while (sum >= u) { @@ -332,7 +330,7 @@ cl_double clprobdistPoissonMean(cl_double lambda, clprobdistStatus* err) { if (err) *err = CLPROBDIST_SUCCESS; //Check Params - if (lambda < 0.0){ + if (lambda < 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): lambda < 0", __func__); return -1; } @@ -344,7 +342,7 @@ cl_double clprobdistPoissonVariance(cl_double lambda, clprobdistStatus* err) { if (err) *err = CLPROBDIST_SUCCESS; //Check Params - if (lambda < 0.0){ + if (lambda < 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): lambda < 0", __func__); return -1; } @@ -356,7 +354,7 @@ cl_double clprobdistPoissonStdDeviation(cl_double lambda, clprobdistStatus* err) if (err) *err = CLPROBDIST_SUCCESS; //Check Params - if (lambda < 0.0){ + if (lambda < 0.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): lambda < 0", __func__); return -1; } @@ -369,7 +367,7 @@ cl_double clprobdistPoissonStdDeviation(cl_double lambda, clprobdistStatus* err) cl_double clprobdistPoissonProbWithObject(_CLPROBDIST_POISSON_OBJ_MEM const clprobdistPoisson* dist, cl_int x, clprobdistStatus* err) { if (err) *err = CLPROBDIST_SUCCESS; - if (x < 0) + if (x < 0) return 0.0; if (!dist->pdf) @@ -387,7 +385,7 @@ cl_double clprobdistPoissonCDFWithObject(_CLPROBDIST_POISSON_OBJ_MEM const clpro if (x < 0) return 0.0; - + if (dist->params.lambda == 0.0) return 1.0; @@ -397,7 +395,7 @@ cl_double clprobdistPoissonCDFWithObject(_CLPROBDIST_POISSON_OBJ_MEM const clpro cdf (lambda, x) = 1 - chiSquare (2x + 2, 2*lambda) which equals also 1 - gamma (x + 1, dist->params.lambda) */ - + if (!dist->cdf) return clprobdistGammaComplCDF_1(x + 1.0, 15, dist->params.lambda, err); @@ -409,7 +407,7 @@ cl_double clprobdistPoissonCDFWithObject(_CLPROBDIST_POISSON_OBJ_MEM const clpro // could also call clprobdistGamma.barF instead. cl_int RMAX = 20; cl_int i; - cl_double term; + cl_double term; Sum = term = clprobdistPoissonProb(dist->params.lambda, x, err); i = x; while (i > 0 && i >= x - RMAX) { @@ -433,7 +431,7 @@ cl_double clprobdistPoissonComplCDFWithObject(_CLPROBDIST_POISSON_OBJ_MEM const */ if (err) *err = CLPROBDIST_SUCCESS; - if (x <= 0) + if (x <= 0) return 1.0; /* For large dist->params.lambda, we use the Chi2 distribution according to the exact @@ -452,30 +450,29 @@ cl_double clprobdistPoissonComplCDFWithObject(_CLPROBDIST_POISSON_OBJ_MEM const if (x <= dist->params.xmin) return 1.0; - if (x > dist->params.xmed){ + if (x > dist->params.xmed) { // We keep the complementary distribution in the upper part of cdf return clprobdistPoisson_CDF(x - dist->params.xmin, dist); - } - else + } else return 1.0 - clprobdistPoisson_CDF(x - 1 - dist->params.xmin, dist); } cl_int clprobdistDiscreteInverseCDF(_CLPROBDIST_POISSON_OBJ_MEM const clprobdistPoisson* dist, cl_double u, clprobdistStatus* err) { if (err) *err = CLPROBDIST_SUCCESS; - cl_int i=0, j=0, k=0; + cl_int i = 0, j = 0, k = 0; - if (u < 0.0 || u > 1.0){ + if (u < 0.0 || u > 1.0) { if (err) *err = clprobdistSetErrorString(CLPROBDIST_INVALID_VALUE, "%s(): u is not in [0,1]", __func__); return -1; } - + if (u <= 0.0) return dist->params.supportA; - + if (u >= 1.0) return dist->params.supportB; - + // Remember: the upper part of cdf contains the complementary distribu- // tion for xmed < s <= xmax, and the lower part of cdf the @@ -485,7 +482,7 @@ cl_int clprobdistDiscreteInverseCDF(_CLPROBDIST_POISSON_OBJ_MEM const clprobdist // In the lower part of cdf if (u <= clprobdistPoisson_CDF(0, dist)) return dist->params.xmin; - + i = 0; j = dist->params.xmed - dist->params.xmin; while (i < j) { @@ -495,24 +492,23 @@ cl_int clprobdistDiscreteInverseCDF(_CLPROBDIST_POISSON_OBJ_MEM const clprobdist else j = k; } - } - else { + } else { // In the upper part of cdf u = 1 - u; if (u < clprobdistPoisson_CDF(dist->params.xmax - dist->params.xmin, dist)) return dist->params.xmax; - i = dist->params.xmed - dist->params.xmin + 1; - j = dist->params.xmax - dist->params.xmin; - while (i < j) { - k = (i + j) / 2; - if (u < clprobdistPoisson_CDF(k, dist)) - i = k + 1; - else - j = k; - } - i--; - + i = dist->params.xmed - dist->params.xmin + 1; + j = dist->params.xmax - dist->params.xmin; + while (i < j) { + k = (i + j) / 2; + if (u < clprobdistPoisson_CDF(k, dist)) + i = k + 1; + else + j = k; + } + i--; + } return i + dist->params.xmin; } @@ -537,11 +533,13 @@ cl_double clprobdistPoissonStdDeviationWithObject(_CLPROBDIST_POISSON_OBJ_MEM co return clprobdistPoissonStdDeviation(dist->params.lambda, err); } -cl_double clprobdistPoissonGetLambda(_CLPROBDIST_POISSON_OBJ_MEM const clprobdistPoisson* dist, clprobdistStatus* err) -{ +cl_double clprobdistPoissonGetLambda(_CLPROBDIST_POISSON_OBJ_MEM const clprobdistPoisson* dist, clprobdistStatus* err) { if (err) *err = CLPROBDIST_SUCCESS; return dist->params.lambda; } +#ifdef __cplusplus +} +#endif #endif // PRIVATE_POISSONDIST_CH