123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355 |
- // 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;
- }
|