Boost logo

Boost-Commit :

From: john_at_[hidden]
Date: 2007-07-28 11:13:53


Author: johnmaddock
Date: 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
New Revision: 7572
URL: http://svn.boost.org/trac/boost/changeset/7572

Log:
Updated function prototypes to include reference to policies.
Added TODO's where more substantive changes are required.
Removed the draft warnings.

Text files modified:
   sandbox/math_toolkit/libs/math/doc/bessel_ik.qbk | 8 ++
   sandbox/math_toolkit/libs/math/doc/bessel_jy.qbk | 8 ++
   sandbox/math_toolkit/libs/math/doc/bessel_spherical.qbk | 7 ++
   sandbox/math_toolkit/libs/math/doc/beta.qbk | 7 +
   sandbox/math_toolkit/libs/math/doc/beta_derivative.qbk | 7 +
   sandbox/math_toolkit/libs/math/doc/compilers.qbk | 6 +
   sandbox/math_toolkit/libs/math/doc/concepts.qbk | 2
   sandbox/math_toolkit/libs/math/doc/digamma.qbk | 7 +
   sandbox/math_toolkit/libs/math/doc/ellint_carlson.qbk | 34 +++++++++++
   sandbox/math_toolkit/libs/math/doc/ellint_legendre.qbk | 62 +++++++++++++++++++---
   sandbox/math_toolkit/libs/math/doc/erf.qbk | 16 +++++
   sandbox/math_toolkit/libs/math/doc/erf_inv.qbk | 16 +++++
   sandbox/math_toolkit/libs/math/doc/error.qbk | 2
   sandbox/math_toolkit/libs/math/doc/error_handling.qbk | 4 +
   sandbox/math_toolkit/libs/math/doc/factorials.qbk | 28 ++++++++++
   sandbox/math_toolkit/libs/math/doc/fraction.qbk | 2
   sandbox/math_toolkit/libs/math/doc/gamma_derivatives.qbk | 7 +
   sandbox/math_toolkit/libs/math/doc/gamma_ratios.qbk | 20 +++++-
   sandbox/math_toolkit/libs/math/doc/hermite.qbk | 10 ++
   sandbox/math_toolkit/libs/math/doc/ibeta.qbk | 28 +++++++++
   sandbox/math_toolkit/libs/math/doc/ibeta_inv.qbk | 77 ++++++++++++++++++++++++++--
   sandbox/math_toolkit/libs/math/doc/igamma.qbk | 28 +++++++++
   sandbox/math_toolkit/libs/math/doc/igamma_inv.qbk | 28 +++++++++
   sandbox/math_toolkit/libs/math/doc/implementation.qbk | 2
   sandbox/math_toolkit/libs/math/doc/inv_hyper.qbk | 21 ++++++-
   sandbox/math_toolkit/libs/math/doc/laguerre.qbk | 16 +++++
   sandbox/math_toolkit/libs/math/doc/lanczos.qbk | 2
   sandbox/math_toolkit/libs/math/doc/legendre.qbk | 22 +++++++
   sandbox/math_toolkit/libs/math/doc/lgamma.qbk | 10 ++
   sandbox/math_toolkit/libs/math/doc/math.qbk | 6 ++
   sandbox/math_toolkit/libs/math/doc/minima.qbk | 2
   sandbox/math_toolkit/libs/math/doc/polynomial.qbk | 2
   sandbox/math_toolkit/libs/math/doc/powers.qbk | 32 +++++++++++
   sandbox/math_toolkit/libs/math/doc/rational.qbk | 5 +
   sandbox/math_toolkit/libs/math/doc/relative_error.qbk | 2
   sandbox/math_toolkit/libs/math/doc/roadmap.qbk | 18 ++++++
   sandbox/math_toolkit/libs/math/doc/roots.qbk | 2
   sandbox/math_toolkit/libs/math/doc/roots_without_derivatives.qbk | 107 ++++++++++++++++++++++++++++++++++++++-
   sandbox/math_toolkit/libs/math/doc/series.qbk | 2
   sandbox/math_toolkit/libs/math/doc/sinc.qbk | 21 ++++++-
   sandbox/math_toolkit/libs/math/doc/spherical_harmonic.qbk | 22 +++++++
   sandbox/math_toolkit/libs/math/doc/structure.qbk | 3 +
   sandbox/math_toolkit/libs/math/doc/tgamma.qbk | 18 +++++
   43 files changed, 646 insertions(+), 83 deletions(-)

Modified: sandbox/math_toolkit/libs/math/doc/bessel_ik.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/bessel_ik.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/bessel_ik.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -6,9 +6,15 @@
    template <class T1, class T2>
    ``__sf_result`` cyl_bessel_i(T1 v, T2 x);
 
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` cyl_bessel_i(T1 v, T2 x, const ``__Policy``&);
+
    template <class T1, class T2>
    ``__sf_result`` cyl_bessel_k(T1 v, T2 x);
    
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` cyl_bessel_k(T1 v, T2 x, const ``__Policy``&);
+
    
 [h4 Description]
 
@@ -29,6 +35,8 @@
 when T1 and T2 are different types. The functions are also optimised for the
 relatively common case that T1 is an integer.
 
+[optional_policy]
+
 The functions return the result of __domain_error whenever the result is
 undefined or complex. For __cyl_bessel_j this occurs when `x < 0` and v is not
 an integer, or when `x == 0` and `v != 0`. For __cyl_neumann this occurs

Modified: sandbox/math_toolkit/libs/math/doc/bessel_jy.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/bessel_jy.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/bessel_jy.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -6,9 +6,15 @@
    template <class T1, class T2>
    ``__sf_result`` cyl_bessel_j(T1 v, T2 x);
 
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` cyl_bessel_j(T1 v, T2 x, const ``__Policy``&);
+
    template <class T1, class T2>
    ``__sf_result`` cyl_neumann(T1 v, T2 x);
    
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` cyl_neumann(T1 v, T2 x, const ``__Policy``&);
+
    
 [h4 Description]
 
@@ -29,6 +35,8 @@
 when T1 and T2 are different types. The functions are also optimised for the
 relatively common case that T1 is an integer.
 
+[optional_policy]
+
 The functions return the result of __domain_error whenever the result is
 undefined or complex. For __cyl_bessel_j this occurs when `x < 0` and v is not
 an integer, or when `x == 0` and `v != 0`. For __cyl_neumann this occurs

Modified: sandbox/math_toolkit/libs/math/doc/bessel_spherical.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/bessel_spherical.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/bessel_spherical.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -6,9 +6,14 @@
    template <class T1, class T2>
    ``__sf_result`` sph_bessel(unsigned v, T2 x);
 
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` sph_bessel(unsigned v, T2 x, const ``__Policy``&);
+
    template <class T1, class T2>
    ``__sf_result`` sph_neumann(unsigned v, T2 x);
    
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` sph_neumann(unsigned v, T2 x, const ``__Policy``&);
    
 [h4 Description]
 
@@ -26,6 +31,8 @@
 The return type of these functions is computed using the __arg_pomotion_rules
 for the single argument type T.
 
+[optional_policy]
+
 The functions return the result of __domain_error whenever the result is
 undefined or complex: this occurs when `x < 0`.
 

Modified: sandbox/math_toolkit/libs/math/doc/beta.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/beta.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/beta.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:beta_function Beta]
 
-[caution __caution ]
-
 [h4 Synopsis]
 
 ``
@@ -13,6 +11,9 @@
    template <class T1, class T2>
    ``__sf_result`` beta(T1 a, T2 b);
    
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` beta(T1 a, T2 b, const ``__Policy``&);
+
    }} // namespaces
 
 [h4 Description]
@@ -23,6 +24,8 @@
 
 [$../graphs/beta.png]
 
+[optional_policy]
+
 There are effectively two versions of this function internally: a fully
 generic version that is slow, but reasonably accurate, and a much more
 efficient approximation that is used where the number of digits in the significand

Modified: sandbox/math_toolkit/libs/math/doc/beta_derivative.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/beta_derivative.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/beta_derivative.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:beta_derivative Derivative of the Incomplete Beta Function]
 
-[caution __caution ]
-
 [h4 Synopsis]
 
 ``
@@ -13,6 +11,9 @@
    template <class T1, class T2, class T3>
    ``__sf_result`` ibeta_derivative(T1 a, T2 b, T3 x);
    
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` ibeta_derivative(T1 a, T2 b, T3 x, const ``__Policy``&);
+
    }} // namespaces
    
 [h4 Description]
@@ -26,6 +27,8 @@
 The return type of this function is computed using the __arg_pomotion_rules
 when T1, T2 and T3 are different types.
 
+[optional_policy]
+
 [h4 Accuracy]
 
 Almost identical to the incomplete beta function __ibeta.

Modified: sandbox/math_toolkit/libs/math/doc/compilers.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/compilers.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/compilers.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -15,7 +15,11 @@
       with that option), otherwise our headers should be warning free.]]
 [[HP aCC][HP-UX][ Unfortunately this emits quite a few warnings from libraries
       upon which we depend (TR1, Array etc).]]
-[[Borland 5.8.2][Windows][Almost works (doesn't compile the test suite, but does cope with *our* code I believe)]]
+[[Borland 5.8.2][Windows][Almost works: some effort has been put into porting to this compiler.
+ However, during testing a number of instances were encountered where this compiler
+ generated incorrect code: completely omitting a function call seemingly at random.
+ For this reason, we cannot recoment using this library with this compiler, as the
+ correct operation of the code cannot be guarenteed.]]
 [[MSVC 8.0][Windows][Warning free at level 4]]
 ]
 

Modified: sandbox/math_toolkit/libs/math/doc/concepts.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/concepts.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/concepts.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -224,6 +224,8 @@
 [[Expression][Result Type][Notes]]
 [[DistributionType::value_type][RealType]
       [The real-number type /RealType/ upon which the distribution operates.]]
+[[DistributionType::policy_type][RealType]
+ [The __Policy to use when evaluating functions that depend on this distribution.]]
 [[d = cd][Distribution&][Distribution types are assignable.]]
 [[Distribution(cd)][Distribution][Distribution types are copy constructible.]]
 [[pdf(cd, cr)][RealType][Returns the PDF of the distribution.]]

Modified: sandbox/math_toolkit/libs/math/doc/digamma.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/digamma.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/digamma.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:digamma Digamma]
 
-[caution __caution ]
-
 [h4 Synopsis]
 
 ``
@@ -13,6 +11,9 @@
   template <class T>
   ``__sf_result`` digamma(T z);
   
+ template <class T, class ``__Policy``>
+ ``__sf_result`` digamma(T z, const ``__Policy``&);
+
   }} // namespaces
   
 [h4 Description]
@@ -24,6 +25,8 @@
 
 [$../graphs/digamma.png]
 
+[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.

Modified: sandbox/math_toolkit/libs/math/doc/ellint_carlson.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/ellint_carlson.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/ellint_carlson.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -8,8 +8,6 @@
 
 [section:ellint_carlson Elliptic Integrals - Carlson Form]
 
-[caution __caution ]
-
 [heading Synopsis]
 
 ``
@@ -21,6 +19,9 @@
   template <class T1, class T2, class T3>
   ``__sf_result`` ellint_rf(T1 x, T2 y, T3 z)
 
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` ellint_rf(T1 x, T2 y, T3 z, const ``__Policy``&)
+
   }} // namespaces
 
 
@@ -33,6 +34,9 @@
   template <class T1, class T2, class T3>
   ``__sf_result`` ellint_rd(T1 x, T2 y, T3 z)
 
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` ellint_rd(T1 x, T2 y, T3 z, const ``__Policy``&)
+
   }} // namespaces
 
 
@@ -45,6 +49,9 @@
   template <class T1, class T2, class T3, class T4>
   ``__sf_result`` ellint_rj(T1 x, T2 y, T3 z, T4 p)
 
+ template <class T1, class T2, class T3, class T4, class ``__Policy``>
+ ``__sf_result`` ellint_rj(T1 x, T2 y, T3 z, T4 p, const ``__Policy``&)
+
   }} // namespaces
 
 
@@ -57,6 +64,9 @@
   template <class T1, class T2>
   ``__sf_result`` ellint_rc(T1 x, T2 y)
 
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` ellint_rc(T1 x, T2 y, const ``__Policy``&)
+
   }} // namespaces
 
 
@@ -75,6 +85,9 @@
   template <class T1, class T2, class T3>
   ``__sf_result`` ellint_rf(T1 x, T2 y, T3 z)
   
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` ellint_rf(T1 x, T2 y, T3 z, const ``__Policy``&)
+
 Returns Carlson's Elliptic Integral R[sub F]:
 
 [$../equations/ellint9.png]
@@ -82,9 +95,14 @@
 Requires that all of the arguments are non-negative, and at most
 one may be zero. Otherwise returns the result of __domain_error.
 
+[optional_policy]
+
   template <class T1, class T2, class T3>
   ``__sf_result`` ellint_rd(T1 x, T2 y, T3 z)
   
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` ellint_rd(T1 x, T2 y, T3 z, const ``__Policy``&)
+
 Returns Carlson's elliptic integral R[sub D]:
 
 [$../equations/ellint10.png]
@@ -92,9 +110,14 @@
 Requires that x and y are non-negative, with at most one of them
 zero, and that z >= 0. Otherwise returns the result of __domain_error.
 
+[optional_policy]
+
   template <class T1, class T2, class T3, class T4>
   ``__sf_result`` ellint_rj(T1 x, T2 y, T3 z, T4 p)
 
+ template <class T1, class T2, class T3, class T4, class ``__Policy``>
+ ``__sf_result`` ellint_rj(T1 x, T2 y, T3 z, T4 p, const ``__Policy``&)
+
 Returns Carlson's elliptic integral R[sub J]:
   
 [$../equations/ellint11.png]
@@ -102,6 +125,8 @@
 Requires that x, y and z are non-negative, with at most one of them
 zero, and that ['p != 0]. Otherwise returns the result of __domain_error.
 
+[optional_policy]
+
 When ['p < 0] the function returns the
 [@http://en.wikipedia.org/wiki/Cauchy_principal_value Cauchy principal value]
 using the relation:
@@ -111,6 +136,9 @@
   template <class T1, class T2>
   ``__sf_result`` ellint_rc(T1 x, T2 y)
 
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` ellint_rc(T1 x, T2 y, const ``__Policy``&)
+
 Returns Carlson's elliptic integral R[sub C]:
   
 [$../equations/ellint12.png]
@@ -118,6 +146,8 @@
 Requires that ['x > 0] and that ['y != 0].
 Otherwise returns the result of __domain_error.
 
+[optional_policy]
+
 When ['y < 0] the function returns the
 [@http://mathworld.wolfram.com/CauchyPrincipalValue.html Cauchy principal value]
 using the relation:

Modified: sandbox/math_toolkit/libs/math/doc/ellint_legendre.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/ellint_legendre.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/ellint_legendre.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -8,8 +8,6 @@
 
 [section:ellint_1 Elliptic Integrals of the First Kind - Legendre Form]
 
-[caution __caution]
-
 [heading Synopsis]
 
 ``
@@ -21,9 +19,15 @@
   template <class T1, class T2>
   ``__sf_result`` ellint_1(T1 k, T2 phi);
 
- template <typename T>
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` ellint_1(T1 k, T2 phi, const ``__Policy``&);
+
+ template <class T>
   ``__sf_result`` ellint_1(T k);
 
+ template <class T, class ``__Policy``>
+ ``__sf_result`` ellint_1(T k, const ``__Policy``&);
+
   }} // namespaces
   
 [heading Description]
@@ -40,21 +44,31 @@
   template <class T1, class T2>
   ``__sf_result`` ellint_1(T1 k, T2 phi);
   
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` ellint_1(T1 k, T2 phi, const ``__Policy``&);
+
 Returns the incomplete elliptic integral of the first kind ['F([phi], k)]:
 
 [$../equations/ellint2.png]
 
 Requires -1 <= k <= 1, otherwise returns the result of __domain_error.
 
- template <typename T>
+[optional_policy]
+
+ template <class T>
   ``__sf_result`` ellint_1(T k);
   
+ template <class T>
+ ``__sf_result`` ellint_1(T k, const ``__Policy``&);
+
 Returns the complete elliptic integral of the first kind ['K(k)]:
 
 [$../equations/ellint6.png]
 
 Requires -1 <= k <= 1, otherwise returns the result of __domain_error.
 
+[optional_policy]
+
 [heading Accuracy]
 
 These functions are computed using only basic arithmetic operations, so
@@ -94,8 +108,6 @@
 
 [section:ellint_2 Elliptic Integrals of the Second Kind - Legendre Form]
 
-[caution __caution]
-
 [heading Synopsis]
 
 ``
@@ -107,9 +119,15 @@
   template <class T1, class T2>
   ``__sf_result`` ellint_2(T1 k, T2 phi);
 
- template <typename T>
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` ellint_2(T1 k, T2 phi, const ``__Policy``&);
+
+ template <class T>
   ``__sf_result`` ellint_2(T k);
 
+ template <class T, class ``__Policy``>
+ ``__sf_result`` ellint_2(T k, const ``__Policy``&);
+
   }} // namespaces
   
 [heading Description]
@@ -126,21 +144,31 @@
   template <class T1, class T2>
   ``__sf_result`` ellint_2(T1 k, T2 phi);
   
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` ellint_2(T1 k, T2 phi, const ``__Policy``&);
+
 Returns the incomplete elliptic integral of the second kind ['E([phi], k)]:
 
 [$../equations/ellint3.png]
 
 Requires -1 <= k <= 1, otherwise returns the result of __domain_error.
 
- template <typename T>
+[optional_policy]
+
+ template <class T>
   ``__sf_result`` ellint_2(T k);
   
+ template <class T>
+ ``__sf_result`` ellint_2(T k, const ``__Policy``&);
+
 Returns the complete elliptic integral of the first kind ['E(k)]:
 
 [$../equations/ellint7.png]
 
 Requires -1 <= k <= 1, otherwise returns the result of __domain_error.
 
+[optional_policy]
+
 [heading Accuracy]
 
 These functions are computed using only basic arithmetic operations, so
@@ -180,8 +208,6 @@
 
 [section:ellint_3 Elliptic Integrals of the Third Kind - Legendre Form]
 
-[caution __caution]
-
 [heading Synopsis]
 
 ``
@@ -193,9 +219,15 @@
   template <class T1, class T2, class T3>
   ``__sf_result`` ellint_3(T1 k, T2 n, T3 phi);
 
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` ellint_3(T1 k, T2 n, T3 phi, const ``__Policy``&);
+
   template <class T1, class T2>
   ``__sf_result`` ellint_3(T1 k, T2 n);
 
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` ellint_3(T1 k, T2 n, const ``__Policy``&);
+
   }} // namespaces
   
 [heading Description]
@@ -212,6 +244,9 @@
   template <class T1, class T2, class T3>
   ``__sf_result`` ellint_3(T1 k, T2 n, T3 phi);
   
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` ellint_3(T1 k, T2 n, T3 phi, const ``__Policy``&);
+
 Returns the incomplete elliptic integral of the third kind ['[Pi](n, [phi], k)]:
 
 [$../equations/ellint4.png]
@@ -220,6 +255,8 @@
 returns the result of __domain_error (outside this range the result
 would be complex).
 
+[optional_policy]
+
 [caution In addition, the region where ['n > 1] and [phi][space] ['is not in the range]
 \[0, [pi]\/2\] is currently unsupported and returns the result of __domain_error.
 For this reason it is recomended that you keep [phi][space] inside its "natural" range
@@ -228,6 +265,9 @@
   template <class T1, class T2>
   ``__sf_result`` ellint_3(T1 k, T2 n);
   
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` ellint_3(T1 k, T2 n, const ``__Policy``&);
+
 Returns the complete elliptic integral of the first kind ['[Pi](n, k)]:
 
 [$../equations/ellint8.png]
@@ -235,6 +275,8 @@
 Requires ['-1 <= k <= 1] and ['n < 1], otherwise returns the
 result of __domain_error (outside this range the result would be complex).
 
+[opitonal_policy]
+
 [heading Accuracy]
 
 These functions are computed using only basic arithmetic operations, so

Modified: sandbox/math_toolkit/libs/math/doc/erf.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/erf.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/erf.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:error_function Error Functions]
 
-[caution __caution ]
-
 [h4 Synopsis]
 
 ``
@@ -13,19 +11,30 @@
    template <class T>
    ``__sf_result`` erf(T z);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` erf(T z, const ``__Policy``&);
+
    template <class T>
    ``__sf_result`` erfc(T z);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` erfc(T z, const ``__Policy``&);
+
    }} // namespaces
    
 The return type of these functions is computed using the __arg_pomotion_rules:
 the return type is `double` if T is an integer type, and T otherwise.
 
+[optional_policy]
+
 [h4 Description]
 
    template <class T>
    ``__sf_result`` erf(T z);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` erf(T z, const ``__Policy``&);
+
 Returns the [@http://en.wikipedia.org/wiki/Error_function error function]
 [@http://functions.wolfram.com/GammaBetaErf/Erf/ erf] of z:
 
@@ -36,6 +45,9 @@
    template <class T>
    ``__sf_result`` erfc(T z);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` erfc(T z, const ``__Policy``&);
+
 Returns the complement of the [@http://functions.wolfram.com/GammaBetaErf/Erfc/ error function] of z:
 
 [$../equations/erf2.png]

Modified: sandbox/math_toolkit/libs/math/doc/erf_inv.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/erf_inv.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/erf_inv.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:error_inv Error Function Inverses]
 
-[caution __caution ]
-
 [h4 Synopsis]
 
 ``
@@ -13,19 +11,30 @@
    template <class T>
    ``__sf_result`` erf_inv(T p);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` erf_inv(T p, const ``__Policy``&);
+
    template <class T>
    ``__sf_result`` erfc_inv(T p);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` erfc_inv(T p, const ``__Policy``&);
+
    }} // namespaces
    
 The return type of these functions is computed using the __arg_pomotion_rules:
 the return type is `double` if T is an integer type, and T otherwise.
 
+[optional_policy]
+
 [h4 Description]
 
    template <class T>
    ``__sf_result`` erf_inv(T z);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` erf_inv(T z, const ``__Policy``&);
+
 Returns the [@http://functions.wolfram.com/GammaBetaErf/InverseErf/ inverse error function]
 of z, that is a value x such that:
 
@@ -36,6 +45,9 @@
    template <class T>
    ``__sf_result`` erfc_inv(T z);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` erfc_inv(T z, const ``__Policy``&);
+
 Returns the inverse of the complement of the error function of z, that is a
 value x such that:
 

Modified: sandbox/math_toolkit/libs/math/doc/error.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/error.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/error.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:relative_error Relative Error]
 
-[caution __caution ]
-
 Given an actual value /a/ and a found value /v/ the relative error can be
 calculated from:
 

Modified: sandbox/math_toolkit/libs/math/doc/error_handling.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/error_handling.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/error_handling.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,6 +1,8 @@
 [section:error_handling Error Handling]
 
-[caution __caution ]
+[caution This documentation is out of date.
+
+TODO: rewrite me please!! ]
 
 [def __format [@../../libs/format/index.html Boost.Format]]
 

Modified: sandbox/math_toolkit/libs/math/doc/factorials.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/factorials.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/factorials.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -13,6 +13,9 @@
    template <class T>
    T factorial(unsigned i);
    
+ template <class T, class ``__Policy``>
+ T factorial(unsigned i, const ``__Policy``&);
+
    template <class T>
    T unchecked_factorial(unsigned i);
    
@@ -26,8 +29,13 @@
    template <class T>
    T factorial(unsigned i);
 
+ template <class T, class ``__Policy``>
+ T factorial(unsigned i, const ``__Policy``&);
+
 Returns [^i!].
 
+[optional_policy]
+
 For [^i <= max_factorial<T>::value] this is implemented by table lookup,
 for larger values of [^i], this function is implemented in terms of __tgamma.
 
@@ -86,10 +94,15 @@
    template <class T>
    T double_factorial(unsigned i);
    
+ template <class T, class ``__Policy``>
+ T double_factorial(unsigned i, const ``__Policy``&);
+
    }} // namespaces
 
 Returns [^i!!].
 
+[optional_policy]
+
 May return the result of __overflow_error if the result is too large
 to represent in type T. The implementation is designed to be optimised
 for small /i/ where table lookup of i! is possible.
@@ -131,6 +144,9 @@
    template <class T>
    ``__sf_result`` rising_factorial(T x, int i);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` rising_factorial(T x, int i, const ``__Policy``&);
+
    }} // namespaces
 
 Returns the rising factorial of /x/ and /i/:
@@ -143,6 +159,8 @@
                           
 Note that both /x/ and /i/ can be negative as well as positive.
 
+[optional_policy]
+
 May return the result of __overflow_error if the result is too large
 to represent in type T.
 
@@ -179,6 +197,9 @@
    template <class T>
    ``__sf_result`` falling_factorial(T x, unsigned i);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` falling_factorial(T x, unsigned i, const ``__Policy``&);
+
    }} // namespaces
 
 Returns the falling factorial of /x/ and /i/:
@@ -189,6 +210,8 @@
 `unsigned` second argument. Argument /x/ can be either positive or
 negative however.
 
+[optional_policy]
+
 May return the result of __overflow_error if the result is too large
 to represent in type T.
 
@@ -225,12 +248,17 @@
    template <class T>
    T binomial_coefficient(unsigned n, unsigned k);
 
+ template <class T, class ``__Policy``>
+ T binomial_coefficient(unsigned n, unsigned k, const ``__Policy``&);
+
    }} // namespaces
 
 Returns the binomial coefficient: [sub n]C[sub k].
 
 Requires k <= n.
 
+[optional_policy]
+
 May return the result of __overflow_error if the result is too large
 to represent in type T.
    

Modified: sandbox/math_toolkit/libs/math/doc/fraction.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/fraction.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/fraction.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:cf Continued Fraction Evaluation]
 
-[caution __caution ]
-
 [h4 Synopsis]
 
 ``

Modified: sandbox/math_toolkit/libs/math/doc/gamma_derivatives.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/gamma_derivatives.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/gamma_derivatives.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:gamma_derivatives Derivative of the Incomplete Gamma Function]
 
-[caution __caution ]
-
 [h4 Synopsis]
 
 ``
@@ -13,6 +11,9 @@
    template <class T1, class T2>
    ``__sf_result`` gamma_p_derivative(T1 a, T2 x);
    
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` gamma_p_derivative(T1 a, T2 x, const ``__Policy``&);
+
    }} // namespaces
    
 [h4 Description]
@@ -23,6 +24,8 @@
 
 [$../equations/derivative1.png]
 
+[optional_policy]
+
 Note that the derivative of the function __gamma_q can be obtained by negating
 the result of this function.
 

Modified: sandbox/math_toolkit/libs/math/doc/gamma_ratios.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/gamma_ratios.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/gamma_ratios.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,9 +1,5 @@
 [section:gamma_ratios Ratios of Gamma Functions]
 
-[caution __caution ]
-
-[h4 Synopsis]
-
 ``
 #include <boost/math/special_functions/gamma.hpp>
 ``
@@ -13,9 +9,15 @@
    template <class T1, class T2>
    ``__sf_result`` tgamma_ratio(T1 a, T2 b);
    
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` tgamma_ratio(T1 a, T2 b, const ``__Policy``&);
+
    template <class T1, class T2>
    ``__sf_result`` tgamma_delta_ratio(T1 a, T2 delta);
    
+ template <class T1, class T2, class Policy>
+ ``__sf_result`` tgamma_delta_ratio(T1 a, T2 delta, const ``__Policy``&);
+
    }} // namespaces
    
 [h4 Description]
@@ -23,19 +25,29 @@
    template <class T1, class T2>
    ``__sf_result`` tgamma_ratio(T1 a, T2 b);
    
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` tgamma_ratio(T1 a, T2 b, const ``__Policy``&);
+
 Returns the ratio of gamma functions:
 
 [$../equations/gamma_ratio0.png]
 
+[optional_policy]
+
 Internally this just calls `tgamma_delta_ratio(a, b-a)`.
    
    template <class T1, class T2>
    ``__sf_result`` tgamma_delta_ratio(T1 a, T2 delta);
    
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` tgamma_delta_ratio(T1 a, T2 delta, const ``__Policy``&);
+
 Returns the ratio of gamma functions:
 
 [$../equations/gamma_ratio1.png]
 
+[optional_policy]
+
 Note that the result is calculated accurately even when /delta/ is
 small compared to /a/: indeed even if /a+delta ~ a/. The function is
 typically used when /a/ is large and /delta/ is very small.

Modified: sandbox/math_toolkit/libs/math/doc/hermite.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/hermite.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/hermite.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:hermite Hermite Polynomials]
 
-[caution __caution ]
-
 [h4 Synopsis]
 
 ``
@@ -13,6 +11,9 @@
    template <class T>
    ``__sf_result`` hermite(unsigned n, T x);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` hermite(unsigned n, T x, const ``__Policy``&);
+
    template <class T1, class T2, class T3>
    ``__sf_result`` hermite_next(unsigned n, T1 x, T2 Hn, T3 Hnm1);
       
@@ -27,10 +28,15 @@
    template <class T>
    ``__sf_result`` hermite(unsigned n, T x);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` hermite(unsigned n, T x, const ``__Policy``&);
+
 Returns the value of the Hermite Polynomial of order /n/ at point /x/:
 
 [$../equations/hermite_0.png]
 
+[optional_policy]
+
 The following graph illustrates the behaviour of the first few
 Hermite Polynomials:
 

Modified: sandbox/math_toolkit/libs/math/doc/ibeta.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/ibeta.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/ibeta.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:ibeta_function Incomplete Beta Functions]
 
-[caution __caution ]
-
 [h4 Synopsis]
 
 ``
@@ -13,15 +11,27 @@
    template <class T1, class T2, class T3>
    ``__sf_result`` ibeta(T1 a, T2 b, T3 x);
    
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` ibeta(T1 a, T2 b, T3 x, const ``__Policy``&);
+
    template <class T1, class T2, class T3>
    ``__sf_result`` ibetac(T1 a, T2 b, T3 x);
    
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` ibetac(T1 a, T2 b, T3 x, const ``__Policy``&);
+
    template <class T1, class T2, class T3>
    ``__sf_result`` beta(T1 a, T2 b, T3 x);
    
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` beta(T1 a, T2 b, T3 x, const ``__Policy``&);
+
    template <class T1, class T2, class T3>
    ``__sf_result`` betac(T1 a, T2 b, T3 x);
    
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` betac(T1 a, T2 b, T3 x, const ``__Policy``&);
+
    }} // namespaces
    
 [h4 Description]
@@ -39,9 +49,14 @@
 The return type of these functions is computed using the __arg_pomotion_rules
 when T1, T2 and T3 are different types.
 
+[optional_policy]
+
    template <class T1, class T2, class T3>
    ``__sf_result`` ibeta(T1 a, T2 b, T3 x);
 
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` ibeta(T1 a, T2 b, T3 x, const ``__Policy``&);
+
 Returns the normalised incomplete beta function of a, b and x:
 
 [$../equations/ibeta3.png]
@@ -51,6 +66,9 @@
    template <class T1, class T2, class T3>
    ``__sf_result`` ibetac(T1 a, T2 b, T3 x);
    
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` ibetac(T1 a, T2 b, T3 x, const ``__Policy``&);
+
 Returns the normalised complement of the incomplete beta function of a, b and x:
 
 [$../equations/ibeta4.png]
@@ -58,6 +76,9 @@
    template <class T1, class T2, class T3>
    ``__sf_result`` beta(T1 a, T2 b, T3 x);
 
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` beta(T1 a, T2 b, T3 x, const ``__Policy``&);
+
 Returns the full (non-normalised) incomplete beta function of a, b and x:
 
 [$../equations/ibeta1.png]
@@ -65,6 +86,9 @@
    template <class T1, class T2, class T3>
    ``__sf_result`` betac(T1 a, T2 b, T3 x);
 
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` betac(T1 a, T2 b, T3 x, const ``__Policy``&);
+
 Returns the full (non-normalised) complement of the incomplete beta function of a, b and x:
 
 [$../equations/ibeta2.png]

Modified: sandbox/math_toolkit/libs/math/doc/ibeta_inv.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/ibeta_inv.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/ibeta_inv.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,9 +1,5 @@
 [section:ibeta_inv_function The Incomplete Beta Function Inverses]
 
-[caution __caution ]
-
-[h4 Synopsis]
-
 ``
 #include <boost/math/special_functions/beta.hpp>
 ``
@@ -13,27 +9,51 @@
    template <class T1, class T2, class T3>
    ``__sf_result`` ibeta_inv(T1 a, T2 b, T3 p);
    
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` ibeta_inv(T1 a, T2 b, T3 p, const ``__Policy``&);
+
    template <class T1, class T2, class T3, class T4>
    ``__sf_result`` ibeta_inv(T1 a, T2 b, T3 p, T4* py);
    
+ template <class T1, class T2, class T3, class T4, class ``__Policy``>
+ ``__sf_result`` ibeta_inv(T1 a, T2 b, T3 p, T4* py, const ``__Policy``&);
+
    template <class T1, class T2, class T3>
    ``__sf_result`` ibetac_inv(T1 a, T2 b, T3 q);
    
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` ibetac_inv(T1 a, T2 b, T3 q, const ``__Policy``&);
+
    template <class T1, class T2, class T3, class T4>
    ``__sf_result`` ibetac_inv(T1 a, T2 b, T3 q, T4* py);
    
+ template <class T1, class T2, class T3, class T4, class ``__Policy``>
+ ``__sf_result`` ibetac_inv(T1 a, T2 b, T3 q, T4* py, const ``__Policy``&);
+
    template <class T1, class T2, class T3>
    ``__sf_result`` ibeta_inva(T1 b, T2 x, T3 p);
    
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` ibeta_inva(T1 b, T2 x, T3 p, const ``__Policy``&);
+
    template <class T1, class T2, class T3>
    ``__sf_result`` ibetac_inva(T1 b, T2 x, T3 q);
    
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` ibetac_inva(T1 b, T2 x, T3 q, const ``__Policy``&);
+
    template <class T1, class T2, class T3>
    ``__sf_result`` ibeta_invb(T1 a, T2 x, T3 p);
    
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` ibeta_invb(T1 a, T2 x, T3 p, const ``__Policy``&);
+
    template <class T1, class T2, class T3>
    ``__sf_result`` ibetac_invb(T1 a, T2 x, T3 q);
    
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` ibetac_invb(T1 a, T2 x, T3 q, const ``__Policy``&);
+
    }} // namespaces
    
 [h4 Description]
@@ -44,6 +64,8 @@
 any of the three parameters to the incomplete beta, starting from either
 the result of the incomplete beta (p) or its complement (q).
 
+[optional_policy]
+
 [tip When people normally talk about the inverse of the incomplete
 beta function, they are talking about inverting on parameter /x/.
 These are implemented here as ibeta_inv and ibeta_inv, and are by
@@ -61,9 +83,15 @@
    template <class T1, class T2, class T3>
    ``__sf_result`` ibeta_inv(T1 a, T2 b, T3 p);
    
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` ibeta_inv(T1 a, T2 b, T3 p, const ``__Policy``&);
+
    template <class T1, class T2, class T3, class T4>
    ``__sf_result`` ibeta_inv(T1 a, T2 b, T3 p, T4* py);
    
+ template <class T1, class T2, class T3, class T4, class ``__Policy``>
+ ``__sf_result`` ibeta_inv(T1 a, T2 b, T3 p, T4* py, const ``__Policy``&);
+
 Returns a value /x/ such that: `p = ibeta(a, b, x);`
 and sets `*py = 1 - x` when the `py` parameter is provided and is non-null.
 Note that internally this function computes whichever is the smaller of
@@ -73,12 +101,20 @@
 
 Requires: /a,b > 0/ and /0 <= p <= 1/.
 
+[optional_policy]
+
    template <class T1, class T2, class T3>
    ``__sf_result`` ibetac_inv(T1 a, T2 b, T3 q);
    
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` ibetac_inv(T1 a, T2 b, T3 q, const ``__Policy``&);
+
    template <class T1, class T2, class T3, class T4>
    ``__sf_result`` ibetac_inv(T1 a, T2 b, T3 q, T4* py);
    
+ template <class T1, class T2, class T3, class T4, class ``__Policy``>
+ ``__sf_result`` ibetac_inv(T1 a, T2 b, T3 q, T4* py, const ``__Policy``&);
+
 Returns a value /x/ such that: `q = ibetac(a, b, x);`
 and sets `*py = 1 - x` when the `py` parameter is provided and is non-null.
 Note that internally this function computes whichever is the smaller of
@@ -88,34 +124,56 @@
 
 Requires: /a,b > 0/ and /0 <= q <= 1/.
 
+[optional_policy]
+
    template <class T1, class T2, class T3>
    ``__sf_result`` ibeta_inva(T1 b, T2 x, T3 p);
 
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` ibeta_inva(T1 b, T2 x, T3 p, const ``__Policy``&);
+
 Returns a value /a/ such that: `p = ibeta(a, b, x);`
 
 Requires: /b > 0/, /0 < x < 1/ and /0 <= p <= 1/.
 
+[optional_policy]
+
    template <class T1, class T2, class T3>
    ``__sf_result`` ibetac_inva(T1 b, T2 x, T3 p);
    
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` ibetac_inva(T1 b, T2 x, T3 p, const ``__Policy``&);
+
 Returns a value /a/ such that: `q = ibetac(a, b, x);`
 
 Requires: /b > 0/, /0 < x < 1/ and /0 <= q <= 1/.
 
+[optional_policy]
+
    template <class T1, class T2, class T3>
    ``__sf_result`` ibeta_invb(T1 b, T2 x, T3 p);
+
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` ibeta_invb(T1 b, T2 x, T3 p, const ``__Policy``&);
 
 Returns a value /b/ such that: `p = ibeta(a, b, x);`
 
 Requires: /a > 0/, /0 < x < 1/ and /0 <= p <= 1/.
 
+[optional_policy]
+
    template <class T1, class T2, class T3>
    ``__sf_result`` ibetac_invb(T1 b, T2 x, T3 p);
    
+ template <class T1, class T2, class T3, class ``__Policy``>
+ ``__sf_result`` ibetac_invb(T1 b, T2 x, T3 p, const ``__Policy``&);
+
 Returns a value /b/ such that: `q = ibetac(a, b, x);`
 
 Requires: /a > 0/, /0 < x < 1/ and /0 <= q <= 1/.
 
+[optional_policy]
+
 [h4 Accuracy]
 
 The accuracy of these functions should closely follow that
@@ -268,8 +326,15 @@
 
 These four functions share a common implementation.
 
-First an initial approximation is computed for /a/ or /b/ so that
-I[sub x](a, b) is near the middle of the range [0,1]. This is then
+First an initial approximation is computed for /a/ or /b/:
+where possible this uses a Cornish-Fisher expansion for the
+negative binomial distribution to get within around 1 of the
+result. However, when /a/ or /b/ are very small the Cornish Fisher
+expansion is not usable, in this case the initial approximation
+is chosen so that
+I[sub x](a, b) is near the middle of the range [0,1].
+
+This initial guess is then
 used as a starting value for a generic root finding algorithm. The
 algorithm converges rapidly on the root once it has been
 bracketed, but bracketing the root may take several iterations.

Modified: sandbox/math_toolkit/libs/math/doc/igamma.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/igamma.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/igamma.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:igamma Incomplete Gamma Functions]
 
-[caution __caution ]
-
 [h4 Synopsis]
 
 ``
@@ -13,15 +11,27 @@
    template <class T1, class T2>
    ``__sf_result`` gamma_p(T1 a, T2 z);
    
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` gamma_p(T1 a, T2 z, const ``__Policy``&);
+
    template <class T1, class T2>
    ``__sf_result`` gamma_q(T1 a, T2 z);
    
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` gamma_q(T1 a, T2 z, const ``__Policy``&);
+
    template <class T1, class T2>
    ``__sf_result`` tgamma_lower(T1 a, T2 z);
    
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` tgamma_lower(T1 a, T2 z, const ``__Policy``&);
+
    template <class T1, class T2>
    ``__sf_result`` tgamma(T1 a, T2 z);
    
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` tgamma(T1 a, T2 z, const ``__Policy``&);
+
    }} // namespaces
    
 [h4 Description]
@@ -37,12 +47,17 @@
 All of these functions require /a > 0/ and /z >= 0/, otherwise they return
 the result of __domain_error.
 
+[optional_policy]
+
 The return type of these functions is computed using the __arg_pomotion_rules
 when T1 and T2 are different types, otherwise the return type is simply T1.
 
    template <class T1, class T2>
    ``__sf_result`` gamma_p(T1 a, T2 z);
    
+ template <class T1, class T2, class Policy>
+ ``__sf_result`` gamma_p(T1 a, T2 z, const ``__Policy``&);
+
 Returns the normalised lower incomplete gamma function of a and z:
 
 [$../equations/igamma4.png]
@@ -54,6 +69,9 @@
    template <class T1, class T2>
    ``__sf_result`` gamma_q(T1 a, T2 z);
 
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` gamma_q(T1 a, T2 z, const ``__Policy``&);
+
 Returns the normalised upper incomplete gamma function of a and z:
 
 [$../equations/igamma3.png]
@@ -65,6 +83,9 @@
    template <class T1, class T2>
    ``__sf_result`` tgamma_lower(T1 a, T2 z);
 
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` tgamma_lower(T1 a, T2 z, const ``__Policy``&);
+
 Returns the full (non-normalised) lower incomplete gamma function of a and z:
 
 [$../equations/igamma2.png]
@@ -72,6 +93,9 @@
    template <class T1, class T2>
    ``__sf_result`` tgamma(T1 a, T2 z);
 
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` tgamma(T1 a, T2 z, const ``__Policy``&);
+
 Returns the full (non-normalised) upper incomplete gamma function of a and z:
 
 [$../equations/igamma1.png]

Modified: sandbox/math_toolkit/libs/math/doc/igamma_inv.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/igamma_inv.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/igamma_inv.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:igamma_inv Incomplete Gamma Function Inverses]
 
-[caution __caution ]
-
 [h4 Synopsis]
 
 ``
@@ -13,15 +11,27 @@
    template <class T1, class T2>
    ``__sf_result`` gamma_q_inv(T1 a, T2 q);
    
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` gamma_q_inv(T1 a, T2 q, const ``__Policy``&);
+
    template <class T1, class T2>
    ``__sf_result`` gamma_p_inv(T1 a, T2 p);
    
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` gamma_p_inv(T1 a, T2 p, const ``__Policy``&);
+
    template <class T1, class T2>
    ``__sf_result`` gamma_q_inva(T1 x, T2 q);
    
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` gamma_q_inva(T1 x, T2 q, const ``__Policy``&);
+
    template <class T1, class T2>
    ``__sf_result`` gamma_p_inva(T1 x, T2 p);
    
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` gamma_p_inva(T1 x, T2 p, const ``__Policy``&);
+
    }} // namespaces
    
 [h4 Description]
@@ -34,6 +44,8 @@
 The return type of these functions is computed using the __arg_pomotion_rules
 when T1 and T2 are different types, otherwise the return type is simply T1.
 
+[optional_policy]
+
 [tip When people normally talk about the inverse of the incomplete
 gamma function, they are talking about inverting on parameter /x/.
 These are implemented here as gamma_p_inv and gamma_q_inv, and are by
@@ -48,6 +60,9 @@
    template <class T1, class T2>
    ``__sf_result`` gamma_q_inv(T1 a, T2 q);
 
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` gamma_q_inv(T1 a, T2 q, const ``__Policy``&);
+
 Returns a value x such that: `q = gamma_q(a, x);`
 
 Requires: /a > 0/ and /1 >= p,q >= 0/.
@@ -55,6 +70,9 @@
    template <class T1, class T2>
    ``__sf_result`` gamma_p_inv(T1 a, T2 p);
    
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` gamma_p_inv(T1 a, T2 p, const ``__Policy``&);
+
 Returns a value x such that: `p = gamma_p(a, x);`
 
 Requires: /a > 0/ and /1 >= p,q >= 0/.
@@ -62,6 +80,9 @@
    template <class T1, class T2>
    ``__sf_result`` gamma_q_inva(T1 x, T2 q);
 
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` gamma_q_inva(T1 x, T2 q, const ``__Policy``&);
+
 Returns a value a such that: `q = gamma_q(a, x);`
 
 Requires: /x > 0/ and /1 >= p,q >= 0/.
@@ -69,6 +90,9 @@
    template <class T1, class T2>
    ``__sf_result`` gamma_p_inva(T1 x, T2 p);
    
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` gamma_p_inva(T1 x, T2 p, const ``__Policy``&);
+
 Returns a value a such that: `p = gamma_p(a, x);`
 
 Requires: /x > 0/ and /1 >= p,q >= 0/.

Modified: sandbox/math_toolkit/libs/math/doc/implementation.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/implementation.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/implementation.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,5 +1,7 @@
 [section:implementation Additional Implementation Notes]
 
+TODO: the error handling in this section is out of date and needs a re-write!!
+
 The majority of the implementation notes are included with the documentation
 of each function or distribution. The notes here are of a more general nature,
 and reflect more the general implementation philosophy used.

Modified: sandbox/math_toolkit/libs/math/doc/inv_hyper.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/inv_hyper.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/inv_hyper.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -77,9 +77,12 @@
 #include <boost/math/special_functions/acosh.hpp>
 ``
 
- template<typename T>
+ template<class T>
    ``__sf_result`` acosh(const T x);
 
+ template<class T, class ``__Policy``>
+ ``__sf_result`` acosh(const T x, const ``__Policy``&);
+
 Computes the reciprocal of (the restriction to the range of __form1)
 [link math_toolkit.special.inv_hyper.inv_hyper_over
 the hyperbolic cosine function], at x. Values returned are positive. Generalised
@@ -91,6 +94,8 @@
 The return type of this function is computed using the __arg_pomotion_rules:
 the return type is `double` when T is an integer type, and T otherwise.
 
+[optional_policy]
+
 [endsect]
 
 [section:asinh asinh]
@@ -99,9 +104,12 @@
 #include <boost/math/special_functions/asinh.hpp>
 ``
 
- template<typename T>
+ template<class T>
    ``__sf_result`` asinh(const T x);
 
+ template<class T, class ``__Policy``>
+ ``__sf_result`` asinh(const T x, const ``__Policy``&);
+
 Computes the reciprocal of
 [link math_toolkit.special.inv_hyper.inv_hyper_over
 the hyperbolic sine function].
@@ -111,6 +119,8 @@
 The return type of this function is computed using the __arg_pomotion_rules:
 the return type is `double` when T is an integer type, and T otherwise.
 
+[optional_policy]
+
 [endsect]
 
 [section:atanh atanh]
@@ -119,14 +129,19 @@
 #include <boost/math/special_functions/atanh.hpp>
 ``
 
- template<typename T>
+ template<class T>
    ``__sf_result`` atanh(const T x);
 
+ template<class T, class ``__Policy``>
+ ``__sf_result`` atanh(const T x, const ``__Policy``&);
+
 Computes the reciprocal of
 [link math_toolkit.special.inv_hyper.inv_hyper_over
 the hyperbolic tangent function], at x.
 Taylor series are used at the origin to ensure accuracy.
 
+[optional_policy]
+
 If x is in the range
 __form3
 or in the range

Modified: sandbox/math_toolkit/libs/math/doc/laguerre.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/laguerre.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/laguerre.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:laguerre Laguerre (and Associated) Polynomials]
 
-[caution __caution ]
-
 [h4 Synopsis]
 
 ``
@@ -13,9 +11,15 @@
    template <class T>
    ``__sf_result`` laguerre(unsigned n, T x);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` laguerre(unsigned n, T x, const ``__Policy``&);
+
    template <class T>
    ``__sf_result`` laguerre(unsigned n, unsigned m, T x);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` laguerre(unsigned n, unsigned m, T x, const ``__Policy``&);
+
    template <class T1, class T2, class T3>
    ``__sf_result`` laguerre_next(unsigned n, T1 x, T2 Ln, T3 Lnm1);
    
@@ -31,9 +35,14 @@
 note than when there is a single template argument the result is the same type
 as that argument or `double` if the template argument is an integer type.
 
+[optional_policy]
+
    template <class T>
    ``__sf_result`` laguerre(unsigned n, T x);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` laguerre(unsigned n, T x, const ``__Policy``&);
+
 Returns the value of the Laguerre Polynomial of order /n/ at point /x/:
 
 [$../equations/laguerre_0.png]
@@ -46,6 +55,9 @@
    template <class T>
    ``__sf_result`` laguerre(unsigned n, unsigned m, T x);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` laguerre(unsigned n, unsigned m, T x, const ``__Policy``&);
+
 Returns the Associated Laguerre polynomial of degree /n/
 and order /m/ at point /x/:
 

Modified: sandbox/math_toolkit/libs/math/doc/lanczos.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/lanczos.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/lanczos.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:lanczos The Lanczos Approximation]
 
-[caution __caution ]
-
 [h4 Motivation]
 
 ['Why base gamma and gamma-like functions on the Lanczos approximation?]

Modified: sandbox/math_toolkit/libs/math/doc/legendre.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/legendre.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/legendre.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:legendre Legendre (and Associated) Polynomials]
 
-[caution __caution ]
-
 [h4 Synopsis]
 
 ``
@@ -13,12 +11,21 @@
    template <class T>
    ``__sf_result`` legendre_p(int n, T x);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` legendre_p(int n, T x, const ``__Policy``&);
+
    template <class T>
    ``__sf_result`` legendre_p(int n, int m, T x);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` legendre_p(int n, int m, T x, const ``__Policy``&);
+
    template <class T>
    ``__sf_result`` legendre_q(int n, T x);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` legendre_q(int n, T x, const ``__Policy``&);
+
    template <class T1, class T2, class T3>
    ``__sf_result`` legendre_next(unsigned l, T1 x, T2 Pl, T3 Plm1);
    
@@ -32,11 +39,16 @@
 note than when there is a single template argument the result is the same type
 as that argument or `double` if the template argument is an integer type.
 
+[optional_policy]
+
 [h4 Description]
 
    template <class T>
    ``__sf_result`` legendre_p(int l, T x);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` legendre_p(int l, T x, const ``__Policy``&);
+
 Returns the Legendre Polynomial of the first kind:
 
 [$../equations/legendre_0.png]
@@ -55,6 +67,9 @@
    template <class T>
    ``__sf_result`` legendre_p(int l, int m, T x);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` legendre_p(int l, int m, T x, const ``__Policy``&);
+
 Returns the associated Legendre polynomial of the first kind:
 
 [$../equations/legendre_1.png]
@@ -91,6 +106,9 @@
    template <class T>
    ``__sf_result`` legendre_q(int n, T x);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` legendre_q(int n, T x, const ``__Policy``&);
+
 Returns the value of the Legendre polynomial that is the second solution
 to the Legendre differential equation, for example:
 

Modified: sandbox/math_toolkit/libs/math/doc/lgamma.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/lgamma.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/lgamma.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:lgamma Log Gamma]
 
-[caution __caution ]
-
 [h4 Synopsis]
 
 ``
@@ -13,9 +11,15 @@
    template <class T>
    ``__sf_result`` lgamma(T z);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` lgamma(T z, const ``__Policy``&);
+
    template <class T>
    ``__sf_result`` lgamma(T z, int* sign);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` lgamma(T z, int* sign, const ``__Policy``&);
+
    }} // namespaces
 
 [h4 Description]
@@ -27,6 +31,8 @@
 The second form of the function takes a pointer to an integer,
 which if non-null is set on output to the sign of tgamma(z).
 
+[optional_policy]
+
 [$../graphs/lgamma.png]
 
 There are effectively two versions of this function internally: a fully

Modified: sandbox/math_toolkit/libs/math/doc/math.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/math.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/math.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -195,6 +195,12 @@
 [template floorlr[x][lfloor][x][rfloor]]
 [template ceil[x] '''&#x2308;'''[x]'''&#x2309;''']
 
+[template optional_policy[]
+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 the
+[link math_toolkit.policy policy documentation for more details].]
+
 [section:intro About the Math Toolkit]
 
 This library is divided into three interconnected parts:

Modified: sandbox/math_toolkit/libs/math/doc/minima.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/minima.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/minima.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:minima Locating Function Minima]
 
-[caution __caution ]
-
 [h4 synopsis]
 
 ``

Modified: sandbox/math_toolkit/libs/math/doc/polynomial.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/polynomial.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/polynomial.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:polynomials Polynomials]
 
-[caution __caution ]
-
 [h4 Synopsis]
 
 ``

Modified: sandbox/math_toolkit/libs/math/doc/powers.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/powers.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/powers.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:powers Logs, Powers, Roots and Exponentials]
 
-[caution __caution ]
-
 [section:log1p log1p]
 
 ``
@@ -13,6 +11,9 @@
    template <class T>
    ``__sf_result`` log1p(T x);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` log1p(T x, const ``__Policy``&);
+
    }} // namespaces
    
 Returns the natural logarithm of `x+1`.
@@ -20,6 +21,8 @@
 The return type of this function is computed using the __arg_pomotion_rules:
 the return is `double` when /x/ is an integer type and T otherwise.
 
+[optional_policy]
+
 There are many situations where it is desirable to compute `log(x+1)`.
 However, for small `x` then `x+1` suffers from catastrophic cancellation errors
 so that `x+1 == 1` and `log(x+1) == 0`, when in fact for very small x, the
@@ -64,6 +67,9 @@
    template <class T>
    ``__sf_result`` expm1(T x);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` expm1(T x, const ``__Policy``&);
+
    }} // namespaces
    
 Returns e[super x] - 1.
@@ -71,6 +77,8 @@
 The return type of this function is computed using the __arg_pomotion_rules:
 the return is `double` when /x/ is an integer type and T otherwise.
 
+[optional_policy]
+
 For small x, then __ex is very close to 1, as a result calculating __exm1 results
 in catastrophic cancellation errors when x is small. `expm1` calculates __exm1 using
 rational approximations (for up to 128-bit long doubles), otherwise via
@@ -103,6 +111,9 @@
    template <class T>
    ``__sf_result`` cbrt(T x);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` cbrt(T x, const ``__Policy``&);
+
    }} // namespaces
    
 Returns the cubed root of x: x[super 1/3].
@@ -110,6 +121,8 @@
 The return type of this function is computed using the __arg_pomotion_rules:
 the return is `double` when /x/ is an integer type and T otherwise.
 
+[optional_policy]
+
 Implemented using Halley iteration.
    
 [h4 Accuracy]
@@ -135,6 +148,9 @@
    template <class T>
    ``__sf_result`` sqrt1pm1(T x);
    
+ template <class T, class ``__Policy``>
+ ``__sf_result`` sqrt1pm1(T x, const ``__Policy``&);
+
    }} // namespaces
    
 Returns `sqrt(1+x) - 1`.
@@ -142,6 +158,8 @@
 The return type of this function is computed using the __arg_pomotion_rules:
 the return is `double` when /x/ is an integer type and T otherwise.
 
+[optional_policy]
+
 This function is useful when you need the difference between sqrt(x) and 1, when
 x is itself close to 1.
 
@@ -170,6 +188,9 @@
    template <class T1, class T2>
    ``__sf_result`` powm1(T1 x, T2 y);
    
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` powm1(T1 x, T2 y, const ``__Policy``&);
+
    }} // namespaces
    
 Returns x[super y ] - 1.
@@ -177,6 +198,8 @@
 The return type of this function is computed using the __arg_pomotion_rules
 when T1 and T2 are dufferent types.
 
+[optional_policy]
+
 There are two domains where this is useful: when y is very small, or when
 x is close to 1.
 
@@ -198,12 +221,17 @@
    template <class T1, class T2>
    ``__sf_result`` hypot(T1 x, T2 y);
    
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` hypot(T1 x, T2 y, const ``__Policy``&);
+
 __effects computes [$../equations/hypot.png]
 in such a way as to avoid undue underflow and overflow.
 
 The return type of this function is computed using the __arg_pomotion_rules
 when T1 and T2 are of different types.
 
+[optional_policy]
+
 When calculating [$../equations/hypot.png] it's quite easy for the intermediate terms to either
 overflow or underflow, even though the result is in fact perfectly representable.
 

Modified: sandbox/math_toolkit/libs/math/doc/rational.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/rational.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/rational.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:rational Polynomial and Rational Function Evaluation]
 
-[caution __caution ]
-
 [h4 synopsis]
 
 ``
@@ -157,6 +155,9 @@
 order as polynomials in ['1\/v]: this avoids unnecessary numerical overflow when the
 coefficients are large.
 
+TODO: mention second order Horner's method and other evaluation options,
+plus config, and performance test app.
+
 [endsect][/section:rational Polynomial and Rational Function Evaluation]
 
 [/

Modified: sandbox/math_toolkit/libs/math/doc/relative_error.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/relative_error.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/relative_error.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:error_test Relative Error and Testing]
 
-[caution __caution ]
-
 [h4 Synopsis]
 
 ``

Modified: sandbox/math_toolkit/libs/math/doc/roadmap.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/roadmap.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/roadmap.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -42,6 +42,24 @@
 and brought them into line with the rest of the code.
 * Added C# "Distribution Explorer" demo application.
 
+[h4 Milestone 5: Post Review First Official Release (in Progress)]
+
+['The documentation for this first release is still being worked upon.]
+
+* Added Policy based framework that allows fine grained control
+over function behaviour.
+* [*Breaking change:] Changed default behaviour for domain, pole and overflow errors
+to throw an exception (based on review feedback), this
+behaviour can be customised using __Policy's.
+* [*Breaking change:] Changed exception thrown when an internal evaluation error
+occurs to boost::math::evaluation_error.
+* [*Breaking change:] Changed discrete quantiles to return an integer result:
+this is anything up to 20 times faster than finding the true root, this
+behaviour can be customised using __Policy's.
+* Polynomial/rational function evaluation is now customisable and hopefully
+faster than before.
+* Added performance test program.
+
 [h4 Future]
 
 * Some TR1 Special functions are still needed.

Modified: sandbox/math_toolkit/libs/math/doc/roots.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/roots.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/roots.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:roots Root Finding With Derivatives]
 
-[caution __caution ]
-
 [h4 Synopsis]
 
 ``

Modified: sandbox/math_toolkit/libs/math/doc/roots_without_derivatives.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/roots_without_derivatives.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/roots_without_derivatives.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:roots2 Root Finding Without Derivatives]
 
-[caution __caution ]
-
 [h4 Synopsis]
 
 ``
@@ -27,6 +25,16 @@
          T max,
          Tol tol);
 
+ template <class F, class T, class Tol, class ``__Policy``>
+ std::pair<T, T>
+ bisect(
+ F f,
+ T min,
+ T max,
+ Tol tol,
+ boost::uintmax_t& max_iter,
+ const ``__Policy``&);
+
    template <class F, class T, class Tol>
    std::pair<T, T>
       bracket_and_solve_root(
@@ -37,6 +45,17 @@
          Tol tol,
          boost::uintmax_t& max_iter);
 
+ template <class F, class T, class Tol, class ``__Policy``>
+ std::pair<T, T>
+ bracket_and_solve_root(
+ F f,
+ const T& guess,
+ const T& factor,
+ bool rising,
+ Tol tol,
+ boost::uintmax_t& max_iter,
+ const ``__Policy``&);
+
    template <class F, class T, class Tol>
    std::pair<T, T>
       toms748_solve(
@@ -46,6 +65,16 @@
          Tol tol,
          boost::uintmax_t& max_iter);
 
+ template <class F, class T, class Tol, class ``__Policy``>
+ std::pair<T, T>
+ toms748_solve(
+ F f,
+ const T& a,
+ const T& b,
+ Tol tol,
+ boost::uintmax_t& max_iter,
+ const ``__Policy``&);
+
    template <class F, class T, class Tol>
    std::pair<T, T>
       toms748_solve(
@@ -57,12 +86,25 @@
          Tol tol,
          boost::uintmax_t& max_iter);
 
+ template <class F, class T, class Tol, class ``__Policy``>
+ std::pair<T, T>
+ toms748_solve(
+ F f,
+ const T& a,
+ const T& b,
+ const T& fa,
+ const T& fb,
+ Tol tol,
+ boost::uintmax_t& max_iter,
+ const ``__Policy``&);
+
    // Termination conditions:
    template <class T>
    struct eps_tolerance;
 
    struct equal_floor;
    struct equal_ceil;
+ struct equal_nearest_integer;
    
    }}} // namespaces
    
@@ -111,7 +153,17 @@
          T max,
          Tol tol);
 
-These two functions locate the root using bisection, function arguments are:
+ template <class F, class T, class Tol, class ``__Policy``>
+ std::pair<T, T>
+ bisect(
+ F f,
+ T min,
+ T max,
+ Tol tol,
+ boost::uintmax_t& max_iter,
+ const ``__Policy``&);
+
+These functions locate the root using bisection, function arguments are:
 
 [variablelist
 [[f] [A unary functor which is the function whose root is to be found.]]
@@ -124,6 +176,8 @@
 [[max_iter][The maximum number of invocations of /f(x)/ to make while searching for the root.]]
 ]
 
+[optional_policy]
+
 Returns: a pair of values /r/ that bracket the root so that:
 
    f(r.first) * f(r.second) <= 0
@@ -153,6 +207,17 @@
          Tol tol,
          boost::uintmax_t& max_iter);
 
+ template <class F, class T, class Tol, class ``__Policy``>
+ std::pair<T, T>
+ bracket_and_solve_root(
+ F f,
+ const T& guess,
+ const T& factor,
+ bool rising,
+ Tol tol,
+ boost::uintmax_t& max_iter,
+ const ``__Policy``&);
+
 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
@@ -177,6 +242,8 @@
 [[max_iter] [The maximum number of function invocations to perform in the search
             for the root.]]
 ]
+
+[optional_policy]
          
 Returns: a pair of values /r/ that bracket the root so that:
 
@@ -206,6 +273,16 @@
          Tol tol,
          boost::uintmax_t& max_iter);
 
+ template <class F, class T, class Tol, class ``__Policy``>
+ std::pair<T, T>
+ toms748_solve(
+ F f,
+ const T& a,
+ const T& b,
+ Tol tol,
+ boost::uintmax_t& max_iter,
+ const ``__Policy``&);
+
    template <class F, class T, class Tol>
    std::pair<T, T>
       toms748_solve(
@@ -217,6 +294,18 @@
          Tol tol,
          boost::uintmax_t& max_iter);
          
+ template <class F, class T, class Tol, class ``__Policy``>
+ std::pair<T, T>
+ toms748_solve(
+ F f,
+ const T& a,
+ const T& b,
+ const T& fa,
+ const T& fb,
+ Tol tol,
+ boost::uintmax_t& max_iter,
+ const ``__Policy``&);
+
 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
@@ -240,6 +329,8 @@
             invocations used.]]
 ]
 
+[optional_policy]
+
 Returns: a pair of values /r/ that bracket the root so that:
 
    f(r.first) * f(r.second) <= 0
@@ -294,6 +385,16 @@
 that is the /ceil/ of the true root. It will terminate as soon as both ends
 of the interval have the same /ceil/.
 
+ struct equal_nearest_integer
+ {
+ equal_nearest_integer();
+ template <class T> 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 /closest/ to the true root. It will terminate as soon as both ends
+of the interval round to the same nearest integer.
+
 [h4 Implementation]
 
 The implementation of the bisection algorithm is extremely straightforward

Modified: sandbox/math_toolkit/libs/math/doc/series.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/series.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/series.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:series_evaluation Series Evaluation]
 
-[caution __caution ]
-
 [h4 Synopsis]
 
 ``

Modified: sandbox/math_toolkit/libs/math/doc/sinc.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/sinc.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/sinc.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -39,12 +39,18 @@
 #include <boost/math/special_functions/sinc.hpp>
 ``
 
- template<typename T>
+ template<class T>
    ``__sf_result`` sinc_pi(const T x);
 
- template<typename T, template<typename> class U>
+ template<class T, class ``__Policy``>
+ ``__sf_result`` sinc_pi(const T x, const ``__Policy``&);
+
+ template<class T, template<typename> class U>
    U<T> sinc_pi(const U<T> x);
 
+ template<class T, template<typename> class U, class ``__Policy``>
+ U<T> sinc_pi(const U<T> x, const ``__Policy``&);
+
 Computes
 [link math_toolkit.special.sinc.sinc_overview
 the Sinus Cardinal] of x:
@@ -55,6 +61,7 @@
 quaternions, octonions etc. Taylor series are used at the origin
 to ensure accuracy.
 
+[optional_policy]
 
 [endsect]
 
@@ -64,12 +71,18 @@
 #include <boost/math/special_functions/sinhc.hpp>
 ``
 
- template<typename T>
+ template<class T>
    ``__sf_result`` sinhc_pi(const T x);
 
+ template<class T, class ``__Policy``>
+ ``__sf_result`` sinhc_pi(const T x, const ``__Policy``&);
+
    template<typename T, template<typename> class U>
    U<T> sinhc_pi(const U<T> x);
 
+ template<class T, template<typename> class U, class ``__Policy``>
+ U<T> sinhc_pi(const U<T> x, const ``__Policy``&);
+
 Computes http://mathworld.wolfram.com/SinhcFunction.html
 [link math_toolkit.special.sinc.sinc_overview
 the Hyperbolic Sinus Cardinal] of x:
@@ -83,6 +96,8 @@
 The return type of the first form is computed using the __arg_pomotion_rules
 when T is an integer type.
 
+[optional_policy]
+
 [endsect]
 
 [endsect]

Modified: sandbox/math_toolkit/libs/math/doc/spherical_harmonic.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/spherical_harmonic.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/spherical_harmonic.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:sph_harm Spherical Harmonics]
 
-[caution __caution ]
-
 [h4 Synopsis]
 
 ``
@@ -13,12 +11,21 @@
    template <class T1, class T2>
    std::complex<``__sf_result``> spherical_harmonic(unsigned n, int m, T1 theta, T2 phi);
 
+ template <class T1, class T2, class ``__Policy``>
+ std::complex<``__sf_result``> spherical_harmonic(unsigned n, int m, T1 theta, T2 phi, const ``__Policy``&);
+
    template <class T1, class T2>
    ``__sf_result`` spherical_harmonic_r(unsigned n, int m, T1 theta, T2 phi);
 
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` spherical_harmonic_r(unsigned n, int m, T1 theta, T2 phi, const ``__Policy``&);
+
    template <class T1, class T2>
    ``__sf_result`` spherical_harmonic_i(unsigned n, int m, T1 theta, T2 phi);
       
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` spherical_harmonic_i(unsigned n, int m, T1 theta, T2 phi, const ``__Policy``&);
+
    }} // namespaces
 
 [h4 Description]
@@ -26,9 +33,14 @@
 The return type of these functions is computed using the __arg_pomotion_rules
 when T1 and T2 are different types.
 
+[optional_policy]
+
    template <class T1, class T2>
    std::complex<``__sf_result``> spherical_harmonic(unsigned n, int m, T1 theta, T2 phi);
    
+ template <class T1, class T2, class ``__Policy``>
+ std::complex<``__sf_result``> spherical_harmonic(unsigned n, int m, T1 theta, T2 phi, const ``__Policy``&);
+
 Returns the value of the Spherical Harmonic Y[sub n][super m](theta, phi):
 
 [$../equations/spherical_0.png]
@@ -67,6 +79,9 @@
    template <class T1, class T2>
    ``__sf_result`` spherical_harmonic_r(unsigned n, int m, T1 theta, T2 phi);
    
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` spherical_harmonic_r(unsigned n, int m, T1 theta, T2 phi, const ``__Policy``&);
+
 Returns the real part of Y[sub n][super m](theta, phi):
 
 [$../equations/spherical_1.png]
@@ -74,6 +89,9 @@
    template <class T1, class T2>
    ``__sf_result`` spherical_harmonic_i(unsigned n, int m, T1 theta, T2 phi);
    
+ template <class T1, class T2, class ``__Policy``>
+ ``__sf_result`` spherical_harmonic_i(unsigned n, int m, T1 theta, T2 phi, const ``__Policy``&);
+
 Returns the imaginary part of Y[sub n][super m](theta, phi):
 
 [$../equations/spherical_2.png]

Modified: sandbox/math_toolkit/libs/math/doc/structure.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/structure.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/structure.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -64,8 +64,11 @@
 See [link math_toolkit.dist.stat_tut Statistical Distributions Tutorial] and
 [link math_toolkit.dist.dist_ref Statistical Distributions Reference]
 and [link math_toolkit.special Special Functions]
+
 [h4 Configuration Policy Macros - Controlling Handling Arguments, Undefined Functions, Throwing of Exceptions, Accuracy, ]
 
+TODO: this is out of date, rewrite me!!
+
 * *BOOST_MATH_DOMAIN_ERROR_POLICY*
 
 See [link math_toolkit.special.error_handling Error handling]

Modified: sandbox/math_toolkit/libs/math/doc/tgamma.qbk
==============================================================================
--- sandbox/math_toolkit/libs/math/doc/tgamma.qbk (original)
+++ sandbox/math_toolkit/libs/math/doc/tgamma.qbk 2007-07-28 11:13:49 EDT (Sat, 28 Jul 2007)
@@ -1,7 +1,5 @@
 [section:tgamma Gamma]
 
-[caution __caution ]
-
 [h4 Synopsis]
 
 ``
@@ -13,9 +11,15 @@
   template <class T>
   ``__sf_result`` tgamma(T z);
   
+ template <class T, class ``__Policy``>
+ ``__sf_result`` tgamma(T z, const ``__Policy``&);
+
   template <class T>
   ``__sf_result`` tgamma1pm1(T dz);
   
+ template <class T, class ``__Policy``>
+ ``__sf_result`` tgamma1pm1(T dz, const ``__Policy``&);
+
   }} // namespaces
   
 [h4 Description]
@@ -23,12 +27,17 @@
   template <class T>
   ``__sf_result`` tgamma(T z);
   
+ template <class T, class ``__Policy``>
+ ``__sf_result`` tgamma(T z, const ``__Policy``&);
+
 Returns the "true gamma" (hence name tgamma) of value z:
 
 [$../equations/gamm1.png]
 
 [$../graphs/gamma.png]
 
+[optional_policy]
+
 There are effectively two versions of the [@http://en.wikipedia.org/wiki/Gamma_function tgamma]
 function internally: a fully
 generic version that is slow, but reasonably accurate, and a much more
@@ -44,6 +53,9 @@
   template <class T>
   ``__sf_result`` tgamma1pm1(T dz);
   
+ template <class T, class ``__Policy``>
+ ``__sf_result`` tgamma1pm1(T dz, const ``__Policy``&);
+
 Returns `tgamma(dz + 1) - 1`. Internally the implementation does not make
 use of the addition and subtraction implied by the definition, leading to
 accurate results even for very small `dz`. However, the implementation is
@@ -53,6 +65,8 @@
 The return type of this function is computed using the __arg_pomotion_rules:
 the result is `double` when T is an integer type, and T otherwise.
 
+[optional_policy]
+
 [h4 Accuracy]
 
 The following table shows the peak errors (in units of epsilon)


Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk