Imported from /home/lorry/working-area/delta_boost-tarball/boost_1_58_0.tar.bz2.HEADboost_1_58_0master

9 files changed, 726 insertions, 9 deletions

diff --git a/libs/math/doc/sf/bernoulli_numbers.qbk b/libs/math/doc/sf/bernoulli_numbers.qbk
new file mode 100644
index 000000000..d97cad6e2
--- /dev/null
+++ b/libs/math/doc/sf/bernoulli_numbers.qbk
@@ -0,0 +1,223 @@
+[section:bernoulli_numbers Bernoulli Numbers]
+
+[@https://en.wikipedia.org/wiki/Bernoulli_number Bernoulli numbers]
+are a sequence of rational numbers useful for the Taylor series expansion,
+Euler-Maclaurin formula, and the Riemann zeta function.
+
+Bernoulli numbers are used in evaluation of some Boost.Math functions,
+including the __tgamma, __lgamma and polygamma functions.
+
+[h4 Single Bernoulli number]
+
+[h4 Synopsis]
+
+``
+#include <boost/math/special_functions/bernoulli.hpp>
+``
+
+  namespace boost { namespace math {
+
+  template <class T>
+  T bernoulli_b2n(const int n);  // Single Bernoulli number (default policy).
+
+  template <class T, class Policy>
+  T bernoulli_b2n(const int n, const Policy &pol); // User policy for errors etc.
+
+  }} // namespaces
+
+[h4 Description]
+
+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.
+
+Refer to __policy_section for more details.
+
+[h4 Examples]
+
+[import ../../example/bernoulli_example.cpp]
+[bernoulli_example_1]
+
+[bernoulli_output_1]
+
+[h4 Single (unchecked) Bernoulli number]
+
+[h4 Synopsis]
+``
+#include <boost/math/special_functions/bernoulli.hpp>
+
+``
+
+  template <>
+  struct max_bernoulli_b2n<T>;
+
+  template<class T>
+  inline T unchecked_bernoulli_b2n(unsigned n);
+
+`unchecked_bernoulli_b2n` provides access to Bernoulli numbers [*without any checks for overflow or invalid parameters].
+It is implemented as a direct (and very fast) table lookup, and while not recommended for general use it can be useful
+inside inner loops where the ultimate performance is required, and error checking is moved outside the loop.
+
+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 value of `boost::math::max_bernoulli_b2n<T>::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<double>::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<boost::multiprecision::cpp_dec_float_50>::value` is just 17.  Of course
+larger indexes can be passed to `bernoulli_b2n<T>(n)`, but then you lose fast table lookup (i.e. values may need to be calculated).
+
+[bernoulli_example_4]
+[bernoulli_output_4]
+
+[h4 Multiple Bernoulli Numbers]
+
+[h4 Synopsis]
+
+``
+#include <boost/math/special_functions/bernoulli.hpp>
+``
+
+  namespace boost { namespace math {
+
+  // Multiple Bernoulli numbers (default policy).
+  template <class T, class OutputIterator>
+  OutputIterator bernoulli_b2n(
+    int start_index,
+    unsigned number_of_bernoullis_b2n,
+    OutputIterator out_it);
+
+  // Multiple Bernoulli numbers (user policy).
+  template <class T, class OutputIterator, class Policy>
+  OutputIterator bernoulli_b2n(
+    int start_index,
+    unsigned number_of_bernoullis_b2n,
+    OutputIterator out_it,
+    const Policy& pol);
+  }} // namespaces
+
+[h4 Description]
+
+Two versions of the Bernoulli number function are provided to compute multiple Bernoulli numbers
+with one call (one with default policy and the other allowing a user-defined policy).
+
+These return a series of Bernoulli numbers:
+
+[: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]
+
+The source of this example is at [@../../example/bernoulli_example.cpp bernoulli_example.cpp]
+
+[h4 Accuracy]
+
+All the functions usually return values within one ULP (unit in the last place) for the floating-point type.
+
+[h4 Implementation]
+
+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<T>::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).
+
+[@http://mathworld.wolfram.com/TangentNumber.html Tangent (or Zag) numbers]
+(an even alternating permutation number) are defined
+and their generating function is also given therein.
+
+The relation of Tangent numbers with Bernoulli numbers  ['B[sub i]]
+is given by Brent and Harvey's equation 14:
+
+__spaces[equation tangent_numbers]
+
+Their relation with Bernoulli numbers ['B[sub i]] are defined by
+
+if i > 0 and i is even then
+__spaces[equation bernoulli_numbers] [br]
+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
+
+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 otherwise be expected - 
+typically for multiprecision types the cost of thread synchronisation is negligible, 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]
+
+
+[section:tangent_numbers Tangent Numbers]
+
+[@http://en.wikipedia.org/wiki/Tangent_numbers Tangent numbers],
+also called a zag function.  See also
+[@http://mathworld.wolfram.com/TangentNumber.html Tangent number].
+
+From the number, An, of alternating permutations of the set {1, ..., n},
+the numbers A2n+1 with odd indices are called tangent numbers or zag numbers.
+The first few values are 1, 2, 16, 272, 7936, 353792, 22368256, 1903757312 ...
+(sequence [@http://oeis.org/A000182 A000182 in OEIS]).
+They are called tangent numbers because they appear as
+numerators in the Maclaurin series of tan x.
+
+Tangent numbers are used in the computation of Bernoulli numbers,
+but are also made available here.
+
+[h4 Synopsis]
+``
+#include <boost/math/special_functions/detail/bernoulli.hpp>
+``
+
+  template <class T>
+  T tangent_t2n(const int i);  // Single tangent number (default policy).
+
+  template <class T, class Policy>
+  T tangent_t2n(const int i, const Policy &pol); // Single tangent number (user policy).
+
+  // Multiple tangent numbers (default policy).
+  template <class T, class OutputIterator>
+  OutputIterator tangent_t2n(const int start_index,
+                                      const unsigned number_of_tangent_t2n,
+                                      OutputIterator out_it);
+
+  // Multiple tangent numbers (user policy).
+  template <class T, class OutputIterator, class Policy>
+  OutputIterator tangent_t2n(const int start_index,
+                                      const unsigned number_of_tangent_t2n,
+                                      OutputIterator out_it,
+                                      const Policy& pol);
+
+[h4 Examples]
+
+[tangent_example_1]
+
+The output is:
+[tangent_output_1]
+
+The source of this example is at [../../example/bernoulli_example.cpp bernoulli_example.cpp]
+
+[endsect] [/section:tangent_numbers Tangent Numbers]
+
+[/
+  Copyright 2013, 2014 Nikhar Agrawal, Christopher Kormanyos, John Maddock, Paul A. Bristow.
+  Distributed under 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).
+]
diff --git a/libs/math/doc/sf/digamma.qbk b/libs/math/doc/sf/digamma.qbk
index cf3b9c232..d62466ac3 100644
--- a/libs/math/doc/sf/digamma.qbk
+++ b/libs/math/doc/sf/digamma.qbk
@@ -27,10 +27,6 @@ derivative of the gamma function:
 
 [optional_policy]
 
-There is no fully generic version of this function:  all the implementations
-are tuned to specific accuracy levels, the most precise of which delivers
-34-digits of precision.
-
 The return type of this function is computed using the __arg_pomotion_rules:
 the result is of type `double` when T is an integer type, and type T otherwise.
 
@@ -97,6 +93,14 @@ Choosing BIG=10 for up to 80-bit reals, and BIG=20 for 128-bit reals allows
 the series to truncated after a suitably small number of terms and evaluated
 as a polynomial in `1/(x*x)`.
 
+The arbitrary precision version of this function uses recurrence relations until
+x > BIG, and then evaluation via the asymptotic expansion above.  As special cases
+integer and half integer arguments are handled via:
+
+[equation digamma4]
+
+[equation digamma5]
+
 The rational approximation [jm_rationals] in the range [1,2] is derived as follows.
 
 First a high precision approximation to digamma was constructed using a 60-term
diff --git a/libs/math/doc/sf/ellint_carlson.qbk b/libs/math/doc/sf/ellint_carlson.qbk
index 512c4d416..b6f2aa42f 100644
--- a/libs/math/doc/sf/ellint_carlson.qbk
+++ b/libs/math/doc/sf/ellint_carlson.qbk
@@ -69,6 +69,21 @@ LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
 
   }} // namespaces
 
+``
+  #include <boost/math/special_functions/ellint_rg.hpp>
+``
+
+  namespace boost { namespace math {
+
+  template <class T1, class T2, class T3>
+  ``__sf_result`` ellint_rg(T1 x, T2 y, T3 z)
+
+  template <class T1, class T2, class T3, class ``__Policy``>
+  ``__sf_result`` ellint_rg(T1 x, T2 y, T3 z, const ``__Policy``&)
+
+  }} // namespaces
+
+
 
 [heading Description]
 
@@ -154,6 +169,20 @@ using the relation:
 
 [equation ellint18]
 
+  template <class T1, class T2, class T3>
+  ``__sf_result`` ellint_rg(T1 x, T2 y, T3 z)
+  
+  template <class T1, class T2, class T3, class ``__Policy``>
+  ``__sf_result`` ellint_rg(T1 x, T2 y, T3 z, const ``__Policy``&)
+  
+Returns Carlson's elliptic integral R[sub G]:
+
+[equation ellint27]
+
+Requires that x and y are non-negative, otherwise returns the result of __domain_error.
+
+[optional_policy]
+
 [heading Testing]
 
 There are two sets of tests.
@@ -199,7 +228,7 @@ closer and closer. When they are nearly equal, the special case equation
 
 is used. More specifically, ['[R F]] is evaluated from a Taylor series
 expansion to the fifth order. The calculations of the other three integrals
-are analogous.
+are analogous, except for R[sub C] which can be computed from elementary functions.
 
 For ['p < 0] in ['R[sub J](x, y, z, p)] and ['y < 0] in ['R[sub C](x, y)],
 the integrals are singular and their
diff --git a/libs/math/doc/sf/ellint_introduction.qbk b/libs/math/doc/sf/ellint_introduction.qbk
index 12a375a5f..80c11f659 100644
--- a/libs/math/doc/sf/ellint_introduction.qbk
+++ b/libs/math/doc/sf/ellint_introduction.qbk
@@ -103,6 +103,14 @@ Complete Elliptic Integral of the Third Kind (Legendre form)
 
 [equation ellint8]
 
+Legendre also defined a forth integral D([phi],k) which is a combination of the other three:
+
+[equation ellint_d]
+
+Like the other Legendre integrals this comes in both complete and incomplete forms.
+
+[h4 Carlson Elliptic Integrals]
+
 Carlson [[link ellint_ref_carlson77 Carlson77]] [[link ellint_ref_carlson78  Carlson78]] gives an alternative definition of
 elliptic integral's canonical forms:
 
@@ -137,6 +145,9 @@ where ['x] is nonnegative and ['y] is nonzero.
 
 ['R[sub D](x, y, z) = R[sub J](x, y, z, z)]]
 
+Carlson's Symmetric Integral
+
+[equation ellint27]
 
 [h4 Duplication Theorem]
 
@@ -155,6 +166,20 @@ In particular,
 
 [equation ellint15]
 
+[h4 Miscellaneous Elliptic Integrals]
+
+There are two functions related to the elliptic integrals which otherwise
+defy categorisation, these are the Jacobi Zeta function:
+
+[equation jacobi_zeta]
+
+and the Heuman Lambda function:
+
+[equation heuman_lambda]
+
+Both of these functions are easily implemented in terms of Carlson's integrals, and are
+provided in this library as __jacobi_zeta and __heuman_lambda.
+
 [h4 Numerical Algorithms]
 
 The conventional methods for computing elliptic integrals are Gauss
diff --git a/libs/math/doc/sf/ellint_legendre.qbk b/libs/math/doc/sf/ellint_legendre.qbk
index 0808058df..a366db948 100644
--- a/libs/math/doc/sf/ellint_legendre.qbk
+++ b/libs/math/doc/sf/ellint_legendre.qbk
@@ -337,3 +337,228 @@ and
 [equation ellint26]
 
 [endsect]
+
+[section:ellint_d Elliptic Integral D - Legendre Form]
+
+[heading Synopsis]
+
+``
+  #include <boost/math/special_functions/ellint_d.hpp>
+``
+
+  namespace boost { namespace math {
+
+  template <class T1, class T2>
+  ``__sf_result`` ellint_d(T1 k, T2 phi);
+
+  template <class T1, class T2, class ``__Policy``>
+  ``__sf_result`` ellint_d(T1 k, T2 phi, const ``__Policy``&);
+
+  template <class T1>
+  ``__sf_result`` ellint_d(T1 k);
+
+  template <class T1, class ``__Policy``>
+  ``__sf_result`` ellint_d(T1 k, const ``__Policy``&);
+
+  }} // namespaces
+  
+[heading Description]
+
+These two functions evaluate the incomplete elliptic integral 
+['D([phi], k)] and its complete counterpart ['D(k) = D([pi]/2, k)].
+
+The return type of these functions is computed using the __arg_pomotion_rules
+when the arguments are of different types: when they are the same type then the result
+is the same type as the arguments.
+
+  template <class T1, class T2>
+  ``__sf_result`` ellint_d(T1 k, T2 phi);
+  
+  template <class T1, class T2, class ``__Policy``>
+  ``__sf_result`` ellint_3(T1 k, T2 phi, const ``__Policy``&);
+  
+Returns the incomplete elliptic integral:
+
+[equation ellint_d]
+
+Requires ['-1 <= k <= 1], otherwise 
+returns the result of __domain_error (outside this range the result 
+would be complex).
+
+[optional_policy]
+
+  template <class T1>
+  ``__sf_result`` ellint_d(T1 k);
+  
+  template <class T1, class ``__Policy``>
+  ``__sf_result`` ellint_d(T1 k, const ``__Policy``&);
+  
+Returns the complete elliptic integral ['D(k) = D([pi]/2, k)]
+
+Requires ['-1 <= k <= 1] otherwise returns the 
+result of __domain_error (outside this range the result would be complex).
+
+[optional_policy]
+
+[heading Accuracy]
+
+These functions are trivially computed in terms of other elliptic integrals
+and generally have very low error rates (a few epsilon) unless parameter [phi]
+is very large, in which case the usual trigonometric function argument-reduction issues apply.
+
+[heading Testing]
+
+The tests use a mixture of spot test values calculated using 
+values calculated at wolframalpha.com, and random test data generated using
+MPFR at 1000-bit precision and a deliberately naive implementation in terms of
+the Legendre integrals.
+
+[heading Implementation]
+
+The implementation for D([phi], k) first performs argument reduction using the relations:
+
+['D(-[phi], k) = -D([phi], k)]
+
+and
+
+['D(n[pi]+[phi], k) = 2nD(k) + D([phi], k)]
+
+to move [phi][space] to the range \[0, [pi]\/2\].
+
+The functions are then implemented in terms of Carlson's integral R[sub D]
+using the relation:
+
+[equation ellint_d]
+
+[endsect]
+
+[section:jacobi_zeta Jacobi Zeta Function]
+
+[heading Synopsis]
+
+``
+  #include <boost/math/special_functions/jacobi_zeta.hpp>
+``
+
+  namespace boost { namespace math {
+
+  template <class T1, class T2>
+  ``__sf_result`` jacobi_zeta(T1 k, T2 phi);
+
+  template <class T1, class T2, class ``__Policy``>
+  ``__sf_result`` jacobi_zeta(T1 k, T2 phi, const ``__Policy``&);
+
+  }} // namespaces
+  
+[heading Description]
+
+This function evaluates the Jacobi Zeta Function ['Z([phi], k)]
+
+[equation jacobi_zeta]
+
+The return type of this function is computed using the __arg_pomotion_rules
+when the arguments are of different types: when they are the same type then the result
+is the same type as the arguments.
+
+Requires ['-1 <= k <= 1], otherwise 
+returns the result of __domain_error (outside this range the result 
+would be complex).
+
+[optional_policy]
+
+Note that there is no complete analogue of this function (where [phi] = [pi] / 2)
+as this takes the value 0 for all ['k].
+
+[heading Accuracy]
+
+These functions are trivially computed in terms of other elliptic integrals
+and generally have very low error rates (a few epsilon) unless parameter [phi]
+is very large, in which case the usual trigonometric function argument-reduction issues apply.
+
+[heading Testing]
+
+The tests use a mixture of spot test values calculated using 
+values calculated at wolframalpha.com, and random test data generated using
+MPFR at 1000-bit precision and a deliberately naive implementation in terms of
+the Legendre integrals.
+
+[heading Implementation]
+
+The implementation for Z([phi], k) first makes the argument [phi] positive using:
+
+['Z(-[phi], k) = -Z([phi], k)]
+
+The function is then implemented in terms of Carlson's integral R[sub J]
+using the relation:
+
+[equation jacobi_zeta]
+
+There is one special case where the above relation fails: when ['k = 1], in that case
+the function simplifies to
+
+['Z([phi], 1) = sign(cos([phi])) sin([phi])]
+
+[endsect]
+
+[section:heuman_lambda Heuman Lambda Function]
+
+[heading Synopsis]
+
+``
+  #include <boost/math/special_functions/heuman_lambda.hpp>
+``
+
+  namespace boost { namespace math {
+
+  template <class T1, class T2>
+  ``__sf_result`` heuman_lambda(T1 k, T2 phi);
+
+  template <class T1, class T2, class ``__Policy``>
+  ``__sf_result`` heuman_lambda(T1 k, T2 phi, const ``__Policy``&);
+
+  }} // namespaces
+  
+[heading Description]
+
+This function evaluates the Heuman Lambda Function ['[Lambda][sub 0]([phi], k)]
+
+[equation heuman_lambda]
+
+The return type of this function is computed using the __arg_pomotion_rules
+when the arguments are of different types: when they are the same type then the result
+is the same type as the arguments.
+
+Requires ['-1 <= k <= 1], otherwise 
+returns the result of __domain_error (outside this range the result 
+would be complex).
+
+[optional_policy]
+
+Note that there is no complete analogue of this function (where [phi] = [pi] / 2)
+as this takes the value 1 for all ['k].
+
+[heading Accuracy]
+
+These functions are trivially computed in terms of other elliptic integrals
+and generally have very low error rates (a few epsilon) unless parameter [phi]
+is very large, in which case the usual trigonometric function argument-reduction issues apply.
+
+[heading Testing]
+
+The tests use a mixture of spot test values calculated using 
+values calculated at wolframalpha.com, and random test data generated using
+MPFR at 1000-bit precision and a deliberately naive implementation in terms of
+the Legendre integrals.
+
+[heading Implementation]
+
+The function is then implemented in terms of Carlson's integrals R[sub J] and R[sub F]
+using the relation:
+
+[equation heuman_lambda]
+
+This relation fails for ['|[phi]| >= [pi]/2] in which case the definition in terms of the
+Jacobi Zeta is used.
+
+[endsect]
+
diff --git a/libs/math/doc/sf/number_series.qbk b/libs/math/doc/sf/number_series.qbk
index 1e1019247..631cc294d 100644
--- a/libs/math/doc/sf/number_series.qbk
+++ b/libs/math/doc/sf/number_series.qbk
@@ -65,7 +65,7 @@ Refer to __policy_section for more details.
   inline T unchecked_bernoulli_b2n(unsigned n);
 
 `unchecked_bernoulli_b2n` provides access to Bernoulli numbers [*without any checks for overflow or invalid parameters].
-It is implemented as a direct (and very fast) table lookup, and while not recomended for general use it can be useful
+It is implemented as a direct (and very fast) table lookup, and while not recommended for general use it can be useful
 inside inner loops where the ultimate performance is required, and error checking is moved outside the loop.
 
 The largest value you can pass to `unchecked_bernoulli_b2n<>` is `max_bernoulli_b2n<>::value`: passing values greater than
@@ -76,7 +76,7 @@ The value of `boost::math::max_bernoulli_b2n<T>::value` varies by the type T, fo
 it's the largest value which doesn't overflow the target type: for example, `boost::math::max_bernoulli_b2n<double>::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<boost::multiprecision::cpp_dec_float_50>::value` is just 17.  Of course
-larger indexes can be passed to `bernoulli_b2n<T>(n)`, but then then you loose fast table lookup (i.e. values may need to be calculated).
+larger indexes can be passed to `bernoulli_b2n<T>(n)`, but then you lose fast table lookup (i.e. values may need to be calculated).
 
 [bernoulli_example_4]
 [bernoulli_output_4]
@@ -156,8 +156,8 @@ elseif i == 1 then  ['B[sub i]] = -1/2 [br]
 elseif i < 0 or i is odd then ['B[sub i]] = 0
 
 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
+free programming), this imparts a much lower overhead on access to cached values than might otherwise be expected - 
+typically for multiprecision types the cost of thread synchronisation is negligible, 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:
 
diff --git a/libs/math/doc/sf/polygamma.qbk b/libs/math/doc/sf/polygamma.qbk
new file mode 100644
index 000000000..19ff6a712
--- /dev/null
+++ b/libs/math/doc/sf/polygamma.qbk
@@ -0,0 +1,115 @@
+[section:polygamma Polygamma]
+
+[h4 Synopsis]
+
+``
+#include <boost/math/special_functions/polygamma.hpp>
+``
+
+  namespace boost{ namespace math{
+  
+  template <class T>
+  ``__sf_result`` polygamma(int n, T z);
+  
+  template <class T, class ``__Policy``>
+  ``__sf_result`` polygamma(int n, T z, const ``__Policy``&);
+  
+  }} // namespaces
+  
+[h4 Description]
+
+Returns the polygamma function of /x/. Polygamma is defined as the n'th
+derivative of the digamma function:
+
+[equation polygamma1]
+
+The following graphs illustrate the behaviour of the function for odd and even order:
+
+[graph polygamma2]
+[graph polygamma3]
+
+[optional_policy]
+
+The return type of this function is computed using the __arg_pomotion_rules:
+the result is of type `double` when T is an integer type, and type T otherwise.
+
+[h4 Accuracy]
+
+The following table shows the peak errors (in units of epsilon) 
+found on various platforms with various floating point types.
+Unless otherwise specified any floating point type that is narrower
+than the one shown will have __zero_error.
+
+[table
+[[Significand Size] [Platform and Compiler] [Small-medium positive arguments] [Small-medium negative x] ]
+[[53]            [Win32 Visual C++ 12]  [Peak=5.0 Mean=1] [Peak=1200 Mean=65]]
+[[64]            [Win64 Mingw GCC]  [Peak=16 Mean=3] [Peak=33 Mean=3] ]
+[[113]           [Win64 Mingw GCC __float128]    [Peak=6.5 Mean=1][Peak=30 Mean=4] ]
+]
+
+As shown above, error rates are generally very acceptable for moderately sized
+arguments.  Error rates should stay low for exact inputs, however, please note that the
+function becomes exceptionally sensitive to small changes in input for large n and negative x,
+indeed for cases where ['n!] would overflow, the function changes directly from -[infin] to
++[infin] somewhere between each negative integer - ['these cases are not handled correctly].
+
+[*For these reasons results should be treated with extreme caution when /n/ is large and x negative].
+
+[h4 Testing]
+
+Testing is against Mathematica generated spot values to 35 digit precision.
+
+[h4 Implementation]
+
+For x < 0 the following reflection formula is used:
+
+[equation polygamma2]
+
+The n'th derivative of ['cot(x)] is tabulated for small /n/, and for larger n
+has the general form:
+
+[equation polygamma3]
+
+The coefficients of the cosine terms can be calculated iteratively starting
+from ['C[sub 1,0] = -1] and then using
+
+[equation polygamma7]
+
+to generate coefficients for n+1.
+
+Note that every other coefficient is zero, and therefore what we have are
+even or odd polynomials depending on whether n is even or odd.
+
+Once x is positive then we have two methods available to us, for small x
+we use the series expansion:
+
+[equation polygamma4]
+
+Note that the evaluation of zeta functions at integer values is essentially a table lookup
+as __zeta is optimized for those cases.
+
+For large x we use the asymptotic expansion:
+
+[equation polygamma5]
+
+For x in-between the two extremes we use the relation:
+
+[equation polygamma6]
+
+to make x large enough for the asymptotic expansion to be used.
+
+There are also two special cases:
+
+[equation polygamma8]
+
+[equation polygamma9]
+
+[endsect][/section:polygamma The Polygamma Function]
+
+[/ 
+  Copyright 2014 John Maddock and Paul A. Bristow.
+  Distributed under 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).
+]
+
diff --git a/libs/math/doc/sf/trigamma.qbk b/libs/math/doc/sf/trigamma.qbk
new file mode 100644
index 000000000..827c8fb52
--- /dev/null
+++ b/libs/math/doc/sf/trigamma.qbk
@@ -0,0 +1,84 @@
+[section:trigamma Trigamma]
+
+[h4 Synopsis]
+
+``
+#include <boost/math/special_functions/trigamma.hpp>
+``
+
+  namespace boost{ namespace math{
+  
+  template <class T>
+  ``__sf_result`` trigamma(T z);
+  
+  template <class T, class ``__Policy``>
+  ``__sf_result`` trigamma(T z, const ``__Policy``&);
+  
+  }} // namespaces
+  
+[h4 Description]
+
+Returns the trigamma function of /x/. Trigamma is defined as the 
+derivative of the digamma function:
+
+[equation trigamma1]
+
+[graph trigamma]
+
+[optional_policy]
+
+The return type of this function is computed using the __arg_pomotion_rules:
+the result is of type `double` when T is an integer type, and type T otherwise.
+
+[h4 Accuracy]
+
+The following table shows the peak errors (in units of epsilon) 
+found on various platforms with various floating point types.
+Unless otherwise specified any floating point type that is narrower
+than the one shown will have __zero_error.
+
+[table
+[[Significand Size] [Platform and Compiler] [Random Values] ]
+[[53]            [Win32 Visual C++ 12]  [Peak=1.0 Mean=0.4] ]
+[[64]            [Win64 Mingw GCC]  [Peak=1.4 Mean=0.4]  ]
+[[113]           [Win64 Mingw GCC __float128]    [Peak=1.0 Mean=0.5] ]
+]
+
+As shown above, error rates are generally very low for built in types.
+For multiprecision types, error rates are typically in the order of a
+few epsilon.
+
+[h4 Testing]
+
+Testing is against Mathematica generated spot values to 35 digit precision.
+
+[h4 Implementation]
+
+The arbitrary precision version of this function simply calls __polygamma.
+
+For built in fixed precision types, negative arguments are first made positive via:
+
+[equation trigamma2]
+
+Then arguments in the range \[0, 1) are shifted to >= 1 via:
+
+[equation trigamma3]
+
+Then evaluation is via one of a number of rational approximations, for small x these are
+of the form:
+
+[equation trigamma4]
+
+and for large x of the form:
+
+[equation trigamma5]
+
+[endsect][/section:digamma The Trigamma Function]
+
+[/ 
+  Copyright 2014 John Maddock and Paul A. Bristow.
+  Distributed under 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).
+]
+
diff --git a/libs/math/doc/sf/zeta.qbk b/libs/math/doc/sf/zeta.qbk
index 08186c6a4..b8520685c 100644
--- a/libs/math/doc/sf/zeta.qbk
+++ b/libs/math/doc/sf/zeta.qbk
@@ -112,6 +112,18 @@ required for R(z-n) is not full machine precision, but an absolute error
 of: [epsilon]/R(0).  This saves us quite a few digits when dealing with large 
 z, especially when [epsilon] is small.
 
+Finally, there are some special cases for integer arguments, there are
+closed forms for negative or even integers:
+
+[equation zeta7]
+
+[equation zeta8]
+
+[equation zeta9]
+
+and for positive odd integers we simply cache pre-computed values as these are of great
+benefit to some infinite series calculations.
+
 [endsect]
 [/ :error_function The Error Functions]