text: //  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 boost::math::ntl::RR 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 "ntl_rr_lanczos.hpp"
  
  #include <boost/math/special_functions/bessel.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;
  
     parameter_info<boost::math::ntl::RR> arg1, arg2;
     test_data<boost::math::ntl::RR> data;
  
     boost::math::ntl::RR::SetPrecision(1000); 
     boost::math::ntl::RR::SetOutputPrecision(40);
  
     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
           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");
        boost::math::ntl::RR (*fp)(boost::math::ntl::RR, boost::math::ntl::RR);
        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
           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;
     write_code(ofs, data, line.c_str());
  
     return 0;
  }
  
  
  
  
  