1 /* Foundation, miscellanea for the statistics modules. 2 */ 3 #ifndef eslSTATS_INCLUDED 4 #define eslSTATS_INCLUDED 5 #include "esl_config.h" 6 7 #include "easel.h" 8 9 /***************************************************************** 10 * Splitting IEEE754 double-precision float into two uint32_t 11 ***************************************************************** 12 * 13 * Currently we only need these macros for one function, 14 * esl_stats_erfc(). The Sun Microsystems erfc() code that we've 15 * borrowed splits an IEEE754 double into two unsigned 32-bit 16 * integers. It uses arcane trickery to deal with endianness at 17 * runtime, using incantations like these: 18 * n0 = ((*(int*)&one)>>29)^1 0|1 = bigendian | littleendian 19 * hx = *(n0+(int*)&x); get high word 20 * (1-n0+(int*)&z) = 0; set low word 21 * 22 * Not only is this arcane and dubious, static code checking (using 23 * the clang/llvm checker) doesn't like it. I found an improvement 24 * in a library called zenilib, at: 25 * http://www-personal.umich.edu/~bazald/l/api/math__private_8h_source.html 26 * 27 * Here we do the same thing in an ANSI-respecting way using unions, 28 * with endianness detected at compile time. 29 * 30 * The zenilib code also appears to derive from (C) Sun Microsystems 31 * code. 32 * 33 * SRE TODO: insert license/copyright info here 34 */ 35 36 #ifdef WORDS_BIGENDIAN 37 typedef union 38 { 39 double val; 40 struct { 41 uint32_t msw; 42 uint32_t lsw; 43 } parts; 44 } esl_double_split_t; 45 #else /* else we're littleendian, such as Intel */ 46 typedef union 47 { 48 double val; 49 struct { 50 uint32_t lsw; 51 uint32_t msw; 52 } parts; 53 } esl_double_split_t; 54 #endif /*WORDS_BIGENDIAN*/ 55 56 #define ESL_GET_WORDS(ix0, ix1, d) \ 57 do { \ 58 esl_double_split_t esltmp_ds; \ 59 esltmp_ds.val = (d); \ 60 (ix0) = esltmp_ds.parts.msw; \ 61 (ix1) = esltmp_ds.parts.lsw; \ 62 } while (0) 63 64 #define ESL_GET_HIGHWORD(ix0, d) \ 65 do { \ 66 esl_double_split_t esltmp_ds; \ 67 esltmp_ds.val = (d); \ 68 (ix0) = esltmp_ds.parts.msw; \ 69 } while (0) 70 71 #define ESL_GET_LOWWORD(ix0, d) \ 72 do { \ 73 esl_double_split_t esltmp_ds; \ 74 esltmp_ds.val = (d); \ 75 (ix0) = esltmp_ds.parts.lsw; \ 76 } while (0) 77 78 #define ESL_SET_WORDS(d, ix0, ix1) \ 79 do { \ 80 esl_double_split_t esltmp_ds; \ 81 esltmp_ds.parts.msw = (ix0); \ 82 esltmp_ds.parts.lsw = (ix1); \ 83 (d) = esltmp_ds.val; \ 84 } while (0) 85 86 #define ESL_SET_HIGHWORD(d, ix0) \ 87 do { \ 88 esl_double_split_t esltmp_ds; \ 89 esltmp_ds.val = (d); \ 90 esltmp_ds.parts.msw = (ix0); \ 91 (d) = esltmp_ds.val; \ 92 } while (0) 93 94 #define ESL_SET_LOWWORD(d, ix1) \ 95 do { \ 96 esl_double_split_t esltmp_ds; \ 97 esltmp_ds.val = (d); \ 98 esltmp_ds.parts.lsw = (ix1); \ 99 (d) = esltmp_ds.val; \ 100 } while (0) 101 102 103 /***************************************************************** 104 * Function declarations 105 *****************************************************************/ 106 107 /* 1. Summary statistics calculations */ 108 extern int esl_stats_DMean(const double *x, int n, double *opt_mean, double *opt_var); 109 extern int esl_stats_FMean(const float *x, int n, double *opt_mean, double *opt_var); 110 extern int esl_stats_IMean(const int *x, int n, double *opt_mean, double *opt_var); 111 112 /* 2. Special functions */ 113 extern int esl_stats_LogGamma(double x, double *ret_answer); 114 extern int esl_stats_Psi(double x, double *ret_answer); 115 extern int esl_stats_IncompleteGamma(double a, double x, double *ret_pax, double *ret_qax); 116 extern double esl_stats_erfc(double x); 117 118 /* 3. Standard statistical tests */ 119 extern int esl_stats_GTest(int ca, int na, int cb, int nb, double *ret_G, double *ret_P); 120 extern int esl_stats_ChiSquaredTest(int v, double x, double *ret_answer); 121 122 /* 4. Data fitting */ 123 extern int esl_stats_LinearRegression(const double *x, const double *y, const double *sigma, int n, 124 double *opt_a, double *opt_b, 125 double *opt_sigma_a, double *opt_sigma_b, double *opt_cov_ab, 126 double *opt_cc, double *opt_Q); 127 128 129 /***************************************************************** 130 * Portability 131 *****************************************************************/ 132 133 #ifndef HAVE_ERFC 134 #define erfc(x) esl_stats_erfc(x) 135 #endif 136 137 #endif /*eslSTATS_INCLUDED*/ 138 139