FFmpeg
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/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 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  double x = u - 0.5;
82 
83  // Beasley-Springer
84  if (fabs(x) < 0.42) {
85 
86  double y = x * x;
87  r = x * (((a[3]*y+a[2])*y+a[1])*y+a[0]) /
88  ((((b[3]*y+b[2])*y+b[1])*y+b[0])*y+1.0);
89  }
90  else {// Moro
91  r = u;
92  if (x > 0.0)
93  r = 1.0 - u;
94  r = log(-log(r));
95  r = c[0] + r*(c[1]+r*(c[2]+r*(c[3]+r*(c[4]+r*(c[5]+r*(c[6]+
96  r*(c[7]+r*c[8])))))));
97  if (x < 0.0)
98  r = -r;
99  }
100 
101  return r;
102 }
103 int main(void)
104 {
105  int x = 0;
106  int i, j;
107  AVLFG state;
108  av_lfg_init(&state, 0xdeadbeef);
109  for (j = 0; j < 10000; j++) {
110  for (i = 0; i < 624; i++) {
111  //av_log(NULL, AV_LOG_ERROR, "%X\n", av_lfg_get(&state));
112  x += av_lfg_get(&state);
113  }
114  }
115  av_log(NULL, AV_LOG_ERROR, "final value:%X\n", x);
116 
117  /* BMG usage example */
118  {
119  double mean = 1000;
120  double stddev = 53;
121  double samp_mean = 0.0, samp_stddev = 0.0, QH = 0;
122  double Z, p_value = -1, tot_samp = 1000;
123  double *PRN_arr = av_malloc_array(tot_samp, sizeof(double));
124 
125  if (!PRN_arr) {
126  fprintf(stderr, "failed to allocate memory!\n");
127  return 1;
128  }
129 
130  av_lfg_init(&state, 42);
131  for (i = 0; i < tot_samp; i += 2) {
132  double bmg_out[2];
133  av_bmg_get(&state, bmg_out);
134  PRN_arr[i ] = bmg_out[0] * stddev + mean;
135  PRN_arr[i+1] = bmg_out[1] * stddev + mean;
136  samp_mean += PRN_arr[i] + PRN_arr[i+1];
137  samp_stddev += PRN_arr[i] * PRN_arr[i] + PRN_arr[i+1] * PRN_arr[i+1];
138  printf("PRN%d : %f\n"
139  "PRN%d : %f\n",
140  i, PRN_arr[i], i+1, PRN_arr[i+1]);
141  }
142  samp_mean /= tot_samp;
143  samp_stddev /= (tot_samp - 1);
144  samp_stddev -= (tot_samp * 1.0 / (tot_samp - 1))*samp_mean*samp_mean;
145  samp_stddev = sqrt(samp_stddev);
146  Z = (mean - samp_mean) / (stddev / sqrt(tot_samp));
147  {
148  int x, y, a, b, flag = 0;
149 
150  if (Z < 0.0) {
151  flag = !flag;
152  Z = Z * -1.0;
153  }
154 
155  a = (int)(Z * 100);
156  b = ((int)Z * 100);
157  x = Z * 10;
158  y = (b > 0) ? a % b : a;
159  y = y % 10;
160  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  p_value = flag ? 1 - Z_TABLE[x][y] : Z_TABLE[x][y];
166  }
167 
168 SKIP: for (i = 0; i < tot_samp; ++i) {
169 
170  if ( i < (tot_samp - 1)) {
171  double H_diff;
172  H_diff = inv_cdf((i + 2.0 - (3.0/8.0)) / (tot_samp + (1.0/4.0)));
173  H_diff -= inv_cdf((i + 1.0 - (3.0/8.0)) / (tot_samp + (1.0/4.0)));
174 
175  QH += ((PRN_arr[i + 1] - PRN_arr[i]) / H_diff);
176  }
177  }
178  QH = 1.0 - QH / ((tot_samp - 1.0) * samp_stddev);
179 
180  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  av_freep(&PRN_arr);
190  }
191  return 0;
192 }
r
const char * r
Definition: vf_curves.c:127
av_lfg_init
av_cold void av_lfg_init(AVLFG *c, unsigned int seed)
Definition: lfg.c:32
u
#define u(width, name, range_min, range_max)
Definition: cbs_h2645.c:251
normalize.log
log
Definition: normalize.py:21
b
#define b
Definition: input.c:41
inv_cdf
static double inv_cdf(double u)
Definition: lfg.c:58
main
int main(void)
Definition: lfg.c:103
AV_LOG_ERROR
#define AV_LOG_ERROR
Something went wrong and cannot losslessly be recovered.
Definition: log.h:209
av_lfg_get
static unsigned int av_lfg_get(AVLFG *c)
Get the next random unsigned 32-bit number using an ALFG.
Definition: lfg.h:53
lfg.h
av_bmg_get
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
fabs
static __device__ float fabs(float a)
Definition: cuda_runtime.h:182
NULL
#define NULL
Definition: coverity.c:32
timer.h
c
Undefined Behavior In the C some operations are like signed integer dereferencing freed accessing outside allocated Undefined Behavior must not occur in a C it is not safe even if the output of undefined operations is unused The unsafety may seem nit picking but Optimizing compilers have in fact optimized code on the assumption that no undefined Behavior occurs Optimizing code based on wrong assumptions can and has in some cases lead to effects beyond the output of computations The signed integer overflow problem in speed critical code Code which is highly optimized and works with signed integers sometimes has the problem that often the output of the computation does not c
Definition: undefined.txt:32
AVLFG
Context structure for the Lagged Fibonacci PRNG.
Definition: lfg.h:33
printf
printf("static const uint8_t my_array[100] = {\n")
a
The reader does not expect b to be semantically here and if the code is changed by maybe adding a a division or other the signedness will almost certainly be mistaken To avoid this confusion a new type was SUINT is the C unsigned type but it holds a signed int to use the same example SUINT a
Definition: undefined.txt:41
AV_LOG_INFO
#define AV_LOG_INFO
Standard information.
Definition: log.h:220
Z_TABLE
static const double Z_TABLE[31][10]
Definition: lfg.c:24
flag
#define flag(name)
Definition: cbs_av1.c:474
log.h
i
#define i(width, name, range_min, range_max)
Definition: cbs_h2645.c:256
av_malloc_array
#define av_malloc_array(a, b)
Definition: tableprint_vlc.h:31
state
static struct @465 state
mean
static float mean(const float *input, int size)
Definition: vf_nnedi.c:866
SKIP
#define SKIP
mem.h
av_freep
#define av_freep(p)
Definition: tableprint_vlc.h:34
av_log
#define av_log(a,...)
Definition: tableprint_vlc.h:27
Generated on Thu Nov 28 2024 19:22:49 for FFmpeg by   doxygen 1.8.17