]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Copyright John Maddock 2007. |
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_NTL_RR_HPP | |
7 | #define BOOST_MATH_NTL_RR_HPP | |
8 | ||
9 | #include <boost/config.hpp> | |
10 | #include <boost/limits.hpp> | |
11 | #include <boost/math/tools/real_cast.hpp> | |
12 | #include <boost/math/tools/precision.hpp> | |
13 | #include <boost/math/constants/constants.hpp> | |
14 | #include <boost/math/tools/roots.hpp> | |
15 | #include <boost/math/special_functions/fpclassify.hpp> | |
16 | #include <boost/math/bindings/detail/big_digamma.hpp> | |
17 | #include <boost/math/bindings/detail/big_lanczos.hpp> | |
18 | ||
19 | #include <ostream> | |
20 | #include <istream> | |
21 | #include <boost/config/no_tr1/cmath.hpp> | |
22 | #include <NTL/RR.h> | |
23 | ||
24 | namespace boost{ namespace math{ | |
25 | ||
26 | namespace ntl | |
27 | { | |
28 | ||
29 | class RR; | |
30 | ||
31 | RR ldexp(RR r, int exp); | |
32 | RR frexp(RR r, int* exp); | |
33 | ||
34 | class RR | |
35 | { | |
36 | public: | |
37 | // Constructors: | |
38 | RR() {} | |
39 | RR(const ::NTL::RR& c) : m_value(c){} | |
40 | RR(char c) | |
41 | { | |
42 | m_value = c; | |
43 | } | |
44 | #ifndef BOOST_NO_INTRINSIC_WCHAR_T | |
45 | RR(wchar_t c) | |
46 | { | |
47 | m_value = c; | |
48 | } | |
49 | #endif | |
50 | RR(unsigned char c) | |
51 | { | |
52 | m_value = c; | |
53 | } | |
54 | RR(signed char c) | |
55 | { | |
56 | m_value = c; | |
57 | } | |
58 | RR(unsigned short c) | |
59 | { | |
60 | m_value = c; | |
61 | } | |
62 | RR(short c) | |
63 | { | |
64 | m_value = c; | |
65 | } | |
66 | RR(unsigned int c) | |
67 | { | |
68 | assign_large_int(c); | |
69 | } | |
70 | RR(int c) | |
71 | { | |
72 | assign_large_int(c); | |
73 | } | |
74 | RR(unsigned long c) | |
75 | { | |
76 | assign_large_int(c); | |
77 | } | |
78 | RR(long c) | |
79 | { | |
80 | assign_large_int(c); | |
81 | } | |
82 | #ifdef BOOST_HAS_LONG_LONG | |
83 | RR(boost::ulong_long_type c) | |
84 | { | |
85 | assign_large_int(c); | |
86 | } | |
87 | RR(boost::long_long_type c) | |
88 | { | |
89 | assign_large_int(c); | |
90 | } | |
91 | #endif | |
92 | RR(float c) | |
93 | { | |
94 | m_value = c; | |
95 | } | |
96 | RR(double c) | |
97 | { | |
98 | m_value = c; | |
99 | } | |
100 | RR(long double c) | |
101 | { | |
102 | assign_large_real(c); | |
103 | } | |
104 | ||
105 | // Assignment: | |
106 | RR& operator=(char c) { m_value = c; return *this; } | |
107 | RR& operator=(unsigned char c) { m_value = c; return *this; } | |
108 | RR& operator=(signed char c) { m_value = c; return *this; } | |
109 | #ifndef BOOST_NO_INTRINSIC_WCHAR_T | |
110 | RR& operator=(wchar_t c) { m_value = c; return *this; } | |
111 | #endif | |
112 | RR& operator=(short c) { m_value = c; return *this; } | |
113 | RR& operator=(unsigned short c) { m_value = c; return *this; } | |
114 | RR& operator=(int c) { assign_large_int(c); return *this; } | |
115 | RR& operator=(unsigned int c) { assign_large_int(c); return *this; } | |
116 | RR& operator=(long c) { assign_large_int(c); return *this; } | |
117 | RR& operator=(unsigned long c) { assign_large_int(c); return *this; } | |
118 | #ifdef BOOST_HAS_LONG_LONG | |
119 | RR& operator=(boost::long_long_type c) { assign_large_int(c); return *this; } | |
120 | RR& operator=(boost::ulong_long_type c) { assign_large_int(c); return *this; } | |
121 | #endif | |
122 | RR& operator=(float c) { m_value = c; return *this; } | |
123 | RR& operator=(double c) { m_value = c; return *this; } | |
124 | RR& operator=(long double c) { assign_large_real(c); return *this; } | |
125 | ||
126 | // Access: | |
127 | NTL::RR& value(){ return m_value; } | |
128 | NTL::RR const& value()const{ return m_value; } | |
129 | ||
130 | // Member arithmetic: | |
131 | RR& operator+=(const RR& other) | |
132 | { m_value += other.value(); return *this; } | |
133 | RR& operator-=(const RR& other) | |
134 | { m_value -= other.value(); return *this; } | |
135 | RR& operator*=(const RR& other) | |
136 | { m_value *= other.value(); return *this; } | |
137 | RR& operator/=(const RR& other) | |
138 | { m_value /= other.value(); return *this; } | |
139 | RR operator-()const | |
140 | { return -m_value; } | |
141 | RR const& operator+()const | |
142 | { return *this; } | |
143 | ||
144 | // RR compatibity: | |
145 | const ::NTL::ZZ& mantissa() const | |
146 | { return m_value.mantissa(); } | |
147 | long exponent() const | |
148 | { return m_value.exponent(); } | |
149 | ||
150 | static void SetPrecision(long p) | |
151 | { ::NTL::RR::SetPrecision(p); } | |
152 | ||
153 | static long precision() | |
154 | { return ::NTL::RR::precision(); } | |
155 | ||
156 | static void SetOutputPrecision(long p) | |
157 | { ::NTL::RR::SetOutputPrecision(p); } | |
158 | static long OutputPrecision() | |
159 | { return ::NTL::RR::OutputPrecision(); } | |
160 | ||
161 | ||
162 | private: | |
163 | ::NTL::RR m_value; | |
164 | ||
165 | template <class V> | |
166 | void assign_large_real(const V& a) | |
167 | { | |
168 | using std::frexp; | |
169 | using std::ldexp; | |
170 | using std::floor; | |
171 | if (a == 0) { | |
172 | clear(m_value); | |
173 | return; | |
174 | } | |
175 | ||
176 | if (a == 1) { | |
177 | NTL::set(m_value); | |
178 | return; | |
179 | } | |
180 | ||
181 | if (!(boost::math::isfinite)(a)) | |
182 | { | |
183 | throw std::overflow_error("Cannot construct an instance of NTL::RR with an infinite value."); | |
184 | } | |
185 | ||
186 | int e; | |
187 | long double f, term; | |
188 | ::NTL::RR t; | |
189 | clear(m_value); | |
190 | ||
191 | f = frexp(a, &e); | |
192 | ||
193 | while(f) | |
194 | { | |
195 | // extract 30 bits from f: | |
196 | f = ldexp(f, 30); | |
197 | term = floor(f); | |
198 | e -= 30; | |
199 | conv(t.x, (int)term); | |
200 | t.e = e; | |
201 | m_value += t; | |
202 | f -= term; | |
203 | } | |
204 | } | |
205 | ||
206 | template <class V> | |
207 | void assign_large_int(V a) | |
208 | { | |
209 | #ifdef BOOST_MSVC | |
210 | #pragma warning(push) | |
211 | #pragma warning(disable:4146) | |
212 | #endif | |
213 | clear(m_value); | |
214 | int exp = 0; | |
215 | NTL::RR t; | |
216 | bool neg = a < V(0) ? true : false; | |
217 | if(neg) | |
218 | a = -a; | |
219 | while(a) | |
220 | { | |
221 | t = static_cast<double>(a & 0xffff); | |
222 | m_value += ldexp(RR(t), exp).value(); | |
223 | a >>= 16; | |
224 | exp += 16; | |
225 | } | |
226 | if(neg) | |
227 | m_value = -m_value; | |
228 | #ifdef BOOST_MSVC | |
229 | #pragma warning(pop) | |
230 | #endif | |
231 | } | |
232 | }; | |
233 | ||
234 | // Non-member arithmetic: | |
235 | inline RR operator+(const RR& a, const RR& b) | |
236 | { | |
237 | RR result(a); | |
238 | result += b; | |
239 | return result; | |
240 | } | |
241 | inline RR operator-(const RR& a, const RR& b) | |
242 | { | |
243 | RR result(a); | |
244 | result -= b; | |
245 | return result; | |
246 | } | |
247 | inline RR operator*(const RR& a, const RR& b) | |
248 | { | |
249 | RR result(a); | |
250 | result *= b; | |
251 | return result; | |
252 | } | |
253 | inline RR operator/(const RR& a, const RR& b) | |
254 | { | |
255 | RR result(a); | |
256 | result /= b; | |
257 | return result; | |
258 | } | |
259 | ||
260 | // Comparison: | |
261 | inline bool operator == (const RR& a, const RR& b) | |
262 | { return a.value() == b.value() ? true : false; } | |
263 | inline bool operator != (const RR& a, const RR& b) | |
264 | { return a.value() != b.value() ? true : false;} | |
265 | inline bool operator < (const RR& a, const RR& b) | |
266 | { return a.value() < b.value() ? true : false; } | |
267 | inline bool operator <= (const RR& a, const RR& b) | |
268 | { return a.value() <= b.value() ? true : false; } | |
269 | inline bool operator > (const RR& a, const RR& b) | |
270 | { return a.value() > b.value() ? true : false; } | |
271 | inline bool operator >= (const RR& a, const RR& b) | |
272 | { return a.value() >= b.value() ? true : false; } | |
273 | ||
274 | #if 0 | |
275 | // Non-member mixed compare: | |
276 | template <class T> | |
277 | inline bool operator == (const T& a, const RR& b) | |
278 | { | |
279 | return a == b.value(); | |
280 | } | |
281 | template <class T> | |
282 | inline bool operator != (const T& a, const RR& b) | |
283 | { | |
284 | return a != b.value(); | |
285 | } | |
286 | template <class T> | |
287 | inline bool operator < (const T& a, const RR& b) | |
288 | { | |
289 | return a < b.value(); | |
290 | } | |
291 | template <class T> | |
292 | inline bool operator > (const T& a, const RR& b) | |
293 | { | |
294 | return a > b.value(); | |
295 | } | |
296 | template <class T> | |
297 | inline bool operator <= (const T& a, const RR& b) | |
298 | { | |
299 | return a <= b.value(); | |
300 | } | |
301 | template <class T> | |
302 | inline bool operator >= (const T& a, const RR& b) | |
303 | { | |
304 | return a >= b.value(); | |
305 | } | |
306 | #endif // Non-member mixed compare: | |
307 | ||
308 | // Non-member functions: | |
309 | /* | |
310 | inline RR acos(RR a) | |
311 | { return ::NTL::acos(a.value()); } | |
312 | */ | |
313 | inline RR cos(RR a) | |
314 | { return ::NTL::cos(a.value()); } | |
315 | /* | |
316 | inline RR asin(RR a) | |
317 | { return ::NTL::asin(a.value()); } | |
318 | inline RR atan(RR a) | |
319 | { return ::NTL::atan(a.value()); } | |
320 | inline RR atan2(RR a, RR b) | |
321 | { return ::NTL::atan2(a.value(), b.value()); } | |
322 | */ | |
323 | inline RR ceil(RR a) | |
324 | { return ::NTL::ceil(a.value()); } | |
325 | /* | |
326 | inline RR fmod(RR a, RR b) | |
327 | { return ::NTL::fmod(a.value(), b.value()); } | |
328 | inline RR cosh(RR a) | |
329 | { return ::NTL::cosh(a.value()); } | |
330 | */ | |
331 | inline RR exp(RR a) | |
332 | { return ::NTL::exp(a.value()); } | |
333 | inline RR fabs(RR a) | |
334 | { return ::NTL::fabs(a.value()); } | |
335 | inline RR abs(RR a) | |
336 | { return ::NTL::abs(a.value()); } | |
337 | inline RR floor(RR a) | |
338 | { return ::NTL::floor(a.value()); } | |
339 | /* | |
340 | inline RR modf(RR a, RR* ipart) | |
341 | { | |
342 | ::NTL::RR ip; | |
343 | RR result = modf(a.value(), &ip); | |
344 | *ipart = ip; | |
345 | return result; | |
346 | } | |
347 | inline RR frexp(RR a, int* expon) | |
348 | { return ::NTL::frexp(a.value(), expon); } | |
349 | inline RR ldexp(RR a, int expon) | |
350 | { return ::NTL::ldexp(a.value(), expon); } | |
351 | */ | |
352 | inline RR log(RR a) | |
353 | { return ::NTL::log(a.value()); } | |
354 | inline RR log10(RR a) | |
355 | { return ::NTL::log10(a.value()); } | |
356 | /* | |
357 | inline RR tan(RR a) | |
358 | { return ::NTL::tan(a.value()); } | |
359 | */ | |
360 | inline RR pow(RR a, RR b) | |
361 | { return ::NTL::pow(a.value(), b.value()); } | |
362 | inline RR pow(RR a, int b) | |
363 | { return ::NTL::power(a.value(), b); } | |
364 | inline RR sin(RR a) | |
365 | { return ::NTL::sin(a.value()); } | |
366 | /* | |
367 | inline RR sinh(RR a) | |
368 | { return ::NTL::sinh(a.value()); } | |
369 | */ | |
370 | inline RR sqrt(RR a) | |
371 | { return ::NTL::sqrt(a.value()); } | |
372 | /* | |
373 | inline RR tanh(RR a) | |
374 | { return ::NTL::tanh(a.value()); } | |
375 | */ | |
376 | inline RR pow(const RR& r, long l) | |
377 | { | |
378 | return ::NTL::power(r.value(), l); | |
379 | } | |
380 | inline RR tan(const RR& a) | |
381 | { | |
382 | return sin(a)/cos(a); | |
383 | } | |
384 | inline RR frexp(RR r, int* exp) | |
385 | { | |
386 | *exp = r.value().e; | |
387 | r.value().e = 0; | |
388 | while(r >= 1) | |
389 | { | |
390 | *exp += 1; | |
391 | r.value().e -= 1; | |
392 | } | |
393 | while(r < 0.5) | |
394 | { | |
395 | *exp -= 1; | |
396 | r.value().e += 1; | |
397 | } | |
398 | BOOST_ASSERT(r < 1); | |
399 | BOOST_ASSERT(r >= 0.5); | |
400 | return r; | |
401 | } | |
402 | inline RR ldexp(RR r, int exp) | |
403 | { | |
404 | r.value().e += exp; | |
405 | return r; | |
406 | } | |
407 | ||
408 | // Streaming: | |
409 | template <class charT, class traits> | |
410 | inline std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& os, const RR& a) | |
411 | { | |
412 | return os << a.value(); | |
413 | } | |
414 | template <class charT, class traits> | |
415 | inline std::basic_istream<charT, traits>& operator>>(std::basic_istream<charT, traits>& is, RR& a) | |
416 | { | |
417 | ::NTL::RR v; | |
418 | is >> v; | |
419 | a = v; | |
420 | return is; | |
421 | } | |
422 | ||
423 | } // namespace ntl | |
424 | ||
425 | namespace lanczos{ | |
426 | ||
427 | struct ntl_lanczos | |
428 | { | |
429 | static ntl::RR lanczos_sum(const ntl::RR& z) | |
430 | { | |
431 | unsigned long p = ntl::RR::precision(); | |
432 | if(p <= 72) | |
433 | return lanczos13UDT::lanczos_sum(z); | |
434 | else if(p <= 120) | |
435 | return lanczos22UDT::lanczos_sum(z); | |
436 | else if(p <= 170) | |
437 | return lanczos31UDT::lanczos_sum(z); | |
438 | else //if(p <= 370) approx 100 digit precision: | |
439 | return lanczos61UDT::lanczos_sum(z); | |
440 | } | |
441 | static ntl::RR lanczos_sum_expG_scaled(const ntl::RR& z) | |
442 | { | |
443 | unsigned long p = ntl::RR::precision(); | |
444 | if(p <= 72) | |
445 | return lanczos13UDT::lanczos_sum_expG_scaled(z); | |
446 | else if(p <= 120) | |
447 | return lanczos22UDT::lanczos_sum_expG_scaled(z); | |
448 | else if(p <= 170) | |
449 | return lanczos31UDT::lanczos_sum_expG_scaled(z); | |
450 | else //if(p <= 370) approx 100 digit precision: | |
451 | return lanczos61UDT::lanczos_sum_expG_scaled(z); | |
452 | } | |
453 | static ntl::RR lanczos_sum_near_1(const ntl::RR& z) | |
454 | { | |
455 | unsigned long p = ntl::RR::precision(); | |
456 | if(p <= 72) | |
457 | return lanczos13UDT::lanczos_sum_near_1(z); | |
458 | else if(p <= 120) | |
459 | return lanczos22UDT::lanczos_sum_near_1(z); | |
460 | else if(p <= 170) | |
461 | return lanczos31UDT::lanczos_sum_near_1(z); | |
462 | else //if(p <= 370) approx 100 digit precision: | |
463 | return lanczos61UDT::lanczos_sum_near_1(z); | |
464 | } | |
465 | static ntl::RR lanczos_sum_near_2(const ntl::RR& z) | |
466 | { | |
467 | unsigned long p = ntl::RR::precision(); | |
468 | if(p <= 72) | |
469 | return lanczos13UDT::lanczos_sum_near_2(z); | |
470 | else if(p <= 120) | |
471 | return lanczos22UDT::lanczos_sum_near_2(z); | |
472 | else if(p <= 170) | |
473 | return lanczos31UDT::lanczos_sum_near_2(z); | |
474 | else //if(p <= 370) approx 100 digit precision: | |
475 | return lanczos61UDT::lanczos_sum_near_2(z); | |
476 | } | |
477 | static ntl::RR g() | |
478 | { | |
479 | unsigned long p = ntl::RR::precision(); | |
480 | if(p <= 72) | |
481 | return lanczos13UDT::g(); | |
482 | else if(p <= 120) | |
483 | return lanczos22UDT::g(); | |
484 | else if(p <= 170) | |
485 | return lanczos31UDT::g(); | |
486 | else //if(p <= 370) approx 100 digit precision: | |
487 | return lanczos61UDT::g(); | |
488 | } | |
489 | }; | |
490 | ||
491 | template<class Policy> | |
492 | struct lanczos<ntl::RR, Policy> | |
493 | { | |
494 | typedef ntl_lanczos type; | |
495 | }; | |
496 | ||
497 | } // namespace lanczos | |
498 | ||
499 | namespace tools | |
500 | { | |
501 | ||
502 | template<> | |
503 | inline int digits<boost::math::ntl::RR>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(boost::math::ntl::RR)) | |
504 | { | |
505 | return ::NTL::RR::precision(); | |
506 | } | |
507 | ||
508 | template <> | |
509 | inline float real_cast<float, boost::math::ntl::RR>(boost::math::ntl::RR t) | |
510 | { | |
511 | double r; | |
512 | conv(r, t.value()); | |
513 | return static_cast<float>(r); | |
514 | } | |
515 | template <> | |
516 | inline double real_cast<double, boost::math::ntl::RR>(boost::math::ntl::RR t) | |
517 | { | |
518 | double r; | |
519 | conv(r, t.value()); | |
520 | return r; | |
521 | } | |
522 | ||
523 | namespace detail{ | |
524 | ||
525 | template<class I> | |
526 | void convert_to_long_result(NTL::RR const& r, I& result) | |
527 | { | |
528 | result = 0; | |
529 | I last_result(0); | |
530 | NTL::RR t(r); | |
531 | double term; | |
532 | do | |
533 | { | |
534 | conv(term, t); | |
535 | last_result = result; | |
536 | result += static_cast<I>(term); | |
537 | t -= term; | |
538 | }while(result != last_result); | |
539 | } | |
540 | ||
541 | } | |
542 | ||
543 | template <> | |
544 | inline long double real_cast<long double, boost::math::ntl::RR>(boost::math::ntl::RR t) | |
545 | { | |
546 | long double result(0); | |
547 | detail::convert_to_long_result(t.value(), result); | |
548 | return result; | |
549 | } | |
550 | template <> | |
551 | inline boost::math::ntl::RR real_cast<boost::math::ntl::RR, boost::math::ntl::RR>(boost::math::ntl::RR t) | |
552 | { | |
553 | return t; | |
554 | } | |
555 | template <> | |
556 | inline unsigned real_cast<unsigned, boost::math::ntl::RR>(boost::math::ntl::RR t) | |
557 | { | |
558 | unsigned result; | |
559 | detail::convert_to_long_result(t.value(), result); | |
560 | return result; | |
561 | } | |
562 | template <> | |
563 | inline int real_cast<int, boost::math::ntl::RR>(boost::math::ntl::RR t) | |
564 | { | |
565 | int result; | |
566 | detail::convert_to_long_result(t.value(), result); | |
567 | return result; | |
568 | } | |
569 | template <> | |
570 | inline long real_cast<long, boost::math::ntl::RR>(boost::math::ntl::RR t) | |
571 | { | |
572 | long result; | |
573 | detail::convert_to_long_result(t.value(), result); | |
574 | return result; | |
575 | } | |
576 | template <> | |
577 | inline long long real_cast<long long, boost::math::ntl::RR>(boost::math::ntl::RR t) | |
578 | { | |
579 | long long result; | |
580 | detail::convert_to_long_result(t.value(), result); | |
581 | return result; | |
582 | } | |
583 | ||
584 | template <> | |
585 | inline boost::math::ntl::RR max_value<boost::math::ntl::RR>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(boost::math::ntl::RR)) | |
586 | { | |
587 | static bool has_init = false; | |
588 | static NTL::RR val; | |
589 | if(!has_init) | |
590 | { | |
591 | val = 1; | |
592 | val.e = NTL_OVFBND-20; | |
593 | has_init = true; | |
594 | } | |
595 | return val; | |
596 | } | |
597 | ||
598 | template <> | |
599 | inline boost::math::ntl::RR min_value<boost::math::ntl::RR>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(boost::math::ntl::RR)) | |
600 | { | |
601 | static bool has_init = false; | |
602 | static NTL::RR val; | |
603 | if(!has_init) | |
604 | { | |
605 | val = 1; | |
606 | val.e = -NTL_OVFBND+20; | |
607 | has_init = true; | |
608 | } | |
609 | return val; | |
610 | } | |
611 | ||
612 | template <> | |
613 | inline boost::math::ntl::RR log_max_value<boost::math::ntl::RR>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(boost::math::ntl::RR)) | |
614 | { | |
615 | static bool has_init = false; | |
616 | static NTL::RR val; | |
617 | if(!has_init) | |
618 | { | |
619 | val = 1; | |
620 | val.e = NTL_OVFBND-20; | |
621 | val = log(val); | |
622 | has_init = true; | |
623 | } | |
624 | return val; | |
625 | } | |
626 | ||
627 | template <> | |
628 | inline boost::math::ntl::RR log_min_value<boost::math::ntl::RR>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(boost::math::ntl::RR)) | |
629 | { | |
630 | static bool has_init = false; | |
631 | static NTL::RR val; | |
632 | if(!has_init) | |
633 | { | |
634 | val = 1; | |
635 | val.e = -NTL_OVFBND+20; | |
636 | val = log(val); | |
637 | has_init = true; | |
638 | } | |
639 | return val; | |
640 | } | |
641 | ||
642 | template <> | |
643 | inline boost::math::ntl::RR epsilon<boost::math::ntl::RR>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(boost::math::ntl::RR)) | |
644 | { | |
645 | return ldexp(boost::math::ntl::RR(1), 1-boost::math::policies::digits<boost::math::ntl::RR, boost::math::policies::policy<> >()); | |
646 | } | |
647 | ||
648 | } // namespace tools | |
649 | ||
650 | // | |
651 | // The number of digits precision in RR can vary with each call | |
652 | // so we need to recalculate these with each call: | |
653 | // | |
654 | namespace constants{ | |
655 | ||
656 | template<> inline boost::math::ntl::RR pi<boost::math::ntl::RR>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(boost::math::ntl::RR)) | |
657 | { | |
658 | NTL::RR result; | |
659 | ComputePi(result); | |
660 | return result; | |
661 | } | |
662 | template<> inline boost::math::ntl::RR e<boost::math::ntl::RR>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(boost::math::ntl::RR)) | |
663 | { | |
664 | NTL::RR result; | |
665 | result = 1; | |
666 | return exp(result); | |
667 | } | |
668 | ||
669 | } // namespace constants | |
670 | ||
671 | namespace ntl{ | |
672 | // | |
673 | // These are some fairly brain-dead versions of the math | |
674 | // functions that NTL fails to provide. | |
675 | // | |
676 | ||
677 | ||
678 | // | |
679 | // Inverse trig functions: | |
680 | // | |
681 | struct asin_root | |
682 | { | |
683 | asin_root(RR const& target) : t(target){} | |
684 | ||
685 | boost::math::tuple<RR, RR, RR> operator()(RR const& p) | |
686 | { | |
687 | RR f0 = sin(p); | |
688 | RR f1 = cos(p); | |
689 | RR f2 = -f0; | |
690 | f0 -= t; | |
691 | return boost::math::make_tuple(f0, f1, f2); | |
692 | } | |
693 | private: | |
694 | RR t; | |
695 | }; | |
696 | ||
697 | inline RR asin(RR z) | |
698 | { | |
699 | double r; | |
700 | conv(r, z.value()); | |
701 | return boost::math::tools::halley_iterate( | |
702 | asin_root(z), | |
703 | RR(std::asin(r)), | |
704 | RR(-boost::math::constants::pi<RR>()/2), | |
705 | RR(boost::math::constants::pi<RR>()/2), | |
706 | NTL::RR::precision()); | |
707 | } | |
708 | ||
709 | struct acos_root | |
710 | { | |
711 | acos_root(RR const& target) : t(target){} | |
712 | ||
713 | boost::math::tuple<RR, RR, RR> operator()(RR const& p) | |
714 | { | |
715 | RR f0 = cos(p); | |
716 | RR f1 = -sin(p); | |
717 | RR f2 = -f0; | |
718 | f0 -= t; | |
719 | return boost::math::make_tuple(f0, f1, f2); | |
720 | } | |
721 | private: | |
722 | RR t; | |
723 | }; | |
724 | ||
725 | inline RR acos(RR z) | |
726 | { | |
727 | double r; | |
728 | conv(r, z.value()); | |
729 | return boost::math::tools::halley_iterate( | |
730 | acos_root(z), | |
731 | RR(std::acos(r)), | |
732 | RR(-boost::math::constants::pi<RR>()/2), | |
733 | RR(boost::math::constants::pi<RR>()/2), | |
734 | NTL::RR::precision()); | |
735 | } | |
736 | ||
737 | struct atan_root | |
738 | { | |
739 | atan_root(RR const& target) : t(target){} | |
740 | ||
741 | boost::math::tuple<RR, RR, RR> operator()(RR const& p) | |
742 | { | |
743 | RR c = cos(p); | |
744 | RR ta = tan(p); | |
745 | RR f0 = ta - t; | |
746 | RR f1 = 1 / (c * c); | |
747 | RR f2 = 2 * ta / (c * c); | |
748 | return boost::math::make_tuple(f0, f1, f2); | |
749 | } | |
750 | private: | |
751 | RR t; | |
752 | }; | |
753 | ||
754 | inline RR atan(RR z) | |
755 | { | |
756 | double r; | |
757 | conv(r, z.value()); | |
758 | return boost::math::tools::halley_iterate( | |
759 | atan_root(z), | |
760 | RR(std::atan(r)), | |
761 | -boost::math::constants::pi<RR>()/2, | |
762 | boost::math::constants::pi<RR>()/2, | |
763 | NTL::RR::precision()); | |
764 | } | |
765 | ||
766 | inline RR atan2(RR y, RR x) | |
767 | { | |
768 | if(x > 0) | |
769 | return atan(y / x); | |
770 | if(x < 0) | |
771 | { | |
772 | return y < 0 ? atan(y / x) - boost::math::constants::pi<RR>() : atan(y / x) + boost::math::constants::pi<RR>(); | |
773 | } | |
774 | return y < 0 ? -boost::math::constants::half_pi<RR>() : boost::math::constants::half_pi<RR>() ; | |
775 | } | |
776 | ||
777 | inline RR sinh(RR z) | |
778 | { | |
779 | return (expm1(z.value()) - expm1(-z.value())) / 2; | |
780 | } | |
781 | ||
782 | inline RR cosh(RR z) | |
783 | { | |
784 | return (exp(z) + exp(-z)) / 2; | |
785 | } | |
786 | ||
787 | inline RR tanh(RR z) | |
788 | { | |
789 | return sinh(z) / cosh(z); | |
790 | } | |
791 | ||
792 | inline RR fmod(RR x, RR y) | |
793 | { | |
794 | // This is a really crummy version of fmod, we rely on lots | |
795 | // of digits to get us out of trouble... | |
796 | RR factor = floor(x/y); | |
797 | return x - factor * y; | |
798 | } | |
799 | ||
800 | template <class Policy> | |
801 | inline int iround(RR const& x, const Policy& pol) | |
802 | { | |
803 | return tools::real_cast<int>(round(x, pol)); | |
804 | } | |
805 | ||
806 | template <class Policy> | |
807 | inline long lround(RR const& x, const Policy& pol) | |
808 | { | |
809 | return tools::real_cast<long>(round(x, pol)); | |
810 | } | |
811 | ||
812 | template <class Policy> | |
813 | inline long long llround(RR const& x, const Policy& pol) | |
814 | { | |
815 | return tools::real_cast<long long>(round(x, pol)); | |
816 | } | |
817 | ||
818 | template <class Policy> | |
819 | inline int itrunc(RR const& x, const Policy& pol) | |
820 | { | |
821 | return tools::real_cast<int>(trunc(x, pol)); | |
822 | } | |
823 | ||
824 | template <class Policy> | |
825 | inline long ltrunc(RR const& x, const Policy& pol) | |
826 | { | |
827 | return tools::real_cast<long>(trunc(x, pol)); | |
828 | } | |
829 | ||
830 | template <class Policy> | |
831 | inline long long lltrunc(RR const& x, const Policy& pol) | |
832 | { | |
833 | return tools::real_cast<long long>(trunc(x, pol)); | |
834 | } | |
835 | ||
836 | } // namespace ntl | |
837 | ||
838 | namespace detail{ | |
839 | ||
840 | template <class Policy> | |
841 | ntl::RR digamma_imp(ntl::RR x, const mpl::int_<0>* , const Policy& pol) | |
842 | { | |
843 | // | |
844 | // This handles reflection of negative arguments, and all our | |
845 | // error handling, then forwards to the T-specific approximation. | |
846 | // | |
847 | BOOST_MATH_STD_USING // ADL of std functions. | |
848 | ||
849 | ntl::RR result = 0; | |
850 | // | |
851 | // Check for negative arguments and use reflection: | |
852 | // | |
853 | if(x < 0) | |
854 | { | |
855 | // Reflect: | |
856 | x = 1 - x; | |
857 | // Argument reduction for tan: | |
858 | ntl::RR remainder = x - floor(x); | |
859 | // Shift to negative if > 0.5: | |
860 | if(remainder > 0.5) | |
861 | { | |
862 | remainder -= 1; | |
863 | } | |
864 | // | |
865 | // check for evaluation at a negative pole: | |
866 | // | |
867 | if(remainder == 0) | |
868 | { | |
869 | return policies::raise_pole_error<ntl::RR>("boost::math::digamma<%1%>(%1%)", 0, (1-x), pol); | |
870 | } | |
871 | result = constants::pi<ntl::RR>() / tan(constants::pi<ntl::RR>() * remainder); | |
872 | } | |
873 | result += big_digamma(x); | |
874 | return result; | |
875 | } | |
876 | ||
877 | } // namespace detail | |
878 | ||
879 | } // namespace math | |
880 | } // namespace boost | |
881 | ||
882 | #endif // BOOST_MATH_REAL_CONCEPT_HPP | |
883 | ||
884 |