From 27f5e4d59ea4137fb2541dfe231c7ced5029c8a6 Mon Sep 17 00:00:00 2001 From: John Maddock Date: Thu, 4 Jan 2007 16:20:49 +0000 Subject: [PATCH] Updated hypot with simpler formula. Added sinh and cosh to ntl.hpp. Fixed FreeBSD test failures (no real long long support). [SVN r3598] --- doc/equations/hypot.png | Bin 283 -> 426 bytes doc/equations/hypot2.png | Bin 415 -> 532 bytes doc/powers.qbk | 11 ++--- .../boost/math/special_functions/hypot.hpp | 45 +++--------------- include/boost/math/tools/ntl.hpp | 10 ++++ test/test_binomial.cpp | 7 +++ test/test_carlson.cpp | 7 +++ test/test_ellint_1.cpp | 7 +++ test/test_ellint_2.cpp | 7 +++ test/test_ellint_3.cpp | 7 +++ test/test_ibeta.cpp | 8 ++-- test/test_laguerre.cpp | 2 + test/test_legendre.cpp | 2 + test/test_negative_binomial.cpp | 7 +++ test/test_spherical_harmonic.cpp | 2 + test/test_tgamma_ratio.cpp | 2 +- 16 files changed, 73 insertions(+), 51 deletions(-) diff --git a/doc/equations/hypot.png b/doc/equations/hypot.png index fa88d6b216613ba0ef0ee28f7b0becf74049020e..1f46aa80f580c36d0ad3baa9d42480a6e4869631 100644 GIT binary patch delta 376 zcmbQuw2E1=Gr-TCmrII^fq{Y7)59f*fq~&90|SFN2Qvc$!!(`Q850$4>uqdo{{R2a z;p;V(fq{XkB*-uLKf}}Q23`ydj53}sjv*Dd-puIcJFLLr%K!iWWpDZ5K*5beHqg(8i^CHoJ5t z*GhXn3Jli`c;i<*sU@ydQR(SL#e+(YY${rlPhAnrJfM}|Hlg6jM?TKj!gsFPk*l+& zi}+b?ow_4y;fl%rWyQCS#ZTP*eMeu}=C{W_vunRoNq8|g?k7L<$H@v2>KU&U&VGEy lC_B&i#Q934kAevd9`*cq&ffq delta 232 zcmZ3*Jex_eGr-TCmrII^fq{Y7)59f*fq}t=fq_AigPDPWVgKEVl@k?h>lqjr{{R2K z^yd7#3=9lRB|(0{{~4ZcH}GO$VCeO9aSW-rb>_-zt_B4jhlh`U?=KNie79)!vI+HE zG4n6$YD%;HcbOq4w@Gy}Z$|~cLx*6~mJ1i>-kos&aAeEY!y#0|SE=2Qvc$!_3lhn~93H^)@y(|NsBz z@b#L?z`(#%666>BpW*3t11|;!#x_qE$B>F!Z)P<19ai9Iv%mL$seZJ1cBn}DU9IbC z(R$Ng{a~nfIP&*!?~mf&$L6UD{C5<&tE^<6^j5&yJth8)W#DuNm%JliM7iqM_w5iZ zKNOz0miq{YKHKMsD>57w*6y}gG_&5o>5TYJ$1DB!E_Jx-Mni|;=;Jzg2-`T(|t+UpP(3uGlU0 zC~U&Cbbb(49Qp?}Zz+v}FUJ}7*cl}U7B za=q086;IKROa0Q7w+U|*S!o$)q7ZaLQgYu0mF=}cNsDAxJXl(MBw*IEZOY4)_ANDe zQx_0&Gk>c2U0>JqITK>?KOMFdbG_R5YWJBv${Rbcr*}kLf0le^=G6Nu&d>Pg{nDqg z_Pplits76&HB>CsJY}n7c6h~>`f&Nq-va$bZ@4!WSgNVWH?`rHk?<+TR$hr7RWbhgc@Z?l-dd2(yXr~He3$pv2W zCrV4dpL1BTmfP@tm^IhlG}gjHIhzfHpIpx8ZkOg%2Zak`)V!3;I`3)Zj>nnZ=7Q%P#1^nSu|5yA;Fo9v0 YiooraS!eGsFfcH9y85}Sb4q9e0HCXy4*&oF diff --git a/doc/powers.qbk b/doc/powers.qbk index 252826783..ad756990a 100644 --- a/doc/powers.qbk +++ b/doc/powers.qbk @@ -193,16 +193,11 @@ overflow or underflow, even though the result is in fact perfectly representable The function is even and symmetric in x and y, so first take assume ['x,y > 0] and ['x > y] (we can permute the arguments if this is not the case). -Then if ['x > HIGH] return [$../equations/hypot2.png] for ['y > 1] or /x/ otherwise. +Then if ['x * [epsilon][space] >= y] we can simply return /x/. -Otherwise if ['y < LOW] return x if ['x > 1 || y == 0] otherwise [$../equations/hypot2.png]. +Otherwise the result is given by: -All other cases are free from overflow or underflow problems and so can be dealt -with via the usual formula [$../equations/hypot.png]. - -The constant HIGH is the square root of half the largest -representable value. The constant LOW is the square root of the -smallest representable value. +[$../equations/hypot2.png] [endsect] diff --git a/include/boost/math/special_functions/hypot.hpp b/include/boost/math/special_functions/hypot.hpp index cc7577e0e..d251e98fd 100644 --- a/include/boost/math/special_functions/hypot.hpp +++ b/include/boost/math/special_functions/hypot.hpp @@ -43,45 +43,12 @@ T hypot(T x, T y) if(y > x) (std::swap)(x, y); - // - // Figure out overflow and underflow limits, - // we could make these constants static to save - // a few cycles, but the code would then not be - // thread safe :-( - // - T safe_upper = sqrt(tools::max_value()) / 2; - T safe_lower = sqrt(tools::min_value()); - static const T one = 1; - // - // Now handle special cases: - // - if(x >= safe_upper) - { - if(y <= one) - { - // y is negligible: - return x; - } - T a = sqrt(x) * sqrt(y); - T b = sqrt(x/y + y/x); - // We may have overflow: - if(tools::max_value() /a < b) - return tools::overflow_error(BOOST_CURRENT_FUNCTION, 0); - return a * b; - } - else if(y <= safe_lower) - { - if((x >= one) || (y == 0)) - { - // y is negligible: - return x; - } - return sqrt(x) * sqrt(y) * sqrt(x/y + y/x); - } - // - // If we get here then x^2+y^2 will not overflow or underflow: - // - return sqrt(x*x + y*y); + + if(x * tools::epsilon() >= y) + return x; + + T rat = y / x; + return x * sqrt(1 + rat*rat); } // template T hypot(T x, T y) diff --git a/include/boost/math/tools/ntl.hpp b/include/boost/math/tools/ntl.hpp index d57e37072..601d26da6 100644 --- a/include/boost/math/tools/ntl.hpp +++ b/include/boost/math/tools/ntl.hpp @@ -325,6 +325,16 @@ namespace NTL{ boost::math::tools::digits()); } + inline NTL::RR sinh(NTL::RR z) + { + return (expm1(z) - expm1(-z)) / 2; + } + + inline NTL::RR cosh(NTL::RR z) + { + return (exp(z) + exp(-z)) / 2; + } + inline NTL::RR fmod(NTL::RR x, NTL::RR y) { // This is a really crummy version of fmod, we rely on lots diff --git a/test/test_binomial.cpp b/test/test_binomial.cpp index df90b5357..df577b2d1 100644 --- a/test/test_binomial.cpp +++ b/test/test_binomial.cpp @@ -640,10 +640,17 @@ int test_main(int, char* []) // (Parameter value, arbitrarily zero, only communicates the floating point type). test_spots(0.0F); // Test float. test_spots(0.0); // Test double. +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS test_spots(0.0L); // Test long double. #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) test_spots(boost::math::concepts::real_concept(0.)); // Test real concept. #endif +#else + std::cout << "The long double tests have been disabled on this platform " + "either because the long double overloads of the usual math functions are " + "not available at all, or because they are too inaccurate for these tests " + "to pass." << std::cout; +#endif return 0; } // int test_main(int, char* []) diff --git a/test/test_carlson.cpp b/test/test_carlson.cpp index 056e09a7f..f35be894d 100644 --- a/test/test_carlson.cpp +++ b/test/test_carlson.cpp @@ -321,8 +321,15 @@ int test_main(int, char* []) expected_results(); test_spots(0.0F, "float"); test_spots(0.0, "double"); +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS test_spots(0.0L, "long double"); test_spots(boost::math::concepts::real_concept(0), "real_concept"); +#else + std::cout << "The long double tests have been disabled on this platform " + "either because the long double overloads of the usual math functions are " + "not available at all, or because they are too inaccurate for these tests " + "to pass." << std::cout; +#endif return 0; } diff --git a/test/test_ellint_1.cpp b/test/test_ellint_1.cpp index 1f62d6ccb..cc3afbcf9 100644 --- a/test/test_ellint_1.cpp +++ b/test/test_ellint_1.cpp @@ -193,8 +193,15 @@ int test_main(int, char* []) test_spots(0.0F, "float"); test_spots(0.0, "double"); +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS test_spots(0.0L, "long double"); test_spots(boost::math::concepts::real_concept(0), "real_concept"); +#else + std::cout << "The long double tests have been disabled on this platform " + "either because the long double overloads of the usual math functions are " + "not available at all, or because they are too inaccurate for these tests " + "to pass." << std::cout; +#endif return 0; } diff --git a/test/test_ellint_2.cpp b/test/test_ellint_2.cpp index 741cfaac3..837b7b2b4 100644 --- a/test/test_ellint_2.cpp +++ b/test/test_ellint_2.cpp @@ -182,8 +182,15 @@ int test_main(int, char* []) expected_results(); test_spots(0.0F, "float"); test_spots(0.0, "double"); +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS test_spots(0.0L, "long double"); test_spots(boost::math::concepts::real_concept(0), "real_concept"); +#else + std::cout << "The long double tests have been disabled on this platform " + "either because the long double overloads of the usual math functions are " + "not available at all, or because they are too inaccurate for these tests " + "to pass." << std::cout; +#endif return 0; } diff --git a/test/test_ellint_3.cpp b/test/test_ellint_3.cpp index 98eb2ef8b..67ccdf677 100644 --- a/test/test_ellint_3.cpp +++ b/test/test_ellint_3.cpp @@ -224,8 +224,15 @@ int test_main(int, char* []) expected_results(); test_spots(0.0F, "float"); test_spots(0.0, "double"); +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS test_spots(0.0L, "long double"); test_spots(boost::math::concepts::real_concept(0), "real_concept"); +#else + std::cout << "The long double tests have been disabled on this platform " + "either because the long double overloads of the usual math functions are " + "not available at all, or because they are too inaccurate for these tests " + "to pass." << std::cout; +#endif return 0; } diff --git a/test/test_ibeta.cpp b/test/test_ibeta.cpp index c54d0fd26..10b92eb0b 100644 --- a/test/test_ibeta.cpp +++ b/test/test_ibeta.cpp @@ -122,16 +122,18 @@ void expected_results() largest_type, // test type(s) "(?i).*large.*", // test data group ".*", 200000, 10000); // test function +#ifdef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS // - // Cygwin: + // No long doubles: // add_expected_result( "[^|]*", // compiler "[^|]*", // stdlib - "cygwin", // platform + BOOST_PLATFORM, // platform largest_type, // test type(s) "(?i).*large.*", // test data group - ".*", 10000, 500); // test function + ".*", 13000, 500); // test function +#endif // // Catch all cases come last: // diff --git a/test/test_laguerre.cpp b/test/test_laguerre.cpp index 4838a051f..4299c59e4 100644 --- a/test/test_laguerre.cpp +++ b/test/test_laguerre.cpp @@ -66,6 +66,7 @@ void expected_results() if((std::numeric_limits::digits <= 64) && (std::numeric_limits::digits != std::numeric_limits::digits)) { +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS add_expected_result( ".*", // compiler ".*", // stdlib @@ -73,6 +74,7 @@ void expected_results() "double", // test type(s) ".*", // test data group ".*", 10, 5); // test function +#endif } add_expected_result( ".*", // compiler diff --git a/test/test_legendre.cpp b/test/test_legendre.cpp index 48e3f533b..64115d113 100644 --- a/test/test_legendre.cpp +++ b/test/test_legendre.cpp @@ -64,6 +64,7 @@ void expected_results() if((std::numeric_limits::digits <= 64) && (std::numeric_limits::digits != std::numeric_limits::digits)) { +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS add_expected_result( ".*", // compiler ".*", // stdlib @@ -71,6 +72,7 @@ void expected_results() "double", // test type(s) ".*", // test data group ".*", 10, 5); // test function +#endif } add_expected_result( ".*", // compiler diff --git a/test/test_negative_binomial.cpp b/test/test_negative_binomial.cpp index 225d7e67c..8442138e3 100644 --- a/test/test_negative_binomial.cpp +++ b/test/test_negative_binomial.cpp @@ -758,10 +758,17 @@ int test_main(int, char* []) // (Parameter value, arbitrarily zero, only communicates the floating point type). test_spots(0.0F); // Test float. test_spots(0.0); // Test double. +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS test_spots(0.0L); // Test long double. #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) test_spots(boost::math::concepts::real_concept(0.)); // Test real concept. #endif +#else + std::cout << "The long double tests have been disabled on this platform " + "either because the long double overloads of the usual math functions are " + "not available at all, or because they are too inaccurate for these tests " + "to pass." << std::cout; +#endif return 0; } // int test_main(int, char* []) diff --git a/test/test_spherical_harmonic.cpp b/test/test_spherical_harmonic.cpp index baccb12b6..7f838999b 100644 --- a/test/test_spherical_harmonic.cpp +++ b/test/test_spherical_harmonic.cpp @@ -57,6 +57,7 @@ void expected_results() #else largest_type = "(long\\s+)?double"; #endif +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS if((std::numeric_limits::digits <= 64) && (std::numeric_limits::digits != std::numeric_limits::digits)) { @@ -68,6 +69,7 @@ void expected_results() ".*", // test data group ".*", 10, 5); // test function } +#endif // // Catch all cases come last: // diff --git a/test/test_tgamma_ratio.cpp b/test/test_tgamma_ratio.cpp index f96ac3535..0b81b896f 100644 --- a/test/test_tgamma_ratio.cpp +++ b/test/test_tgamma_ratio.cpp @@ -63,7 +63,7 @@ void expected_results() add_expected_result( "[^|]*", // compiler "[^|]*", // stdlib - "HP-UX|linux", // platform + "HP-UX|linux|.*(bsd|BSD).*", // platform "float", // test type(s) "[^|]*", // test data group "boost::math::tgamma_ratio[^|]*", 35, 8); // test function