Re: random number generators for calc

> > just happened to reimplement RAND() with a B.A Wichmann
> > and I.D Hill generator.. see
>
>       As they are entitled to. No doubt we'll come up with a comprehensive,
> nicely architected and beautiful solution of our own :-)

Ok, I'll go ahead with a suggestion then. :)

Attached patch is just a simple wrapper around boost, to be used like

sc::rng::seed() to replace libc srand()
sc::rng::rand() to replace libc rand()

and a few more distributions.

I've also modified ScInterpreter::ScRandom()
to simply call sc::rng::uniform()
which should solve bug 33365.

I've not modified RANDBETWEEN(a,b) yet but this could simply call
sc::rng::uniform_int(a,b).

Also, the rand in Basic could be changed in the future.

Small problems:
 - compiler warnings from within boost
 - few asserts might need to be replaced

From 569a0e05f1c2320041786be051842a552aca885c Mon Sep 17 00:00:00 2001
From: tino <ttk448@gmail.com>
Date: Thu, 6 Dec 2012 14:05:11 +0000
Subject: [PATCH] fdo#33365 added wrapper for boost random, use that in RAND()

Change-Id: Iafc524d12c76423f74dc16b42595e52fbc5a1e54
---
 sc/Library_sc.mk                 |   1 +
 sc/source/core/data/global.cxx   |   2 +
 sc/source/core/inc/random.hxx    |  35 ++++++
 sc/source/core/tool/interpr1.cxx |   4 +-
 sc/source/core/tool/random.cxx   | 255 +++++++++++++++++++++++++++++++++++++++
 5 files changed, 296 insertions(+), 1 deletion(-)
 create mode 100644 sc/source/core/inc/random.hxx
 create mode 100644 sc/source/core/tool/random.cxx

diff --git a/sc/Library_sc.mk b/sc/Library_sc.mk
index 348300f..10918e4 100644
--- a/sc/Library_sc.mk
+++ b/sc/Library_sc.mk
@@ -213,6 +213,7 @@ $(eval $(call gb_Library_add_exception_objects,sc,\
        sc/source/core/tool/progress \
        sc/source/core/tool/queryentry \
        sc/source/core/tool/queryparam \
+       sc/source/core/tool/random \
        sc/source/core/tool/rangelst \
        sc/source/core/tool/rangenam \
        sc/source/core/tool/rangeseq \
diff --git a/sc/source/core/data/global.cxx b/sc/source/core/data/global.cxx
index 4e449c5..e18d241 100644
--- a/sc/source/core/data/global.cxx
+++ b/sc/source/core/data/global.cxx
@@ -75,6 +75,7 @@
 #include "sc.hrc"
 #include "scmod.hxx"
 #include "appoptio.hxx"
+#include "random.hxx"
 
 // -----------------------------------------------------------------------
 
@@ -557,6 +558,7 @@ void ScGlobal::Init()
     // names from the compiler.
     ScParameterClassification::Init();
     srand( (unsigned) time( NULL ) );       // Random Seed Init fuer Interpreter
+    sc::rng::seed( time( NULL ) );          // seed for libc rand() replacement
 
     InitAddIns();
 
diff --git a/sc/source/core/inc/random.hxx b/sc/source/core/inc/random.hxx
new file mode 100644
index 0000000..dc068b1
--- /dev/null
+++ b/sc/source/core/inc/random.hxx
@@ -0,0 +1,35 @@
+/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+/*
+ * This file is part of the LibreOffice project.
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ */
+
+#ifndef SC_RANDOM_HXX
+#define SC_RANDOM_HXX
+
+namespace sc
+{
+namespace rng
+{
+
+void seed(int i);       // set initial seed (equivalent of libc srand())
+uint32_t max();         // max number by rand() (equivalent of RAND_MAX)
+uint32_t min();         // min number by rand() (most likely 0)
+uint32_t rand();        // replacement for libc's rand()
+
+int uniform_int(int from, int to);  // uniform distribution given int range
+double uniform();                   // uniform distribution in [0,1)
+double normal();                    // normal distribution N(0,1)
+double normal(double mu, double sigma); // N(mu, sigma^2)
+double lognormal();                 // lognormal distribution: exp(N(0,1))
+double lognormal(double mu, double sigma);
+
+} // namespace
+} // namespace
+
+#endif
+
+/* vim:set shiftwidth=4 softtabstop=4 expandtab: */
diff --git a/sc/source/core/tool/interpr1.cxx b/sc/source/core/tool/interpr1.cxx
index 43a21eb..adecc75 100644
--- a/sc/source/core/tool/interpr1.cxx
+++ b/sc/source/core/tool/interpr1.cxx
@@ -44,6 +44,7 @@
 #include "globstr.hrc"
 #include "attrib.hxx"
 #include "jumpmatrix.hxx"
+#include "random.hxx"
 
 #include <comphelper/processfactory.hxx>
 #include <comphelper/string.hxx>
@@ -1711,7 +1712,8 @@ void ScInterpreter::ScPi()
 void ScInterpreter::ScRandom()
 {
     RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScRandom" );
-    PushDouble((double)rand() / ((double)RAND_MAX+1.0));
+    //PushDouble((double)rand() / ((double)RAND_MAX+1.0));
+    PushDouble(sc::rng::uniform());
 }
 
 
diff --git a/sc/source/core/tool/random.cxx b/sc/source/core/tool/random.cxx
new file mode 100644
index 0000000..4195e26
--- /dev/null
+++ b/sc/source/core/tool/random.cxx
@@ -0,0 +1,255 @@
+/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+/*
+ * This file is part of the LibreOffice project.
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ *
+ * Contributor(s):
+ *   Copyright (C) 2012 Tino Kluge <tino.kluge@hrz.tu-chemnitz.de>
+ *
+ */
+
+#include <boost/random.hpp>
+
+// this is nothing but a simple wrapper around
+// the boost random generators
+
+namespace sc
+{
+namespace rng
+{
+
+// underlying random number generator
+// boost::mt19937 implements the Mersenne twister algorithm which
+// is fast and has good statistical properties, it produces integers
+// in the range of [0, 2^32-1] internally
+// memory requirement: 625*sizeof(uint32_t)
+// http://en.wikipedia.org/wiki/Mersenne_twister
+#define BOOST_RNG_ALGO  boost::mt19937
+BOOST_RNG_ALGO global_rng;
+
+// initialises the state of the global random number generator
+// should only be called once at the start of the main programme
+// (note, a few boost::variate_generator<> (like normal) have their
+// own state which would need a reset as well to guarantee identical
+// sequence of numbers, e.g. via myrand.distribution().reset())
+void seed(int i)
+{
+    global_rng.seed(i);
+}
+
+// uniform distributed integers, for boost::mt19937 [0, 2^32-1]
+uint32_t rand()
+{
+    return global_rng();
+}
+
+// biggest possible number generated by global_rng (like RAND_MAX)
+uint32_t max()
+{
+    return global_rng.max();
+}
+
+// smallest possible number generated by global_rng (0 for boost::mt19937)
+uint32_t min()
+{
+    return global_rng.min();
+}
+
+
+// all rand function below use global_rng and apply certain
+// transformations to obtain specific non-uniform distributions
+// if the inverse cdf (cumulative distribution function) F_inv(x)
+// is easily computable then one would simply do:
+//
+//  rand_nonuniform() = F_inv(rand_uniform01())
+//
+// even if F_inv(x) is not easily computable there may be other
+// simple ways, all of this intelligence is hidden within the
+// abstraction of boost::variate_generator<> so we don't have
+// to think too much
+
+
+// uniform integers in given range [from ... to]
+int uniform_int(int from, int to)
+{
+    assert(from<=to);
+    return from + rng::rand() % (to-from+1);
+}
+
+// uniform [0,1) or [a,b) distribution
+double uniform()
+{
+    static boost::uniform_01<BOOST_RNG_ALGO&> myrand(global_rng);
+    return myrand();
+}
+double uniform(double a, double b)
+{
+    assert(a<=b);
+    return a + (b-a)*uniform();
+}
+
+// normal N(0,1) or N(mu, sigma^2) distribution
+double normal()
+{
+    static boost::variate_generator<BOOST_RNG_ALGO&,
+           boost::normal_distribution<> >
+           myrand(global_rng, boost::normal_distribution<>());
+    return myrand();
+}
+double normal(double mu, double sigma)
+{
+    return mu + sigma*normal();
+}
+
+// Lognormal distribution
+double lognormal()
+{
+    return exp(normal());
+}
+double lognormal(double mu, double sigma)
+{
+    return exp(normal(mu,sigma));
+}
+
+
+
+
+
+// requires newer version of boost (at least 1.47 or so)
+// a list of all available distributions see for example
+// 
http://www.boost.org/doc/libs/1_52_0/doc/html/boost_random/reference.html#boost_random.reference.distributions
+#if 0
+
+
+// uniform integers in given range [from ... to]
+// should replace the above uniform_int() once newer boost is available
+/*
+int uniform_int(int from, int to){
+   assert(from<=to);
+   // set up random generator based on global_rng and specified distribution
+   static boost::variate_generator<BOOST_RNG_ALGO&,
+               boost::random::uniform_int_distribution<> >
+               myrand(global_rng,boost::random::uniform_int_distribution<>());
+   // change distribution parameters according to input values
+   myrand.distribution().param(
+         boost::random::uniform_int_distribution<>::param_type(from,to));
+   // return a generated random value
+   return myrand();
+}
+*/
+
+double normal2(double mu, double sigma)
+{
+    static boost::variate_generator<BOOST_RNG_ALGO&,
+           boost::random::normal_distribution<> >
+           myrand(global_rng, boost::random::normal_distribution<>());
+    myrand.distribution().param(
+        boost::random::normal_distribution<>::param_type(mu,sigma));
+    return myrand();
+}
+
+// Cauchy distribution
+double cauchy(double x0, double gamma)
+{
+    static boost::variate_generator<BOOST_RNG_ALGO&,
+           boost::random::cauchy_distribution<> >
+           myrand(global_rng, boost::random::cauchy_distribution<>());
+    myrand.distribution().param(
+        boost::random::cauchy_distribution<>::param_type(x0,gamma));
+    return myrand();
+}
+
+// Poisson distribution
+double poisson(double lambda)
+{
+    static boost::variate_generator<BOOST_RNG_ALGO&,
+           boost::random::poisson_distribution<> >
+           myrand(global_rng, boost::random::poisson_distribution<>());
+    myrand.distribution().param(
+        boost::random::poisson_distribution<>::param_type(lambda));
+    return myrand();
+}
+
+// Triangular distribution, density function looks like a triangle
+double triangle(double a, double b, double c)
+{
+    assert(a<=c && c<=b);
+    static boost::variate_generator<BOOST_RNG_ALGO&,
+           boost::random::triangle_distribution<> >
+           myrand(global_rng, boost::random::triangle_distribution<>());
+    myrand.distribution().param(
+        boost::random::triangle_distribution<>::param_type(a,b,c));
+    return myrand();
+}
+
+// exponential distribution
+double exponential(double lambda)
+{
+    assert(lambda>=0.0);
+    static boost::variate_generator<BOOST_RNG_ALGO&,
+           boost::random::exponential_distribution<> >
+           myrand(global_rng,boost::random::exponential_distribution<>());
+    myrand.distribution().param(
+        boost::random::exponential_distribution<>::param_type(lambda));
+    return myrand();
+}
+
+
+
+// Bernoulli distribution, returns true with probability p
+bool bernoulli(double p)
+{
+    static boost::variate_generator<BOOST_RNG_ALGO&,
+           boost::random::bernoulli_distribution<> >
+           myrand(global_rng, boost::random::bernoulli_distribution<>());
+    myrand.distribution().param(
+        boost::random::bernoulli_distribution<>::param_type(p) );
+    return myrand();
+}
+
+// binomial distribution, returns the sum of n Bernoulli trials
+int binomial(double p, int n)
+{
+    static boost::variate_generator<BOOST_RNG_ALGO&,
+           boost::random::binomial_distribution<> >
+           myrand(global_rng, boost::random::binomial_distribution<>());
+    myrand.distribution().param(
+        boost::random::binomial_distribution<>::param_type(p,n) );
+    return myrand();
+}
+
+// geometric distribution, number of Bernoulli trials till success
+int geometric(double p)
+{
+    static boost::variate_generator<BOOST_RNG_ALGO&,
+           boost::random::geometric_distribution<> >
+           myrand(global_rng, boost::random::geometric_distribution<>());
+    myrand.distribution().param(
+        boost::random::geometric_distribution<>::param_type(p));
+    return myrand();
+}
+
+
+
+// Chi-squared Chi^2(k) distribution
+double chisq(int k)
+{
+    assert(k>0);
+    static boost::variate_generator<BOOST_RNG_ALGO&,
+           boost::random::chi_squared_distribution<> >
+           myrand(global_rng,boost::random::chi_squared_distribution<>());
+    myrand.distribution().param(
+        boost::random::chi_squared_distribution<>::param_type(k));
+    return myrand();
+}
+
+
+#endif // newer boost version
+
+} // namespace
+} // namespace
+
+/* vim:set shiftwidth=4 softtabstop=4 expandtab: */
-- 
1.7.11.2