view release on metacpan or  search on metacpan

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

      k += 1;
      sq += 2;
      T mult = (sq * sq - mu) / (k * z8);
      ok = fabs(mult) < 0.5f;
      term *= mult;
      *p += term;
      k += 1;
      sq += 2;
   }
   while((fabs(term) > tolerance * *p) && ok);
   return ok;
}

// Calculate Y(v, x) and Y(v+1, x) by Temme's method, see
// Temme, Journal of Computational Physics, vol 21, 343 (1976)
template <typename T, typename Policy>
int temme_jy(T v, T x, T* Y, T* Y1, const Policy& pol)
{
    T g, h, p, q, f, coef, sum, sum1, tolerance;
    T a, d, e, sigma;
    unsigned long k;

    BOOST_MATH_STD_USING
    using namespace boost::math::tools;
    using namespace boost::math::constants;

    BOOST_ASSERT(fabs(v) <= 0.5f);  // precondition for using this routine

    T gp = boost::math::tgamma1pm1(v, pol);
    T gm = boost::math::tgamma1pm1(-v, pol);
    T spv = boost::math::sin_pi(v, pol);
    T spv2 = boost::math::sin_pi(v/2, pol);
    T xp = pow(x/2, v);

    a = log(x / 2);
    sigma = -a * v;
    d = abs(sigma) < tools::epsilon<T>() ?
        T(1) : sinh(sigma) / sigma;
    e = abs(v) < tools::epsilon<T>() ? T(v*pi<T>()*pi<T>() / 2)
        : T(2 * spv2 * spv2 / v);

    T g1 = (v == 0) ? T(-euler<T>()) : T((gp - gm) / ((1 + gp) * (1 + gm) * 2 * v));
    T g2 = (2 + gp + gm) / ((1 + gp) * (1 + gm) * 2);
    T vspv = (fabs(v) < tools::epsilon<T>()) ? T(1/constants::pi<T>()) : T(v / spv);
    f = (g1 * cosh(sigma) - g2 * a * d) * 2 * vspv;

    p = vspv / (xp * (1 + gm));
    q = vspv * xp / (1 + gp);

    g = f + e * q;
    h = p;
    coef = 1;
    sum = coef * g;
    sum1 = coef * h;

    T v2 = v * v;
    T coef_mult = -x * x / 4;

    // series summation
    tolerance = policies::get_epsilon<T, Policy>();
    for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++)
    {
        f = (k * f + p + q) / (k*k - v2);
        p /= k - v;
        q /= k + v;
        g = f + e * q;
        h = p - k * g;
        coef *= coef_mult / k;
        sum += coef * g;
        sum1 += coef * h;
        if (abs(coef * g) < abs(sum) * tolerance) 
        { 
           break; 
        }
    }
    policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in temme_jy", k, pol);
    *Y = -sum;
    *Y1 = -2 * sum1 / x;

    return 0;
}

// Evaluate continued fraction fv = J_(v+1) / J_v, see
// Abramowitz and Stegun, Handbook of Mathematical Functions, 1972, 9.1.73
template <typename T, typename Policy>
int CF1_jy(T v, T x, T* fv, int* sign, const Policy& pol)
{
    T C, D, f, a, b, delta, tiny, tolerance;
    unsigned long k;
    int s = 1;

    BOOST_MATH_STD_USING

    // |x| <= |v|, CF1_jy converges rapidly
    // |x| > |v|, CF1_jy needs O(|x|) iterations to converge

    // modified Lentz's method, see
    // Lentz, Applied Optics, vol 15, 668 (1976)
    tolerance = 2 * policies::get_epsilon<T, Policy>();;
    tiny = sqrt(tools::min_value<T>());
    C = f = tiny;                           // b0 = 0, replace with tiny
    D = 0;
    for (k = 1; k < policies::get_max_series_iterations<Policy>() * 100; k++)
    {
        a = -1;
        b = 2 * (v + k) / x;
        C = b + a / C;
        D = b + a * D;
        if (C == 0) { C = tiny; }
        if (D == 0) { D = tiny; }
        D = 1 / D;
        delta = C * D;
        f *= delta;
        if (D < 0) { s = -s; }
        if (abs(delta - 1) < tolerance) 
        { break; }
    }
    policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF1_jy", k / 100, pol);
    *fv = -f;
    *sign = s;                              // sign of denominator

    return 0;
}
//
// This algorithm was originally written by Xiaogang Zhang
// using std::complex to perform the complex arithmetic.
// However, that turns out to 10x or more slower than using
// all real-valued arithmetic, so it's been rewritten using
// real values only.
//
template <typename T, typename Policy>
int CF2_jy(T v, T x, T* p, T* q, const Policy& pol)
{
   BOOST_MATH_STD_USING

   T Cr, Ci, Dr, Di, fr, fi, a, br, bi, delta_r, delta_i, temp;
   T tiny;
   unsigned long k;

   // |x| >= |v|, CF2_jy converges rapidly
   // |x| -> 0, CF2_jy fails to converge
   BOOST_ASSERT(fabs(x) > 1);

   // modified Lentz's method, complex numbers involved, see
   // Lentz, Applied Optics, vol 15, 668 (1976)
   T tolerance = 2 * policies::get_epsilon<T, Policy>();
   tiny = sqrt(tools::min_value<T>());
   Cr = fr = -0.5f / x;
   Ci = fi = 1;
   //Dr = Di = 0;
   T v2 = v * v;
   a = (0.25f - v2) / x; // Note complex this one time only!
   br = 2 * x;
   bi = 2;
   temp = Cr * Cr + 1;
   Ci = bi + a * Cr / temp;
   Cr = br + a / temp;
   Dr = br;
   Di = bi;
   if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; }
   if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; }
   temp = Dr * Dr + Di * Di;
   Dr = Dr / temp;
   Di = -Di / temp;
   delta_r = Cr * Dr - Ci * Di;
   delta_i = Ci * Dr + Cr * Di;
   temp = fr;
   fr = temp * delta_r - fi * delta_i;
   fi = temp * delta_i + fi * delta_r;
   for (k = 2; k < policies::get_max_series_iterations<Policy>(); k++)
   {
      a = k - 0.5f;
      a *= a;
      a -= v2;
      bi += 2;
      temp = Cr * Cr + Ci * Ci;
      Cr = br + a * Cr / temp;
      Ci = bi - a * Ci / temp;
      Dr = br + a * Dr;
      Di = bi + a * Di;
      if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; }
      if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; }
      temp = Dr * Dr + Di * Di;
      Dr = Dr / temp;
      Di = -Di / temp;
      delta_r = Cr * Dr - Ci * Di;
      delta_i = Ci * Dr + Cr * Di;
      temp = fr;
      fr = temp * delta_r - fi * delta_i;
      fi = temp * delta_i + fi * delta_r;
      if (fabs(delta_r - 1) + fabs(delta_i) < tolerance)
         break;
   }
   policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF2_jy", k, pol);
   *p = fr;
   *q = fi;

   return 0;
}

static const int need_j = 1;
static const int need_y = 2;

// Compute J(v, x) and Y(v, x) simultaneously by Steed's method, see
// Barnett et al, Computer Physics Communications, vol 8, 377 (1974)
template <typename T, typename Policy>
int bessel_jy(T v, T x, T* J, T* Y, int kind, const Policy& pol)
{
    BOOST_ASSERT(x >= 0);

    T u, Jv, Ju, Yv, Yv1, Yu, Yu1(0), fv, fu;
    T W, p, q, gamma, current, prev, next;
    bool reflect = false;
    unsigned n, k;
    int s;
    int org_kind = kind;
    T cp = 0;
    T sp = 0;

    static const char* function = "boost::math::bessel_jy<%1%>(%1%,%1%)";

    BOOST_MATH_STD_USING
    using namespace boost::math::tools;
    using namespace boost::math::constants;

    if (v < 0)
    {
        reflect = true;
        v = -v;                             // v is non-negative from here
        kind = need_j|need_y;               // need both for reflection formula
    }
    n = iround(v, pol);
    u = v - n;                              // -1/2 <= u < 1/2

    if(reflect)
    {
        T z = (u + n % 2);
        cp = boost::math::cos_pi(z, pol);
        sp = boost::math::sin_pi(z, pol);
    }

    if (x == 0)
    {
       *J = *Y = policies::raise_overflow_error<T>(
          function, 0, pol);
       return 1;
    }

    // x is positive until reflection
    W = T(2) / (x * pi<T>());               // Wronskian
    T Yv_scale = 1;
    if((x > 8) && (x < 1000) && hankel_PQ(v, x, &p, &q, pol))
    {
       //