1 files changed, 355 insertions, 0 deletions
diff --git a/src/boost/libs/math/tools/bessel_data.cpp b/src/boost/libs/math/tools/bessel_data.cpp
new file mode 100644
index 000000000..6c9503587
--- /dev/null
+++ b/src/boost/libs/math/tools/bessel_data.cpp
@@ -0,0 +1,355 @@
+//  Copyright (c) 2007 John Maddock
+//  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)
+//
+// Computes test data for the various bessel functions using
+// archived - deliberately naive - version of the code.
+// We'll rely on the high precision of mp_t to get us out of
+// trouble and not worry about how long the calculations take.
+// This provides a reasonably independent set of test data to
+// compare against newly added asymptotic expansions etc.
+//
+#include <fstream>
+
+#include <boost/math/tools/test_data.hpp>
+#include <boost/math/special_functions/bessel.hpp>
+#include "mp_t.hpp"
+
+using namespace boost::math::tools;
+using namespace boost::math;
+using namespace boost::math::detail;
+using namespace std;
+
+// 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>
+int bessel_jy_bare(T v, T x, T* J, T* Y, int kind = need_j|need_y)
+{
+    // Jv1 = J_(v+1), Yv1 = Y_(v+1), fv = J_(v+1) / J_v
+    // Ju1 = J_(u+1), Yu1 = Y_(u+1), fu = J_(u+1) / J_u
+    T u, Jv, Ju, Yv, Yv1, Yu, Yu1, fv, fu;
+    T W, p, q, gamma, current, prev, next;
+    bool reflect = false;
+    int n, k, s;
+
+    using namespace std;
+    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 = real_cast<int>(v + 0.5L);
+    u = v - n;                              // -1/2 <= u < 1/2
+
+    if (x < 0)
+    {
+       *J = *Y = policies::raise_domain_error<T>("",
+          "Real argument x=%1% must be non-negative, complex number result not supported", x, policies::policy<>());
+        return 1;
+    }
+    if (x == 0)
+    {
+       *J = *Y = policies::raise_overflow_error<T>(
+          "", 0, policies::policy<>());
+       return 1;
+    }
+
+    // x is positive until reflection
+    W = T(2) / (x * pi<T>());               // Wronskian
+    if (x <= 2)                           // x in (0, 2]
+    {
+       if(temme_jy(u, x, &Yu, &Yu1, policies::policy<>()))             // 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;
+        CF1_jy(v, x, &fv, &s, policies::policy<>());                 // continued fraction CF1
+        Jv = W / (Yv * fv - Yv1);           // Wronskian relation
+    }
+    else                                    // x in (2, \infty)
+    {
+        // Get Y(u, x):
+        CF1_jy(v, x, &fv, &s, policies::policy<>());
+        // 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() to get fu, not much difference in precision
+        fu = prev / current;
+        CF2_jy(u, x, &p, &q, policies::policy<>());                  // continued fraction CF2
+        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);
+
+        // 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;
+    }
+
+    if (reflect)
+    {
+        T z = (u + n % 2) * pi<T>();
+        *J = cos(z) * Jv - sin(z) * Yv;     // reflection formula
+        *Y = sin(z) * Jv + cos(z) * Yv;
+    }
+    else
+    {
+        *J = Jv;
+        *Y = Yv;
+    }
+
+    return 0;
+}
+
+int progress = 0;
+
+template <class T>
+T cyl_bessel_j_bare(T v, T x)
+{
+   T j, y;
+   bessel_jy_bare(v, x, &j, &y);
+
+   std::cout << progress++ << ":   J(" << v << ", " << x << ") = " << j << std::endl;
+
+   if(fabs(j) > 1e30)
+      throw std::domain_error("");
+
+   return j;
+}
+
+template <class T>
+T cyl_bessel_i_bare(T v, T x)
+{
+   using namespace std;
+   if(x < 0)
+   {
+      // better have integer v:
+      if(floor(v) == v)
+      {
+         T r = cyl_bessel_i_bare(v, -x);
+         if(tools::real_cast<int>(v) & 1)
+            r = -r;
+         return r;
+      }
+      else
+         return policies::raise_domain_error<T>(
+            "",
+            "Got x = %1%, but we need x >= 0", x, policies::policy<>());
+   }
+   if(x == 0)
+   {
+      return (v == 0) ? 1 : 0;
+   }
+   T I, K;
+   boost::math::detail::bessel_ik(v, x, &I, &K, 0xffff, policies::policy<>());
+
+   std::cout << progress++ << ":   I(" << v << ", " << x << ") = " << I << std::endl;
+
+   if(fabs(I) > 1e30)
+      throw std::domain_error("");
+
+   return I;
+}
+
+template <class T>
+T cyl_bessel_k_bare(T v, T x)
+{
+   using namespace std;
+   if(x < 0)
+   {
+      return policies::raise_domain_error<T>(
+         "",
+         "Got x = %1%, but we need x > 0", x, policies::policy<>());
+   }
+   if(x == 0)
+   {
+      return (v == 0) ? policies::raise_overflow_error<T>("", 0, policies::policy<>())
+         : policies::raise_domain_error<T>(
+         "",
+         "Got x = %1%, but we need x > 0", x, policies::policy<>());
+   }
+   T I, K;
+   bessel_ik(v, x, &I, &K, 0xFFFF, policies::policy<>());
+
+   std::cout << progress++ << ":   K(" << v << ", " << x << ") = " << K << std::endl;
+
+   if(fabs(K) > 1e30)
+      throw std::domain_error("");
+
+   return K;
+}
+
+template <class T>
+T cyl_neumann_bare(T v, T x)
+{
+   T j, y;
+   bessel_jy(v, x, &j, &y, 0xFFFF, policies::policy<>());
+
+   std::cout << progress++ << ":   Y(" << v << ", " << x << ") = " << y << std::endl;
+
+   if(fabs(y) > 1e30)
+      throw std::domain_error("");
+
+   return y;
+}
+
+template <class T>
+T sph_bessel_j_bare(T v, T x)
+{
+   std::cout << progress++ << ":   j(" << v << ", " << x << ") = ";
+   if((v < 0) || (floor(v) != v))
+      throw std::domain_error("");
+   T r = sqrt(constants::pi<T>() / (2 * x)) * cyl_bessel_j_bare(v+0.5, x);
+   std::cout << r << std::endl;
+   return r;
+}
+
+template <class T>
+T sph_bessel_y_bare(T v, T x)
+{
+   std::cout << progress++ << ":   y(" << v << ", " << x << ") = ";
+   if((v < 0) || (floor(v) != v))
+      throw std::domain_error("");
+   T r = sqrt(constants::pi<T>() / (2 * x)) * cyl_neumann_bare(v+0.5, x);
+   std::cout << r << std::endl;
+   return r;
+}
+
+enum
+{
+   func_J = 0,
+   func_Y,
+   func_I,
+   func_K,
+   func_j,
+   func_y
+};
+
+int main(int argc, char* argv[])
+{
+   std::cout << std::setprecision(17) << std::scientific;
+   std::cout << sph_bessel_j_bare(0., 0.1185395751953125e4) << std::endl;
+   std::cout << sph_bessel_j_bare(22., 0.6540834903717041015625) << std::endl;
+
+   std::cout << std::setprecision(40) << std::scientific;
+
+   parameter_info<mp_t> arg1, arg2;
+   test_data<mp_t> data;
+
+   int functype = 0;
+   std::string letter = "J";
+
+   if(argc == 2)
+   {
+      if(std::strcmp(argv[1], "--Y") == 0)
+      {
+         functype = func_Y;
+         letter = "Y";
+      }
+      else if(std::strcmp(argv[1], "--I") == 0)
+      {
+         functype = func_I;
+         letter = "I";
+      }
+      else if(std::strcmp(argv[1], "--K") == 0)
+      {
+         functype = func_K;
+         letter = "K";
+      }
+      else if(std::strcmp(argv[1], "--j") == 0)
+      {
+         functype = func_j;
+         letter = "j";
+      }
+      else if(std::strcmp(argv[1], "--y") == 0)
+      {
+         functype = func_y;
+         letter = "y";
+      }
+      else
+         BOOST_ASSERT(0);
+   }
+
+   bool cont;
+   std::string line;
+
+   std::cout << "Welcome.\n"
+      "This program will generate spot tests for the Bessel " << letter << " function\n\n";
+   do{
+      get_user_parameter_info(arg1, "v");
+      get_user_parameter_info(arg2, "x");
+      mp_t (*fp)(mp_t, mp_t);
+      if(functype == func_J)
+         fp = cyl_bessel_j_bare;
+      else if(functype == func_I)
+         fp = cyl_bessel_i_bare;
+      else if(functype == func_K)
+         fp = cyl_bessel_k_bare;
+      else if(functype == func_Y)
+         fp = cyl_neumann_bare;
+      else if(functype == func_j)
+         fp = sph_bessel_j_bare;
+      else if(functype == func_y)
+         fp = sph_bessel_y_bare;
+      else
+         BOOST_ASSERT(0);
+
+      data.insert(fp, arg1, arg2);
+
+      std::cout << "Any more data [y/n]?";
+      std::getline(std::cin, line);
+      boost::algorithm::trim(line);
+      cont = (line == "y");
+   }while(cont);
+
+   std::cout << "Enter name of test data file [default=bessel_j_data.ipp]";
+   std::getline(std::cin, line);
+   boost::algorithm::trim(line);
+   if(line == "")
+      line = "bessel_j_data.ipp";
+   std::ofstream ofs(line.c_str());
+   line.erase(line.find('.'));
+   ofs << std::scientific << std::setprecision(40);
+   write_code(ofs, data, line.c_str());
+
+   return 0;
+}
+
+
+
+