blitz Version 0.9
|
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