]>
git.proxmox.com Git - ceph.git/blob - ceph/src/boost/libs/math/minimax/f.cpp
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)
7 //#include "../tools/ntl_rr_lanczos.hpp"
8 //#include "../tools/ntl_rr_digamma.hpp"
9 #include "multiprecision.hpp"
10 #include <boost/math/tools/polynomial.hpp>
11 #include <boost/math/special_functions.hpp>
12 #include <boost/math/special_functions/zeta.hpp>
13 #include <boost/math/special_functions/expint.hpp>
14 #include <boost/math/special_functions/lambert_w.hpp>
19 mp_type
f(const mp_type
& x
, int variant
)
21 static const mp_type tiny
= boost::math::tools::min_value
<mp_type
>() * 64;
26 mp_type x_
= sqrt(x
== 0 ? 1e-80 : x
);
27 return boost::math::erf(x_
) / x_
;
32 return boost::math::erfc(x_
) * x_
/ exp(-x_
* x_
);
36 return boost::math::erfc(x
) * x
/ exp(-x
* x
);
43 return boost::math::lgamma(y
+2) / y
- 0.5;
47 // lgamma in the range [2,3], use:
49 // lgamma(x) = (x-2) * (x + 1) * (c + R(x - 2))
51 // Works well at 80-bit long double precision, but doesn't
52 // stretch to 128-bit precision.
56 return boost::lexical_cast
<mp_type
>("0.42278433509846713939348790991759756895784066406008") / 3;
58 return boost::math::lgamma(x
+2) / (x
* (x
+3));
62 // lgamma in the range [1,2], use:
64 // lgamma(x) = (x - 1) * (x - 2) * (c + R(x - 1))
66 // works well over [1, 1.5] but not near 2 :-(
68 mp_type r1
= boost::lexical_cast
<mp_type
>("0.57721566490153286060651209008240243104215933593992");
69 mp_type r2
= boost::lexical_cast
<mp_type
>("0.42278433509846713939348790991759756895784066406008");
78 return boost::math::lgamma(x
+1) / (x
* (x
- 1));
83 // lgamma in the range [1.5,2], use:
85 // lgamma(x) = (2 - x) * (1 - x) * (c + R(2 - x))
87 // works well over [1.5, 2] but not near 1 :-(
89 mp_type r1
= boost::lexical_cast
<mp_type
>("0.57721566490153286060651209008240243104215933593992");
90 mp_type r2
= boost::lexical_cast
<mp_type
>("0.42278433509846713939348790991759756895784066406008");
99 return boost::math::lgamma(2-x
) / (x
* (x
- 1));
104 // erf_inv in range [0, 0.5]
108 y
= boost::math::tools::epsilon
<mp_type
>() / 64;
109 return boost::math::erf_inv(y
) / (y
* (y
+10));
114 // erfc_inv in range [0.25, 0.5]
115 // Use an y-offset of 0.25, and range [0, 0.25]
116 // abs error, auto y-offset.
120 y
= boost::lexical_cast
<mp_type
>("1e-5000");
121 return sqrt(-2 * log(y
)) / boost::math::erfc_inv(y
);
127 x2
= boost::lexical_cast
<mp_type
>("1e-5000");
128 mp_type y
= exp(-x2
*x2
); // sqrt(-log(x2)) - 5;
129 return boost::math::erfc_inv(y
) / x2
;
134 // Digamma over the interval [1,2], set x-offset to 1
135 // and optimise for absolute error over [0,1].
137 int current_precision
= get_working_precision();
138 if(current_precision
< 1000)
139 set_working_precision(1000);
141 // This value for the root of digamma is calculated using our
142 // differentiated lanczos approximation. It agrees with Cody
143 // to ~ 25 digits and to Morris to 35 digits. See:
144 // TOMS ALGORITHM 708 (Didonato and Morris).
145 // and Math. Comp. 27, 123-127 (1973) by Cody, Strecok and Thacher.
147 //mp_type root = boost::lexical_cast<mp_type>("1.4616321449683623412626595423257213234331845807102825466429633351908372838889871");
149 // Actually better to calculate the root on the fly, it appears to be more
150 // accurate: convergence is easier with the 1000-bit value, the approximation
151 // produced agrees with functions.mathworld.com values to 35 digits even quite
154 static boost::math::tools::eps_tolerance
<mp_type
> tol(1000);
155 static boost::uintmax_t max_iter
= 1000;
156 mp_type (*pdg
)(mp_type
) = &boost::math::digamma
;
157 static const mp_type root
= boost::math::tools::bracket_and_solve_root(pdg
, mp_type(1.4), mp_type(1.5), true, tol
, max_iter
).first
;
161 if(fabs(x2
- root
) < lim
)
164 // This is a problem area:
165 // x2-root suffers cancellation error, so does digamma.
166 // That gets compounded again when Remez calculates the error
167 // function. This cludge seems to stop the worst of the problems:
169 static const mp_type a
= boost::math::digamma(root
- lim
) / -lim
;
170 static const mp_type b
= boost::math::digamma(root
+ lim
) / lim
;
171 mp_type fract
= (x2
- root
+ lim
) / (2*lim
);
172 mp_type r
= (1-fract
) * a
+ fract
* b
;
173 std::cout
<< "In root area: " << r
;
176 mp_type result
= boost::math::digamma(x2
) / (x2
- root
);
177 if(current_precision
< 1000)
178 set_working_precision(current_precision
);
185 static mp_type lim
= 1e-80;
186 static mp_type a
= boost::math::expm1(-lim
);
187 static mp_type b
= boost::math::expm1(lim
);
188 static mp_type l
= (b
-a
) / (2 * lim
);
191 return boost::math::expm1(x
) / x
;
193 // demo, and test case:
198 return boost::math::ellint_1(x
);
203 return boost::math::ellint_1(1-x
) / log(x
);
209 mp_type z
= 1 - x
* log(x
);
210 return boost::math::ellint_2(sqrt(1-x
)) / z
;
213 // Bessel I0(x) over [0,16]:
215 return boost::math::cyl_bessel_i(0, sqrt(x
));
218 // Bessel I0(x) over [16,INF]
220 mp_type z
= 1 / (mp_type(1)/16 - x
);
221 return boost::math::cyl_bessel_i(0, z
) * sqrt(z
) / exp(z
);
226 return boost::math::zeta(1 - x
) * x
- x
;
231 return boost::math::zeta(x
) - 1 / (x
- 1);
234 // Zeta over [a, b] : a >> 1
236 return log(boost::math::zeta(x
) - 1);
239 // expint[1] over [0,1]:
241 mp_type tiny
= boost::lexical_cast
<mp_type
>("1e-5000");
242 mp_type z
= (x
<= tiny
) ? tiny
: x
;
243 return boost::math::expint(1, z
) - z
+ log(z
);
246 // expint[1] over [1,N],
247 // Note that x varies from [0,1]:
250 return boost::math::expint(1, z
) * exp(z
) * z
;
253 // expin Ei over [0,R]
255 static const mp_type root
=
256 boost::lexical_cast
<mp_type
>("0.372507410781366634461991866580119133535689497771654051555657435242200120636201854384926049951548942392");
257 mp_type z
= x
< (std::numeric_limits
<long double>::min
)() ? (std::numeric_limits
<long double>::min
)() : x
;
258 return (boost::math::expint(z
) - log(z
/ root
)) / (z
- root
);
261 // Expint Ei for large x:
263 static const mp_type root
=
264 boost::lexical_cast
<mp_type
>("0.372507410781366634461991866580119133535689497771654051555657435242200120636201854384926049951548942392");
265 mp_type z
= x
< (std::numeric_limits
<long double>::min
)() ? (std::numeric_limits
<long double>::max
)() : mp_type(1 / x
);
266 return (boost::math::expint(z
) - z
) * z
* exp(-z
);
267 //return (boost::math::expint(z) - log(z)) * z * exp(-z);
270 // Expint Ei for large x:
272 return (boost::math::expint(x
) - x
) * x
* exp(-x
);
277 // erf_inv in range [0, 0.5]
281 y
= boost::math::tools::epsilon
<mp_type
>() / 64;
283 return boost::math::erf_inv(y
) / (y
);
287 // log1p over [-0.5,0.5]
290 y
= (y
== 0) ? 1e-100 : boost::math::sign(y
) * 1e-100;
291 return (boost::math::log1p(y
) - y
+ y
* y
/ 2) / (y
);
295 // cbrt over [0.5, 1]
296 return boost::math::cbrt(x
);
300 // trigamma over [x,y]
303 return boost::math::trigamma(x
) * (x
* x
);
307 // trigamma over [x, INF]
309 mp_type y
= (x
== 0) ? (std::numeric_limits
<double>::max
)() / 2 : mp_type(1/x
);
310 return boost::math::trigamma(y
) * y
;
315 // Don't need to go past x = 1/1000 = 1e-3 for double, or
316 // 1/15000 = 0.0006 for long double, start at 1/7.75=0.13
318 return sqrt(arg
) * exp(-arg
) * boost::math::cyl_bessel_i(0, arg
);
323 mp_type xx
= sqrt(x
) * 2;
324 return (boost::math::cyl_bessel_i(0, xx
) - 1) / x
;
329 mp_type xx
= sqrt(x
) * 2;
330 return (boost::math::cyl_bessel_i(1, xx
) * 2 / xx
- 1 - x
/ 2) / (x
* x
);
336 return boost::math::cyl_bessel_i(1, xx
) * sqrt(xx
) * exp(-xx
);
341 mp_type xx
= sqrt(x
);
342 return boost::math::cyl_bessel_k(0, xx
) + log(xx
) * boost::math::cyl_bessel_i(0, xx
);
348 return boost::math::cyl_bessel_k(0, xx
) * exp(xx
) * sqrt(xx
);
353 mp_type xx
= sqrt(x
);
354 return (boost::math::cyl_bessel_k(1, xx
) - log(xx
) * boost::math::cyl_bessel_i(1, xx
) - 1 / xx
) / xx
;
360 return boost::math::cyl_bessel_k(1, xx
) * sqrt(xx
) * exp(xx
);
364 return boost::math::lambert_w0(x
);
369 return boost::math::lambert_w0(x
) / x
;
373 static const mp_type e1
= exp(mp_type(-1));
374 return x
/ -boost::math::lambert_w0(-e1
+ x
);
379 return 1 / boost::math::lambert_w0(xx
);
384 return boost::math::lambert_w0(ex
) - x
;
391 const boost::math::tools::polynomial
<mp_type
>& n
,
392 const boost::math::tools::polynomial
<mp_type
>& d
,
393 const mp_type
& x_offset
,
394 const mp_type
& y_offset
,
400 // do nothing here...