source: tags/ms_r17q3/ARBDB/admath.cxx

Last change on this file was 15587, checked in by westram, 8 years ago
  • reintegrates 'fix' into 'trunk'
    • fix OSX unittest for jenkins project names (up to max. 16 chars long)
    • optimize GB_log_fak
    • forward RegExpr-compile-errors (or at least assert on failure)
  • adds: log:branches/fix@15582:15586
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.2 KB
Line 
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
21double GB_log_fak(int n) {
22    // returns log(n!)
23    static int     cached = 0;
24    static double *res    = 0;
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
45class RandomNumberGenerator : virtual Noncopyable {
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
55public:
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
78static RandomNumberGenerator gen(time(0));
79
80void GB_random_seed(unsigned seed) {
81    /*! normally you should not use GB_random_seed.
82     * Use it only to reproduce some pseudo-random result (e.g. in unittests)
83     */
84    gen.reseed(seed);
85}
86
87int GB_random(int range) {
88    // produces a random number in range [0 .. range-1]
89    return int(gen.get_rand()*((double)range));
90}
91
92// --------------------------------------------------------------------------------
93
94#ifdef UNIT_TESTS
95#ifndef TEST_UNIT_H
96#include <test_unit.h>
97#endif
98
99void TEST_random() {
100    // test random numbers are portable (=reproducable on all tested systems)
101    GB_random_seed(42);
102
103    TEST_EXPECT_EQUAL(GB_random(1000), 571);
104    TEST_EXPECT_EQUAL(GB_random(1000), 256);
105    TEST_EXPECT_EQUAL(GB_random(1000), 447);
106    TEST_EXPECT_EQUAL(GB_random(1000), 654);
107
108    int rmax = -1;
109    int rmin = 99999;
110    for (int i = 0; i<2000; ++i) {
111        int r = GB_random(1000);
112        rmax  = std::max(r, rmax);
113        rmin  = std::min(r, rmin);
114    }
115
116    TEST_EXPECT_EQUAL(rmin, 0);
117    TEST_EXPECT_EQUAL(rmax, 999);
118}
119
120void TEST_log_fak() {
121    const double EPS = 0.000001;
122
123    for (int loop = 0; loop<=1; ++loop) {
124        TEST_ANNOTATE(GBS_global_string("loop=%i", loop));
125
126        TEST_EXPECT_SIMILAR(GB_log_fak(0),    0.0,      EPS);
127        TEST_EXPECT_SIMILAR(GB_log_fak(1),    0.0,      EPS);
128        TEST_EXPECT_SIMILAR(GB_log_fak(2),    0.693147, EPS);
129        TEST_EXPECT_SIMILAR(GB_log_fak(3),    1.791759, EPS);
130        TEST_EXPECT_SIMILAR(GB_log_fak(10),  15.104413, EPS);
131        TEST_EXPECT_SIMILAR(GB_log_fak(30),  74.658236, EPS);
132        TEST_EXPECT_SIMILAR(GB_log_fak(70), 230.439044, EPS);
133        TEST_EXPECT_SIMILAR(GB_log_fak(99), 359.134205, EPS);
134    }
135}
136
137#endif // UNIT_TESTS
138
139// --------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.