trng 4.16

AUTHORS

@@ -1,2 +1,2 @@
 Heiko Bauke
-heiko.bauke@snafu.de
+trng@mail.de

TODO

@@ -1,18 +1,7 @@
 Implement:
 
-* more distributions 
-
-  * Continous distributions 
-    - Kolmogorov-Smirnov distribution
-      http://en.wikipedia.org/wiki/Kolmogorov-Smirnov
-    - Levy distribution
-      http://en.wikipedia.org/wiki/L%C3%A9vy_distribution
-
-  * and other distributions, see
-    http://en.wikipedia.org/wiki/List_of_probability_distributions
-
-  * multinomial distribution, see
-    - http://en.wikipedia.org/wiki/Multinomial_distribution
+* more distributions, see
+  http://en.wikipedia.org/wiki/List_of_probability_distributions
 
 * ISAAC generator
   http://burtleburtle.net/bob/rand/isaacafa.html
@@ -29,5 +18,4 @@ Implement:
   http://www.thesalmons.org/john/random123/papers/random123sc11.pdf
   Parallel Random Numbers: As Easy as 1, 2, 3
 
-* Intel hardware random number generators 
-  http://software.intel.com/en-us/articles/user-manual-for-the-rdrand-library-linux-version/?cid=sw:ISNnews_87_Eng_7264
+* Intel hardware random number generators (RDRAND)

configure

@@ -1,6 +1,6 @@
 #! /bin/sh
 # Guess values for system-dependent variables and create Makefiles.
-# Generated by GNU Autoconf 2.69 for trng 4.15.
+# Generated by GNU Autoconf 2.69 for trng 4.16.
 #
 # Report bugs to <[email protected]>.
 #
@@ -589,8 +589,8 @@ MAKEFLAGS=
 # Identity of this package.
 PACKAGE_NAME='trng'
 PACKAGE_TARNAME='trng'
-PACKAGE_VERSION='4.15'
-PACKAGE_STRING='trng 4.15'
+PACKAGE_VERSION='4.16'
+PACKAGE_STRING='trng 4.16'
 PACKAGE_BUGREPORT='[email protected]'
 PACKAGE_URL=''
 
@@ -1326,7 +1326,7 @@ if test "$ac_init_help" = "long"; then
   # Omit some internal or obsolete options to make the list less imposing.
   # This message is too long to be a string in the A/UX 3.1 sh.
   cat <<_ACEOF
-\`configure' configures trng 4.15 to adapt to many kinds of systems.
+\`configure' configures trng 4.16 to adapt to many kinds of systems.
 
 Usage: $0 [OPTION]... [VAR=VALUE]...
 
@@ -1396,7 +1396,7 @@ fi
 
 if test -n "$ac_init_help"; then
   case $ac_init_help in
-     short | recursive ) echo "Configuration of trng 4.15:";;
+     short | recursive ) echo "Configuration of trng 4.16:";;
    esac
   cat <<\_ACEOF
 
@@ -1507,7 +1507,7 @@ fi
 test -n "$ac_init_help" && exit $ac_status
 if $ac_init_version; then
   cat <<\_ACEOF
-trng configure 4.15
+trng configure 4.16
 generated by GNU Autoconf 2.69
 
 Copyright (C) 2012 Free Software Foundation, Inc.
@@ -2064,7 +2064,7 @@ cat >config.log <<_ACEOF
 This file contains any messages produced by compilers while
 running configure, to aid debugging if configure makes a mistake.
 
-It was created by trng $as_me 4.15, which was
+It was created by trng $as_me 4.16, which was
 generated by GNU Autoconf 2.69.  Invocation command line was
 
   $ $0 $@
@@ -2927,7 +2927,7 @@ fi
 
 # Define the identity of the package.
  PACKAGE='trng'
- VERSION='4.15'
+ VERSION='4.16'
 
 
 cat >>confdefs.h <<_ACEOF
@@ -3021,7 +3021,7 @@ END
 fi
 
 
-$as_echo "#define TRNG_VERSION 4.15" >>confdefs.h
+$as_echo "#define TRNG_VERSION 4.16" >>confdefs.h
 
 
 ac_ext=cpp
@@ -16445,7 +16445,7 @@ cat >>$CONFIG_STATUS <<\_ACEOF || ac_write_fail=1
 # report actual input values of CONFIG_FILES etc. instead of their
 # values after options handling.
 ac_log="
-This file was extended by trng $as_me 4.15, which was
+This file was extended by trng $as_me 4.16, which was
 generated by GNU Autoconf 2.69.  Invocation command line was
 
   CONFIG_FILES    = $CONFIG_FILES
@@ -16511,7 +16511,7 @@ _ACEOF
 cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1
 ac_cs_config="`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`"
 ac_cs_version="\\
-trng config.status 4.15
+trng config.status 4.16
 configured by $0, generated by GNU Autoconf 2.69,
   with options \\"\$ac_cs_config\\"
 

configure.ac

@@ -1,7 +1,7 @@
-AC_INIT([trng], [4.15], [[email protected]])
+AC_INIT([trng], [4.16], [[email protected]])
 AM_INIT_AUTOMAKE()
 
-AC_DEFINE([TRNG_VERSION], [4.15], [TRNG version number.])
+AC_DEFINE([TRNG_VERSION], [4.16], [TRNG version number.])
 
 AC_LANG_CPLUSPLUS
 AC_LANG(C)

examples/distributions.cc

@@ -61,6 +61,7 @@
 #include <trng/rayleigh_dist.hpp>
 #include <trng/bernoulli_dist.hpp>
 #include <trng/binomial_dist.hpp>
+#include <trng/negative_binomial_dist.hpp>
 #include <trng/hypergeometric_dist.hpp>
 #include <trng/geometric_dist.hpp>
 #include <trng/poisson_dist.hpp>
@@ -294,6 +295,10 @@ int main() {
     trng::binomial_dist d(0.4, 20);
     discrete_test(d, "binomial distribution");
   }
+  {
+    trng::negative_binomial_dist d(0.4, 20);
+    discrete_test(d, "negative binomial distribution");
+  }
   {
     trng::hypergeometric_dist d(10, 5, 5);
     discrete_test(d, "hypergeometric distribution");

src/Makefile.am

@@ -27,7 +27,8 @@ pkginclude_HEADERS = config.hpp cuda.hpp utility.hpp uniformxx.hpp static_assert
 		     powerlaw_dist.hpp tent_dist.hpp weibull_dist.hpp \
 		     extreme_value_dist.hpp gamma_dist.hpp beta_dist.hpp chi_square_dist.hpp student_t_dist.hpp snedecor_f_dist.hpp\
 		     rayleigh_dist.hpp \
-		     bernoulli_dist.hpp binomial_dist.hpp hypergeometric_dist.hpp geometric_dist.hpp \
+		     bernoulli_dist.hpp binomial_dist.hpp negative_binomial_dist.hpp \
+		     hypergeometric_dist.hpp geometric_dist.hpp \
 		     poisson_dist.hpp discrete_dist.hpp fast_discrete_dist.hpp
 
 libtrng4_la_LDFLAGS = -version-info 0:0:0

src/Makefile.in

@@ -347,7 +347,8 @@ pkginclude_HEADERS = config.hpp cuda.hpp utility.hpp uniformxx.hpp static_assert
 		     powerlaw_dist.hpp tent_dist.hpp weibull_dist.hpp \
 		     extreme_value_dist.hpp gamma_dist.hpp beta_dist.hpp chi_square_dist.hpp student_t_dist.hpp snedecor_f_dist.hpp\
 		     rayleigh_dist.hpp \
-		     bernoulli_dist.hpp binomial_dist.hpp hypergeometric_dist.hpp geometric_dist.hpp \
+		     bernoulli_dist.hpp binomial_dist.hpp negative_binomial_dist.hpp \
+		     hypergeometric_dist.hpp geometric_dist.hpp \
 		     poisson_dist.hpp discrete_dist.hpp fast_discrete_dist.hpp
 
 libtrng4_la_LDFLAGS = -version-info 0:0:0

src/beta_dist.hpp

@@ -66,11 +66,11 @@ namespace trng {
       result_type norm() const { return norm_; }
     public:
       TRNG_CUDA_ENABLE
-      void alpha(result_type alpha_new) { alpha_=alpha_new; norm_=math::Beta(alpha, beta); }
+      void alpha(result_type alpha_new) { alpha_=alpha_new; norm_=math::Beta(alpha_, beta_); }
       TRNG_CUDA_ENABLE
-      void beta(result_type beta_new) { beta_=beta_new; norm_=math::Beta(alpha, beta); }
+      void beta(result_type beta_new) { beta_=beta_new; norm_=math::Beta(alpha_, beta_); }
       param_type() : 
-	alpha_(1), beta_(1), norm_(math::Beta(alpha, beta)) {
+	alpha_(1), beta_(1), norm_(math::Beta(alpha_, beta_)) {
       }
       TRNG_CUDA_ENABLE
       param_type(result_type alpha, result_type beta) : 

src/config.hpp

@@ -41,7 +41,7 @@
 #define TRNG_HAVE_OPENMP 1
 
 /* Define to 1 if Intel Threading Building Blocks are available. */
-#define TRNG_HAVE_TBB 1
+/* #undef TRNG_HAVE_TBB */
 
 /* Define to 1 if you have the `tgamma' function. */
 #define TRNG_HAVE_TGAMMA 1
@@ -53,4 +53,4 @@
 #define TRNG_HAVE_TGAMMAL 1
 
 /* TRNG version number. */
-#define TRNG_VERSION 4.15
+#define TRNG_VERSION 4.16

src/constants.hpp

@@ -61,7 +61,7 @@ namespace trng {
       TRNG_NEW_CONSTANT(float, 2.50662827463100050241577f, sqrt_2pi);
       TRNG_NEW_CONSTANT(float, 0.3183098861837906715377676f, one_over_pi);
       TRNG_NEW_CONSTANT(float, .707106781186547524400845f, one_over_sqrt_2);
-      TRNG_NEW_CONSTANT(float, .398942280401432677939946f, one_over_sqrt_2pi)
+      TRNG_NEW_CONSTANT(float, .398942280401432677939946f, one_over_sqrt_2pi);
       TRNG_NEW_CONSTANT(float, .797884560802865355879892f, sqrt_2_over_pi);
     };
 
@@ -77,7 +77,7 @@ namespace trng {
       TRNG_NEW_CONSTANT(double, 2.50662827463100050241577, sqrt_2pi);
       TRNG_NEW_CONSTANT(double, 0.3183098861837906715377676, one_over_pi);
       TRNG_NEW_CONSTANT(double, .707106781186547524400845, one_over_sqrt_2);
-      TRNG_NEW_CONSTANT(double, .398942280401432677939946, one_over_sqrt_2pi)
+      TRNG_NEW_CONSTANT(double, .398942280401432677939946, one_over_sqrt_2pi);
       TRNG_NEW_CONSTANT(double, .797884560802865355879892, sqrt_2_over_pi);
     };
 
@@ -93,7 +93,7 @@ namespace trng {
       TRNG_NEW_CONSTANT(long double, 2.50662827463100050241577l, sqrt_2pi);
       TRNG_NEW_CONSTANT(long double, .707106781186547524400845l, one_over_sqrt_2);
       TRNG_NEW_CONSTANT(long double, 0.3183098861837906715377676l, one_over_pi);
-      TRNG_NEW_CONSTANT(long double, .398942280401432677939946l, one_over_sqrt_2pi)
+      TRNG_NEW_CONSTANT(long double, .398942280401432677939946l, one_over_sqrt_2pi);
       TRNG_NEW_CONSTANT(long double, .797884560802865355879892l, sqrt_2_over_pi);
     };
 

src/discrete_dist.hpp

@@ -59,6 +59,8 @@ namespace trng {
       size_type N, offset, layers;
 
     public:
+      param_type() : P(), N(0), offset(0), layers(0) {
+      }
       template<typename iter>
       param_type(iter first, iter last) :
 	P(first, last) {
@@ -67,7 +69,7 @@ namespace trng {
 	offset=math::pow2(layers)-1;
 	P.resize(N+offset);
 	std::copy_backward(P.begin(), P.begin()+N, P.end());
-	std::fill(P.begin(), P.begin()+N, 0);
+	std::fill(P.begin(), P.begin()+offset, 0);
 	update_all_layers();
       }
       explicit param_type(int n) :
@@ -77,7 +79,7 @@ namespace trng {
 	offset=math::pow2(layers)-1;
 	P.resize(N+offset);
 	std::copy_backward(P.begin(), P.begin()+N, P.end());
-	std::fill(P.begin(), P.begin()+N, 0);
+	std::fill(P.begin(), P.begin()+offset, 0);
 	update_all_layers();
       }
     private:
@@ -130,8 +132,10 @@ namespace trng {
     // random numbers
     template<typename R>
     int operator()(R &r) {
-      double u=utility::uniformco<double>(r)*P.P[0];
-      param_type::size_type x=0;
+      if (P.N==0)
+	return -1;
+      double u(utility::uniformco<double>(r)*P.P[0]);
+      param_type::size_type x(0);
       while (x<P.offset) {
 	if (u<P.P[2*x+1]) {
 	  x=2*x+1;

src/gamma_dist.hpp

@@ -125,7 +125,7 @@ namespace trng {
 	++num_iterations;
 	y_old=y;
 	result_type f0=math::GammaP(p.kappa(), y)-x;
-	result_type f1=math::pow(y, p.kappa()-1)*math::exp(-y-ln_Gamma_kappa);
+	result_type f1=math::exp((p.kappa()-1)*math::ln(y)-y-ln_Gamma_kappa);
 	result_type f2=f1*(p.kappa()-1-y)/y;
 	y-=f0/f1*(1+f0*f2/(2*f1*f1));
       } while (num_iterations<16 &&
@@ -180,8 +180,7 @@ namespace trng {
 	return 0;
       else {
 	x/=p.theta();
-	return math::pow(x, p.kappa()-1)/
-	  (math::exp(x+math::ln_Gamma(p.kappa()))*p.theta());
+	return math::exp((p.kappa()-1)*math::ln(x)-x-math::ln_Gamma(p.kappa())) /(p.theta());
       }
     }
     // cumulative density function 

src/poisson_dist.hpp

@@ -60,11 +60,12 @@ namespace trng {
 	P_=std::vector<double>();
 	int x=0;
 	double p=0.0;
-	while (p<1.0-1.0/4096.0) {
+	while (p<1.0-512.0/4096.0) {
 	  p=math::GammaQ(x+1.0, mu_);
 	  P_.push_back(p);
 	  ++x;
 	}
+	P_.push_back(1);
       }
 
     public:
@@ -92,9 +93,13 @@ namespace trng {
     int operator()(R &r) {
       double p(utility::uniformco<double>(r));
       int x(utility::discrete(p, P.P_.begin(), P.P_.end()));
-      if (x+1==P.P_.size())
-        while (p>cdf(x))
+      if (x+1==P.P_.size()) {
+	p-=cdf(x);
+        while (p>0) {
           ++x;
+	  p-=pdf(x);
+	}
+      }
       return x;
     }
     template<typename R>