1 [section:special_tut Tutorial: How to Write a New Special Function]
3 [section:special_tut_impl Implementation]
5 In this section, we'll provide a "recipe" for adding a new special function to this library to make life easier for
6 future authors wishing to contribute. We'll assume the function returns a single floating-point result, and takes
7 two floating-point arguments. For the sake of exposition we'll give the function the name [~my_special].
9 Normally, the implementation of such a function is split into two layers - a public user layer, and an internal
10 implementation layer that does the actual work.
11 The implementation layer is declared inside a `detail` namespace and has a simple signature:
13 namespace boost { namespace math { namespace detail {
15 template <class T, class Policy>
16 T my_special_imp(const T& a, const T&b, const Policy& pol)
18 /* Implementation goes here */
23 We'll come back to what can go inside the implementation later, but first lets look at the user layer.
24 This consists of two overloads of the function, with and without a __Policy argument:
26 namespace boost{ namespace math{
28 template <class T, class U>
29 typename tools::promote_args<T, U>::type my_special(const T& a, const U& b);
31 template <class T, class U, class Policy>
32 typename tools::promote_args<T, U>::type my_special(const T& a, const U& b, const Policy& pol);
36 Note how each argument has a different template type - this allows for mixed type arguments - the return
37 type is computed from a traits class and is the "common type" of all the arguments after any integer
38 arguments have been promoted to type `double`.
40 The implementation of the non-policy overload is trivial:
42 namespace boost{ namespace math{
44 template <class T, class U>
45 inline typename tools::promote_args<T, U>::type my_special(const T& a, const U& b)
47 // Simply forward with a default policy:
48 return my_special(a, b, policies::policy<>();
53 The implementation of the other overload is somewhat more complex, as there's some meta-programming to do,
54 but from a runtime perspective is still a one-line forwarding function. Here it is with comments explaining
57 namespace boost{ namespace math{
59 template <class T, class U, class Policy>
60 inline typename tools::promote_args<T, U>::type my_special(const T& a, const U& b, const Policy& pol)
63 // We've found some standard library functions to misbehave if any FPU exception flags
64 // are set prior to their call, this code will clear those flags, then reset them
67 BOOST_FPU_EXCEPTION_GUARD
69 // The type of the result - the common type of T and U after
70 // any integer types have been promoted to double:
72 typedef typename tools::promote_args<T, U>::type result_type;
74 // The type used for the calculation. This may be a wider type than
75 // the result in order to ensure full precision:
77 typedef typename policies::evaluation<result_type, Policy>::type value_type;
79 // The type of the policy to forward to the actual implementation.
80 // We disable promotion of float and double as that's [possibly]
81 // happened already in the line above. Also reset to the default
82 // any policies we don't use (reduces code bloat if we're called
83 // multiple times with differing policies we don't actually use).
84 // Also normalise the type, again to reduce code bloat in case we're
85 // called multiple times with functionally identical policies that happen
86 // to be different types.
88 typedef typename policies::normalise<
90 policies::promote_float<false>,
91 policies::promote_double<false>,
92 policies::discrete_quantile<>,
93 policies::assert_undefined<> >::type forwarding_policy;
95 // Whew. Now we can make the actual call to the implementation.
96 // Arguments are explicitly cast to the evaluation type, and the result
97 // passed through checked_narrowing_cast which handles things like overflow
98 // according to the policy passed:
100 return policies::checked_narrowing_cast<result_type, forwarding_policy>(
101 detail::my_special_imp(
102 static_cast<value_type>(a),
103 static_cast<value_type>(x),
104 forwarding_policy()),
105 "boost::math::my_special<%1%>(%1%, %1%)");
110 We're now almost there, we just need to flesh out the details of the implementation layer:
112 namespace boost { namespace math { namespace detail {
114 template <class T, class Policy>
115 T my_special_imp(const T& a, const T&b, const Policy& pol)
117 /* Implementation goes here */
122 The following guidelines indicate what (other than basic arithmetic) can go in the implementation:
124 * Error conditions (for example bad arguments) should be handled by calling one of the
125 [link math_toolkit.error_handling.finding_more_information policy based error handlers].
126 * Calls to standard library functions should be made unqualified (this allows argument
127 dependent lookup to find standard library functions for user-defined floating point
128 types such as those from __multiprecision). In addition, the macro `BOOST_MATH_STD_USING`
129 should appear at the start of the function (note no semi-colon afterwards!) so that
130 all the math functions in `namespace std` are visible in the current scope.
131 * Calls to other special functions should be made as fully qualified calls, and include the
132 policy parameter as the last argument, for example `boost::math::tgamma(a, pol)`.
133 * Where possible, evaluation of series, continued fractions, polynomials, or root
134 finding should use one of the [link math_toolkit.internals_overview boiler-plate functions]. In any case, after
135 any iterative method, you should verify that the number of iterations did not exceed the
136 maximum specified in the __Policy type, and if it did terminate as a result of exceeding the
137 maximum, then the appropriate error handler should be called (see existing code for examples).
138 * Numeric constants such as [pi] etc should be obtained via a call to the [link math_toolkit.constants appropriate function],
139 for example: `constants::pi<T>()`.
140 * Where tables of coefficients are used (for example for rational approximations), care should be taken
141 to ensure these are initialized at program startup to ensure thread safety when using user-defined number types.
142 See for example the use of `erf_initializer` in [@../../include/boost/math/special_functions/erf.hpp erf.hpp].
144 Here are some other useful internal functions:
147 [[function][Meaning]]
148 [[`policies::digits<T, Policy>()`][Returns number of binary digits in T (possible overridden by the policy).]]
149 [[`policies::get_max_series_iterations<Policy>()`][Maximum number of iterations for series evaluation.]]
150 [[`policies::get_max_root_iterations<Policy>()`][Maximum number of iterations for root finding.]]
151 [[`polices::get_epsilon<T, Policy>()`][Epsilon for type T, possibly overridden by the Policy.]]
152 [[`tools::digits<T>()`][Returns the number of binary digits in T.]]
153 [[`tools::max_value<T>()`][Equivalent to `std::numeric_limits<T>::max()`]]
154 [[`tools::min_value<T>()`][Equivalent to `std::numeric_limits<T>::min()`]]
155 [[`tools::log_max_value<T>()`][Equivalent to the natural logarithm of `std::numeric_limits<T>::max()`]]
156 [[`tools::log_min_value<T>()`][Equivalent to the natural logarithm of `std::numeric_limits<T>::min()`]]
157 [[`tools::epsilon<T>()`][Equivalent to `std::numeric_limits<T>::epsilon()`.]]
158 [[`tools::root_epsilon<T>()`][Equivalent to the square root of `std::numeric_limits<T>::epsilon()`.]]
159 [[`tools::forth_root_epsilon<T>()`][Equivalent to the forth root of `std::numeric_limits<T>::epsilon()`.]]
164 [section:special_tut_test Testing]
166 We work under the assumption that untested code doesn't work, so some tests for your new special function are in order,
167 we'll divide these up in to 3 main categories:
171 Spot tests consist of checking that the expected exception is generated when the inputs are in error (or
172 otherwise generate undefined values), and checking any special values. We can check for expected exceptions
173 with `BOOST_CHECK_THROW`, so for example if it's a domain error for the last parameter to be outside the range
174 `[0,1]` then we might have:
176 BOOST_CHECK_THROW(my_special(0, -0.1), std::domain_error);
177 BOOST_CHECK_THROW(my_special(0, 1.1), std::domain_error);
179 When the function has known exact values (typically integer values) we can use `BOOST_CHECK_EQUAL`:
181 BOOST_CHECK_EQUAL(my_special(1.0, 0.0), 0);
182 BOOST_CHECK_EQUAL(my_special(1.0, 1.0), 1);
184 When the function has known values which are not exact (from a floating point perspective) then we can use
185 `BOOST_CHECK_CLOSE_FRACTION`:
187 // Assumes 4 epsilon is as close as we can get to a true value of 2Pi:
188 BOOST_CHECK_CLOSE_FRACTION(my_special(0.5, 0.5), 2 * constants::pi<double>(), std::numeric_limits<double>::epsilon() * 4);
190 [h4 Independent Test Values]
192 If the function is implemented by some other known good source (for example Mathematica or it's online versions
193 [@http://functions.wolfram.com functions.wolfram.com] or [@http://www.wolframalpha.com www.wolframalpha.com]
194 then it's a good idea to sanity check our implementation by having at least one independendly generated value
195 for each code branch our implementation may take. To slot these in nicely with our testing framework it's best to
196 tabulate these like this:
198 // function values calculated on http://functions.wolfram.com/
199 static const boost::array<boost::array<T, 3>, 10> my_special_data = {{
200 {{ SC_(0), SC_(0), SC_(1) }},
201 {{ SC_(0), SC_(1), SC_(1.26606587775200833559824462521471753760767031135496220680814) }},
202 /* More values here... */
205 We'll see how to use this table and the meaning of the `SC_` macro later. One important point
206 is to make sure that the input values have exact binary representations: so choose values such as
207 1.5, 1.25, 1.125 etc. This ensures that if `my_special` is unusually sensitive in one area, that
208 we don't get apparently large errors just because the inputs are 0.5 ulp in error.
210 [h4 Random Test Values]
212 We can generate a large number of test values to check both for future regressions, and for
213 accumulated rounding or cancellation error in our implementation. Ideally we would use an
214 independent implementation for this (for example my_special may be defined in directly terms
215 of other special functions but not implemented that way for performance or accuracy reasons).
216 Alternatively we may use our own implementation directly, but with any special cases (asymptotic
217 expansions etc) disabled. We have a set of [link math_toolkit.internals.test_data tools]
218 to generate test data directly, here's a typical example:
220 [import ../../example/special_data.cpp]
221 [special_data_example]
223 Typically several sets of data will be generated this way, including random values in some "normal"
224 range, extreme values (very large or very small), and values close to any "interesting" behaviour
225 of the function (singularities etc).
227 [h4 The Test File Header]
229 We split the actual test file into 2 distinct parts: a header that contains the testing code
230 as a series of function templates, and the actual .cpp test driver that decides which types
231 are tested, and sets the "expected" error rates for those types. It's done this way because:
233 * We want to test with both built in floating point types, and with multiprecision types.
234 However, both compile and runtimes with the latter can be too long for the folks who run
235 the tests to realistically cope with, so it makes sense to split the test into (at least)
237 * The definition of the SC_ macro used in our tables of data may differ depending on what type
238 we're testing (see below). Again this is largely a matter of managing compile times as large tables
239 of user-defined-types can take a crazy amount of time to compile with some compilers.
241 The test header contains 2 functions:
243 template <class Real, class T>
244 void do_test(const T& data, const char* type_name, const char* test_name);
247 void test(T, const char* type_name);
249 Before implementing those, we'll include the headers we'll need, and provide a default
250 definition for the SC_ macro:
252 // A couple of Boost.Test headers in case we need any BOOST_CHECK_* macros:
253 #include <boost/test/unit_test.hpp>
254 #include <boost/test/floating_point_comparison.hpp>
255 // Our function to test:
256 #include <boost/math/special_functions/my_special.hpp>
257 // We need boost::array for our test data, plus a few headers from
258 // libs/math/test that contain our testing machinary:
259 #include <boost/array.hpp>
260 #include "functor.hpp"
261 #include "handle_test_result.hpp"
262 #include "table_type.hpp"
265 #define SC_(x) static_cast<typename table_type<T>::type>(BOOST_JOIN(x, L))
268 The easiest function to implement is the "test" function which is what we'll be calling
269 from the test-driver program. It simply includes the files containing the tabular
270 test data and calls `do_test` function for each table, along with a description of what's
274 void test(T, const char* type_name)
277 // The actual test data is rather verbose, so it's in a separate file
279 // The contents are as follows, each row of data contains
280 // three items, input value a, input value b and my_special(a, b):
282 # include "my_special_1.ipp"
284 do_test<T>(my_special_1, name, "MySpecial Function: Mathematica Values");
286 # include "my_special_2.ipp"
288 do_test<T>(my_special_2, name, "MySpecial Function: Random Values");
290 # include "my_special_3.ipp"
292 do_test<T>(my_special_3, name, "MySpecial Function: Very Small Values");
295 The function `do_test` takes each table of data and calculates values for each row
296 of data, along with statistics for max and mean error etc, most of this is handled
297 by some boilerplate code:
299 template <class Real, class T>
300 void do_test(const T& data, const char* type_name, const char* test_name)
302 // Get the type of each row and each element in the rows:
303 typedef typename T::value_type row_type;
304 typedef Real value_type;
306 // Get a pointer to our function, we have to use a workaround here
307 // as some compilers require the template types to be explicitly
308 // specified, while others don't much like it if it is!
309 typedef value_type (*pg)(value_type, value_type);
310 #if defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
311 pg funcp = boost::math::my_special<value_type, value_type>;
313 pg funcp = boost::math::my_special;
316 // Somewhere to hold our results:
317 boost::math::tools::test_result<value_type> result;
318 // And some pretty printing:
319 std::cout << "Testing " << test_name << " with type " << type_name
320 << "\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
323 // Test my_special against data:
325 result = boost::math::tools::test_hetero<Real>(
326 /* First argument is the table */
328 /* Next comes our function pointer, plus the indexes of it's arguments in the table */
329 bind_func<Real>(funcp, 0, 1),
330 /* Then the index of the result in the table - potentially we can test several
331 related functions this way, each having the same input arguments, and different
332 output values in different indexes in the table */
333 extract_result<Real>(2));
335 // Finish off with some boilerplate to check the results were within the expected errors,
336 // and pretty print the results:
338 handle_test_result(result, data[result.worst()], result.worst(), type_name, "boost::math::my_special", test_name);
341 Now we just need to write the test driver program, at it's most basic it looks something like this:
343 #include <boost/math/special_functions/math_fwd.hpp>
344 #include <boost/math/tools/test.hpp>
345 #include <boost/math/tools/stats.hpp>
346 #include <boost/type_traits.hpp>
347 #include <boost/array.hpp>
348 #include "functor.hpp"
350 #include "handle_test_result.hpp"
351 #include "test_my_special.hpp"
353 BOOST_AUTO_TEST_CASE( test_main )
356 // Test each floating point type, plus real_concept.
357 // We specify the name of each type by hand as typeid(T).name()
358 // often gives an unreadable mangled name.
363 // Testing of long double and real_concept is protected
364 // by some logic to disable these for unsupported
365 // or problem compilers.
367 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
368 test(0.1L, "long double");
369 #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
370 #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
371 test(boost::math::concepts::real_concept(0.1), "real_concept");
375 std::cout << "<note>The long double tests have been disabled on this platform "
376 "either because the long double overloads of the usual math functions are "
377 "not available at all, or because they are too inaccurate for these tests "
378 "to pass.</note>" << std::cout;
382 That's almost all there is too it - except that if the above program is run it's very likely that
383 all the tests will fail as the default maximum allowable error is 1 epsilon. So we'll
384 define a function (don't forget to call it from the start of the `test_main` above) to
385 up the limits to something sensible, based both on the function we're calling and on
386 the particular tests plus the platform and compiler:
388 void expected_results()
391 // Define the max and mean errors expected for
392 // various compilers and platforms.
394 const char* largest_type;
395 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
396 if(boost::math::policies::digits<double, boost::math::policies::policy<> >() == boost::math::policies::digits<long double, boost::math::policies::policy<> >())
398 largest_type = "(long\\s+)?double|real_concept";
402 largest_type = "long double|real_concept";
405 largest_type = "(long\\s+)?double";
408 // We call add_expected_result for each error rate we wish to adjust, these tell
409 // handle_test_result what level of error is acceptable. We can have as many calls
410 // to add_expected_result as we need, each one establishes a rule for acceptable error
411 // with rules set first given preference.
414 /* First argument is a regular expression to match against the name of the compiler
415 set in BOOST_COMPILER */
417 /* Second argument is a regular expression to match against the name of the
418 C++ standard library as set in BOOST_STDLIB */
420 /* Third argument is a regular expression to match against the name of the
421 platform as set in BOOST_PLATFORM */
423 /* Forth argument is the name of the type being tested, normally we will
424 only need to up the acceptable error rate for the widest floating
425 point type being tested */
427 /* Fifth argument is a regular expression to match against
428 the name of the group of data being tested */
429 "MySpecial Function:.*Small.*",
430 /* Sixth argument is a regular expression to match against the name
431 of the function being tested */
432 "boost::math::my_special",
433 /* Seventh argument is the maximum allowable error expressed in units
434 of machine epsilon passed as a long integer value */
436 /* Eighth argument is the maximum allowable mean error expressed in units
437 of machine epsilon passed as a long integer value */
441 [h4 Testing Multiprecision Types]
443 Testing of multiprecision types is handled by the test drivers in libs/multiprecision/test/math,
444 please refer to these for examples. Note that these tests are run only occationally as they take
445 a lot of CPU cycles to build and run.
447 [h4 Improving Compile Times]
449 As noted above, these test programs can take a while to build as we're instantiating a lot of templates
450 for several different types, and our test runners are already stretched to the limit, and probably
451 using outdated "spare" hardware. There are two things we can do to speed things up:
453 * Use a precompiled header.
454 * Use separate compilation of our special function templates.
456 We can make these changes by changing the list of includes from:
458 #include <boost/math/special_functions/math_fwd.hpp>
459 #include <boost/math/tools/test.hpp>
460 #include <boost/math/tools/stats.hpp>
461 #include <boost/type_traits.hpp>
462 #include <boost/array.hpp>
463 #include "functor.hpp"
465 #include "handle_test_result.hpp"
469 #include <pch_light.hpp>
473 #include <boost/math/special_functions/my_special.hpp>
477 #include <boost/math/special_functions/math_fwd.hpp>
479 The Jamfile target that builds the test program will need the targets
481 test_instances//test_instances pch_light
483 adding to it's list of source dependencies (see the Jamfile for examples).
485 Finally the project in libs/math/test/test_instances will need modifying
486 to instantiate function `my_special`.
488 These changes should be made last, when `my_special` is stable and the code is in Trunk.
492 Our concept checks verify that your function's implementation makes no assumptions that aren't
493 required by our [link math_toolkit.real_concepts Real number conceptual requirements]. They also
494 check for various common bugs and programming traps that we've fallen into over time. To
495 add your function to these tests, edit libs/math/test/compile_test/instantiate.hpp to add
496 calls to your function: there are 7 calls to each function, each with a different purpose.
497 Search for something like "ibeta" or "gamm_p" and follow their examples.
504 Copyright 2013 John Maddock.
505 Distributed under the Boost Software License, Version 1.0.
506 (See accompanying file LICENSE_1_0.txt or copy at
507 http://www.boost.org/LICENSE_1_0.txt).