diff --git a/libavutil/tests/lfg.c b/libavutil/tests/lfg.c
index 1425e02..645d39d 100644
--- a/libavutil/tests/lfg.c
+++ b/libavutil/tests/lfg.c
@@ -20,6 +20,59 @@
 #include "libavutil/timer.h"
 #include "libavutil/lfg.h"
 
+#define TEST    1
+
+// Inverse cumulative distribution function
+static double inv_cdf(double u)
+{
+    const double a[4] = { 2.50662823884,
+                         -18.61500062529,
+                          41.39119773534,
+                         -25.44106049637};
+
+    const double b[4] = {-8.47351093090,
+                          23.08336743743,
+                         -21.06224101826,
+                           3.13082909833};
+
+    const double c[9] = {0.3374754822726147,
+                          0.9761690190917186,
+                          0.1607979714918209,
+                          0.0276438810333863,
+                          0.0038405729373609,
+                          0.0003951896511919,
+                          0.0000321767881768,
+                          0.0000002888167364,
+                          0.0000003960315187};
+
+    double r;
+    double x = u - 0.5;
+
+    // Beasley-Springer
+#if TEST == 0
+    if (fabs(x) < 0.42) {
+#endif
+#if TEST == 1
+    if (u > 0.08 && u < 0.92) {
+#endif
+        double y = x * x;
+        r        = x * (((a[3]*y+a[2])*y+a[1])*y+a[0]) /
+                        ((((b[3]*y+b[2])*y+b[1])*y+b[0])*y+1.0);
+    }
+    else {// Moro
+        r = u;
+        if (x > 0.0)
+            r = 1.0 - u;
+        r = log(-log(r));
+        r = c[0] + r*(c[1]+r*(c[2]+r*(c[3]+r*(c[4]+r*(c[5]+r*(c[6]+
+                                                        r*(c[7]+r*c[8])))))));
+        if (x < 0.0)
+            r = -r;
+    }
+
+    return r;
+}
+
 int main(void)
 {
     int x = 0;
@@ -37,38 +90,58 @@ int main(void)
     }
     av_log(NULL, AV_LOG_ERROR, "final value:%X\n", x);
 
-    /* BMG usage example */
+    /*  Chen-Shapiro Normality Test */
     {
-        double mean   = 1000;
-        double stddev = 53;
-        double samp_mean = 0.0, samp_stddev = 0.0;
-        double samp0, samp1;
+        double QH           = 0.0;
+        const double stddev = 53.0;
+        double samp_stddev  = 0.0;
+        const double mean   = 1000;
+        double samp_mean    = 0.0;
+        const int tot_samp  = 2000;
+        double *PRN_arr     = av_malloc_array(tot_samp, sizeof(double));
 
-        av_lfg_init(&state, 42);
+        av_lfg_init(&state, 0x7a69);
+        if (!PRN_arr) {
+            fprintf(stderr, "failed to allocate memory!\n");
+            return 1;
+        }
 
-        for (i = 0; i < 1000; i += 2) {
+        for (i = 0; i < tot_samp; i += 2) {
             double bmg_out[2];
             av_bmg_get(&state, bmg_out);
-            samp0 = bmg_out[0] * stddev + mean;
-            samp1 = bmg_out[1] * stddev + mean;
-            samp_mean += samp0 + samp1;
-            samp_stddev += samp0 * samp0 + samp1 * samp1;
-            av_log(NULL, AV_LOG_INFO,
-                   "%f\n%f\n",
-                   samp0,
-                   samp1);
+            PRN_arr[i  ] = bmg_out[0] * stddev + mean;
+            PRN_arr[i+1] = bmg_out[1] * stddev + mean;
+            samp_mean   += PRN_arr[i] + PRN_arr[i+1];
+            av_log(NULL, AV_LOG_INFO, "PRN%d : %f\n"
+                                      "PRN%d : %f\n",
+                                       i, PRN_arr[i], i+1, PRN_arr[i+1]);
         }
-        /* TODO: add proper normality test */
-        samp_mean /= 1000;
-        samp_stddev /= 999;
-        samp_stddev -= (1000.0/999.0)*samp_mean*samp_mean;
-        samp_stddev = sqrt(samp_stddev);
-        av_log(NULL, AV_LOG_INFO, "sample mean  : %f\n"
-                                  "true mean    : %f\n"
-                                  "sample stddev: %f\n"
-                                  "true stddev  : %f\n",
-                                   samp_mean, mean, samp_stddev, stddev);
-    }
+        samp_mean /= tot_samp;
+
+        for (i = 0; i < tot_samp; ++i) {
+            double temp;
+            temp         = PRN_arr[i] - samp_mean;
+            temp        *= temp;
+            samp_stddev += temp;
 
+            if ( i < (tot_samp - 1)) {
+                double H_diff;
+                H_diff  = inv_cdf((i + 2.0 - (3.0/8.0)) / (tot_samp + (1.0/4.0)));
+                H_diff -= inv_cdf((i + 1.0 - (3.0/8.0)) / (tot_samp + (1.0/4.0)));
+
+                QH     += ((PRN_arr[i + 1] - PRN_arr[i]) / H_diff);
+            }
+        }
+        samp_stddev  = sqrt(samp_stddev / (tot_samp - 1.0));
+        QH           = 1.0 - QH / ((tot_samp - 1.0) * samp_stddev);
+
+        av_log(NULL, AV_LOG_INFO, "sample mean      : %f\n"
+                                  "true mean        : %f\n"
+                                  "sample stddev    : %f\n"
+                                  "true stddev      : %f\n"
+                                  "Chen-Shapiro [QH]: %f\n",
+                                   samp_mean, mean, samp_stddev, stddev, QH);
+        av_freep(&PRN_arr);
+    }
     return 0;
 }
diff --git a/tests/fate/libavutil.mak b/tests/fate/libavutil.mak
index a7bf739..c3b9091 100644
--- a/tests/fate/libavutil.mak
+++ b/tests/fate/libavutil.mak
@@ -101,6 +101,11 @@ FATE_LIBAVUTIL += fate-imgutils
 fate-imgutils: libavutil/tests/imgutils$(EXESUF)
 fate-imgutils: CMD = run libavutil/tests/imgutils
 
+FATE_LIBAVUTIL += fate-lfg
+fate-lfg: libavutil/tests/lfg$(EXESUF)
+fate-lfg: CMD = run libavutil/tests/lfg
+fate-lfg: REF = /dev/null
+
 FATE_LIBAVUTIL += fate-md5
 fate-md5: libavutil/tests/md5$(EXESUF)
 fate-md5: CMD = run libavutil/tests/md5
