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