libs/math/tools/bessel_data.cpp

   1 //  Copyright (c) 2007 John Maddock
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5 //
   6 // Computes test data for the various bessel functions using
   7 // archived - deliberately naive - version of the code.
   8 // We'll rely on the high precision of mp_t to get us out of
   9 // trouble and not worry about how long the calculations take.
  10 // This provides a reasonably independent set of test data to
  11 // compare against newly added asymptotic expansions etc.
  12 //
  13 #include <fstream>
  14
  15 #include <boost/math/tools/test_data.hpp>
  16 #include <boost/math/special_functions/bessel.hpp>
  17 #include "mp_t.hpp"
  18
  19 using namespace boost::math::tools;
  20 using namespace boost::math;
  21 using namespace boost::math::detail;
  22 using namespace std;
  23
  24 // Compute J(v, x) and Y(v, x) simultaneously by Steed's method, see
  25 // Barnett et al, Computer Physics Communications, vol 8, 377 (1974)
  26 template <typename T>
  27 int bessel_jy_bare(T v, T x, T* J, T* Y, int kind = need_j|need_y)
  28 {
  29     // Jv1 = J_(v+1), Yv1 = Y_(v+1), fv = J_(v+1) / J_v
  30     // Ju1 = J_(u+1), Yu1 = Y_(u+1), fu = J_(u+1) / J_u
  31     T u, Jv, Ju, Yv, Yv1, Yu, Yu1, fv, fu;
  32     T W, p, q, gamma, current, prev, next;
  33     bool reflect = false;
  34     int n, k, s;
  35
  36     using namespace std;
  37     using namespace boost::math::tools;
  38     using namespace boost::math::constants;
  39
  40     if (v < 0)
  41     {
  42         reflect = true;
  43         v = -v;                             // v is non-negative from here
  44         kind = need_j|need_y;               // need both for reflection formula
  45     }
  46     n = real_cast<int>(v + 0.5L);
  47     u = v - n;                              // -1/2 <= u < 1/2
  48
  49     if (x < 0)
  50     {
  51        *J = *Y = policies::raise_domain_error<T>("",
  52           "Real argument x=%1% must be non-negative, complex number result not supported", x, policies::policy<>());
  53         return 1;
  54     }
  55     if (x == 0)
  56     {
  57        *J = *Y = policies::raise_overflow_error<T>(
  58           "", 0, policies::policy<>());
  59        return 1;
  60     }
  61
  62     // x is positive until reflection
  63     W = T(2) / (x * pi<T>());               // Wronskian
  64     if (x <= 2)                           // x in (0, 2]
  65     {
  66        if(temme_jy(u, x, &Yu, &Yu1, policies::policy<>()))             // Temme series
  67         {
  68            // domain error:
  69            *J = *Y = Yu;
  70            return 1;
  71         }
  72         prev = Yu;
  73         current = Yu1;
  74         for (k = 1; k <= n; k++)            // forward recurrence for Y
  75         {
  76             next = 2 * (u + k) * current / x - prev;
  77             prev = current;
  78             current = next;
  79         }
  80         Yv = prev;
  81         Yv1 = current;
  82         CF1_jy(v, x, &fv, &s, policies::policy<>());                 // continued fraction CF1
  83         Jv = W / (Yv * fv - Yv1);           // Wronskian relation
  84     }
  85     else                                    // x in (2, \infty)
  86     {
  87         // Get Y(u, x):
  88         CF1_jy(v, x, &fv, &s, policies::policy<>());
  89         // tiny initial value to prevent overflow
  90         T init = sqrt(tools::min_value<T>());
  91         prev = fv * s * init;
  92         current = s * init;
  93         for (k = n; k > 0; k--)             // backward recurrence for J
  94         {
  95             next = 2 * (u + k) * current / x - prev;
  96             prev = current;
  97             current = next;
  98         }
  99         T ratio = (s * init) / current;     // scaling ratio
 100         // can also call CF1() to get fu, not much difference in precision
 101         fu = prev / current;
 102         CF2_jy(u, x, &p, &q, policies::policy<>());                  // continued fraction CF2
 103         T t = u / x - fu;                   // t = J'/J
 104         gamma = (p - t) / q;
 105         Ju = sign(current) * sqrt(W / (q + gamma * (p - t)));
 106
 107         Jv = Ju * ratio;                    // normalization
 108
 109         Yu = gamma * Ju;
 110         Yu1 = Yu * (u/x - p - q/gamma);
 111
 112         // compute Y:
 113         prev = Yu;
 114         current = Yu1;
 115         for (k = 1; k <= n; k++)            // forward recurrence for Y
 116         {
 117             next = 2 * (u + k) * current / x - prev;
 118             prev = current;
 119             current = next;
 120         }
 121         Yv = prev;
 122     }
 123
 124     if (reflect)
 125     {
 126         T z = (u + n % 2) * pi<T>();
 127         *J = cos(z) * Jv - sin(z) * Yv;     // reflection formula
 128         *Y = sin(z) * Jv + cos(z) * Yv;
 129     }
 130     else
 131     {
 132         *J = Jv;
 133         *Y = Yv;
 134     }
 135
 136     return 0;
 137 }
 138
 139 int progress = 0;
 140
 141 template <class T>
 142 T cyl_bessel_j_bare(T v, T x)
 143 {
 144    T j, y;
 145    bessel_jy_bare(v, x, &j, &y);
 146
 147    std::cout << progress++ << ":   J(" << v << ", " << x << ") = " << j << std::endl;
 148
 149    if(fabs(j) > 1e30)
 150       throw std::domain_error("");
 151
 152    return j;
 153 }
 154
 155 template <class T>
 156 T cyl_bessel_i_bare(T v, T x)
 157 {
 158    using namespace std;
 159    if(x < 0)
 160    {
 161       // better have integer v:
 162       if(floor(v) == v)
 163       {
 164          T r = cyl_bessel_i_bare(v, -x);
 165          if(tools::real_cast<int>(v) & 1)
 166             r = -r;
 167          return r;
 168       }
 169       else
 170          return policies::raise_domain_error<T>(
 171             "",
 172             "Got x = %1%, but we need x >= 0", x, policies::policy<>());
 173    }
 174    if(x == 0)
 175    {
 176       return (v == 0) ? 1 : 0;
 177    }
 178    T I, K;
 179    boost::math::detail::bessel_ik(v, x, &I, &K, 0xffff, policies::policy<>());
 180
 181    std::cout << progress++ << ":   I(" << v << ", " << x << ") = " << I << std::endl;
 182
 183    if(fabs(I) > 1e30)
 184       throw std::domain_error("");
 185
 186    return I;
 187 }
 188
 189 template <class T>
 190 T cyl_bessel_k_bare(T v, T x)
 191 {
 192    using namespace std;
 193    if(x < 0)
 194    {
 195       return policies::raise_domain_error<T>(
 196          "",
 197          "Got x = %1%, but we need x > 0", x, policies::policy<>());
 198    }
 199    if(x == 0)
 200    {
 201       return (v == 0) ? policies::raise_overflow_error<T>("", 0, policies::policy<>())
 202          : policies::raise_domain_error<T>(
 203          "",
 204          "Got x = %1%, but we need x > 0", x, policies::policy<>());
 205    }
 206    T I, K;
 207    bessel_ik(v, x, &I, &K, 0xFFFF, policies::policy<>());
 208
 209    std::cout << progress++ << ":   K(" << v << ", " << x << ") = " << K << std::endl;
 210
 211    if(fabs(K) > 1e30)
 212       throw std::domain_error("");
 213
 214    return K;
 215 }
 216
 217 template <class T>
 218 T cyl_neumann_bare(T v, T x)
 219 {
 220    T j, y;
 221    bessel_jy(v, x, &j, &y, 0xFFFF, policies::policy<>());
 222
 223    std::cout << progress++ << ":   Y(" << v << ", " << x << ") = " << y << std::endl;
 224
 225    if(fabs(y) > 1e30)
 226       throw std::domain_error("");
 227
 228    return y;
 229 }
 230
 231 template <class T>
 232 T sph_bessel_j_bare(T v, T x)
 233 {
 234    std::cout << progress++ << ":   j(" << v << ", " << x << ") = ";
 235    if((v < 0) || (floor(v) != v))
 236       throw std::domain_error("");
 237    T r = sqrt(constants::pi<T>() / (2 * x)) * cyl_bessel_j_bare(v+0.5, x);
 238    std::cout << r << std::endl;
 239    return r;
 240 }
 241
 242 template <class T>
 243 T sph_bessel_y_bare(T v, T x)
 244 {
 245    std::cout << progress++ << ":   y(" << v << ", " << x << ") = ";
 246    if((v < 0) || (floor(v) != v))
 247       throw std::domain_error("");
 248    T r = sqrt(constants::pi<T>() / (2 * x)) * cyl_neumann_bare(v+0.5, x);
 249    std::cout << r << std::endl;
 250    return r;
 251 }
 252
 253 enum
 254 {
 255    func_J = 0,
 256    func_Y,
 257    func_I,
 258    func_K,
 259    func_j,
 260    func_y
 261 };
 262
 263 int main(int argc, char* argv[])
 264 {
 265    std::cout << std::setprecision(17) << std::scientific;
 266    std::cout << sph_bessel_j_bare(0., 0.1185395751953125e4) << std::endl;
 267    std::cout << sph_bessel_j_bare(22., 0.6540834903717041015625) << std::endl;
 268
 269    std::cout << std::setprecision(40) << std::scientific;
 270
 271    parameter_info<mp_t> arg1, arg2;
 272    test_data<mp_t> data;
 273
 274    int functype = 0;
 275    std::string letter = "J";
 276
 277    if(argc == 2)
 278    {
 279       if(std::strcmp(argv[1], "--Y") == 0)
 280       {
 281          functype = func_Y;
 282          letter = "Y";
 283       }
 284       else if(std::strcmp(argv[1], "--I") == 0)
 285       {
 286          functype = func_I;
 287          letter = "I";
 288       }
 289       else if(std::strcmp(argv[1], "--K") == 0)
 290       {
 291          functype = func_K;
 292          letter = "K";
 293       }
 294       else if(std::strcmp(argv[1], "--j") == 0)
 295       {
 296          functype = func_j;
 297          letter = "j";
 298       }
 299       else if(std::strcmp(argv[1], "--y") == 0)
 300       {
 301          functype = func_y;
 302          letter = "y";
 303       }
 304       else
 305          assert(0);
 306    }
 307
 308    bool cont;
 309    std::string line;
 310
 311    std::cout << "Welcome.\n"
 312       "This program will generate spot tests for the Bessel " << letter << " function\n\n";
 313    do{
 314       get_user_parameter_info(arg1, "v");
 315       get_user_parameter_info(arg2, "x");
 316       mp_t (*fp)(mp_t, mp_t);
 317       if(functype == func_J)
 318          fp = cyl_bessel_j_bare;
 319       else if(functype == func_I)
 320          fp = cyl_bessel_i_bare;
 321       else if(functype == func_K)
 322          fp = cyl_bessel_k_bare;
 323       else if(functype == func_Y)
 324          fp = cyl_neumann_bare;
 325       else if(functype == func_j)
 326          fp = sph_bessel_j_bare;
 327       else if(functype == func_y)
 328          fp = sph_bessel_y_bare;
 329       else
 330          assert(0);
 331
 332       data.insert(fp, arg1, arg2);
 333
 334       std::cout << "Any more data [y/n]?";
 335       std::getline(std::cin, line);
 336       boost::algorithm::trim(line);
 337       cont = (line == "y");
 338    }while(cont);
 339
 340    std::cout << "Enter name of test data file [default=bessel_j_data.ipp]";
 341    std::getline(std::cin, line);
 342    boost::algorithm::trim(line);
 343    if(line == "")
 344       line = "bessel_j_data.ipp";
 345    std::ofstream ofs(line.c_str());
 346    line.erase(line.find('.'));
 347    ofs << std::scientific << std::setprecision(40);
 348    write_code(ofs, data, line.c_str());
 349
 350    return 0;
 351 }
 352
 353
 354
 355