From 3c8349e7eca33ff2b2cd22d14922ee0d5fbd01ab Mon Sep 17 00:00:00 2001 From: Jens Maurer Date: Wed, 26 Sep 2001 21:29:40 +0000 Subject: [PATCH] update from Hubert Holin now works with significantly more compilers [SVN r11284] --- Jamfile | 14 + include/boost/math/octonion.hpp | 11 + include/boost/math/quaternion.hpp | 11 + .../boost/math/special_functions/acosh.hpp | 129 +++ .../boost/math/special_functions/asinh.hpp | 98 ++ .../boost/math/special_functions/atanh.hpp | 89 +- include/boost/math/special_functions/sinc.hpp | 41 +- .../boost/math/special_functions/sinhc.hpp | 41 +- octonion/octonion_test.cpp | 139 +-- octonion/output.txt | Bin 3802 -> 3965 bytes quaternion/HSO3.hpp | 883 +++++++++--------- quaternion/HSO3SO4.cpp | 810 ++++++++-------- quaternion/HSO4.hpp | 270 +++--- quaternion/output.txt | Bin 2203 -> 2528 bytes quaternion/quaternion_test.cpp | 381 ++++---- special_functions/index.html | 7 +- special_functions/output.txt | 330 +------ special_functions/special_functions.html | 59 +- special_functions/special_functions_test.cpp | 342 ++++--- 19 files changed, 1891 insertions(+), 1764 deletions(-) create mode 100644 Jamfile create mode 100644 include/boost/math/special_functions/acosh.hpp create mode 100644 include/boost/math/special_functions/asinh.hpp diff --git a/Jamfile b/Jamfile new file mode 100644 index 000000000..4241544a5 --- /dev/null +++ b/Jamfile @@ -0,0 +1,14 @@ +subproject libs/math ; + +exe octonion_test : octonion/octonion_test.cpp + : $(BOOST_ROOT) : ; + +exe quaternion_test : quaternion/quaternion_test.cpp + : $(BOOST_ROOT) : ; + +exe HSO3SO4 : quaternion/HSO3SO4.cpp + : $(BOOST_ROOT) : ; + +exe special_functions_test : special_functions/special_functions_test.cpp + : $(BOOST_ROOT) : ; + diff --git a/include/boost/math/octonion.hpp b/include/boost/math/octonion.hpp index 7d456d468..0788b99b8 100644 --- a/include/boost/math/octonion.hpp +++ b/include/boost/math/octonion.hpp @@ -18,6 +18,17 @@ namespace boost { namespace math { +#if defined(__GNUC__) && (__GNUC__ < 3) + // gcc 2.x ignores function scope using declarations, + // put them in the scope of the enclosing namespace instead: + using ::std::valarray; + using ::std::sqrt; + using ::std::cos; + using ::std::sin; + using ::std::exp; + using ::std::cosh; +#endif + #define BOOST_OCTONION_ACCESSOR_GENERATOR(type) \ type real() const \ { \ diff --git a/include/boost/math/quaternion.hpp b/include/boost/math/quaternion.hpp index a7ee5389f..817d7e8e2 100644 --- a/include/boost/math/quaternion.hpp +++ b/include/boost/math/quaternion.hpp @@ -23,6 +23,17 @@ namespace boost { namespace math { +#if defined(__GNUC__) && (__GNUC__ < 3) + // gcc 2.x ignores function scope using declarations, + // put them in the scope of the enclosing namespace instead: + using ::std::valarray; + using ::std::sqrt; + using ::std::cos; + using ::std::sin; + using ::std::exp; + using ::std::cosh; +#endif + #define BOOST_QUATERNION_ACCESSOR_GENERATOR(type) \ type real() const \ { \ diff --git a/include/boost/math/special_functions/acosh.hpp b/include/boost/math/special_functions/acosh.hpp new file mode 100644 index 000000000..848f4ac57 --- /dev/null +++ b/include/boost/math/special_functions/acosh.hpp @@ -0,0 +1,129 @@ +// boost asinh.hpp header file + +// (C) Copyright Eric Ford 2001 & Hubert Holin. Permission to copy, use, modify, sell and +// distribute this software is granted provided this copyright notice appears +// in all copies. This software is provided "as is" without express or implied +// warranty, and with no claim as to its suitability for any purpose. + +// See http://www.boost.org for updates, documentation, and revision history. + +#ifndef BOOST_ACOSH_HPP +#define BOOST_ACOSH_HPP + + +#include +#include +#include +#include + + +// This is the inverse of the hyperbolic cosine function. + +namespace boost +{ + namespace math + { +#if defined(__GNUC__) && (__GNUC__ < 3) + // gcc 2.x ignores function scope using declarations, + // put them in the scope of the enclosing namespace instead: + + using ::std::abs; + using ::std::sqrt; + using ::std::log; + + using ::std::numeric_limits; +#endif + + // These are implementation details (for main fare see below) + + namespace detail + { + template < + typename T, + bool QuietNanSupported + > + struct acosh_helper2_t + { + static T get_NaN() + { + return(::std::numeric_limits::quiet_NaN()); + } + }; // boost::detail::acosh_helper2_t + + + template + struct acosh_helper2_t + { + static T get_NaN() + { + ::std::string error_reporting("Argument to acosh is greater than or equal to +1!"); + ::std::domain_error bad_argument(error_reporting); + + throw(bad_argument); + } + }; // boost::detail::acosh_helper2_t + + } // boost::detail + + + // This is the main fare + + template + inline T acosh(const T x) + { + using ::std::abs; + using ::std::sqrt; + using ::std::log; + + using ::std::numeric_limits; + + typedef detail::acosh_helper2_t::has_quiet_NaN> helper2_type; + + + T const one = static_cast(1); + T const two = static_cast(2); + + static T const taylor_2_bound = sqrt(numeric_limits::epsilon()); + static T const taylor_n_bound = sqrt(taylor_2_bound); + static T const upper_taylor_2_bound = one/taylor_2_bound; + + if (x < one) + { + return(helper2_type::get_NaN()); + } + else if (x >= taylor_n_bound) + { + if (x > upper_taylor_2_bound) + { + // approximation by laurent series in 1/x at 0+ order from -1 to 0 + return( log( x*two) ); + } + else + { + return( log( x + sqrt(x*x-one) ) ); + } + } + else + { + T y = sqrt(x-one); + + // approximation by taylor series in y at 0 up to order 2 + T result = y; + + if (y >= taylor_2_bound) + { + T y3 = y*y*y; + + // approximation by taylor series in y at 0 up to order 4 + result -= y3/static_cast(12); + } + + return(sqrt(static_cast(2))*result); + } + } + } +} + +#endif /* BOOST_ACOSH_HPP */ + + diff --git a/include/boost/math/special_functions/asinh.hpp b/include/boost/math/special_functions/asinh.hpp new file mode 100644 index 000000000..77fed9182 --- /dev/null +++ b/include/boost/math/special_functions/asinh.hpp @@ -0,0 +1,98 @@ +// boost asinh.hpp header file + +// (C) Copyright Eric Ford & Hubert Holin 2001. Permission to copy, use, modify, sell and +// distribute this software is granted provided this copyright notice appears +// in all copies. This software is provided "as is" without express or implied +// warranty, and with no claim as to its suitability for any purpose. + +// See http://www.boost.org for updates, documentation, and revision history. + +#ifndef BOOST_ASINH_HPP +#define BOOST_ASINH_HPP + + +#include +#include +#include +#include + + +// This is the inverse of the hyperbolic sine function. + +namespace boost +{ + namespace math + { +#if defined(__GNUC__) && (__GNUC__ < 3) + // gcc 2.x ignores function scope using declarations, + // put them in the scope of the enclosing namespace instead: + + using ::std::abs; + using ::std::sqrt; + using ::std::log; + + using ::std::numeric_limits; +#endif + + template + inline T asinh(const T x) + { + using ::std::abs; + using ::std::sqrt; + using ::std::log; + + using ::std::numeric_limits; + + + T const one = static_cast(1); + T const two = static_cast(2); + + static T const taylor_2_bound = sqrt(numeric_limits::epsilon()); + static T const taylor_n_bound = sqrt(taylor_2_bound); + static T const upper_taylor_2_bound = one/taylor_2_bound; + static T const upper_taylor_n_bound = one/taylor_n_bound; + + if (x >= +taylor_n_bound) + { + if (x > upper_taylor_n_bound) + { + if (x > upper_taylor_2_bound) + { + // approximation by laurent series in 1/x at 0+ order from -1 to 0 + return( log( x * two) ); + } + else + { + // approximation by laurent series in 1/x at 0+ order from -1 to 1 + return( log( x*two + (one/(x*two)) ) ); + } + } + else + { + return( log( x + sqrt(x*x+one) ) ); + } + } + else if (x <= -taylor_n_bound) + { + return(-asinh(-x)); + } + else + { + // approximation by taylor series in x at 0 up to order 2 + T result = x; + + if (abs(x) >= taylor_2_bound) + { + T x3 = x*x*x; + + // approximation by taylor series in x at 0 up to order 4 + result -= x3/static_cast(6); + } + + return(result); + } + } + } +} + +#endif /* BOOST_ASINH_HPP */ diff --git a/include/boost/math/special_functions/atanh.hpp b/include/boost/math/special_functions/atanh.hpp index fe5567a99..ec23ed89e 100644 --- a/include/boost/math/special_functions/atanh.hpp +++ b/include/boost/math/special_functions/atanh.hpp @@ -23,6 +23,17 @@ namespace boost { namespace math { +#if defined(__GNUC__) && (__GNUC__ < 3) + // gcc 2.x ignores function scope using declarations, + // put them in the scope of the enclosing namespace instead: + + using ::std::abs; + using ::std::sqrt; + using ::std::log; + + using ::std::numeric_limits; +#endif + // These are implementation details (for main fare see below) namespace detail @@ -35,12 +46,12 @@ namespace boost { static T get_pos_infinity() { - return(+std::numeric_limits::infinity()); + return(+::std::numeric_limits::infinity()); } static T get_neg_infinity() { - return(-std::numeric_limits::infinity()); + return(-::std::numeric_limits::infinity()); } }; // boost::math::detail::atanh_helper1_t @@ -51,7 +62,7 @@ namespace boost static T get_pos_infinity() { ::std::string error_reporting("Argument to atanh is +1 (result: +Infinity)!"); - ::std::out_of_range bad_argument(error_reporting); + ::std::out_of_range bad_argument(error_reporting); throw(bad_argument); } @@ -59,7 +70,7 @@ namespace boost static T get_neg_infinity() { ::std::string error_reporting("Argument to atanh is -1 (result: -Infinity)!"); - ::std::out_of_range bad_argument(error_reporting); + ::std::out_of_range bad_argument(error_reporting); throw(bad_argument); } @@ -72,14 +83,9 @@ namespace boost > struct atanh_helper2_t { - static T get_pos_NaN() + static T get_NaN() { - return(+std::numeric_limits::quiet_NaN()); - } - - static T get_neg_NaN() - { - return(-std::numeric_limits::quiet_NaN()); + return(::std::numeric_limits::quiet_NaN()); } }; // boost::detail::atanh_helper2_t @@ -87,18 +93,10 @@ namespace boost template struct atanh_helper2_t { - static T get_pos_NaN() + static T get_NaN() { - ::std::string error_reporting("Argument to atanh is strictly greater than +1!"); - ::std::domain_error bad_argument(error_reporting); - - throw(bad_argument); - } - - static T get_neg_NaN() - { - ::std::string error_reporting("Argument to atanh is strictly smaller than -1!"); - ::std::domain_error bad_argument(error_reporting); + ::std::string error_reporting("Argument to atanh is strictly greater than +1 or strictly smaller than -1!"); + ::std::domain_error bad_argument(error_reporting); throw(bad_argument); } @@ -113,44 +111,55 @@ namespace boost { using ::std::abs; using ::std::sqrt; + using ::std::log; using ::std::numeric_limits; - typedef detail::atanh_helper1_t::has_infinity> helper1_type; - typedef detail::atanh_helper2_t::has_quiet_NaN> helper2_type; + typedef detail::atanh_helper1_t::has_infinity> helper1_type; + typedef detail::atanh_helper2_t::has_quiet_NaN> helper2_type; - static T const e0 = sqrt(numeric_limits::epsilon()); + T const one = static_cast(1); + T const two = static_cast(2); - T const one = static_cast(1); - T const two = static_cast(2); + static T const taylor_2_bound = sqrt(numeric_limits::epsilon()); + static T const taylor_n_bound = sqrt(taylor_2_bound); - if (x < -one) + if (x < -one) { - return(helper2_type::get_neg_NaN()); + return(helper2_type::get_NaN()); } - else if (x == -one) + else if (x < -one+numeric_limits::epsilon()) { return(helper1_type::get_neg_infinity()); } - else if (x == +one) + else if (x > +one-numeric_limits::epsilon()) { return(helper1_type::get_pos_infinity()); } - else if (x > +one) + else if (x > +one) { - return(helper2_type::get_pos_NaN()); + return(helper2_type::get_NaN()); } - else if (abs(x) < e0) + else if (abs(x) >= taylor_n_bound) { - return(x); - } - else // -one < x <= -e0 or +e0 <= x < +one - { - using ::std::log; - return(log( (one + x) / (one - x) ) / two); } + else + { + // approximation by taylor series in x at 0 up to order 2 + T result = x; + + if (abs(x) >= taylor_2_bound) + { + T x3 = x*x*x; + + // approximation by taylor series in x at 0 up to order 4 + result += x3/static_cast(3); + } + + return(result); + } } } } diff --git a/include/boost/math/special_functions/sinc.hpp b/include/boost/math/special_functions/sinc.hpp index 18d7bfa5f..538cf7c5a 100644 --- a/include/boost/math/special_functions/sinc.hpp +++ b/include/boost/math/special_functions/sinc.hpp @@ -23,6 +23,17 @@ namespace boost { namespace math { +#if defined(__GNUC__) && (__GNUC__ < 3) + // gcc 2.x ignores function scope using declarations, + // put them in the scope of the enclosing namespace instead: + + using ::std::abs; + using ::std::sqrt; + using ::std::sin; + + using ::std::numeric_limits; +#endif + // This is the "Sinus Cardinal" of index Pi. template @@ -34,26 +45,29 @@ namespace boost using ::std::numeric_limits; - static T const e1 = numeric_limits::epsilon(); - static T const e2 = sqrt(e1); - static T const e3 = sqrt(e2); + static T const taylor_0_bound = numeric_limits::epsilon(); + static T const taylor_2_bound = sqrt(taylor_0_bound); + static T const taylor_n_bound = sqrt(taylor_2_bound); - if (abs(x) > e3) + if (abs(x) >= taylor_n_bound) { return(sin(x)/x); } else { + // approximation by taylor series in x at 0 up to order 0 T result = static_cast(1); - if (abs(x) > e1) + if (abs(x) >= taylor_0_bound) { T x2 = x*x; + // approximation by taylor series in x at 0 up to order 2 result -= x2/static_cast(6); - if (abs(x) > e2) + if (abs(x) >= taylor_2_bound) { + // approximation by taylor series in x at 0 up to order 4 result += (x2*x2)/static_cast(120); } } @@ -72,26 +86,29 @@ namespace boost using ::std::numeric_limits; - static T const e1 = numeric_limits::epsilon(); - static T const e2 = sqrt(e1); - static T const e3 = sqrt(e2); + static T const taylor_0_bound = numeric_limits::epsilon(); + static T const taylor_2_bound = sqrt(taylor_0_bound); + static T const taylor_n_bound = sqrt(taylor_2_bound); - if (abs(x) > e3) + if (abs(x) >= taylor_n_bound) { return(sin(x)/x); } else { + // approximation by taylor series in x at 0 up to order 0 U result = static_cast< U >(1); - if (abs(x) > e1) + if (abs(x) >= taylor_0_bound) { U x2 = x*x; + // approximation by taylor series in x at 0 up to order 2 result -= x2/static_cast(6); - if (abs(x) > e2) + if (abs(x) >= taylor_2_bound) { + // approximation by taylor series in x at 0 up to order 4 result += (x2*x2)/static_cast(120); } } diff --git a/include/boost/math/special_functions/sinhc.hpp b/include/boost/math/special_functions/sinhc.hpp index bf5363464..4ad863852 100644 --- a/include/boost/math/special_functions/sinhc.hpp +++ b/include/boost/math/special_functions/sinhc.hpp @@ -23,6 +23,17 @@ namespace boost { namespace math { +#if defined(__GNUC__) && (__GNUC__ < 3) + // gcc 2.x ignores function scope using declarations, + // put them in the scope of the enclosing namespace instead: + + using ::std::abs; + using ::std::sqrt; + using ::std::sinh; + + using ::std::numeric_limits; +#endif + // This is the "Hyperbolic Sinus Cardinal" of index Pi. template @@ -34,26 +45,29 @@ namespace boost using ::std::numeric_limits; - static T const e1 = numeric_limits::epsilon(); - static T const e2 = sqrt(e1); - static T const e3 = sqrt(e2); + static T const taylor_0_bound = numeric_limits::epsilon(); + static T const taylor_2_bound = sqrt(taylor_0_bound); + static T const taylor_n_bound = sqrt(taylor_2_bound); - if (abs(x) > e3) + if (abs(x) >= taylor_n_bound) { return(sinh(x)/x); } else { + // approximation by taylor series in x at 0 up to order 0 T result = static_cast(1); - if (abs(x) > e1) + if (abs(x) >= taylor_0_bound) { T x2 = x*x; + // approximation by taylor series in x at 0 up to order 2 result += x2/static_cast(6); - if (abs(x) > e2) + if (abs(x) >= taylor_2_bound) { + // approximation by taylor series in x at 0 up to order 4 result += (x2*x2)/static_cast(120); } } @@ -72,26 +86,29 @@ namespace boost using ::std::numeric_limits; - static T const e1 = numeric_limits::epsilon(); - static T const e2 = sqrt(e1); - static T const e3 = sqrt(e2); + static T const taylor_0_bound = numeric_limits::epsilon(); + static T const taylor_2_bound = sqrt(taylor_0_bound); + static T const taylor_n_bound = sqrt(taylor_2_bound); - if (abs(x) > e3) + if (abs(x) >= taylor_n_bound) { return(sinh(x)/x); } else { + // approximation by taylor series in x at 0 up to order 0 U result = static_cast< U >(1); - if (abs(x) > e1) + if (abs(x) >= taylor_0_bound) { U x2 = x*x; + // approximation by taylor series in x at 0 up to order 2 result += x2/static_cast(6); - if (abs(x) > e2) + if (abs(x) >= taylor_2_bound) { + // approximation by taylor series in x at 0 up to order 4 result += (x2*x2)/static_cast(120); } } diff --git a/octonion/octonion_test.cpp b/octonion/octonion_test.cpp index e4363eae9..3f179a09a 100644 --- a/octonion/octonion_test.cpp +++ b/octonion/octonion_test.cpp @@ -13,13 +13,6 @@ #include -#include -#ifdef BOOST_NO_STDC_NAMESPACE -namespace std { - using ::sqrt; - using ::atan; -} -#endif // explicit (if ludicrous) instanciation #ifndef __GNUC__ @@ -321,7 +314,7 @@ int test_main(int, char *[]) << o0 << " ." << ::std::endl; } #else - ::std::istringstream bogus("(1,2,3,4,5,6,7,8)"); + ::std::istringstream bogus("(1,2,3,4,5,6,7,8)"); bogus >> o0; @@ -404,14 +397,14 @@ int test_main(int, char *[]) ::std::cout << ::std::endl; - float rho = ::std::sqrt(4096.0f); - float theta = ::std::atan(1.0f); - float phi1 = ::std::atan(1.0f); - float phi2 = ::std::atan(1.0f); - float phi3 = ::std::atan(1.0f); - float phi4 = ::std::atan(1.0f); - float phi5 = ::std::atan(1.0f); - float phi6 = ::std::atan(1.0f); + float rho = ::std::sqrt(4096.0f); + float theta = ::std::atan(1.0f); + float phi1 = ::std::atan(1.0f); + float phi2 = ::std::atan(1.0f); + float phi3 = ::std::atan(1.0f); + float phi4 = ::std::atan(1.0f); + float phi5 = ::std::atan(1.0f); + float phi6 = ::std::atan(1.0f); ::std::cout << "The value of the octonion represented " << "in spherical form by " @@ -424,14 +417,14 @@ int test_main(int, char *[]) phi1, phi2, phi3, phi4, phi5, phi6) << ::std::endl; - float rho1 = 1; - float rho2 = 2; - float rho3 = ::std::sqrt(2.0f); - float rho4 = ::std::sqrt(8.0f); - float theta1 = 0; - float theta2 = ::std::atan(1.0f)*2; - float theta3 = ::std::atan(1.0f); - float theta4 = ::std::atan(::std::sqrt(3.0f)); + float rho1 = 1; + float rho2 = 2; + float rho3 = ::std::sqrt(2.0f); + float rho4 = ::std::sqrt(8.0f); + float theta1 = 0; + float theta2 = ::std::atan(1.0f)*2; + float theta3 = ::std::atan(1.0f); + float theta4 = ::std::atan(::std::sqrt(3.0f)); ::std::cout << "The value of the octonion represented " << "in multipolar form by " @@ -444,14 +437,14 @@ int test_main(int, char *[]) rho3, theta3, rho4, theta4) << ::std::endl; - float r = ::std::sqrt(2.0f); - float angle = ::std::atan(1.0f); - float h1 = 3; - float h2 = 4; - float h3 = 5; - float h4 = 6; - float h5 = 7; - float h6 = 8; + float r = ::std::sqrt(2.0f); + float angle = ::std::atan(1.0f); + float h1 = 3; + float h2 = 4; + float h3 = 5; + float h4 = 6; + float h5 = 7; + float h6 = 8; ::std::cout << "The value of the octonion represented " << "in cylindrical form by " @@ -612,40 +605,56 @@ int test_main(int, char *[]) using ::boost::math::abs; - #define BOOST_OCTONION_MULTIPLICATION_TEST(type) \ - \ - BOOST_CRITICAL_TEST(abs( \ - ::boost::math::octonion(1,0,0,0,0,0,0,0)* \ - ::boost::math::octonion(1,0,0,0,0,0,0,0)- \ - static_cast(1)) <= \ - numeric_limits::epsilon()); \ - \ - for (int idx = 1; idx < 8; ++idx) \ - { \ - ::boost::math::octonion toto = \ - index_i_element(idx); \ - \ - BOOST_CRITICAL_TEST(abs(toto*toto+static_cast(1)) <= \ - numeric_limits::epsilon()); \ - } +#define BOOST_OCTONION_MULTIPLICATION_TEST(type) \ + \ + ::std::cout << "Testing multiplication." << std::endl; \ + \ + BOOST_CRITICAL_TEST(abs( \ + ::boost::math::octonion(1,0,0,0,0,0,0,0)* \ + ::boost::math::octonion(1,0,0,0,0,0,0,0)- \ + static_cast(1)) <= \ + numeric_limits::epsilon()); \ + \ + for (int idx = 1; idx < 8; ++idx) \ + { \ + ::boost::math::octonion toto = \ + index_i_element(idx); \ + \ + BOOST_CRITICAL_TEST(abs(toto*toto+static_cast(1)) <= \ + numeric_limits::epsilon()); \ + } - #define BOOST_OCTONION_TRENSCENDENTALS_TEST(type) \ - \ - for (int idx = 1; idx < 8; ++idx) \ - { \ - ::boost::math::octonion toto = \ - static_cast(4)*::std::atan(static_cast(1))* \ - index_i_element(idx); \ - \ - BOOST_TEST(abs(exp(toto)+static_cast(1)) <= \ - 2*numeric_limits::epsilon()); \ - } +#define BOOST_OCTONION_TRENSCENDENTALS_TEST(type) \ + \ + ::std::cout << "Testing exp." << std::endl; \ + \ + for (int idx = 1; idx < 8; ++idx) \ + { \ + ::boost::math::octonion toto = \ + static_cast(4)*::std::atan(static_cast(1))* \ + index_i_element(idx); \ + \ + BOOST_TEST(abs(exp(toto)+static_cast(1)) <= \ + 2*numeric_limits::epsilon()); \ + } - #define BOOST_OCTONION_TEST(type) \ - BOOST_OCTONION_MULTIPLICATION_TEST(type) \ - BOOST_OCTONION_TRENSCENDENTALS_TEST(type) +#if defined(__GNUC__) || defined(__COMO__) || \ + (defined(__MWERKS__) && (__MWERKS__ <= 0x2301)) +#define BOOST_OCTONION_TEST(type) \ + \ + ::std::cout << "Testing " << #type << "." << std::endl; \ + \ + BOOST_OCTONION_MULTIPLICATION_TEST(type) +#else +#define BOOST_OCTONION_TEST(type) \ + \ + ::std::cout << "Testing " << #type << "." << std::endl; \ + \ + BOOST_OCTONION_MULTIPLICATION_TEST(type) \ + BOOST_OCTONION_TRENSCENDENTALS_TEST(type) +#endif BOOST_OCTONION_TEST(float) @@ -653,10 +662,10 @@ int test_main(int, char *[]) BOOST_OCTONION_TEST(long double) - #undef BOOST_OCTONION_TEST +#undef BOOST_OCTONION_TEST - #undef BOOST_OCTONION_MULTIPLICATION_TEST - #undef BOOST_OCTONION_TRENSCENDENTALS_TEST +#undef BOOST_OCTONION_MULTIPLICATION_TEST +#undef BOOST_OCTONION_TRENSCENDENTALS_TEST return(::boost::exit_success); } diff --git a/octonion/output.txt b/octonion/output.txt index 1f7a02086d3f31dcb489bfdff20f58ec1914e3f1..181fdd7898cf34fd1bb5a06f93394a2e471d5912 100644 GIT binary patch delta 191 zcmca5`&VwmalZPH)Z&uNymW=MoczQRJuVnGw=}0DvmhrkIk6-&KMx_ET2X-Dq~w<- q<)k9`L}|*&&r4T;>okX2k5_A6zCvnIQGQXeLP}~$YH~?x3Ksw$ct*zn delta 27 icmew>cT0A|aX#U^e1+7aqWq#_g_P8i)Z~)X6fOXv*b1=# diff --git a/quaternion/HSO3.hpp b/quaternion/HSO3.hpp index 1e904e297..f100d3ac2 100644 --- a/quaternion/HSO3.hpp +++ b/quaternion/HSO3.hpp @@ -25,473 +25,480 @@ #include + +#if defined(__GNUC__) && (__GNUC__ < 3) +// gcc 2.x ignores function scope using declarations, put them here instead: +using namespace ::std; +using namespace ::boost::math; +#endif + template struct R3_matrix { - TYPE_FLOAT a11, a12, a13; - TYPE_FLOAT a21, a22, a23; - TYPE_FLOAT a31, a32, a33; + TYPE_FLOAT a11, a12, a13; + TYPE_FLOAT a21, a22, a23; + TYPE_FLOAT a31, a32, a33; }; -// Note: the input quaternion need not be of norm 1 for the following function +// Note: the input quaternion need not be of norm 1 for the following function template -R3_matrix quaternion_to_R3_rotation(::boost::quaternion const & q) +R3_matrix quaternion_to_R3_rotation(::boost::math::quaternion const & q) { - using ::std::numeric_limits; - - TYPE_FLOAT a = q.R_component_1(); - TYPE_FLOAT b = q.R_component_2(); - TYPE_FLOAT c = q.R_component_3(); - TYPE_FLOAT d = q.R_component_4(); - - TYPE_FLOAT aa = a*a; - TYPE_FLOAT ab = a*b; - TYPE_FLOAT ac = a*c; - TYPE_FLOAT ad = a*d; - TYPE_FLOAT bb = b*b; - TYPE_FLOAT bc = b*c; - TYPE_FLOAT bd = b*d; - TYPE_FLOAT cc = c*c; - TYPE_FLOAT cd = c*d; - TYPE_FLOAT dd = d*d; - - TYPE_FLOAT norme_carre = aa+bb+cc+dd; - - if (norme_carre <= numeric_limits::epsilon()) - { - ::std::string error_reporting("Argument to quaternion_to_R3_rotation is too small!"); - ::std::underflow_error bad_argument(error_reporting); - - throw(bad_argument); - } - - R3_matrix out_matrix; - - out_matrix.a11 = (aa+bb-cc-dd)/norme_carre; - out_matrix.a12 = 2*(-ad+bc)/norme_carre; - out_matrix.a13 = 2*(ac+bd)/norme_carre; - out_matrix.a21 = 2*(ad+bc)/norme_carre; - out_matrix.a22 = (aa-bb+cc-dd)/norme_carre; - out_matrix.a23 = 2*(-ab+cd)/norme_carre; - out_matrix.a31 = 2*(-ac+bd)/norme_carre; - out_matrix.a32 = 2*(ab+cd)/norme_carre; - out_matrix.a33 = (aa-bb-cc+dd)/norme_carre; - - return(out_matrix); + using ::std::numeric_limits; + + TYPE_FLOAT a = q.R_component_1(); + TYPE_FLOAT b = q.R_component_2(); + TYPE_FLOAT c = q.R_component_3(); + TYPE_FLOAT d = q.R_component_4(); + + TYPE_FLOAT aa = a*a; + TYPE_FLOAT ab = a*b; + TYPE_FLOAT ac = a*c; + TYPE_FLOAT ad = a*d; + TYPE_FLOAT bb = b*b; + TYPE_FLOAT bc = b*c; + TYPE_FLOAT bd = b*d; + TYPE_FLOAT cc = c*c; + TYPE_FLOAT cd = c*d; + TYPE_FLOAT dd = d*d; + + TYPE_FLOAT norme_carre = aa+bb+cc+dd; + + if (norme_carre <= numeric_limits::epsilon()) + { + ::std::string error_reporting("Argument to quaternion_to_R3_rotation is too small!"); + ::std::underflow_error bad_argument(error_reporting); + + throw(bad_argument); + } + + R3_matrix out_matrix; + + out_matrix.a11 = (aa+bb-cc-dd)/norme_carre; + out_matrix.a12 = 2*(-ad+bc)/norme_carre; + out_matrix.a13 = 2*(ac+bd)/norme_carre; + out_matrix.a21 = 2*(ad+bc)/norme_carre; + out_matrix.a22 = (aa-bb+cc-dd)/norme_carre; + out_matrix.a23 = 2*(-ab+cd)/norme_carre; + out_matrix.a31 = 2*(-ac+bd)/norme_carre; + out_matrix.a32 = 2*(ab+cd)/norme_carre; + out_matrix.a33 = (aa-bb-cc+dd)/norme_carre; + + return(out_matrix); } namespace { - template - void find_invariant_vector(R3_matrix const & rot, - TYPE_FLOAT & x, - TYPE_FLOAT & y, - TYPE_FLOAT & z) - { - using ::std::sqrt; - - using ::std::numeric_limits; - - TYPE_FLOAT b11 = rot.a11 - static_cast(1); - TYPE_FLOAT b12 = rot.a12; - TYPE_FLOAT b13 = rot.a13; - TYPE_FLOAT b21 = rot.a21; - TYPE_FLOAT b22 = rot.a22 - static_cast(1); - TYPE_FLOAT b23 = rot.a23; - TYPE_FLOAT b31 = rot.a31; - TYPE_FLOAT b32 = rot.a32; - TYPE_FLOAT b33 = rot.a33 - static_cast(1); - - TYPE_FLOAT minors[9] = - { - b11*b22-b12*b21, - b11*b23-b13*b21, - b12*b23-b13*b22, - b11*b32-b12*b31, - b11*b33-b13*b31, - b12*b33-b13*b32, - b21*b32-b22*b31, - b21*b33-b23*b31, - b22*b33-b23*b32 - }; - - TYPE_FLOAT * where = ::std::max_element(minors, minors+9); - - TYPE_FLOAT det = *where; - - if (det <= numeric_limits::epsilon()) - { - ::std::string error_reporting("Underflow error in find_invariant_vector!"); - ::std::underflow_error processing_error(error_reporting); - - throw(processing_error); - } - - switch (where-minors) - { - case 0: - - z = static_cast(1); - - x = (-b13*b22+b12*b23)/det; - y = (-b11*b23+b13*b21)/det; - - break; - - case 1: - - y = static_cast(1); - - x = (-b12*b23+b13*b22)/det; - z = (-b11*b22+b12*b21)/det; - - break; - - case 2: - - x = static_cast(1); - - y = (-b11*b23+b13*b21)/det; - z = (-b12*b21+b11*b22)/det; - - break; - - case 3: - - z = static_cast(1); - - x = (-b13*b32+b12*b33)/det; - y = (-b11*b33+b13*b31)/det; - - break; - - case 4: - - y = static_cast(1); - - x = (-b12*b33+b13*b32)/det; - z = (-b11*b32+b12*b31)/det; - - break; - - case 5: - - x = static_cast(1); - - y = (-b11*b33+b13*b31)/det; - z = (-b12*b31+b11*b32)/det; - - break; - - case 6: - - z = static_cast(1); - - x = (-b23*b32+b22*b33)/det; - y = (-b21*b33+b23*b31)/det; - - break; - - case 7: - - y = static_cast(1); - - x = (-b22*b33+b23*b32)/det; - z = (-b21*b32+b22*b31)/det; - - break; - - case 8: - - x = static_cast(1); - - y = (-b21*b33+b23*b31)/det; - z = (-b22*b31+b21*b32)/det; - - break; - - default: - - ::std::string error_reporting("Impossible condition in find_invariant_vector"); - ::std::logic_error processing_error(error_reporting); + template + void find_invariant_vector( R3_matrix const & rot, + TYPE_FLOAT & x, + TYPE_FLOAT & y, + TYPE_FLOAT & z) + { + using ::std::sqrt; + + using ::std::numeric_limits; + + TYPE_FLOAT b11 = rot.a11 - static_cast(1); + TYPE_FLOAT b12 = rot.a12; + TYPE_FLOAT b13 = rot.a13; + TYPE_FLOAT b21 = rot.a21; + TYPE_FLOAT b22 = rot.a22 - static_cast(1); + TYPE_FLOAT b23 = rot.a23; + TYPE_FLOAT b31 = rot.a31; + TYPE_FLOAT b32 = rot.a32; + TYPE_FLOAT b33 = rot.a33 - static_cast(1); + + TYPE_FLOAT minors[9] = + { + b11*b22-b12*b21, + b11*b23-b13*b21, + b12*b23-b13*b22, + b11*b32-b12*b31, + b11*b33-b13*b31, + b12*b33-b13*b32, + b21*b32-b22*b31, + b21*b33-b23*b31, + b22*b33-b23*b32 + }; + + TYPE_FLOAT * where = ::std::max_element(minors, minors+9); + + TYPE_FLOAT det = *where; + + if (det <= numeric_limits::epsilon()) + { + ::std::string error_reporting("Underflow error in find_invariant_vector!"); + ::std::underflow_error processing_error(error_reporting); throw(processing_error); - - break; - } - - TYPE_FLOAT vecnorm = sqrt(x*x+y*y+z*z); - - if (vecnorm <= numeric_limits::epsilon()) - { - ::std::string error_reporting("Overflow error in find_invariant_vector!"); - ::std::overflow_error processing_error(error_reporting); - - throw(processing_error); - } - - x /= vecnorm; - y /= vecnorm; - z /= vecnorm; - } - - - template - void find_orthogonal_vector(TYPE_FLOAT x, - TYPE_FLOAT y, - TYPE_FLOAT z, - TYPE_FLOAT & u, - TYPE_FLOAT & v, - TYPE_FLOAT & w) - { - using ::std::abs; - using ::std::sqrt; - - using ::std::numeric_limits; - - TYPE_FLOAT vecnormsqr = x*x+y*y+z*z; - - if (vecnormsqr <= numeric_limits::epsilon()) - { - ::std::string error_reporting("Underflow error in find_orthogonal_vector!"); - ::std::underflow_error processing_error(error_reporting); - - throw(processing_error); - } - - TYPE_FLOAT lambda; - - TYPE_FLOAT components[3] = - { - abs(x), - abs(y), - abs(z) - }; - - TYPE_FLOAT * where = ::std::min_element(components, components+3); - - switch (where-components) - { - case 0: - - if (*where <= numeric_limits::epsilon()) - { - v = - w = static_cast(0); - u = static_cast(1); - } - else - { - lambda = -x/vecnormsqr; - - u = static_cast(1) + lambda*x; - v = lambda*y; - w = lambda*z; - } - - break; - - case 1: - - if (*where <= numeric_limits::epsilon()) - { - u = - w = static_cast(0); - v = static_cast(1); - } - else - { - lambda = -y/vecnormsqr; - - u = lambda*x; - v = static_cast(1) + lambda*y; - w = lambda*z; - } - - break; - - case 2: - - if (*where <= numeric_limits::epsilon()) - { - u = - v = static_cast(0); - w = static_cast(1); - } - else - { - lambda = -z/vecnormsqr; - - u = lambda*x; - v = lambda*y; - w = static_cast(1) + lambda*z; - } - - break; - - default: - - ::std::string error_reporting("Impossible condition in find_invariant_vector"); - ::std::logic_error processing_error(error_reporting); + } + + switch (where-minors) + { + case 0: + + z = static_cast(1); + + x = (-b13*b22+b12*b23)/det; + y = (-b11*b23+b13*b21)/det; + + break; + + case 1: + + y = static_cast(1); + + x = (-b12*b23+b13*b22)/det; + z = (-b11*b22+b12*b21)/det; + + break; + + case 2: + + x = static_cast(1); + + y = (-b11*b23+b13*b21)/det; + z = (-b12*b21+b11*b22)/det; + + break; + + case 3: + + z = static_cast(1); + + x = (-b13*b32+b12*b33)/det; + y = (-b11*b33+b13*b31)/det; + + break; + + case 4: + + y = static_cast(1); + + x = (-b12*b33+b13*b32)/det; + z = (-b11*b32+b12*b31)/det; + + break; + + case 5: + + x = static_cast(1); + + y = (-b11*b33+b13*b31)/det; + z = (-b12*b31+b11*b32)/det; + + break; + + case 6: + + z = static_cast(1); + + x = (-b23*b32+b22*b33)/det; + y = (-b21*b33+b23*b31)/det; + + break; + + case 7: + + y = static_cast(1); + + x = (-b22*b33+b23*b32)/det; + z = (-b21*b32+b22*b31)/det; + + break; + + case 8: + + x = static_cast(1); + + y = (-b21*b33+b23*b31)/det; + z = (-b22*b31+b21*b32)/det; + + break; + + default: + + ::std::string error_reporting("Impossible condition in find_invariant_vector"); + ::std::logic_error processing_error(error_reporting); + + throw(processing_error); + + break; + } + + TYPE_FLOAT vecnorm = sqrt(x*x+y*y+z*z); + + if (vecnorm <= numeric_limits::epsilon()) + { + ::std::string error_reporting("Overflow error in find_invariant_vector!"); + ::std::overflow_error processing_error(error_reporting); throw(processing_error); + } + + x /= vecnorm; + y /= vecnorm; + z /= vecnorm; + } + + + template + void find_orthogonal_vector( TYPE_FLOAT x, + TYPE_FLOAT y, + TYPE_FLOAT z, + TYPE_FLOAT & u, + TYPE_FLOAT & v, + TYPE_FLOAT & w) + { + using ::std::abs; + using ::std::sqrt; + + using ::std::numeric_limits; + + TYPE_FLOAT vecnormsqr = x*x+y*y+z*z; + + if (vecnormsqr <= numeric_limits::epsilon()) + { + ::std::string error_reporting("Underflow error in find_orthogonal_vector!"); + ::std::underflow_error processing_error(error_reporting); - break; - } - - TYPE_FLOAT vecnorm = sqrt(u*u+v*v+w*w); - - if (vecnorm <= numeric_limits::epsilon()) - { - ::std::string error_reporting("Underflow error in find_orthogonal_vector!"); - ::std::underflow_error processing_error(error_reporting); - - throw(processing_error); - } - - u /= vecnorm; - v /= vecnorm; - w /= vecnorm; - } - - - // Note: we want [[v, v, w], [r, s, t], [x, y, z]] to be a direct orthogonal basis - // of R^3. It might not be orthonormal, however, and we do not check if the - // two input vectors are colinear or not. - - template - void find_vector_for_BOD(TYPE_FLOAT x, - TYPE_FLOAT y, - TYPE_FLOAT z, - TYPE_FLOAT u, - TYPE_FLOAT v, - TYPE_FLOAT w, - TYPE_FLOAT & r, - TYPE_FLOAT & s, - TYPE_FLOAT & t) - { - r = +y*w-z*v; - s = -x*w+z*u; - t = +x*v-y*u; - } + throw(processing_error); + } + + TYPE_FLOAT lambda; + + TYPE_FLOAT components[3] = + { + abs(x), + abs(y), + abs(z) + }; + + TYPE_FLOAT * where = ::std::min_element(components, components+3); + + switch (where-components) + { + case 0: + + if (*where <= numeric_limits::epsilon()) + { + v = + w = static_cast(0); + u = static_cast(1); + } + else + { + lambda = -x/vecnormsqr; + + u = static_cast(1) + lambda*x; + v = lambda*y; + w = lambda*z; + } + + break; + + case 1: + + if (*where <= numeric_limits::epsilon()) + { + u = + w = static_cast(0); + v = static_cast(1); + } + else + { + lambda = -y/vecnormsqr; + + u = lambda*x; + v = static_cast(1) + lambda*y; + w = lambda*z; + } + + break; + + case 2: + + if (*where <= numeric_limits::epsilon()) + { + u = + v = static_cast(0); + w = static_cast(1); + } + else + { + lambda = -z/vecnormsqr; + + u = lambda*x; + v = lambda*y; + w = static_cast(1) + lambda*z; + } + + break; + + default: + + ::std::string error_reporting("Impossible condition in find_invariant_vector"); + ::std::logic_error processing_error(error_reporting); + + throw(processing_error); + + break; + } + + TYPE_FLOAT vecnorm = sqrt(u*u+v*v+w*w); + + if (vecnorm <= numeric_limits::epsilon()) + { + ::std::string error_reporting("Underflow error in find_orthogonal_vector!"); + ::std::underflow_error processing_error(error_reporting); + + throw(processing_error); + } + + u /= vecnorm; + v /= vecnorm; + w /= vecnorm; + } + + + // Note: we want [[v, v, w], [r, s, t], [x, y, z]] to be a direct orthogonal basis + // of R^3. It might not be orthonormal, however, and we do not check if the + // two input vectors are colinear or not. + + template + void find_vector_for_BOD(TYPE_FLOAT x, + TYPE_FLOAT y, + TYPE_FLOAT z, + TYPE_FLOAT u, + TYPE_FLOAT v, + TYPE_FLOAT w, + TYPE_FLOAT & r, + TYPE_FLOAT & s, + TYPE_FLOAT & t) + { + r = +y*w-z*v; + s = -x*w+z*u; + t = +x*v-y*u; + } } template -inline bool is_R3_rotation_matrix(R3_matrix const & mat) +inline bool is_R3_rotation_matrix(R3_matrix const & mat) { - using ::std::abs; - - using ::std::numeric_limits; - - return ( - !( - (abs(mat.a11*mat.a11+mat.a21*mat.a21+mat.a31*mat.a31 - static_cast(1)) > static_cast(10)*numeric_limits::epsilon())|| - (abs(mat.a11*mat.a12+mat.a21*mat.a22+mat.a31*mat.a32 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| - (abs(mat.a11*mat.a13+mat.a21*mat.a23+mat.a31*mat.a33 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| - //(abs(mat.a11*mat.a12+mat.a21*mat.a22+mat.a31*mat.a32 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| - (abs(mat.a12*mat.a12+mat.a22*mat.a22+mat.a32*mat.a32 - static_cast(1)) > static_cast(10)*numeric_limits::epsilon())|| - (abs(mat.a12*mat.a13+mat.a22*mat.a23+mat.a32*mat.a33 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| - //(abs(mat.a11*mat.a13+mat.a21*mat.a23+mat.a31*mat.a33 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| - //(abs(mat.a12*mat.a13+mat.a22*mat.a23+mat.a32*mat.a33 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| - (abs(mat.a13*mat.a13+mat.a23*mat.a23+mat.a33*mat.a33 - static_cast(1)) > static_cast(10)*numeric_limits::epsilon()) - ) - ); + using ::std::abs; + + using ::std::numeric_limits; + + return ( + !( + (abs(mat.a11*mat.a11+mat.a21*mat.a21+mat.a31*mat.a31 - static_cast(1)) > static_cast(10)*numeric_limits::epsilon())|| + (abs(mat.a11*mat.a12+mat.a21*mat.a22+mat.a31*mat.a32 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| + (abs(mat.a11*mat.a13+mat.a21*mat.a23+mat.a31*mat.a33 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| + //(abs(mat.a11*mat.a12+mat.a21*mat.a22+mat.a31*mat.a32 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| + (abs(mat.a12*mat.a12+mat.a22*mat.a22+mat.a32*mat.a32 - static_cast(1)) > static_cast(10)*numeric_limits::epsilon())|| + (abs(mat.a12*mat.a13+mat.a22*mat.a23+mat.a32*mat.a33 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| + //(abs(mat.a11*mat.a13+mat.a21*mat.a23+mat.a31*mat.a33 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| + //(abs(mat.a12*mat.a13+mat.a22*mat.a23+mat.a32*mat.a33 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| + (abs(mat.a13*mat.a13+mat.a23*mat.a23+mat.a33*mat.a33 - static_cast(1)) > static_cast(10)*numeric_limits::epsilon()) + ) + ); } template -::boost::quaternion R3_rotation_to_quaternion(R3_matrix const & rot, - ::boost::quaternion const * hint = 0) +::boost::math::quaternion R3_rotation_to_quaternion( R3_matrix const & rot, + ::boost::math::quaternion const * hint = 0) { - using ::std::abs; - using ::std::sqrt; - - using ::std::numeric_limits; - - using ::boost::abs; - - if (!is_R3_rotation_matrix(rot)) - { - ::std::string error_reporting("Argument to R3_rotation_to_quaternion is not an R^3 rotation matrix!"); - ::std::range_error bad_argument(error_reporting); - - throw(bad_argument); - } - - ::boost::quaternion q; - - if ( - (abs(rot.a11 - static_cast(1)) <= numeric_limits::epsilon())&& - (abs(rot.a22 - static_cast(1)) <= numeric_limits::epsilon())&& - (abs(rot.a33 - static_cast(1)) <= numeric_limits::epsilon()) - ) - { - q = ::boost::quaternion(1); - } - else - { - TYPE_FLOAT cos_theta = (rot.a11+rot.a22+rot.a33-static_cast(1))/static_cast(2); - TYPE_FLOAT stuff = (cos_theta+static_cast(1))/static_cast(2); - TYPE_FLOAT cos_theta_sur_2 = sqrt(stuff); - TYPE_FLOAT sin_theta_sur_2 = sqrt(1-stuff); - - TYPE_FLOAT x; - TYPE_FLOAT y; - TYPE_FLOAT z; - - find_invariant_vector(rot, x, y, z); - - TYPE_FLOAT u; - TYPE_FLOAT v; - TYPE_FLOAT w; - - find_orthogonal_vector(x, y, z, u, v, w); - - TYPE_FLOAT r; - TYPE_FLOAT s; - TYPE_FLOAT t; - - find_vector_for_BOD(x, y, z, u, v, w, r, s, t); - - TYPE_FLOAT ru = rot.a11*u+rot.a12*v+rot.a13*w; - TYPE_FLOAT rv = rot.a21*u+rot.a22*v+rot.a23*w; - TYPE_FLOAT rw = rot.a31*u+rot.a32*v+rot.a33*w; - - TYPE_FLOAT angle_sign_determinator = r*ru+s*rv+t*rw; - - if (angle_sign_determinator > +numeric_limits::epsilon()) - { - q = ::boost::quaternion(cos_theta_sur_2, +x*sin_theta_sur_2, +y*sin_theta_sur_2, +z*sin_theta_sur_2); - } - else if (angle_sign_determinator < -numeric_limits::epsilon()) - { - q = ::boost::quaternion(cos_theta_sur_2, -x*sin_theta_sur_2, -y*sin_theta_sur_2, -z*sin_theta_sur_2); - } - else - { - TYPE_FLOAT desambiguator = u*ru+v*rv+w*rw; - - if (desambiguator >= static_cast(1)) - { - q = ::boost::quaternion(0, +x, +y, +z); - } - else - { - q = ::boost::quaternion(0, -x, -y, -z); - } - } - } - - if ((hint != 0) && (abs(*hint+q) < abs(*hint-q))) - { - return(-q); - } - - return(q); + using ::std::abs; + using ::std::sqrt; + + using ::std::numeric_limits; + + using ::boost::math::abs; + + if (!is_R3_rotation_matrix(rot)) + { + ::std::string error_reporting("Argument to R3_rotation_to_quaternion is not an R^3 rotation matrix!"); + ::std::range_error bad_argument(error_reporting); + + throw(bad_argument); + } + + ::boost::math::quaternion q; + + if ( + (abs(rot.a11 - static_cast(1)) <= numeric_limits::epsilon())&& + (abs(rot.a22 - static_cast(1)) <= numeric_limits::epsilon())&& + (abs(rot.a33 - static_cast(1)) <= numeric_limits::epsilon()) + ) + { + q = ::boost::math::quaternion(1); + } + else + { + TYPE_FLOAT cos_theta = (rot.a11+rot.a22+rot.a33-static_cast(1))/static_cast(2); + TYPE_FLOAT stuff = (cos_theta+static_cast(1))/static_cast(2); + TYPE_FLOAT cos_theta_sur_2 = sqrt(stuff); + TYPE_FLOAT sin_theta_sur_2 = sqrt(1-stuff); + + TYPE_FLOAT x; + TYPE_FLOAT y; + TYPE_FLOAT z; + + find_invariant_vector(rot, x, y, z); + + TYPE_FLOAT u; + TYPE_FLOAT v; + TYPE_FLOAT w; + + find_orthogonal_vector(x, y, z, u, v, w); + + TYPE_FLOAT r; + TYPE_FLOAT s; + TYPE_FLOAT t; + + find_vector_for_BOD(x, y, z, u, v, w, r, s, t); + + TYPE_FLOAT ru = rot.a11*u+rot.a12*v+rot.a13*w; + TYPE_FLOAT rv = rot.a21*u+rot.a22*v+rot.a23*w; + TYPE_FLOAT rw = rot.a31*u+rot.a32*v+rot.a33*w; + + TYPE_FLOAT angle_sign_determinator = r*ru+s*rv+t*rw; + + if (angle_sign_determinator > +numeric_limits::epsilon()) + { + q = ::boost::math::quaternion(cos_theta_sur_2, +x*sin_theta_sur_2, +y*sin_theta_sur_2, +z*sin_theta_sur_2); + } + else if (angle_sign_determinator < -numeric_limits::epsilon()) + { + q = ::boost::math::quaternion(cos_theta_sur_2, -x*sin_theta_sur_2, -y*sin_theta_sur_2, -z*sin_theta_sur_2); + } + else + { + TYPE_FLOAT desambiguator = u*ru+v*rv+w*rw; + + if (desambiguator >= static_cast(1)) + { + q = ::boost::math::quaternion(0, +x, +y, +z); + } + else + { + q = ::boost::math::quaternion(0, -x, -y, -z); + } + } + } + + if ((hint != 0) && (abs(*hint+q) < abs(*hint-q))) + { + return(-q); + } + + return(q); } #endif /* TEST_HSO3_HPP */ diff --git a/quaternion/HSO3SO4.cpp b/quaternion/HSO3SO4.cpp index ae19b46c6..960f33460 100644 --- a/quaternion/HSO3SO4.cpp +++ b/quaternion/HSO3SO4.cpp @@ -14,431 +14,431 @@ #include "HSO4.hpp" -const int number_of_intervals = 5; +const int number_of_intervals = 5; -const float pi = ::std::atan(1.0f)*4; +const float pi = ::std::atan(1.0f)*4; -void test_SO3(); - -void test_SO4(); +void test_SO3(); + +void test_SO4(); -int main() +int main() { - test_SO3(); - - test_SO4(); - - ::std::cout << "That's all folks!" << ::std::endl; + test_SO3(); + + test_SO4(); + + ::std::cout << "That's all folks!" << ::std::endl; } // -// Test of quaternion and R^3 rotation relationship +// Test of quaternion and R^3 rotation relationship // -void test_SO3_spherical() +void test_SO3_spherical() { - ::std::cout << "Testing spherical:" << ::std::endl; - ::std::cout << ::std::endl; - - const float rho = 1.0f; - - float theta; - float phi1; - float phi2; - - for (int idxphi2 = 0; idxphi2 <= number_of_intervals; idxphi2++) - { - phi2 = (-pi/2)+(idxphi2*pi)/number_of_intervals; - - for (int idxphi1 = 0; idxphi1 <= number_of_intervals; idxphi1++) - { - phi1 = (-pi/2)+(idxphi1*pi)/number_of_intervals; - - for (int idxtheta = 0; idxtheta <= number_of_intervals; idxtheta++) - { - theta = -pi+(idxtheta*(2*pi))/number_of_intervals; + ::std::cout << "Testing spherical:" << ::std::endl; + ::std::cout << ::std::endl; + + const float rho = 1.0f; + + float theta; + float phi1; + float phi2; + + for (int idxphi2 = 0; idxphi2 <= number_of_intervals; idxphi2++) + { + phi2 = (-pi/2)+(idxphi2*pi)/number_of_intervals; + + for (int idxphi1 = 0; idxphi1 <= number_of_intervals; idxphi1++) + { + phi1 = (-pi/2)+(idxphi1*pi)/number_of_intervals; - //::std::cout << "theta = " << theta << " ; "; - //::std::cout << "phi1 = " << phi1 << " ; "; - //::std::cout << "phi2 = " << phi2; - //::std::cout << ::std::endl; - - ::boost::quaternion q = ::boost::spherical(rho, theta, phi1, phi2); - - //::std::cout << "q = " << q << ::std::endl; - - R3_matrix rot = quaternion_to_R3_rotation(q); - - //::std::cout << "rot = "; - //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << ::std::endl; - //::std::cout << "\t"; - //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << ::std::endl; - //::std::cout << "\t"; - //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << ::std::endl; - - ::boost::quaternion p = R3_rotation_to_quaternion(rot, &q); - - //::std::cout << "p = " << p << ::std::endl; - - //::std::cout << "round trip discrepancy: " << ::boost::abs(q-p) << ::std::endl; - - //::std::cout << ::std::endl; - } - } - } - - ::std::cout << ::std::endl; -} - - -void test_SO3_semipolar() -{ - ::std::cout << "Testing semipolar:" << ::std::endl; - ::std::cout << ::std::endl; - - const float rho = 1.0f; - - float alpha; - float theta1; - float theta2; - - for (int idxalpha = 0; idxalpha <= number_of_intervals; idxalpha++) - { - alpha = (idxalpha*(pi/2))/number_of_intervals; - - for (int idxtheta1 = 0; idxtheta1 <= number_of_intervals; idxtheta1++) - { - theta1 = -pi+(idxtheta1*(2*pi))/number_of_intervals; - - for (int idxtheta2 = 0; idxtheta2 <= number_of_intervals; idxtheta2++) - { - theta2 = -pi+(idxtheta2*(2*pi))/number_of_intervals; - - //::std::cout << "alpha = " << alpha << " ; "; - //::std::cout << "theta1 = " << theta1 << " ; "; - //::std::cout << "theta2 = " << theta2; - //::std::cout << ::std::endl; - - ::boost::quaternion q = ::boost::semipolar(rho, alpha, theta1, theta2); - - //::std::cout << "q = " << q << ::std::endl; - - R3_matrix rot = quaternion_to_R3_rotation(q); - - //::std::cout << "rot = "; - //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << ::std::endl; - //::std::cout << "\t"; - //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << ::std::endl; - //::std::cout << "\t"; - //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << ::std::endl; - - ::boost::quaternion p = R3_rotation_to_quaternion(rot, &q); - - //::std::cout << "p = " << p << ::std::endl; - - //::std::cout << "round trip discrepancy: " << ::boost::abs(q-p) << ::std::endl; - - //::std::cout << ::std::endl; - } - } - } - - ::std::cout << ::std::endl; -} - - -void test_SO3_multipolar() -{ - ::std::cout << "Testing multipolar:" << ::std::endl; - ::std::cout << ::std::endl; - - float rho1; - float rho2; - - float theta1; - float theta2; - - for (int idxrho = 0; idxrho <= number_of_intervals; idxrho++) - { - rho1 = (idxrho*1.0f)/number_of_intervals; - rho2 = ::std::sqrt(1.0f-rho1*rho1); - - for (int idxtheta1 = 0; idxtheta1 <= number_of_intervals; idxtheta1++) - { - theta1 = -pi+(idxtheta1*(2*pi))/number_of_intervals; - - for (int idxtheta2 = 0; idxtheta2 <= number_of_intervals; idxtheta2++) - { - theta2 = -pi+(idxtheta2*(2*pi))/number_of_intervals; - - //::std::cout << "rho1 = " << rho1 << " ; "; - //::std::cout << "theta1 = " << theta1 << " ; "; - //::std::cout << "theta2 = " << theta2; - //::std::cout << ::std::endl; - - ::boost::quaternion q = ::boost::multipolar(rho1, theta1, rho2, theta2); - - //::std::cout << "q = " << q << ::std::endl; - - R3_matrix rot = quaternion_to_R3_rotation(q); - - //::std::cout << "rot = "; - //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << ::std::endl; - //::std::cout << "\t"; - //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << ::std::endl; - //::std::cout << "\t"; - //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << ::std::endl; - - ::boost::quaternion p = R3_rotation_to_quaternion(rot, &q); - - //::std::cout << "p = " << p << ::std::endl; - - //::std::cout << "round trip discrepancy: " << ::boost::abs(q-p) << ::std::endl; - - //::std::cout << ::std::endl; - } - } - } - - ::std::cout << ::std::endl; -} - - -void test_SO3_cylindrospherical() -{ - ::std::cout << "Testing cylindrospherical:" << ::std::endl; - ::std::cout << ::std::endl; - - float t; - - float radius; - float longitude; - float latitude; - - for (int idxt = 0; idxt <= number_of_intervals; idxt++) - { - t = -1.0f+(idxt*2.0f)/number_of_intervals; - radius = ::std::sqrt(1.0f-t*t); - - for (int idxlatitude = 0; idxlatitude <= number_of_intervals; idxlatitude++) - { - latitude = (-pi/2)+(idxlatitude*pi)/number_of_intervals; - - for (int idxlongitude = 0; idxlongitude <= number_of_intervals; idxlongitude++) - { - longitude = -pi+(idxlongitude*(2*pi))/number_of_intervals; - - //::std::cout << "t = " << t << " ; "; - //::std::cout << "longitude = " << longitude; - //::std::cout << "latitude = " << latitude; - //::std::cout << ::std::endl; - - ::boost::quaternion q = ::boost::cylindrospherical(t, radius, longitude, latitude); - - //::std::cout << "q = " << q << ::std::endl; - - R3_matrix rot = quaternion_to_R3_rotation(q); - - //::std::cout << "rot = "; - //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << ::std::endl; - //::std::cout << "\t"; - //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << ::std::endl; - //::std::cout << "\t"; - //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << ::std::endl; - - ::boost::quaternion p = R3_rotation_to_quaternion(rot, &q); - - //::std::cout << "p = " << p << ::std::endl; - - //::std::cout << "round trip discrepancy: " << ::boost::abs(q-p) << ::std::endl; - - //::std::cout << ::std::endl; - } - } - } - - ::std::cout << ::std::endl; -} - - -void test_SO3_cylindrical() -{ - ::std::cout << "Testing cylindrical:" << ::std::endl; - ::std::cout << ::std::endl; - - float r; - float angle; - - float h1; - float h2; - - for (int idxh2 = 0; idxh2 <= number_of_intervals; idxh2++) - { - h2 = -1.0f+(idxh2*2.0f)/number_of_intervals; - - for (int idxh1 = 0; idxh1 <= number_of_intervals; idxh1++) - { - h1 = ::std::sqrt(1.0f-h2*h2)*(-1.0f+(idxh2*2.0f)/number_of_intervals); - r = ::std::sqrt(1.0f-h1*h1-h2*h2); - - for (int idxangle = 0; idxangle <= number_of_intervals; idxangle++) - { - angle = -pi+(idxangle*(2*pi))/number_of_intervals; - - //::std::cout << "angle = " << angle << " ; "; - //::std::cout << "h1 = " << h1; - //::std::cout << "h2 = " << h2; - //::std::cout << ::std::endl; - - ::boost::quaternion q = ::boost::cylindrical(r, angle, h1, h2); - - //::std::cout << "q = " << q << ::std::endl; - - R3_matrix rot = quaternion_to_R3_rotation(q); - - //::std::cout << "rot = "; - //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << ::std::endl; - //::std::cout << "\t"; - //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << ::std::endl; - //::std::cout << "\t"; - //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << ::std::endl; - - ::boost::quaternion p = R3_rotation_to_quaternion(rot, &q); - - //::std::cout << "p = " << p << ::std::endl; - - //::std::cout << "round trip discrepancy: " << ::boost::abs(q-p) << ::std::endl; - - //::std::cout << ::std::endl; - } - } - } - - ::std::cout << ::std::endl; -} - - -void test_SO3() -{ - ::std::cout << "Testing SO3:" << ::std::endl; - ::std::cout << ::std::endl; - - test_SO3_spherical(); - - test_SO3_semipolar(); - - test_SO3_multipolar(); - - test_SO3_cylindrospherical(); - - test_SO3_cylindrical(); -} - - -// -// Test of quaternion and R^4 rotation relationship -// - -void test_SO4_spherical() -{ - ::std::cout << "Testing spherical:" << ::std::endl; - ::std::cout << ::std::endl; - - const float rho1 = 1.0f; - const float rho2 = 1.0f; - - float theta1; - float phi11; - float phi21; - - float theta2; - float phi12; - float phi22; - - for (int idxphi21 = 0; idxphi21 <= number_of_intervals; idxphi21++) - { - phi21 = (-pi/2)+(idxphi21*pi)/number_of_intervals; - - for (int idxphi22 = 0; idxphi22 <= number_of_intervals; idxphi22++) - { - phi22 = (-pi/2)+(idxphi22*pi)/number_of_intervals; - - for (int idxphi11 = 0; idxphi11 <= number_of_intervals; idxphi11++) - { - phi11 = (-pi/2)+(idxphi11*pi)/number_of_intervals; - - for (int idxphi12 = 0; idxphi12 <= number_of_intervals; idxphi12++) + for (int idxtheta = 0; idxtheta <= number_of_intervals; idxtheta++) { - phi12 = (-pi/2)+(idxphi12*pi)/number_of_intervals; - - for (int idxtheta1 = 0; idxtheta1 <= number_of_intervals; idxtheta1++) - { - theta1 = -pi+(idxtheta1*(2*pi))/number_of_intervals; - - for (int idxtheta2 = 0; idxtheta2 <= number_of_intervals; idxtheta2++) - { - theta2 = -pi+(idxtheta2*(2*pi))/number_of_intervals; - - //::std::cout << "theta1 = " << theta1 << " ; "; - //::std::cout << "phi11 = " << phi11 << " ; "; - //::std::cout << "phi21 = " << phi21; - //::std::cout << "theta2 = " << theta2 << " ; "; - //::std::cout << "phi12 = " << phi12 << " ; "; - //::std::cout << "phi22 = " << phi22; - //::std::cout << ::std::endl; - - ::boost::quaternion p1 = ::boost::spherical(rho1, theta1, phi11, phi21); - - //::std::cout << "p1 = " << p1 << ::std::endl; - - ::boost::quaternion q1 = ::boost::spherical(rho2, theta2, phi12, phi22); - - //::std::cout << "q1 = " << q1 << ::std::endl; - - ::std::pair< ::boost::quaternion , ::boost::quaternion > pq1 = - ::std::make_pair(p1,q1); - - R4_matrix rot = quaternions_to_R4_rotation(pq1); - - //::std::cout << "rot = "; - //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << "\t" << rot.a14 << ::std::endl; - //::std::cout << "\t"; - //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << "\t" << rot.a24 << ::std::endl; - //::std::cout << "\t"; - //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << "\t" << rot.a34 << ::std::endl; - //::std::cout << "\t"; - //::std::cout << "\t" << rot.a41 << "\t" << rot.a42 << "\t" << rot.a43 << "\t" << rot.a44 << ::std::endl; - - ::std::pair< ::boost::quaternion , ::boost::quaternion > pq2 = - R4_rotation_to_quaternions(rot, &pq1); - - //::std::cout << "p1 = " << pq.first << ::std::endl; - //::std::cout << "p2 = " << pq.second << ::std::endl; - - //::std::cout << "round trip discrepancy: " << ::std::sqrt(::boost::norm(pq1.first-pq2.first)+::boost::norm(pq1.second-pq2.second)) << ::std::endl; - - //::std::cout << ::std::endl; - } - } + theta = -pi+(idxtheta*(2*pi))/number_of_intervals; + + //::std::cout << "theta = " << theta << " ; "; + //::std::cout << "phi1 = " << phi1 << " ; "; + //::std::cout << "phi2 = " << phi2; + //::std::cout << ::std::endl; + + ::boost::math::quaternion q = ::boost::math::spherical(rho, theta, phi1, phi2); + + //::std::cout << "q = " << q << ::std::endl; + + R3_matrix rot = quaternion_to_R3_rotation(q); + + //::std::cout << "rot = "; + //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << ::std::endl; + //::std::cout << "\t"; + //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << ::std::endl; + //::std::cout << "\t"; + //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << ::std::endl; + + ::boost::math::quaternion p = R3_rotation_to_quaternion(rot, &q); + + //::std::cout << "p = " << p << ::std::endl; + + //::std::cout << "round trip discrepancy: " << ::boost::math::abs(q-p) << ::std::endl; + + //::std::cout << ::std::endl; } - } - } - } - - ::std::cout << ::std::endl; + } + } + + ::std::cout << ::std::endl; } - -void test_SO4() + +void test_SO3_semipolar() { - ::std::cout << "Testing SO4:" << ::std::endl; - ::std::cout << ::std::endl; - - test_SO4_spherical(); + ::std::cout << "Testing semipolar:" << ::std::endl; + ::std::cout << ::std::endl; + + const float rho = 1.0f; + + float alpha; + float theta1; + float theta2; + + for (int idxalpha = 0; idxalpha <= number_of_intervals; idxalpha++) + { + alpha = (idxalpha*(pi/2))/number_of_intervals; + + for (int idxtheta1 = 0; idxtheta1 <= number_of_intervals; idxtheta1++) + { + theta1 = -pi+(idxtheta1*(2*pi))/number_of_intervals; + + for (int idxtheta2 = 0; idxtheta2 <= number_of_intervals; idxtheta2++) + { + theta2 = -pi+(idxtheta2*(2*pi))/number_of_intervals; + + //::std::cout << "alpha = " << alpha << " ; "; + //::std::cout << "theta1 = " << theta1 << " ; "; + //::std::cout << "theta2 = " << theta2; + //::std::cout << ::std::endl; + + ::boost::math::quaternion q = ::boost::math::semipolar(rho, alpha, theta1, theta2); + + //::std::cout << "q = " << q << ::std::endl; + + R3_matrix rot = quaternion_to_R3_rotation(q); + + //::std::cout << "rot = "; + //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << ::std::endl; + //::std::cout << "\t"; + //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << ::std::endl; + //::std::cout << "\t"; + //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << ::std::endl; + + ::boost::math::quaternion p = R3_rotation_to_quaternion(rot, &q); + + //::std::cout << "p = " << p << ::std::endl; + + //::std::cout << "round trip discrepancy: " << ::boost::math::abs(q-p) << ::std::endl; + + //::std::cout << ::std::endl; + } + } + } + + ::std::cout << ::std::endl; +} + + +void test_SO3_multipolar() +{ + ::std::cout << "Testing multipolar:" << ::std::endl; + ::std::cout << ::std::endl; + + float rho1; + float rho2; + + float theta1; + float theta2; + + for (int idxrho = 0; idxrho <= number_of_intervals; idxrho++) + { + rho1 = (idxrho*1.0f)/number_of_intervals; + rho2 = ::std::sqrt(1.0f-rho1*rho1); + + for (int idxtheta1 = 0; idxtheta1 <= number_of_intervals; idxtheta1++) + { + theta1 = -pi+(idxtheta1*(2*pi))/number_of_intervals; + + for (int idxtheta2 = 0; idxtheta2 <= number_of_intervals; idxtheta2++) + { + theta2 = -pi+(idxtheta2*(2*pi))/number_of_intervals; + + //::std::cout << "rho1 = " << rho1 << " ; "; + //::std::cout << "theta1 = " << theta1 << " ; "; + //::std::cout << "theta2 = " << theta2; + //::std::cout << ::std::endl; + + ::boost::math::quaternion q = ::boost::math::multipolar(rho1, theta1, rho2, theta2); + + //::std::cout << "q = " << q << ::std::endl; + + R3_matrix rot = quaternion_to_R3_rotation(q); + + //::std::cout << "rot = "; + //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << ::std::endl; + //::std::cout << "\t"; + //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << ::std::endl; + //::std::cout << "\t"; + //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << ::std::endl; + + ::boost::math::quaternion p = R3_rotation_to_quaternion(rot, &q); + + //::std::cout << "p = " << p << ::std::endl; + + //::std::cout << "round trip discrepancy: " << ::boost::math::abs(q-p) << ::std::endl; + + //::std::cout << ::std::endl; + } + } + } + + ::std::cout << ::std::endl; +} + + +void test_SO3_cylindrospherical() +{ + ::std::cout << "Testing cylindrospherical:" << ::std::endl; + ::std::cout << ::std::endl; + + float t; + + float radius; + float longitude; + float latitude; + + for (int idxt = 0; idxt <= number_of_intervals; idxt++) + { + t = -1.0f+(idxt*2.0f)/number_of_intervals; + radius = ::std::sqrt(1.0f-t*t); + + for (int idxlatitude = 0; idxlatitude <= number_of_intervals; idxlatitude++) + { + latitude = (-pi/2)+(idxlatitude*pi)/number_of_intervals; + + for (int idxlongitude = 0; idxlongitude <= number_of_intervals; idxlongitude++) + { + longitude = -pi+(idxlongitude*(2*pi))/number_of_intervals; + + //::std::cout << "t = " << t << " ; "; + //::std::cout << "longitude = " << longitude; + //::std::cout << "latitude = " << latitude; + //::std::cout << ::std::endl; + + ::boost::math::quaternion q = ::boost::math::cylindrospherical(t, radius, longitude, latitude); + + //::std::cout << "q = " << q << ::std::endl; + + R3_matrix rot = quaternion_to_R3_rotation(q); + + //::std::cout << "rot = "; + //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << ::std::endl; + //::std::cout << "\t"; + //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << ::std::endl; + //::std::cout << "\t"; + //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << ::std::endl; + + ::boost::math::quaternion p = R3_rotation_to_quaternion(rot, &q); + + //::std::cout << "p = " << p << ::std::endl; + + //::std::cout << "round trip discrepancy: " << ::boost::math::abs(q-p) << ::std::endl; + + //::std::cout << ::std::endl; + } + } + } + + ::std::cout << ::std::endl; +} + + +void test_SO3_cylindrical() +{ + ::std::cout << "Testing cylindrical:" << ::std::endl; + ::std::cout << ::std::endl; + + float r; + float angle; + + float h1; + float h2; + + for (int idxh2 = 0; idxh2 <= number_of_intervals; idxh2++) + { + h2 = -1.0f+(idxh2*2.0f)/number_of_intervals; + + for (int idxh1 = 0; idxh1 <= number_of_intervals; idxh1++) + { + h1 = ::std::sqrt(1.0f-h2*h2)*(-1.0f+(idxh2*2.0f)/number_of_intervals); + r = ::std::sqrt(1.0f-h1*h1-h2*h2); + + for (int idxangle = 0; idxangle <= number_of_intervals; idxangle++) + { + angle = -pi+(idxangle*(2*pi))/number_of_intervals; + + //::std::cout << "angle = " << angle << " ; "; + //::std::cout << "h1 = " << h1; + //::std::cout << "h2 = " << h2; + //::std::cout << ::std::endl; + + ::boost::math::quaternion q = ::boost::math::cylindrical(r, angle, h1, h2); + + //::std::cout << "q = " << q << ::std::endl; + + R3_matrix rot = quaternion_to_R3_rotation(q); + + //::std::cout << "rot = "; + //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << ::std::endl; + //::std::cout << "\t"; + //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << ::std::endl; + //::std::cout << "\t"; + //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << ::std::endl; + + ::boost::math::quaternion p = R3_rotation_to_quaternion(rot, &q); + + //::std::cout << "p = " << p << ::std::endl; + + //::std::cout << "round trip discrepancy: " << ::boost::math::abs(q-p) << ::std::endl; + + //::std::cout << ::std::endl; + } + } + } + + ::std::cout << ::std::endl; +} + + +void test_SO3() +{ + ::std::cout << "Testing SO3:" << ::std::endl; + ::std::cout << ::std::endl; + + test_SO3_spherical(); + + test_SO3_semipolar(); + + test_SO3_multipolar(); + + test_SO3_cylindrospherical(); + + test_SO3_cylindrical(); +} + + +// +// Test of quaternion and R^4 rotation relationship +// + +void test_SO4_spherical() +{ + ::std::cout << "Testing spherical:" << ::std::endl; + ::std::cout << ::std::endl; + + const float rho1 = 1.0f; + const float rho2 = 1.0f; + + float theta1; + float phi11; + float phi21; + + float theta2; + float phi12; + float phi22; + + for (int idxphi21 = 0; idxphi21 <= number_of_intervals; idxphi21++) + { + phi21 = (-pi/2)+(idxphi21*pi)/number_of_intervals; + + for (int idxphi22 = 0; idxphi22 <= number_of_intervals; idxphi22++) + { + phi22 = (-pi/2)+(idxphi22*pi)/number_of_intervals; + + for (int idxphi11 = 0; idxphi11 <= number_of_intervals; idxphi11++) + { + phi11 = (-pi/2)+(idxphi11*pi)/number_of_intervals; + + for (int idxphi12 = 0; idxphi12 <= number_of_intervals; idxphi12++) + { + phi12 = (-pi/2)+(idxphi12*pi)/number_of_intervals; + + for (int idxtheta1 = 0; idxtheta1 <= number_of_intervals; idxtheta1++) + { + theta1 = -pi+(idxtheta1*(2*pi))/number_of_intervals; + + for (int idxtheta2 = 0; idxtheta2 <= number_of_intervals; idxtheta2++) + { + theta2 = -pi+(idxtheta2*(2*pi))/number_of_intervals; + + //::std::cout << "theta1 = " << theta1 << " ; "; + //::std::cout << "phi11 = " << phi11 << " ; "; + //::std::cout << "phi21 = " << phi21; + //::std::cout << "theta2 = " << theta2 << " ; "; + //::std::cout << "phi12 = " << phi12 << " ; "; + //::std::cout << "phi22 = " << phi22; + //::std::cout << ::std::endl; + + ::boost::math::quaternion p1 = ::boost::math::spherical(rho1, theta1, phi11, phi21); + + //::std::cout << "p1 = " << p1 << ::std::endl; + + ::boost::math::quaternion q1 = ::boost::math::spherical(rho2, theta2, phi12, phi22); + + //::std::cout << "q1 = " << q1 << ::std::endl; + + ::std::pair< ::boost::math::quaternion , ::boost::math::quaternion > pq1 = + ::std::make_pair(p1,q1); + + R4_matrix rot = quaternions_to_R4_rotation(pq1); + + //::std::cout << "rot = "; + //::std::cout << "\t" << rot.a11 << "\t" << rot.a12 << "\t" << rot.a13 << "\t" << rot.a14 << ::std::endl; + //::std::cout << "\t"; + //::std::cout << "\t" << rot.a21 << "\t" << rot.a22 << "\t" << rot.a23 << "\t" << rot.a24 << ::std::endl; + //::std::cout << "\t"; + //::std::cout << "\t" << rot.a31 << "\t" << rot.a32 << "\t" << rot.a33 << "\t" << rot.a34 << ::std::endl; + //::std::cout << "\t"; + //::std::cout << "\t" << rot.a41 << "\t" << rot.a42 << "\t" << rot.a43 << "\t" << rot.a44 << ::std::endl; + + ::std::pair< ::boost::math::quaternion , ::boost::math::quaternion > pq2 = + R4_rotation_to_quaternions(rot, &pq1); + + //::std::cout << "p1 = " << pq.first << ::std::endl; + //::std::cout << "p2 = " << pq.second << ::std::endl; + + //::std::cout << "round trip discrepancy: " << ::std::sqrt(::boost::math::norm(pq1.first-pq2.first)+::boost::math::norm(pq1.second-pq2.second)) << ::std::endl; + + //::std::cout << ::std::endl; + } + } + } + } + } + } + + ::std::cout << ::std::endl; +} + + +void test_SO4() +{ + ::std::cout << "Testing SO4:" << ::std::endl; + ::std::cout << ::std::endl; + + test_SO4_spherical(); } diff --git a/quaternion/HSO4.hpp b/quaternion/HSO4.hpp index 0045c4f03..909772cf0 100644 --- a/quaternion/HSO4.hpp +++ b/quaternion/HSO4.hpp @@ -26,157 +26,157 @@ template struct R4_matrix { - TYPE_FLOAT a11, a12, a13, a14; - TYPE_FLOAT a21, a22, a23, a24; - TYPE_FLOAT a31, a32, a33, a34; - TYPE_FLOAT a41, a42, a43, a44; + TYPE_FLOAT a11, a12, a13, a14; + TYPE_FLOAT a21, a22, a23, a24; + TYPE_FLOAT a31, a32, a33, a34; + TYPE_FLOAT a41, a42, a43, a44; }; -// Note: the input quaternions need not be of norm 1 for the following function +// Note: the input quaternions need not be of norm 1 for the following function template -R4_matrix quaternions_to_R4_rotation(::std::pair< ::boost::quaternion , ::boost::quaternion > const & pq) +R4_matrix quaternions_to_R4_rotation(::std::pair< ::boost::math::quaternion , ::boost::math::quaternion > const & pq) { - using ::std::numeric_limits; - - TYPE_FLOAT a0 = pq.first.R_component_1(); - TYPE_FLOAT b0 = pq.first.R_component_2(); - TYPE_FLOAT c0 = pq.first.R_component_3(); - TYPE_FLOAT d0 = pq.first.R_component_4(); - - TYPE_FLOAT norme_carre0 = a0*a0+b0*b0+c0*c0+d0*d0; - - if (norme_carre0 <= numeric_limits::epsilon()) - { - ::std::string error_reporting("Argument to quaternions_to_R4_rotation is too small!"); - ::std::underflow_error bad_argument(error_reporting); - - throw(bad_argument); - } - - TYPE_FLOAT a1 = pq.second.R_component_1(); - TYPE_FLOAT b1 = pq.second.R_component_2(); - TYPE_FLOAT c1 = pq.second.R_component_3(); - TYPE_FLOAT d1 = pq.second.R_component_4(); - - TYPE_FLOAT norme_carre1 = a1*a1+b1*b1+c1*c1+d1*d1; - - if (norme_carre1 <= numeric_limits::epsilon()) - { - ::std::string error_reporting("Argument to quaternions_to_R4_rotation is too small!"); - ::std::underflow_error bad_argument(error_reporting); - - throw(bad_argument); - } - - TYPE_FLOAT prod_norm = norme_carre0*norme_carre1; - - TYPE_FLOAT a0a1 = a0*a1; - TYPE_FLOAT a0b1 = a0*b1; - TYPE_FLOAT a0c1 = a0*c1; - TYPE_FLOAT a0d1 = a0*d1; - TYPE_FLOAT b0a1 = b0*a1; - TYPE_FLOAT b0b1 = b0*b1; - TYPE_FLOAT b0c1 = b0*c1; - TYPE_FLOAT b0d1 = b0*d1; - TYPE_FLOAT c0a1 = c0*a1; - TYPE_FLOAT c0b1 = c0*b1; - TYPE_FLOAT c0c1 = c0*c1; - TYPE_FLOAT c0d1 = c0*d1; - TYPE_FLOAT d0a1 = d0*a1; - TYPE_FLOAT d0b1 = d0*b1; - TYPE_FLOAT d0c1 = d0*c1; - TYPE_FLOAT d0d1 = d0*d1; - - R4_matrix out_matrix; - - out_matrix.a11 = (+a0a1+b0b1+c0c1+d0d1)/prod_norm; - out_matrix.a12 = (+a0b1-b0a1-c0d1+d0c1)/prod_norm; - out_matrix.a13 = (+a0c1+b0d1-c0a1-d0b1)/prod_norm; - out_matrix.a14 = (+a0d1-b0c1+c0b1-d0a1)/prod_norm; - out_matrix.a21 = (-a0b1+b0a1-c0d1+d0c1)/prod_norm; - out_matrix.a22 = (+a0a1+b0b1-c0c1-d0d1)/prod_norm; - out_matrix.a23 = (-a0d1+b0c1+c0b1-d0a1)/prod_norm; - out_matrix.a24 = (+a0c1+b0d1+c0a1+d0b1)/prod_norm; - out_matrix.a31 = (-a0c1+b0d1+c0a1-d0b1)/prod_norm; - out_matrix.a32 = (+a0d1+b0c1+c0b1+d0a1)/prod_norm; - out_matrix.a33 = (+a0a1-b0b1+c0c1-d0d1)/prod_norm; - out_matrix.a34 = (-a0b1-b0a1+c0d1+d0c1)/prod_norm; - out_matrix.a41 = (-a0d1-b0c1+c0b1+d0a1)/prod_norm; - out_matrix.a42 = (-a0c1+b0d1-c0a1+d0b1)/prod_norm; - out_matrix.a43 = (+a0b1+b0a1+c0d1+d0c1)/prod_norm; - out_matrix.a44 = (+a0a1-b0b1-c0c1+d0d1)/prod_norm; - - return(out_matrix); + using ::std::numeric_limits; + + TYPE_FLOAT a0 = pq.first.R_component_1(); + TYPE_FLOAT b0 = pq.first.R_component_2(); + TYPE_FLOAT c0 = pq.first.R_component_3(); + TYPE_FLOAT d0 = pq.first.R_component_4(); + + TYPE_FLOAT norme_carre0 = a0*a0+b0*b0+c0*c0+d0*d0; + + if (norme_carre0 <= numeric_limits::epsilon()) + { + ::std::string error_reporting("Argument to quaternions_to_R4_rotation is too small!"); + ::std::underflow_error bad_argument(error_reporting); + + throw(bad_argument); + } + + TYPE_FLOAT a1 = pq.second.R_component_1(); + TYPE_FLOAT b1 = pq.second.R_component_2(); + TYPE_FLOAT c1 = pq.second.R_component_3(); + TYPE_FLOAT d1 = pq.second.R_component_4(); + + TYPE_FLOAT norme_carre1 = a1*a1+b1*b1+c1*c1+d1*d1; + + if (norme_carre1 <= numeric_limits::epsilon()) + { + ::std::string error_reporting("Argument to quaternions_to_R4_rotation is too small!"); + ::std::underflow_error bad_argument(error_reporting); + + throw(bad_argument); + } + + TYPE_FLOAT prod_norm = norme_carre0*norme_carre1; + + TYPE_FLOAT a0a1 = a0*a1; + TYPE_FLOAT a0b1 = a0*b1; + TYPE_FLOAT a0c1 = a0*c1; + TYPE_FLOAT a0d1 = a0*d1; + TYPE_FLOAT b0a1 = b0*a1; + TYPE_FLOAT b0b1 = b0*b1; + TYPE_FLOAT b0c1 = b0*c1; + TYPE_FLOAT b0d1 = b0*d1; + TYPE_FLOAT c0a1 = c0*a1; + TYPE_FLOAT c0b1 = c0*b1; + TYPE_FLOAT c0c1 = c0*c1; + TYPE_FLOAT c0d1 = c0*d1; + TYPE_FLOAT d0a1 = d0*a1; + TYPE_FLOAT d0b1 = d0*b1; + TYPE_FLOAT d0c1 = d0*c1; + TYPE_FLOAT d0d1 = d0*d1; + + R4_matrix out_matrix; + + out_matrix.a11 = (+a0a1+b0b1+c0c1+d0d1)/prod_norm; + out_matrix.a12 = (+a0b1-b0a1-c0d1+d0c1)/prod_norm; + out_matrix.a13 = (+a0c1+b0d1-c0a1-d0b1)/prod_norm; + out_matrix.a14 = (+a0d1-b0c1+c0b1-d0a1)/prod_norm; + out_matrix.a21 = (-a0b1+b0a1-c0d1+d0c1)/prod_norm; + out_matrix.a22 = (+a0a1+b0b1-c0c1-d0d1)/prod_norm; + out_matrix.a23 = (-a0d1+b0c1+c0b1-d0a1)/prod_norm; + out_matrix.a24 = (+a0c1+b0d1+c0a1+d0b1)/prod_norm; + out_matrix.a31 = (-a0c1+b0d1+c0a1-d0b1)/prod_norm; + out_matrix.a32 = (+a0d1+b0c1+c0b1+d0a1)/prod_norm; + out_matrix.a33 = (+a0a1-b0b1+c0c1-d0d1)/prod_norm; + out_matrix.a34 = (-a0b1-b0a1+c0d1+d0c1)/prod_norm; + out_matrix.a41 = (-a0d1-b0c1+c0b1+d0a1)/prod_norm; + out_matrix.a42 = (-a0c1+b0d1-c0a1+d0b1)/prod_norm; + out_matrix.a43 = (+a0b1+b0a1+c0d1+d0c1)/prod_norm; + out_matrix.a44 = (+a0a1-b0b1-c0c1+d0d1)/prod_norm; + + return(out_matrix); } template -inline bool is_R4_rotation_matrix(R4_matrix const & mat) +inline bool is_R4_rotation_matrix(R4_matrix const & mat) { - using ::std::abs; - - using ::std::numeric_limits; - - return ( - !( - (abs(mat.a11*mat.a11+mat.a21*mat.a21+mat.a31*mat.a31+mat.a41*mat.a41 - static_cast(1)) > static_cast(10)*numeric_limits::epsilon())|| - (abs(mat.a11*mat.a12+mat.a21*mat.a22+mat.a31*mat.a32+mat.a41*mat.a42 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| - (abs(mat.a11*mat.a13+mat.a21*mat.a23+mat.a31*mat.a33+mat.a41*mat.a43 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| - (abs(mat.a11*mat.a14+mat.a21*mat.a24+mat.a31*mat.a34+mat.a41*mat.a44 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| - //(abs(mat.a11*mat.a12+mat.a21*mat.a22+mat.a31*mat.a32+mat.a41*mat.a42 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| - (abs(mat.a12*mat.a12+mat.a22*mat.a22+mat.a32*mat.a32+mat.a42*mat.a42 - static_cast(1)) > static_cast(10)*numeric_limits::epsilon())|| - (abs(mat.a12*mat.a13+mat.a22*mat.a23+mat.a32*mat.a33+mat.a42*mat.a43 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| - (abs(mat.a12*mat.a14+mat.a22*mat.a24+mat.a32*mat.a34+mat.a42*mat.a44 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| - //(abs(mat.a11*mat.a13+mat.a21*mat.a23+mat.a31*mat.a33+mat.a41*mat.a43 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| - //(abs(mat.a12*mat.a13+mat.a22*mat.a23+mat.a32*mat.a33+mat.a42*mat.a43 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| - (abs(mat.a13*mat.a13+mat.a23*mat.a23+mat.a33*mat.a33+mat.a43*mat.a43 - static_cast(1)) > static_cast(10)*numeric_limits::epsilon())|| - (abs(mat.a13*mat.a14+mat.a23*mat.a24+mat.a33*mat.a34+mat.a43*mat.a44 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| - //(abs(mat.a11*mat.a14+mat.a21*mat.a24+mat.a31*mat.a34+mat.a41*mat.a44 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| - //(abs(mat.a12*mat.a14+mat.a22*mat.a24+mat.a32*mat.a34+mat.a42*mat.a44 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| - //(abs(mat.a13*mat.a14+mat.a23*mat.a24+mat.a33*mat.a34+mat.a43*mat.a44 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| - (abs(mat.a14*mat.a14+mat.a24*mat.a24+mat.a34*mat.a34+mat.a44*mat.a44 - static_cast(1)) > static_cast(10)*numeric_limits::epsilon()) - ) - ); + using ::std::abs; + + using ::std::numeric_limits; + + return ( + !( + (abs(mat.a11*mat.a11+mat.a21*mat.a21+mat.a31*mat.a31+mat.a41*mat.a41 - static_cast(1)) > static_cast(10)*numeric_limits::epsilon())|| + (abs(mat.a11*mat.a12+mat.a21*mat.a22+mat.a31*mat.a32+mat.a41*mat.a42 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| + (abs(mat.a11*mat.a13+mat.a21*mat.a23+mat.a31*mat.a33+mat.a41*mat.a43 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| + (abs(mat.a11*mat.a14+mat.a21*mat.a24+mat.a31*mat.a34+mat.a41*mat.a44 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| + //(abs(mat.a11*mat.a12+mat.a21*mat.a22+mat.a31*mat.a32+mat.a41*mat.a42 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| + (abs(mat.a12*mat.a12+mat.a22*mat.a22+mat.a32*mat.a32+mat.a42*mat.a42 - static_cast(1)) > static_cast(10)*numeric_limits::epsilon())|| + (abs(mat.a12*mat.a13+mat.a22*mat.a23+mat.a32*mat.a33+mat.a42*mat.a43 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| + (abs(mat.a12*mat.a14+mat.a22*mat.a24+mat.a32*mat.a34+mat.a42*mat.a44 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| + //(abs(mat.a11*mat.a13+mat.a21*mat.a23+mat.a31*mat.a33+mat.a41*mat.a43 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| + //(abs(mat.a12*mat.a13+mat.a22*mat.a23+mat.a32*mat.a33+mat.a42*mat.a43 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| + (abs(mat.a13*mat.a13+mat.a23*mat.a23+mat.a33*mat.a33+mat.a43*mat.a43 - static_cast(1)) > static_cast(10)*numeric_limits::epsilon())|| + (abs(mat.a13*mat.a14+mat.a23*mat.a24+mat.a33*mat.a34+mat.a43*mat.a44 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| + //(abs(mat.a11*mat.a14+mat.a21*mat.a24+mat.a31*mat.a34+mat.a41*mat.a44 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| + //(abs(mat.a12*mat.a14+mat.a22*mat.a24+mat.a32*mat.a34+mat.a42*mat.a44 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| + //(abs(mat.a13*mat.a14+mat.a23*mat.a24+mat.a33*mat.a34+mat.a43*mat.a44 - static_cast(0)) > static_cast(10)*numeric_limits::epsilon())|| + (abs(mat.a14*mat.a14+mat.a24*mat.a24+mat.a34*mat.a34+mat.a44*mat.a44 - static_cast(1)) > static_cast(10)*numeric_limits::epsilon()) + ) + ); } template -::std::pair< ::boost::quaternion , ::boost::quaternion > R4_rotation_to_quaternions(R4_matrix const & rot, - ::std::pair< ::boost::quaternion , ::boost::quaternion > const * hint = 0) +::std::pair< ::boost::math::quaternion , ::boost::math::quaternion > R4_rotation_to_quaternions( R4_matrix const & rot, + ::std::pair< ::boost::math::quaternion , ::boost::math::quaternion > const * hint = 0) { - if (!is_R4_rotation_matrix(rot)) - { - ::std::string error_reporting("Argument to R4_rotation_to_quaternions is not an R^4 rotation matrix!"); - ::std::range_error bad_argument(error_reporting); - - throw(bad_argument); - } - - R3_matrix mat; - - mat.a11 = -rot.a31*rot.a42+rot.a32*rot.a41+rot.a22*rot.a11-rot.a21*rot.a12; - mat.a12 = -rot.a31*rot.a43+rot.a33*rot.a41+rot.a23*rot.a11-rot.a21*rot.a13; - mat.a13 = -rot.a31*rot.a44+rot.a34*rot.a41+rot.a24*rot.a11-rot.a21*rot.a14; - mat.a21 = -rot.a31*rot.a12-rot.a22*rot.a41+rot.a32*rot.a11+rot.a21*rot.a42; - mat.a22 = -rot.a31*rot.a13-rot.a23*rot.a41+rot.a33*rot.a11+rot.a21*rot.a43; - mat.a23 = -rot.a31*rot.a14-rot.a24*rot.a41+rot.a34*rot.a11+rot.a21*rot.a44; - mat.a31 = +rot.a31*rot.a22-rot.a12*rot.a41+rot.a42*rot.a11-rot.a21*rot.a32; - mat.a32 = +rot.a31*rot.a23-rot.a13*rot.a41+rot.a43*rot.a11-rot.a21*rot.a33; - mat.a33 = +rot.a31*rot.a24-rot.a14*rot.a41+rot.a44*rot.a11-rot.a21*rot.a34; - - ::boost::quaternion q = R3_rotation_to_quaternion(mat); - - ::boost::quaternion p = - ::boost::quaternion(rot.a11,rot.a12,rot.a13,rot.a14)*q; - - if ((hint != 0) && (abs(hint->second+q) < abs(hint->second-q))) - { - return(::std::make_pair(-p,-q)); - } - - return(::std::make_pair(p,q)); + if (!is_R4_rotation_matrix(rot)) + { + ::std::string error_reporting("Argument to R4_rotation_to_quaternions is not an R^4 rotation matrix!"); + ::std::range_error bad_argument(error_reporting); + + throw(bad_argument); + } + + R3_matrix mat; + + mat.a11 = -rot.a31*rot.a42+rot.a32*rot.a41+rot.a22*rot.a11-rot.a21*rot.a12; + mat.a12 = -rot.a31*rot.a43+rot.a33*rot.a41+rot.a23*rot.a11-rot.a21*rot.a13; + mat.a13 = -rot.a31*rot.a44+rot.a34*rot.a41+rot.a24*rot.a11-rot.a21*rot.a14; + mat.a21 = -rot.a31*rot.a12-rot.a22*rot.a41+rot.a32*rot.a11+rot.a21*rot.a42; + mat.a22 = -rot.a31*rot.a13-rot.a23*rot.a41+rot.a33*rot.a11+rot.a21*rot.a43; + mat.a23 = -rot.a31*rot.a14-rot.a24*rot.a41+rot.a34*rot.a11+rot.a21*rot.a44; + mat.a31 = +rot.a31*rot.a22-rot.a12*rot.a41+rot.a42*rot.a11-rot.a21*rot.a32; + mat.a32 = +rot.a31*rot.a23-rot.a13*rot.a41+rot.a43*rot.a11-rot.a21*rot.a33; + mat.a33 = +rot.a31*rot.a24-rot.a14*rot.a41+rot.a44*rot.a11-rot.a21*rot.a34; + + ::boost::math::quaternion q = R3_rotation_to_quaternion(mat); + + ::boost::math::quaternion p = + ::boost::math::quaternion(rot.a11,rot.a12,rot.a13,rot.a14)*q; + + if ((hint != 0) && (abs(hint->second+q) < abs(hint->second-q))) + { + return(::std::make_pair(-p,-q)); + } + + return(::std::make_pair(p,q)); } #endif /* TEST_HSO4_HPP */ diff --git a/quaternion/output.txt b/quaternion/output.txt index 39f44275dbd22760dd58c75b95d66b4b00ee021d..8688bc6a98b7f6c6dbb61de6b0fc4a5181761cef 100644 GIT binary patch delta 354 zcmbO&_&|6=2S+sUb;eBPJUvE9v6(8TbfgnS&);NoLG{XpN9}nttdcnlJkoZ zoZ?I*Wy$%)87NXnyp;UXq?}Z!Y<*t30uAlT$%>2ZwN8zCvnIQGQXeLP}~$YH~?x3Ksy7d -#include -#ifdef BOOST_NO_STDC_NAMESPACE -namespace std { - using ::sqrt; - using ::atan; - using ::log; - using ::exp; - using ::cos; - using ::sin; - using ::tan; - using ::cosh; - using ::sinh; - using ::tanh; -} -#endif // explicit (if ludicrous) instanciation #ifndef __GNUC__ @@ -236,7 +221,7 @@ int test_main(int, char *[]) << q0 << " ." << ::std::endl; } #else - ::std::istringstream bogus("(1,2,3,4)"); + ::std::istringstream bogus("(1,2,3,4)"); bogus >> q0; @@ -319,10 +304,10 @@ int test_main(int, char *[]) ::std::cout << ::std::endl; - float rho = ::std::sqrt(8.0f); - float theta = ::std::atan(1.0f); - float phi1 = ::std::atan(1.0f); - float phi2 = ::std::atan(1.0f); + float rho = ::std::sqrt(8.0f); + float theta = ::std::atan(1.0f); + float phi1 = ::std::atan(1.0f); + float phi2 = ::std::atan(1.0f); ::std::cout << "The value of the quaternion represented " << "in spherical form by " @@ -332,7 +317,7 @@ int test_main(int, char *[]) << ::boost::math::spherical(rho, theta, phi1, phi2) << ::std::endl; - float alpha = ::std::atan(1.0f); + float alpha = ::std::atan(1.0f); ::std::cout << "The value of the quaternion represented " << "in semipolar form by " @@ -342,10 +327,10 @@ int test_main(int, char *[]) << ::boost::math::semipolar(rho, alpha, phi1, phi2) << ::std::endl; - float rho1 = 1; - float rho2 = 2; - float theta1 = 0; - float theta2 = ::std::atan(1.0f)*2; + float rho1 = 1; + float rho2 = 2; + float theta1 = 0; + float theta2 = ::std::atan(1.0f)*2; ::std::cout << "The value of the quaternion represented " << "in multipolar form by " @@ -355,10 +340,10 @@ int test_main(int, char *[]) << ::boost::math::multipolar(rho1, theta1, rho2, theta2) << ::std::endl; - float t = 5; - float radius = ::std::sqrt(2.0f); - float longitude = ::std::atan(1.0f); - float lattitude = ::std::atan(::std::sqrt(3.0f)); + float t = 5; + float radius = ::std::sqrt(2.0f); + float longitude = ::std::atan(1.0f); + float lattitude = ::std::atan(::std::sqrt(3.0f)); ::std::cout << "The value of the quaternion represented " << "in cylindrospherical form by " @@ -366,13 +351,13 @@ int test_main(int, char *[]) << " , longitude = " << longitude << " , latitude = " << lattitude << " is " << ::boost::math::cylindrospherical(t, radius, - longitude, lattitude) + longitude, lattitude) << ::std::endl; - float r = ::std::sqrt(2.0f); - float angle = ::std::atan(1.0f); - float h1 = 3; - float h2 = 4; + float r = ::std::sqrt(2.0f); + float angle = ::std::atan(1.0f); + float h1 = 3; + float h2 = 4; ::std::cout << "The value of the quaternion represented " << "in cylindrical form by " @@ -431,164 +416,188 @@ int test_main(int, char *[]) using ::boost::math::abs; - #define BOOST_QUATERNION_MULTIPLICATION_TEST(type) \ - \ - BOOST_CRITICAL_TEST(abs(::boost::math::quaternion(1,0,0,0)* \ - ::boost::math::quaternion(1,0,0,0)- \ - static_cast(1)) <= \ - numeric_limits::epsilon()); \ - \ - BOOST_CRITICAL_TEST(abs(::boost::math::quaternion(0,1,0,0)* \ - ::boost::math::quaternion(0,1,0,0)+ \ - static_cast(1)) <= \ - numeric_limits::epsilon()); \ - \ - BOOST_CRITICAL_TEST(abs(::boost::math::quaternion(0,0,1,0)* \ - ::boost::math::quaternion(0,0,1,0)+ \ - static_cast(1)) <= \ - numeric_limits::epsilon()); \ - \ - BOOST_CRITICAL_TEST(abs(::boost::math::quaternion(0,0,0,1)* \ - ::boost::math::quaternion(0,0,0,1)+ \ - static_cast(1)) <= \ - numeric_limits::epsilon()); \ - \ - BOOST_CRITICAL_TEST(abs(::boost::math::quaternion(0,1,0,0)* \ - ::boost::math::quaternion(0,0,1,0)- \ - ::boost::math::quaternion(0,0,0,1)) <= \ - numeric_limits::epsilon()); \ - \ - BOOST_CRITICAL_TEST(abs(::boost::math::quaternion(0,0,1,0)* \ - ::boost::math::quaternion(0,1,0,0)+ \ - ::boost::math::quaternion(0,0,0,1)) <= \ - numeric_limits::epsilon()); \ - \ - BOOST_CRITICAL_TEST(abs(::boost::math::quaternion(0,0,1,0)* \ - ::boost::math::quaternion(0,0,0,1)- \ - ::boost::math::quaternion(0,1,0,0)) <= \ - numeric_limits::epsilon()); \ - \ - BOOST_CRITICAL_TEST(abs(::boost::math::quaternion(0,0,0,1)* \ - ::boost::math::quaternion(0,0,1,0)+ \ - ::boost::math::quaternion(0,1,0,0)) <= \ - numeric_limits::epsilon()); \ - \ - BOOST_CRITICAL_TEST(abs(::boost::math::quaternion(0,0,0,1)* \ - ::boost::math::quaternion(0,1,0,0)- \ - ::boost::math::quaternion(0,0,1,0)) <= \ - numeric_limits::epsilon()); \ - \ - BOOST_CRITICAL_TEST(abs(::boost::math::quaternion(0,1,0,0)* \ - ::boost::math::quaternion(0,0,0,1)+ \ - ::boost::math::quaternion(0,0,1,0)) <= \ - numeric_limits::epsilon()); +#define BOOST_QUATERNION_MULTIPLICATION_TEST(type) \ + \ + ::std::cout << "Testing multiplication." << std::endl; \ + \ + BOOST_CRITICAL_TEST(abs(::boost::math::quaternion(1,0,0,0)* \ + ::boost::math::quaternion(1,0,0,0)- \ + static_cast(1)) <= \ + numeric_limits::epsilon()); \ + \ + BOOST_CRITICAL_TEST(abs(::boost::math::quaternion(0,1,0,0)* \ + ::boost::math::quaternion(0,1,0,0)+ \ + static_cast(1)) <= \ + numeric_limits::epsilon()); \ + \ + BOOST_CRITICAL_TEST(abs(::boost::math::quaternion(0,0,1,0)* \ + ::boost::math::quaternion(0,0,1,0)+ \ + static_cast(1)) <= \ + numeric_limits::epsilon()); \ + \ + BOOST_CRITICAL_TEST(abs(::boost::math::quaternion(0,0,0,1)* \ + ::boost::math::quaternion(0,0,0,1)+ \ + static_cast(1)) <= \ + numeric_limits::epsilon()); \ + \ + BOOST_CRITICAL_TEST(abs(::boost::math::quaternion(0,1,0,0)* \ + ::boost::math::quaternion(0,0,1,0)- \ + ::boost::math::quaternion(0,0,0,1)) <= \ + numeric_limits::epsilon()); \ + \ + BOOST_CRITICAL_TEST(abs(::boost::math::quaternion(0,0,1,0)* \ + ::boost::math::quaternion(0,1,0,0)+ \ + ::boost::math::quaternion(0,0,0,1)) <= \ + numeric_limits::epsilon()); \ + \ + BOOST_CRITICAL_TEST(abs(::boost::math::quaternion(0,0,1,0)* \ + ::boost::math::quaternion(0,0,0,1)- \ + ::boost::math::quaternion(0,1,0,0)) <= \ + numeric_limits::epsilon()); \ + \ + BOOST_CRITICAL_TEST(abs(::boost::math::quaternion(0,0,0,1)* \ + ::boost::math::quaternion(0,0,1,0)+ \ + ::boost::math::quaternion(0,1,0,0)) <= \ + numeric_limits::epsilon()); \ + \ + BOOST_CRITICAL_TEST(abs(::boost::math::quaternion(0,0,0,1)* \ + ::boost::math::quaternion(0,1,0,0)- \ + ::boost::math::quaternion(0,0,1,0)) <= \ + numeric_limits::epsilon()); \ + \ + BOOST_CRITICAL_TEST(abs(::boost::math::quaternion(0,1,0,0)* \ + ::boost::math::quaternion(0,0,0,1)+ \ + ::boost::math::quaternion(0,0,1,0)) <= \ + numeric_limits::epsilon()); - #define BOOST_QUATERNION_EXP_TEST(type) \ - \ - BOOST_TEST(abs(exp(::boost::math::quaternion \ - (0,4*::std::atan(static_cast(1)),0,0) \ - )+static_cast(1)) <= 2*numeric_limits::epsilon()); \ - \ - BOOST_TEST(abs(exp(::boost::math::quaternion \ - (0,0,4*::std::atan(static_cast(1)),0) \ - )+static_cast(1)) <= 2*numeric_limits::epsilon()); \ - \ - BOOST_TEST(abs(exp(::boost::math::quaternion \ - (0,0,0,4*::std::atan(static_cast(1))) \ - )+static_cast(1)) <= 2*numeric_limits::epsilon()); +#define BOOST_QUATERNION_EXP_TEST(type) \ + \ + ::std::cout << "Testing exp." << std::endl; \ + \ + BOOST_TEST(abs(exp(::boost::math::quaternion \ + (0,4*::std::atan(static_cast(1)),0,0) \ + )+static_cast(1)) <= 2*numeric_limits::epsilon()); \ + \ + BOOST_TEST(abs(exp(::boost::math::quaternion \ + (0,0,4*::std::atan(static_cast(1)),0) \ + )+static_cast(1)) <= 2*numeric_limits::epsilon()); \ + \ + BOOST_TEST(abs(exp(::boost::math::quaternion \ + (0,0,0,4*::std::atan(static_cast(1))) \ + )+static_cast(1)) <= 2*numeric_limits::epsilon()); - #define BOOST_QUATERNION_COS_TEST(type) \ - \ - BOOST_TEST(abs(static_cast(4)* \ - cos(::boost::math::quaternion \ - (0,::std::log(static_cast(2)),0,0) \ - )-static_cast(5)) <= 4*numeric_limits::epsilon()); \ - \ - BOOST_TEST(abs(static_cast(4)* \ - cos(::boost::math::quaternion \ - (0,0,::std::log(static_cast(2)),0) \ - )-static_cast(5)) <= 4*numeric_limits::epsilon()); \ - \ - BOOST_TEST(abs(static_cast(4)* \ - cos(::boost::math::quaternion \ - (0,0,0,::std::log(static_cast(2))) \ - )-static_cast(5)) <= 4*numeric_limits::epsilon()); +#define BOOST_QUATERNION_COS_TEST(type) \ + \ + ::std::cout << "Testing cos." << std::endl; \ + \ + BOOST_TEST(abs(static_cast(4)* \ + cos(::boost::math::quaternion \ + (0,::std::log(static_cast(2)),0,0) \ + )-static_cast(5)) <= 4*numeric_limits::epsilon()); \ + \ + BOOST_TEST(abs(static_cast(4)* \ + cos(::boost::math::quaternion \ + (0,0,::std::log(static_cast(2)),0) \ + )-static_cast(5)) <= 4*numeric_limits::epsilon()); \ + \ + BOOST_TEST(abs(static_cast(4)* \ + cos(::boost::math::quaternion \ + (0,0,0,::std::log(static_cast(2))) \ + )-static_cast(5)) <= 4*numeric_limits::epsilon()); - #define BOOST_QUATERNION_SIN_TEST(type) \ - \ - BOOST_TEST(abs(static_cast(4)* \ - sin(::boost::math::quaternion \ - (0,::std::log(static_cast(2)),0,0) \ - )-::boost::math::quaternion(0,3,0,0)) <= \ - 4*numeric_limits::epsilon()); \ - \ - BOOST_TEST(abs(static_cast(4)* \ - sin(::boost::math::quaternion \ - (0,0,::std::log(static_cast(2)),0) \ - )-::boost::math::quaternion(0,0,3,0)) <= \ - 4*numeric_limits::epsilon()); \ - \ - BOOST_TEST(abs(static_cast(4)* \ - sin(::boost::math::quaternion \ - (0,0,0,::std::log(static_cast(2))) \ - )-::boost::math::quaternion(0,0,0,3)) <= \ - 4*numeric_limits::epsilon()); +#define BOOST_QUATERNION_SIN_TEST(type) \ + \ + ::std::cout << "Testing sin." << std::endl; \ + \ + BOOST_TEST(abs(static_cast(4)* \ + sin(::boost::math::quaternion \ + (0,::std::log(static_cast(2)),0,0) \ + )-::boost::math::quaternion(0,3,0,0)) <= \ + 4*numeric_limits::epsilon()); \ + \ + BOOST_TEST(abs(static_cast(4)* \ + sin(::boost::math::quaternion \ + (0,0,::std::log(static_cast(2)),0) \ + )-::boost::math::quaternion(0,0,3,0)) <= \ + 4*numeric_limits::epsilon()); \ + \ + BOOST_TEST(abs(static_cast(4)* \ + sin(::boost::math::quaternion \ + (0,0,0,::std::log(static_cast(2))) \ + )-::boost::math::quaternion(0,0,0,3)) <= \ + 4*numeric_limits::epsilon()); - #define BOOST_QUATERNION_COSH_TEST(type) \ - \ - BOOST_TEST(abs( \ - cosh(::boost::math::quaternion \ - (0,4*::std::atan(static_cast(1)),0,0) \ - )+static_cast(1)) <= 4*numeric_limits::epsilon()); \ - \ - BOOST_TEST(abs( \ - cosh(::boost::math::quaternion \ - (0,0,4*::std::atan(static_cast(1)),0) \ - )+static_cast(1)) <= 4*numeric_limits::epsilon()); \ - \ - BOOST_TEST(abs( \ - cosh(::boost::math::quaternion \ - (0,0,0,4*::std::atan(static_cast(1))) \ - )+static_cast(1)) <= 4*numeric_limits::epsilon()); +#define BOOST_QUATERNION_COSH_TEST(type) \ + \ + ::std::cout << "Testing cosh." << std::endl; \ + \ + BOOST_TEST(abs( \ + cosh(::boost::math::quaternion \ + (0,4*::std::atan(static_cast(1)),0,0) \ + )+static_cast(1)) <= 4*numeric_limits::epsilon()); \ + \ + BOOST_TEST(abs( \ + cosh(::boost::math::quaternion \ + (0,0,4*::std::atan(static_cast(1)),0) \ + )+static_cast(1)) <= 4*numeric_limits::epsilon()); \ + \ + BOOST_TEST(abs( \ + cosh(::boost::math::quaternion \ + (0,0,0,4*::std::atan(static_cast(1))) \ + )+static_cast(1)) <= 4*numeric_limits::epsilon()); - #define BOOST_QUATERNION_SINH_TEST(type) \ - \ - BOOST_TEST(abs( \ - sinh(::boost::math::quaternion \ - (0,2*::std::atan(static_cast(1)),0,0) \ - )-::boost::math::quaternion(0,1,0,0)) <= \ - 4*numeric_limits::epsilon()); \ - \ - BOOST_TEST(abs( \ - sinh(::boost::math::quaternion \ - (0,0,2*::std::atan(static_cast(1)),0) \ - )-::boost::math::quaternion(0,0,1,0)) <= \ - 4*numeric_limits::epsilon()); \ - \ - BOOST_TEST(abs( \ - sinh(::boost::math::quaternion \ - (0,0,0,2*::std::atan(static_cast(1))) \ - )-::boost::math::quaternion(0,0,0,1)) <= \ - 4*numeric_limits::epsilon()); +#define BOOST_QUATERNION_SINH_TEST(type) \ + \ + ::std::cout << "Testing sinh." << std::endl; \ + \ + BOOST_TEST(abs( \ + sinh(::boost::math::quaternion \ + (0,2*::std::atan(static_cast(1)),0,0) \ + )-::boost::math::quaternion(0,1,0,0)) <= \ + 4*numeric_limits::epsilon()); \ + \ + BOOST_TEST(abs( \ + sinh(::boost::math::quaternion \ + (0,0,2*::std::atan(static_cast(1)),0) \ + )-::boost::math::quaternion(0,0,1,0)) <= \ + 4*numeric_limits::epsilon()); \ + \ + BOOST_TEST(abs( \ + sinh(::boost::math::quaternion \ + (0,0,0,2*::std::atan(static_cast(1))) \ + )-::boost::math::quaternion(0,0,0,1)) <= \ + 4*numeric_limits::epsilon()); - #define BOOST_QUATERNION_TRENSCENDENTALS_TEST(type) \ - \ - BOOST_QUATERNION_EXP_TEST(type) \ - BOOST_QUATERNION_COS_TEST(type) \ - BOOST_QUATERNION_SIN_TEST(type) \ - BOOST_QUATERNION_COSH_TEST(type) \ - BOOST_QUATERNION_SINH_TEST(type) +#define BOOST_QUATERNION_TRENSCENDENTALS_TEST(type) \ + \ + BOOST_QUATERNION_EXP_TEST(type) \ + BOOST_QUATERNION_COS_TEST(type) \ + BOOST_QUATERNION_SIN_TEST(type) \ + BOOST_QUATERNION_COSH_TEST(type) \ + BOOST_QUATERNION_SINH_TEST(type) - #define BOOST_QUATERNION_TEST(type) \ - BOOST_QUATERNION_MULTIPLICATION_TEST(type) \ - BOOST_QUATERNION_TRENSCENDENTALS_TEST(type) +#if defined(__GNUC__) || defined(__COMO__) || \ + (defined(__MWERKS__) && (__MWERKS__ <= 0x2301)) +#define BOOST_QUATERNION_TEST(type) \ + \ + ::std::cout << "Testing " << #type << "." << std::endl; \ + \ + BOOST_QUATERNION_MULTIPLICATION_TEST(type) +#else +#define BOOST_QUATERNION_TEST(type) \ + \ + ::std::cout << "Testing " << #type << "." << std::endl; \ + \ + BOOST_QUATERNION_MULTIPLICATION_TEST(type) \ + BOOST_QUATERNION_TRENSCENDENTALS_TEST(type) +#endif BOOST_QUATERNION_TEST(float) @@ -596,16 +605,16 @@ int test_main(int, char *[]) BOOST_QUATERNION_TEST(long double) - #undef BOOST_QUATERNION_TEST +#undef BOOST_QUATERNION_TEST - #undef BOOST_QUATERNION_MULTIPLICATION_TEST - #undef BOOST_QUATERNION_TRENSCENDENTALS_TEST +#undef BOOST_QUATERNION_MULTIPLICATION_TEST +#undef BOOST_QUATERNION_TRENSCENDENTALS_TEST - #undef BOOST_QUATERNION_EXP_TEST - #undef BOOST_QUATERNION_COS_TEST - #undef BOOST_QUATERNION_SIN_TEST - #undef BOOST_QUATERNION_COSH_TEST - #undef BOOST_QUATERNION_SINH_TEST +#undef BOOST_QUATERNION_EXP_TEST +#undef BOOST_QUATERNION_COS_TEST +#undef BOOST_QUATERNION_SIN_TEST +#undef BOOST_QUATERNION_COSH_TEST +#undef BOOST_QUATERNION_SINH_TEST return(::boost::exit_success); diff --git a/special_functions/index.html b/special_functions/index.html index 50e68573f..46e7d41b6 100644 --- a/special_functions/index.html +++ b/special_functions/index.html @@ -41,13 +41,14 @@

The Special Functions library is a stop-gap collection of numerical functions, some of which are needed for our implementation of quaternions and octonions.

Comments

diff --git a/special_functions/output.txt b/special_functions/output.txt index fdbc7e419..f1a07117a 100644 --- a/special_functions/output.txt +++ b/special_functions/output.txt @@ -1,307 +1,25 @@ - -inf -10 -10 - -inf -9.8 -9.8 - -inf -9.6 -9.6 - -inf -9.4 -9.4 - -inf -9.2 -9.2 - -8.66434 -9 -9 - -8.66434 -8.8 -8.8 - -8.66434 -8.6 -8.6 - -8.31777 -8.4 -8.4 - -8.11503 -8.2 -8.2 - -7.97119 -8 -8 - -7.76846 -7.8 -7.8 - -7.62462 -7.6 -7.6 - -7.38187 -7.4 -7.4 - -7.19212 -7.2 -7.2 - -6.99824 -7 -7 - -6.7955 -6.8 -6.8 - -6.60077 -6.6 -6.6 - -6.39804 -6.4 -6.4 - -6.20071 -6.2 -6.2 - -6.0004 -6 -6 - -5.79929 -5.8 -5.8 - -5.59981 -5.6 -5.6 - -5.40035 -5.4 -5.4 - -5.20006 -5.2 -5.2 - -5.0001 -5 -5 - -4.8001 -4.8 -4.8 - -4.6 -4.6 -4.6 - -4.4 -4.4 -4.4 - -4.19997 -4.2 -4.2 - -4.00002 -4 -4 - -3.8 -3.8 -3.8 - -3.59999 -3.6 -3.6 - -3.39999 -3.4 -3.4 - -3.2 -3.2 -3.2 - -3 -3 -3 - -2.8 -2.8 -2.8 - -2.6 -2.6 -2.6 - -2.4 -2.4 -2.4 - -2.2 -2.2 -2.2 - -2 -2 -2 - -1.8 -1.8 -1.8 - -1.6 -1.6 -1.6 - -1.4 -1.4 -1.4 - -1.2 -1.2 -1.2 - -1 -1 -1 - -0.8 -0.8 -0.8 - -0.6 -0.6 -0.6 - -0.4 -0.4 -0.4 - -0.2 -0.2 -0.2 - 0 0 0 - 0.2 0.2 0.2 - 0.4 0.4 0.4 - 0.6 0.6 0.6 - 0.8 0.8 0.8 - 1 1 1 - 1.2 1.2 1.2 - 1.4 1.4 1.4 - 1.6 1.6 1.6 - 1.8 1.8 1.8 - 2 2 2 - 2.2 2.2 2.2 - 2.4 2.4 2.4 - 2.6 2.6 2.6 - 2.8 2.8 2.8 - 3 3 3 - 3.2 3.2 3.2 - 3.39999 3.4 3.4 - 3.59999 3.6 3.6 - 3.8 3.8 3.8 - 4.00002 4 4 - 4.19997 4.2 4.2 - 4.4 4.4 4.4 - 4.6 4.6 4.6 - 4.8001 4.8 4.8 - 5.0001 5 5 - 5.20006 5.2 5.2 - 5.40035 5.4 5.4 - 5.59981 5.6 5.6 - 5.79929 5.8 5.8 - 6.0004 6 6 - 6.20071 6.2 6.2 - 6.39804 6.4 6.4 - 6.60077 6.6 6.6 - 6.7955 6.8 6.8 - 6.99824 7 7 - 7.19212 7.2 7.2 - 7.38187 7.4 7.4 - 7.62462 7.6 7.6 - 7.76846 7.8 7.8 - 7.97119 8 8 - 8.11503 8.2 8.2 - 8.31777 8.4 8.4 - 8.66434 8.6 8.6 - 8.66434 8.8 8.8 - 8.66434 9 9 - inf 9.2 9.2 - inf 9.4 9.4 - inf 9.6 9.6 - inf 9.8 9.8 - inf 10 10 - - 0.841471 0.841471 0.841471 - 0.847446 0.847446 0.847446 - 0.853325 0.853325 0.853325 - 0.859104 0.859104 0.859104 - 0.864784 0.864784 0.864784 - 0.870363 0.870363 0.870363 - 0.87584 0.87584 0.87584 - 0.881212 0.881212 0.881212 - 0.88648 0.88648 0.88648 - 0.891641 0.891641 0.891641 - 0.896695 0.896695 0.896695 - 0.90164 0.90164 0.90164 - 0.906476 0.906476 0.906476 - 0.9112 0.9112 0.9112 - 0.915812 0.915812 0.915812 - 0.920311 0.920311 0.920311 - 0.924696 0.924696 0.924696 - 0.928965 0.928965 0.928965 - 0.933118 0.933118 0.933118 - 0.937153 0.937153 0.937153 - 0.941071 0.941071 0.941071 - 0.944869 0.944869 0.944869 - 0.948547 0.948547 0.948547 - 0.952104 0.952104 0.952104 - 0.955539 0.955539 0.955539 - 0.958851 0.958851 0.958851 - 0.96204 0.96204 0.96204 - 0.965105 0.965105 0.965105 - 0.968044 0.968044 0.968044 - 0.970858 0.970858 0.970858 - 0.973546 0.973546 0.973546 - 0.976107 0.976106 0.976106 - 0.97854 0.97854 0.97854 - 0.980844 0.980844 0.980844 - 0.98302 0.983021 0.983021 - 0.985067 0.985067 0.985067 - 0.986984 0.986984 0.986984 - 0.988771 0.988771 0.988771 - 0.990428 0.990428 0.990428 - 0.991953 0.991953 0.991953 - 0.993347 0.993347 0.993347 - 0.994609 0.994609 0.994609 - 0.995739 0.995739 0.995739 - 0.996737 0.996737 0.996737 - 0.997602 0.997602 0.997602 - 0.998334 0.998334 0.998334 - 0.998934 0.998934 0.998934 - 0.9994 0.9994 0.9994 - 0.999733 0.999733 0.999733 - 0.999933 0.999933 0.999933 - 1 1 1 - 0.999933 0.999933 0.999933 - 0.999733 0.999733 0.999733 - 0.9994 0.9994 0.9994 - 0.998934 0.998934 0.998934 - 0.998334 0.998334 0.998334 - 0.997602 0.997602 0.997602 - 0.996737 0.996737 0.996737 - 0.995739 0.995739 0.995739 - 0.994609 0.994609 0.994609 - 0.993347 0.993347 0.993347 - 0.991953 0.991953 0.991953 - 0.990428 0.990428 0.990428 - 0.988771 0.988771 0.988771 - 0.986984 0.986984 0.986984 - 0.985067 0.985067 0.985067 - 0.98302 0.983021 0.983021 - 0.980844 0.980844 0.980844 - 0.97854 0.97854 0.97854 - 0.976107 0.976106 0.976106 - 0.973546 0.973546 0.973546 - 0.970858 0.970858 0.970858 - 0.968044 0.968044 0.968044 - 0.965105 0.965105 0.965105 - 0.96204 0.96204 0.96204 - 0.958851 0.958851 0.958851 - 0.955539 0.955539 0.955539 - 0.952104 0.952104 0.952104 - 0.948547 0.948547 0.948547 - 0.944869 0.944869 0.944869 - 0.941071 0.941071 0.941071 - 0.937153 0.937153 0.937153 - 0.933118 0.933118 0.933118 - 0.928965 0.928965 0.928965 - 0.924696 0.924696 0.924696 - 0.920311 0.920311 0.920311 - 0.915812 0.915812 0.915812 - 0.9112 0.9112 0.9112 - 0.906476 0.906476 0.906476 - 0.90164 0.90164 0.90164 - 0.896695 0.896695 0.896695 - 0.891641 0.891641 0.891641 - 0.88648 0.88648 0.88648 - 0.881212 0.881212 0.881212 - 0.87584 0.87584 0.87584 - 0.870363 0.870363 0.870363 - 0.864784 0.864784 0.864784 - 0.859104 0.859104 0.859104 - 0.853325 0.853325 0.853325 - 0.847446 0.847446 0.847446 - 0.841471 0.841471 0.841471 - - 1.1752 1.1752 1.1752 - 1.16793 1.16793 1.16793 - 1.16084 1.16084 1.16084 - 1.15391 1.15391 1.15391 - 1.14716 1.14716 1.14716 - 1.14057 1.14057 1.14057 - 1.13416 1.13416 1.13416 - 1.12791 1.12791 1.12791 - 1.12182 1.12182 1.12182 - 1.1159 1.1159 1.1159 - 1.11013 1.11013 1.11013 - 1.10453 1.10453 1.10453 - 1.09909 1.09909 1.09909 - 1.0938 1.0938 1.0938 - 1.08867 1.08867 1.08867 - 1.08369 1.08369 1.08369 - 1.07887 1.07887 1.07887 - 1.0742 1.0742 1.0742 - 1.06968 1.06968 1.06968 - 1.06531 1.06531 1.06531 - 1.06109 1.06109 1.06109 - 1.05702 1.05702 1.05702 - 1.05309 1.05309 1.05309 - 1.04931 1.04931 1.04931 - 1.04568 1.04568 1.04568 - 1.04219 1.04219 1.04219 - 1.03884 1.03884 1.03884 - 1.03564 1.03564 1.03564 - 1.03258 1.03258 1.03258 - 1.02966 1.02966 1.02966 - 1.02688 1.02688 1.02688 - 1.02424 1.02424 1.02424 - 1.02174 1.02174 1.02174 - 1.01938 1.01938 1.01938 - 1.01715 1.01715 1.01715 - 1.01507 1.01507 1.01507 - 1.01312 1.01312 1.01312 - 1.0113 1.0113 1.0113 - 1.00963 1.00963 1.00963 - 1.00809 1.00809 1.00809 - 1.00668 1.00668 1.00668 - 1.00541 1.00541 1.00541 - 1.00427 1.00427 1.00427 - 1.00327 1.00327 1.00327 - 1.0024 1.0024 1.0024 - 1.00167 1.00167 1.00167 - 1.00107 1.00107 1.00107 - 1.0006 1.0006 1.0006 - 1.00027 1.00027 1.00027 - 1.00007 1.00007 1.00007 - 1 1 1 - 1.00007 1.00007 1.00007 - 1.00027 1.00027 1.00027 - 1.0006 1.0006 1.0006 - 1.00107 1.00107 1.00107 - 1.00167 1.00167 1.00167 - 1.0024 1.0024 1.0024 - 1.00327 1.00327 1.00327 - 1.00427 1.00427 1.00427 - 1.00541 1.00541 1.00541 - 1.00668 1.00668 1.00668 - 1.00809 1.00809 1.00809 - 1.00963 1.00963 1.00963 - 1.0113 1.0113 1.0113 - 1.01312 1.01312 1.01312 - 1.01507 1.01507 1.01507 - 1.01715 1.01715 1.01715 - 1.01938 1.01938 1.01938 - 1.02174 1.02174 1.02174 - 1.02424 1.02424 1.02424 - 1.02688 1.02688 1.02688 - 1.02966 1.02966 1.02966 - 1.03258 1.03258 1.03258 - 1.03564 1.03564 1.03564 - 1.03884 1.03884 1.03884 - 1.04219 1.04219 1.04219 - 1.04568 1.04568 1.04568 - 1.04931 1.04931 1.04931 - 1.05309 1.05309 1.05309 - 1.05702 1.05702 1.05702 - 1.06109 1.06109 1.06109 - 1.06531 1.06531 1.06531 - 1.06968 1.06968 1.06968 - 1.0742 1.0742 1.0742 - 1.07887 1.07887 1.07887 - 1.08369 1.08369 1.08369 - 1.08867 1.08867 1.08867 - 1.0938 1.0938 1.0938 - 1.09909 1.09909 1.09909 - 1.10453 1.10453 1.10453 - 1.11013 1.11013 1.11013 - 1.1159 1.1159 1.1159 - 1.12182 1.12182 1.12182 - 1.12791 1.12791 1.12791 - 1.13416 1.13416 1.13416 - 1.14057 1.14057 1.14057 - 1.14716 1.14716 1.14716 - 1.15391 1.15391 1.15391 - 1.16084 1.16084 1.16084 - 1.16793 1.16793 1.16793 - 1.1752 1.1752 1.1752 - +Testing float. +Testing atanh in the real domain. +Testing asinh in the real domain. +Testing acosh in the real domain. +Testing sinc_pi in the real domain. +Testing sinhc_pi in the real domain. +Testing sinc_pi in the complex domain. +Testing sinhc_pi in the complex domain. +Testing double. +Testing atanh in the real domain. +Testing asinh in the real domain. +Testing acosh in the real domain. +Testing sinc_pi in the real domain. +Testing sinhc_pi in the real domain. +Testing sinc_pi in the complex domain. +Testing sinhc_pi in the complex domain. +Testing long double. +Testing atanh in the real domain. +Testing asinh in the real domain. +Testing acosh in the real domain. +Testing sinc_pi in the real domain. +Testing sinhc_pi in the real domain. +Testing sinc_pi in the complex domain. +Testing sinhc_pi in the complex domain. no errors detected diff --git a/special_functions/special_functions.html b/special_functions/special_functions.html index 98cd74dcb..69d28c4b7 100644 --- a/special_functions/special_functions.html +++ b/special_functions/special_functions.html @@ -16,9 +16,9 @@
  • Functions
  • History -

    The Special Functions library currently provides three templated special functions, in namespace boost. Two of these (sinc_pi and sinhc_pi) are needed by our implementation of quaternions and octonions.

    -

    The function atanh is entirely classical, the function sinc_pi sees heavy use in signal processing tasks, and the function sinhc_pi is an ad'hoc function whose naming is modelled on sinc_pi and hyperbolic functions.

    -

    The exponential function is defined, for all object for which this makes sense, as the power series , with (and by definition) being the factorial of . In particular, the exponential function is well defined for real numbers, complex number, quaternions, octonions, and matrices of complex numbers, among others.

    +

    The Special Functions library currently provides five templated special functions, in namespace boost. Two of these (sinc_pi and sinhc_pi) are needed by our implementation of quaternions and octonions.

    +

    The functions acosh, asinh and atanh are entirely classical, the function sinc_pi sees heavy use in signal processing tasks, and the function sinhc_pi is an ad'hoc function whose naming is modelled on sinc_pi and hyperbolic functions.

    +

    The exponential funtion is defined, for all object for which this makes sense, as the power series , with (and by definition) being the factorial of . In particular, the exponential function is well defined for real numbers, complex number, quaternions, octonions, and matrices of complex numbers, among others.

    Graph of exp on R

    @@ -28,11 +28,12 @@

    Real and Imaginary parts of exp on C

    +

    The hyperbolic functions are defined as power series which can be computed (for reals, complex, quaternions and octonions) as:

     

    -

    Hyperbolic cosinus:

    -

    Hyperbolic sinus:

    -

    Hyperbolic tangent:

    +

    Hyperbolic cosine:

    +

    Hyperbolic sine:

    +

    Hyperbolic tangent:

     

    @@ -40,10 +41,13 @@

    Hyperbolic functions on r (cosh: purple; sinh: red; tanh: blue)

    -

    The hyperbolic sinus is one to one on the set of real numbers, with range the full set of reals, while the hyperbolic tangent is also one to one on the set of real numbers but with range , and therefore both have inverses.

    -

    The inverse of the hyperbolic tangent is called the Argument hyperbolic tangent, and can be computed as .

    -

    The Sinus Cardinal family of functions is defined by the family of indices by . It is also an entire function (sum of a power series).

    -

    We define, by analogy, the Hyperbolic Sinus Cardinal family of functions is defined by the family of indices by . It, too, is an entire function.

    +

    +

    The hyperbolic sine is one to one on the set of real numbers, with range the full set of reals, while the hyperbolic tangent is also one to one on the set of real numbers but with range , and therefore both have inverses. The hyperbolic cosine is one to one from onto (and from onto ); the inverse function we use here is defined on with range .

    +

    The inverse of the hyperbolic tangent is called the Argument hyperbolic tangent, and can be computed as .

    +

    The inverse of the hyperbolic sine is called the Argument hyperbolic sine, and can be computed (for ) as .

    +

    The inverse of the hyperbolic cosine is called the Argument hyperbolic cosine, and can be computed as .

    +

    The sine Cardinal family of functions is defined by the family of indices by .

    +

    We define, by analogy, the Hyperbolic sine Cardinal family of functions is defined by the family of indices by . It, too, is an entire function.

    Sinus Cardinal of index pi (purple) and Hyperbolic Sinus Cardinal of index pi (red) on R

    @@ -53,6 +57,8 @@

    Header Files

    The interface and implementation for each function (or forms of a function) are both supplied by one header file:

      +
    • acosh.hpp +
    • asinh.hpp
    • atahn.hpp
    • sinc.hpp
    • sinhc.hpp @@ -72,6 +78,10 @@ namespace math { + template<typename T> inline T acosh(const T x); + + template<typename T> inline T asinh(const T x); + template<typename T> inline T atanh(const T x); template<typename T> inline T sinc_pi(const T x); @@ -90,29 +100,42 @@ }

      Functions

      The functions implemented here can throw standard exceptions, but no exception specification has been made.

      +

      acosh

      + 
      template<typename T> inline T acosh(const T x);
      +
      +

      Computes the reciprocal of (the restriction to the range of) the hyperbolic cosine function, at x. Values returned are positive. Generalised Taylor series are used near 1 and Laurent series are used near the infinity to ensure accuracy.

      +

      If x is in the range a quiet NaN is returned (if the system allows, otherwise a domain_error exception is generated).

      +
      +

      asinh

      + 
      template<typename T> inline T asinh(const T x);
      +
      +

      Computes the reciprocal of the hyperbolic sine function. Taylor series are used at the origin and Laurent series are used near the infinity to ensure accuracy.

      +

      atanh

      template<typename T> inline T atanh(const T x);
      -

      Computes the reciprocal of the hyperbolic tangent function, at x.

      -

      If x is outside the range ] -1 ; +1 [ , an out_of_range exception is generated.

      -

      For the value -1 , minus infinity is returned (if the system allow, otherwise an out_of_range exception is generated).

      -

      For the value +1 , plus infinity is returned (if the system allow, otherwise an out_of_range exception is generated).

      +

      Computes the reciprocal of the hyperbolic tangent function, at x. Taylor series are used at the origin to ensure accuracy.

      +

      If x is in the range or in the range a quiet NaN is returned (if the system allows, otherwise a domain_error exception is generated).

      +

      If x is in the range , minus infinity is returned (if the system allows, otherwise an out_of_range exception is generated), with denoting numeric_limits<T>::epsilon().

      +

      If x is in the range , plus infinity is returned (if the system allows, otherwise an out_of_range exception is generated), with denoting numeric_limits<T>::epsilon().

      sinc_pi

      template<typename T> inline T 																															sinc_pi(const T x);
      template<typename T, template<typename> class U> inline U<T> sinc_pi(const U<T> x);
      -

      Computes the Sinus Cardinal of x. The second form is for complexes, quaternions, octonions...

      +

      Computes the Sinus Cardinal of x. The second form is for complexes, quaternions, octonions... Taylor series are used at the origin to ensure accuracy.

      sinhc_pi

      template<typename T> inline T																								 							sinhc_pi(const T x);
      template<typename T, template<typename> class U> inline U<T> sinhc_pi(const U<T> x);
      -

      Computes the Hyperbolic Sinus Cardinal of x. The second form is for complexes, quaternions, octonions...

      +

      Computes the Hyperbolic Sinus Cardinal of x. The second form is for complexes, quaternions, octonions... Taylor series are used at the origin to ensure accuracy.

      History

        -
      • 1.3.2- 07/07/2001: introduced namespace math. +
      • 1.4.0 - 14/09/2001: added (after rewrite) acosh and asinh, which were submited by Eric Ford; applied changes for Gcc 2.9.x suggested by John Maddock; improved accuracy; sanity check for test file, related to accuracy. +
      • 1.3.2 - 07/07/2001: introduced namespace math. +
      • 1.3.1 - 07/06/2001:(end of Boost review) split special_functions.hpp into atanh.hpp, sinc.hpp and sinhc.hpp; improved efficiency of atanh with compile-time technique (Daryle Walker); improved accuracy of all functions near zero (Peter Schmitteckert).
      • 1.3.0 - 26/03/2001: support for complexes & all, for cardinal functions. @@ -122,7 +145,7 @@
      • 1.0.0 - 10/08/1999: first public version.
      -

      Revised 11 July 2001

      +

      Revised 14 Sep 2001

      © Copyright Hubert Holin 2001. Permission to copy, use, modify, sell and distribute this document is granted provided this copyright notice appears in all copies. This software is provided "as is" without express or implied  warranty, and with no claim as to its suitability for any purpose.

      diff --git a/special_functions/special_functions_test.cpp b/special_functions/special_functions_test.cpp index 8b0ac0aa6..a644929ec 100644 --- a/special_functions/special_functions_test.cpp +++ b/special_functions/special_functions_test.cpp @@ -6,183 +6,233 @@ // warranty, and with no claim as to its suitability for any purpose. +#include #include #include #include #include +#include +#include #include #include -#include -#ifdef BOOST_NO_STDC_NAMESPACE -namespace std { - using ::abs; - using ::tanh; - using ::log; - using ::sinh; - using ::sin; - using ::sqrt; -} -#endif #define BOOST_INCLUDE_MAIN // for testing, include rather than link #include +template +T atanh_error_evaluator(T x) +{ + using ::std::abs; + using ::std::tanh; + using ::std::cosh; + + using ::std::numeric_limits; + + using ::boost::math::atanh; + + + static T const epsilon = numeric_limits::epsilon(); + + T y = tanh(x); + T z = atanh(y); + + T absolute_error = abs(z-x); + T relative_error = absolute_error/(cosh(x)*cosh(x)); + T scaled_error = relative_error/epsilon; + + return(scaled_error); +} + + +template +T asinh_error_evaluator(T x) +{ + using ::std::abs; + using ::std::sinh; + using ::std::cosh; + + using ::std::numeric_limits; + + using ::boost::math::asinh; + + + static T const epsilon = numeric_limits::epsilon(); + + T y = sinh(x); + T z = asinh(y); + + T absolute_error = abs(z-x); + T relative_error = absolute_error*cosh(x); + T scaled_error = relative_error/epsilon; + + return(scaled_error); +} + + +template +T acosh_error_evaluator(T x) +{ + using ::std::abs; + using ::std::sinh; + using ::std::cosh; + + using ::std::numeric_limits; + + using ::boost::math::acosh; + + + static T const epsilon = numeric_limits::epsilon(); + + T y = cosh(x); + T z = acosh(y); + + T absolute_error = abs(z-abs(x)); + T relative_error = absolute_error*abs(sinh(x)); + T scaled_error = relative_error/epsilon; + + return(scaled_error); +} + + int test_main(int, char *[]) { - using ::std::abs; - using ::std::tanh; - using ::std::log; - using ::std::sinh; - using ::std::sin; + using ::std::abs; + using ::std::sin; + using ::std::log; + using ::std::tanh; + using ::std::sinh; + using ::std::cosh; - using ::std::numeric_limits; + using ::std::numeric_limits; - using ::boost::math::atanh; - using ::boost::math::sinc_pi; - using ::boost::math::sinhc_pi; + using ::boost::math::atanh; + using ::boost::math::asinh; + using ::boost::math::acosh; + using ::boost::math::sinc_pi; + using ::boost::math::sinhc_pi; // tests for evaluation by humans - - for (int i = 0; i <= 100; i++) - { - float xf = - tanh(static_cast(i-50)/static_cast(5)); - - double xd = - tanh(static_cast(i-50)/static_cast(5)); - - long double xl = - tanh(static_cast(i-50)/static_cast(5)); - - if ( - std::numeric_limits::has_infinity && - std::numeric_limits::has_infinity && - std::numeric_limits::has_infinity - ) - { - ::std::cout << ::std::setw(15) - << atanh(xf) - << ::std::setw(15) - << atanh(xd) - << ::std::setw(15) - << atanh(xl) - << ::std::endl; - } - else - { - if ( - (abs(xf-static_cast(1)) < - numeric_limits::epsilon())|| - (abs(xf+static_cast(1)) < - numeric_limits::epsilon())|| - (abs(xf-static_cast(1)) < - numeric_limits::epsilon())|| - (abs(xf+static_cast(1)) < - numeric_limits::epsilon())|| - (abs(xf-static_cast(1)) < - numeric_limits::epsilon())|| - (abs(xf+static_cast(1)) < - numeric_limits::epsilon()) - ) - { - ::std::cout << "Platform's numerics may lack precision." - << ::std::endl; - } - else - { - ::std::cout << ::std::setw(15) - << atanh(xf) - << ::std::setw(15) - << atanh(xd) - << ::std::setw(15) - << atanh(xl) - << ::std::endl; - } - } - } - ::std::cout << ::std::endl; - - for (int i = 0; i <= 100; i++) - { - ::std::cout << ::std::setw(15) - << sinc_pi(static_cast(i-50)/ - static_cast(50)) - << ::std::setw(15) - << sinc_pi(static_cast(i-50)/ - static_cast(50)) - << ::std::setw(15) - << sinc_pi(static_cast(i-50)/ - static_cast(50)) - << ::std::endl; - } - - ::std::cout << ::std::endl; - - for (int i = 0; i <= 100; i++) - { - ::std::cout << ::std::setw(15) - << sinhc_pi(static_cast(i-50)/ - static_cast(50)) - << ::std::setw(15) - << sinhc_pi(static_cast(i-50)/ - static_cast(50)) - << ::std::setw(15) - << sinhc_pi(static_cast(i-50)/ - static_cast(50)) - << ::std::endl; - } - - ::std::cout << ::std::endl; // tests for evaluation by scripts -#define BOOST_SPECIAL_FUNCTIONS_TEST_ATANH(type) \ - \ - BOOST_TEST(abs(atanh(static_cast(0))) <= \ - numeric_limits::epsilon()); \ - \ - BOOST_TEST(abs(atanh(static_cast(3)/5)- \ - log(static_cast(2))) <= \ - numeric_limits::epsilon()); \ - \ - BOOST_TEST(abs(atanh(static_cast(-3)/5)+ \ - log(static_cast(2))) <= \ +#define BOOST_SPECIAL_FUNCTIONS_TEST_ATANH(type) \ + \ + ::std::cout << "Testing atanh in the real domain." << ::std::endl; \ + \ + BOOST_TEST(abs(atanh(static_cast(0))) <= \ + numeric_limits::epsilon()); \ + \ + BOOST_TEST(abs(atanh(static_cast(3)/5)- \ + log(static_cast(2))) <= \ + numeric_limits::epsilon()); \ + \ + BOOST_TEST(abs(atanh(static_cast(-3)/5)+ \ + log(static_cast(2))) <= \ + numeric_limits::epsilon()); \ + \ + for (int i = 0; i <= 100; i++) \ + { \ + type x = static_cast(i-50)/static_cast(5); \ + type y = tanh(x); \ + \ + if ( \ + (abs(y-static_cast(1)) >= \ + numeric_limits::epsilon())&& \ + (abs(y+static_cast(1)) >= \ + numeric_limits::epsilon()) \ + ) \ + { \ + BOOST_TEST(atanh_error_evaluator(x) <= static_cast(4)); \ + } \ + } + + + +#define BOOST_SPECIAL_FUNCTIONS_TEST_ASINH(type) \ + \ + ::std::cout << "Testing asinh in the real domain." << ::std::endl; \ + \ + for (int i = 0; i <= 100; i++) \ + { \ + type x = static_cast(i-50)/static_cast(5); \ + \ + BOOST_TEST(asinh_error_evaluator(x) <= static_cast(4)); \ + } + + +#define BOOST_SPECIAL_FUNCTIONS_TEST_ACOSH(type) \ + \ + ::std::cout << "Testing acosh in the real domain." << ::std::endl; \ + \ + for (int i = 0; i <= 100; i++) \ + { \ + type x = static_cast(i-50)/static_cast(5); \ + \ + BOOST_TEST(acosh_error_evaluator(x) <= static_cast(4)); \ + } + + +#define BOOST_SPECIAL_FUNCTIONS_TEST_SINC_PI(type) \ + \ + ::std::cout << "Testing sinc_pi in the real domain." << ::std::endl; \ + \ + BOOST_TEST(abs(sinc_pi(static_cast(0))- \ + static_cast(1)) <= \ + numeric_limits::epsilon()); + +#define BOOST_SPECIAL_FUNCTIONS_TEST_SINC_PI_FOR_C(type) \ + \ + ::std::cout << "Testing sinc_pi in the complex domain." << ::std::endl; \ + \ + BOOST_TEST(abs(sinc_pi(::std::complex(0, 1))- \ + ::std::complex(sinh(static_cast(1)))) <= \ numeric_limits::epsilon()); -#define BOOST_SPECIAL_FUNCTIONS_TEST_SINC_PI(type) \ - \ - BOOST_TEST(abs(sinc_pi(static_cast(0))- \ - static_cast(1)) <= \ - numeric_limits::epsilon()); \ - \ - BOOST_TEST(abs(sinc_pi(::std::complex(0, 1))- \ - ::std::complex(sinh(static_cast(1)))) <= \ +#define BOOST_SPECIAL_FUNCTIONS_TEST_SINHC_PI(type) \ + \ + ::std::cout << "Testing sinhc_pi in the real domain." << ::std::endl; \ + \ + BOOST_TEST(abs(sinhc_pi(static_cast(0))- \ + static_cast(1)) <= \ + numeric_limits::epsilon()); + +#define BOOST_SPECIAL_FUNCTIONS_TEST_SINHC_PI_FOR_C(type) \ + \ + ::std::cout << "Testing sinhc_pi in the complex domain." << ::std::endl;\ + \ + BOOST_TEST(abs(sinhc_pi(::std::complex(0, 1))- \ + ::std::complex(sin(static_cast(1)))) <= \ numeric_limits::epsilon()); -#define BOOST_SPECIAL_FUNCTIONS_TEST_SINHC_PI(type) \ - \ - BOOST_TEST(abs(sinhc_pi(static_cast(0))- \ - static_cast(1)) <= \ - numeric_limits::epsilon()); \ - \ - BOOST_TEST(abs(sinhc_pi(::std::complex(0, 1))- \ - ::std::complex(sin(static_cast(1)))) <= \ - numeric_limits::epsilon()); - - -#define BOOST_SPECIAL_FUNCTIONS_TEST(type) \ - BOOST_SPECIAL_FUNCTIONS_TEST_ATANH(type) \ - BOOST_SPECIAL_FUNCTIONS_TEST_SINC_PI(type) \ +#if defined(__GNUC__) || defined(__COMO__) || \ + (defined(__MWERKS__) && (__MWERKS__ <= 0x2301)) +#define BOOST_SPECIAL_FUNCTIONS_TEST(type) \ + ::std::cout << "Testing " << #type << "." << std::endl; \ + BOOST_SPECIAL_FUNCTIONS_TEST_ATANH(type) \ + BOOST_SPECIAL_FUNCTIONS_TEST_ASINH(type) \ + BOOST_SPECIAL_FUNCTIONS_TEST_ACOSH(type) \ + BOOST_SPECIAL_FUNCTIONS_TEST_SINC_PI(type) \ BOOST_SPECIAL_FUNCTIONS_TEST_SINHC_PI(type) +#else +#define BOOST_SPECIAL_FUNCTIONS_TEST(type) \ + ::std::cout << "Testing " << #type << "." << std::endl; \ + BOOST_SPECIAL_FUNCTIONS_TEST_ATANH(type) \ + BOOST_SPECIAL_FUNCTIONS_TEST_ASINH(type) \ + BOOST_SPECIAL_FUNCTIONS_TEST_ACOSH(type) \ + BOOST_SPECIAL_FUNCTIONS_TEST_SINC_PI(type) \ + BOOST_SPECIAL_FUNCTIONS_TEST_SINHC_PI(type) \ + BOOST_SPECIAL_FUNCTIONS_TEST_SINC_PI_FOR_C(type) \ + BOOST_SPECIAL_FUNCTIONS_TEST_SINHC_PI_FOR_C(type) +#endif BOOST_SPECIAL_FUNCTIONS_TEST(float) @@ -193,8 +243,12 @@ int test_main(int, char *[]) #undef BOOST_SPECIAL_FUNCTIONS_TEST #undef BOOST_SPECIAL_FUNCTIONS_TEST_ATANH +#undef BOOST_SPECIAL_FUNCTIONS_TEST_ASINH +#undef BOOST_SPECIAL_FUNCTIONS_TEST_ACOSH #undef BOOST_SPECIAL_FUNCTIONS_TEST_SINC_PI +#undef BOOST_SPECIAL_FUNCTIONS_TEST_SINC_PI_FOR_C #undef BOOST_SPECIAL_FUNCTIONS_TEST_SINHC_PI +#undef BOOST_SPECIAL_FUNCTIONS_TEST_SINHC_PI_FOR_C return(::boost::exit_success); }