1 // (C) Copyright John Maddock 2006.
2 // Use, modification and distribution are subject to the
3 // Boost Software License, Version 1.0. (See accompanying file
4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5
6 #include <boost/math/special_functions/gamma.hpp>
7 #include <boost/math/special_functions/beta.hpp>
8 #include <boost/math/constants/constants.hpp>
9 #include <boost/lexical_cast.hpp>
10 #include <fstream>
11 #include <map>
12 #include <boost/math/tools/test_data.hpp>
13 #include <boost/random.hpp>
14 #include "mp_t.hpp"
15
16 using namespace boost::math::tools;
17 using namespace boost::math;
18 using namespace std;
19
20 template <class T>
21 struct ibeta_fraction1_t
22 {
23 typedef std::pair<T, T> result_type;
24
ibeta_fraction1_tibeta_fraction1_t25 ibeta_fraction1_t(T a_, T b_, T x_) : a(a_), b(b_), x(x_), k(1) {}
26
operator ()ibeta_fraction1_t27 result_type operator()()
28 {
29 T aN;
30 if(k & 1)
31 {
32 int m = (k - 1) / 2;
33 aN = -(a + m) * (a + b + m) * x;
34 aN /= a + 2*m;
35 aN /= a + 2*m + 1;
36 }
37 else
38 {
39 int m = k / 2;
40 aN = m * (b - m) *x;
41 aN /= a + 2*m - 1;
42 aN /= a + 2*m;
43 }
44 ++k;
45 return std::make_pair(aN, T(1));
46 }
47
48 private:
49 T a, b, x;
50 int k;
51 };
52
53 //
54 // This function caches previous calls to beta
55 // just so we can speed things up a bit:
56 //
57 template <class T>
get_beta(T a,T b)58 T get_beta(T a, T b)
59 {
60 static std::map<std::pair<T, T>, T> m;
61
62 if(a < b)
63 std::swap(a, b);
64
65 std::pair<T, T> p(a, b);
66 typename std::map<std::pair<T, T>, T>::const_iterator i = m.find(p);
67 if(i != m.end())
68 return i->second;
69
70 T r = beta(a, b);
71 p.first = a;
72 p.second = b;
73 m[p] = r;
74
75 return r;
76 }
77
78 //
79 // compute the continued fraction:
80 //
81 template <class T>
get_ibeta_fraction1(T a,T b,T x)82 T get_ibeta_fraction1(T a, T b, T x)
83 {
84 ibeta_fraction1_t<T> f(a, b, x);
85 T fract = boost::math::tools::continued_fraction_a(f, boost::math::policies::digits<T, boost::math::policies::policy<> >());
86 T denom = (a * (fract + 1));
87 T num = pow(x, a) * pow(1 - x, b);
88 if(num == 0)
89 return 0;
90 else if(denom == 0)
91 return -1;
92 return num / denom;
93 }
94 //
95 // calculate the incomplete beta from the fraction:
96 //
97 template <class T>
ibeta_fraction1(T a,T b,T x)98 std::pair<T,T> ibeta_fraction1(T a, T b, T x)
99 {
100 T bet = get_beta(a, b);
101 if(x > ((a+1)/(a+b+2)))
102 {
103 T fract = get_ibeta_fraction1(b, a, 1-x);
104 if(fract/bet > 0.75)
105 {
106 fract = get_ibeta_fraction1(a, b, x);
107 return std::make_pair(fract, bet - fract);
108 }
109 return std::make_pair(bet - fract, fract);
110 }
111 T fract = get_ibeta_fraction1(a, b, x);
112 if(fract/bet > 0.75)
113 {
114 fract = get_ibeta_fraction1(b, a, 1-x);
115 return std::make_pair(bet - fract, fract);
116 }
117 return std::make_pair(fract, bet - fract);
118
119 }
120 //
121 // calculate the regularised incomplete beta from the fraction:
122 //
123 template <class T>
ibeta_fraction1_regular(T a,T b,T x)124 std::pair<T,T> ibeta_fraction1_regular(T a, T b, T x)
125 {
126 T bet = get_beta(a, b);
127 if(x > ((a+1)/(a+b+2)))
128 {
129 T fract = get_ibeta_fraction1(b, a, 1-x);
130 if(fract == 0)
131 bet = 1; // normalise so we don't get 0/0
132 else if(bet == 0)
133 return std::make_pair(T(-1), T(-1)); // Yikes!!
134 if(fract / bet > 0.75)
135 {
136 fract = get_ibeta_fraction1(a, b, x);
137 return std::make_pair(fract / bet, 1 - (fract / bet));
138 }
139 return std::make_pair(1 - (fract / bet), fract / bet);
140 }
141 T fract = get_ibeta_fraction1(a, b, x);
142 if(fract / bet > 0.75)
143 {
144 fract = get_ibeta_fraction1(b, a, 1-x);
145 return std::make_pair(1 - (fract / bet), fract / bet);
146 }
147 return std::make_pair(fract / bet, 1 - (fract / bet));
148 }
149
150 //
151 // we absolutely must truncate the input values to float
152 // precision: we have to be certain that the input values
153 // can be represented exactly in whatever width floating
154 // point type we are testing, otherwise the output will
155 // necessarily be off.
156 //
157 float external_f;
force_truncate(const float * f)158 float force_truncate(const float* f)
159 {
160 external_f = *f;
161 return external_f;
162 }
163
truncate_to_float(mp_t r)164 float truncate_to_float(mp_t r)
165 {
166 float f = boost::math::tools::real_cast<float>(r);
167 return force_truncate(&f);
168 }
169
170 boost::mt19937 rnd;
171 boost::uniform_real<float> ur_a(1.0F, 5.0F);
172 boost::variate_generator<boost::mt19937, boost::uniform_real<float> > gen(rnd, ur_a);
173 boost::uniform_real<float> ur_a2(0.0F, 100.0F);
174 boost::variate_generator<boost::mt19937, boost::uniform_real<float> > gen2(rnd, ur_a2);
175
176 struct beta_data_generator
177 {
operator ()beta_data_generator178 boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t, mp_t, mp_t> operator()(mp_t ap, mp_t bp, mp_t x_)
179 {
180 float a = truncate_to_float(real_cast<float>(gen() * pow(mp_t(10), ap)));
181 float b = truncate_to_float(real_cast<float>(gen() * pow(mp_t(10), bp)));
182 float x = truncate_to_float(real_cast<float>(x_));
183 std::cout << a << " " << b << " " << x << std::endl;
184 std::pair<mp_t, mp_t> ib_full = ibeta_fraction1(mp_t(a), mp_t(b), mp_t(x));
185 std::pair<mp_t, mp_t> ib_reg = ibeta_fraction1_regular(mp_t(a), mp_t(b), mp_t(x));
186 return boost::math::make_tuple(a, b, x, ib_full.first, ib_full.second, ib_reg.first, ib_reg.second);
187 }
188 };
189
190 // medium sized values:
191 struct beta_data_generator_medium
192 {
operator ()beta_data_generator_medium193 boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t, mp_t, mp_t> operator()(mp_t x_)
194 {
195 mp_t a = gen2();
196 mp_t b = gen2();
197 mp_t x = x_;
198 a = ConvPrec(a, 22);
199 b = ConvPrec(b, 22);
200 x = ConvPrec(x, 22);
201 std::cout << a << " " << b << " " << x << std::endl;
202 //mp_t exp_beta = boost::math::beta(a, b, x);
203 std::pair<mp_t, mp_t> ib_full = ibeta_fraction1(mp_t(a), mp_t(b), mp_t(x));
204 /*exp_beta = boost::math::tools::relative_error(ib_full.first, exp_beta);
205 if(exp_beta > 1e-40)
206 {
207 std::cout << exp_beta << std::endl;
208 }*/
209 std::pair<mp_t, mp_t> ib_reg = ibeta_fraction1_regular(mp_t(a), mp_t(b), mp_t(x));
210 return boost::math::make_tuple(a, b, x, ib_full.first, ib_full.second, ib_reg.first, ib_reg.second);
211 }
212 };
213
214 struct beta_data_generator_small
215 {
operator ()beta_data_generator_small216 boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t, mp_t, mp_t> operator()(mp_t x_)
217 {
218 float a = truncate_to_float(gen2()/10);
219 float b = truncate_to_float(gen2()/10);
220 float x = truncate_to_float(real_cast<float>(x_));
221 std::cout << a << " " << b << " " << x << std::endl;
222 std::pair<mp_t, mp_t> ib_full = ibeta_fraction1(mp_t(a), mp_t(b), mp_t(x));
223 std::pair<mp_t, mp_t> ib_reg = ibeta_fraction1_regular(mp_t(a), mp_t(b), mp_t(x));
224 return boost::math::make_tuple(a, b, x, ib_full.first, ib_full.second, ib_reg.first, ib_reg.second);
225 }
226 };
227
228 struct beta_data_generator_int
229 {
operator ()beta_data_generator_int230 boost::math::tuple<mp_t, mp_t, mp_t, mp_t, mp_t, mp_t, mp_t> operator()(mp_t a, mp_t b, mp_t x_)
231 {
232 float x = truncate_to_float(real_cast<float>(x_));
233 std::cout << a << " " << b << " " << x << std::endl;
234 std::pair<mp_t, mp_t> ib_full = ibeta_fraction1(a, b, mp_t(x));
235 std::pair<mp_t, mp_t> ib_reg = ibeta_fraction1_regular(a, b, mp_t(x));
236 return boost::math::make_tuple(a, b, x, ib_full.first, ib_full.second, ib_reg.first, ib_reg.second);
237 }
238 };
239
main(int,char * [])240 int main(int, char* [])
241 {
242 parameter_info<mp_t> arg1, arg2, arg3, arg4, arg5;
243 test_data<mp_t> data;
244
245 std::cout << "Welcome.\n"
246 "This program will generate spot tests for the incomplete beta functions:\n"
247 " beta(a, b, x) and ibeta(a, b, x)\n\n"
248 "This is not an interactive program be prepared for a long wait!!!\n\n";
249
250 arg1 = make_periodic_param(mp_t(-5), mp_t(6), 11);
251 arg2 = make_periodic_param(mp_t(-5), mp_t(6), 11);
252 arg3 = make_random_param(mp_t(0.0001), mp_t(1), 10);
253 arg4 = make_random_param(mp_t(0.0001), mp_t(1), 100 /*500*/);
254 arg5 = make_periodic_param(mp_t(1), mp_t(41), 10);
255
256 arg1.type |= dummy_param;
257 arg2.type |= dummy_param;
258 arg3.type |= dummy_param;
259 arg4.type |= dummy_param;
260 arg5.type |= dummy_param;
261
262 // comment out all but one of the following when running
263 // or this program will take forever to complete!
264 //data.insert(beta_data_generator(), arg1, arg2, arg3);
265 //data.insert(beta_data_generator_medium(), arg4);
266 //data.insert(beta_data_generator_small(), arg4);
267 data.insert(beta_data_generator_int(), arg5, arg5, arg3);
268
269 test_data<mp_t>::const_iterator i, j;
270 i = data.begin();
271 j = data.end();
272 while(i != j)
273 {
274 mp_t v1 = beta((*i)[0], (*i)[1], (*i)[2]);
275 mp_t v2 = relative_error(v1, (*i)[3]);
276 std::string s = boost::lexical_cast<std::string>((*i)[3]);
277 mp_t v3 = boost::lexical_cast<mp_t>(s);
278 mp_t v4 = relative_error(v3, (*i)[3]);
279 if(v2 > 1e-40)
280 {
281 std::cout << v2 << std::endl;
282 }
283 if(v4 > 1e-60)
284 {
285 std::cout << v4 << std::endl;
286 }
287 ++ i;
288 }
289
290 std::cout << "Enter name of test data file [default=ibeta_data.ipp]";
291 std::string line;
292 std::getline(std::cin, line);
293 boost::algorithm::trim(line);
294 if(line == "")
295 line = "ibeta_data.ipp";
296 std::ofstream ofs(line.c_str());
297 ofs << std::scientific << std::setprecision(40);
298 write_code(ofs, data, "ibeta_data");
299
300 return 0;
301 }
302
303
304