1    	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2    	/*                                                                           */
3    	/*                  This file is part of the class library                   */
4    	/*       SoPlex --- the Sequential object-oriented simPlex.                  */
5    	/*                                                                           */
6    	/*  Copyright (c) 1996-2023 Zuse Institute Berlin (ZIB)                      */
7    	/*                                                                           */
8    	/*  Licensed under the Apache License, Version 2.0 (the "License");          */
9    	/*  you may not use this file except in compliance with the License.         */
10   	/*  You may obtain a copy of the License at                                  */
11   	/*                                                                           */
12   	/*      http://www.apache.org/licenses/LICENSE-2.0                           */
13   	/*                                                                           */
14   	/*  Unless required by applicable law or agreed to in writing, software      */
15   	/*  distributed under the License is distributed on an "AS IS" BASIS,        */
16   	/*  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
17   	/*  See the License for the specific language governing permissions and      */
18   	/*  limitations under the License.                                           */
19   	/*                                                                           */
20   	/*  You should have received a copy of the Apache-2.0 license                */
21   	/*  along with SoPlex; see the file LICENSE. If not email to soplex@zib.de.  */
22   	/*                                                                           */
23   	/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24   	
25   	/**@file  random.h
26   	 * @brief Random numbers.
27   	 */
28   	#ifndef _RANDOM_H_
29   	#define _RANDOM_H_
30   	
31   	#include <assert.h>
32   	#include <stdint.h>
33   	#include <algorithm>
34   	
35   	namespace soplex
36   	{
37   	
38   	/* initial seeds for KISS random number generator */
39   	#define SOPLEX_DEFAULT_LIN  UINT32_C(123456789)
40   	#define SOPLEX_DEFAULT_XOR  UINT32_C(362436000)
41   	#define SOPLEX_DEFAULT_MWC  UINT32_C(521288629)
42   	#define SOPLEX_DEFAULT_CST  UINT32_C(7654321)
43   	
44   	/* defines for linear congruential generator */
45   	#define SOPLEX_RSTEP    UINT64_C(1103515245)
46   	#define SOPLEX_RADD     UINT64_C(12345)
47   	
48   	/**@brief   Random numbers.
49   	   @ingroup Elementary
50   	
51   	   Class Random provides random Real variables, i.e. a value variable that
52   	   gives another value each time it is accessed. It may be used just like an
53   	   ordinary Real by means of an overloaded cast operator Real()%.
54   	
55   	   This is an implementation of KISS random number generator developed by George Marsaglia.
56   	   KISS is combination of three different random number generators:
57   	    - Linear congruential generator
58   	    - Xorshift
59   	    - Lag-1 Multiply-with-carry
60   	
61   	   KISS has a period of 2^123 and passes all statistical test part of BigCrush-Test of TestU01 [1].
62   	
63   	   [1] http://dl.acm.org/citation.cfm?doid=1268776.1268777
64   	*/
65   	class Random
66   	{
67   	private:
68   	
69   	   //--------------------------------------
70   	   /**@name Data */
71   	   ///@{
72   	   uint32_t seedshift;  ///< initial shift for random seeds.
73   	   uint32_t lin_seed;   ///< random seed for linear congruential RNS.
74   	   uint32_t xor_seed;   ///< random seed for XOR-shift RNS.
75   	   uint32_t mwc_seed;   ///< random seed Multiple-with-carry RNS.
76   	   uint32_t cst_seed;   ///< random seed shifted for mwc_seed.
77   	   ///@}
78   	
79   	   //--------------------------------------
80   	   /**@name Helpers */
81   	   ///@{
82   	   /// executes KISS random number generator and returns a pseudo random Real value in [0,1].
83   	   Real next_random()
84   	   {
85   	      uint64_t t;
86   	
87   	      /* linear congruential */
88   	      lin_seed = (uint32_t)(lin_seed * SOPLEX_RSTEP + SOPLEX_RADD);
89   	
90   	      /* Xorshift */
91   	      xor_seed ^= (xor_seed << 13);
92   	      xor_seed ^= (xor_seed >> 17);
93   	      xor_seed ^= (xor_seed << 5);
94   	
95   	      /* Multiply-with-carry */
96   	      t = UINT64_C(698769069) * mwc_seed + cst_seed;
97   	      cst_seed = (uint32_t)(t >> 32);
98   	      mwc_seed = (uint32_t) t;
99   	
100  	      return (lin_seed + xor_seed + mwc_seed) / (Real)UINT32_MAX;
101  	   }
102  	
103  	   ///@}
104  	
105  	public:
106  	
107  	   //--------------------------------------
108  	   /**@name Access */
109  	   ///@{
110  	   /// returns next random number.
111  	   Real next(Real minimum = 0.0, Real maximum = 1.0)
112  	   {
113  	      Real randnumber = next_random();
114  	
115  	      /* we multiply minimum and maximum separately by randnumber in order to avoid overflow if they are more than
116  	       * std::numeric_limits<Real>::max() apart
117  	       */
118  	      return minimum * (1.0 - randnumber) + maximum * randnumber;
119  	   }
120  	
121  	   /// returns the initial seed shift
122  	   uint32_t getSeed() const
123  	   {
124  	      return seedshift;
125  	   }
126  	
127  	   ///@}
128  	
129  	   //--------------------------------------
130  	   /**@name Modification */
131  	   ///@{
132  	   /// initialize all seeds of the random number generator.
133  	   void setSeed(uint32_t initshift)
134  	   {
135  	      seedshift = initshift;
136  	
137  	      /* use SOPLEX_MAX to avoid zero after over flowing */
138  	      lin_seed = SOPLEX_MAX(SOPLEX_DEFAULT_LIN + initshift, 1u);
139  	      xor_seed = SOPLEX_MAX(SOPLEX_DEFAULT_XOR + initshift, 1u);
140  	      mwc_seed = SOPLEX_MAX(SOPLEX_DEFAULT_MWC + initshift, 1u);
141  	      cst_seed = SOPLEX_DEFAULT_CST + initshift;
142  	
143  	      assert(lin_seed > 0);
144  	      assert(xor_seed > 0);
145  	      assert(mwc_seed > 0);
146  	
147  	      /* advance state once to have more random values */
148  	      (void) next_random();
149  	   }
150  	
151  	   ///@}
152  	
153  	
154  	   //--------------------------------------
155  	   /**@name Constructors / destructors */
156  	   ///@{
157  	   /// default constructor.
158  	   /** Constructs a new (pseudo) Random variable using \p randomseed as seed for the random
159  	       variable's sequence.
160  	   */
161  	   explicit
162  	   Random(uint32_t randomseed = 0)
163  	   {
164  	      setSeed(randomseed);
165  	   }
166  	
167  	   /// destructor
168  	   ~Random()
169  	   {}
170  	   ///@}
171  	};
172  	
173  	} // namespace soplex
174  	#endif // _RANDOM_H_
175