Line | Branch | Exec | Source |
---|---|---|---|
1 | /* | ||
2 | * This file is part of FFmpeg. | ||
3 | * | ||
4 | * FFmpeg is free software; you can redistribute it and/or | ||
5 | * modify it under the terms of the GNU Lesser General Public | ||
6 | * License as published by the Free Software Foundation; either | ||
7 | * version 2.1 of the License, or (at your option) any later version. | ||
8 | * | ||
9 | * FFmpeg is distributed in the hope that it will be useful, | ||
10 | * but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
11 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | ||
12 | * Lesser General Public License for more details. | ||
13 | * | ||
14 | * You should have received a copy of the GNU Lesser General Public | ||
15 | * License along with FFmpeg; if not, write to the Free Software | ||
16 | * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA | ||
17 | */ | ||
18 | |||
19 | #include "libavutil/log.h" | ||
20 | #include "libavutil/mem.h" | ||
21 | #include "libavutil/timer.h" | ||
22 | #include "libavutil/lfg.h" | ||
23 | |||
24 | static const double Z_TABLE[31][10] = { | ||
25 | {0.5000, 0.5040, 0.5080, 0.5120, 0.5160, 0.5199, 0.5239, 0.5279, 0.5319, 0.5359}, | ||
26 | {0.5398, 0.5438, 0.5478, 0.5517, 0.5557, 0.5596, 0.5636, 0.5675, 0.5714, 0.5753}, | ||
27 | {0.5793, 0.5832, 0.5871, 0.5910, 0.5948, 0.5987, 0.6026, 0.6064, 0.6103, 0.6141}, | ||
28 | {0.6179, 0.6217, 0.6255, 0.6293, 0.6331, 0.6368, 0.6406, 0.6443, 0.6480, 0.6517}, | ||
29 | {0.6554, 0.6591, 0.6628, 0.6664, 0.6700, 0.6736, 0.6772, 0.6808, 0.6844, 0.6879}, | ||
30 | {0.6915, 0.6950, 0.6985, 0.7019, 0.7054, 0.7088, 0.7123, 0.7157, 0.7190, 0.7224}, | ||
31 | {0.7257, 0.7291, 0.7324, 0.7357, 0.7389, 0.7422, 0.7454, 0.7486, 0.7517, 0.7549}, | ||
32 | {0.7580, 0.7611, 0.7642, 0.7673, 0.7704, 0.7734, 0.7764, 0.7794, 0.7823, 0.7852}, | ||
33 | {0.7881, 0.7910, 0.7939, 0.7967, 0.7995, 0.8023, 0.8051, 0.8078, 0.8106, 0.8133}, | ||
34 | {0.8159, 0.8186, 0.8212, 0.8238, 0.8264, 0.8289, 0.8315, 0.8340, 0.8365, 0.8389}, | ||
35 | {0.8413, 0.8438, 0.8461, 0.8485, 0.8508, 0.8531, 0.8554, 0.8577, 0.8599, 0.8621}, | ||
36 | {0.8643, 0.8665, 0.8686, 0.8708, 0.8729, 0.8749, 0.8770, 0.8790, 0.8810, 0.8830}, | ||
37 | {0.8849, 0.8869, 0.8888, 0.8907, 0.8925, 0.8944, 0.8962, 0.8980, 0.8997, 0.9015}, | ||
38 | {0.9032, 0.9049, 0.9066, 0.9082, 0.9099, 0.9115, 0.9131, 0.9147, 0.9162, 0.9177}, | ||
39 | {0.9192, 0.9207, 0.9222, 0.9236, 0.9251, 0.9265, 0.9279, 0.9292, 0.9306, 0.9319}, | ||
40 | {0.9332, 0.9345, 0.9357, 0.9370, 0.9382, 0.9394, 0.9406, 0.9418, 0.9429, 0.9441}, | ||
41 | {0.9452, 0.9463, 0.9474, 0.9484, 0.9495, 0.9505, 0.9515, 0.9525, 0.9535, 0.9545}, | ||
42 | {0.9554, 0.9564, 0.9573, 0.9582, 0.9591, 0.9599, 0.9608, 0.9616, 0.9625, 0.9633}, | ||
43 | {0.9641, 0.9649, 0.9656, 0.9664, 0.9671, 0.9678, 0.9686, 0.9693, 0.9699, 0.9706}, | ||
44 | {0.9713, 0.9719, 0.9726, 0.9732, 0.9738, 0.9744, 0.9750, 0.9756, 0.9761, 0.9767}, | ||
45 | {0.9772, 0.9778, 0.9783, 0.9788, 0.9793, 0.9798, 0.9803, 0.9808, 0.9812, 0.9817}, | ||
46 | {0.9821, 0.9826, 0.9830, 0.9834, 0.9838, 0.9842, 0.9846, 0.9850, 0.9854, 0.9857}, | ||
47 | {0.9861, 0.9864, 0.9868, 0.9871, 0.9875, 0.9878, 0.9881, 0.9884, 0.9887, 0.9890}, | ||
48 | {0.9893, 0.9896, 0.9898, 0.9901, 0.9904, 0.9906, 0.9909, 0.9911, 0.9913, 0.9916}, | ||
49 | {0.9918, 0.9920, 0.9922, 0.9925, 0.9927, 0.9929, 0.9931, 0.9932, 0.9934, 0.9936}, | ||
50 | {0.9938, 0.9940, 0.9941, 0.9943, 0.9945, 0.9946, 0.9948, 0.9949, 0.9951, 0.9952}, | ||
51 | {0.9953, 0.9955, 0.9956, 0.9957, 0.9959, 0.9960, 0.9961, 0.9962, 0.9963, 0.9964}, | ||
52 | {0.9965, 0.9966, 0.9967, 0.9968, 0.9969, 0.9970, 0.9971, 0.9972, 0.9973, 0.9974}, | ||
53 | {0.9974, 0.9975, 0.9976, 0.9977, 0.9977, 0.9978, 0.9979, 0.9979, 0.9980, 0.9981}, | ||
54 | {0.9981, 0.9982, 0.9982, 0.9983, 0.9984, 0.9984, 0.9985, 0.9985, 0.9986, 0.9986}, | ||
55 | {0.9987, 0.9987, 0.9987, 0.9988, 0.9988, 0.9989, 0.9989, 0.9989, 0.9990, 0.9990} }; | ||
56 | |||
57 | // Inverse cumulative distribution function | ||
58 | 1998 | static double inv_cdf(double u) | |
59 | { | ||
60 | const double a[4] = { 2.50662823884, | ||
61 | -18.61500062529, | ||
62 | 41.39119773534, | ||
63 | -25.44106049637}; | ||
64 | |||
65 | const double b[4] = {-8.47351093090, | ||
66 | 23.08336743743, | ||
67 | -21.06224101826, | ||
68 | 3.13082909833}; | ||
69 | |||
70 | const double c[9] = {0.3374754822726147, | ||
71 | 0.9761690190917186, | ||
72 | 0.1607979714918209, | ||
73 | 0.0276438810333863, | ||
74 | 0.0038405729373609, | ||
75 | 0.0003951896511919, | ||
76 | 0.0000321767881768, | ||
77 | 0.0000002888167364, | ||
78 | 0.0000003960315187}; | ||
79 | |||
80 | double r; | ||
81 | 1998 | double x = u - 0.5; | |
82 | |||
83 | // Beasley-Springer | ||
84 |
2/2✓ Branch 0 taken 1680 times.
✓ Branch 1 taken 318 times.
|
1998 | if (fabs(x) < 0.42) { |
85 | |||
86 | 1680 | double y = x * x; | |
87 | 1680 | r = x * (((a[3]*y+a[2])*y+a[1])*y+a[0]) / | |
88 | 1680 | ((((b[3]*y+b[2])*y+b[1])*y+b[0])*y+1.0); | |
89 | } | ||
90 | else {// Moro | ||
91 | 318 | r = u; | |
92 |
2/2✓ Branch 0 taken 159 times.
✓ Branch 1 taken 159 times.
|
318 | if (x > 0.0) |
93 | 159 | r = 1.0 - u; | |
94 | 318 | r = log(-log(r)); | |
95 | 318 | r = c[0] + r*(c[1]+r*(c[2]+r*(c[3]+r*(c[4]+r*(c[5]+r*(c[6]+ | |
96 | 318 | r*(c[7]+r*c[8]))))))); | |
97 |
2/2✓ Branch 0 taken 159 times.
✓ Branch 1 taken 159 times.
|
318 | if (x < 0.0) |
98 | 159 | r = -r; | |
99 | } | ||
100 | |||
101 | 1998 | return r; | |
102 | } | ||
103 | 1 | int main(void) | |
104 | { | ||
105 | 1 | int x = 0; | |
106 | int i, j; | ||
107 | AVLFG state; | ||
108 | 1 | av_lfg_init(&state, 0xdeadbeef); | |
109 |
2/2✓ Branch 0 taken 10000 times.
✓ Branch 1 taken 1 times.
|
10001 | for (j = 0; j < 10000; j++) { |
110 |
2/2✓ Branch 0 taken 6240000 times.
✓ Branch 1 taken 10000 times.
|
6250000 | for (i = 0; i < 624; i++) { |
111 | //av_log(NULL, AV_LOG_ERROR, "%X\n", av_lfg_get(&state)); | ||
112 | 6240000 | x += av_lfg_get(&state); | |
113 | } | ||
114 | } | ||
115 | 1 | av_log(NULL, AV_LOG_ERROR, "final value:%X\n", x); | |
116 | |||
117 | /* BMG usage example */ | ||
118 | { | ||
119 | 1 | double mean = 1000; | |
120 | 1 | double stddev = 53; | |
121 | 1 | double samp_mean = 0.0, samp_stddev = 0.0, QH = 0; | |
122 | 1 | double Z, p_value = -1, tot_samp = 1000; | |
123 | 1 | double *PRN_arr = av_malloc_array(tot_samp, sizeof(double)); | |
124 | |||
125 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
1 | if (!PRN_arr) { |
126 | ✗ | fprintf(stderr, "failed to allocate memory!\n"); | |
127 | ✗ | return 1; | |
128 | } | ||
129 | |||
130 | 1 | av_lfg_init(&state, 42); | |
131 |
2/2✓ Branch 0 taken 500 times.
✓ Branch 1 taken 1 times.
|
501 | for (i = 0; i < tot_samp; i += 2) { |
132 | double bmg_out[2]; | ||
133 | 500 | av_bmg_get(&state, bmg_out); | |
134 | 500 | PRN_arr[i ] = bmg_out[0] * stddev + mean; | |
135 | 500 | PRN_arr[i+1] = bmg_out[1] * stddev + mean; | |
136 | 500 | samp_mean += PRN_arr[i] + PRN_arr[i+1]; | |
137 | 500 | samp_stddev += PRN_arr[i] * PRN_arr[i] + PRN_arr[i+1] * PRN_arr[i+1]; | |
138 | 500 | printf("PRN%d : %f\n" | |
139 | "PRN%d : %f\n", | ||
140 | 500 | i, PRN_arr[i], i+1, PRN_arr[i+1]); | |
141 | } | ||
142 | 1 | samp_mean /= tot_samp; | |
143 | 1 | samp_stddev /= (tot_samp - 1); | |
144 | 1 | samp_stddev -= (tot_samp * 1.0 / (tot_samp - 1))*samp_mean*samp_mean; | |
145 | 1 | samp_stddev = sqrt(samp_stddev); | |
146 | 1 | Z = (mean - samp_mean) / (stddev / sqrt(tot_samp)); | |
147 | { | ||
148 | 1 | int x, y, a, b, flag = 0; | |
149 | |||
150 |
1/2✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
|
1 | if (Z < 0.0) { |
151 | 1 | flag = !flag; | |
152 | 1 | Z = Z * -1.0; | |
153 | } | ||
154 | |||
155 | 1 | a = (int)(Z * 100); | |
156 | 1 | b = ((int)Z * 100); | |
157 | 1 | x = Z * 10; | |
158 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
1 | y = (b > 0) ? a % b : a; |
159 | 1 | y = y % 10; | |
160 |
2/4✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
|
1 | if (x > 30 || y > 9) { |
161 | ✗ | av_log(NULL, AV_LOG_INFO, "error: out of bounds! tried to access" | |
162 | "Z_TABLE[%d][%d]\n", x, y); | ||
163 | ✗ | goto SKIP; | |
164 | } | ||
165 |
1/2✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
|
1 | p_value = flag ? 1 - Z_TABLE[x][y] : Z_TABLE[x][y]; |
166 | } | ||
167 | |||
168 |
2/2✓ Branch 0 taken 1000 times.
✓ Branch 1 taken 1 times.
|
1001 | SKIP: for (i = 0; i < tot_samp; ++i) { |
169 | |||
170 |
2/2✓ Branch 0 taken 999 times.
✓ Branch 1 taken 1 times.
|
1000 | if ( i < (tot_samp - 1)) { |
171 | double H_diff; | ||
172 | 999 | H_diff = inv_cdf((i + 2.0 - (3.0/8.0)) / (tot_samp + (1.0/4.0))); | |
173 | 999 | H_diff -= inv_cdf((i + 1.0 - (3.0/8.0)) / (tot_samp + (1.0/4.0))); | |
174 | |||
175 | 999 | QH += ((PRN_arr[i + 1] - PRN_arr[i]) / H_diff); | |
176 | } | ||
177 | } | ||
178 | 1 | QH = 1.0 - QH / ((tot_samp - 1.0) * samp_stddev); | |
179 | |||
180 | 1 | printf("sample mean : %f\n" | |
181 | "true mean : %f\n" | ||
182 | "sample stddev: %f\n" | ||
183 | "true stddev : %f\n" | ||
184 | "z-score : %f\n" | ||
185 | "p-value : %f\n" | ||
186 | "QH[normality]: %f\n", | ||
187 | samp_mean, mean, samp_stddev, stddev, Z, p_value, QH); | ||
188 | |||
189 | 1 | av_freep(&PRN_arr); | |
190 | } | ||
191 | 1 | return 0; | |
192 | } | ||
193 |