00001 // -*- C++ -*- 00002 /*************************************************************************** 00003 * blitz/rand-normal.h Random Gaussian (Normal) generator 00004 * 00005 * $Id: rand-normal.h,v 1.5 2005/10/19 21:17:00 julianc Exp $ 00006 * 00007 * Copyright (C) 1997-2001 Todd Veldhuizen <tveldhui@oonumerics.org> 00008 * 00009 * This program is free software; you can redistribute it and/or 00010 * modify it under the terms of the GNU General Public License 00011 * as published by the Free Software Foundation; either version 2 00012 * of the License, or (at your option) any later version. 00013 * 00014 * This program is distributed in the hope that it will be useful, 00015 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00017 * GNU General Public License for more details. 00018 * 00019 * Suggestions: blitz-dev@oonumerics.org 00020 * Bugs: blitz-bugs@oonumerics.org 00021 * 00022 * For more information, please see the Blitz++ Home Page: 00023 * http://oonumerics.org/blitz/ 00024 * 00025 *************************************************************************** 00026 * 00027 * This generator transforms a (0,1] uniform distribution into 00028 * a Normal distribution. Let u,v be (0,1] random variables. Then 00029 * 00030 * x = sqrt(-2 ln v) cos(pi (2u-1)) 00031 * 00032 * is N(0,1) distributed. 00033 * 00034 * Reference: Athanasios Papoulis, "Probability, random variables, 00035 * and stochastic processes," McGraw-Hill : Toronto, 1991. 00036 * 00037 ***************************************************************************/ 00038 00039 #ifndef BZ_RAND_NORMAL_H 00040 #define BZ_RAND_NORMAL_H 00041 00042 #include <blitz/random.h> 00043 #include <blitz/rand-uniform.h> 00044 00045 #if defined(BZ_HAVE_STD) 00046 #include <cmath> 00047 #else 00048 #include <math.h> 00049 #endif 00050 00051 #ifndef M_PI 00052 #define M_PI 3.14159265358979323846 /* pi */ 00053 #endif 00054 00055 BZ_NAMESPACE(blitz) 00056 00057 template<typename P_uniform BZ_TEMPLATE_DEFAULT(Uniform)> 00058 class Normal { 00059 00060 public: 00061 typedef double T_numtype; 00062 00063 Normal(double mean = 0.0, double variance = 1.0, double = 0.0) 00064 : mean_(mean), sigma_(BZ_MATHFN_SCOPE(sqrt)(variance)) 00065 { 00066 } 00067 00068 void randomize() 00069 { 00070 uniform_.randomize(); 00071 } 00072 00073 double random() 00074 { 00075 double u, v; 00076 00077 do { 00078 u = uniform_.random(); 00079 v = uniform_.random(); 00080 } while (v == 0); 00081 00082 return mean_ + sigma_ * BZ_MATHFN_SCOPE(sqrt)(-2*BZ_MATHFN_SCOPE(log)(v)) * BZ_MATHFN_SCOPE(cos)(M_PI * (2*u - 1)); 00083 } 00084 00085 private: 00086 double mean_, sigma_; 00087 P_uniform uniform_; 00088 }; 00089 00090 BZ_NAMESPACE_END 00091 00092 #endif // BZ_RAND_NORMAL_H 00093