1 // Copyright John Maddock 2008.
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)
6 // Wrapper that works with mpfr_class defined in gmpfrxx.h
7 // See http://math.berkeley.edu/~wilken/code/gmpfrxx/
8 // Also requires the gmp and mpfr libraries.
11 #ifndef BOOST_MATH_E_FLOAT_BINDINGS_HPP
12 #define BOOST_MATH_E_FLOAT_BINDINGS_HPP
14 #include <boost/config.hpp>
17 #include <e_float/e_float.h>
18 #include <functions/functions.h>
20 #include <boost/math/tools/precision.hpp>
21 #include <boost/math/tools/real_cast.hpp>
22 #include <boost/math/policies/policy.hpp>
23 #include <boost/math/distributions/fwd.hpp>
24 #include <boost/math/special_functions/math_fwd.hpp>
25 #include <boost/math/special_functions/fpclassify.hpp>
26 #include <boost/math/bindings/detail/big_digamma.hpp>
27 #include <boost/math/bindings/detail/big_lanczos.hpp>
28 #include <boost/lexical_cast.hpp>
31 namespace boost{ namespace math{ namespace ef{
38 e_float(const ::e_float& c) : m_value(c){}
41 m_value = ::e_float(c);
43 #ifndef BOOST_NO_INTRINSIC_WCHAR_T
46 m_value = ::e_float(c);
49 e_float(unsigned char c)
51 m_value = ::e_float(c);
53 e_float(signed char c)
55 m_value = ::e_float(c);
57 e_float(unsigned short c)
59 m_value = ::e_float(c);
63 m_value = ::e_float(c);
65 e_float(unsigned int c)
67 m_value = ::e_float(c);
71 m_value = ::e_float(c);
73 e_float(unsigned long c)
75 m_value = ::e_float((UINT64)c);
79 m_value = ::e_float((INT64)c);
81 #ifdef BOOST_HAS_LONG_LONG
82 e_float(boost::ulong_long_type c)
84 m_value = ::e_float(c);
86 e_float(boost::long_long_type c)
88 m_value = ::e_float(c);
99 e_float(long double c)
101 assign_large_real(c);
105 e_float& operator=(char c) { m_value = ::e_float(c); return *this; }
106 e_float& operator=(unsigned char c) { m_value = ::e_float(c); return *this; }
107 e_float& operator=(signed char c) { m_value = ::e_float(c); return *this; }
108 #ifndef BOOST_NO_INTRINSIC_WCHAR_T
109 e_float& operator=(wchar_t c) { m_value = ::e_float(c); return *this; }
111 e_float& operator=(short c) { m_value = ::e_float(c); return *this; }
112 e_float& operator=(unsigned short c) { m_value = ::e_float(c); return *this; }
113 e_float& operator=(int c) { m_value = ::e_float(c); return *this; }
114 e_float& operator=(unsigned int c) { m_value = ::e_float(c); return *this; }
115 e_float& operator=(long c) { m_value = ::e_float((INT64)c); return *this; }
116 e_float& operator=(unsigned long c) { m_value = ::e_float((UINT64)c); return *this; }
117 #ifdef BOOST_HAS_LONG_LONG
118 e_float& operator=(boost::long_long_type c) { m_value = ::e_float(c); return *this; }
119 e_float& operator=(boost::ulong_long_type c) { m_value = ::e_float(c); return *this; }
121 e_float& operator=(float c) { assign_large_real(c); return *this; }
122 e_float& operator=(double c) { assign_large_real(c); return *this; }
123 e_float& operator=(long double c) { assign_large_real(c); return *this; }
126 ::e_float& value(){ return m_value; }
127 ::e_float const& value()const{ return m_value; }
129 // Member arithmetic:
130 e_float& operator+=(const e_float& other)
131 { m_value += other.value(); return *this; }
132 e_float& operator-=(const e_float& other)
133 { m_value -= other.value(); return *this; }
134 e_float& operator*=(const e_float& other)
135 { m_value *= other.value(); return *this; }
136 e_float& operator/=(const e_float& other)
137 { m_value /= other.value(); return *this; }
138 e_float operator-()const
140 e_float const& operator+()const
147 void assign_large_real(const V& a)
153 m_value = ::ef::zero();
158 m_value = ::ef::one();
162 if ((boost::math::isinf)(a))
164 m_value = a > 0 ? m_value.my_value_inf() : -m_value.my_value_inf();
167 if((boost::math::isnan)(a))
169 m_value = m_value.my_value_nan();
176 m_value = ::ef::zero();
180 ::e_float shift = ::ef::pow2(30);
184 // extract 30 bits from f:
189 m_value += ::e_float(static_cast<INT64>(term));
192 m_value *= ::ef::pow2(e);
197 // Non-member arithmetic:
198 inline e_float operator+(const e_float& a, const e_float& b)
204 inline e_float operator-(const e_float& a, const e_float& b)
210 inline e_float operator*(const e_float& a, const e_float& b)
216 inline e_float operator/(const e_float& a, const e_float& b)
224 inline bool operator == (const e_float& a, const e_float& b)
225 { return a.value() == b.value() ? true : false; }
226 inline bool operator != (const e_float& a, const e_float& b)
227 { return a.value() != b.value() ? true : false;}
228 inline bool operator < (const e_float& a, const e_float& b)
229 { return a.value() < b.value() ? true : false; }
230 inline bool operator <= (const e_float& a, const e_float& b)
231 { return a.value() <= b.value() ? true : false; }
232 inline bool operator > (const e_float& a, const e_float& b)
233 { return a.value() > b.value() ? true : false; }
234 inline bool operator >= (const e_float& a, const e_float& b)
235 { return a.value() >= b.value() ? true : false; }
237 std::istream& operator >> (std::istream& is, e_float& f)
239 return is >> f.value();
242 std::ostream& operator << (std::ostream& os, const e_float& f)
244 return os << f.value();
247 inline e_float fabs(const e_float& v)
249 return ::ef::fabs(v.value());
252 inline e_float abs(const e_float& v)
254 return ::ef::fabs(v.value());
257 inline e_float floor(const e_float& v)
259 return ::ef::floor(v.value());
262 inline e_float ceil(const e_float& v)
264 return ::ef::ceil(v.value());
267 inline e_float pow(const e_float& v, const e_float& w)
269 return ::ef::pow(v.value(), w.value());
272 inline e_float pow(const e_float& v, int i)
274 return ::ef::pow(v.value(), ::e_float(i));
277 inline e_float exp(const e_float& v)
279 return ::ef::exp(v.value());
282 inline e_float log(const e_float& v)
284 return ::ef::log(v.value());
287 inline e_float sqrt(const e_float& v)
289 return ::ef::sqrt(v.value());
292 inline e_float sin(const e_float& v)
294 return ::ef::sin(v.value());
297 inline e_float cos(const e_float& v)
299 return ::ef::cos(v.value());
302 inline e_float tan(const e_float& v)
304 return ::ef::tan(v.value());
307 inline e_float acos(const e_float& v)
309 return ::ef::acos(v.value());
312 inline e_float asin(const e_float& v)
314 return ::ef::asin(v.value());
317 inline e_float atan(const e_float& v)
319 return ::ef::atan(v.value());
322 inline e_float atan2(const e_float& v, const e_float& u)
324 return ::ef::atan2(v.value(), u.value());
327 inline e_float ldexp(const e_float& v, int e)
329 return v.value() * ::ef::pow2(e);
332 inline e_float frexp(const e_float& v, int* expon)
336 v.value().extract_parts(d, i);
337 *expon = static_cast<int>(i);
338 return v.value() * ::ef::pow2(-i);
341 inline e_float sinh (const e_float& x)
343 return ::ef::sinh(x.value());
346 inline e_float cosh (const e_float& x)
348 return ::ef::cosh(x.value());
351 inline e_float tanh (const e_float& x)
353 return ::ef::tanh(x.value());
356 inline e_float asinh (const e_float& x)
358 return ::ef::asinh(x.value());
361 inline e_float acosh (const e_float& x)
363 return ::ef::acosh(x.value());
366 inline e_float atanh (const e_float& x)
368 return ::ef::atanh(x.value());
371 e_float fmod(const e_float& v1, const e_float& v2)
384 inline int fpclassify_imp< boost::math::ef::e_float> BOOST_NO_MACRO_EXPAND(boost::math::ef::e_float x, const generic_tag<true>&)
386 if(x.value().isnan())
388 if(x.value().isinf())
397 template <class Policy>
398 inline int itrunc(const e_float& v, const Policy& pol)
401 e_float r = boost::math::trunc(v, pol);
402 if(fabs(r) > (std::numeric_limits<int>::max)())
403 return static_cast<int>(policies::raise_rounding_error("boost::math::itrunc<%1%>(%1%)", 0, 0, v, pol));
404 return static_cast<int>(r.value().extract_int64());
407 template <class Policy>
408 inline long ltrunc(const e_float& v, const Policy& pol)
411 e_float r = boost::math::trunc(v, pol);
412 if(fabs(r) > (std::numeric_limits<long>::max)())
413 return static_cast<long>(policies::raise_rounding_error("boost::math::ltrunc<%1%>(%1%)", 0, 0L, v, pol));
414 return static_cast<long>(r.value().extract_int64());
417 #ifdef BOOST_HAS_LONG_LONG
418 template <class Policy>
419 inline boost::long_long_type lltrunc(const e_float& v, const Policy& pol)
422 e_float r = boost::math::trunc(v, pol);
423 if(fabs(r) > (std::numeric_limits<boost::long_long_type>::max)())
424 return static_cast<boost::long_long_type>(policies::raise_rounding_error("boost::math::lltrunc<%1%>(%1%)", 0, v, 0LL, pol).value().extract_int64());
425 return static_cast<boost::long_long_type>(r.value().extract_int64());
429 template <class Policy>
430 inline int iround(const e_float& v, const Policy& pol)
433 e_float r = boost::math::round(v, pol);
434 if(fabs(r) > (std::numeric_limits<int>::max)())
435 return static_cast<int>(policies::raise_rounding_error("boost::math::iround<%1%>(%1%)", 0, v, 0, pol).value().extract_int64());
436 return static_cast<int>(r.value().extract_int64());
439 template <class Policy>
440 inline long lround(const e_float& v, const Policy& pol)
443 e_float r = boost::math::round(v, pol);
444 if(fabs(r) > (std::numeric_limits<long>::max)())
445 return static_cast<long int>(policies::raise_rounding_error("boost::math::lround<%1%>(%1%)", 0, v, 0L, pol).value().extract_int64());
446 return static_cast<long int>(r.value().extract_int64());
449 #ifdef BOOST_HAS_LONG_LONG
450 template <class Policy>
451 inline boost::long_long_type llround(const e_float& v, const Policy& pol)
454 e_float r = boost::math::round(v, pol);
455 if(fabs(r) > (std::numeric_limits<boost::long_long_type>::max)())
456 return static_cast<boost::long_long_type>(policies::raise_rounding_error("boost::math::llround<%1%>(%1%)", 0, v, 0LL, pol).value().extract_int64());
457 return static_cast<boost::long_long_type>(r.value().extract_int64());
466 class numeric_limits< ::boost::math::ef::e_float> : public numeric_limits< ::e_float>
469 static const ::boost::math::ef::e_float (min) (void)
471 return (numeric_limits< ::e_float>::min)();
473 static const ::boost::math::ef::e_float (max) (void)
475 return (numeric_limits< ::e_float>::max)();
477 static const ::boost::math::ef::e_float epsilon (void)
479 return (numeric_limits< ::e_float>::epsilon)();
481 static const ::boost::math::ef::e_float round_error(void)
483 return (numeric_limits< ::e_float>::round_error)();
485 static const ::boost::math::ef::e_float infinity (void)
487 return (numeric_limits< ::e_float>::infinity)();
489 static const ::boost::math::ef::e_float quiet_NaN (void)
491 return (numeric_limits< ::e_float>::quiet_NaN)();
494 // e_float's supplied digits member is wrong
495 // - it should be same the same as digits 10
496 // - given that radix is 10.
498 static const int digits = digits10;
503 namespace boost{ namespace math{
507 template <class Policy>
508 struct precision< ::boost::math::ef::e_float, Policy>
510 typedef typename Policy::precision_type precision_type;
511 typedef digits2<((::std::numeric_limits< ::boost::math::ef::e_float>::digits10 + 1) * 1000L) / 301L> digits_2;
512 typedef typename mpl::if_c<
513 ((digits_2::value <= precision_type::value)
514 || (Policy::precision_type::value <= 0)),
515 // Default case, full precision for RealType:
517 // User customised precision:
527 inline int digits< ::boost::math::ef::e_float>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC( ::boost::math::ef::e_float))
529 return ((::std::numeric_limits< ::boost::math::ef::e_float>::digits10 + 1) * 1000L) / 301L;
533 inline ::boost::math::ef::e_float root_epsilon< ::boost::math::ef::e_float>()
535 return detail::root_epsilon_imp(static_cast< ::boost::math::ef::e_float const*>(0), mpl::int_<0>());
539 inline ::boost::math::ef::e_float forth_root_epsilon< ::boost::math::ef::e_float>()
541 return detail::forth_root_epsilon_imp(static_cast< ::boost::math::ef::e_float const*>(0), mpl::int_<0>());
548 template<class Policy>
549 struct lanczos<boost::math::ef::e_float, Policy>
551 typedef typename mpl::if_c<
552 std::numeric_limits< ::e_float>::digits10 < 22,
555 std::numeric_limits< ::e_float>::digits10 < 36,
558 std::numeric_limits< ::e_float>::digits10 < 50,
561 std::numeric_limits< ::e_float>::digits10 < 110,
570 } // namespace lanczos
572 template <class Policy>
573 inline boost::math::ef::e_float skewness(const extreme_value_distribution<boost::math::ef::e_float, Policy>& /*dist*/)
576 // This is 12 * sqrt(6) * zeta(3) / pi^3:
577 // See http://mathworld.wolfram.com/ExtremeValueDistribution.html
579 return boost::lexical_cast<boost::math::ef::e_float>("1.1395470994046486574927930193898461120875997958366");
582 template <class Policy>
583 inline boost::math::ef::e_float skewness(const rayleigh_distribution<boost::math::ef::e_float, Policy>& /*dist*/)
585 // using namespace boost::math::constants;
586 return boost::lexical_cast<boost::math::ef::e_float>("0.63111065781893713819189935154422777984404221106391");
587 // Computed using NTL at 150 bit, about 50 decimal digits.
588 // return 2 * root_pi<RealType>() * pi_minus_three<RealType>() / pow23_four_minus_pi<RealType>();
591 template <class Policy>
592 inline boost::math::ef::e_float kurtosis(const rayleigh_distribution<boost::math::ef::e_float, Policy>& /*dist*/)
594 // using namespace boost::math::constants;
595 return boost::lexical_cast<boost::math::ef::e_float>("3.2450893006876380628486604106197544154170667057995");
596 // Computed using NTL at 150 bit, about 50 decimal digits.
597 // return 3 - (6 * pi<RealType>() * pi<RealType>() - 24 * pi<RealType>() + 16) /
598 // (four_minus_pi<RealType>() * four_minus_pi<RealType>());
601 template <class Policy>
602 inline boost::math::ef::e_float kurtosis_excess(const rayleigh_distribution<boost::math::ef::e_float, Policy>& /*dist*/)
604 //using namespace boost::math::constants;
605 // Computed using NTL at 150 bit, about 50 decimal digits.
606 return boost::lexical_cast<boost::math::ef::e_float>("0.2450893006876380628486604106197544154170667057995");
607 // return -(6 * pi<RealType>() * pi<RealType>() - 24 * pi<RealType>() + 16) /
608 // (four_minus_pi<RealType>() * four_minus_pi<RealType>());
614 // Version of Digamma accurate to ~100 decimal digits.
616 template <class Policy>
617 boost::math::ef::e_float digamma_imp(boost::math::ef::e_float x, const mpl::int_<0>* , const Policy& pol)
620 // This handles reflection of negative arguments, and all our
621 // eboost::math::ef::e_floator handling, then forwards to the T-specific approximation.
623 BOOST_MATH_STD_USING // ADL of std functions.
625 boost::math::ef::e_float result = 0;
627 // Check for negative arguments and use reflection:
633 // Argument reduction for tan:
634 boost::math::ef::e_float remainder = x - floor(x);
635 // Shift to negative if > 0.5:
641 // check for evaluation at a negative pole:
645 return policies::raise_pole_error<boost::math::ef::e_float>("boost::math::digamma<%1%>(%1%)", 0, (1-x), pol);
647 result = constants::pi<boost::math::ef::e_float>() / tan(constants::pi<boost::math::ef::e_float>() * remainder);
649 result += big_digamma(x);
652 boost::math::ef::e_float bessel_i0(boost::math::ef::e_float x)
654 static const boost::math::ef::e_float P1[] = {
655 boost::lexical_cast<boost::math::ef::e_float>("-2.2335582639474375249e+15"),
656 boost::lexical_cast<boost::math::ef::e_float>("-5.5050369673018427753e+14"),
657 boost::lexical_cast<boost::math::ef::e_float>("-3.2940087627407749166e+13"),
658 boost::lexical_cast<boost::math::ef::e_float>("-8.4925101247114157499e+11"),
659 boost::lexical_cast<boost::math::ef::e_float>("-1.1912746104985237192e+10"),
660 boost::lexical_cast<boost::math::ef::e_float>("-1.0313066708737980747e+08"),
661 boost::lexical_cast<boost::math::ef::e_float>("-5.9545626019847898221e+05"),
662 boost::lexical_cast<boost::math::ef::e_float>("-2.4125195876041896775e+03"),
663 boost::lexical_cast<boost::math::ef::e_float>("-7.0935347449210549190e+00"),
664 boost::lexical_cast<boost::math::ef::e_float>("-1.5453977791786851041e-02"),
665 boost::lexical_cast<boost::math::ef::e_float>("-2.5172644670688975051e-05"),
666 boost::lexical_cast<boost::math::ef::e_float>("-3.0517226450451067446e-08"),
667 boost::lexical_cast<boost::math::ef::e_float>("-2.6843448573468483278e-11"),
668 boost::lexical_cast<boost::math::ef::e_float>("-1.5982226675653184646e-14"),
669 boost::lexical_cast<boost::math::ef::e_float>("-5.2487866627945699800e-18"),
671 static const boost::math::ef::e_float Q1[] = {
672 boost::lexical_cast<boost::math::ef::e_float>("-2.2335582639474375245e+15"),
673 boost::lexical_cast<boost::math::ef::e_float>("7.8858692566751002988e+12"),
674 boost::lexical_cast<boost::math::ef::e_float>("-1.2207067397808979846e+10"),
675 boost::lexical_cast<boost::math::ef::e_float>("1.0377081058062166144e+07"),
676 boost::lexical_cast<boost::math::ef::e_float>("-4.8527560179962773045e+03"),
677 boost::lexical_cast<boost::math::ef::e_float>("1.0"),
679 static const boost::math::ef::e_float P2[] = {
680 boost::lexical_cast<boost::math::ef::e_float>("-2.2210262233306573296e-04"),
681 boost::lexical_cast<boost::math::ef::e_float>("1.3067392038106924055e-02"),
682 boost::lexical_cast<boost::math::ef::e_float>("-4.4700805721174453923e-01"),
683 boost::lexical_cast<boost::math::ef::e_float>("5.5674518371240761397e+00"),
684 boost::lexical_cast<boost::math::ef::e_float>("-2.3517945679239481621e+01"),
685 boost::lexical_cast<boost::math::ef::e_float>("3.1611322818701131207e+01"),
686 boost::lexical_cast<boost::math::ef::e_float>("-9.6090021968656180000e+00"),
688 static const boost::math::ef::e_float Q2[] = {
689 boost::lexical_cast<boost::math::ef::e_float>("-5.5194330231005480228e-04"),
690 boost::lexical_cast<boost::math::ef::e_float>("3.2547697594819615062e-02"),
691 boost::lexical_cast<boost::math::ef::e_float>("-1.1151759188741312645e+00"),
692 boost::lexical_cast<boost::math::ef::e_float>("1.3982595353892851542e+01"),
693 boost::lexical_cast<boost::math::ef::e_float>("-6.0228002066743340583e+01"),
694 boost::lexical_cast<boost::math::ef::e_float>("8.5539563258012929600e+01"),
695 boost::lexical_cast<boost::math::ef::e_float>("-3.1446690275135491500e+01"),
696 boost::lexical_cast<boost::math::ef::e_float>("1.0"),
698 boost::math::ef::e_float value, factor, r;
701 using namespace boost::math::tools;
705 x = -x; // even function
709 return static_cast<boost::math::ef::e_float>(1);
711 if (x <= 15) // x in (0, 15]
713 boost::math::ef::e_float y = x * x;
714 value = evaluate_polynomial(P1, y) / evaluate_polynomial(Q1, y);
716 else // x in (15, \infty)
718 boost::math::ef::e_float y = 1 / x - boost::math::ef::e_float(1) / 15;
719 r = evaluate_polynomial(P2, y) / evaluate_polynomial(Q2, y);
720 factor = exp(x) / sqrt(x);
727 boost::math::ef::e_float bessel_i1(boost::math::ef::e_float x)
729 static const boost::math::ef::e_float P1[] = {
730 lexical_cast<boost::math::ef::e_float>("-1.4577180278143463643e+15"),
731 lexical_cast<boost::math::ef::e_float>("-1.7732037840791591320e+14"),
732 lexical_cast<boost::math::ef::e_float>("-6.9876779648010090070e+12"),
733 lexical_cast<boost::math::ef::e_float>("-1.3357437682275493024e+11"),
734 lexical_cast<boost::math::ef::e_float>("-1.4828267606612366099e+09"),
735 lexical_cast<boost::math::ef::e_float>("-1.0588550724769347106e+07"),
736 lexical_cast<boost::math::ef::e_float>("-5.1894091982308017540e+04"),
737 lexical_cast<boost::math::ef::e_float>("-1.8225946631657315931e+02"),
738 lexical_cast<boost::math::ef::e_float>("-4.7207090827310162436e-01"),
739 lexical_cast<boost::math::ef::e_float>("-9.1746443287817501309e-04"),
740 lexical_cast<boost::math::ef::e_float>("-1.3466829827635152875e-06"),
741 lexical_cast<boost::math::ef::e_float>("-1.4831904935994647675e-09"),
742 lexical_cast<boost::math::ef::e_float>("-1.1928788903603238754e-12"),
743 lexical_cast<boost::math::ef::e_float>("-6.5245515583151902910e-16"),
744 lexical_cast<boost::math::ef::e_float>("-1.9705291802535139930e-19"),
746 static const boost::math::ef::e_float Q1[] = {
747 lexical_cast<boost::math::ef::e_float>("-2.9154360556286927285e+15"),
748 lexical_cast<boost::math::ef::e_float>("9.7887501377547640438e+12"),
749 lexical_cast<boost::math::ef::e_float>("-1.4386907088588283434e+10"),
750 lexical_cast<boost::math::ef::e_float>("1.1594225856856884006e+07"),
751 lexical_cast<boost::math::ef::e_float>("-5.1326864679904189920e+03"),
752 lexical_cast<boost::math::ef::e_float>("1.0"),
754 static const boost::math::ef::e_float P2[] = {
755 lexical_cast<boost::math::ef::e_float>("1.4582087408985668208e-05"),
756 lexical_cast<boost::math::ef::e_float>("-8.9359825138577646443e-04"),
757 lexical_cast<boost::math::ef::e_float>("2.9204895411257790122e-02"),
758 lexical_cast<boost::math::ef::e_float>("-3.4198728018058047439e-01"),
759 lexical_cast<boost::math::ef::e_float>("1.3960118277609544334e+00"),
760 lexical_cast<boost::math::ef::e_float>("-1.9746376087200685843e+00"),
761 lexical_cast<boost::math::ef::e_float>("8.5591872901933459000e-01"),
762 lexical_cast<boost::math::ef::e_float>("-6.0437159056137599999e-02"),
764 static const boost::math::ef::e_float Q2[] = {
765 lexical_cast<boost::math::ef::e_float>("3.7510433111922824643e-05"),
766 lexical_cast<boost::math::ef::e_float>("-2.2835624489492512649e-03"),
767 lexical_cast<boost::math::ef::e_float>("7.4212010813186530069e-02"),
768 lexical_cast<boost::math::ef::e_float>("-8.5017476463217924408e-01"),
769 lexical_cast<boost::math::ef::e_float>("3.2593714889036996297e+00"),
770 lexical_cast<boost::math::ef::e_float>("-3.8806586721556593450e+00"),
771 lexical_cast<boost::math::ef::e_float>("1.0"),
773 boost::math::ef::e_float value, factor, r, w;
776 using namespace boost::math::tools;
781 return static_cast<boost::math::ef::e_float>(0);
783 if (w <= 15) // w in (0, 15]
785 boost::math::ef::e_float y = x * x;
786 r = evaluate_polynomial(P1, y) / evaluate_polynomial(Q1, y);
790 else // w in (15, \infty)
792 boost::math::ef::e_float y = 1 / w - boost::math::ef::e_float(1) / 15;
793 r = evaluate_polynomial(P2, y) / evaluate_polynomial(Q2, y);
794 factor = exp(w) / sqrt(w);
800 value *= -value; // odd function
805 } // namespace detail
808 #endif // BOOST_MATH_E_FLOAT_BINDINGS_HPP