]>
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 | #ifndef BOOST_MATH_SPECIAL_BETA_HPP | |
7 | #define BOOST_MATH_SPECIAL_BETA_HPP | |
8 | ||
9 | #ifdef _MSC_VER | |
10 | #pragma once | |
11 | #endif | |
12 | ||
13 | #include <boost/math/special_functions/math_fwd.hpp> | |
14 | #include <boost/math/tools/config.hpp> | |
15 | #include <boost/math/special_functions/gamma.hpp> | |
16 | #include <boost/math/special_functions/binomial.hpp> | |
17 | #include <boost/math/special_functions/factorials.hpp> | |
18 | #include <boost/math/special_functions/erf.hpp> | |
19 | #include <boost/math/special_functions/log1p.hpp> | |
20 | #include <boost/math/special_functions/expm1.hpp> | |
21 | #include <boost/math/special_functions/trunc.hpp> | |
22 | #include <boost/math/tools/roots.hpp> | |
23 | #include <boost/static_assert.hpp> | |
24 | #include <boost/config/no_tr1/cmath.hpp> | |
25 | ||
26 | namespace boost{ namespace math{ | |
27 | ||
28 | namespace detail{ | |
29 | ||
30 | // | |
31 | // Implementation of Beta(a,b) using the Lanczos approximation: | |
32 | // | |
33 | template <class T, class Lanczos, class Policy> | |
34 | T beta_imp(T a, T b, const Lanczos&, const Policy& pol) | |
35 | { | |
36 | BOOST_MATH_STD_USING // for ADL of std names | |
37 | ||
38 | if(a <= 0) | |
39 | return policies::raise_domain_error<T>("boost::math::beta<%1%>(%1%,%1%)", "The arguments to the beta function must be greater than zero (got a=%1%).", a, pol); | |
40 | if(b <= 0) | |
41 | return policies::raise_domain_error<T>("boost::math::beta<%1%>(%1%,%1%)", "The arguments to the beta function must be greater than zero (got b=%1%).", b, pol); | |
42 | ||
43 | T result; | |
44 | ||
45 | T prefix = 1; | |
46 | T c = a + b; | |
47 | ||
48 | // Special cases: | |
49 | if((c == a) && (b < tools::epsilon<T>())) | |
50 | return 1 / b; | |
51 | else if((c == b) && (a < tools::epsilon<T>())) | |
52 | return 1 / a; | |
53 | if(b == 1) | |
54 | return 1/a; | |
55 | else if(a == 1) | |
56 | return 1/b; | |
57 | else if(c < tools::epsilon<T>()) | |
58 | { | |
59 | result = c / a; | |
60 | result /= b; | |
61 | return result; | |
62 | } | |
63 | ||
64 | /* | |
65 | // | |
66 | // This code appears to be no longer necessary: it was | |
67 | // used to offset errors introduced from the Lanczos | |
68 | // approximation, but the current Lanczos approximations | |
69 | // are sufficiently accurate for all z that we can ditch | |
70 | // this. It remains in the file for future reference... | |
71 | // | |
72 | // If a or b are less than 1, shift to greater than 1: | |
73 | if(a < 1) | |
74 | { | |
75 | prefix *= c / a; | |
76 | c += 1; | |
77 | a += 1; | |
78 | } | |
79 | if(b < 1) | |
80 | { | |
81 | prefix *= c / b; | |
82 | c += 1; | |
83 | b += 1; | |
84 | } | |
85 | */ | |
86 | ||
87 | if(a < b) | |
88 | std::swap(a, b); | |
89 | ||
90 | // Lanczos calculation: | |
91 | T agh = static_cast<T>(a + Lanczos::g() - 0.5f); | |
92 | T bgh = static_cast<T>(b + Lanczos::g() - 0.5f); | |
93 | T cgh = static_cast<T>(c + Lanczos::g() - 0.5f); | |
94 | result = Lanczos::lanczos_sum_expG_scaled(a) * (Lanczos::lanczos_sum_expG_scaled(b) / Lanczos::lanczos_sum_expG_scaled(c)); | |
95 | T ambh = a - 0.5f - b; | |
96 | if((fabs(b * ambh) < (cgh * 100)) && (a > 100)) | |
97 | { | |
98 | // Special case where the base of the power term is close to 1 | |
99 | // compute (1+x)^y instead: | |
100 | result *= exp(ambh * boost::math::log1p(-b / cgh, pol)); | |
101 | } | |
102 | else | |
103 | { | |
104 | result *= pow(agh / cgh, a - T(0.5) - b); | |
105 | } | |
106 | if(cgh > 1e10f) | |
107 | // this avoids possible overflow, but appears to be marginally less accurate: | |
108 | result *= pow((agh / cgh) * (bgh / cgh), b); | |
109 | else | |
110 | result *= pow((agh * bgh) / (cgh * cgh), b); | |
111 | result *= sqrt(boost::math::constants::e<T>() / bgh); | |
112 | ||
113 | // If a and b were originally less than 1 we need to scale the result: | |
114 | result *= prefix; | |
115 | ||
116 | return result; | |
117 | } // template <class T, class Lanczos> beta_imp(T a, T b, const Lanczos&) | |
118 | ||
119 | // | |
120 | // Generic implementation of Beta(a,b) without Lanczos approximation support | |
121 | // (Caution this is slow!!!): | |
122 | // | |
123 | template <class T, class Policy> | |
92f5a8d4 | 124 | T beta_imp(T a, T b, const lanczos::undefined_lanczos& l, const Policy& pol) |
7c673cae FG |
125 | { |
126 | BOOST_MATH_STD_USING | |
127 | ||
128 | if(a <= 0) | |
129 | return policies::raise_domain_error<T>("boost::math::beta<%1%>(%1%,%1%)", "The arguments to the beta function must be greater than zero (got a=%1%).", a, pol); | |
130 | if(b <= 0) | |
131 | return policies::raise_domain_error<T>("boost::math::beta<%1%>(%1%,%1%)", "The arguments to the beta function must be greater than zero (got b=%1%).", b, pol); | |
132 | ||
92f5a8d4 | 133 | const T c = a + b; |
7c673cae | 134 | |
92f5a8d4 TL |
135 | // Special cases: |
136 | if ((c == a) && (b < tools::epsilon<T>())) | |
137 | return 1 / b; | |
138 | else if ((c == b) && (a < tools::epsilon<T>())) | |
139 | return 1 / a; | |
140 | if (b == 1) | |
141 | return 1 / a; | |
142 | else if (a == 1) | |
143 | return 1 / b; | |
144 | else if (c < tools::epsilon<T>()) | |
7c673cae | 145 | { |
92f5a8d4 TL |
146 | T result = c / a; |
147 | result /= b; | |
148 | return result; | |
7c673cae | 149 | } |
7c673cae | 150 | |
92f5a8d4 TL |
151 | // Regular cases start here: |
152 | const T min_sterling = minimum_argument_for_bernoulli_recursion<T>(); | |
7c673cae | 153 | |
92f5a8d4 TL |
154 | long shift_a = 0; |
155 | long shift_b = 0; | |
7c673cae | 156 | |
92f5a8d4 TL |
157 | if(a < min_sterling) |
158 | shift_a = 1 + ltrunc(min_sterling - a); | |
159 | if(b < min_sterling) | |
160 | shift_b = 1 + ltrunc(min_sterling - b); | |
161 | long shift_c = shift_a + shift_b; | |
7c673cae | 162 | |
92f5a8d4 TL |
163 | if ((shift_a == 0) && (shift_b == 0)) |
164 | { | |
165 | return pow(a / c, a) * pow(b / c, b) * scaled_tgamma_no_lanczos(a, pol) * scaled_tgamma_no_lanczos(b, pol) / scaled_tgamma_no_lanczos(c, pol); | |
166 | } | |
167 | else if ((a < 1) && (b < 1)) | |
168 | { | |
169 | return boost::math::tgamma(a, pol) * (boost::math::tgamma(b, pol) / boost::math::tgamma(c)); | |
170 | } | |
171 | else if(a < 1) | |
172 | return boost::math::tgamma(a, pol) * boost::math::tgamma_delta_ratio(b, a, pol); | |
173 | else if(b < 1) | |
174 | return boost::math::tgamma(b, pol) * boost::math::tgamma_delta_ratio(a, b, pol); | |
175 | else | |
176 | { | |
177 | T result = beta_imp(T(a + shift_a), T(b + shift_b), l, pol); | |
178 | // | |
179 | // Recursion: | |
180 | // | |
181 | for (long i = 0; i < shift_c; ++i) | |
182 | { | |
183 | result *= c + i; | |
184 | if (i < shift_a) | |
185 | result /= a + i; | |
186 | if (i < shift_b) | |
187 | result /= b + i; | |
188 | } | |
189 | return result; | |
190 | } | |
7c673cae | 191 | |
7c673cae FG |
192 | } // template <class T>T beta_imp(T a, T b, const lanczos::undefined_lanczos& l) |
193 | ||
194 | ||
195 | // | |
196 | // Compute the leading power terms in the incomplete Beta: | |
197 | // | |
198 | // (x^a)(y^b)/Beta(a,b) when normalised, and | |
199 | // (x^a)(y^b) otherwise. | |
200 | // | |
201 | // Almost all of the error in the incomplete beta comes from this | |
202 | // function: particularly when a and b are large. Computing large | |
203 | // powers are *hard* though, and using logarithms just leads to | |
204 | // horrendous cancellation errors. | |
205 | // | |
206 | template <class T, class Lanczos, class Policy> | |
207 | T ibeta_power_terms(T a, | |
208 | T b, | |
209 | T x, | |
210 | T y, | |
211 | const Lanczos&, | |
212 | bool normalised, | |
213 | const Policy& pol, | |
214 | T prefix = 1, | |
215 | const char* function = "boost::math::ibeta<%1%>(%1%, %1%, %1%)") | |
216 | { | |
217 | BOOST_MATH_STD_USING | |
218 | ||
219 | if(!normalised) | |
220 | { | |
221 | // can we do better here? | |
222 | return pow(x, a) * pow(y, b); | |
223 | } | |
224 | ||
225 | T result; | |
226 | ||
227 | T c = a + b; | |
228 | ||
229 | // combine power terms with Lanczos approximation: | |
230 | T agh = static_cast<T>(a + Lanczos::g() - 0.5f); | |
231 | T bgh = static_cast<T>(b + Lanczos::g() - 0.5f); | |
232 | T cgh = static_cast<T>(c + Lanczos::g() - 0.5f); | |
233 | result = Lanczos::lanczos_sum_expG_scaled(c) / (Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b)); | |
234 | result *= prefix; | |
235 | // combine with the leftover terms from the Lanczos approximation: | |
236 | result *= sqrt(bgh / boost::math::constants::e<T>()); | |
237 | result *= sqrt(agh / cgh); | |
238 | ||
239 | // l1 and l2 are the base of the exponents minus one: | |
240 | T l1 = (x * b - y * agh) / agh; | |
241 | T l2 = (y * a - x * bgh) / bgh; | |
242 | if(((std::min)(fabs(l1), fabs(l2)) < 0.2)) | |
243 | { | |
244 | // when the base of the exponent is very near 1 we get really | |
245 | // gross errors unless extra care is taken: | |
246 | if((l1 * l2 > 0) || ((std::min)(a, b) < 1)) | |
247 | { | |
248 | // | |
92f5a8d4 | 249 | // This first branch handles the simple cases where either: |
7c673cae | 250 | // |
92f5a8d4 TL |
251 | // * The two power terms both go in the same direction |
252 | // (towards zero or towards infinity). In this case if either | |
253 | // term overflows or underflows, then the product of the two must | |
254 | // do so also. | |
255 | // *Alternatively if one exponent is less than one, then we | |
256 | // can't productively use it to eliminate overflow or underflow | |
257 | // from the other term. Problems with spurious overflow/underflow | |
258 | // can't be ruled out in this case, but it is *very* unlikely | |
7c673cae FG |
259 | // since one of the power terms will evaluate to a number close to 1. |
260 | // | |
261 | if(fabs(l1) < 0.1) | |
262 | { | |
263 | result *= exp(a * boost::math::log1p(l1, pol)); | |
264 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
265 | } | |
266 | else | |
267 | { | |
268 | result *= pow((x * cgh) / agh, a); | |
269 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
270 | } | |
271 | if(fabs(l2) < 0.1) | |
272 | { | |
273 | result *= exp(b * boost::math::log1p(l2, pol)); | |
274 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
275 | } | |
276 | else | |
277 | { | |
278 | result *= pow((y * cgh) / bgh, b); | |
279 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
280 | } | |
281 | } | |
282 | else if((std::max)(fabs(l1), fabs(l2)) < 0.5) | |
283 | { | |
284 | // | |
92f5a8d4 TL |
285 | // Both exponents are near one and both the exponents are |
286 | // greater than one and further these two | |
287 | // power terms tend in opposite directions (one towards zero, | |
288 | // the other towards infinity), so we have to combine the terms | |
7c673cae FG |
289 | // to avoid any risk of overflow or underflow. |
290 | // | |
291 | // We do this by moving one power term inside the other, we have: | |
292 | // | |
293 | // (1 + l1)^a * (1 + l2)^b | |
294 | // = ((1 + l1)*(1 + l2)^(b/a))^a | |
295 | // = (1 + l1 + l3 + l1*l3)^a ; l3 = (1 + l2)^(b/a) - 1 | |
296 | // = exp((b/a) * log(1 + l2)) - 1 | |
297 | // | |
298 | // The tricky bit is deciding which term to move inside :-) | |
299 | // By preference we move the larger term inside, so that the | |
300 | // size of the largest exponent is reduced. However, that can | |
301 | // only be done as long as l3 (see above) is also small. | |
302 | // | |
303 | bool small_a = a < b; | |
304 | T ratio = b / a; | |
305 | if((small_a && (ratio * l2 < 0.1)) || (!small_a && (l1 / ratio > 0.1))) | |
306 | { | |
307 | T l3 = boost::math::expm1(ratio * boost::math::log1p(l2, pol), pol); | |
308 | l3 = l1 + l3 + l3 * l1; | |
309 | l3 = a * boost::math::log1p(l3, pol); | |
310 | result *= exp(l3); | |
311 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
312 | } | |
313 | else | |
314 | { | |
315 | T l3 = boost::math::expm1(boost::math::log1p(l1, pol) / ratio, pol); | |
316 | l3 = l2 + l3 + l3 * l2; | |
317 | l3 = b * boost::math::log1p(l3, pol); | |
318 | result *= exp(l3); | |
319 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
320 | } | |
321 | } | |
322 | else if(fabs(l1) < fabs(l2)) | |
323 | { | |
324 | // First base near 1 only: | |
325 | T l = a * boost::math::log1p(l1, pol) | |
326 | + b * log((y * cgh) / bgh); | |
327 | if((l <= tools::log_min_value<T>()) || (l >= tools::log_max_value<T>())) | |
328 | { | |
329 | l += log(result); | |
330 | if(l >= tools::log_max_value<T>()) | |
331 | return policies::raise_overflow_error<T>(function, 0, pol); | |
332 | result = exp(l); | |
333 | } | |
334 | else | |
335 | result *= exp(l); | |
336 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
337 | } | |
338 | else | |
339 | { | |
340 | // Second base near 1 only: | |
341 | T l = b * boost::math::log1p(l2, pol) | |
342 | + a * log((x * cgh) / agh); | |
343 | if((l <= tools::log_min_value<T>()) || (l >= tools::log_max_value<T>())) | |
344 | { | |
345 | l += log(result); | |
346 | if(l >= tools::log_max_value<T>()) | |
347 | return policies::raise_overflow_error<T>(function, 0, pol); | |
348 | result = exp(l); | |
349 | } | |
350 | else | |
351 | result *= exp(l); | |
352 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
353 | } | |
354 | } | |
355 | else | |
356 | { | |
357 | // general case: | |
358 | T b1 = (x * cgh) / agh; | |
359 | T b2 = (y * cgh) / bgh; | |
360 | l1 = a * log(b1); | |
361 | l2 = b * log(b2); | |
362 | BOOST_MATH_INSTRUMENT_VARIABLE(b1); | |
363 | BOOST_MATH_INSTRUMENT_VARIABLE(b2); | |
364 | BOOST_MATH_INSTRUMENT_VARIABLE(l1); | |
365 | BOOST_MATH_INSTRUMENT_VARIABLE(l2); | |
366 | if((l1 >= tools::log_max_value<T>()) | |
367 | || (l1 <= tools::log_min_value<T>()) | |
368 | || (l2 >= tools::log_max_value<T>()) | |
369 | || (l2 <= tools::log_min_value<T>()) | |
370 | ) | |
371 | { | |
372 | // Oops, under/overflow, sidestep if we can: | |
373 | if(a < b) | |
374 | { | |
375 | T p1 = pow(b2, b / a); | |
376 | T l3 = a * (log(b1) + log(p1)); | |
377 | if((l3 < tools::log_max_value<T>()) | |
378 | && (l3 > tools::log_min_value<T>())) | |
379 | { | |
380 | result *= pow(p1 * b1, a); | |
381 | } | |
382 | else | |
383 | { | |
384 | l2 += l1 + log(result); | |
385 | if(l2 >= tools::log_max_value<T>()) | |
386 | return policies::raise_overflow_error<T>(function, 0, pol); | |
387 | result = exp(l2); | |
388 | } | |
389 | } | |
390 | else | |
391 | { | |
392 | T p1 = pow(b1, a / b); | |
393 | T l3 = (log(p1) + log(b2)) * b; | |
394 | if((l3 < tools::log_max_value<T>()) | |
395 | && (l3 > tools::log_min_value<T>())) | |
396 | { | |
397 | result *= pow(p1 * b2, b); | |
398 | } | |
399 | else | |
400 | { | |
401 | l2 += l1 + log(result); | |
402 | if(l2 >= tools::log_max_value<T>()) | |
403 | return policies::raise_overflow_error<T>(function, 0, pol); | |
404 | result = exp(l2); | |
405 | } | |
406 | } | |
407 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
408 | } | |
409 | else | |
410 | { | |
411 | // finally the normal case: | |
412 | result *= pow(b1, a) * pow(b2, b); | |
413 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
414 | } | |
415 | } | |
416 | ||
417 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
418 | ||
419 | return result; | |
420 | } | |
421 | // | |
422 | // Compute the leading power terms in the incomplete Beta: | |
423 | // | |
424 | // (x^a)(y^b)/Beta(a,b) when normalised, and | |
425 | // (x^a)(y^b) otherwise. | |
426 | // | |
427 | // Almost all of the error in the incomplete beta comes from this | |
428 | // function: particularly when a and b are large. Computing large | |
429 | // powers are *hard* though, and using logarithms just leads to | |
430 | // horrendous cancellation errors. | |
431 | // | |
432 | // This version is generic, slow, and does not use the Lanczos approximation. | |
433 | // | |
434 | template <class T, class Policy> | |
435 | T ibeta_power_terms(T a, | |
436 | T b, | |
437 | T x, | |
438 | T y, | |
92f5a8d4 | 439 | const boost::math::lanczos::undefined_lanczos& l, |
7c673cae | 440 | bool normalised, |
92f5a8d4 | 441 | const Policy& pol, |
7c673cae FG |
442 | T prefix = 1, |
443 | const char* = "boost::math::ibeta<%1%>(%1%, %1%, %1%)") | |
444 | { | |
445 | BOOST_MATH_STD_USING | |
446 | ||
447 | if(!normalised) | |
448 | { | |
92f5a8d4 | 449 | return prefix * pow(x, a) * pow(y, b); |
7c673cae FG |
450 | } |
451 | ||
7c673cae FG |
452 | T c = a + b; |
453 | ||
92f5a8d4 TL |
454 | const T min_sterling = minimum_argument_for_bernoulli_recursion<T>(); |
455 | ||
456 | long shift_a = 0; | |
457 | long shift_b = 0; | |
458 | ||
459 | if (a < min_sterling) | |
460 | shift_a = 1 + ltrunc(min_sterling - a); | |
461 | if (b < min_sterling) | |
462 | shift_b = 1 + ltrunc(min_sterling - b); | |
463 | ||
464 | if ((shift_a == 0) && (shift_b == 0)) | |
7c673cae | 465 | { |
92f5a8d4 TL |
466 | T power1, power2; |
467 | if (a < b) | |
468 | { | |
469 | power1 = pow((x * y * c * c) / (a * b), a); | |
470 | power2 = pow((y * c) / b, b - a); | |
471 | } | |
472 | else | |
473 | { | |
474 | power1 = pow((x * y * c * c) / (a * b), b); | |
475 | power2 = pow((x * c) / a, a - b); | |
476 | } | |
477 | if (!(boost::math::isnormal)(power1) || !(boost::math::isnormal)(power2)) | |
478 | { | |
479 | // We have to use logs :( | |
480 | return prefix * exp(a * log(x * c / a) + b * log(y * c / b)) * scaled_tgamma_no_lanczos(c, pol) / (scaled_tgamma_no_lanczos(a, pol) * scaled_tgamma_no_lanczos(b, pol)); | |
481 | } | |
482 | return prefix * power1 * power2 * scaled_tgamma_no_lanczos(c, pol) / (scaled_tgamma_no_lanczos(a, pol) * scaled_tgamma_no_lanczos(b, pol)); | |
7c673cae | 483 | } |
92f5a8d4 TL |
484 | |
485 | T power1 = pow(x, a); | |
486 | T power2 = pow(y, b); | |
487 | T bet = beta_imp(a, b, l, pol); | |
488 | ||
489 | if(!(boost::math::isnormal)(power1) || !(boost::math::isnormal)(power2) || !(boost::math::isnormal)(bet)) | |
7c673cae | 490 | { |
92f5a8d4 TL |
491 | int shift_c = shift_a + shift_b; |
492 | T result = ibeta_power_terms(T(a + shift_a), T(b + shift_b), x, y, l, normalised, pol, prefix); | |
493 | if ((boost::math::isnormal)(result)) | |
494 | { | |
495 | for (int i = 0; i < shift_c; ++i) | |
496 | { | |
497 | result /= c + i; | |
498 | if (i < shift_a) | |
499 | { | |
500 | result *= a + i; | |
501 | result /= x; | |
502 | } | |
503 | if (i < shift_b) | |
504 | { | |
505 | result *= b + i; | |
506 | result /= y; | |
507 | } | |
508 | } | |
509 | return prefix * result; | |
510 | } | |
7c673cae | 511 | else |
92f5a8d4 TL |
512 | { |
513 | T log_result = log(x) * a + log(y) * b + log(prefix); | |
514 | if ((boost::math::isnormal)(bet)) | |
515 | log_result -= log(bet); | |
516 | else | |
517 | log_result += boost::math::lgamma(c, pol) - boost::math::lgamma(a) - boost::math::lgamma(c, pol); | |
518 | return exp(log_result); | |
519 | } | |
7c673cae | 520 | } |
92f5a8d4 | 521 | return prefix * power1 * (power2 / bet); |
7c673cae FG |
522 | } |
523 | // | |
524 | // Series approximation to the incomplete beta: | |
525 | // | |
526 | template <class T> | |
527 | struct ibeta_series_t | |
528 | { | |
529 | typedef T result_type; | |
530 | ibeta_series_t(T a_, T b_, T x_, T mult) : result(mult), x(x_), apn(a_), poch(1-b_), n(1) {} | |
531 | T operator()() | |
532 | { | |
533 | T r = result / apn; | |
534 | apn += 1; | |
535 | result *= poch * x / n; | |
536 | ++n; | |
537 | poch += 1; | |
538 | return r; | |
539 | } | |
540 | private: | |
541 | T result, x, apn, poch; | |
542 | int n; | |
543 | }; | |
544 | ||
545 | template <class T, class Lanczos, class Policy> | |
546 | T ibeta_series(T a, T b, T x, T s0, const Lanczos&, bool normalised, T* p_derivative, T y, const Policy& pol) | |
547 | { | |
548 | BOOST_MATH_STD_USING | |
549 | ||
550 | T result; | |
551 | ||
552 | BOOST_ASSERT((p_derivative == 0) || normalised); | |
553 | ||
554 | if(normalised) | |
555 | { | |
556 | T c = a + b; | |
557 | ||
558 | // incomplete beta power term, combined with the Lanczos approximation: | |
559 | T agh = static_cast<T>(a + Lanczos::g() - 0.5f); | |
560 | T bgh = static_cast<T>(b + Lanczos::g() - 0.5f); | |
561 | T cgh = static_cast<T>(c + Lanczos::g() - 0.5f); | |
562 | result = Lanczos::lanczos_sum_expG_scaled(c) / (Lanczos::lanczos_sum_expG_scaled(a) * Lanczos::lanczos_sum_expG_scaled(b)); | |
563 | ||
564 | T l1 = log(cgh / bgh) * (b - 0.5f); | |
565 | T l2 = log(x * cgh / agh) * a; | |
566 | // | |
567 | // Check for over/underflow in the power terms: | |
568 | // | |
569 | if((l1 > tools::log_min_value<T>()) | |
570 | && (l1 < tools::log_max_value<T>()) | |
571 | && (l2 > tools::log_min_value<T>()) | |
572 | && (l2 < tools::log_max_value<T>())) | |
573 | { | |
574 | if(a * b < bgh * 10) | |
575 | result *= exp((b - 0.5f) * boost::math::log1p(a / bgh, pol)); | |
576 | else | |
577 | result *= pow(cgh / bgh, b - 0.5f); | |
578 | result *= pow(x * cgh / agh, a); | |
579 | result *= sqrt(agh / boost::math::constants::e<T>()); | |
580 | ||
581 | if(p_derivative) | |
582 | { | |
583 | *p_derivative = result * pow(y, b); | |
584 | BOOST_ASSERT(*p_derivative >= 0); | |
585 | } | |
586 | } | |
587 | else | |
588 | { | |
589 | // | |
590 | // Oh dear, we need logs, and this *will* cancel: | |
591 | // | |
592 | result = log(result) + l1 + l2 + (log(agh) - 1) / 2; | |
593 | if(p_derivative) | |
594 | *p_derivative = exp(result + b * log(y)); | |
595 | result = exp(result); | |
596 | } | |
597 | } | |
598 | else | |
599 | { | |
600 | // Non-normalised, just compute the power: | |
601 | result = pow(x, a); | |
602 | } | |
603 | if(result < tools::min_value<T>()) | |
604 | return s0; // Safeguard: series can't cope with denorms. | |
605 | ibeta_series_t<T> s(a, b, x, result); | |
606 | boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); | |
607 | result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, s0); | |
608 | policies::check_series_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%) in ibeta_series (with lanczos)", max_iter, pol); | |
609 | return result; | |
610 | } | |
611 | // | |
612 | // Incomplete Beta series again, this time without Lanczos support: | |
613 | // | |
614 | template <class T, class Policy> | |
92f5a8d4 | 615 | T ibeta_series(T a, T b, T x, T s0, const boost::math::lanczos::undefined_lanczos& l, bool normalised, T* p_derivative, T y, const Policy& pol) |
7c673cae FG |
616 | { |
617 | BOOST_MATH_STD_USING | |
618 | ||
619 | T result; | |
620 | BOOST_ASSERT((p_derivative == 0) || normalised); | |
621 | ||
622 | if(normalised) | |
623 | { | |
92f5a8d4 TL |
624 | const T min_sterling = minimum_argument_for_bernoulli_recursion<T>(); |
625 | ||
626 | long shift_a = 0; | |
627 | long shift_b = 0; | |
628 | ||
629 | if (a < min_sterling) | |
630 | shift_a = 1 + ltrunc(min_sterling - a); | |
631 | if (b < min_sterling) | |
632 | shift_b = 1 + ltrunc(min_sterling - b); | |
633 | ||
7c673cae FG |
634 | T c = a + b; |
635 | ||
92f5a8d4 | 636 | if ((shift_a == 0) && (shift_b == 0)) |
7c673cae | 637 | { |
92f5a8d4 | 638 | result = pow(x * c / a, a) * pow(c / b, b) * scaled_tgamma_no_lanczos(c, pol) / (scaled_tgamma_no_lanczos(a, pol) * scaled_tgamma_no_lanczos(b, pol)); |
7c673cae | 639 | } |
92f5a8d4 TL |
640 | else if ((a < 1) && (b > 1)) |
641 | result = pow(x, a) / (boost::math::tgamma(a, pol) * boost::math::tgamma_delta_ratio(b, a, pol)); | |
7c673cae FG |
642 | else |
643 | { | |
92f5a8d4 TL |
644 | T power = pow(x, a); |
645 | T bet = beta_imp(a, b, l, pol); | |
646 | if (!(boost::math::isnormal)(power) || !(boost::math::isnormal)(bet)) | |
647 | { | |
648 | result = exp(a * log(x) + boost::math::lgamma(c, pol) - boost::math::lgamma(a, pol) - boost::math::lgamma(b, pol)); | |
649 | } | |
7c673cae | 650 | else |
92f5a8d4 | 651 | result = power / bet; |
7c673cae | 652 | } |
7c673cae FG |
653 | if(p_derivative) |
654 | { | |
655 | *p_derivative = result * pow(y, b); | |
656 | BOOST_ASSERT(*p_derivative >= 0); | |
657 | } | |
658 | } | |
659 | else | |
660 | { | |
661 | // Non-normalised, just compute the power: | |
662 | result = pow(x, a); | |
663 | } | |
664 | if(result < tools::min_value<T>()) | |
665 | return s0; // Safeguard: series can't cope with denorms. | |
666 | ibeta_series_t<T> s(a, b, x, result); | |
667 | boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); | |
668 | result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, s0); | |
669 | policies::check_series_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%) in ibeta_series (without lanczos)", max_iter, pol); | |
670 | return result; | |
671 | } | |
672 | ||
673 | // | |
674 | // Continued fraction for the incomplete beta: | |
675 | // | |
676 | template <class T> | |
677 | struct ibeta_fraction2_t | |
678 | { | |
679 | typedef std::pair<T, T> result_type; | |
680 | ||
681 | ibeta_fraction2_t(T a_, T b_, T x_, T y_) : a(a_), b(b_), x(x_), y(y_), m(0) {} | |
682 | ||
683 | result_type operator()() | |
684 | { | |
685 | T aN = (a + m - 1) * (a + b + m - 1) * m * (b - m) * x * x; | |
686 | T denom = (a + 2 * m - 1); | |
687 | aN /= denom * denom; | |
688 | ||
689 | T bN = static_cast<T>(m); | |
690 | bN += (m * (b - m) * x) / (a + 2*m - 1); | |
691 | bN += ((a + m) * (a * y - b * x + 1 + m *(2 - x))) / (a + 2*m + 1); | |
692 | ||
693 | ++m; | |
694 | ||
695 | return std::make_pair(aN, bN); | |
696 | } | |
697 | ||
698 | private: | |
699 | T a, b, x, y; | |
700 | int m; | |
701 | }; | |
702 | // | |
703 | // Evaluate the incomplete beta via the continued fraction representation: | |
704 | // | |
705 | template <class T, class Policy> | |
706 | inline T ibeta_fraction2(T a, T b, T x, T y, const Policy& pol, bool normalised, T* p_derivative) | |
707 | { | |
708 | typedef typename lanczos::lanczos<T, Policy>::type lanczos_type; | |
709 | BOOST_MATH_STD_USING | |
710 | T result = ibeta_power_terms(a, b, x, y, lanczos_type(), normalised, pol); | |
711 | if(p_derivative) | |
712 | { | |
713 | *p_derivative = result; | |
714 | BOOST_ASSERT(*p_derivative >= 0); | |
715 | } | |
716 | if(result == 0) | |
717 | return result; | |
718 | ||
719 | ibeta_fraction2_t<T> f(a, b, x, y); | |
720 | T fract = boost::math::tools::continued_fraction_b(f, boost::math::policies::get_epsilon<T, Policy>()); | |
721 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
722 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
723 | return result / fract; | |
724 | } | |
725 | // | |
726 | // Computes the difference between ibeta(a,b,x) and ibeta(a+k,b,x): | |
727 | // | |
728 | template <class T, class Policy> | |
729 | T ibeta_a_step(T a, T b, T x, T y, int k, const Policy& pol, bool normalised, T* p_derivative) | |
730 | { | |
731 | typedef typename lanczos::lanczos<T, Policy>::type lanczos_type; | |
732 | ||
733 | BOOST_MATH_INSTRUMENT_VARIABLE(k); | |
734 | ||
735 | T prefix = ibeta_power_terms(a, b, x, y, lanczos_type(), normalised, pol); | |
736 | if(p_derivative) | |
737 | { | |
738 | *p_derivative = prefix; | |
739 | BOOST_ASSERT(*p_derivative >= 0); | |
740 | } | |
741 | prefix /= a; | |
742 | if(prefix == 0) | |
743 | return prefix; | |
744 | T sum = 1; | |
745 | T term = 1; | |
746 | // series summation from 0 to k-1: | |
747 | for(int i = 0; i < k-1; ++i) | |
748 | { | |
749 | term *= (a+b+i) * x / (a+i+1); | |
750 | sum += term; | |
751 | } | |
752 | prefix *= sum; | |
753 | ||
754 | return prefix; | |
755 | } | |
756 | // | |
757 | // This function is only needed for the non-regular incomplete beta, | |
758 | // it computes the delta in: | |
759 | // beta(a,b,x) = prefix + delta * beta(a+k,b,x) | |
760 | // it is currently only called for small k. | |
761 | // | |
762 | template <class T> | |
763 | inline T rising_factorial_ratio(T a, T b, int k) | |
764 | { | |
765 | // calculate: | |
766 | // (a)(a+1)(a+2)...(a+k-1) | |
767 | // _______________________ | |
768 | // (b)(b+1)(b+2)...(b+k-1) | |
769 | ||
770 | // This is only called with small k, for large k | |
771 | // it is grossly inefficient, do not use outside it's | |
772 | // intended purpose!!! | |
773 | BOOST_MATH_INSTRUMENT_VARIABLE(k); | |
774 | if(k == 0) | |
775 | return 1; | |
776 | T result = 1; | |
777 | for(int i = 0; i < k; ++i) | |
778 | result *= (a+i) / (b+i); | |
779 | return result; | |
780 | } | |
781 | // | |
782 | // Routine for a > 15, b < 1 | |
783 | // | |
784 | // Begin by figuring out how large our table of Pn's should be, | |
f67539c2 | 785 | // quoted accuracies are "guesstimates" based on empirical observation. |
7c673cae FG |
786 | // Note that the table size should never exceed the size of our |
787 | // tables of factorials. | |
788 | // | |
789 | template <class T> | |
790 | struct Pn_size | |
791 | { | |
792 | // This is likely to be enough for ~35-50 digit accuracy | |
793 | // but it's hard to quantify exactly: | |
92f5a8d4 TL |
794 | BOOST_STATIC_CONSTANT(unsigned, value = |
795 | ::boost::math::max_factorial<T>::value >= 100 ? 50 | |
796 | : ::boost::math::max_factorial<T>::value >= ::boost::math::max_factorial<double>::value ? 30 | |
797 | : ::boost::math::max_factorial<T>::value >= ::boost::math::max_factorial<float>::value ? 15 : 1); | |
798 | BOOST_STATIC_ASSERT(::boost::math::max_factorial<T>::value >= ::boost::math::max_factorial<float>::value); | |
7c673cae FG |
799 | }; |
800 | template <> | |
801 | struct Pn_size<float> | |
802 | { | |
803 | BOOST_STATIC_CONSTANT(unsigned, value = 15); // ~8-15 digit accuracy | |
804 | BOOST_STATIC_ASSERT(::boost::math::max_factorial<float>::value >= 30); | |
805 | }; | |
806 | template <> | |
807 | struct Pn_size<double> | |
808 | { | |
809 | BOOST_STATIC_CONSTANT(unsigned, value = 30); // 16-20 digit accuracy | |
810 | BOOST_STATIC_ASSERT(::boost::math::max_factorial<double>::value >= 60); | |
811 | }; | |
812 | template <> | |
813 | struct Pn_size<long double> | |
814 | { | |
815 | BOOST_STATIC_CONSTANT(unsigned, value = 50); // ~35-50 digit accuracy | |
816 | BOOST_STATIC_ASSERT(::boost::math::max_factorial<long double>::value >= 100); | |
817 | }; | |
818 | ||
819 | template <class T, class Policy> | |
820 | T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const Policy& pol, bool normalised) | |
821 | { | |
822 | typedef typename lanczos::lanczos<T, Policy>::type lanczos_type; | |
823 | BOOST_MATH_STD_USING | |
824 | // | |
825 | // This is DiDonato and Morris's BGRAT routine, see Eq's 9 through 9.6. | |
826 | // | |
827 | // Some values we'll need later, these are Eq 9.1: | |
828 | // | |
829 | T bm1 = b - 1; | |
830 | T t = a + bm1 / 2; | |
831 | T lx, u; | |
832 | if(y < 0.35) | |
833 | lx = boost::math::log1p(-y, pol); | |
834 | else | |
835 | lx = log(x); | |
836 | u = -t * lx; | |
837 | // and from from 9.2: | |
838 | T prefix; | |
839 | T h = regularised_gamma_prefix(b, u, pol, lanczos_type()); | |
840 | if(h <= tools::min_value<T>()) | |
841 | return s0; | |
842 | if(normalised) | |
843 | { | |
844 | prefix = h / boost::math::tgamma_delta_ratio(a, b, pol); | |
845 | prefix /= pow(t, b); | |
846 | } | |
847 | else | |
848 | { | |
849 | prefix = full_igamma_prefix(b, u, pol) / pow(t, b); | |
850 | } | |
851 | prefix *= mult; | |
852 | // | |
f67539c2 | 853 | // now we need the quantity Pn, unfortunately this is computed |
7c673cae FG |
854 | // recursively, and requires a full history of all the previous values |
855 | // so no choice but to declare a big table and hope it's big enough... | |
856 | // | |
857 | T p[ ::boost::math::detail::Pn_size<T>::value ] = { 1 }; // see 9.3. | |
858 | // | |
859 | // Now an initial value for J, see 9.6: | |
860 | // | |
861 | T j = boost::math::gamma_q(b, u, pol) / h; | |
862 | // | |
863 | // Now we can start to pull things together and evaluate the sum in Eq 9: | |
864 | // | |
865 | T sum = s0 + prefix * j; // Value at N = 0 | |
866 | // some variables we'll need: | |
867 | unsigned tnp1 = 1; // 2*N+1 | |
868 | T lx2 = lx / 2; | |
869 | lx2 *= lx2; | |
870 | T lxp = 1; | |
871 | T t4 = 4 * t * t; | |
872 | T b2n = b; | |
873 | ||
874 | for(unsigned n = 1; n < sizeof(p)/sizeof(p[0]); ++n) | |
875 | { | |
876 | /* | |
877 | // debugging code, enable this if you want to determine whether | |
878 | // the table of Pn's is large enough... | |
879 | // | |
880 | static int max_count = 2; | |
881 | if(n > max_count) | |
882 | { | |
883 | max_count = n; | |
884 | std::cerr << "Max iterations in BGRAT was " << n << std::endl; | |
885 | } | |
886 | */ | |
887 | // | |
888 | // begin by evaluating the next Pn from Eq 9.4: | |
889 | // | |
890 | tnp1 += 2; | |
891 | p[n] = 0; | |
892 | T mbn = b - n; | |
893 | unsigned tmp1 = 3; | |
894 | for(unsigned m = 1; m < n; ++m) | |
895 | { | |
896 | mbn = m * b - n; | |
897 | p[n] += mbn * p[n-m] / boost::math::unchecked_factorial<T>(tmp1); | |
898 | tmp1 += 2; | |
899 | } | |
900 | p[n] /= n; | |
901 | p[n] += bm1 / boost::math::unchecked_factorial<T>(tnp1); | |
902 | // | |
903 | // Now we want Jn from Jn-1 using Eq 9.6: | |
904 | // | |
905 | j = (b2n * (b2n + 1) * j + (u + b2n + 1) * lxp) / t4; | |
906 | lxp *= lx2; | |
907 | b2n += 2; | |
908 | // | |
909 | // pull it together with Eq 9: | |
910 | // | |
911 | T r = prefix * p[n] * j; | |
912 | sum += r; | |
913 | if(r > 1) | |
914 | { | |
915 | if(fabs(r) < fabs(tools::epsilon<T>() * sum)) | |
916 | break; | |
917 | } | |
918 | else | |
919 | { | |
920 | if(fabs(r / tools::epsilon<T>()) < fabs(sum)) | |
921 | break; | |
922 | } | |
923 | } | |
924 | return sum; | |
925 | } // template <class T, class Lanczos>T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const Lanczos& l, bool normalised) | |
926 | ||
927 | // | |
928 | // For integer arguments we can relate the incomplete beta to the | |
929 | // complement of the binomial distribution cdf and use this finite sum. | |
930 | // | |
931 | template <class T> | |
932 | T binomial_ccdf(T n, T k, T x, T y) | |
933 | { | |
934 | BOOST_MATH_STD_USING // ADL of std names | |
935 | ||
936 | T result = pow(x, n); | |
937 | ||
938 | if(result > tools::min_value<T>()) | |
939 | { | |
940 | T term = result; | |
941 | for(unsigned i = itrunc(T(n - 1)); i > k; --i) | |
942 | { | |
943 | term *= ((i + 1) * y) / ((n - i) * x); | |
944 | result += term; | |
945 | } | |
946 | } | |
947 | else | |
948 | { | |
949 | // First term underflows so we need to start at the mode of the | |
950 | // distribution and work outwards: | |
951 | int start = itrunc(n * x); | |
952 | if(start <= k + 1) | |
953 | start = itrunc(k + 2); | |
954 | result = pow(x, start) * pow(y, n - start) * boost::math::binomial_coefficient<T>(itrunc(n), itrunc(start)); | |
955 | if(result == 0) | |
956 | { | |
92f5a8d4 | 957 | // OK, starting slightly above the mode didn't work, |
7c673cae FG |
958 | // we'll have to sum the terms the old fashioned way: |
959 | for(unsigned i = start - 1; i > k; --i) | |
960 | { | |
961 | result += pow(x, (int)i) * pow(y, n - i) * boost::math::binomial_coefficient<T>(itrunc(n), itrunc(i)); | |
962 | } | |
963 | } | |
964 | else | |
965 | { | |
966 | T term = result; | |
967 | T start_term = result; | |
968 | for(unsigned i = start - 1; i > k; --i) | |
969 | { | |
970 | term *= ((i + 1) * y) / ((n - i) * x); | |
971 | result += term; | |
972 | } | |
973 | term = start_term; | |
974 | for(unsigned i = start + 1; i <= n; ++i) | |
975 | { | |
976 | term *= (n - i + 1) * x / (i * y); | |
977 | result += term; | |
978 | } | |
979 | } | |
980 | } | |
981 | ||
982 | return result; | |
983 | } | |
984 | ||
985 | ||
986 | // | |
987 | // The incomplete beta function implementation: | |
f67539c2 | 988 | // This is just a big bunch of spaghetti code to divide up the |
7c673cae FG |
989 | // input range and select the right implementation method for |
990 | // each domain: | |
991 | // | |
992 | template <class T, class Policy> | |
993 | T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_derivative) | |
994 | { | |
995 | static const char* function = "boost::math::ibeta<%1%>(%1%, %1%, %1%)"; | |
996 | typedef typename lanczos::lanczos<T, Policy>::type lanczos_type; | |
997 | BOOST_MATH_STD_USING // for ADL of std math functions. | |
998 | ||
999 | BOOST_MATH_INSTRUMENT_VARIABLE(a); | |
1000 | BOOST_MATH_INSTRUMENT_VARIABLE(b); | |
1001 | BOOST_MATH_INSTRUMENT_VARIABLE(x); | |
1002 | BOOST_MATH_INSTRUMENT_VARIABLE(inv); | |
1003 | BOOST_MATH_INSTRUMENT_VARIABLE(normalised); | |
1004 | ||
1005 | bool invert = inv; | |
1006 | T fract; | |
1007 | T y = 1 - x; | |
1008 | ||
1009 | BOOST_ASSERT((p_derivative == 0) || normalised); | |
1010 | ||
1011 | if(p_derivative) | |
1012 | *p_derivative = -1; // value not set. | |
1013 | ||
1014 | if((x < 0) || (x > 1)) | |
1015 | return policies::raise_domain_error<T>(function, "Parameter x outside the range [0,1] in the incomplete beta function (got x=%1%).", x, pol); | |
1016 | ||
1017 | if(normalised) | |
1018 | { | |
1019 | if(a < 0) | |
1020 | return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be >= zero (got a=%1%).", a, pol); | |
1021 | if(b < 0) | |
1022 | return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be >= zero (got b=%1%).", b, pol); | |
1023 | // extend to a few very special cases: | |
1024 | if(a == 0) | |
1025 | { | |
1026 | if(b == 0) | |
1027 | return policies::raise_domain_error<T>(function, "The arguments a and b to the incomplete beta function cannot both be zero, with x=%1%.", x, pol); | |
1028 | if(b > 0) | |
1029 | return static_cast<T>(inv ? 0 : 1); | |
1030 | } | |
1031 | else if(b == 0) | |
1032 | { | |
1033 | if(a > 0) | |
1034 | return static_cast<T>(inv ? 1 : 0); | |
1035 | } | |
1036 | } | |
1037 | else | |
1038 | { | |
1039 | if(a <= 0) | |
1040 | return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be greater than zero (got a=%1%).", a, pol); | |
1041 | if(b <= 0) | |
1042 | return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be greater than zero (got b=%1%).", b, pol); | |
1043 | } | |
1044 | ||
1045 | if(x == 0) | |
1046 | { | |
1047 | if(p_derivative) | |
1048 | { | |
1049 | *p_derivative = (a == 1) ? (T)1 : (a < 1) ? T(tools::max_value<T>() / 2) : T(tools::min_value<T>() * 2); | |
1050 | } | |
1051 | return (invert ? (normalised ? T(1) : boost::math::beta(a, b, pol)) : T(0)); | |
1052 | } | |
1053 | if(x == 1) | |
1054 | { | |
1055 | if(p_derivative) | |
1056 | { | |
1057 | *p_derivative = (b == 1) ? T(1) : (b < 1) ? T(tools::max_value<T>() / 2) : T(tools::min_value<T>() * 2); | |
1058 | } | |
1059 | return (invert == 0 ? (normalised ? 1 : boost::math::beta(a, b, pol)) : 0); | |
1060 | } | |
1061 | if((a == 0.5f) && (b == 0.5f)) | |
1062 | { | |
1063 | // We have an arcsine distribution: | |
1064 | if(p_derivative) | |
1065 | { | |
1066 | *p_derivative = 1 / constants::pi<T>() * sqrt(y * x); | |
1067 | } | |
1068 | T p = invert ? asin(sqrt(y)) / constants::half_pi<T>() : asin(sqrt(x)) / constants::half_pi<T>(); | |
1069 | if(!normalised) | |
1070 | p *= constants::pi<T>(); | |
1071 | return p; | |
1072 | } | |
1073 | if(a == 1) | |
1074 | { | |
1075 | std::swap(a, b); | |
1076 | std::swap(x, y); | |
1077 | invert = !invert; | |
1078 | } | |
1079 | if(b == 1) | |
1080 | { | |
1081 | // | |
1082 | // Special case see: http://functions.wolfram.com/GammaBetaErf/BetaRegularized/03/01/01/ | |
1083 | // | |
1084 | if(a == 1) | |
1085 | { | |
1086 | if(p_derivative) | |
1087 | *p_derivative = 1; | |
1088 | return invert ? y : x; | |
1089 | } | |
92f5a8d4 | 1090 | |
7c673cae FG |
1091 | if(p_derivative) |
1092 | { | |
1093 | *p_derivative = a * pow(x, a - 1); | |
1094 | } | |
1095 | T p; | |
1096 | if(y < 0.5) | |
1097 | p = invert ? T(-boost::math::expm1(a * boost::math::log1p(-y, pol), pol)) : T(exp(a * boost::math::log1p(-y, pol))); | |
1098 | else | |
1099 | p = invert ? T(-boost::math::powm1(x, a, pol)) : T(pow(x, a)); | |
1100 | if(!normalised) | |
1101 | p /= a; | |
1102 | return p; | |
1103 | } | |
1104 | ||
1105 | if((std::min)(a, b) <= 1) | |
1106 | { | |
1107 | if(x > 0.5) | |
1108 | { | |
1109 | std::swap(a, b); | |
1110 | std::swap(x, y); | |
1111 | invert = !invert; | |
1112 | BOOST_MATH_INSTRUMENT_VARIABLE(invert); | |
1113 | } | |
1114 | if((std::max)(a, b) <= 1) | |
1115 | { | |
1116 | // Both a,b < 1: | |
1117 | if((a >= (std::min)(T(0.2), b)) || (pow(x, a) <= 0.9)) | |
1118 | { | |
1119 | if(!invert) | |
1120 | { | |
1121 | fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol); | |
1122 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1123 | } | |
1124 | else | |
1125 | { | |
1126 | fract = -(normalised ? 1 : boost::math::beta(a, b, pol)); | |
1127 | invert = false; | |
1128 | fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol); | |
1129 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1130 | } | |
1131 | } | |
1132 | else | |
1133 | { | |
1134 | std::swap(a, b); | |
1135 | std::swap(x, y); | |
1136 | invert = !invert; | |
1137 | if(y >= 0.3) | |
1138 | { | |
1139 | if(!invert) | |
1140 | { | |
1141 | fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol); | |
1142 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1143 | } | |
1144 | else | |
1145 | { | |
1146 | fract = -(normalised ? 1 : boost::math::beta(a, b, pol)); | |
1147 | invert = false; | |
1148 | fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol); | |
1149 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1150 | } | |
1151 | } | |
1152 | else | |
1153 | { | |
1154 | // Sidestep on a, and then use the series representation: | |
1155 | T prefix; | |
1156 | if(!normalised) | |
1157 | { | |
1158 | prefix = rising_factorial_ratio(T(a+b), a, 20); | |
1159 | } | |
1160 | else | |
1161 | { | |
1162 | prefix = 1; | |
1163 | } | |
1164 | fract = ibeta_a_step(a, b, x, y, 20, pol, normalised, p_derivative); | |
1165 | if(!invert) | |
1166 | { | |
1167 | fract = beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised); | |
1168 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1169 | } | |
1170 | else | |
1171 | { | |
1172 | fract -= (normalised ? 1 : boost::math::beta(a, b, pol)); | |
1173 | invert = false; | |
1174 | fract = -beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised); | |
1175 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1176 | } | |
1177 | } | |
1178 | } | |
1179 | } | |
1180 | else | |
1181 | { | |
1182 | // One of a, b < 1 only: | |
1183 | if((b <= 1) || ((x < 0.1) && (pow(b * x, a) <= 0.7))) | |
1184 | { | |
1185 | if(!invert) | |
1186 | { | |
1187 | fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol); | |
1188 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1189 | } | |
1190 | else | |
1191 | { | |
1192 | fract = -(normalised ? 1 : boost::math::beta(a, b, pol)); | |
1193 | invert = false; | |
1194 | fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol); | |
1195 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1196 | } | |
1197 | } | |
1198 | else | |
1199 | { | |
1200 | std::swap(a, b); | |
1201 | std::swap(x, y); | |
1202 | invert = !invert; | |
1203 | ||
1204 | if(y >= 0.3) | |
1205 | { | |
1206 | if(!invert) | |
1207 | { | |
1208 | fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol); | |
1209 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1210 | } | |
1211 | else | |
1212 | { | |
1213 | fract = -(normalised ? 1 : boost::math::beta(a, b, pol)); | |
1214 | invert = false; | |
1215 | fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol); | |
1216 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1217 | } | |
1218 | } | |
1219 | else if(a >= 15) | |
1220 | { | |
1221 | if(!invert) | |
1222 | { | |
1223 | fract = beta_small_b_large_a_series(a, b, x, y, T(0), T(1), pol, normalised); | |
1224 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1225 | } | |
1226 | else | |
1227 | { | |
1228 | fract = -(normalised ? 1 : boost::math::beta(a, b, pol)); | |
1229 | invert = false; | |
1230 | fract = -beta_small_b_large_a_series(a, b, x, y, fract, T(1), pol, normalised); | |
1231 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1232 | } | |
1233 | } | |
1234 | else | |
1235 | { | |
1236 | // Sidestep to improve errors: | |
1237 | T prefix; | |
1238 | if(!normalised) | |
1239 | { | |
1240 | prefix = rising_factorial_ratio(T(a+b), a, 20); | |
1241 | } | |
1242 | else | |
1243 | { | |
1244 | prefix = 1; | |
1245 | } | |
1246 | fract = ibeta_a_step(a, b, x, y, 20, pol, normalised, p_derivative); | |
1247 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1248 | if(!invert) | |
1249 | { | |
1250 | fract = beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised); | |
1251 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1252 | } | |
1253 | else | |
1254 | { | |
1255 | fract -= (normalised ? 1 : boost::math::beta(a, b, pol)); | |
1256 | invert = false; | |
1257 | fract = -beta_small_b_large_a_series(T(a + 20), b, x, y, fract, prefix, pol, normalised); | |
1258 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1259 | } | |
1260 | } | |
1261 | } | |
1262 | } | |
1263 | } | |
1264 | else | |
1265 | { | |
1266 | // Both a,b >= 1: | |
1267 | T lambda; | |
1268 | if(a < b) | |
1269 | { | |
1270 | lambda = a - (a + b) * x; | |
1271 | } | |
1272 | else | |
1273 | { | |
1274 | lambda = (a + b) * y - b; | |
1275 | } | |
1276 | if(lambda < 0) | |
1277 | { | |
1278 | std::swap(a, b); | |
1279 | std::swap(x, y); | |
1280 | invert = !invert; | |
1281 | BOOST_MATH_INSTRUMENT_VARIABLE(invert); | |
1282 | } | |
92f5a8d4 | 1283 | |
7c673cae FG |
1284 | if(b < 40) |
1285 | { | |
1286 | if((floor(a) == a) && (floor(b) == b) && (a < (std::numeric_limits<int>::max)() - 100) && (y != 1)) | |
1287 | { | |
1288 | // relate to the binomial distribution and use a finite sum: | |
1289 | T k = a - 1; | |
1290 | T n = b + k; | |
1291 | fract = binomial_ccdf(n, k, x, y); | |
1292 | if(!normalised) | |
1293 | fract *= boost::math::beta(a, b, pol); | |
1294 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1295 | } | |
1296 | else if(b * x <= 0.7) | |
1297 | { | |
1298 | if(!invert) | |
1299 | { | |
1300 | fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol); | |
1301 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1302 | } | |
1303 | else | |
1304 | { | |
1305 | fract = -(normalised ? 1 : boost::math::beta(a, b, pol)); | |
1306 | invert = false; | |
1307 | fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol); | |
1308 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1309 | } | |
1310 | } | |
1311 | else if(a > 15) | |
1312 | { | |
1313 | // sidestep so we can use the series representation: | |
1314 | int n = itrunc(T(floor(b)), pol); | |
1315 | if(n == b) | |
1316 | --n; | |
1317 | T bbar = b - n; | |
1318 | T prefix; | |
1319 | if(!normalised) | |
1320 | { | |
1321 | prefix = rising_factorial_ratio(T(a+bbar), bbar, n); | |
1322 | } | |
1323 | else | |
1324 | { | |
1325 | prefix = 1; | |
1326 | } | |
1327 | fract = ibeta_a_step(bbar, a, y, x, n, pol, normalised, static_cast<T*>(0)); | |
1328 | fract = beta_small_b_large_a_series(a, bbar, x, y, fract, T(1), pol, normalised); | |
1329 | fract /= prefix; | |
1330 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1331 | } | |
1332 | else if(normalised) | |
1333 | { | |
1334 | // The formula here for the non-normalised case is tricky to figure | |
1335 | // out (for me!!), and requires two pochhammer calculations rather | |
1336 | // than one, so leave it for now and only use this in the normalized case.... | |
1337 | int n = itrunc(T(floor(b)), pol); | |
1338 | T bbar = b - n; | |
1339 | if(bbar <= 0) | |
1340 | { | |
1341 | --n; | |
1342 | bbar += 1; | |
1343 | } | |
1344 | fract = ibeta_a_step(bbar, a, y, x, n, pol, normalised, static_cast<T*>(0)); | |
1345 | fract += ibeta_a_step(a, bbar, x, y, 20, pol, normalised, static_cast<T*>(0)); | |
1346 | if(invert) | |
1347 | fract -= 1; // Note this line would need changing if we ever enable this branch in non-normalized case | |
1348 | fract = beta_small_b_large_a_series(T(a+20), bbar, x, y, fract, T(1), pol, normalised); | |
1349 | if(invert) | |
1350 | { | |
1351 | fract = -fract; | |
1352 | invert = false; | |
1353 | } | |
1354 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1355 | } | |
1356 | else | |
1357 | { | |
1358 | fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative); | |
1359 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1360 | } | |
1361 | } | |
1362 | else | |
1363 | { | |
1364 | fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative); | |
1365 | BOOST_MATH_INSTRUMENT_VARIABLE(fract); | |
1366 | } | |
1367 | } | |
1368 | if(p_derivative) | |
1369 | { | |
1370 | if(*p_derivative < 0) | |
1371 | { | |
1372 | *p_derivative = ibeta_power_terms(a, b, x, y, lanczos_type(), true, pol); | |
1373 | } | |
1374 | T div = y * x; | |
1375 | ||
1376 | if(*p_derivative != 0) | |
1377 | { | |
1378 | if((tools::max_value<T>() * div < *p_derivative)) | |
1379 | { | |
f67539c2 | 1380 | // overflow, return an arbitrarily large value: |
7c673cae FG |
1381 | *p_derivative = tools::max_value<T>() / 2; |
1382 | } | |
1383 | else | |
1384 | { | |
1385 | *p_derivative /= div; | |
1386 | } | |
1387 | } | |
1388 | } | |
1389 | return invert ? (normalised ? 1 : boost::math::beta(a, b, pol)) - fract : fract; | |
1390 | } // template <class T, class Lanczos>T ibeta_imp(T a, T b, T x, const Lanczos& l, bool inv, bool normalised) | |
1391 | ||
1392 | template <class T, class Policy> | |
1393 | inline T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised) | |
1394 | { | |
1395 | return ibeta_imp(a, b, x, pol, inv, normalised, static_cast<T*>(0)); | |
1396 | } | |
1397 | ||
1398 | template <class T, class Policy> | |
1399 | T ibeta_derivative_imp(T a, T b, T x, const Policy& pol) | |
1400 | { | |
1401 | static const char* function = "ibeta_derivative<%1%>(%1%,%1%,%1%)"; | |
1402 | // | |
1403 | // start with the usual error checks: | |
1404 | // | |
1405 | if(a <= 0) | |
1406 | return policies::raise_domain_error<T>(function, "The argument a to the incomplete beta function must be greater than zero (got a=%1%).", a, pol); | |
1407 | if(b <= 0) | |
1408 | return policies::raise_domain_error<T>(function, "The argument b to the incomplete beta function must be greater than zero (got b=%1%).", b, pol); | |
1409 | if((x < 0) || (x > 1)) | |
1410 | return policies::raise_domain_error<T>(function, "Parameter x outside the range [0,1] in the incomplete beta function (got x=%1%).", x, pol); | |
1411 | // | |
1412 | // Now the corner cases: | |
1413 | // | |
1414 | if(x == 0) | |
1415 | { | |
92f5a8d4 | 1416 | return (a > 1) ? 0 : |
7c673cae FG |
1417 | (a == 1) ? 1 / boost::math::beta(a, b, pol) : policies::raise_overflow_error<T>(function, 0, pol); |
1418 | } | |
1419 | else if(x == 1) | |
1420 | { | |
1421 | return (b > 1) ? 0 : | |
1422 | (b == 1) ? 1 / boost::math::beta(a, b, pol) : policies::raise_overflow_error<T>(function, 0, pol); | |
1423 | } | |
1424 | // | |
1425 | // Now the regular cases: | |
1426 | // | |
1427 | typedef typename lanczos::lanczos<T, Policy>::type lanczos_type; | |
1428 | T y = (1 - x) * x; | |
1429 | T f1 = ibeta_power_terms<T>(a, b, x, 1 - x, lanczos_type(), true, pol, 1 / y, function); | |
1430 | return f1; | |
1431 | } | |
1432 | // | |
1433 | // Some forwarding functions that dis-ambiguate the third argument type: | |
1434 | // | |
1435 | template <class RT1, class RT2, class Policy> | |
92f5a8d4 | 1436 | inline typename tools::promote_args<RT1, RT2>::type |
f67539c2 | 1437 | beta(RT1 a, RT2 b, const Policy&, const boost::true_type*) |
7c673cae FG |
1438 | { |
1439 | BOOST_FPU_EXCEPTION_GUARD | |
1440 | typedef typename tools::promote_args<RT1, RT2>::type result_type; | |
1441 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
1442 | typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type; | |
1443 | typedef typename policies::normalise< | |
92f5a8d4 TL |
1444 | Policy, |
1445 | policies::promote_float<false>, | |
1446 | policies::promote_double<false>, | |
7c673cae FG |
1447 | policies::discrete_quantile<>, |
1448 | policies::assert_undefined<> >::type forwarding_policy; | |
1449 | ||
1450 | return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::beta_imp(static_cast<value_type>(a), static_cast<value_type>(b), evaluation_type(), forwarding_policy()), "boost::math::beta<%1%>(%1%,%1%)"); | |
1451 | } | |
1452 | template <class RT1, class RT2, class RT3> | |
92f5a8d4 | 1453 | inline typename tools::promote_args<RT1, RT2, RT3>::type |
f67539c2 | 1454 | beta(RT1 a, RT2 b, RT3 x, const boost::false_type*) |
7c673cae FG |
1455 | { |
1456 | return boost::math::beta(a, b, x, policies::policy<>()); | |
1457 | } | |
1458 | } // namespace detail | |
1459 | ||
1460 | // | |
1461 | // The actual function entry-points now follow, these just figure out | |
1462 | // which Lanczos approximation to use | |
1463 | // and forward to the implementation functions: | |
1464 | // | |
1465 | template <class RT1, class RT2, class A> | |
92f5a8d4 | 1466 | inline typename tools::promote_args<RT1, RT2, A>::type |
7c673cae FG |
1467 | beta(RT1 a, RT2 b, A arg) |
1468 | { | |
1469 | typedef typename policies::is_policy<A>::type tag; | |
1470 | return boost::math::detail::beta(a, b, arg, static_cast<tag*>(0)); | |
1471 | } | |
1472 | ||
1473 | template <class RT1, class RT2> | |
92f5a8d4 | 1474 | inline typename tools::promote_args<RT1, RT2>::type |
7c673cae FG |
1475 | beta(RT1 a, RT2 b) |
1476 | { | |
1477 | return boost::math::beta(a, b, policies::policy<>()); | |
1478 | } | |
1479 | ||
1480 | template <class RT1, class RT2, class RT3, class Policy> | |
92f5a8d4 | 1481 | inline typename tools::promote_args<RT1, RT2, RT3>::type |
7c673cae FG |
1482 | beta(RT1 a, RT2 b, RT3 x, const Policy&) |
1483 | { | |
1484 | BOOST_FPU_EXCEPTION_GUARD | |
1485 | typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type; | |
1486 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
1487 | typedef typename policies::normalise< | |
92f5a8d4 TL |
1488 | Policy, |
1489 | policies::promote_float<false>, | |
1490 | policies::promote_double<false>, | |
7c673cae FG |
1491 | policies::discrete_quantile<>, |
1492 | policies::assert_undefined<> >::type forwarding_policy; | |
1493 | ||
1494 | return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), false, false), "boost::math::beta<%1%>(%1%,%1%,%1%)"); | |
1495 | } | |
1496 | ||
1497 | template <class RT1, class RT2, class RT3, class Policy> | |
92f5a8d4 | 1498 | inline typename tools::promote_args<RT1, RT2, RT3>::type |
7c673cae FG |
1499 | betac(RT1 a, RT2 b, RT3 x, const Policy&) |
1500 | { | |
1501 | BOOST_FPU_EXCEPTION_GUARD | |
1502 | typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type; | |
1503 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
1504 | typedef typename policies::normalise< | |
92f5a8d4 TL |
1505 | Policy, |
1506 | policies::promote_float<false>, | |
1507 | policies::promote_double<false>, | |
7c673cae FG |
1508 | policies::discrete_quantile<>, |
1509 | policies::assert_undefined<> >::type forwarding_policy; | |
1510 | ||
1511 | return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), true, false), "boost::math::betac<%1%>(%1%,%1%,%1%)"); | |
1512 | } | |
1513 | template <class RT1, class RT2, class RT3> | |
92f5a8d4 | 1514 | inline typename tools::promote_args<RT1, RT2, RT3>::type |
7c673cae FG |
1515 | betac(RT1 a, RT2 b, RT3 x) |
1516 | { | |
1517 | return boost::math::betac(a, b, x, policies::policy<>()); | |
1518 | } | |
1519 | ||
1520 | template <class RT1, class RT2, class RT3, class Policy> | |
92f5a8d4 | 1521 | inline typename tools::promote_args<RT1, RT2, RT3>::type |
7c673cae FG |
1522 | ibeta(RT1 a, RT2 b, RT3 x, const Policy&) |
1523 | { | |
1524 | BOOST_FPU_EXCEPTION_GUARD | |
1525 | typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type; | |
1526 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
1527 | typedef typename policies::normalise< | |
92f5a8d4 TL |
1528 | Policy, |
1529 | policies::promote_float<false>, | |
1530 | policies::promote_double<false>, | |
7c673cae FG |
1531 | policies::discrete_quantile<>, |
1532 | policies::assert_undefined<> >::type forwarding_policy; | |
1533 | ||
1534 | return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), false, true), "boost::math::ibeta<%1%>(%1%,%1%,%1%)"); | |
1535 | } | |
1536 | template <class RT1, class RT2, class RT3> | |
92f5a8d4 | 1537 | inline typename tools::promote_args<RT1, RT2, RT3>::type |
7c673cae FG |
1538 | ibeta(RT1 a, RT2 b, RT3 x) |
1539 | { | |
1540 | return boost::math::ibeta(a, b, x, policies::policy<>()); | |
1541 | } | |
1542 | ||
1543 | template <class RT1, class RT2, class RT3, class Policy> | |
92f5a8d4 | 1544 | inline typename tools::promote_args<RT1, RT2, RT3>::type |
7c673cae FG |
1545 | ibetac(RT1 a, RT2 b, RT3 x, const Policy&) |
1546 | { | |
1547 | BOOST_FPU_EXCEPTION_GUARD | |
1548 | typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type; | |
1549 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
1550 | typedef typename policies::normalise< | |
92f5a8d4 TL |
1551 | Policy, |
1552 | policies::promote_float<false>, | |
1553 | policies::promote_double<false>, | |
7c673cae FG |
1554 | policies::discrete_quantile<>, |
1555 | policies::assert_undefined<> >::type forwarding_policy; | |
1556 | ||
1557 | return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy(), true, true), "boost::math::ibetac<%1%>(%1%,%1%,%1%)"); | |
1558 | } | |
1559 | template <class RT1, class RT2, class RT3> | |
92f5a8d4 | 1560 | inline typename tools::promote_args<RT1, RT2, RT3>::type |
7c673cae FG |
1561 | ibetac(RT1 a, RT2 b, RT3 x) |
1562 | { | |
1563 | return boost::math::ibetac(a, b, x, policies::policy<>()); | |
1564 | } | |
1565 | ||
1566 | template <class RT1, class RT2, class RT3, class Policy> | |
92f5a8d4 | 1567 | inline typename tools::promote_args<RT1, RT2, RT3>::type |
7c673cae FG |
1568 | ibeta_derivative(RT1 a, RT2 b, RT3 x, const Policy&) |
1569 | { | |
1570 | BOOST_FPU_EXCEPTION_GUARD | |
1571 | typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type; | |
1572 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
1573 | typedef typename policies::normalise< | |
92f5a8d4 TL |
1574 | Policy, |
1575 | policies::promote_float<false>, | |
1576 | policies::promote_double<false>, | |
7c673cae FG |
1577 | policies::discrete_quantile<>, |
1578 | policies::assert_undefined<> >::type forwarding_policy; | |
1579 | ||
1580 | return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::ibeta_derivative_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), forwarding_policy()), "boost::math::ibeta_derivative<%1%>(%1%,%1%,%1%)"); | |
1581 | } | |
1582 | template <class RT1, class RT2, class RT3> | |
92f5a8d4 | 1583 | inline typename tools::promote_args<RT1, RT2, RT3>::type |
7c673cae FG |
1584 | ibeta_derivative(RT1 a, RT2 b, RT3 x) |
1585 | { | |
1586 | return boost::math::ibeta_derivative(a, b, x, policies::policy<>()); | |
1587 | } | |
1588 | ||
1589 | } // namespace math | |
1590 | } // namespace boost | |
1591 | ||
1592 | #include <boost/math/special_functions/detail/ibeta_inverse.hpp> | |
1593 | #include <boost/math/special_functions/detail/ibeta_inv_ab.hpp> | |
1594 | ||
1595 | #endif // BOOST_MATH_SPECIAL_BETA_HPP |