view release on metacpan or  search on metacpan

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

         lambda = (a + b) * y - b;
      }
      if(lambda < 0)
      {
         std::swap(a, b);
         std::swap(x, y);
         invert = !invert;
         BOOST_MATH_INSTRUMENT_VARIABLE(invert);
      }
      
      if(b < 40)
      {
         if((floor(a) == a) && (floor(b) == b) && (a < (std::numeric_limits<int>::max)() - 100))
         {
            // relate to the binomial distribution and use a finite sum:
            T k = a - 1;
            T n = b + k;
            fract = binomial_ccdf(n, k, x, y);
            if(!normalised)
               fract *= boost::math::beta(a, b, pol);
            BOOST_MATH_INSTRUMENT_VARIABLE(fract);
         }
         else if(b * x <= 0.7)
         {
            if(!invert)
            {
               fract = ibeta_series(a, b, x, T(0), lanczos_type(), normalised, p_derivative, y, pol);
               BOOST_MATH_INSTRUMENT_VARIABLE(fract);
            }
            else
            {
               fract = -(normalised ? 1 : boost::math::beta(a, b, pol));
               invert = false;
               fract = -ibeta_series(a, b, x, fract, lanczos_type(), normalised, p_derivative, y, pol);
               BOOST_MATH_INSTRUMENT_VARIABLE(fract);
            }
         }
         else if(a > 15)
         {
            // sidestep so we can use the series representation:
            int n = itrunc(T(floor(b)), pol);
            if(n == b)
               --n;
            T bbar = b - n;
            T prefix;
            if(!normalised)
            {
               prefix = rising_factorial_ratio(T(a+bbar), bbar, n);
            }
            else
            {
               prefix = 1;
            }
            fract = ibeta_a_step(bbar, a, y, x, n, pol, normalised, static_cast<T*>(0));
            fract = beta_small_b_large_a_series(a,  bbar, x, y, fract, T(1), pol, normalised);
            fract /= prefix;
            BOOST_MATH_INSTRUMENT_VARIABLE(fract);
         }
         else if(normalised)
         {
            // the formula here for the non-normalised case is tricky to figure
            // out (for me!!), and requires two pochhammer calculations rather
            // than one, so leave it for now....
            int n = itrunc(T(floor(b)), pol);
            T bbar = b - n;
            if(bbar <= 0)
            {
               --n;
               bbar += 1;
            }
            fract = ibeta_a_step(bbar, a, y, x, n, pol, normalised, static_cast<T*>(0));
            fract += ibeta_a_step(a, bbar, x, y, 20, pol, normalised, static_cast<T*>(0));
            if(invert)
               fract -= (normalised ? 1 : boost::math::beta(a, b, pol));
            //fract = ibeta_series(a+20, bbar, x, fract, l, normalised, p_derivative, y);
            fract = beta_small_b_large_a_series(T(a+20),  bbar, x, y, fract, T(1), pol, normalised);
            if(invert)
            {
               fract = -fract;
               invert = false;
            }
            BOOST_MATH_INSTRUMENT_VARIABLE(fract);
         }
         else
         {
            fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative);
            BOOST_MATH_INSTRUMENT_VARIABLE(fract);
         }
      }
      else
      {
         fract = ibeta_fraction2(a, b, x, y, pol, normalised, p_derivative);
         BOOST_MATH_INSTRUMENT_VARIABLE(fract);
      }
   }
   if(p_derivative)
   {
      if(*p_derivative < 0)
      {
         *p_derivative = ibeta_power_terms(a, b, x, y, lanczos_type(), true, pol);
      }
      T div = y * x;

      if(*p_derivative != 0)
      {
         if((tools::max_value<T>() * div < *p_derivative))
         {
            // overflow, return an arbitarily large value:
            *p_derivative = tools::max_value<T>() / 2;
         }
         else
         {
            *p_derivative /= div;
         }
      }
   }
   return invert ? (normalised ? 1 : boost::math::beta(a, b, pol)) - fract : fract;
} // template <class T, class Lanczos>T ibeta_imp(T a, T b, T x, const Lanczos& l, bool inv, bool normalised)

template <class T, class Policy>
inline T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised)