34 * |
34 * |
35 * 2. Redistributions in binary form must reproduce the above copyright |
35 * 2. Redistributions in binary form must reproduce the above copyright |
36 * notice, this list of conditions and the following disclaimer in the |
36 * notice, this list of conditions and the following disclaimer in the |
37 * documentation and/or other materials provided with the distribution. |
37 * documentation and/or other materials provided with the distribution. |
38 * |
38 * |
39 * 3. The names of its contributors may not be used to endorse or promote |
39 * 3. The names of its contributors may not be used to endorse or promote |
40 * products derived from this software without specific prior written |
40 * products derived from this software without specific prior written |
41 * permission. |
41 * permission. |
42 * |
42 * |
43 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
43 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
44 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |
44 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |
45 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS |
45 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS |
165 |
165 |
166 void initState() { |
166 void initState() { |
167 static const Word seedArray[4] = { |
167 static const Word seedArray[4] = { |
168 0x12345u, 0x23456u, 0x34567u, 0x45678u |
168 0x12345u, 0x23456u, 0x34567u, 0x45678u |
169 }; |
169 }; |
170 |
170 |
171 initState(seedArray, seedArray + 4); |
171 initState(seedArray, seedArray + 4); |
172 } |
172 } |
173 |
173 |
174 void initState(Word seed) { |
174 void initState(Word seed) { |
175 |
175 |
176 static const Word mul = RandomTraits<Word>::mul; |
176 static const Word mul = RandomTraits<Word>::mul; |
177 |
177 |
178 current = state; |
178 current = state; |
179 |
179 |
180 Word *curr = state + length - 1; |
180 Word *curr = state + length - 1; |
181 curr[0] = seed; --curr; |
181 curr[0] = seed; --curr; |
182 for (int i = 1; i < length; ++i) { |
182 for (int i = 1; i < length; ++i) { |
183 curr[0] = (mul * ( curr[1] ^ (curr[1] >> (bits - 2)) ) + i); |
183 curr[0] = (mul * ( curr[1] ^ (curr[1] >> (bits - 2)) ) + i); |
267 state[0] = (((state[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^ |
267 state[0] = (((state[0] & hiMask) | (curr[length - 1] & loMask)) >> 1) ^ |
268 curr[length - shift] ^ mask[curr[length - 1] & 1ul]; |
268 curr[length - shift] ^ mask[curr[length - 1] & 1ul]; |
269 |
269 |
270 } |
270 } |
271 |
271 |
272 |
272 |
273 Word *current; |
273 Word *current; |
274 Word state[length]; |
274 Word state[length]; |
275 |
275 |
276 }; |
276 }; |
277 |
277 |
278 |
278 |
279 template <typename Result, |
279 template <typename Result, |
280 int shift = (std::numeric_limits<Result>::digits + 1) / 2> |
280 int shift = (std::numeric_limits<Result>::digits + 1) / 2> |
281 struct Masker { |
281 struct Masker { |
282 static Result mask(const Result& result) { |
282 static Result mask(const Result& result) { |
283 return Masker<Result, (shift + 1) / 2>:: |
283 return Masker<Result, (shift + 1) / 2>:: |
284 mask(static_cast<Result>(result | (result >> shift))); |
284 mask(static_cast<Result>(result | (result >> shift))); |
285 } |
285 } |
286 }; |
286 }; |
287 |
287 |
288 template <typename Result> |
288 template <typename Result> |
289 struct Masker<Result, 1> { |
289 struct Masker<Result, 1> { |
290 static Result mask(const Result& result) { |
290 static Result mask(const Result& result) { |
291 return static_cast<Result>(result | (result >> 1)); |
291 return static_cast<Result>(result | (result >> 1)); |
292 } |
292 } |
293 }; |
293 }; |
294 |
294 |
295 template <typename Result, typename Word, |
295 template <typename Result, typename Word, |
296 int rest = std::numeric_limits<Result>::digits, int shift = 0, |
296 int rest = std::numeric_limits<Result>::digits, int shift = 0, |
297 bool last = rest <= std::numeric_limits<Word>::digits> |
297 bool last = rest <= std::numeric_limits<Word>::digits> |
298 struct IntConversion { |
298 struct IntConversion { |
299 static const int bits = std::numeric_limits<Word>::digits; |
299 static const int bits = std::numeric_limits<Word>::digits; |
300 |
300 |
301 static Result convert(RandomCore<Word>& rnd) { |
301 static Result convert(RandomCore<Word>& rnd) { |
302 return static_cast<Result>(rnd() >> (bits - rest)) << shift; |
302 return static_cast<Result>(rnd() >> (bits - rest)) << shift; |
303 } |
303 } |
304 |
304 |
305 }; |
305 }; |
306 |
306 |
307 template <typename Result, typename Word, int rest, int shift> |
307 template <typename Result, typename Word, int rest, int shift> |
308 struct IntConversion<Result, Word, rest, shift, false> { |
308 struct IntConversion<Result, Word, rest, shift, false> { |
309 static const int bits = std::numeric_limits<Word>::digits; |
309 static const int bits = std::numeric_limits<Word>::digits; |
310 |
310 |
311 static Result convert(RandomCore<Word>& rnd) { |
311 static Result convert(RandomCore<Word>& rnd) { |
312 return (static_cast<Result>(rnd()) << shift) | |
312 return (static_cast<Result>(rnd()) << shift) | |
313 IntConversion<Result, Word, rest - bits, shift + bits>::convert(rnd); |
313 IntConversion<Result, Word, rest - bits, shift + bits>::convert(rnd); |
314 } |
314 } |
315 }; |
315 }; |
316 |
316 |
317 |
317 |
318 template <typename Result, typename Word, |
318 template <typename Result, typename Word, |
319 bool one_word = (std::numeric_limits<Word>::digits < |
319 bool one_word = (std::numeric_limits<Word>::digits < |
320 std::numeric_limits<Result>::digits) > |
320 std::numeric_limits<Result>::digits) > |
321 struct Mapping { |
321 struct Mapping { |
322 static Result map(RandomCore<Word>& rnd, const Result& bound) { |
322 static Result map(RandomCore<Word>& rnd, const Result& bound) { |
323 Word max = Word(bound - 1); |
323 Word max = Word(bound - 1); |
324 Result mask = Masker<Result>::mask(bound - 1); |
324 Result mask = Masker<Result>::mask(bound - 1); |
325 Result num; |
325 Result num; |
326 do { |
326 do { |
327 num = IntConversion<Result, Word>::convert(rnd) & mask; |
327 num = IntConversion<Result, Word>::convert(rnd) & mask; |
328 } while (num > max); |
328 } while (num > max); |
329 return num; |
329 return num; |
330 } |
330 } |
331 }; |
331 }; |
332 |
332 |
348 struct ShiftMultiplier { |
348 struct ShiftMultiplier { |
349 static const Result multiplier() { |
349 static const Result multiplier() { |
350 Result res = ShiftMultiplier<Result, exp / 2>::multiplier(); |
350 Result res = ShiftMultiplier<Result, exp / 2>::multiplier(); |
351 res *= res; |
351 res *= res; |
352 if ((exp & 1) == 1) res *= static_cast<Result>(2.0); |
352 if ((exp & 1) == 1) res *= static_cast<Result>(2.0); |
353 return res; |
353 return res; |
354 } |
354 } |
355 }; |
355 }; |
356 |
356 |
357 template <typename Result, int exp> |
357 template <typename Result, int exp> |
358 struct ShiftMultiplier<Result, exp, false> { |
358 struct ShiftMultiplier<Result, exp, false> { |
359 static const Result multiplier() { |
359 static const Result multiplier() { |
360 Result res = ShiftMultiplier<Result, exp / 2>::multiplier(); |
360 Result res = ShiftMultiplier<Result, exp / 2>::multiplier(); |
361 res *= res; |
361 res *= res; |
362 if ((exp & 1) == 1) res *= static_cast<Result>(0.5); |
362 if ((exp & 1) == 1) res *= static_cast<Result>(0.5); |
363 return res; |
363 return res; |
364 } |
364 } |
365 }; |
365 }; |
366 |
366 |
367 template <typename Result> |
367 template <typename Result> |
368 struct ShiftMultiplier<Result, 0, true> { |
368 struct ShiftMultiplier<Result, 0, true> { |
369 static const Result multiplier() { |
369 static const Result multiplier() { |
370 return static_cast<Result>(1.0); |
370 return static_cast<Result>(1.0); |
371 } |
371 } |
372 }; |
372 }; |
373 |
373 |
374 template <typename Result> |
374 template <typename Result> |
375 struct ShiftMultiplier<Result, -20, true> { |
375 struct ShiftMultiplier<Result, -20, true> { |
376 static const Result multiplier() { |
376 static const Result multiplier() { |
377 return static_cast<Result>(1.0/1048576.0); |
377 return static_cast<Result>(1.0/1048576.0); |
378 } |
378 } |
379 }; |
379 }; |
380 |
380 |
381 template <typename Result> |
381 template <typename Result> |
382 struct ShiftMultiplier<Result, -32, true> { |
382 struct ShiftMultiplier<Result, -32, true> { |
383 static const Result multiplier() { |
383 static const Result multiplier() { |
384 return static_cast<Result>(1.0/424967296.0); |
384 return static_cast<Result>(1.0/424967296.0); |
385 } |
385 } |
386 }; |
386 }; |
387 |
387 |
388 template <typename Result> |
388 template <typename Result> |
389 struct ShiftMultiplier<Result, -53, true> { |
389 struct ShiftMultiplier<Result, -53, true> { |
390 static const Result multiplier() { |
390 static const Result multiplier() { |
391 return static_cast<Result>(1.0/9007199254740992.0); |
391 return static_cast<Result>(1.0/9007199254740992.0); |
392 } |
392 } |
393 }; |
393 }; |
394 |
394 |
395 template <typename Result> |
395 template <typename Result> |
396 struct ShiftMultiplier<Result, -64, true> { |
396 struct ShiftMultiplier<Result, -64, true> { |
397 static const Result multiplier() { |
397 static const Result multiplier() { |
398 return static_cast<Result>(1.0/18446744073709551616.0); |
398 return static_cast<Result>(1.0/18446744073709551616.0); |
399 } |
399 } |
400 }; |
400 }; |
401 |
401 |
402 template <typename Result, int exp> |
402 template <typename Result, int exp> |
403 struct Shifting { |
403 struct Shifting { |
405 return result * ShiftMultiplier<Result, exp>::multiplier(); |
405 return result * ShiftMultiplier<Result, exp>::multiplier(); |
406 } |
406 } |
407 }; |
407 }; |
408 |
408 |
409 template <typename Result, typename Word, |
409 template <typename Result, typename Word, |
410 int rest = std::numeric_limits<Result>::digits, int shift = 0, |
410 int rest = std::numeric_limits<Result>::digits, int shift = 0, |
411 bool last = rest <= std::numeric_limits<Word>::digits> |
411 bool last = rest <= std::numeric_limits<Word>::digits> |
412 struct RealConversion{ |
412 struct RealConversion{ |
413 static const int bits = std::numeric_limits<Word>::digits; |
413 static const int bits = std::numeric_limits<Word>::digits; |
414 |
414 |
415 static Result convert(RandomCore<Word>& rnd) { |
415 static Result convert(RandomCore<Word>& rnd) { |
416 return Shifting<Result, - shift - rest>:: |
416 return Shifting<Result, - shift - rest>:: |
417 shift(static_cast<Result>(rnd() >> (bits - rest))); |
417 shift(static_cast<Result>(rnd() >> (bits - rest))); |
418 } |
418 } |
419 }; |
419 }; |
420 |
420 |
421 template <typename Result, typename Word, int rest, int shift> |
421 template <typename Result, typename Word, int rest, int shift> |
422 struct RealConversion<Result, Word, rest, shift, false> { |
422 struct RealConversion<Result, Word, rest, shift, false> { |
423 static const int bits = std::numeric_limits<Word>::digits; |
423 static const int bits = std::numeric_limits<Word>::digits; |
424 |
424 |
425 static Result convert(RandomCore<Word>& rnd) { |
425 static Result convert(RandomCore<Word>& rnd) { |
426 return Shifting<Result, - shift - bits>:: |
426 return Shifting<Result, - shift - bits>:: |
427 shift(static_cast<Result>(rnd())) + |
427 shift(static_cast<Result>(rnd())) + |
552 /// \brief Constructor with seed |
552 /// \brief Constructor with seed |
553 /// |
553 /// |
554 /// Constructor with seed. The current number type will be converted |
554 /// Constructor with seed. The current number type will be converted |
555 /// to the architecture word type. |
555 /// to the architecture word type. |
556 template <typename Number> |
556 template <typename Number> |
557 Random(Number seed) { |
557 Random(Number seed) { |
558 _random_bits::Initializer<Number, Word>::init(core, seed); |
558 _random_bits::Initializer<Number, Word>::init(core, seed); |
559 } |
559 } |
560 |
560 |
561 /// \brief Constructor with array seeding |
561 /// \brief Constructor with array seeding |
562 /// |
562 /// |
563 /// Constructor with array seeding. The given range should contain |
563 /// Constructor with array seeding. The given range should contain |
564 /// any number type and the numbers will be converted to the |
564 /// any number type and the numbers will be converted to the |
565 /// architecture word type. |
565 /// architecture word type. |
566 template <typename Iterator> |
566 template <typename Iterator> |
567 Random(Iterator begin, Iterator end) { |
567 Random(Iterator begin, Iterator end) { |
568 typedef typename std::iterator_traits<Iterator>::value_type Number; |
568 typedef typename std::iterator_traits<Iterator>::value_type Number; |
569 _random_bits::Initializer<Number, Word>::init(core, begin, end); |
569 _random_bits::Initializer<Number, Word>::init(core, begin, end); |
570 } |
570 } |
571 |
571 |
572 /// \brief Copy constructor |
572 /// \brief Copy constructor |
595 /// \brief Seeding random sequence |
595 /// \brief Seeding random sequence |
596 /// |
596 /// |
597 /// Seeding the random sequence. The current number type will be |
597 /// Seeding the random sequence. The current number type will be |
598 /// converted to the architecture word type. |
598 /// converted to the architecture word type. |
599 template <typename Number> |
599 template <typename Number> |
600 void seed(Number seed) { |
600 void seed(Number seed) { |
601 _random_bits::Initializer<Number, Word>::init(core, seed); |
601 _random_bits::Initializer<Number, Word>::init(core, seed); |
602 } |
602 } |
603 |
603 |
604 /// \brief Seeding random sequence |
604 /// \brief Seeding random sequence |
605 /// |
605 /// |
606 /// Seeding the random sequence. The given range should contain |
606 /// Seeding the random sequence. The given range should contain |
607 /// any number type and the numbers will be converted to the |
607 /// any number type and the numbers will be converted to the |
608 /// architecture word type. |
608 /// architecture word type. |
609 template <typename Iterator> |
609 template <typename Iterator> |
610 void seed(Iterator begin, Iterator end) { |
610 void seed(Iterator begin, Iterator end) { |
611 typedef typename std::iterator_traits<Iterator>::value_type Number; |
611 typedef typename std::iterator_traits<Iterator>::value_type Number; |
612 _random_bits::Initializer<Number, Word>::init(core, begin, end); |
612 _random_bits::Initializer<Number, Word>::init(core, begin, end); |
613 } |
613 } |
614 |
614 |
615 /// \brief Seeding from file or from process id and time |
615 /// \brief Seeding from file or from process id and time |
694 |
694 |
695 /// \brief Returns a random real number the range [0, b) |
695 /// \brief Returns a random real number the range [0, b) |
696 /// |
696 /// |
697 /// It returns a random real number from the range [0, b). |
697 /// It returns a random real number from the range [0, b). |
698 template <typename Number> |
698 template <typename Number> |
699 Number real(Number b) { |
699 Number real(Number b) { |
700 return real<Number>() * b; |
700 return real<Number>() * b; |
701 } |
701 } |
702 |
702 |
703 /// \brief Returns a random real number from the range [a, b) |
703 /// \brief Returns a random real number from the range [a, b) |
704 /// |
704 /// |
705 /// It returns a random real number from the range [a, b). |
705 /// It returns a random real number from the range [a, b). |
706 template <typename Number> |
706 template <typename Number> |
707 Number real(Number a, Number b) { |
707 Number real(Number a, Number b) { |
708 return real<Number>() * (b - a) + a; |
708 return real<Number>() * (b - a) + a; |
709 } |
709 } |
710 |
710 |
711 /// @} |
711 /// @} |
712 |
712 |
713 ///\name Uniform distributions |
713 ///\name Uniform distributions |
723 |
723 |
724 /// \brief Returns a random real number from the range [0, b) |
724 /// \brief Returns a random real number from the range [0, b) |
725 /// |
725 /// |
726 /// It returns a random real number from the range [0, b). |
726 /// It returns a random real number from the range [0, b). |
727 template <typename Number> |
727 template <typename Number> |
728 Number operator()(Number b) { |
728 Number operator()(Number b) { |
729 return real<Number>() * b; |
729 return real<Number>() * b; |
730 } |
730 } |
731 |
731 |
732 /// \brief Returns a random real number from the range [a, b) |
732 /// \brief Returns a random real number from the range [a, b) |
733 /// |
733 /// |
734 /// It returns a random real number from the range [a, b). |
734 /// It returns a random real number from the range [a, b). |
735 template <typename Number> |
735 template <typename Number> |
736 Number operator()(Number a, Number b) { |
736 Number operator()(Number a, Number b) { |
737 return real<Number>() * (b - a) + a; |
737 return real<Number>() * (b - a) + a; |
738 } |
738 } |
739 |
739 |
740 /// \brief Returns a random integer from a range |
740 /// \brief Returns a random integer from a range |
741 /// |
741 /// |
742 /// It returns a random integer from the range {0, 1, ..., b - 1}. |
742 /// It returns a random integer from the range {0, 1, ..., b - 1}. |
782 /// It returns a random integer uniformly from the whole range of |
782 /// It returns a random integer uniformly from the whole range of |
783 /// the current \c Number type. The default result type of this |
783 /// the current \c Number type. The default result type of this |
784 /// function is \c int. |
784 /// function is \c int. |
785 template <typename Number> |
785 template <typename Number> |
786 Number integer() { |
786 Number integer() { |
787 static const int nb = std::numeric_limits<Number>::digits + |
787 static const int nb = std::numeric_limits<Number>::digits + |
788 (std::numeric_limits<Number>::is_signed ? 1 : 0); |
788 (std::numeric_limits<Number>::is_signed ? 1 : 0); |
789 return _random_bits::IntConversion<Number, Word, nb>::convert(core); |
789 return _random_bits::IntConversion<Number, Word, nb>::convert(core); |
790 } |
790 } |
791 |
791 |
792 int integer() { |
792 int integer() { |
793 return integer<int>(); |
793 return integer<int>(); |
794 } |
794 } |
795 |
795 |
796 /// \brief Returns a random bool |
796 /// \brief Returns a random bool |
797 /// |
797 /// |
798 /// It returns a random bool. The generator holds a buffer for |
798 /// It returns a random bool. The generator holds a buffer for |
799 /// random bits. Every time when it become empty the generator makes |
799 /// random bits. Every time when it become empty the generator makes |
800 /// a new random word and fill the buffer up. |
800 /// a new random word and fill the buffer up. |
852 } |
852 } |
853 |
853 |
854 /// Gamma distribution with given integer shape |
854 /// Gamma distribution with given integer shape |
855 |
855 |
856 /// This function generates a gamma distribution random number. |
856 /// This function generates a gamma distribution random number. |
857 /// |
857 /// |
858 ///\param k shape parameter (<tt>k>0</tt> integer) |
858 ///\param k shape parameter (<tt>k>0</tt> integer) |
859 double gamma(int k) |
859 double gamma(int k) |
860 { |
860 { |
861 double s = 0; |
861 double s = 0; |
862 for(int i=0;i<k;i++) s-=std::log(1.0-real<double>()); |
862 for(int i=0;i<k;i++) s-=std::log(1.0-real<double>()); |
863 return s; |
863 return s; |
864 } |
864 } |
865 |
865 |
866 /// Gamma distribution with given shape and scale parameter |
866 /// Gamma distribution with given shape and scale parameter |
867 |
867 |
868 /// This function generates a gamma distribution random number. |
868 /// This function generates a gamma distribution random number. |
869 /// |
869 /// |
870 ///\param k shape parameter (<tt>k>0</tt>) |
870 ///\param k shape parameter (<tt>k>0</tt>) |
871 ///\param theta scale parameter |
871 ///\param theta scale parameter |
872 /// |
872 /// |
873 double gamma(double k,double theta=1.0) |
873 double gamma(double k,double theta=1.0) |
874 { |
874 { |
875 double xi,nu; |
875 double xi,nu; |
876 const double delta = k-std::floor(k); |
876 const double delta = k-std::floor(k); |
877 const double v0=E/(E-delta); |
877 const double v0=E/(E-delta); |
878 do { |
878 do { |
879 double V0=1.0-real<double>(); |
879 double V0=1.0-real<double>(); |
880 double V1=1.0-real<double>(); |
880 double V1=1.0-real<double>(); |
881 double V2=1.0-real<double>(); |
881 double V2=1.0-real<double>(); |
882 if(V2<=v0) |
882 if(V2<=v0) |
883 { |
883 { |
884 xi=std::pow(V1,1.0/delta); |
884 xi=std::pow(V1,1.0/delta); |
885 nu=V0*std::pow(xi,delta-1.0); |
885 nu=V0*std::pow(xi,delta-1.0); |
886 } |
886 } |
887 else |
887 else |
888 { |
888 { |
889 xi=1.0-std::log(V1); |
889 xi=1.0-std::log(V1); |
890 nu=V0*std::exp(-xi); |
890 nu=V0*std::exp(-xi); |
891 } |
891 } |
892 } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi)); |
892 } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi)); |
893 return theta*(xi+gamma(int(std::floor(k)))); |
893 return theta*(xi+gamma(int(std::floor(k)))); |
894 } |
894 } |
895 |
895 |
896 /// Weibull distribution |
896 /// Weibull distribution |
897 |
897 |
898 /// This function generates a Weibull distribution random number. |
898 /// This function generates a Weibull distribution random number. |
899 /// |
899 /// |
900 ///\param k shape parameter (<tt>k>0</tt>) |
900 ///\param k shape parameter (<tt>k>0</tt>) |
901 ///\param lambda scale parameter (<tt>lambda>0</tt>) |
901 ///\param lambda scale parameter (<tt>lambda>0</tt>) |
902 /// |
902 /// |
903 double weibull(double k,double lambda) |
903 double weibull(double k,double lambda) |
904 { |
904 { |
905 return lambda*pow(-std::log(1.0-real<double>()),1.0/k); |
905 return lambda*pow(-std::log(1.0-real<double>()),1.0/k); |
906 } |
906 } |
907 |
907 |
908 /// Pareto distribution |
908 /// Pareto distribution |
909 |
909 |
910 /// This function generates a Pareto distribution random number. |
910 /// This function generates a Pareto distribution random number. |
911 /// |
911 /// |
912 ///\param k shape parameter (<tt>k>0</tt>) |
912 ///\param k shape parameter (<tt>k>0</tt>) |
913 ///\param x_min location parameter (<tt>x_min>0</tt>) |
913 ///\param x_min location parameter (<tt>x_min>0</tt>) |
914 /// |
914 /// |
915 double pareto(double k,double x_min) |
915 double pareto(double k,double x_min) |
916 { |
916 { |
917 return exponential(gamma(k,1.0/x_min))+x_min; |
917 return exponential(gamma(k,1.0/x_min))+x_min; |
918 } |
918 } |
919 |
919 |
920 /// Poisson distribution |
920 /// Poisson distribution |
921 |
921 |
922 /// This function generates a Poisson distribution random number with |
922 /// This function generates a Poisson distribution random number with |
923 /// parameter \c lambda. |
923 /// parameter \c lambda. |
924 /// |
924 /// |
925 /// The probability mass function of this distribusion is |
925 /// The probability mass function of this distribusion is |
926 /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f] |
926 /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f] |
927 /// \note The algorithm is taken from the book of Donald E. Knuth titled |
927 /// \note The algorithm is taken from the book of Donald E. Knuth titled |
928 /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the |
928 /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the |
929 /// return value. |
929 /// return value. |
930 |
930 |
931 int poisson(double lambda) |
931 int poisson(double lambda) |
932 { |
932 { |
933 const double l = std::exp(-lambda); |
933 const double l = std::exp(-lambda); |
934 int k=0; |
934 int k=0; |
935 double p = 1.0; |
935 double p = 1.0; |
936 do { |
936 do { |
937 k++; |
937 k++; |
938 p*=real<double>(); |
938 p*=real<double>(); |
939 } while (p>=l); |
939 } while (p>=l); |
940 return k-1; |
940 return k-1; |
941 } |
941 } |
942 |
942 |
943 ///@} |
943 ///@} |
944 |
944 |
945 ///\name Two dimensional distributions |
945 ///\name Two dimensional distributions |
946 /// |
946 /// |
947 |
947 |
948 ///@{ |
948 ///@{ |
949 |
949 |
950 /// Uniform distribution on the full unit circle |
950 /// Uniform distribution on the full unit circle |
951 |
951 |
952 /// Uniform distribution on the full unit circle. |
952 /// Uniform distribution on the full unit circle. |
953 /// |
953 /// |
954 dim2::Point<double> disc() |
954 dim2::Point<double> disc() |
955 { |
955 { |
956 double V1,V2; |
956 double V1,V2; |
957 do { |
957 do { |
958 V1=2*real<double>()-1; |
958 V1=2*real<double>()-1; |
959 V2=2*real<double>()-1; |
959 V2=2*real<double>()-1; |
960 |
960 |
961 } while(V1*V1+V2*V2>=1); |
961 } while(V1*V1+V2*V2>=1); |
962 return dim2::Point<double>(V1,V2); |
962 return dim2::Point<double>(V1,V2); |
963 } |
963 } |
964 /// A kind of two dimensional Gauss distribution |
964 /// A kind of two dimensional Gauss distribution |
965 |
965 |
971 /// the Box-Muller method. |
971 /// the Box-Muller method. |
972 dim2::Point<double> gauss2() |
972 dim2::Point<double> gauss2() |
973 { |
973 { |
974 double V1,V2,S; |
974 double V1,V2,S; |
975 do { |
975 do { |
976 V1=2*real<double>()-1; |
976 V1=2*real<double>()-1; |
977 V2=2*real<double>()-1; |
977 V2=2*real<double>()-1; |
978 S=V1*V1+V2*V2; |
978 S=V1*V1+V2*V2; |
979 } while(S>=1); |
979 } while(S>=1); |
980 double W=std::sqrt(-2*std::log(S)/S); |
980 double W=std::sqrt(-2*std::log(S)/S); |
981 return dim2::Point<double>(W*V1,W*V2); |
981 return dim2::Point<double>(W*V1,W*V2); |
982 } |
982 } |
983 /// A kind of two dimensional exponential distribution |
983 /// A kind of two dimensional exponential distribution |
984 |
984 |
985 /// This function provides a turning symmetric two-dimensional distribution. |
985 /// This function provides a turning symmetric two-dimensional distribution. |
986 /// The x-coordinate is of conditionally exponential distribution |
986 /// The x-coordinate is of conditionally exponential distribution |
987 /// with the condition that x is positive and y=0. If x is negative and |
987 /// with the condition that x is positive and y=0. If x is negative and |
988 /// y=0 then, -x is of exponential distribution. The same is true for the |
988 /// y=0 then, -x is of exponential distribution. The same is true for the |
989 /// y-coordinate. |
989 /// y-coordinate. |
990 dim2::Point<double> exponential2() |
990 dim2::Point<double> exponential2() |
991 { |
991 { |
992 double V1,V2,S; |
992 double V1,V2,S; |
993 do { |
993 do { |
994 V1=2*real<double>()-1; |
994 V1=2*real<double>()-1; |
995 V2=2*real<double>()-1; |
995 V2=2*real<double>()-1; |
996 S=V1*V1+V2*V2; |
996 S=V1*V1+V2*V2; |
997 } while(S>=1); |
997 } while(S>=1); |
998 double W=-std::log(S)/S; |
998 double W=-std::log(S)/S; |
999 return dim2::Point<double>(W*V1,W*V2); |
999 return dim2::Point<double>(W*V1,W*V2); |
1000 } |
1000 } |
1001 |
1001 |
1002 ///@} |
1002 ///@} |
1003 }; |
1003 }; |
1004 |
1004 |
1005 |
1005 |
1006 extern Random rnd; |
1006 extern Random rnd; |
1007 |
1007 |