ARB
admath.cxx
Go to the documentation of this file.
1 // =============================================================== //
2 // //
3 // File : admath.cxx //
4 // Purpose : //
5 // //
6 // Institute of Microbiology (Technical University Munich) //
7 // http://www.arb-home.de/ //
8 // //
9 // =============================================================== //
10 
11 #include <cmath>
12 #include <ctime>
13 #include <algorithm>
14 
15 #include <boost/random/linear_congruential.hpp>
16 #include <boost/random/uniform_real.hpp>
17 #include <boost/random/variate_generator.hpp>
18 
19 #include "gb_local.h"
20 
21 double GB_log_fak(int n) {
22  // returns log(n!)
23  static int cached = 0;
24  static double *res = NULp;
25 
26  if (n<=1) return 0.0; // log 1 = 0
27 
28  if (n > cached) {
29  int new_size = std::max(n + 50, cached*4/3);
30  gb_assert(new_size>(cached+1));
31  ARB_realloc(res, new_size);
32 
33  res[0] = 0.0;
34  for (int i=cached+1; i<new_size; i++) {
35  res[i] = res[i-1]+log((double)i);
36  }
37  cached = new_size-1;
38  }
39  return res[n];
40 }
41 
42 // -------------------------------------------
43 // portable random number generation
44 
46  typedef boost::minstd_rand base_generator_type;
47  typedef boost::uniform_real<> distribution_type;
48  typedef boost::variate_generator<base_generator_type&, distribution_type> generator_type;
49 
50  unsigned usedSeed; // unused (for debugging purposes only)
51  base_generator_type base;
52  generator_type *generator;
53 
54 
55 public:
56  RandomNumberGenerator(unsigned seed) :
57  usedSeed(seed),
58  base(seed),
59  generator(new generator_type(base, distribution_type(0, 1)))
60  {}
61  ~RandomNumberGenerator() { delete generator; }
62 
63  void reseed(unsigned seed) {
64  base.seed(seed);
65 
66  delete generator;
67  generator = new generator_type(base, distribution_type(0, 1));
68  usedSeed = seed;
69 
70  get_rand(); // the 1st random number depends too much on the seed => drop it
71  }
72 
73  double get_rand() {
74  return (*generator)();
75  }
76 };
77 
78 static RandomNumberGenerator gen(time(NULp));
79 
80 void GB_random_seed(unsigned seed) {
84  if (!seed) seed = 1; // fails assertion if seed==0
85  gen.reseed(seed);
86 }
87 
88 int GB_random(int range) {
89  // produces a random number in range [0 .. range-1]
90  return int(gen.get_rand()*((double)range));
91 }
92 
93 // --------------------------------------------------------------------------------
94 
95 #ifdef UNIT_TESTS
96 #ifndef TEST_UNIT_H
97 #include <test_unit.h>
98 #endif
99 
100 void TEST_random() {
101  // test random numbers are portable (=reproducable on all tested systems)
102  GB_random_seed(42);
103 
104  TEST_EXPECT_EQUAL(GB_random(1000), 571);
105  TEST_EXPECT_EQUAL(GB_random(1000), 256);
106  TEST_EXPECT_EQUAL(GB_random(1000), 447);
107  TEST_EXPECT_EQUAL(GB_random(1000), 654);
108 
109  int rmax = -1;
110  int rmin = 99999;
111  for (int i = 0; i<2000; ++i) {
112  int r = GB_random(1000);
113  rmax = std::max(r, rmax);
114  rmin = std::min(r, rmin);
115  }
116 
117  TEST_EXPECT_EQUAL(rmin, 0);
118  TEST_EXPECT_EQUAL(rmax, 999);
119 }
120 
121 void TEST_log_fak() {
122  const double EPS = 0.000001;
123 
124  for (int loop = 0; loop<=1; ++loop) {
125  TEST_ANNOTATE(GBS_global_string("loop=%i", loop));
126 
127  TEST_EXPECT_SIMILAR(GB_log_fak(0), 0.0, EPS);
128  TEST_EXPECT_SIMILAR(GB_log_fak(1), 0.0, EPS);
129  TEST_EXPECT_SIMILAR(GB_log_fak(2), 0.693147, EPS);
130  TEST_EXPECT_SIMILAR(GB_log_fak(3), 1.791759, EPS);
131  TEST_EXPECT_SIMILAR(GB_log_fak(10), 15.104413, EPS);
132  TEST_EXPECT_SIMILAR(GB_log_fak(30), 74.658236, EPS);
133  TEST_EXPECT_SIMILAR(GB_log_fak(70), 230.439044, EPS);
134  TEST_EXPECT_SIMILAR(GB_log_fak(99), 359.134205, EPS);
135  }
136 }
137 
138 #endif // UNIT_TESTS
139 
140 // --------------------------------------------------------------------------------
#define TEST_EXPECT_SIMILAR(expr, want, epsilon)
Definition: test_unit.h:1287
double GB_log_fak(int n)
Definition: admath.cxx:21
int GB_random(int range)
Definition: admath.cxx:88
const char * GBS_global_string(const char *templat,...)
Definition: arb_msg.cxx:204
void GB_random_seed(unsigned seed)
Definition: admath.cxx:80
#define EPS
Definition: awt_canvas.hxx:291
#define gb_assert(cond)
Definition: arbdbt.h:11
void ARB_realloc(TYPE *&tgt, size_t nelem)
Definition: arb_mem.h:43
#define NULp
Definition: cxxforward.h:97
static RandomNumberGenerator gen(time(NULp))
#define min(a, b)
Definition: f2c.h:153
RandomNumberGenerator(unsigned seed)
Definition: admath.cxx:56
void reseed(unsigned seed)
Definition: admath.cxx:63
#define TEST_EXPECT_EQUAL(expr, want)
Definition: test_unit.h:1283
#define max(a, b)
Definition: f2c.h:154