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 #define L22
7 //#include "../tools/ntl_rr_lanczos.hpp"
8 //#include "../tools/ntl_rr_digamma.hpp"
9 #include <boost/math/bindings/rr.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
15 #include <cmath>
16
17
f(const boost::math::ntl::RR & x,int variant)18 boost::math::ntl::RR f(const boost::math::ntl::RR& x, int variant)
19 {
20 static const boost::math::ntl::RR tiny = boost::math::tools::min_value<boost::math::ntl::RR>() * 64;
21 switch(variant)
22 {
23 case 0:
24 {
25 boost::math::ntl::RR x_ = sqrt(x == 0 ? 1e-80 : x);
26 return boost::math::erf(x_) / x_;
27 }
28 case 1:
29 {
30 boost::math::ntl::RR x_ = 1 / x;
31 return boost::math::erfc(x_) * x_ / exp(-x_ * x_);
32 }
33 case 2:
34 {
35 return boost::math::erfc(x) * x / exp(-x * x);
36 }
37 case 3:
38 {
39 boost::math::ntl::RR y(x);
40 if(y == 0)
41 y += tiny;
42 return boost::math::lgamma(y+2) / y - 0.5;
43 }
44 case 4:
45 //
46 // lgamma in the range [2,3], use:
47 //
48 // lgamma(x) = (x-2) * (x + 1) * (c + R(x - 2))
49 //
50 // Works well at 80-bit long double precision, but doesn't
51 // stretch to 128-bit precision.
52 //
53 if(x == 0)
54 {
55 return boost::lexical_cast<boost::math::ntl::RR>("0.42278433509846713939348790991759756895784066406008") / 3;
56 }
57 return boost::math::lgamma(x+2) / (x * (x+3));
58 case 5:
59 {
60 //
61 // lgamma in the range [1,2], use:
62 //
63 // lgamma(x) = (x - 1) * (x - 2) * (c + R(x - 1))
64 //
65 // works well over [1, 1.5] but not near 2 :-(
66 //
67 boost::math::ntl::RR r1 = boost::lexical_cast<boost::math::ntl::RR>("0.57721566490153286060651209008240243104215933593992");
68 boost::math::ntl::RR r2 = boost::lexical_cast<boost::math::ntl::RR>("0.42278433509846713939348790991759756895784066406008");
69 if(x == 0)
70 {
71 return r1;
72 }
73 if(x == 1)
74 {
75 return r2;
76 }
77 return boost::math::lgamma(x+1) / (x * (x - 1));
78 }
79 case 6:
80 {
81 //
82 // lgamma in the range [1.5,2], use:
83 //
84 // lgamma(x) = (2 - x) * (1 - x) * (c + R(2 - x))
85 //
86 // works well over [1.5, 2] but not near 1 :-(
87 //
88 boost::math::ntl::RR r1 = boost::lexical_cast<boost::math::ntl::RR>("0.57721566490153286060651209008240243104215933593992");
89 boost::math::ntl::RR r2 = boost::lexical_cast<boost::math::ntl::RR>("0.42278433509846713939348790991759756895784066406008");
90 if(x == 0)
91 {
92 return r2;
93 }
94 if(x == 1)
95 {
96 return r1;
97 }
98 return boost::math::lgamma(2-x) / (x * (x - 1));
99 }
100 case 7:
101 {
102 //
103 // erf_inv in range [0, 0.5]
104 //
105 boost::math::ntl::RR y = x;
106 if(y == 0)
107 y = boost::math::tools::epsilon<boost::math::ntl::RR>() / 64;
108 return boost::math::erf_inv(y) / (y * (y+10));
109 }
110 case 8:
111 {
112 //
113 // erfc_inv in range [0.25, 0.5]
114 // Use an y-offset of 0.25, and range [0, 0.25]
115 // abs error, auto y-offset.
116 //
117 boost::math::ntl::RR y = x;
118 if(y == 0)
119 y = boost::lexical_cast<boost::math::ntl::RR>("1e-5000");
120 return sqrt(-2 * log(y)) / boost::math::erfc_inv(y);
121 }
122 case 9:
123 {
124 boost::math::ntl::RR x2 = x;
125 if(x2 == 0)
126 x2 = boost::lexical_cast<boost::math::ntl::RR>("1e-5000");
127 boost::math::ntl::RR y = exp(-x2*x2); // sqrt(-log(x2)) - 5;
128 return boost::math::erfc_inv(y) / x2;
129 }
130 case 10:
131 {
132 //
133 // Digamma over the interval [1,2], set x-offset to 1
134 // and optimise for absolute error over [0,1].
135 //
136 int current_precision = boost::math::ntl::RR::precision();
137 if(current_precision < 1000)
138 boost::math::ntl::RR::SetPrecision(1000);
139 //
140 // This value for the root of digamma is calculated using our
141 // differentiated lanczos approximation. It agrees with Cody
142 // to ~ 25 digits and to Morris to 35 digits. See:
143 // TOMS ALGORITHM 708 (Didonato and Morris).
144 // and Math. Comp. 27, 123-127 (1973) by Cody, Strecok and Thacher.
145 //
146 //boost::math::ntl::RR root = boost::lexical_cast<boost::math::ntl::RR>("1.4616321449683623412626595423257213234331845807102825466429633351908372838889871");
147 //
148 // Actually better to calculate the root on the fly, it appears to be more
149 // accurate: convergence is easier with the 1000-bit value, the approximation
150 // produced agrees with functions.mathworld.com values to 35 digits even quite
151 // near the root.
152 //
153 static boost::math::tools::eps_tolerance<boost::math::ntl::RR> tol(1000);
154 static boost::uintmax_t max_iter = 1000;
155 boost::math::ntl::RR (*pdg)(boost::math::ntl::RR) = &boost::math::digamma;
156 static const boost::math::ntl::RR root = boost::math::tools::bracket_and_solve_root(pdg, boost::math::ntl::RR(1.4), boost::math::ntl::RR(1.5), true, tol, max_iter).first;
157
158 boost::math::ntl::RR x2 = x;
159 double lim = 1e-65;
160 if(fabs(x2 - root) < lim)
161 {
162 //
163 // This is a problem area:
164 // x2-root suffers cancellation error, so does digamma.
165 // That gets compounded again when Remez calculates the error
166 // function. This cludge seems to stop the worst of the problems:
167 //
168 static const boost::math::ntl::RR a = boost::math::digamma(root - lim) / -lim;
169 static const boost::math::ntl::RR b = boost::math::digamma(root + lim) / lim;
170 boost::math::ntl::RR fract = (x2 - root + lim) / (2*lim);
171 boost::math::ntl::RR r = (1-fract) * a + fract * b;
172 std::cout << "In root area: " << r;
173 return r;
174 }
175 boost::math::ntl::RR result = boost::math::digamma(x2) / (x2 - root);
176 if(current_precision < 1000)
177 boost::math::ntl::RR::SetPrecision(current_precision);
178 return result;
179 }
180 case 11:
181 // expm1:
182 if(x == 0)
183 {
184 static boost::math::ntl::RR lim = 1e-80;
185 static boost::math::ntl::RR a = boost::math::expm1(-lim);
186 static boost::math::ntl::RR b = boost::math::expm1(lim);
187 static boost::math::ntl::RR l = (b-a) / (2 * lim);
188 return l;
189 }
190 return boost::math::expm1(x) / x;
191 case 12:
192 // demo, and test case:
193 return exp(x);
194 case 13:
195 // K(k):
196 {
197 return boost::math::ellint_1(x);
198 }
199 case 14:
200 // K(k)
201 {
202 return boost::math::ellint_1(1-x) / log(x);
203 }
204 case 15:
205 // E(k)
206 {
207 // x = 1-k^2
208 boost::math::ntl::RR z = 1 - x * log(x);
209 return boost::math::ellint_2(sqrt(1-x)) / z;
210 }
211 case 16:
212 // Bessel I0(x) over [0,16]:
213 {
214 return boost::math::cyl_bessel_i(0, sqrt(x));
215 }
216 case 17:
217 // Bessel I0(x) over [16,INF]
218 {
219 boost::math::ntl::RR z = 1 / (boost::math::ntl::RR(1)/16 - x);
220 return boost::math::cyl_bessel_i(0, z) * sqrt(z) / exp(z);
221 }
222 case 18:
223 // Zeta over [0, 1]
224 {
225 return boost::math::zeta(1 - x) * x - x;
226 }
227 case 19:
228 // Zeta over [1, n]
229 {
230 return boost::math::zeta(x) - 1 / (x - 1);
231 }
232 case 20:
233 // Zeta over [a, b] : a >> 1
234 {
235 return log(boost::math::zeta(x) - 1);
236 }
237 case 21:
238 // expint[1] over [0,1]:
239 {
240 boost::math::ntl::RR tiny = boost::lexical_cast<boost::math::ntl::RR>("1e-5000");
241 boost::math::ntl::RR z = (x <= tiny) ? tiny : x;
242 return boost::math::expint(1, z) - z + log(z);
243 }
244 case 22:
245 // expint[1] over [1,N],
246 // Note that x varies from [0,1]:
247 {
248 boost::math::ntl::RR z = 1 / x;
249 return boost::math::expint(1, z) * exp(z) * z;
250 }
251 case 23:
252 // expin Ei over [0,R]
253 {
254 static const boost::math::ntl::RR root =
255 boost::lexical_cast<boost::math::ntl::RR>("0.372507410781366634461991866580119133535689497771654051555657435242200120636201854384926049951548942392");
256 boost::math::ntl::RR z = x < (std::numeric_limits<long double>::min)() ? (std::numeric_limits<long double>::min)() : x;
257 return (boost::math::expint(z) - log(z / root)) / (z - root);
258 }
259 case 24:
260 // Expint Ei for large x:
261 {
262 static const boost::math::ntl::RR root =
263 boost::lexical_cast<boost::math::ntl::RR>("0.372507410781366634461991866580119133535689497771654051555657435242200120636201854384926049951548942392");
264 boost::math::ntl::RR z = x < (std::numeric_limits<long double>::min)() ? (std::numeric_limits<long double>::max)() : 1 / x;
265 return (boost::math::expint(z) - z) * z * exp(-z);
266 //return (boost::math::expint(z) - log(z)) * z * exp(-z);
267 }
268 case 25:
269 // Expint Ei for large x:
270 {
271 return (boost::math::expint(x) - x) * x * exp(-x);
272 }
273 case 26:
274 {
275 //
276 // erf_inv in range [0, 0.5]
277 //
278 boost::math::ntl::RR y = x;
279 if(y == 0)
280 y = boost::math::tools::epsilon<boost::math::ntl::RR>() / 64;
281 y = sqrt(y);
282 return boost::math::erf_inv(y) / (y);
283 }
284 case 28:
285 {
286 // log1p over [-0.5,0.5]
287 boost::math::ntl::RR y = x;
288 if(fabs(y) < 1e-100)
289 y = (y == 0) ? 1e-100 : boost::math::sign(y) * 1e-100;
290 return (boost::math::log1p(y) - y + y * y / 2) / (y);
291 }
292 case 29:
293 {
294 // cbrt over [0.5, 1]
295 return boost::math::cbrt(x);
296 }
297 case 30:
298 {
299 // trigamma over [x,y]
300 boost::math::ntl::RR y = x;
301 y = sqrt(y);
302 return boost::math::trigamma(x) * (x * x);
303 }
304 case 31:
305 {
306 // trigamma over [x, INF]
307 if(x == 0) return 1;
308 boost::math::ntl::RR y = (x == 0) ? (std::numeric_limits<double>::max)() / 2 : 1/x;
309 return boost::math::trigamma(y) * y;
310 }
311 }
312 return 0;
313 }
314
show_extra(const boost::math::tools::polynomial<boost::math::ntl::RR> & n,const boost::math::tools::polynomial<boost::math::ntl::RR> & d,const boost::math::ntl::RR & x_offset,const boost::math::ntl::RR & y_offset,int variant)315 void show_extra(
316 const boost::math::tools::polynomial<boost::math::ntl::RR>& n,
317 const boost::math::tools::polynomial<boost::math::ntl::RR>& d,
318 const boost::math::ntl::RR& x_offset,
319 const boost::math::ntl::RR& y_offset,
320 int variant)
321 {
322 switch(variant)
323 {
324 default:
325 // do nothing here...
326 ;
327 }
328 }
329
330