//  Copyright (c) 2006 Xiaogang Zhang
//  Use, modification and distribution are subject to the
//  Boost Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

#ifndef BOOST_MATH_BESSEL_JY_HPP
#define BOOST_MATH_BESSEL_JY_HPP

#ifdef _MSC_VER
#pragma once
#endif

#include <boost/math/tools/config.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/special_functions/sign.hpp>
#include <boost/math/special_functions/hypot.hpp>
#include <boost/math/special_functions/sin_pi.hpp>
#include <boost/math/special_functions/cos_pi.hpp>
#include <boost/math/special_functions/detail/bessel_jy_asym.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/mpl/if.hpp>
#include <boost/type_traits/is_floating_point.hpp>
#include <complex>

// Bessel functions of the first and second kind of fractional order

namespace boost { namespace math {

namespace detail {

//
// Simultaneous calculation of A&S 9.2.9 and 9.2.10
// for use in A&S 9.2.5 and 9.2.6.
// This series is quick to evaluate, but divergent unless
// x is very large, in fact it's pretty hard to figure out
// with any degree of precision when this series actually 
// *will* converge!!  Consequently, we may just have to
// try it and see...
//
template <class T, class Policy>
bool hankel_PQ(T v, T x, T* p, T* q, const Policy& )
{
   BOOST_MATH_STD_USING
   T tolerance = 2 * policies::get_epsilon<T, Policy>();
   *p = 1;
   *q = 0;
   T k = 1;
   T z8 = 8 * x;
   T sq = 1;
   T mu = 4 * v * v;
   T term = 1;
   bool ok = true;
   do
   {
      term *= (mu - sq * sq) / (k * z8);
      *q += term;
      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("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("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;
   //std::cout << "C = " << Cr << " " << Ci << std::endl;
   //std::cout << "D = " << Dr << " " << Di << std::endl;
   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;
   //std::cout << fr << " " << fi << std::endl;
   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;
      //std::cout << "C = " << Cr << " " << Ci << std::endl;
      //std::cout << "D = " << Dr << " " << Di << std::endl;
      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;
      //std::cout << fr << " " << fi << std::endl;
   }
   policies::check_series_iterations("boost::math::bessel_jy<%1%>(%1%,%1%) in CF2_jy", k, pol);
   *p = fr;
   *q = fi;

   return 0;
}

enum
{
   need_j = 1, 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;

    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 (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
    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 <= 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;
        for (k = 1; k <= n; k++)            // forward recurrence for Y
        {
            next = 2 * (u + k) * current / x - 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 = W / (Yv * fv - Yv1);           // Wronskian relation
        }
        else
           Jv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
    }
    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;
        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;
           for (k = n; k > 0; k--)             // backward recurrence for J
           {
               next = 2 * (u + k) * current / x - prev;
               prev = current;
               current = next;
           }
           T ratio = (s * init) / current;     // scaling ratio
           // can also call CF1_jy() to get fu, not much difference in precision
           fu = prev / current;
           CF2_jy(u, x, &p, &q, pol);                  // continued fraction CF2_jy
           T t = u / x - fu;                   // t = J'/J
           gamma = (p - t) / q;
           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
           {
               next = 2 * (u + k) * current / x - 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)
    {
        T z = (u + n % 2);
        *J = boost::math::cos_pi(z, pol) * Jv - boost::math::sin_pi(z, pol) * Yv;     // reflection formula
        *Y = boost::math::sin_pi(z, pol) * Jv + boost::math::cos_pi(z, pol) * Yv;
    }
    else
    {
        *J = Jv;
        *Y = Yv;
    }

    return 0;
}

} // namespace detail

}} // namespaces

#endif // BOOST_MATH_BESSEL_JY_HPP