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