blitz Version 0.9
random/F.h
Go to the documentation of this file.
00001 /*
00002  * F distribution
00003  *
00004  * This code has been adapted from RANDLIB.C 1.3, by
00005  * Barry W. Brown, James Lovato, Kathy Russell, and John Venier.
00006  * Code was originally by Ahrens and Dieter (see above).
00007  *
00008  * Adapter's notes:
00009  * BZ_NEEDS_WORK: how to handle seeding for the two gamma RNGs if
00010  * independentState is used?
00011  * BZ_NEEDS_WORK: code for handling possible overflow when xden
00012  * is tiny seems a bit flaky.
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 * tiny(T());
00034     }
00035 
00036     void setDF(T _dfn, T _dfd)
00037     {
00038         BZPRECONDITION(_dfn > 0.0);
00039         BZPRECONDITION(_dfd > 0.0);
00040         dfn = _dfn;
00041         dfd = _dfd;
00042 
00043         ngamma.setMean(dfn/2.0);
00044         dgamma.setMean(dfd/2.0);
00045     }
00046 
00047     T random()
00048     {
00049         T xnum = 2.0 * ngamma.random() / dfn;
00050         T xden = 2.0 * ngamma.random() / dfd;
00051 
00052         // Rare event: Will an overflow probably occur?
00053         if (xden <= mindenom)
00054         {
00055             // Yes, just return huge(T())
00056             return huge(T());
00057         }
00058 
00059         return xnum / xden;
00060     }
00061 
00062     void seed(IRNG_int s)
00063     {
00064         // This is such a bad idea if independentState is used. Ugh.
00065         // If sharedState is used, it is merely inefficient (the
00066         // same RNG is seeded twice).
00067 
00068         ngamma.seed(s);
00069         dgamma.seed(s);
00070     }
00071 
00072 protected:
00073     Gamma<T,IRNG,stateTag> ngamma, dgamma;
00074     T dfn, dfd, mindenom;
00075 };
00076 
00077 BZ_NAMESPACE_END
00078 
00079 #endif // BZ_RANDOM_F
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines