view release on metacpan or  search on metacpan

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

    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))
    {
       //
       // Hankel approximation: note that this method works best when x 
       // is large, but in that case we end up calculating sines and cosines
       // of large values, with horrendous resulting accuracy.  It is fast though
       // when it works....
       //
       T chi = x - fmod(T(v / 2 + 0.25f), T(2)) * boost::math::constants::pi<T>();
       T sc = sin(chi);
       T cc = cos(chi);
       chi = sqrt(2 / (boost::math::constants::pi<T>() * x));
       Yv = chi * (p * sc + q * cc);
       Jv = chi * (p * cc - q * sc);
    }
    else if((x < 1) && (u != 0) && (log(policies::get_epsilon<T, Policy>() / 2) > v * log((x/2) * (x/2) / v)))
    {
       // Evaluate using series representations.
       // This is particularly important for x << v as in this
       // area temme_jy may be slow to converge, if it converges at all.
       // Requires x is not an integer.
       if(kind&need_j)
          Jv = bessel_j_small_z_series(v, x, pol);
       else
          Jv = std::numeric_limits<T>::quiet_NaN();
       if((org_kind&need_y && (!reflect || (cp != 0))) 
          || (org_kind & need_j && (reflect && (sp != 0))))
       {
          // Only calculate if we need it, and if the reflection formula will actually use it:
          Yv = bessel_y_small_z_series(v, x, &Yv_scale, pol);
       }
       else
          Yv = std::numeric_limits<T>::quiet_NaN();
    }
    else if((u == 0) && (x < policies::get_epsilon<T, Policy>()))
    {
       // Truncated series evaluation for small x and v an integer,
       // much quicker in this area than temme_jy below.
       if(kind&need_j)
          Jv = bessel_j_small_z_series(v, x, pol);
       else
          Jv = std::numeric_limits<T>::quiet_NaN();
       if((org_kind&need_y && (!reflect || (cp != 0))) 
          || (org_kind & need_j && (reflect && (sp != 0))))
       {
          // Only calculate if we need it, and if the reflection formula will actually use it:
          Yv = bessel_yn_small_z(n, x, &Yv_scale, pol);
       }
       else
          Yv = std::numeric_limits<T>::quiet_NaN();
    }
    else if (x <= 2)                           // x in (0, 2]
    {
        if(temme_jy(u, x, &Yu, &Yu1, pol))             // Temme series
        {
           // domain error:
           *J = *Y = Yu;
           return 1;
        }
        prev = Yu;
        current = Yu1;
        T scale = 1;
        for (k = 1; k <= n; k++)            // forward recurrence for Y
        {
            T fact = 2 * (u + k) / x;
            if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))
            {
               scale /= current;
               prev /= current;
               current = 1;
            }
            next = fact * current - prev;
            prev = current;
            current = next;
        }
         Yv = prev;
         Yv1 = current;
         if(kind&need_j)
         {
            CF1_jy(v, x, &fv, &s, pol);                 // continued fraction CF1_jy
            Jv = scale * W / (Yv * fv - Yv1);           // Wronskian relation
         }
         else
            Jv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
         Yv_scale = scale;
    }
    else                                    // x in (2, \infty)
    {
        // Get Y(u, x):
        // define tag type that will dispatch to right limits:
        typedef typename bessel_asymptotic_tag<T, Policy>::type tag_type;

        T lim, ratio;
        switch(kind)
        {
        case need_j:
           lim = asymptotic_bessel_j_limit<T>(v, tag_type());
           break;
        case need_y:
           lim = asymptotic_bessel_y_limit<T>(tag_type());
           break;
        default:
           lim = (std::max)(
              asymptotic_bessel_j_limit<T>(v, tag_type()),
              asymptotic_bessel_y_limit<T>(tag_type()));
           break;
        }
        if(x > lim)
        {
           if(kind&need_y)
           {
              Yu = asymptotic_bessel_y_large_x_2(u, x);
              Yu1 = asymptotic_bessel_y_large_x_2(T(u + 1), x);
           }
           else
              Yu = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
           if(kind&need_j)
           {
              Jv = asymptotic_bessel_j_large_x_2(v, x);
           }
           else
              Jv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
        }
        else
        {
           CF1_jy(v, x, &fv, &s, pol);
           // tiny initial value to prevent overflow
           T init = sqrt(tools::min_value<T>());
           prev = fv * s * init;
           current = s * init;
           if(v < max_factorial<T>::value)
           {
              for (k = n; k > 0; k--)             // backward recurrence for J
              {
                  next = 2 * (u + k) * current / x - prev;
                  prev = current;
                  current = next;
              }
              ratio = (s * init) / current;     // scaling ratio
              // can also call CF1_jy() to get fu, not much difference in precision
              fu = prev / current;
           }
           else
           {
              //
              // When v is large we may get overflow in this calculation
              // leading to NaN's and other nasty surprises:
              //
              bool over = false;
              for (k = n; k > 0; k--)             // backward recurrence for J
              {
                  T t = 2 * (u + k) / x;
                  if(tools::max_value<T>() / t < current)
                  {
                     over = true;
                     break;
                  }
                  next = t * current - prev;
                  prev = current;
                  current = next;
              }
              if(!over)
              {
                 ratio = (s * init) / current;     // scaling ratio
                 // can also call CF1_jy() to get fu, not much difference in precision
                 fu = prev / current;
              }
              else
              {
                 ratio = 0;
                 fu = 1;
              }
           }
           CF2_jy(u, x, &p, &q, pol);                  // continued fraction CF2_jy
           T t = u / x - fu;                   // t = J'/J
           gamma = (p - t) / q;
           //
           // We can't allow gamma to cancel out to zero competely as it messes up
           // the subsequent logic.  So pretend that one bit didn't cancel out
           // and set to a suitably small value.  The only test case we've been able to
           // find for this, is when v = 8.5 and x = 4*PI.
           //
           if(gamma == 0)
           {
              gamma = u * tools::epsilon<T>() / x;
           }
           Ju = sign(current) * sqrt(W / (q + gamma * (p - t)));

           Jv = Ju * ratio;                    // normalization

           Yu = gamma * Ju;
           Yu1 = Yu * (u/x - p - q/gamma);
        }
        if(kind&need_y)
        {
           // compute Y:
           prev = Yu;
           current = Yu1;
           for (k = 1; k <= n; k++)            // forward recurrence for Y
           {
               T fact = 2 * (u + k) / x;
               if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))
               {
                  prev /= current;
                  Yv_scale /= current;
                  current = 1;
               }
               next = fact * current - prev;
               prev = current;
               current = next;
           }
           Yv = prev;
        }
        else
           Yv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
    }

    if (reflect)
    {
        if((sp != 0) && (tools::max_value<T>() * fabs(Yv_scale) < fabs(sp * Yv)))
           *J = org_kind & need_j ? T(-sign(sp) * sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
        else
            *J = cp * Jv - (sp == 0 ? T(0) : T((sp * Yv) / Yv_scale));     // reflection formula
        if((cp != 0) && (tools::max_value<T>() * fabs(Yv_scale) < fabs(cp * Yv)))
           *Y = org_kind & need_y ? T(-sign(cp) * sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
        else
           *Y = sp * Jv + (cp == 0 ? T(0) : T((cp * Yv) / Yv_scale));
    }
    else
    {
        *J = Jv;
        if(tools::max_value<T>() * fabs(Yv_scale) < fabs(Yv))
           *Y = org_kind & need_y ? T(sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
        else
         *Y = Yv / Yv_scale;
    }

    return 0;
}

} // namespace detail

}} // namespaces

#endif // BOOST_MATH_BESSEL_JY_HPP