Updated boost to v1.55.0

[rsem.git] / boost / math / special_functions / detail / bessel_jn.hpp

diff --git a/boost/math/special_functions/detail/bessel_jn.hpp b/boost/math/special_functions/detail/bessel_jn.hpp
new file mode 100644 (file)
index 0000000..3f15f9c
--- /dev/null
+++ b/boost/math/special_functions/detail/bessel_jn.hpp
@@ -0,0 +1,134 @@
+//  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_JN_HPP
+#define BOOST_MATH_BESSEL_JN_HPP
+
+#ifdef _MSC_VER
+#pragma once
+#endif
+
+#include <boost/math/special_functions/detail/bessel_j0.hpp>
+#include <boost/math/special_functions/detail/bessel_j1.hpp>
+#include <boost/math/special_functions/detail/bessel_jy.hpp>
+#include <boost/math/special_functions/detail/bessel_jy_asym.hpp>
+#include <boost/math/special_functions/detail/bessel_jy_series.hpp>
+
+// Bessel function of the first kind of integer order
+// J_n(z) is the minimal solution
+// n < abs(z), forward recurrence stable and usable
+// n >= abs(z), forward recurrence unstable, use Miller's algorithm
+
+namespace boost { namespace math { namespace detail{
+
+template <typename T, typename Policy>
+T bessel_jn(int n, T x, const Policy& pol)
+{
+    T value(0), factor, current, prev, next;
+
+    BOOST_MATH_STD_USING
+
+    //
+    // Reflection has to come first:
+    //
+    if (n < 0)
+    {
+        factor = (n & 0x1) ? -1 : 1;  // J_{-n}(z) = (-1)^n J_n(z)
+        n = -n;
+    }
+    else
+    {
+        factor = 1;
+    }
+    if(x < 0)
+    {
+        factor *= (n & 0x1) ? -1 : 1;  // J_{n}(-z) = (-1)^n J_n(z)
+        x = -x;
+    }
+    //
+    // Special cases:
+    //
+    if (n == 0)
+    {
+        return factor * bessel_j0(x);
+    }
+    if (n == 1)
+    {
+        return factor * bessel_j1(x);
+    }
+
+    if (x == 0)                             // n >= 2
+    {
+        return static_cast<T>(0);
+    }
+
+    if(asymptotic_bessel_large_x_limit(T(n), x))
+      return factor * asymptotic_bessel_j_large_x_2<T>(n, x);
+
+    BOOST_ASSERT(n > 1);
+    T scale = 1;
+    if (n < abs(x))                         // forward recurrence
+    {
+        prev = bessel_j0(x);
+        current = bessel_j1(x);
+        policies::check_series_iterations<T>("boost::math::bessel_j_n<%1%>(%1%,%1%)", n, pol);
+        for (int k = 1; k < n; k++)
+        {
+            T fact = 2 * k / x;
+            //
+            // rescale if we would overflow or underflow:
+            //
+            if((fabs(fact) > 1) && ((tools::max_value<T>() - fabs(prev)) / fabs(fact) < fabs(current)))
+            {
+               scale /= current;
+               prev /= current;
+               current = 1;
+            }
+            value = fact * current - prev;
+            prev = current;
+            current = value;
+        }
+    }
+    else if((x < 1) || (n > x * x / 4) || (x < 5))
+    {
+       return factor * bessel_j_small_z_series(T(n), x, pol);
+    }
+    else                                    // backward recurrence
+    {
+        T fn; int s;                        // fn = J_(n+1) / J_n
+        // |x| <= n, fast convergence for continued fraction CF1
+        boost::math::detail::CF1_jy(static_cast<T>(n), x, &fn, &s, pol);
+        prev = fn;
+        current = 1;
+        // Check recursion won't go on too far:
+        policies::check_series_iterations<T>("boost::math::bessel_j_n<%1%>(%1%,%1%)", n, pol);
+        for (int k = n; k > 0; k--)
+        {
+            T fact = 2 * k / x;
+            if((fabs(fact) > 1) && ((tools::max_value<T>() - fabs(prev)) / fabs(fact) < fabs(current)))
+            {
+               prev /= current;
+               scale /= current;
+               current = 1;
+            }
+            next = fact * current - prev;
+            prev = current;
+            current = next;
+        }
+        value = bessel_j0(x) / current;       // normalization
+        scale = 1 / scale;
+    }
+    value *= factor;
+
+    if(tools::max_value<T>() * scale < fabs(value))
+       return policies::raise_overflow_error<T>("boost::math::bessel_jn<%1%>(%1%,%1%)", 0, pol);
+
+    return value / scale;
+}
+
+}}} // namespaces
+
+#endif // BOOST_MATH_BESSEL_JN_HPP
+