]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | ||
25 | ibeta_fraction1_t(T a_, T b_, T x_) : a(a_), b(b_), x(x_), k(1) {} | |
26 | ||
27 | 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> | |
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> | |
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> | |
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> | |
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 trunctate 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; | |
158 | float force_truncate(const float* f) | |
159 | { | |
160 | external_f = *f; | |
161 | return external_f; | |
162 | } | |
163 | ||
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 | { | |
178 | 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 | { | |
193 | 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 | { | |
216 | 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 | { | |
230 | 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 | ||
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 |