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