diff --git a/include/boost/math/special_functions/bessel.hpp b/include/boost/math/special_functions/bessel.hpp
index e7309baab..c458b6c9d 100644
--- a/include/boost/math/special_functions/bessel.hpp
+++ b/include/boost/math/special_functions/bessel.hpp
@@ -363,12 +363,13 @@ inline T cyl_bessel_j_zero_imp(T v, int m, const Policy& pol)
 
    const T half_epsilon(boost::math::tools::epsilon<T>() / 2U);
 
-   // Handle the non-finite order.
+   // Handle non-finite order.
    if (!(boost::math::isfinite)(v) )
    {
      return policies::raise_domain_error<T>(function, "Order argument is %1%, but must be finite >= 0 !", v, pol);
    }
 
+   // Handle negative rank.
    if(m < 0)
    {
       // Zeros of Jv(x) with negative rank are not defined and requesting one raises a domain error.
@@ -415,12 +416,14 @@ inline T cyl_bessel_j_zero_imp(T v, int m, const Policy& pol)
    // Account for the radix of number representations having non-two radix!
    const int my_digits2 = policies::digits<T, Policy>();
 
+   const T delta_lo = ((guess_root > 0.2F) ? T(0.2) : T(guess_root / 2U));
+
    // Perform the root-finding using Newton-Raphson iteration from Boost.Math.
    const T jvm =
       boost::math::tools::newton_raphson_iterate(
          boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::function_object_jv_and_jv_prime<T, Policy>((order_is_integer ? vv : v), order_is_zero, pol),
          guess_root,
-         (std::max)(T(guess_root - 0.2F), T(0)),
+         T(guess_root - delta_lo),
          T(guess_root + 0.2F),
          my_digits2,
          number_of_iterations);
@@ -447,7 +450,7 @@ inline T cyl_neumann_zero_imp(T v, int m, const Policy& pol)
      return policies::raise_domain_error<T>(function, "Order argument is %1%, but must be finite >= 0 !", v, pol);
    }
 
-   // Handle if the zero'th root or a negative root is requested.
+   // Handle negative rank.
    if(m < 0)
    {
       return policies::raise_domain_error<T>(function, "Requested the %1%'th zero, but the rank must be positive !", m, pol);
@@ -459,10 +462,19 @@ inline T cyl_neumann_zero_imp(T v, int m, const Policy& pol)
    const bool order_is_negative = (v < 0);
    const T vv((!order_is_negative) ? v : -v);
 
-   // Check if the order is very close to a negative half-integer.
-   const bool order_is_negative_half_integer = (order_is_negative && (((vv * 2U) - floor(vv * 2U)) < boost::math::tools::epsilon<T>()));
+   const bool order_is_integer = ((vv - floor(vv)) < half_epsilon);
 
-   // The zero'th zero of Jv(x) for v < 0 is not defined
+   // For negative integers, use reflection to positive integer order.
+   if(order_is_negative && order_is_integer)
+      return boost::math::detail::cyl_neumann_zero_imp(vv, m, pol);
+
+   // Check if the order is very close to a negative half-integer.
+   const T delta_half_integer(vv - (floor(vv) + 0.5F));
+
+   const bool order_is_negative_half_integer =
+      (order_is_negative && ((delta_half_integer > -half_epsilon) && (delta_half_integer < +half_epsilon)));
+
+   // The zero'th zero of Yv(x) for v < 0 is not defined
    // unless the order is a negative integer.
    if((m == 0) && (!order_is_negative_half_integer))
    {
@@ -470,7 +482,7 @@ inline T cyl_neumann_zero_imp(T v, int m, const Policy& pol)
       return policies::raise_domain_error<T>(function, "Requested the %1%'th zero of Yv for negative, non-half-integer order, but the rank must be > 0 !", m, pol);
    }
 
-   // For nrgative half-integers, use the corresponding
+   // For negative half-integers, use the corresponding
    // spherical Bessel function of positive half-integer order.
    if(order_is_negative_half_integer)
       return boost::math::detail::cyl_bessel_j_zero_imp(vv, m, pol);
@@ -487,12 +499,14 @@ inline T cyl_neumann_zero_imp(T v, int m, const Policy& pol)
    // Account for the radix of number representations having non-two radix!
    const int my_digits2 = policies::digits<T, Policy>();
 
+   const T delta_lo = ((guess_root > 0.2F) ? T(0.2) : T(guess_root / 2U));
+
    // Perform the root-finding using Newton-Raphson iteration from Boost.Math.
    const T yvm =
       boost::math::tools::newton_raphson_iterate(
          boost::math::detail::bessel_zero::cyl_neumann_zero_detail::function_object_yv_and_yv_prime<T, Policy>(v, pol),
          guess_root,
-         (std::max)(T(guess_root - 0.2F), T(0)),
+         T(guess_root - delta_lo),
          T(guess_root + 0.2F),
          my_digits2,
          number_of_iterations);
diff --git a/include/boost/math/special_functions/detail/bessel_jy_zero.hpp b/include/boost/math/special_functions/detail/bessel_jy_zero.hpp
index ef04afec1..917530173 100644
--- a/include/boost/math/special_functions/detail/bessel_jy_zero.hpp
+++ b/include/boost/math/special_functions/detail/bessel_jy_zero.hpp
@@ -9,9 +9,11 @@
 //
 // This header contains implementation details for estimating the zeros
 // of cylindrical Bessel and Neumann functions on the positive real axis.
+// Support is included for both positive as well as negative order.
 // Various methods are used to estimate the roots. These include
 // empirical curve fitting and McMahon's asymptotic approximation
-// for small order, and uniform asymptotic expansion for large order.
+// for small order, uniform asymptotic expansion for large order,
+// and iteration and root interlacing for negative order.
 //
 #ifndef _BESSEL_JY_ZERO_2013_01_18_HPP_
   #define _BESSEL_JY_ZERO_2013_01_18_HPP_
@@ -283,7 +285,14 @@
                 root_hi = root_lo;
 
                 // Decrease the lower end of the bracket using an adaptive algorithm.
-                (root_lo > 0.5F) ? root_lo -= 0.5F : root_lo *= 0.85F;
+                if(root_lo > 0.5F)
+                {
+                  root_lo -= 0.5F;
+                }
+                else
+                {
+                  root_lo *= 0.85F;
+                }
               }
             }
             else
@@ -437,7 +446,7 @@
         template<class T, class Policy>
         T initial_guess(const T& v, const int m, const Policy& pol)
         {
-          BOOST_MATH_STD_USING // ADL of std names, needed for ceil.
+          BOOST_MATH_STD_USING // ADL of std names, needed for floor.
 
           // Compute an estimate of the m'th root of cyl_neumann.
 
@@ -448,12 +457,17 @@
           {
             // Create the positive order and extract its positive floor and ceiling integer parts.
             const T vv(-v);
-            const T vv_ceil (ceil (vv));
             const T vv_floor(floor(vv));
 
             // The to-be-found root is bracketed by the roots of the
             // Bessel function whose reflected, positive integer order
-            // is greater than, but nearest to vv.
+            // is less than, but nearest to vv.
+
+            // The special case of negative, half-integer order uses
+            // the relation between Yv and spherical Bessel functions
+            // in order to obtain the bracket for the root.
+            // In these special cases, cyl_neumann(-n/2, x) = sph_bessel_j(+n/2, x)
+            // for v = -n/2.
 
             T root_hi;
             T root_lo;
@@ -462,6 +476,8 @@
             {
               // The estimate of the first root for negative order is found using
               // an adaptive range-searching algorithm.
+              // Take special precautions for the discontinuity at negative,
+              // half-integer orders and use different brackets above and below these.
               if(T(vv - vv_floor) < 0.5F)
               {
                 root_hi = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(vv_floor, m, pol);
@@ -487,7 +503,14 @@
                 root_hi = root_lo;
 
                 // Decrease the lower end of the bracket using an adaptive algorithm.
-                (root_lo > 0.5F) ? root_lo -= 0.5F : root_lo *= 0.85F;
+                if(root_lo > 0.5F)
+                {
+                  root_lo -= 0.5F;
+                }
+                else
+                {
+                  root_lo *= 0.85F;
+                }
               }
             }
             else