- メインページ
- モジュール
- ネームスペース
- データ構造
- ファイル

- ファイル一覧
- グローバル

random.h
説明を見る。
1 //
2 // Copyright (c) 2003-2011, MIST Project, Nagoya University
3 // All rights reserved.
4 //
5 // Redistribution and use in source and binary forms, with or without modification,
6 // are permitted provided that the following conditions are met:
7 //
8 // 1. Redistributions of source code must retain the above copyright notice,
9 // this list of conditions and the following disclaimer.
10 //
11 // 2. Redistributions in binary form must reproduce the above copyright notice,
12 // this list of conditions and the following disclaimer in the documentation
13 // and/or other materials provided with the distribution.
14 //
15 // 3. Neither the name of the Nagoya University nor the names of its contributors
16 // may be used to endorse or promote products derived from this software
17 // without specific prior written permission.
18 //
19 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
20 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
21 // FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
22 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
24 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
25 // IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
26 // THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 //
28
33 #ifndef __INCLUDE_MIST_RANDOM__
34 #define __INCLUDE_MIST_RANDOM__
35
36 #include <cmath>
37
38 #ifndef __INCLUDE_MIST_H__
39 #include "mist.h"
40 #endif
41
42 #ifndef __INCLUDE_MIST_MATRIX__
43 #include "matrix.h"
44 #endif
45
46 // mist名前空間の始まり
47 _MIST_BEGIN
48
49
57
58
59 namespace __random__
60 {
61 static const unsigned long n = 624;
62 static const unsigned long m = 397;
63 static const unsigned long vec_a = 0x9908b0dfUL;
64 static const unsigned long upper = 0x80000000UL;
65 static const unsigned long lower = 0x7fffffffUL;
66 static const double pai_timed_by_2 = 6.283185307179586;
67 }
68
69
71 namespace uniform
72 {
73
84 class random
85 {
86 // Period parameters
87
88 // static unsigned long vec_[n_]; /* the array for the state vector */
89 // static int i_=n_+1; /* i_==n_+1 means vec_[n_] is not initialized */
90
91 array< unsigned long > vec_;
92 unsigned long i_;
93
94
95 public:
100 random( ) : vec_( __random__::n ), i_( __random__::n + 1 )
101 {
102 }
103
104
111 explicit random( const unsigned long seed ) : vec_( __random__::n ), i_( __random__::n + 1 )
112 {
113 init( seed );
114 }
115
116
123 explicit random( const array< unsigned long >& seed_array ) : vec_( __random__::n ), i_( __random__::n + 1 )
124 {
125 init( seed_array );
126 }
127
128
135 void init( const unsigned long& seed )
136 {
137 vec_[ 0 ] = seed & 0xffffffffUL;
138
139 for( i_ = 1 ; i_ < __random__::n ; i_ ++ )
140 {
141 vec_[ i_ ] = ( 1812433253UL * ( vec_[ i_ - 1 ] ^ ( vec_[ i_ - 1 ] >> 30 ) ) + i_ );
142 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
143 /* In the previous versions, MSBs of the seed affect */
144 /* only MSBs of the array vec_[ ]. */
145 /* 2002/01/09 modified by Makoto Matsumoto */
146 vec_[ i_ ] &= 0xffffffffUL;
147 /* for >32 bit machines */
148 }
149 }
150
151 //
160 void init( const array< unsigned long >& seed_array )
161 {
162 unsigned long i, j, k;
163 init( 19650218UL );
164 i = 1;
165 j = 0;
166
167 for( k = ( __random__::n > seed_array.size( ) ? __random__::n : static_cast< unsigned long >( seed_array.size( ) ) ) ; k > 0 ; k -- )
168 {
169 vec_[ i ] = ( vec_[ i ] ^ ( ( vec_[ i - 1 ] ^ ( vec_[ i - 1 ] >> 30 ) ) * 1664525UL ) ) + seed_array[ j ] + j; /* non linear */
170 vec_[ i ] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
171 i ++;
172 j ++;
173
174 if( i >= __random__::n )
175 {
176 vec_[ 0 ] = vec_[ __random__::n - 1 ];
177 i = 1;
178 }
179
180 if( j >= seed_array.size( ) )
181 {
182 j = 0;
183 }
184 }
185
186 for( k = __random__::n - 1 ; k > 0 ; k -- )
187 {
188 vec_[ i ] = ( vec_[ i ] ^ ( ( vec_[ i - 1 ] ^ ( vec_[ i - 1 ] >> 30 ) ) * 1566083941UL ) ) - i; /* non linear */
189 vec_[ i ] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
190 i++;
191
192 if( i >= __random__::n )
193 {
194 vec_[ 0 ] = vec_[ __random__::n - 1 ];
195 i = 1;
196 }
197 }
198
199 vec_[ 0 ] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
200 }
201
208 unsigned long int32( )
209 {
210 unsigned long y;
211 static unsigned long mag01[ 2 ] = { 0x0UL, __random__::vec_a };
212 /* mag01[x] = x * __random__::vec_a for x=0,1 */
213
214 if( i_ >= __random__::n ) /* generate __random__::n words at one time */
215 {
216 unsigned long kk;
217
218 if( i_ == __random__::n + 1 ) /* if init_genrand() has not been called, */
219 {
220 init( 5489UL ); /* a default initial seed is used */
221 }
222 for( kk = 0 ; kk < __random__::n - __random__::m ; kk ++ )
223 {
224 y = ( vec_[ kk ] & __random__::upper ) | ( vec_[ kk + 1 ] & __random__::lower );
225 vec_[ kk ] = vec_[ kk + __random__::m ] ^ ( y >> 1 ) ^ mag01[ y & 0x1UL ];
226 }
227
228 for( ; kk < __random__::n - 1 ; kk ++ )
229 {
230 y = ( vec_[ kk ] & __random__::upper ) | ( vec_[ kk + 1 ] & __random__::lower );
231 vec_[ kk ] = vec_[ kk + ( __random__::m - __random__::n ) ] ^ ( y >> 1 ) ^ mag01[ y & 0x1UL ];
232 }
233
234 y = ( vec_[ __random__::n - 1 ] & __random__::upper ) | ( vec_[ 0 ] & __random__::lower );
235 vec_[ __random__::n - 1 ] = vec_[ __random__::m - 1 ] ^ ( y >> 1 ) ^ mag01[ y & 0x1UL ];
236
237 i_ = 0;
238 }
239
240 y = vec_[ i_ ++ ];
241
242 /* Tempering */
243 y ^= ( y >> 11 );
244 y ^= ( y << 7 ) & 0x9d2c5680UL;
245 y ^= ( y << 15 ) & 0xefc60000UL;
246 y ^= ( y >> 18 );
247
248 return y;
249 }
250
251
258 unsigned long int31( )
259 {
260 return ( int32( ) >> 1 );
261 }
262
263
270 double real1( )
271 {
272 return ( int32( ) * ( 1.0 / 4294967295.0 ) );
273 /* divided by 2^32-1 */
274 }
275
276
283 double real2( )
284 {
285 return ( int32( ) * ( 1.0 / 4294967296.0 ) );
286 /* divided by 2^32 */
287 }
288
289
296 double real3( )
297 {
298 return ( ( ( double ) int32( ) ) + 0.5 ) * ( 1.0 / 4294967296.0 );
299 /* divided by 2^32 */
300 }
301
302
309 double res53( )
310 {
311 const unsigned long a = int32( ) >> 5;
312 const unsigned long b = int32( ) >> 6;
313
314 return ( ( a * 67108864.0 + b ) * ( 1.0 / 9007199254740992.0 ) );
315 }
316
317
324 unsigned long generate( const unsigned int n = 0xffffffffUL )
325 {
326 return ( static_cast< unsigned long >( real2( ) * n ) );
327 }
328
329
336 unsigned long operator( )( const unsigned int n = 0xffffffffUL )
337 {
338 return generate( n );
339 }
340
341 };
342
343 } // uniform
344
345
346
348 namespace gauss
349 {
354 class random
355 {
356 uniform::random u_rand_;
357
358 double mean_;
359 double standard_deviation_;
360
361 public:
365 random( ) :
366 u_rand_( ),
367 mean_( 0.0 ),
368 standard_deviation_( 1.0 )
369 {
370 }
371
377 random( const double mean, const double variance ) :
378 u_rand_( ),
379 mean_( mean ),
380 standard_deviation_( std::sqrt( variance ) )
381 {
382 }
383
390 random( const unsigned long seed, const double mean = 0.0, const double variance = 1.0 ) :
391 u_rand_( seed ),
392 mean_( mean ),
393 standard_deviation_( std::sqrt( variance ) )
394 {
395 }
396
403 random( const array< unsigned long >& seed_array, const double mean = 0.0, const double variance = 1.0 ) :
404 u_rand_( seed_array ),
405 mean_( mean ),
406 standard_deviation_( std::sqrt( variance ) )
407 {
408 }
409
416 void init( const unsigned long& seed )
417 {
418 u_rand_.init( seed );
419 }
420
421 //
430 void init( const array< unsigned long >& seed_array )
431 {
432 u_rand_.init( seed_array );
433 }
434
440 void set_param( const double& mean = 0.0, const double& variance = 1.0 )
441 {
442 mean_ = mean;
443 standard_deviation_ = std::sqrt( variance );
444 }
445
450 double generate( )
451 {
452 return( standard_deviation_ * std::sqrt( -2.0 * std::log( 1.0 - u_rand_.real2( ) ) ) * std::cos( __random__::pai_timed_by_2 * ( 1.0 - u_rand_.real2( ) ) ) + mean_ );
453 }
454
459 double operator( )( )
460 {
461 return generate( );
462 }
463
464 };
465
466 } // gauss
467
468
469
471 namespace poisson
472 {
477 class random
478 {
479 uniform::random u_rand_;
480
481 double lambda_exp_;
482
483 public:
487 random( ) :
488 u_rand_( ),
489 lambda_exp_( std::exp( 1.0 ) )
490 {
491 }
492
497 random( const double lambda ) :
498 u_rand_( ),
499 lambda_exp_( std::exp( lambda ) )
500 {
501 }
502
508 random( const unsigned long seed, const double lambda = 1.0 ) :
509 u_rand_( seed ),
510 lambda_exp_( std::exp( lambda ) )
511 {
512 }
513
519 random( const array< unsigned long >& seed_array, const double lambda = 1.0 ) :
520 u_rand_( seed_array ),
521 lambda_exp_( std::exp( lambda ) )
522 {
523 }
524
531 void init( const unsigned long& seed )
532 {
533 u_rand_.init( seed );
534 }
535
536 //
545 void init( const array< unsigned long >& seed_array )
546 {
547 u_rand_.init( seed_array );
548 }
549
554 void set_param( const double lambda = 1.0 )
555 {
556 lambda_exp_ = std::exp( lambda );
557 }
558
563 int generate( )
564 {
565 int k = 0;
566 double lambda = lambda_exp_ * u_rand_.real2( );
567 for( ; lambda > 1.0 ; k++ )
568 {
569 lambda *= u_rand_.real2( );
570 }
571 return( k );
572 }
573
580 int generate( double l )
581 {
582 int k = 0;
583 double lambda = std::exp( l ) * u_rand_.real2( );
584 for( ; lambda > 1.0 ; k++ )
585 {
586 lambda *= u_rand_.real2( );
587 }
588 return( k );
589 }
590
595 int operator( )( )
596 {
597 return( generate( ) );
598 }
599
606 int operator( )( double l )
607 {
608 return( generate( l ) );
609 }
610 };
611
612 } // poisson
613
614
616 namespace multivariate_gauss
617 {
622 class random
623 {
624 private:
625 gauss::random g_rand_;
626
627 matrix< double > mean_;
628 matrix< double > l_triangle_;
629
630 const matrix< double > choleski( const matrix< double > &mat1 )
631 {
632 // coded by h.ishida
633 matrix< double > mat2( mat1.rows( ), mat1.cols( ) );
634 size_t i, j, m, k = mat1.cols( );
635
636 mat2( 0, 0 ) = std::sqrt( mat1( 0, 0 ) );
637
638 for( i = 0 ; i < k ; i ++ )
639 {
640 mat2( i, 0 ) = mat1( i, 0 ) / mat2( 0, 0 );
641 }
642
643 for( i = 1 ; i < k ; i ++ )
644 {
645 double s = 0.0;
646 for( j = 0 ; j < i ; j ++ )
647 {
648 s += mat2( i, j ) * mat2( i, j );
649 }
650
651 mat2( i, i ) = sqrt( mat1( i, i ) - s );
652 for( j = 1 ; j < i ; j ++ )
653 {
654 s = 0.0;
655 for( m = 0 ; m < j ; m ++ )
656 {
657 s += mat2( i, m ) * mat2( j, m );
658 }
659 mat2( i, j ) = ( mat1( i, j ) - s ) / mat2( j, j );
660 mat2( j, i ) = 0.0;
661 }
662 }
663
664 return( mat2 );
665 }
666
667 public:
671 random( ) :
672 g_rand_( ),
673 mean_( matrix< double >( ) ),
674 l_triangle_( matrix< double >( ) )
675 {
676 }
677
683 random( const matrix< double > &mean, const matrix< double > &covariance ) :
684 g_rand_( ),
685 mean_( mean ),
686 l_triangle_( choleski( covariance ) )
687 {
688 }
689
696 random( const unsigned long seed, const matrix< double > &mean, const matrix< double > &covariance ) :
697 g_rand_( seed ),
698 mean_( mean ),
699 l_triangle_( choleski( covariance ) )
700 {
701 }
702
709 random( const array< unsigned long >& seed_array, const matrix< double > &mean, const matrix< double > &covariance ) :
710 g_rand_( seed_array ),
711 mean_( mean ),
712 l_triangle_( choleski( covariance ) )
713 {
714 }
715
722 void init( const unsigned long& seed )
723 {
724 g_rand_.init( seed );
725 }
726
732 void set_param( const matrix< double > &mean, const matrix< double > &covariance )
733 {
734 mean_ = mean,
735 l_triangle_ = choleski( covariance );
736 }
737
742 matrix< double > generate( )
743 {
744 matrix< double > r_vec( mean_.rows( ), mean_.cols( ) );
745 for( size_t i = 0 ; i < r_vec.size( ) ; i ++ )
746 {
747 r_vec[ i ] = g_rand_.generate( );
748 }
749 return( mean_ + l_triangle_ * r_vec );
750 }
751
756 matrix< double > operator( )( )
757 {
758 return generate( );
759 }
760 };
761
762
763 } // multivariate_gauss
764
766 // 擬似乱数の生成の終わり
767
768
769 // mist名前空間の終わり
770 _MIST_END
771
772 /*
773 A C-program for MT19937, with initialization improved 2002/1/26.
774 Coded by Takuji Nishimura and Makoto Matsumoto.
775
776 Before using, initialize the state by using init_genrand(seed)
777 or init_by_array(init_key, key_length).
778
779 Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
780 All rights reserved.
781
782 Redistribution and use in source and binary forms, with or without
783 modification, are permitted provided that the following conditions
784 are met:
785
786 1. Redistributions of source code must retain the above copyright
787 notice, this list of conditions and the following disclaimer.
788
789 2. Redistributions in binary form must reproduce the above copyright
790 notice, this list of conditions and the following disclaimer in the
791 documentation and/or other materials provided with the distribution.
792
793 3. The names of its contributors may not be used to endorse or promote
794 products derived from this software without specific prior written
795 permission.
796
797 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
798 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
799 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
800 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
801 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
802 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
803 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
804 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
805 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
806 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
807 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
808 */
809
810 #endif // __INCLUDE_MIST_RANDOM__
811
812