random.cc 2.78 KB
Newer Older
1
// -*- coding: utf-8 -*-
2
// Copyright (C) 2011, 2012, 2013, 2014, 2015 Laboratoire de Recherche et
3
// Développement de l'Epita (LRDE).
4
5
// Copyright (C) 2004 Laboratoire d'Informatique de Paris 6 (LIP6),
// département Systèmes Répartis Coopératifs (SRC), Université Pierre
6
7
8
9
10
11
// et Marie Curie.
//
// This file is part of Spot, a model checking library.
//
// Spot is free software; you can redistribute it and/or modify it
// under the terms of the GNU General Public License as published by
12
// the Free Software Foundation; either version 3 of the License, or
13
14
15
16
17
18
19
20
// (at your option) any later version.
//
// Spot is distributed in the hope that it will be useful, but WITHOUT
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
// or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
// License for more details.
//
// You should have received a copy of the GNU General Public License
21
// along with this program.  If not, see <http://www.gnu.org/licenses/>.
22

23
#include "config.h"
24
#include "random.hh"
25
#include <random>
26
27
28

namespace spot
{
29
  static std::mt19937 gen;
30

31
32
33
  void
  srand(unsigned int seed)
  {
34
    gen.seed(seed);
35
36
37
38
39
  }

  double
  drand()
  {
40
    return  gen() / (1.0 + gen.max());
41
42
43
44
45
  }

  int
  mrand(int max)
  {
46
    return static_cast<int>(max * drand());
47
48
49
50
51
  }

  int
  rrand(int min, int max)
  {
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
    return min + static_cast<int>((max - min + 1) * drand());
  }

  double
  nrand()
  {
    const double r = drand();

    const double lim = 1.e-20;
    if (r < lim)
      return -1./lim;
    if (r > 1.0 - lim)
      return 1./lim;

    double t;
    if (r < 0.5)
      t = sqrt(-2.0 * log(r));
    else
      t = sqrt(-2.0 * log(1.0 - r));

    const double p0 = 0.322232431088;
    const double p1 = 1.0;
    const double p2 = 0.342242088547;
    const double p3 = 0.204231210245e-1;
    const double p4 = 0.453642210148e-4;
    const double q0 = 0.099348462606;
    const double q1 = 0.588581570495;
    const double q2 = 0.531103462366;
    const double q3 = 0.103537752850;
    const double q4 = 0.385607006340e-2;
    const double p = p0 + t * (p1 + t * (p2 + t * (p3 + t * p4)));
    const double q = q0 + t * (q1 + t * (q2 + t * (q3 + t * q4)));

    if (r < 0.5)
      return (p / q) - t;
    else
      return t - (p / q);
  }

  double
  bmrand()
  {
    static double next;
    static bool has_next = false;

    if (has_next)
      {
	has_next = false;
	return next;
      }

    double x;
    double y;
    double r;
    do
      {
	x = 2.0 * drand() - 1.0;
	y = 2.0 * drand() - 1.0;
	r = x * x + y * y;
      }
    while (r >= 1.0 || r == 0.0);
    r = sqrt(-2 * log(r) / r);
    next = y * r;
    has_next = true;
    return x * r;
117
118
119
  }

  int
120
  prand(double p)
121
  {
122
123
124
125
126
127
128
129
130
    double s = 0.0;
    long x = 0;

    while (s < p)
      {
	s -= log(1.0 - drand());
	++x;
      }
    return x - 1;
131
  }
132
}