From 3ba69a15c680596feeef1d213ce2b85823dde985 Mon Sep 17 00:00:00 2001
From: Stefano Sabatini <stefano.sabatini-lala@poste.it>
Date: Sat, 12 Dec 2009 16:24:37 +0000
Subject: [PATCH] Implement av_bmg_next(), a Box-Muller Gaussian random
 generator.

See the thread:
"[FFmpeg-devel] [PATCH] Box-Muller gaussian generator".

Originally committed as revision 20808 to svn://svn.ffmpeg.org/ffmpeg/trunk
---
 libavutil/avutil.h |  4 ++--
 libavutil/lfg.c    | 33 +++++++++++++++++++++++++++++++++
 libavutil/lfg.h    |  8 ++++++++
 3 files changed, 43 insertions(+), 2 deletions(-)

diff --git a/libavutil/avutil.h b/libavutil/avutil.h
index d799211ffc..0a398680f4 100644
--- a/libavutil/avutil.h
+++ b/libavutil/avutil.h
@@ -35,8 +35,8 @@
 #define AV_VERSION(a, b, c) AV_VERSION_DOT(a, b, c)
 
 #define LIBAVUTIL_VERSION_MAJOR 50
-#define LIBAVUTIL_VERSION_MINOR  5
-#define LIBAVUTIL_VERSION_MICRO  1
+#define LIBAVUTIL_VERSION_MINOR  6
+#define LIBAVUTIL_VERSION_MICRO  0
 
 #define LIBAVUTIL_VERSION_INT   AV_VERSION_INT(LIBAVUTIL_VERSION_MAJOR, \
                                                LIBAVUTIL_VERSION_MINOR, \
diff --git a/libavutil/lfg.c b/libavutil/lfg.c
index e876d18fef..a6828634ae 100644
--- a/libavutil/lfg.c
+++ b/libavutil/lfg.c
@@ -39,6 +39,21 @@ void av_cold av_lfg_init(AVLFG *c, unsigned int seed){
     c->index=0;
 }
 
+void av_bmg_get(AVLFG *lfg, double out[2])
+{
+    double x1, x2, w;
+
+    do {
+        x1 = 2.0/UINT_MAX*av_lfg_get(lfg) - 1.0;
+        x2 = 2.0/UINT_MAX*av_lfg_get(lfg) - 1.0;
+        w = x1*x1 + x2*x2;
+    } while (w >= 1.0);
+
+    w = sqrt((-2.0 * log(w)) / w);
+    out[0] = x1 * w;
+    out[1] = x2 * w;
+}
+
 #ifdef TEST
 #include "log.h"
 #include "common.h"
@@ -59,6 +74,24 @@ int main(void)
         STOP_TIMER("624 calls of av_lfg_get");
     }
     av_log(NULL, AV_LOG_ERROR, "final value:%X\n", x);
+
+    /* BMG usage example */
+    {
+        double mean   = 1000;
+        double stddev = 53;
+
+        av_lfg_init(&state, 42);
+
+        for (i = 0; i < 1000; i += 2) {
+            double bmg_out[2];
+            av_bmg_get(&state, bmg_out);
+            av_log(NULL, AV_LOG_INFO,
+                   "%f\n%f\n",
+                   bmg_out[0] * stddev + mean,
+                   bmg_out[1] * stddev + mean);
+        }
+    }
+
     return 0;
 }
 #endif
diff --git a/libavutil/lfg.h b/libavutil/lfg.h
index 3250c18e79..ac89d120d5 100644
--- a/libavutil/lfg.h
+++ b/libavutil/lfg.h
@@ -51,4 +51,12 @@ static inline unsigned int av_mlfg_get(AVLFG *c){
     return c->state[c->index++ & 63] = 2*a*b+a+b;
 }
 
+/**
+ * Gets the next two numbers generated by a Box-Muller Gaussian
+ * generator using the random numbers issued by lfg.
+ *
+ * @param out[2] array where are placed the two generated numbers
+ */
+void av_bmg_get(AVLFG *lfg, double out[2]);
+
 #endif /* AVUTIL_LFG_H */