diff --git a/doc/equations/bernoulli_numbers.png b/doc/equations/bernoulli_numbers.png
index bb4353612..f8de0ef35 100644
Binary files a/doc/equations/bernoulli_numbers.png and b/doc/equations/bernoulli_numbers.png differ
diff --git a/doc/equations/bernoulli_numbers.svg b/doc/equations/bernoulli_numbers.svg
index e1963536f..dec5eb1d3 100644
--- a/doc/equations/bernoulli_numbers.svg
+++ b/doc/equations/bernoulli_numbers.svg
@@ -1,341 +1,2 @@
-
-
+
+Bj=(−1)j2−1jTj/2(4j−2j)
\ No newline at end of file
diff --git a/doc/equations/bernoulli_numbers2.mml b/doc/equations/bernoulli_numbers2.mml
new file mode 100644
index 000000000..53816d7e5
--- /dev/null
+++ b/doc/equations/bernoulli_numbers2.mml
@@ -0,0 +1,230 @@
+
+
+
\ No newline at end of file
diff --git a/doc/equations/bernoulli_numbers2.png b/doc/equations/bernoulli_numbers2.png
new file mode 100644
index 000000000..c88976447
Binary files /dev/null and b/doc/equations/bernoulli_numbers2.png differ
diff --git a/doc/equations/bernoulli_numbers2.svg b/doc/equations/bernoulli_numbers2.svg
new file mode 100644
index 000000000..06ae9c70a
--- /dev/null
+++ b/doc/equations/bernoulli_numbers2.svg
@@ -0,0 +1,2 @@
+
+ln(Bn)=(12+n)ln(n)+(12−n)ln(π)+(32−n)ln(2)−R(n)R(n)=n(1−112(1−130(1−27n−2)n−2)n−2)
\ No newline at end of file
diff --git a/doc/equations/generate.sh b/doc/equations/generate.sh
index 27d1d5eff..77e0427da 100755
--- a/doc/equations/generate.sh
+++ b/doc/equations/generate.sh
@@ -7,9 +7,9 @@
#
# Paths to tools come first, change these to match your system:
#
-math2svg='m:\download\open\SVGMath-0.3.1\math2svg.py'
-python=/cygdrive/c/Python26/python.exe
-inkscape=/cygdrive/c/progra~1/Inkscape/inkscape
+math2svg='d:\download\open\SVGMath-0.3.1\math2svg.py'
+python=/cygdrive/c/Python27/python.exe
+inkscape=/cygdrive/c/progra~1/Inkscape/inkscape.exe
# Image DPI:
dpi=120
@@ -18,7 +18,7 @@ for mmlfile in $*; do
pngfile=$(basename $svgfile .svg).png
tempfile=temp.mml
# strip html wrappers put in by MathCast:
- cat $mmlfile | tr -d "\r\n" | sed -e 's/.*\(]*>.*<\/math>\).*/\1/' > $tempfile
+ cat $mmlfile | tr -d "\r\n" | sed -e 's/.*\(]*>.*<\/math>\).*/\1/' -e 's///g' -e 's/<\/semantics>//g' > $tempfile
echo Generating $svgfile
$python $math2svg $tempfile > $svgfile
@@ -31,3 +31,4 @@ done
+
diff --git a/doc/sf/bernoulli_numbers.qbk b/doc/sf/bernoulli_numbers.qbk
index 6ef84ae5f..8551216e3 100644
--- a/doc/sf/bernoulli_numbers.qbk
+++ b/doc/sf/bernoulli_numbers.qbk
@@ -18,16 +18,23 @@ including the __tgamma, __lgamma and polygamma functions.
namespace boost { namespace math {
template
- T bernoulli_b2n(const int i); // Single Bernoulli number (default policy).
+ T bernoulli_b2n(const int n); // Single Bernoulli number (default policy).
template
- T bernoulli_b2n(const int i, const Policy &pol); // User policy for errors etc.
+ T bernoulli_b2n(const int n, const Policy &pol); // User policy for errors etc.
}} // namespaces
[h4 Description]
-Both return the (2 * i)[super th] Bernoulli number.
+Both return the (2 * n)[super th] Bernoulli number B[sub 2n].
+
+Note that since all odd numbered Bernoulli numbers are zero (apart from B[sub 1] which is [plusminus][frac12])
+the interface will only return the even numbered Bernoulli numbers.
+
+This function uses fast table lookup for low-indexed Bernoulli numbers, while larger values are calculated
+as needed and then cached. The caching mechanism requires a certain amount of thread safety code, so
+`unchecked_bernoulli_b2n` may provide a better interface for performance critical code.
The final __Policy argument is optional and can be used to control the behaviour of the function:
how it handles errors, what level of precision to use, etc.
@@ -45,23 +52,32 @@ Refer to __policy_section for more details.
[h4 Synopsis]
``
-#include
+#include
``
+ template <>
+ struct max_bernoulli_b2n;
+
template
inline T unchecked_bernoulli_b2n(unsigned n);
`unchecked_bernoulli_b2n` provides access to Bernoulli numbers [*without any checks for overflow or invalid parameters].
-(Not recommended for normal use, but bypassing overflow checks is quicker.)
+It is implemented as a direct (and very fast) table lookup, and while not recomended for general use it can be useful
+inside inner loops where the ultimate performance is required, and error checking is moved outside the loop.
- template <>
- struct max_bernoulli_index::value
+The largest value you can pass to `unchecked_bernoulli_b2n<>` is `max_bernoulli_b2n<>::value`: passing values greater than
+that will result in a buffer overrun error, so it's clearly important to place the error handling in your own code
+when using this direct interface.
-The largest value you can (without overflow) pass to `unchecked_bernoulli_b2n<>`
-is `max_bernoulli_b2n<>::value`.
+The value of `boost::math::max_bernoulli_b2n::value` varies by the type T, for types `float`/`double`/`long double`
+it's the largest value which doesn't overflow the target type: for example, `boost::math::max_bernoulli_b2n::value` is 129.
+However, for multiprecision types, it's the largest value for which the result can be represented as the ratio of two 64-bit
+integers, for example `boost::math::max_bernoulli_b2n::value` is just 17. Of course
+larger indexes can be passed to `bernoulli_b2n(n)`, but then then you loose fast table lookup (i.e. values may need to be calculated).
-For example, `boost::math::max_bernoulli_b2n::value` is 129.
+[bernoulli_example_4]
+[bernoulli_output_4]
[h4 Multiple Bernoulli Numbers]
@@ -96,20 +112,15 @@ with one call (one with default policy and the other allowing a user-defined pol
These return a series of Bernoulli numbers:
- Bernoulli(2*start_index),
- Bernoulli(2*(start_index+1))
- ...
- Bernoulli(2*(number_of_bernoullis_b2n-1))
+[:B[sub 2*start_index],B[sub 2*(start_index+1)],...,B[sub 2*(start_index+number_of_bernoullis_b2n-1)]]
[h4 Examples]
[bernoulli_example_2]
[bernoulli_output_2]
[bernoulli_example_3]
[bernoulli_output_3]
-[bernoulli_example_4]
-[bernoulli_output_4]
-The source of this example is at [../../example/bernoulli_example.cpp bernoulli_example.cpp]
+The source of this example is at [@../../example/bernoulli_example.cpp bernoulli_example.cpp]
[h4 Accuracy]
@@ -117,9 +128,10 @@ All the functions usually return values within one ULP (unit in the last place)
[h4 Implementation]
-The implementation details are in [@../../include/boost/math/special_functions/detail/bernoulli_details.hpp bernoulli_details.hpp].
+The implementation details are in [@../../include/boost/math/special_functions/detail/bernoulli_details.hpp bernoulli_details.hpp]
+and [@../../include/boost/math/special_functions/detail/unchecked_bernoulli.hpp unchecked_bernoulli.hpp].
-For `i <= max_bernoulli_index::value` this is implemented by table lookup;
+For `i <= max_bernoulli_index::value` this is implemented by simple table lookup from a statically initialized table;
for larger values of `i`, this is implemented by the Tangent Numbers algorithm as described in the paper:
Fast Computation of Bernoulli, Tangent and Secant Numbers, Richard P. Brent and David Harvey,
[@http://arxiv.org/pdf/1108.0286v3.pdf] (2011).
@@ -128,8 +140,6 @@ Fast Computation of Bernoulli, Tangent and Secant Numbers, Richard P. Brent and
(an even alternating permutation number) are defined
and their generating function is also given therein.
- [/@http://mathworld.wolfram.com/images/equations/TangentNumber/Inline15.gif]
-
The relation of Tangent numbers with Bernoulli numbers ['B[sub i]]
is given by Brent and Harvey's equation 14:
@@ -143,7 +153,13 @@ elseif i == 0 then ['B[sub i]] = 1 [br]
elseif i == 1 then ['B[sub i]] = -1/2 [br]
elseif i < 0 or i is odd then ['B[sub i]] = 0
-[/@http://s7.postimg.org/mygi2kror/bernoulli_1.png]
+Note that computed values are stored in a fixed-size table, access is thread safe via atomic operations (i.e. lock
+free programming), this imparts a much lower overhead on access to cached values than might overwise be expected -
+typically for multiprecision types the cost of thread synchronisation is negligable, while for built in types
+this code is not normally executed anyway. For very large arguments which cannot be reasonably computed or
+stored in our cache, an asymptotic expansion [@http://www.luschny.de/math/primes/bernincl.html due to Luschny] is used:
+
+[equation bernoulli_numbers2]
[endsect] [/section:bernoulli_numbers Bernoulli Numbers]
diff --git a/example/bernoulli_example.cpp b/example/bernoulli_example.cpp
index ed373c6fe..5a2c58969 100644
--- a/example/bernoulli_example.cpp
+++ b/example/bernoulli_example.cpp
@@ -34,7 +34,8 @@ int main()
{
//[bernoulli_example_1
-/*`A simple example computes the value of `Bernoulli(2)` where the return type is `double`.
+/*`A simple example computes the value of B[sub 4] where the return type is `double`,
+note that the argument to bernoulli_b2n is ['2] not ['4] since it computes B[sub 2N].
*/
@@ -48,7 +49,7 @@ int main()
/*`So B[sub 4] == -1/30 == -0.0333333333333333
If we use Boost.Multiprecision and its 50 decimal digit floating-point type `cpp_dec_float_50`,
-we can calculate the value of much larger numbers like `Bernoulli(100)`
+we can calculate the value of much larger numbers like B[sub 200]
and also obtain much higher precision.
*/
@@ -107,11 +108,7 @@ and we will get a helpful error message (provided try'n'catch blocks are used).
//] //[/bernoulli_example_3]
//[bernoulli_example_4
-/*`Users can determine from `max_bernoulli_b2n<>::value`
-the largest Bernoulli number you can evaluate for any type
-(or safely pass to `unchecked_bernoulli_b2n<>`).
-
-For example:
+/*For example:
*/
std::cout << "boost::math::max_bernoulli_b2n::value = " << boost::math::max_bernoulli_b2n::value << std::endl;
std::cout << "Maximum Bernoulli number using float is " << boost::math::bernoulli_b2n( boost::math::max_bernoulli_b2n::value) << std::endl;