00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015 #ifndef BZ_RANDOM_F
00016 #define BZ_RANDOM_F
00017
00018 #ifndef BZ_RANDOM_GAMMA
00019 #include <random/gamma.h>
00020 #endif
00021
00022 BZ_NAMESPACE(ranlib)
00023
00024 template<typename T = double, typename IRNG = defaultIRNG,
00025 typename stateTag = defaultState>
00026 class F {
00027 public:
00028 typedef T T_numtype;
00029
00030 F(T numeratorDF, T denominatorDF)
00031 {
00032 setDF(numeratorDF, denominatorDF);
00033 mindenom = 0.085 * blitz::tiny(T());
00034 }
00035
00036 F(T numeratorDF, T denominatorDF, unsigned int i) :
00037 ngamma(i), dgamma(i)
00038 {
00039 setDF(numeratorDF, denominatorDF);
00040 mindenom = 0.085 * blitz::tiny(T());
00041 }
00042
00043 void setDF(T _dfn, T _dfd)
00044 {
00045 BZPRECONDITION(_dfn > 0.0);
00046 BZPRECONDITION(_dfd > 0.0);
00047 dfn = _dfn;
00048 dfd = _dfd;
00049
00050 ngamma.setMean(dfn/2.0);
00051 dgamma.setMean(dfd/2.0);
00052 }
00053
00054 T random()
00055 {
00056 T xnum = 2.0 * ngamma.random() / dfn;
00057 T xden = 2.0 * ngamma.random() / dfd;
00058
00059
00060 if (xden <= mindenom)
00061 {
00062
00063 return blitz::huge(T());
00064 }
00065
00066 return xnum / xden;
00067 }
00068
00069 void seed(IRNG_int s, IRNG_int r)
00070 {
00071
00072
00073
00074
00075
00076
00077
00078 ngamma.seed(s);
00079 dgamma.seed(r);
00080 }
00081
00082 protected:
00083 Gamma<T,IRNG,stateTag> ngamma, dgamma;
00084 T dfn, dfd, mindenom;
00085 };
00086
00087 BZ_NAMESPACE_END
00088
00089 #endif // BZ_RANDOM_F