view release on metacpan or  search on metacpan

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

            catch(...)
            {
               n = (std::numeric_limits<int>::max)();
            }
            n = (std::min)(n, 1500);
            T d = pow(3 + sqrt(T(8)), n);
            d = (d + 1 / d) / 2;
            T b = -1;
            T c = -d;
            c = b - c;
            sum *= c;
            b = -n * n * b * 2;
            abs_err = ldexp(fabs(sum), -tools::digits<T>());

            while(j < n)
            {
               a_pow *= aa;
               dj_pow *= half_h_h / j;
               one_minus_dj_sum += dj_pow;
               term = one_minus_dj_sum * a_pow / (2 * j + 1);
               c = b - c;
               sum += c * term;
               abs_err += ldexp((std::max)(T(fabs(sum)), T(fabs(c*term))), -tools::digits<T>());
               b = (j + n) * (j - n) * b / ((j + T(0.5)) * (j + 1));
               ++j;
               //
               // Include an escape route to prevent calculating too many terms:
               //
               if((j > 10) && (fabs(sum * tools::epsilon<T>()) > fabs(c * term)))
                  break;
            }
            abs_err += fabs(c * term);
            if(sum < 0)  // sum must always be positive, if it's negative something really bad has happend:
               policies::raise_evaluation_error(function, 0, T(0), pol);
            return std::pair<T, T>((sum / d) / boost::math::constants::two_pi<T>(), abs_err / sum);
         }

         template<typename RealType, class Policy>
         inline RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah, const Policy&, const mpl::true_&)
         {
            BOOST_MATH_STD_USING
            using namespace boost::math::constants;

            const unsigned short maxii = m+m+1;
            const RealType hs = h*h;
            const RealType as = -a*a;
            const RealType y = static_cast<RealType>(1) / hs;

            unsigned short ii = 1;
            RealType val = 0;
            RealType vi = a * exp( -ah*ah*half<RealType>() ) / root_two_pi<RealType>();
            RealType z = owens_t_znorm1(ah)/h;
            RealType last_z = fabs(z);
            RealType lim = policies::get_epsilon<RealType, Policy>();

            while( true )
            {
               val += z;
               //
               // This series stops converging after a while, so put a limit
               // on how far we go before returning our best guess:
               //
               if((fabs(lim * val) > fabs(z)) || ((ii > maxii) && (fabs(z) > last_z)) || (z == 0))
               {
                  val *= exp( -hs*half<RealType>() ) / root_two_pi<RealType>();
                  break;
               } // if( maxii <= ii )
               last_z = fabs(z);
               z = y * ( vi - static_cast<RealType>(ii) * z );
               vi *= as;
               ii += 2;
            } // while( true )

            return val;
         } // RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah)

         template<typename RealType, class Policy>
         inline std::pair<RealType, RealType> owens_t_T2_accelerated(const RealType h, const RealType a, const RealType ah, const Policy&)
         {
            //
            // This is the same series as T2, but with acceleration applied.
            // Note that we have to be *very* careful to check that nothing bad
            // has happened during evaluation - this series will go divergent
            // and/or fail to alternate at a drop of a hat! :-(
            //
            BOOST_MATH_STD_USING
            using namespace boost::math::constants;

            const RealType hs = h*h;
            const RealType as = -a*a;
            const RealType y = static_cast<RealType>(1) / hs;

            unsigned short ii = 1;
            RealType val = 0;
            RealType vi = a * exp( -ah*ah*half<RealType>() ) / root_two_pi<RealType>();
            RealType z = boost::math::detail::owens_t_znorm1(ah)/h;
            RealType last_z = fabs(z);

            //
            // Normally with this form of series acceleration we can calculate
            // up front how many terms will be required - based on the assumption
            // that each term decreases in size by a factor of 3.  However,
            // that assumption does not apply here, as the underlying T1 series can 
            // go quite strongly divergent in the early terms, before strongly
            // converging later.  Various "guestimates" have been tried to take account
            // of this, but they don't always work.... so instead set "n" to the 
            // largest value that won't cause overflow later, and abort iteration
            // when the last accelerated term was small enough...
            //
            int n;
            try
            {
               n = itrunc(RealType(tools::log_max_value<RealType>() / 6));
            }
            catch(...)
            {
               n = (std::numeric_limits<int>::max)();
            }
            n = (std::min)(n, 1500);
            RealType d = pow(3 + sqrt(RealType(8)), n);
            d = (d + 1 / d) / 2;