From 46d9b3364520e103d4eff6fc2862c36ba702ee9e Mon Sep 17 00:00:00 2001 From: John Maddock Date: Tue, 1 Aug 2006 17:23:10 +0000 Subject: [PATCH] Tightened up root finding tests. Added docs for roots without derivatives. Added derivatives of incomplete gamma and beta. [SVN r3112] --- doc/beta_gamma_derivatives.qbk | 44 +++ doc/equations/derivative1.png | Bin 0 -> 1087 bytes doc/equations/derivative2.png | Bin 0 -> 1242 bytes doc/equations/derivatives.xml | 145 ++++++++ doc/math.qbk | 1 + doc/roots_without_derivatives.qbk | 310 ++++++++++++++++++ include/boost/math/special_functions/beta.hpp | 53 +++ .../boost/math/special_functions/gamma.hpp | 41 +++ test/test_ibeta.cpp | 10 + test/test_igamma.cpp | 7 + test/test_toms748_solve.cpp | 5 +- tools/log1p_expm1_data.cpp | 63 ++++ 12 files changed, 678 insertions(+), 1 deletion(-) create mode 100644 doc/beta_gamma_derivatives.qbk create mode 100644 doc/equations/derivative1.png create mode 100644 doc/equations/derivative2.png create mode 100644 doc/equations/derivatives.xml create mode 100644 doc/roots_without_derivatives.qbk create mode 100644 tools/log1p_expm1_data.cpp diff --git a/doc/beta_gamma_derivatives.qbk b/doc/beta_gamma_derivatives.qbk new file mode 100644 index 000000000..a1ec54472 --- /dev/null +++ b/doc/beta_gamma_derivatives.qbk @@ -0,0 +1,44 @@ +[#beta_gamma_derivatives][section Derivatives of the Incomplete Beta and Gamma Functions] + +[caution __caution ] + +[h4 Synopsis] + +`` +#include +`` + + namespace boost{ namespace math{ + + template + T gamma_P_derivative(T a, T x); + + template + T ibeta_derivative(T a, T b, T x); + + }} // namespaces + +[h4 Description] + +These two functions find some uses in statistical distributions, they +implement the partial derivative with respect to /x/ of the incomplete +gamma and beta functions respectively. + +[$../equations/derivative1.png] + +[$../equations/derivative2.png] + +[h4 Accuracy] + +Almost identical to the incomplete gamma and beta functions, refer to +the documentation for those function for more information. + +[h4 Implementation] + +These two functions just expose some of the internals of the incomplete +gamma and beta functions, refer to the documentation for those functions +for more information. + + +[endsect] + diff --git a/doc/equations/derivative1.png b/doc/equations/derivative1.png new file mode 100644 index 0000000000000000000000000000000000000000..abd91399a3fdaa85bee967572af14c4fc8a066d7 GIT binary patch literal 1087 zcmeAS@N?(olHy`uVBq!ia0y~yV4T6gz@Wv!#=yYP=jS`0fq{V~-O<;Pfnog#bJnhx z3=9mCC9V-A!TD(=<%vb942~)JNvR5+xryniL8*x;m4zo$Z5SAsFM7H-hE&{oGncn- zwSqvycIo$*e#@ty-{86aL3aGLBN7J~(-Ioqd;44tb4_=9=6>Y(zW4DA3=Fqwl$jV9 z4jg1)XyBIXKYl*_Ac(2;d!PNu1O|o$+|o)6AC7PRqR=UpD0})mN-^3UIV-*dNeZjgKF8|AHko7I0Z{wZ}>c#gr{xINkOe)#4qE7^BH zRI_im{9+nUnqAcEP21lFC@~k@EO~?9tV|9ZRm=La&>We{WnKW zmYL7(;)boyKF<7`lzCV20^`k-hHB9zJYhQ3{L9LI@$!pzeEp(mYsIw8b*6QL!*L4< zh69-lA)BANomdzsu<+94CuzawcC9q4PT}1c)EuIA{l3%0Sg&i(Go`kC(|o(cDR)|0 zoZ{c@QgM~yuII0KKA&>tf?>w8)d!!R^?y_SX8-x;8|VDz&Yiva-QABK({8L>(*vZKVje=uexnZl%%(&LrCo_dD-9BzB6#14vtZpJgqB=bNx2oY-jzu z|F`ZD)#iO>wd{(WZd&2ow~Myh#P@GIEmGS*%>-oM{P2PLdFly`6-GA@@T7SuBQlQ_zX2B93tum47 z3r6WIk+1F@zFRe^H2A)E+~1X78&cowmDP+qTGt&rW##9-spnK~roGeHG-tYmCzH_0 z7fyQTKPf0aKily_VwRPJzW<*?6@nX|X2%rry_mu9Xy2tEt^CO^R~Hy8NG&O``TG2J zGM7?`&A+V?zvgev>+(ubt}M%3xu8tmVrJnj$9pd-MAAkn)2+J(9blDXX4lc7Gp*4*@q4?l(mVcttG4F#l9D13=N-;oJ%TQ z#wse?P~`A2e#?n_izfcM>r?F%IzhBi+R$#r3kDJP4e~yB66_S4?RNI3-~O^u{aL!| zKXGez-u++SI=)~Kdm<+h!>}mzn~3H9^XcD7IyPWg(?)57=wBYY{wD7$DzXm_TzlAf>cn$l52vkJX!iQr zi6`NO5oa9_a5?Y%H&eGrQ`+y{bJ2CZhFwNZ>Y_ZRvIm{qv#fue;_JRKSMY&=b@cPT z6Q&%p0wRt74uF>3UJ#ba{(b0@m2pdJSXnZ7#xA~iqbK}O{ETf^Ugk?o z7udUnFUIPx;#-ro%%AMu&X~1rt>B3(hMX3UoBF?T1;3su^VcWT)Ajg+=u6Myci-C| zzb}1Wve$|9w1VATmHyG*KTn_8=a#adb$Q$3>GwqAqtD)5z48tpi&QXA)cmVPy}vo5 zbVD~@d7Qk$;J8}YT&3$z^lh>@wn+ppoEsPw&e8U7UUcEjO^@`qHmxZ;nsVaajV-As zp52&|cJk}ATcOjJ@<-`pp1<@@C-Ug@U-d^e=s#SU>91R_8L8@Pf4Inr{oCf>f>Wm* zDfHuRuv|DP4i*cgyT<7{qw;CyfU+aFe6c`CT}-^6*57ml0}J|?MCy_NrO zrQ7ez=hM|Du2!}1iDlwoXZU{VukN>c2BAfhSMN}nH(UI3Vs`Kh{W~lU>JN^HE&E^) zCh%eF_xs!L+jjSd-526p$Sk2#WE!cf#2$1t`nBE8wflA?O*(c(wL##HGso)ACz=W+ z48l#k*Pk6|(BJ1T7|_C(5b({Xl|8KH)fU}{HzZjd3~dS5`w +]> + + + + + + + + gamma_P_derivative + + + a + , + x + + + + = + + + ∂ + + ∂ + x + + + P + + + a + , + x + + + + = + + + + + e + + − + x + + + + x + + a + − + 1 + + + + + Γ + + + a + + + + + + + + + ibeta_derivative + + + a + , + b + , + x + + + + = + + + ∂ + + ∂ + x + + + + I + x + + + + a + , + b + + + + = + + + + + + + 1 + − + x + + + + b + − + 1 + + + + x + + a + − + 1 + + + + + B + + + a + , + b + + + + + + + + \ No newline at end of file diff --git a/doc/math.qbk b/doc/math.qbk index 2bb57864f..5fef09bcf 100644 --- a/doc/math.qbk +++ b/doc/math.qbk @@ -81,6 +81,7 @@ data and/or data for output to an external graphing application. [include beta.qbk] [include ibeta.qbk] [include ibeta_inv.qbk] +[include beta_gamma_derivatives.qbk] [include fpclassify.qbk] [include error_handling.qbk] diff --git a/doc/roots_without_derivatives.qbk b/doc/roots_without_derivatives.qbk new file mode 100644 index 000000000..3866206a1 --- /dev/null +++ b/doc/roots_without_derivatives.qbk @@ -0,0 +1,310 @@ +[#roots2][section Root Finding Without Derivatives] + +[caution __caution ] + +[h4 Synopsis] + +`` +#include +`` + + namespace boost{ namespace math{ namespace tools{ + + template + std::pair + bisect( + F f, + T min, + T max, + Tol tol, + boost::uintmax_t& max_iter); + + template + std::pair + bisect( + F f, + T min, + T max, + Tol tol); + + template + std::pair + bracket_and_solve_root( + F f, + const T& guess, + const T& factor, + bool rising, + Tol tol, + boost::uintmax_t& max_iter); + + template + std::pair + toms748_solve( + F f, + const T& a, + const T& b, + Tol tol, + boost::uintmax_t& max_iter); + + template + std::pair + toms748_solve( + F f, + const T& a, + const T& b, + const T& fa, + const T& fb, + Tol tol, + boost::uintmax_t& max_iter); + + // Termination conditions: + template + struct eps_tolerance; + + struct equal_floor; + struct equal_ceil; + + }}} // namespaces + +[h4 Description] + +These functions solve the root of some function /f(x)/ without the +need for the derivatives of /f(x)/. The functions here that use TOMS +Algorithm 748 are asymptotically the most efficient known, and have +been shown to be optimal for a certain classes of smooth functions. + +Alternatively there is a simple bisection routine which can be useful +in it's own right in some situations, or alternatively for narrowing +down the range containing the root prior to calling a more advanced +algorithm. + +All the algorithms in this section reduce the diameter of the enclosing +interval with the same asymptotic efficiency with which they locate the +root. This is in contrast to the derivative based methods which may /never/ +significantly reduce the enclosing interval, even though they rapidly approach +the root. This is also in contrast to some other derivate free methods +(for example the methods of Brent or Dekker) which only reduce the enclosing +interval on the final step. Therefore these methods return a std::pair +containing the enclosing +interval found, and accept a function object specifying the termination condition. +Three function objects are provided for ready-made termination conditions: +/eps_tolerance/ causes termination when the relative error in the enclosing +interval is below a certain threshold, while /equal_floor/ and /equal_ceil/ are +useful for certain statistical applications where the result is known to be +an integer. Other user defined termination conditions are likely to be used +only rairly, but may be useful in some specific circumstances. + + template + std::pair + bisect( + F f, + T min, + T max, + Tol tol, + boost::uintmax_t& max_iter); + + template + std::pair + bisect( + F f, + T min, + T max, + Tol tol); + +These two functions locate the root using bisection, function arguments are: + +[variablelist +[[f] [A unary functor which the function whose root is to be found.]] +[[min] [The left bracket of the interval known to contain the root.]] +[[max] [The right bracket of the interval known to contain the root. + It is a precondition that /min < max/ and /f(min)*f(max) <= 0/, + the function calls __logic_error if these preconditions are violated.]] +[[tol] [A binary functor that specifies the termination condition: the function + will return the current brackets enclosing the root when /tol(min,max)/ becomes true.]] +[[max_iter][The maximum number of invocations of /f(x)/ to make while searching for the root.]] +] + +Returns: a pair of values /r/ that bracket the root so that: + + f(r.first) * f(r.second) <= 0 + +and either + + tol(r.first, r.second) == true + +or + + max_iter >= m + +where /m/ is the initial value of /max_iter/ passed to the function. + +In other words it's up to the caller to verify whether termination occured +as a result of exceeding /max_iter/ function invocations (easily done by +checking the value of /max_iter/), rather than because the termination +condition /tol/ was satisfied. + + template + std::pair + bracket_and_solve_root( + F f, + const T& guess, + const T& factor, + bool rising, + Tol tol, + boost::uintmax_t& max_iter); + +This is a convenience function that calls /toms748_solve/ internally +to find the root of /f(x)/. It's usable only when /f(x)/ is a monotonic +function, and the location of the root is known approximately, and in +particular it is known whether the root is occurs for positive or negative +/x/. The parameters are: + +[variablelist +[[f][A unary functor that is the function whose root is to be solved. + f(x) must be uniformly increasing or decreasing on /x/.]] +[[guess][An initial approximation to the root]] +[[factor][A scaling factor that is used to bracket the root: the value + /guess/ is mutiplied (or divided as appropriate) by /factor/ + until two values are found that bracket the root. A value + such as 2 is a typical choice for /factor/.]] +[[rising][Set to /true/ if /f(x)/ is rising on /x/ and /false/ if /f(x)/ + is falling on /x/. This value is used along with the result + of /f(guess)/ to determine if /guess/ is + above or below the root.]] +[[tol] [A binary functor that determines the termination condition for the search + for the root. /tol/ is passed the current brackets at each step, + when it returns true then the current brackets are returned as the result.]] +[[max_iter] [The maximum number of function invocations to perform in the search + for the root.]] +] + +Returns: a pair of values /r/ that bracket the root so that: + + f(r.first) * f(r.second) <= 0 + +and either + + tol(r.first, r.second) == true + +or + + max_iter >= m + +where /m/ is the initial value of /max_iter/ passed to the function. + +In other words it's up to the caller to verify whether termination occured +as a result of exceeding /max_iter/ function invocations (easily done by +checking the value of /max_iter/), rather than because the termination +condition /tol/ was satisfied. + + template + std::pair + toms748_solve( + F f, + const T& a, + const T& b, + Tol tol, + boost::uintmax_t& max_iter); + + template + std::pair + toms748_solve( + F f, + const T& a, + const T& b, + const T& fa, + const T& fb, + Tol tol, + boost::uintmax_t& max_iter); + +These two functions implement TOMS Algorithm 748: it uses a mixture of +cubic, quadratic and linear (secant) interpolation to locate the root of +/f(x)/. The two functions differ only by whether values for /f(a)/ and +/f(b)/ are already available. The parameters are: + +[variablelist +[[f] [A unary functor that is the function whose root is to be solved. + f(x) need not be uniformly increasing or decreasing on /x/ and + may have mutilple roots.]] +[[a] [ The lower bound for the initial bracket of the root.]] +[[b] [The upper bound for the initial bracket of the root. + It is a precondition that /a < b/ and that /a/ and /b/ + bracket the root to find so that /f(a)*f(b) < 0/.]] +[[fa] [Optional: the value of /f(a)/.]] +[[fb] [Optional: the value of /f(b)/.]] +[[tol] [A binary functor that determines the termination condition for the search + for the root. /tol/ is passed the current brackets at each step, + when it returns true then the current brackets are returned as the result.]] +[[max_iter] [The maximum number of function invocations to perform in the search + for the root. On exit /max_iter/ is set to actual number of function + invocations used.]] +] + +Returns: a pair of values /r/ that bracket the root so that: + + f(r.first) * f(r.second) <= 0 + +and either + + tol(r.first, r.second) == true + +or + + max_iter >= m + +where /m/ is the initial value of /max_iter/ passed to the function. + +In other words it's up to the caller to verify whether termination occured +as a result of exceeding /max_iter/ function invocations (easily done by +checking the value of /max_iter/), rather than because the termination +condition /tol/ was satisfied. + + template + struct eps_tolerance + { + eps_tolerance(int bits); + bool operator()(const T& a, const T& b)const; + }; + +This is the usual termination condition used with these root finding functions. +It's operator() will return true when the relative distance between /a/ and /b/ +is less than twice the machine epsilon for T, or 2[super 1-bits], whichever is +the larger. In other words you set /bits/ to the number of bits of precision you +want in the result. The minimal tolerance of twice the machine epsilon of T is +required to ensure that we get back a bracketing interval: since this must clearly +be at least 1 epsilon in size. + + struct equal_floor + { + equal_floor(); + template bool operator()(const T& a, const T& b)const; + }; + +This termination condition is used when you want to find an integer result +that is the /floor/ of the true root. It will terminate as soon as both ends +of the interval have the same /floor/. + + struct equal_ceil + { + equal_ceil(); + template bool operator()(const T& a, const T& b)const; + }; + +This termination condition is used when you want to find an integer result +that is the /ceil/ of the true root. It will terminate as soon as both ends +of the interval have the same /ceil/. + +[h4 Implementation] + +The implementation of the bisection algorithm is extrememly straightforward +and not detailed here. TOMS algorithm 748 is described in detail in: + +['Algorithm 748: Enclosing Zeros of Continuous Functions, +G. E. Alefeld, F. A. Potra and Yixun Shi, +ACM Transactions on Mathematica1 Software, Vol. 21. No. 3. September 1995. +Pages 327-344.] + +The implementation here is a faithful translation of this paper into C++. + +[endsect] + diff --git a/include/boost/math/special_functions/beta.hpp b/include/boost/math/special_functions/beta.hpp index 7ce1532b3..77be924f7 100644 --- a/include/boost/math/special_functions/beta.hpp +++ b/include/boost/math/special_functions/beta.hpp @@ -926,6 +926,51 @@ T ibeta_imp(T a, T b, T x, const L& l, bool inv, bool normalised) return invert ? (normalised ? 1 : beta_imp(a, b, l)) - fract : fract; } // template T ibeta_imp(T a, T b, T x, const L& l, bool inv, bool normalised) +template +T ibeta_derivative_imp(T a, T b, T x, const L& l) +{ + // + // start with the usual error checks: + // + if(a <= 0) + tools::domain_error(BOOST_CURRENT_FUNCTION, "The argument a to the incomplete beta function must be greater than zero (got a=%1%).", a); + if(b <= 0) + tools::domain_error(BOOST_CURRENT_FUNCTION, "The argument b to the incomplete beta function must be greater than zero (got b=%1%).", b); + if((x < 0) || (x > 1)) + tools::domain_error(BOOST_CURRENT_FUNCTION, "Parameter x outside the range [0,1] in the incomplete beta function (got x=%1%).", x); + // + // Now the corner cases: + // + if(x == 0) + { + return (a > 1) ? 0 : + (a == 1) ? 1 : tools::overflow_error(BOOST_CURRENT_FUNCTION, 0); + } + else if(x == 1) + { + return (b > 1) ? 0 : + (b == 1) ? 1 : tools::overflow_error(BOOST_CURRENT_FUNCTION, 0); + } + // + // Now the regular cases: + // + T f1 = ibeta_power_terms(a, b, x, 1 - x, l, true); + T y = (1 - x) * x; + + if(f1 == 0) + return 0; + + if((tools::max_value() * y < f1)) + { + // overflow: + return tools::overflow_error(BOOST_CURRENT_FUNCTION, 0); + } + + f1 /= y; + + return f1; +} + } // namespace detail // @@ -973,6 +1018,14 @@ T ibetac(T a, T b, T x) return tools::checked_narrowing_cast(detail::ibeta_imp(static_cast(a), static_cast(b), static_cast(x), evaluation_type(), true, true), BOOST_CURRENT_FUNCTION); } +template +T ibeta_derivative(T a, T b, T x) +{ + typedef typename lanczos::lanczos_traits::value_type value_type; + typedef typename lanczos::lanczos_traits::evaluation_type evaluation_type; + return tools::checked_narrowing_cast(detail::ibeta_derivative_imp(static_cast(a), static_cast(b), static_cast(x), evaluation_type()), BOOST_CURRENT_FUNCTION); +} + } // namespace math } // namespace boost diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index bf5f9e4e1..3c29a8ce9 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -873,6 +873,39 @@ T tgamma_delta_ratio_imp(T z, T delta, const lanczos::undefined_lanczos&) return sum * prefix; } +template +T gamma_P_derivative_imp(T a, T x, const L& l) +{ + // + // Usual error checks first: + // + if(a <= 0) + tools::domain_error(BOOST_CURRENT_FUNCTION, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a); + if(x < 0) + tools::domain_error(BOOST_CURRENT_FUNCTION, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x); + // + // Now special cases: + // + if(x == 0) + { + return (a > 1) ? 0 : + (a == 1) ? 1 : tools::overflow_error(BOOST_CURRENT_FUNCTION, 0); + } + // + // Normal case: + // + T f1 = detail::regularised_gamma_prefix(a, x, l); + if((x < 1) && (tools::max_value() * x < f1)) + { + // overflow: + return tools::overflow_error(BOOST_CURRENT_FUNCTION, 0); + } + + f1 /= x; + + return f1; +} + } // namespace detail template @@ -974,6 +1007,14 @@ inline T tgamma_ratio(T a, T b) return tools::checked_narrowing_cast(detail::tgamma_delta_ratio_imp(static_cast(a), static_cast(b - a), evaluation_type()), BOOST_CURRENT_FUNCTION); } +template +T gamma_P_derivative(T a, T x) +{ + typedef typename lanczos::lanczos_traits::value_type value_type; + typedef typename lanczos::lanczos_traits::evaluation_type evaluation_type; + return tools::checked_narrowing_cast(detail::gamma_P_derivative_imp(static_cast(a), static_cast(x), evaluation_type()), BOOST_CURRENT_FUNCTION); +} + } // namespace math } // namespace boost diff --git a/test/test_ibeta.cpp b/test/test_ibeta.cpp index 4ade0ea08..418779741 100644 --- a/test/test_ibeta.cpp +++ b/test/test_ibeta.cpp @@ -319,6 +319,16 @@ void test_spots(T) static_cast(1), static_cast(0.32)), static_cast(0.948565954109602496577407403168592262389L), tolerance); + + // very naive check on derivative: + using namespace std; // For ADL of std functions + tolerance = boost::math::tools::epsilon() * 10000; // 100 eps + BOOST_CHECK_CLOSE( + ::boost::math::ibeta_derivative( + static_cast(2), + static_cast(3), + static_cast(0.5)), + pow(static_cast(0.5), static_cast(2)) * pow(static_cast(0.5), static_cast(1)) / boost::math::beta(static_cast(2), static_cast(3)), tolerance); } int test_main(int, char* []) diff --git a/test/test_igamma.cpp b/test/test_igamma.cpp index 0632b9b8c..ec999cff2 100644 --- a/test/test_igamma.cpp +++ b/test/test_igamma.cpp @@ -275,6 +275,13 @@ void test_spots(T) BOOST_CHECK_CLOSE(::boost::math::gamma_P(static_cast(5), static_cast(100)), static_cast(0.9999999999999999999999999999999999998386069466302L), tolerance); BOOST_CHECK_CLOSE(::boost::math::gamma_P(static_cast(1.5), static_cast(2)), static_cast(0.73853587005088937779717792402407879809718939080921L), tolerance); BOOST_CHECK_CLOSE(::boost::math::gamma_P(static_cast(20.5), static_cast(22)), static_cast(0.65424667956532673185028409120341593367429721070928L), tolerance); + + // naive check on derivative function: + using namespace std; // For ADL of std functions + tolerance = boost::math::tools::epsilon() * 5000; // 50 eps + BOOST_CHECK_CLOSE(::boost::math::gamma_P_derivative(static_cast(20.5), static_cast(22)), + exp(static_cast(-22)) * pow(static_cast(22), static_cast(19.5)) / boost::math::tgamma(static_cast(20.5)), tolerance); + } int test_main(int, char* []) diff --git a/test/test_toms748_solve.cpp b/test/test_toms748_solve.cpp index 973c0c667..98f225064 100644 --- a/test/test_toms748_solve.cpp +++ b/test/test_toms748_solve.cpp @@ -120,6 +120,7 @@ void run_test(T a, T b, int id) b, boost::math::tools::eps_tolerance(std::numeric_limits::digits), c); + BOOST_CHECK_EQUAL(c, toms748tester::total_calls()); total += c; invocations += toms748tester::total_calls(); std::cout << "Function " << id << "\nresult={" << r.first << ", " << r.second << "} total calls=" << toms748tester::total_calls() << "\n\n"; @@ -135,6 +136,7 @@ void run_test(T a, T b, int id, int p) b, boost::math::tools::eps_tolerance(std::numeric_limits::digits), c); + BOOST_CHECK_EQUAL(c, toms748tester::total_calls()); total += c; invocations += toms748tester::total_calls(); std::cout << "Function " << id << "\nresult={" << r.first << ", " << r.second << "} total calls=" << toms748tester::total_calls() << "\n\n"; @@ -150,6 +152,7 @@ void run_test(T a, T b, int id, T p1, T p2) b, boost::math::tools::eps_tolerance(std::numeric_limits::digits), c); + BOOST_CHECK_EQUAL(c, toms748tester::total_calls()); total += c; invocations += toms748tester::total_calls(); std::cout << "Function " << id << "\n Result={" << r.first << ", " << r.second << "} total calls=" << toms748tester::total_calls() << "\n\n"; @@ -259,7 +262,7 @@ int test_main(int, char* []) std::cout << std::setprecision(18); std::cout << "Function " << 4 << "\n Result={" << r.first << ", " << r.second << "} total calls=" << toms748tester::total_calls() << "\n\n"; toms748tester::reset(); - BOOST_CHECK(c < 5); + BOOST_CHECK(c < 20); } return 0; diff --git a/tools/log1p_expm1_data.cpp b/tools/log1p_expm1_data.cpp new file mode 100644 index 000000000..7b2ab24cc --- /dev/null +++ b/tools/log1p_expm1_data.cpp @@ -0,0 +1,63 @@ +// (C) Copyright John Maddock 2006. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include + +#include +#include +#include + +using namespace boost::math::tools; +using namespace std; + +struct data_generator +{ + std::tr1::tuple operator()(NTL::RR z) + { + return std::tr1::make_tuple(boost::math::log1p(z), boost::math::expm1(z)); + } +}; + +int main(int argc, char* argv[]) +{ + NTL::RR::SetPrecision(1000); + NTL::RR::SetOutputPrecision(40); + + parameter_info arg1; + test_data data; + + std::cout << "Welcome.\n" + "This program will generate spot tests for the log1p and expm1 functions:\n\n"; + + bool cont; + std::string line; + + do{ + if(0 == get_user_parameter_info(arg1, "z")) + return 1; + data.insert(data_generator(), arg1); + + std::cout << "Any more data [y/n]?"; + std::getline(std::cin, line); + boost::algorithm::trim(line); + cont = (line == "y"); + }while(cont); + + std::cout << "Enter name of test data file [default=log1p_expm1_data.ipp]"; + std::getline(std::cin, line); + boost::algorithm::trim(line); + if(line == "") + line = "log1p_expm1_data.ipp"; + std::ofstream ofs(line.c_str()); + ofs << std::scientific; + write_code(ofs, data, "log1p_expm1_data"); + + return 0; +} + + + +