blitz  Version 0.9
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
F.h
Go to the documentation of this file.
1 /*
2  * F distribution
3  *
4  * This code has been adapted from RANDLIB.C 1.3, by
5  * Barry W. Brown, James Lovato, Kathy Russell, and John Venier.
6  * Code was originally by Ahrens and Dieter (see above).
7  *
8  * Adapter's notes:
9  * BZ_NEEDS_WORK: how to handle seeding for the two gamma RNGs if
10  * independentState is used?
11  * BZ_NEEDS_WORK: code for handling possible overflow when xden
12  * is tiny seems a bit flaky.
13  */
14 
15 #ifndef BZ_RANDOM_F
16 #define BZ_RANDOM_F
17 
18 #ifndef BZ_RANDOM_GAMMA
19  #include <random/gamma.h>
20 #endif
21 
22 BZ_NAMESPACE(ranlib)
23 
24 template<typename T = double, typename IRNG = defaultIRNG,
25  typename stateTag = defaultState>
26 class F {
27 public:
28  typedef T T_numtype;
29 
30  F(T numeratorDF, T denominatorDF)
31  {
32  setDF(numeratorDF, denominatorDF);
33  mindenom = 0.085 * tiny(T());
34  }
35 
36  void setDF(T _dfn, T _dfd)
37  {
38  BZPRECONDITION(_dfn > 0.0);
39  BZPRECONDITION(_dfd > 0.0);
40  dfn = _dfn;
41  dfd = _dfd;
42 
43  ngamma.setMean(dfn/2.0);
44  dgamma.setMean(dfd/2.0);
45  }
46 
47  T random()
48  {
49  T xnum = 2.0 * ngamma.random() / dfn;
50  T xden = 2.0 * ngamma.random() / dfd;
51 
52  // Rare event: Will an overflow probably occur?
53  if (xden <= mindenom)
54  {
55  // Yes, just return huge(T())
56  return huge(T());
57  }
58 
59  return xnum / xden;
60  }
61 
62  void seed(IRNG_int s)
63  {
64  // This is such a bad idea if independentState is used. Ugh.
65  // If sharedState is used, it is merely inefficient (the
66  // same RNG is seeded twice).
67 
68  ngamma.seed(s);
69  dgamma.seed(s);
70  }
71 
72 protected:
74  T dfn, dfd, mindenom;
75 };
76 
78 
79 #endif // BZ_RANDOM_F