view release on metacpan or  search on metacpan

src/boost/math/special_functions/owens_t.hpp  view on Meta::CPAN

            {
               n = (std::numeric_limits<int>::max)();
            }
            n = (std::min)(n, 1500);
            RealType d = pow(3 + sqrt(RealType(8)), n);
            d = (d + 1 / d) / 2;
            RealType b = -1;
            RealType c = -d;
            int s = 1;

            for(int k = 0; k < n; ++k)
            {
               //
               // Check for both convergence and whether the series has gone bad:
               //
               if(
                  (fabs(z) > last_z)     // Series has gone divergent, abort
                  || (fabs(val) * tools::epsilon<RealType>() > fabs(c * s * z))  // Convergence!
                  || (z * s < 0)         // Series has stopped alternating - all bets are off - abort.
                  )
               {
                  break;
               }
               c = b - c;
               val += c * s * z;
               b = (k + n) * (k - n) * b / ((k + RealType(0.5)) * (k + 1));
               last_z = fabs(z);
               s = -s;
               z = y * ( vi - static_cast<RealType>(ii) * z );
               vi *= as;
               ii += 2;
            } // while( true )
            RealType err = fabs(c * z) / val;
            return std::pair<RealType, RealType>(val * exp( -hs*half<RealType>() ) / (d * root_two_pi<RealType>()), err);
         } // RealType owens_t_T2_accelerated(const RealType h, const RealType a, const RealType ah, const Policy&)

         template<typename RealType, typename Policy>
         inline RealType T4_mp(const RealType h, const RealType a, const Policy& pol)
         {
            BOOST_MATH_STD_USING
            
            const RealType hs = h*h;
            const RealType as = -a*a;

            unsigned short ii = 1;
            RealType ai = constants::one_div_two_pi<RealType>() * a * exp( -0.5*hs*(1.0-as) );
            RealType yi = 1.0;
            RealType val = 0.0;

            RealType lim = boost::math::policies::get_epsilon<RealType, Policy>();

            while( true )
            {
               RealType term = ai*yi;
               val += term;
               if((yi != 0) && (fabs(val * lim) > fabs(term)))
                  break;
               ii += 2;
               yi = (1.0-hs*yi) / static_cast<RealType>(ii);
               ai *= as;
               if(ii > (std::min)(1500, (int)policies::get_max_series_iterations<Policy>()))
                  policies::raise_evaluation_error("boost::math::owens_t<%1%>", 0, val, pol);
            } // while( true )

            return val;
         } // arg_type owens_t_T4(const arg_type h, const arg_type a, const unsigned short m)


         // This routine dispatches the call to one of six subroutines, depending on the values
         // of h and a.
         // preconditions: h >= 0, 0<=a<=1, ah=a*h
         //
         // Note there are different versions for different precisions....
         template<typename RealType, typename Policy>
         inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, mpl::int_<64> const&)
         {
            // Simple main case for 64-bit precision or less, this is as per the Patefield-Tandy paper:
            BOOST_MATH_STD_USING
            //
            // Handle some special cases first, these are from
            // page 1077 of Owen's original paper:
            //
            if(h == 0)
            {
               return atan(a) * constants::one_div_two_pi<RealType>();
            }
            if(a == 0)
            {
               return 0;
            }
            if(a == 1)
            {
               return owens_t_znorm2(RealType(-h)) * owens_t_znorm2(h) / 2;
            }
            if(a >= tools::max_value<RealType>())
            {
               return owens_t_znorm2(RealType(fabs(h)));
            }
            RealType val = 0; // avoid compiler warnings, 0 will be overwritten in any case
            const unsigned short icode = owens_t_compute_code(h, a);
            const unsigned short m = owens_t_get_order(icode, val /* just a dummy for the type */, pol);
            static const unsigned short meth[] = {1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 4, 4, 4, 4, 5, 6}; // 18 entries

            // determine the appropriate method, T1 ... T6
            switch( meth[icode] )
            {
            case 1: // T1
               val = owens_t_T1(h,a,m);
               break;
            case 2: // T2
               typedef typename policies::precision<RealType, Policy>::type precision_type;
               typedef mpl::bool_<(precision_type::value == 0) || (precision_type::value > 64)> tag_type;
               val = owens_t_T2(h, a, m, ah, pol, tag_type());
               break;
            case 3: // T3
               val = owens_t_T3(h,a,ah, pol);
               break;
            case 4: // T4
               val = owens_t_T4(h,a,m);
               break;
            case 5: // T5

src/boost/math/special_functions/owens_t.hpp  view on Meta::CPAN

               break;
            default:
               BOOST_THROW_EXCEPTION(std::logic_error("selection routine in Owen's T function failed"));
            }
            return val;
         }

         template<typename RealType, typename Policy>
         inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, const mpl::int_<65>&)
         {
            // Arbitrary precision version:
            BOOST_MATH_STD_USING
            //
            // Handle some special cases first, these are from
            // page 1077 of Owen's original paper:
            //
            if(h == 0)
            {
               return atan(a) * constants::one_div_two_pi<RealType>();
            }
            if(a == 0)
            {
               return 0;
            }
            if(a == 1)
            {
               return owens_t_znorm2(RealType(-h)) * owens_t_znorm2(h) / 2;
            }
            if(a >= tools::max_value<RealType>())
            {
               return owens_t_znorm2(RealType(fabs(h)));
            }
            // Attempt arbitrary precision code, this will throw if it goes wrong:
            typedef typename boost::math::policies::normalise<Policy, boost::math::policies::evaluation_error<> >::type forwarding_policy;
            std::pair<RealType, RealType> p1(0, tools::max_value<RealType>()), p2(0, tools::max_value<RealType>());
            RealType target_precision = policies::get_epsilon<RealType, Policy>() * 1000;
            bool have_t1(false), have_t2(false);
            if(ah < 3)
            {
               try
               {
                  have_t1 = true;
                  p1 = owens_t_T1_accelerated(h, a, forwarding_policy());
                  if(p1.second < target_precision)
                     return p1.first;
               }
               catch(const boost::math::evaluation_error&){}  // T1 may fail and throw, that's OK
            }
            if(ah > 1)
            {
               try
               {
                  have_t2 = true;
                  p2 = owens_t_T2_accelerated(h, a, ah, forwarding_policy());
                  if(p2.second < target_precision)
                     return p2.first;
               }
               catch(const boost::math::evaluation_error&){}  // T2 may fail and throw, that's OK
            }
            //
            // If we haven't tried T1 yet, do it now - sometimes it succeeds and the number of iterations
            // is fairly low compared to T4.
            //
            if(!have_t1)
            {
               try
               {
                  have_t1 = true;
                  p1 = owens_t_T1_accelerated(h, a, forwarding_policy());
                  if(p1.second < target_precision)
                     return p1.first;
               }
               catch(const boost::math::evaluation_error&){}  // T1 may fail and throw, that's OK
            }
            //
            // If we haven't tried T2 yet, do it now - sometimes it succeeds and the number of iterations
            // is fairly low compared to T4.
            //
            if(!have_t2)
            {
               try
               {
                  have_t2 = true;
                  p2 = owens_t_T2_accelerated(h, a, ah, forwarding_policy());
                  if(p2.second < target_precision)
                     return p2.first;
               }
               catch(const boost::math::evaluation_error&){}  // T2 may fail and throw, that's OK
            }
            //
            // OK, nothing left to do but try the most expensive option which is T4,
            // this is often slow to converge, but when it does converge it tends to
            // be accurate:
            try
            {
               return T4_mp(h, a, pol);
            }
            catch(const boost::math::evaluation_error&){}  // T4 may fail and throw, that's OK
            //
            // Now look back at the results from T1 and T2 and see if either gave better
            // results than we could get from the 64-bit precision versions.
            //
            if((std::min)(p1.second, p2.second) < 1e-20)
            {
               return p1.second < p2.second ? p1.first : p2.first;
            }
            //
            // We give up - no arbitrary precision versions succeeded!
            //
            return owens_t_dispatch(h, a, ah, pol, mpl::int_<64>());
         } // RealType owens_t_dispatch(RealType h, RealType a, RealType ah)
         template<typename RealType, typename Policy>
         inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, const mpl::int_<0>&)
         {
            // We don't know what the precision is until runtime:
            if(tools::digits<RealType>() <= 64)
               return owens_t_dispatch(h, a, ah, pol, mpl::int_<64>());
            return owens_t_dispatch(h, a, ah, pol, mpl::int_<65>());
         }
         template<typename RealType, typename Policy>
         inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol)
         {
            // Figure out the precision and forward to the correct version:
            typedef typename policies::precision<RealType, Policy>::type precision_type;
            typedef typename mpl::if_c<
               precision_type::value == 0,
               mpl::int_<0>,
               typename mpl::if_c<
                  precision_type::value <= 64,
                  mpl::int_<64>,
                  mpl::int_<65>
               >::type
            >::type tag_type;
            return owens_t_dispatch(h, a, ah, pol, tag_type());
         }
         // compute Owen's T function, T(h,a), for arbitrary values of h and a