Context Navigation

source: XMLIO_V2/external/include/random/uniform.h @ 80

Last change on this file since 80 was 80, checked in by ymipsl, 14 years ago
ajout lib externe
Property svn:eol-style set to `native`
File size: 11.3 KB

Line
1	#ifndef BZ_RANDOM_UNIFORM_H
2	#define BZ_RANDOM_UNIFORM_H
3
4	#include <random/default.h>
5
6	#ifndef FLT_MANT_DIG
7	#include <float.h>
8	#endif
9
10	BZ_NAMESPACE(ranlib)
11
12	/*****************************************************************************
13	* UniformClosedOpen generator: uniform random numbers in [0,1).
14	*****************************************************************************/
15
16	template<typename T = defaultType, typename IRNG = defaultIRNG,
17	typename stateTag = defaultState>
18	class UniformClosedOpen { };
19
20	// These constants are 1/2^32, 1/2^64, 1/2^96, 1/2^128
21	const long double norm32open = .2328306436538696289062500000000000000000E-9L;
22	const long double norm64open = .5421010862427522170037264004349708557129E-19L;
23	const long double norm96open = .1262177448353618888658765704452457967477E-28L;
24	const long double norm128open = .2938735877055718769921841343055614194547E-38L;
25
26
27	template<typename IRNG, typename stateTag>
28	class UniformClosedOpen<float,IRNG,stateTag>
29	: public IRNGWrapper<IRNG,stateTag>
30	{
31
32	public:
33	typedef float T_numtype;
34
35	float random()
36	{
37	#if FLT_MANT_DIG > 96
38	#if (FLT_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS)
39	#error RNG code assumes float mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning.
40	#endif
41	IRNG_int i1 = this->irng_.random();
42	IRNG_int i2 = this->irng_.random();
43	IRNG_int i3 = this->irng_.random();
44	IRNG_int i4 = this->irng_.random();
45	return i1 * norm128open + i2 * norm96open + i3 * norm64open
46	+ i4 * norm32open;
47	#elif FLT_MANT_DIG > 64
48	IRNG_int i1 = this->irng_.random();
49	IRNG_int i2 = this->irng_.random();
50	IRNG_int i3 = this->irng_.random();
51	return i1 * norm96open + i2 * norm64open + i3 * norm32open;
52	#elif FLT_MANT_DIG > 32
53	IRNG_int i1 = this->irng_.random();
54	IRNG_int i2 = this->irng_.random();
55	return i1 * norm64open + i2 * norm32open;
56	#else
57	IRNG_int i1 = this->irng_.random();
58	return i1 * norm32open;
59	#endif
60	}
61
62	float getUniform()
63	{ return random(); }
64	};
65
66	template<typename IRNG, typename stateTag>
67	class UniformClosedOpen<double,IRNG,stateTag>
68	: public IRNGWrapper<IRNG,stateTag>
69	{
70	public:
71	typedef double T_numtype;
72
73	double random()
74	{
75	#if DBL_MANT_DIG > 96
76	#if (DBL_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS)
77	#error RNG code assumes double mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning.
78	#endif
79
80	IRNG_int i1 = this->irng_.random();
81	IRNG_int i2 = this->irng_.random();
82	IRNG_int i3 = this->irng_.random();
83	IRNG_int i4 = this->irng_.random();
84	return i1 * norm128open + i2 * norm96open + i3 * norm64open
85	+ i4 * norm32open;
86	#elif DBL_MANT_DIG > 64
87	IRNG_int i1 = this->irng_.random();
88	IRNG_int i2 = this->irng_.random();
89	IRNG_int i3 = this->irng_.random();
90	return i1 * norm96open + i2 * norm64open + i3 * norm32open;
91	#elif DBL_MANT_DIG > 32
92	IRNG_int i1 = this->irng_.random();
93	IRNG_int i2 = this->irng_.random();
94	return i1 * norm64open + i2 * norm32open;
95	#else
96	IRNG_int i1 = this->irng_.random();
97	return i1 * norm32open;
98	#endif
99
100	}
101
102	double getUniform() { return random(); }
103	};
104
105	template<typename IRNG, typename stateTag>
106	class UniformClosedOpen<long double,IRNG,stateTag>
107	: public IRNGWrapper<IRNG,stateTag>
108	{
109	public:
110	typedef long double T_numtype;
111
112	long double random()
113	{
114	#if LDBL_MANT_DIG > 96
115	#if (LDBL_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS)
116	#error RNG code assumes long double mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning.
117	#endif
118
119	IRNG_int i1 = this->irng_.random();
120	IRNG_int i2 = this->irng_.random();
121	IRNG_int i3 = this->irng_.random();
122	IRNG_int i4 = this->irng_.random();
123	return i1 * norm128open + i2 * norm96open + i3 * norm64open
124	+ i4 * norm32open;
125	#elif LDBL_MANT_DIG > 64
126	IRNG_int i1 = this->irng_.random();
127	IRNG_int i2 = this->irng_.random();
128	IRNG_int i3 = this->irng_.random();
129	return i1 * norm96open + i2 * norm64open + i3 * norm32open;
130	#elif LDBL_MANT_DIG > 32
131	IRNG_int i1 = this->irng_.random();
132	IRNG_int i2 = this->irng_.random();
133	return i1 * norm64open + i2 * norm32open;
134	#else
135	IRNG_int i1 = this->irng_.random();
136	return i1 * norm32open;
137	#endif
138	}
139
140	long double getUniform() { return random(); }
141	};
142
143	// For people who don't care about open or closed: just give them
144	// ClosedOpen (this is the fastest).
145
146	template<class T, typename IRNG = defaultIRNG,
147	typename stateTag = defaultState>
148	class Uniform : public UniformClosedOpen<T,IRNG,stateTag>
149	{ };
150
151	/*****************************************************************************
152	* UniformClosed generator: uniform random numbers in [0,1].
153	*****************************************************************************/
154
155	// This constant is 1/(2^32-1)
156	const long double norm32closed = .2328306437080797375431469961868475648078E-9L;
157
158	// These constants are 2^32/(2^64-1) and 1/(2^64-1), respectively.
159
160	const long double norm64closed1 =
161	.23283064365386962891887177448353618888727188481031E-9L;
162	const long double norm64closed2 =
163	.54210108624275221703311375920552804341370213034169E-19L;
164
165	// These constants are 2^64/(2^96-1), 2^32/(2^96-1), and 1/(2^96-1)
166	const long double norm96closed1 = .2328306436538696289062500000029387358771E-9L;
167	const long double norm96closed2 =
168	.5421010862427522170037264004418131333707E-19L;
169	const long double norm96closed3 =
170	.1262177448353618888658765704468388886588E-28L;
171
172	// These constants are 2^96/(2^128-1), 2^64/(2^128-1), 2^32/(2^128-1) and
173	// 1/(2^128-1).
174	const long double norm128clos1 = .2328306436538696289062500000000000000007E-9L;
175	const long double norm128clos2 = .5421010862427522170037264004349708557145E-19L;
176	const long double norm128clos3 = .1262177448353618888658765704452457967481E-28L;
177	const long double norm128clos4 = .2938735877055718769921841343055614194555E-38L;
178
179
180	template<typename T = double, typename IRNG = defaultIRNG,
181	typename stateTag = defaultState>
182	class UniformClosed { };
183
184	template<typename IRNG, typename stateTag>
185	class UniformClosed<float,IRNG,stateTag> : public IRNGWrapper<IRNG,stateTag> {
186
187	public:
188	typedef float T_numtype;
189
190	float random()
191	{
192	#if FLTMANT_DIG > 96
193	#if (FLT_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS)
194	#error RNG code assumes float mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning.
195	#endif
196	IRNG_int i1 = this->irng_.random();
197	IRNG_int i2 = this->irng_.random();
198	IRNG_int i3 = this->irng_.random();
199	IRNG_int i4 = this->irng_.random();
200
201	return i1 * norm128clos1 + i2 * norm128clos2
202	+ i3 * norm128clos3 + i4 * norm128clos4;
203	#elif FLT_MANT_DIG > 64
204	IRNG_int i1 = this->irng_.random();
205	IRNG_int i2 = this->irng_.random();
206	IRNG_int i3 = this->irng_.random();
207
208	return i1 * norm96closed1 + i2 * norm96closed2
209	+ i3 * norm96closed3;
210	#elif FLT_MANT_DIG > 32
211	IRNG_int i1 = this->irng_.random();
212	IRNG_int i2 = this->irng_.random();
213
214	return i1 * norm64closed1 + i2 * norm64closed2;
215	#else
216	IRNG_int i = this->irng_.random();
217	return i * norm32closed;
218	#endif
219
220	}
221
222	float getUniform()
223	{ return random(); }
224	};
225
226	template<typename IRNG, typename stateTag>
227	class UniformClosed<double,IRNG,stateTag> : public IRNGWrapper<IRNG,stateTag> {
228
229	public:
230	typedef double T_numtype;
231
232	double random()
233	{
234	#if DBL_MANT_DIG > 96
235	#if (DBL_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS)
236	#error RNG code assumes double mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning.
237	#endif
238	IRNG_int i1 = this->irng_.random();
239	IRNG_int i2 = this->irng_.random();
240	IRNG_int i3 = this->irng_.random();
241	IRNG_int i4 = this->irng_.random();
242
243	return i1 * norm128clos1 + i2 * norm128clos2
244	+ i3 * norm128clos3 + i4 * norm128clos4;
245	#elif DBL_MANT_DIG > 64
246	IRNG_int i1 = this->irng_.random();
247	IRNG_int i2 = this->irng_.random();
248	IRNG_int i3 = this->irng_.random();
249
250	return i1 * norm96closed1 + i2 * norm96closed2
251	+ i3 * norm96closed3;
252	#elif DBL_MANT_DIG > 32
253	IRNG_int i1 = this->irng_.random();
254	IRNG_int i2 = this->irng_.random();
255
256	return i1 * norm64closed1 + i2 * norm64closed2;
257	#else
258	IRNG_int i = this->irng_.random();
259	return i * norm32closed;
260	#endif
261
262	}
263
264	double getUniform()
265	{ return random(); }
266	};
267
268	template<typename IRNG, typename stateTag>
269	class UniformClosed<long double,IRNG,stateTag>
270	: public IRNGWrapper<IRNG,stateTag> {
271
272	public:
273	typedef long double T_numtype;
274
275	long double random()
276	{
277	#if LDBL_MANT_DIG > 96
278	#if (LDBL_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS)
279	#error RNG code assumes long double mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning.
280	#endif
281	IRNG_int i1 = this->irng_.random();
282	IRNG_int i2 = this->irng_.random();
283	IRNG_int i3 = this->irng_.random();
284	IRNG_int i4 = this->irng_.random();
285
286	return i1 * norm128clos1 + i2 * norm128clos2
287	+ i3 * norm128clos3 + i4 * norm128clos4;
288	#elif LDBL_MANT_DIG > 64
289	IRNG_int i1 = this->irng_.random();
290	IRNG_int i2 = this->irng_.random();
291	IRNG_int i3 = this->irng_.random();
292
293	return i1 * norm96closed1 + i2 * norm96closed2
294	+ i3 * norm96closed3;
295	#elif LDBL_MANT_DIG > 32
296	IRNG_int i1 = this->irng_.random();
297	IRNG_int i2 = this->irng_.random();
298
299	return i1 * norm64closed1 + i2 * norm64closed2;
300	#else
301	IRNG_int i = this->irng_.random();
302	return i * norm32closed;
303	#endif
304	}
305
306	long double getUniform()
307	{ return random(); }
308
309	};
310
311	/*****************************************************************************
312	* UniformOpen generator: uniform random numbers in (0,1).
313	*****************************************************************************/
314
315	template<typename T = double, typename IRNG = defaultIRNG,
316	typename stateTag = defaultState>
317	class UniformOpen : public UniformClosedOpen<T,IRNG,stateTag> {
318	public:
319	typedef T T_numtype;
320
321	T random()
322	{
323	// Turn a [0,1) into a (0,1) interval by weeding out
324	// any zeros.
325	T x;
326	do {
327	x = UniformClosedOpen<T,IRNG,stateTag>::random();
328	} while (x == 0.0L);
329
330	return x;
331	}
332
333	T getUniform()
334	{ return random(); }
335
336	};
337
338	/*****************************************************************************
339	* UniformOpenClosed generator: uniform random numbers in (0,1]
340	*****************************************************************************/
341
342	template<typename T = double, typename IRNG = defaultIRNG,
343	typename stateTag = defaultState>
344	class UniformOpenClosed : public UniformClosedOpen<T,IRNG,stateTag> {
345
346	public:
347	typedef T T_numtype;
348
349	T random()
350	{
351	// Antithetic value: taking 1-X where X is [0,1) results
352	// in a (0,1] distribution.
353	return 1.0 - UniformClosedOpen<T,IRNG,stateTag>::random();
354	}
355
356	T getUniform()
357	{ return random(); }
358	};
359
360	BZ_NAMESPACE_END
361
362	#endif // BZ_RANDOM_UNIFORM_H

Note: See TracBrowser for help on using the repository browser.

Download in other formats: