Merge and run CI

.gitignore

@@ -5,3 +5,4 @@ reporting/*/html
 inspect.html
 /example/float128_examples.cpp
 test/cuda
+*.DS_Store

.travis.yml

@@ -546,83 +546,24 @@ matrix:
 
     - os: osx
       env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE=special_fun
-      osx_image: xcode8.1
+      osx_image: xcode9.3
 
     - os: osx
       env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE=distribution_tests
-      osx_image: xcode8.1
+      osx_image: xcode9.3
 
     - os: osx
       env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE=misc
-      osx_image: xcode8.1
+      osx_image: xcode9.3
 
     - os: osx
       env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE=float128_tests
-      osx_image: xcode8.1
+      osx_image: xcode9.3
 
     - os: osx
       env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE=../example//examples
-      osx_image: xcode8.1
+      osx_image: xcode9.3
 
-    - os: osx
-      env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE=special_fun
-      osx_image: xcode8.3
-
-    - os: osx
-      env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE=distribution_tests
-      osx_image: xcode8.3
-
-    - os: osx
-      env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE=misc
-      osx_image: xcode8.3
-
-    - os: osx
-      env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE=float128_tests
-      osx_image: xcode8.3
-
-    - os: osx
-      env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE=../example//examples
-      osx_image: xcode8.3
-
-    - os: osx
-      env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE=special_fun
-      osx_image: xcode7.3
-
-    - os: osx
-      env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE=distribution_tests
-      osx_image: xcode7.3
-
-    - os: osx
-      env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE=misc
-      osx_image: xcode7.3
-
-    - os: osx
-      env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE=float128_tests
-      osx_image: xcode7.3
-
-    - os: osx
-      env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE=../example//examples
-      osx_image: xcode7.3
-
-    - os: osx
-      env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE=special_fun
-      osx_image: xcode6.4
-
-    - os: osx
-      env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE=distribution_tests
-      osx_image: xcode6.4
-
-    - os: osx
-      env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE=misc
-      osx_image: xcode6.4
-
-    - os: osx
-      env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE=float128_tests
-      osx_image: xcode6.4
-
-    - os: osx
-      env: TOOLSET=clang COMPILER=clang++ CXXSTD=c++14 TEST_SUITE=../example//examples
-      osx_image: xcode6.4
 
 install:
   - cd ..
@@ -632,6 +573,7 @@ install:
   - git submodule update --init libs/config
   - git submodule update --init libs/format
   - git submodule update --init libs/numeric
+  - git submodule update --init libs/type_traits
   - git submodule update --init tools/boostdep
   - cp -r $TRAVIS_BUILD_DIR/* libs/math
   - python tools/boostdep/depinst/depinst.py math
@@ -642,7 +584,7 @@ script:
   - |-
     echo "using $TOOLSET : : $COMPILER : <cxxflags>-std=$CXXSTD ;" > ~/user-config.jam
   - (cd libs/config/test && ../../../b2 config_info_travis_install toolset=$TOOLSET && ./config_info_travis)
-  - (cd libs/math/test && ../../../b2 -j3 toolset=$TOOLSET $TEST_SUITE define=CI_SUPPRESS_KNOWN_ISSUES define=SLOW_COMPILER)
+  - (cd libs/math/test && travis_wait 60 ../../../b2 -j3 -d+0 -q toolset=$TOOLSET $TEST_SUITE define=CI_SUPPRESS_KNOWN_ISSUES define=SLOW_COMPILER)
 
 notifications:
   email:

README.md

@@ -0,0 +1,91 @@
+Boost Math Library
+==================
+
+This library is divided into several interconnected parts:
+
+### Floating Point Utilities
+
+Utility functions for dealing with floating point arithmetic, includes functions for floating point classification (fpclassify, isnan, isinf etc), sign manipulation, rounding, comparison, and computing the distance between floating point numbers.
+
+### Specific Width Floating Point Types
+
+A set of typedefs similar to those provided by <cstdint> but for floating point types.
+
+### Mathematical Constants
+
+A wide range of constants ranging from various multiples of π, fractions, through to euler's constant etc.
+
+These are of course usable from template code, or as non-templates with a simplified interface if that is more appropriate.
+
+### Statistical Distributions
+
+Provides a reasonably comprehensive set of statistical distributions, upon which higher level statistical tests can be built.
+
+The initial focus is on the central univariate distributions. Both continuous (like normal & Fisher) and discrete (like binomial & Poisson) distributions are provided.
+
+A comprehensive tutorial is provided, along with a series of worked examples illustrating how the library is used to conduct statistical tests.
+
+### Special Functions
+
+Provides a small number of high quality special functions, initially these were concentrated on functions used in statistical applications along with those in the Technical Report on C++ Library Extensions.
+
+The function families currently implemented are the gamma, beta & erf functions along with the incomplete gamma and beta functions (four variants of each) and all the possible inverses of these, plus digamma, various factorial functions, Bessel functions, elliptic integrals, sinus cardinals (along with their hyperbolic variants), inverse hyperbolic functions, Legrendre/Laguerre/Hermite/Chebyshev polynomials and various special power and logarithmic functions.
+
+All the implementations are fully generic and support the use of arbitrary "real-number" types, including Boost.Multiprecision, although they are optimised for use with types with known-about significand (or mantissa) sizes: typically float, double or long double.
+
+These functions also provide the basis of support for the TR1 special functions.
+
+### Root Finding and Function Minimisation
+
+A comprehensive set of root finding algorithms over the real-line, both with derivatives and derivative free.
+
+Also function minimisation via Brent's Method.
+
+### Polynomials and Rational Functions
+
+Tools for manipulating polynomials and for efficient evaluation of rationals or polynomials.
+
+### Interpolation
+
+Function interpolation via barycentric rational interpolation, compactly supported cubic B-splines, and the Chebyshev transform.
+
+### Numerical Integration and Differentiation
+
+A reasonably comprehensive set of routines for integration (trapezoidal, Gauss-Legendre, Gauss-Kronrod and double-exponential) and differentiation.
+
+The integration routines are usable for functions returning complex results - and hence can be used for computation of  contour integrals.
+
+### Quaternions and Octonions
+
+Quaternion and Octonians as class templates similar to std::complex. 
+
+The full documentation is available on [boost.org](http://www.boost.org/doc/libs/release/libs/math).
+
+## Support, bugs and feature requests ##
+
+Bugs and feature requests can be reported through the [Gitub issue tracker](https://github.com/boostorg/math/issues)
+(see [open issues](https://github.com/boostorg/math/issues) and
+[closed issues](https://github.com/boostorg/math/issues?utf8=%E2%9C%93&q=is%3Aissue+is%3Aclosed)).
+
+You can submit your changes through a [pull request](https://github.com/boostorg/math/pulls).
+
+There is no mailing-list specific to Boost Math, although you can use the general-purpose Boost [mailing-list](http://lists.boost.org/mailman/listinfo.cgi/boost-users) using the tag [math].
+
+
+## Development ##
+
+Clone the whole boost project, which includes the individual Boost projects as submodules ([see boost+git doc](https://github.com/boostorg/boost/wiki/Getting-Started)): 
+
+    git clone https://github.com/boostorg/boost
+    cd boost
+    git submodule update --init
+
+The Boost Math Library is located in `libs/math/`. 
+
+### Running tests ###
+First, make sure you are in `libs/math/test`. 
+You can either run all the tests listed in `Jamfile.v2` or run a single test:
+
+    ../../../b2                        <- run all tests
+    ../../../b2 static_assert_test     <- single test
+

appveyor.yml

@@ -23,16 +23,25 @@ environment:
     - ARGS: --toolset=msvc-14.0 address-model=32
     - ARGS: --toolset=msvc-12.0 address-model=64
     - ARGS: --toolset=msvc-14.0 address-model=64
-    #- ARGS: --toolset=msvc-14.0 address-model=64 cxxflags=-std:c++latest
+      TEST_SUITE: special_fun distribution_tests
+    - ARGS: --toolset=msvc-14.0 address-model=64
+      TEST_SUITE: misc ../example//examples
 
     - APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2017
       ARGS: --toolset=msvc-14.1 address-model=64
+      TEST_SUITE: special_fun distribution_tests
+    
+    - APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2017
+      ARGS: --toolset=msvc-14.1 address-model=64
+      TEST_SUITE: misc ../example//examples
     
     - APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2017
       ARGS: --toolset=msvc-14.1 address-model=32
+      TEST_SUITE: special_fun distribution_tests
     
-    #- APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2017
-    #  ARGS: --toolset=msvc-14.1 address-model=64 cxxflags=-std:c++latest
+    - APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2017
+      ARGS: --toolset=msvc-14.1 address-model=32
+      TEST_SUITE: misc ../example//examples
     
     - ARGS: --toolset=gcc address-model=64
       TEST_SUITE: float128_tests

build/Jamfile.v2

@@ -120,8 +120,7 @@ lib boost_math_c99l : ../src/tr1/$(C99_SOURCES)l.cpp pch
          [ check-target-builds ../config//has_long_double_support "long double support" : : <build>no ]
    ;
 
-boost-install boost_math_c99 boost_math_c99f boost_math_c99l boost_math_tr1 boost_math_tr1f boost_math_tr1l ;
-
-
-
+# Fake 'library' target for dependency resolution purposes
+alias boost_math ;
 
+boost-install boost_math_c99 boost_math_c99f boost_math_c99l boost_math_tr1 boost_math_tr1f boost_math_tr1l boost_math ;

config/has_fftw3.cpp

@@ -7,10 +7,13 @@
 
 int main()
 {
+   fftwq_plan plan;  // early versions don't have this it seems.
+
    fftw_cleanup();
    fftwf_cleanup();
    fftwl_cleanup();
 
+
    return 0;
 }
 

doc/background/remez.qbk

@@ -8,9 +8,9 @@ should consult your favorite textbook.
 
 Imagine that you want to approximate some function f(x) by way of a rational
 function R(x), where R(x) may be either a polynomial P(x) or a ratio of two
-polynomials P(x)/Q(x) (a rational function).  Initially we'll concentrate on the 
-polynomial case, as it's by far the easier to deal with, later we'll extend 
-to the full rational function case.  
+polynomials P(x)/Q(x) (a rational function).  Initially we'll concentrate on the
+polynomial case, as it's by far the easier to deal with, later we'll extend
+to the full rational function case.
 
 We want to find the "best" rational approximation, where
 "best" is defined to be the approximation that has the least deviation
@@ -61,11 +61,11 @@ the range \[-1, 1\].
 
 Before we can begin the Remez method, we must obtain an initial value
 for the location of the extrema of the error function.  We could "guess"
-these, but a much closer first approximation can be obtained by first  
+these, but a much closer first approximation can be obtained by first
 constructing an interpolated polynomial approximation to f(x).
 
 In order to obtain the N+1 coefficients of the interpolated polynomial
-we need N+1 points (x[sub 0]...x[sub N]): with our interpolated form 
+we need N+1 points (x[sub 0]...x[sub N]): with our interpolated form
 passing through each of those points
 that yields N+1 simultaneous equations:
 
@@ -73,12 +73,12 @@ f(x[sub i]) = P(x[sub i]) = c[sub 0] + c[sub 1]x[sub i] ... + c[sub N]x[sub i][s
 
 Which can be solved for the coefficients c[sub 0]...c[sub N] in P(x).
 
-Obviously this is not a minimax solution, indeed our only guarantee is that f(x) and 
+Obviously this is not a minimax solution, indeed our only guarantee is that f(x) and
 P(x) touch at N+1 locations, away from those points the error may be arbitrarily
 large.  However, we would clearly like this initial approximation to be as close to
 f(x) as possible, and it turns out that using the zeros of an orthogonal polynomial
-as the initial interpolation points is a good choice.  In our example we'll use the 
-zeros of a Chebyshev polynomial as these are particularly easy to calculate, 
+as the initial interpolation points is a good choice.  In our example we'll use the
+zeros of a Chebyshev polynomial as these are particularly easy to calculate,
 interpolating for a polynomial of degree 4, and measuring /relative error/
 we get the following error function:
 
@@ -86,7 +86,7 @@ we get the following error function:
 
 Which has a peak relative error of 1.2x10[super -3].
 
-While this is a pretty good approximation already, judging by the 
+While this is a pretty good approximation already, judging by the
 shape of the error function we can clearly do better.  Before starting
 on the Remez method propper, we have one more step to perform: locate
 all the extrema of the error function, and store
@@ -101,7 +101,7 @@ achieve, and typically is very close to the minimax solution.
 
 However, if we want to optimise for /relative error/, or if the approximation
 is a rational function, then the initial Chebyshev solution can be quite far
-from the ideal minimax solution.  
+from the ideal minimax solution.
 
 A more technical discussion of the theory involved can be found in this
 [@http://math.fullerton.edu/mathews/n2003/ChebyshevPolyMod.html online course].]
@@ -118,7 +118,7 @@ To obtain the error term E, and the coefficients of the polynomial P(x).
 
 This gives us a new approximation to f(x) that has the same error /E/ at
 each of the control points, and whose error function ['alternates in sign]
-at the control points.  This is still not necessarily the minimax 
+at the control points.  This is still not necessarily the minimax
 solution though: since the control points may not be at the extrema of the error
 function.  After this first step here's what our approximation's error
 function looks like:
@@ -131,10 +131,10 @@ dropped to 5.6x10[super -4].
 
 [h4 Remez Step 2]
 
-The second step is to locate the extrema of the new approximation, which we do 
+The second step is to locate the extrema of the new approximation, which we do
 in two stages:  first, since the error function changes sign at each
 control point, we must have N+1 roots of the error function located between
-each pair of N+2 control points.  Once these roots are found by standard root finding 
+each pair of N+2 control points.  Once these roots are found by standard root finding
 techniques, we know that N extrema are bracketed between each pair of
 roots, plus two more between the endpoints of the range and the first and last roots.
 The N+2 extrema can then be found using standard function minimisation techniques.
@@ -195,10 +195,10 @@ polynomial alternatives.  For example, if we takes our previous example
 of an approximation to e[super x], we obtained 5x10[super -4] accuracy
 with an order 4 polynomial.  If we move two of the unknowns into the denominator
 to give a pair of order 2 polynomials, and re-minimise, then the peak relative error drops
-to 8.7x10[super -5].  That's a 5 fold increase in accuracy, for the same number 
+to 8.7x10[super -5].  That's a 5 fold increase in accuracy, for the same number
 of terms overall.
 
-[h4 Practical Considerations]
+[h4:remez_practical Practical Considerations]
 
 Most treatises on approximation theory stop at this point.  However, from
 a practical point of view, most of the work involves finding the right
@@ -211,12 +211,12 @@ f(x) = R(x)
 
 But this will converge to a useful approximation only if f(x) is smooth.  In
 addition round-off errors when evaluating the rational form mean that this
-will never get closer than within a few epsilon of machine precision.  
+will never get closer than within a few epsilon of machine precision.
 Therefore this form of direct approximation is often reserved for situations
 where we want efficiency, rather than accuracy.
 
 The first step in improving the situation is generally to split f(x) into
-a dominant part that we can compute accurately by another method, and a 
+a dominant part that we can compute accurately by another method, and a
 slowly changing remainder which can be approximated by a rational approximation.
 We might be tempted to write:
 
@@ -235,7 +235,7 @@ to the constant /c/, that error will effectively get wiped out when R(x) is adde
 
 The difficult part is obviously finding the right g(x) to extract from your
 function: often the asymptotic behaviour of the function will give a clue, so
-for example the function __erfc becomes proportional to 
+for example the function __erfc becomes proportional to
 e[super -x[super 2]]\/x as x becomes large.  Therefore using:
 
 erfc(z) = (C + R(x)) e[super -x[super 2]]/x
@@ -262,15 +262,15 @@ the roots, the approximation is nonetheless useless.
 Assuming that the approximation does not have any fatal errors, and that the only
 issue is converging adequately on the minimax solution, the aim is to
 get as close as possible to the minimax solution before beginning the Remez method.
-Using the zeros of a Chebyshev polynomial for the initial interpolation is a 
+Using the zeros of a Chebyshev polynomial for the initial interpolation is a
 good start, but may not be ideal when dealing with relative errors and\/or
 rational (rather than polynomial) approximations.  One approach is to skew
 the initial interpolation points to one end: for example if we raise the
-roots of the Chebyshev polynomial to a positive power greater than 1 
-then the roots will be skewed towards the middle of the \[-1,1\] interval, 
+roots of the Chebyshev polynomial to a positive power greater than 1
+then the roots will be skewed towards the middle of the \[-1,1\] interval,
 while a positive power less than one
 will skew them towards either end.  More usefully, if we initially rescale the
-points over \[0,1\] and then raise to a positive power, we can skew them to the left 
+points over \[0,1\] and then raise to a positive power, we can skew them to the left
 or right.  Returning to our example of e[super x][space] over \[-1,1\], the initial
 interpolated form was some way from the minimax solution:
 
@@ -287,7 +287,7 @@ our desired minimax solution (5x10[super -4]).
 
 [h4 Remez Method Checklist]
 
-The following lists some of the things to check if the Remez method goes wrong, 
+The following lists some of the things to check if the Remez method goes wrong,
 it is by no means an exhaustive list, but is provided in the hopes that it will
 prove useful.
 
@@ -312,7 +312,7 @@ location, try starting from a different location.
 * If a rational function won't converge, one can minimise a polynomial
 (which presents no problems), then rotate one term from the numerator to
 the denominator and minimise again.  In theory one can continue moving
-terms one at a time from numerator to denominator, and then re-minimising, 
+terms one at a time from numerator to denominator, and then re-minimising,
 retaining the last set of control points at each stage.
 * Try using a smaller interval.  It may also be possible to optimise over
 one (small) interval, rescale the control points over a larger interval,
@@ -325,41 +325,41 @@ over, say \[0, b\], for some smallish value /b/.
 The original references for the Remez Method and it's extension
 to rational functions are unfortunately in Russian:
 
-Remez, E.Ya., ['Fundamentals of numerical methods for Chebyshev approximations], 
+Remez, E.Ya., ['Fundamentals of numerical methods for Chebyshev approximations],
 "Naukova Dumka", Kiev, 1969.
 
-Remez, E.Ya., Gavrilyuk, V.T., ['Computer development of certain approaches 
-to the approximate construction of solutions of Chebyshev problems 
+Remez, E.Ya., Gavrilyuk, V.T., ['Computer development of certain approaches
+to the approximate construction of solutions of Chebyshev problems
 nonlinearly depending on parameters], Ukr. Mat. Zh. 12 (1960), 324-338.
 
-Gavrilyuk, V.T., ['Generalization of the first polynomial algorithm of 
-E.Ya.Remez for the problem of constructing rational-fractional 
+Gavrilyuk, V.T., ['Generalization of the first polynomial algorithm of
+E.Ya.Remez for the problem of constructing rational-fractional
 Chebyshev approximations], Ukr. Mat. Zh. 16 (1961), 575-585.
 
 Some English language sources include:
 
-Fraser, W., Hart, J.F., ['On the computation of rational approximations 
+Fraser, W., Hart, J.F., ['On the computation of rational approximations
 to continuous functions], Comm. of the ACM 5 (1962), 401-403, 414.
 
-Ralston, A., ['Rational Chebyshev approximation by Remes' algorithms], 
+Ralston, A., ['Rational Chebyshev approximation by Remes' algorithms],
 Numer.Math. 7 (1965), no. 4, 322-330.
 
-A. Ralston, ['Rational Chebyshev approximation, Mathematical 
-Methods for Digital Computers v. 2] (Ralston A., Wilf H., eds.), 
+A. Ralston, ['Rational Chebyshev approximation, Mathematical
+Methods for Digital Computers v. 2] (Ralston A., Wilf H., eds.),
 Wiley, New York, 1967, pp. 264-284.
 
 Hart, J.F. e.a., ['Computer approximations], Wiley, New York a.o., 1968.
 
-Cody, W.J., Fraser, W., Hart, J.F., ['Rational Chebyshev approximation 
+Cody, W.J., Fraser, W., Hart, J.F., ['Rational Chebyshev approximation
 using linear equations], Numer.Math. 12 (1968), 242-251.
 
-Cody, W.J., ['A survey of practical rational and polynomial 
+Cody, W.J., ['A survey of practical rational and polynomial
 approximation of functions], SIAM Review 12 (1970), no. 3, 400-423.
 
-Barrar, R.B., Loeb, H.J., ['On the Remez algorithm for non-linear 
+Barrar, R.B., Loeb, H.J., ['On the Remez algorithm for non-linear
 families], Numer.Math. 15 (1970), 382-391.
 
-Dunham, Ch.B., ['Convergence of the Fraser-Hart algorithm for rational 
+Dunham, Ch.B., ['Convergence of the Fraser-Hart algorithm for rational
 Chebyshev approximation], Math. Comp. 29 (1975), no. 132, 1078-1082.
 
 G. L. Litvinov, ['Approximate construction of rational
@@ -368,7 +368,7 @@ Russian Journal of Mathematical Physics, vol.1, No. 3, 1994.
 
 [endsect][/section:remez The Remez Method]
 
-[/ 
+[/
   Copyright 2006 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

doc/differentiation/numerical_differentiation.qbk

@@ -0,0 +1,108 @@
+[/
+Copyright (c) 2018 Nick Thompson
+Use, modification and distribution are subject to the
+Boost Software License, Version 1.0. (See accompanying file
+LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+]
+
+[section:diff Numerical Differentiation]
+
+[heading Synopsis]
+
+``
+#include <boost/math/tools/numerical_differentiation.hpp>
+
+
+namespace boost::math::tools {
+
+    template <class F, class Real>
+    Real complex_step_derivative(const F f, Real x);
+
+    template <class F, class Real, size_t order = 6>
+    Real finite_difference_derivative(const F f, Real x, Real* error = nullptr);
+
+} // namespaces
+``
+
+[heading Description]
+
+The function `finite_difference_derivative` calculates a finite-difference approximation to the derivative of of a function /f/ at point /x/.
+A basic usage is
+
+    auto f = [](double x) { return std::exp(x); };
+    double x = 1.7;
+    double dfdx = finite_difference_derivative(f, x);
+
+Finite differencing is complicated, as finite-difference approximations to the derivative are /infinitely/ ill-conditioned.
+In addition, for any function implemented in finite-precision arithmetic, the "true" derivative is /zero/ almost everywhere, and undefined at representables.
+However, some tricks allow for reasonable results to be obtained in many cases.
+
+There are two sources of error from finite differences: One, the truncation error arising from using a finite number of samples to cancel out higher order terms in the Taylor series.
+The second is the roundoff error involved in evaluating the function.
+The truncation error goes to zero as /h/ \u2192 0, but the roundoff error becomes unbounded.
+By balancing these two sources of error, we can choose a value of /h/ that minimizes the maximum total error.
+For this reason boost's `finite_difference_derivative` does not require the user to input a stepsize.
+For more details about the theoretical error analysis involved in finite-difference approximations to the derivative, see [@http://web.archive.org/web/20150420195907/http://www.uio.no/studier/emner/matnat/math/MAT-INF1100/h08/kompendiet/diffint.pdf here].
+
+Despite the effort that has went into choosing a reasonable value of /h/, the problem is still fundamentally ill-conditioned, and hence an error estimate is essential.
+It can be queried as follows
+
+    double error_estimate;
+    double d = finite_difference_derivative(f, x, &error_estimate);
+
+N.B.: Producing an error estimate requires additional function evaluations and as such is slower than simple evaluation of the derivative.
+It also expands the domain over which the function must be differentiable and requires the function to have two more continuous derivatives.
+The error estimate is computed under the assumption that /f/ is evaluated to 1ULP.
+This might seem an extreme assumption, but it is the only sensible one, as the routine cannot know the functions rounding error.
+If the function cannot be evaluated with very great accuracy, Lanczos's smoothing differentiation is recommended as an alternative.
+
+The default order of accuracy is 6, which reflects that fact that people tend to be interested in functions with many continuous derivatives.
+If your function does not have 7 continuous derivatives, is may be of interest to use a lower order method, which can be achieved via (say)
+
+    double d = finite_difference_derivative<decltype(f), Real, 2>(f, x);
+
+This requests a second-order accurate derivative be computed.
+
+It is emphatically /not/ the case that higher order methods always give higher accuracy for smooth functions.
+Higher order methods require more addition of positive and negative terms, which can lead to catastrophic cancellation.
+A function which is very good at making a mockery of finite-difference differentiation is exp(x)/(cos(x)[super 3] + sin(x)[super 3]).
+Differentiating this function by `finite_difference_derivative` in double precision at /x=5.5/ gives zero correct digits at order 4, 6, and 8, but recovers 5 correct digits at order 2.
+These are dangerous waters; use the error estimates to tread carefully.
+
+For a finite-difference method of order /k/, the error is /C/ \u03B5[super k/k+1].
+In the limit /k/ \u2192 \u221E, we see that the error tends to \u03B5, recovering the full precision for the type.
+However, this ignores the fact that higher-order methods require subtracting more nearly-equal (perhaps noisy) terms, so the constant /C/ grows with /k/.
+Since /C/ grows quickly and \u03B5[super k/k+1] approaches \u03B5 slowly, we can see there is a compromise between high-order accuracy and conditioning of the difference quotient.
+In practice we have found that /k=6/ seems to be a good compromise between the two (and have made this the default), but users are encouraged to examine the error estimates to choose an optimal order of accuracy for the given problem.
+
+[table:id Cost of Finite-Difference Numerical Differentiation
+    [[Order of Accuracy] [Function Evaluations] [Error]            [Continuous Derivatives Required for Error Estimate to Hold] [Additional Function Evaluations to Produce Error Estimates]]
+    [[1]                 [2]                    [\u03B5[super 1/2]]  [2]                                                     [1]]
+    [[2]                 [2]                    [\u03B5[super 2/3]]  [3]                                                     [2]]
+    [[4]                 [4]                    [\u03B5[super 4/5]]  [5]                                                     [2]]
+    [[6]                 [6]                    [\u03B5[super 6/7]]  [7]                                                     [2]]
+    [[8]                 [8]                    [\u03B5[super 8/9]]  [9]                                                     [2]]
+]
+
+
+Given all the caveats which must be kept in mind for successful use of finite-difference differentiation, it is reasonable to try to avoid it if possible.
+Boost provides two possibilities: The Chebyshev transform (see [link math_toolkit.sf_poly.chebyshev here]) and the complex step derivative.
+If your function is the restriction to the real line of a holomorphic function which takes real values at real argument, then the *complex step derivative* can be used.
+The idea is very simple: Since /f/ is complex-differentiable, ['f(x+\u2148 h) = f(x) + \u2148 hf'(x) - h[super 2]f''(x) + [bigo](h[super 3])].
+As long as /f(x)/ \u2208 \u211D, then ['f'(x) = \u2111 f(x+\u2148 h)/h + [bigo](h[super 2])].
+This method requires a single complex function evaluation and is not subject to the catastrophic subtractive cancellation that plagues finite-difference calculations.
+
+An example usage:
+
+    double x = 7.2;
+    double e_prime = complex_step_derivative(std::exp<std::complex<double>>, x);
+
+References:
+
+Squire, William, and George Trapp. ['Using complex variables to estimate derivatives of real functions.] Siam Review 40.1 (1998): 110-112.
+
+Fornberg, Bengt. ['Generation of finite difference formulas on arbitrarily spaced grids.] Mathematics of computation 51.184 (1988): 699-706.
+
+Corless, Robert M., and Nicolas Fillion. ['A graduate introduction to numerical methods.] AMC 10 (2013): 12.
+
+[endsect]

doc/distributions/gamma.qbk

@@ -101,7 +101,7 @@ The domain of the random variable is \[0,+[infin]\].
 
 [h4 Accuracy]
 
-The lognormal distribution is implemented in terms of the 
+The gamma distribution is implemented in terms of the 
 incomplete gamma functions __gamma_p and __gamma_q and their
 inverses __gamma_p_inv and __gamma_q_inv: refer to the accuracy
 data for those functions for more information.

doc/distributions/non_members.qbk

@@ -1,11 +1,11 @@
 [section:nmp Non-Member Properties]
 
-Properties that are common to all distributions are accessed via non-member 
+Properties that are common to all distributions are accessed via non-member
 getter functions: non-membership allows more of these functions to be added over time,
 as the need arises.  Unfortunately the literature uses many different and
 confusing names to refer to a rather small number of actual concepts; refer
-to the [link math_toolkit.dist_ref.nmp.concept_index concept index] to find the property you 
-want by the name you are most familiar with. 
+to the [link math_toolkit.dist_ref.nmp.concept_index concept index] to find the property you
+want by the name you are most familiar with.
 Or use the [link math_toolkit.dist_ref.nmp.function_index function index]
 to go straight to the function you want if you already know its name.
 
@@ -62,8 +62,8 @@ to go straight to the function you want if you already know its name.
 
    template <class RealType, class ``__Policy``>
    RealType cdf(const ``['Distribution-Type]``<RealType, ``__Policy``>& dist, const RealType& x);
-   
-The __cdf is the probability that 
+
+The __cdf is the probability that
 the variable takes a value less than or equal to x.  It is equivalent
 to the integral from -infinity to x of the __pdf.
 
@@ -79,11 +79,11 @@ normal distribution:
 
    template <class Distribution, class RealType>
    RealType cdf(const ``['Unspecified-Complement-Type]``<Distribution, RealType>& comp);
-   
-The complement of the __cdf 
-is the probability that 
+
+The complement of the __cdf
+is the probability that
 the variable takes a value greater than x.  It is equivalent
-to the integral from x to infinity of the __pdf, or 1 minus the __cdf of x. 
+to the integral from x to infinity of the __pdf, or 1 minus the __cdf of x.
 
 This is also known as the survival function.
 
@@ -118,7 +118,7 @@ the defined range for the distribution.
 [equation hazard]
 
 [caution
-Some authors refer to this as the conditional failure 
+Some authors refer to this as the conditional failure
 density function rather than the hazard function.]
 
 [h4:chf Cumulative Hazard Function]
@@ -133,14 +133,14 @@ the defined range for the distribution.
 
 [equation chf]
 
-[caution 
+[caution
 Some authors refer to this as simply the "Hazard Function".]
 
 [h4:mean mean]
 
    template<class RealType, class ``__Policy``>
    RealType mean(const ``['Distribution-Type]``<RealType, ``__Policy``>& dist);
-   
+
 Returns the mean of the distribution /dist/.
 
 This function may return a __domain_error if the distribution does not have
@@ -150,14 +150,14 @@ a defined mean (for example the Cauchy distribution).
 
    template<class RealType, class ``__Policy``>
    RealType median(const ``['Distribution-Type]``<RealType, ``__Policy``>& dist);
-   
+
 Returns the median of the distribution /dist/.
 
 [h4:mode mode]
 
    template<class RealType, ``__Policy``>
    RealType mode(const ``['Distribution-Type]``<RealType, ``__Policy``>& dist);
-   
+
 Returns the mode of the distribution /dist/.
 
 This function may return a __domain_error if the distribution does not have
@@ -167,14 +167,14 @@ a defined mode.
 
    template <class RealType, class ``__Policy``>
    RealType pdf(const ``['Distribution-Type]``<RealType, ``__Policy``>& dist, const RealType& x);
-   
-For a continuous function, the probability density function (pdf) returns 
-the probability that the variate has the value x. 
-Since for continuous distributions the probability at a single point is actually zero, 
+
+For a continuous function, the probability density function (pdf) returns
+the probability that the variate has the value x.
+Since for continuous distributions the probability at a single point is actually zero,
 the probability is better expressed as the integral of the pdf between two points:
 see the __cdf.
 
-For a discrete distribution, the pdf is the probability that the 
+For a discrete distribution, the pdf is the probability that the
 variate takes the value x.
 
 This function may return a __domain_error if the random variable is outside
@@ -188,14 +188,14 @@ For example, for a standard normal distribution the pdf looks like this:
 
    template<class RealType, class ``__Policy``>
    std::pair<RealType, RealType> range(const ``['Distribution-Type]``<RealType, ``__Policy``>& dist);
-   
+
 Returns the valid range of the random variable over distribution /dist/.
 
 [h4:quantile Quantile]
 
    template <class RealType, class ``__Policy``>
    RealType quantile(const ``['Distribution-Type]``<RealType, ``__Policy``>& dist, const RealType& p);
-   
+
 The quantile is best viewed as the inverse of the __cdf, it returns
 a value /x/ such that `cdf(dist, x) == p`.
 
@@ -217,7 +217,7 @@ See also [link math_toolkit.stat_tut.overview.complements complements].
 
    template <class Distribution, class RealType>
    RealType quantile(const ``['Unspecified-Complement-Type]``<Distribution, RealType>& comp);
-   
+
 This is the inverse of the __ccdf.  It is calculated by wrapping
 the arguments in a call to the quantile function in a call to
 /complement/.  For example:
@@ -250,8 +250,8 @@ distribution:
 
    template <class RealType, class ``__Policy``>
    RealType standard_deviation(const ``['Distribution-Type]``<RealType, ``__Policy``>& dist);
-   
-Returns the standard deviation of distribution /dist/.   
+
+Returns the standard deviation of distribution /dist/.
 
 This function may return a __domain_error if the distribution does not have
 a defined standard deviation.
@@ -260,7 +260,7 @@ a defined standard deviation.
 
    template<class RealType, class ``__Policy``>
    std::pair<RealType, RealType> support(const ``['Distribution-Type]``<RealType, ``__Policy``>& dist);
-   
+
 Returns the supported range of random variable over the distribution /dist/.
 
 The distribution is said to be 'supported' over a range that is
@@ -274,7 +274,7 @@ Outside are uninteresting zones where the pdf is zero, and the cdf zero or unity
 
    template <class RealType, class ``__Policy``>
    RealType variance(const ``['Distribution-Type]``<RealType, ``__Policy``>& dist);
-   
+
 Returns the variance of the distribution /dist/.
 
 This function may return a __domain_error if the distribution does not have
@@ -284,7 +284,7 @@ a defined variance.
 
    template <class RealType, class ``__Policy``>
    RealType skewness(const ``['Distribution-Type]``<RealType, ``__Policy``>& dist);
-   
+
 Returns the skewness of the distribution /dist/.
 
 This function may return a __domain_error if the distribution does not have
@@ -294,7 +294,7 @@ a defined skewness.
 
    template <class RealType, class ``__Policy``>
    RealType kurtosis(const ``['Distribution-Type]``<RealType, ``__Policy``>& dist);
-   
+
 Returns the 'proper' kurtosis (normalized fourth moment) of the distribution /dist/.
 
 kertosis = [beta][sub 2][space]= [mu][sub 4][space] / [mu][sub 2][super 2]
@@ -326,7 +326,7 @@ a defined kurtosis.
 
    template <class RealType, ``__Policy``>
    RealType kurtosis_excess(const ``['Distribution-Type]``<RealType, ``__Policy``>& dist);
-   
+
 Returns the kurtosis excess of the distribution /dist/.
 
 kurtosis excess = [gamma][sub 2][space]= [mu][sub 4][space] / [mu][sub 2][super 2][space]- 3 = kurtosis - 3
@@ -334,7 +334,7 @@ kurtosis excess = [gamma][sub 2][space]= [mu][sub 4][space] / [mu][sub 2][super
 Where [mu][sub i][space] is the i'th central moment of the distribution, and
 in particular [mu][sub 2][space] is the variance of the distribution.
 
-The kurtosis excess is a measure of the "peakedness" of a distribution, and 
+The kurtosis excess is a measure of the "peakedness" of a distribution, and
 is more widely used than the "kurtosis proper".  It is defined so that
 the kurtosis excess of a normal distribution is zero.
 
@@ -344,7 +344,7 @@ a defined kurtosis excess.
 Kurtosis excess can have a value from -2 to + infinity.
 
   kurtosis = kurtosis_excess +3;
-  
+
 The kurtosis excess of a normal distribution is zero.
 
 [h4:cdfPQ P and Q]
@@ -361,12 +361,12 @@ the __quantile.
 
 [h4:cdf_inv Inverse CDF Function.]
 
-The inverse of the cumulative distribution function, is the same as the 
+The inverse of the cumulative distribution function, is the same as the
 __quantile.
 
 [h4:survival_inv Inverse Survival Function.]
 
-The inverse of the survival function, is the same as computing the 
+The inverse of the survival function, is the same as computing the
 [link math_toolkit.dist_ref.nmp.quantile_c quantile
 from the complement of the probability].
 
@@ -380,13 +380,13 @@ while the term __pdf applies to continuous distributions.
 [h4:lower_critical Lower Critical Value.]
 
 The lower critical value calculates the value of the random variable
-given the area under the left tail of the distribution.  
+given the area under the left tail of the distribution.
 It is equivalent to calculating the __quantile.
 
-[h4: upper_critical Upper Critical Value.]
+[h4:upper_critical Upper Critical Value.]
 
 The upper critical value calculates the value of the random variable
-given the area under the right tail of the distribution.  It is equivalent to 
+given the area under the right tail of the distribution.  It is equivalent to
 calculating the [link math_toolkit.dist_ref.nmp.quantile_c quantile from the complement of the
 probability].
 
@@ -394,8 +394,7 @@ probability].
 
 Refer to the __ccdf.
 
-[endsect][/section:nmp Non-Member Properties]
-
+[endsect] [/section:nmp Non-Member Properties]
 
 [/ non_members.qbk
   Copyright 2006 John Maddock and Paul A. Bristow.

doc/distributions/students_t_examples.qbk

@@ -611,7 +611,7 @@ all the possible tests are summed up in the following table:
 *difference* in means]
 [Reject if complement of CDF for |t| > significance level / 2:
 
-`cdf(complement(dist, fabs(t))) < alpha / 2`]]
+`cdf(complement(dist, fabs(t))) > alpha / 2`]]
 
 [[The Alternative-hypothesis: Sample 1 Mean is *less* than
 Sample 2 Mean.]

doc/equations/complex_bessel_k_integral.mml

@@ -0,0 +1,51 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.1 plus MathML 2.0//EN" "http://www.w3.org/TR/MathML2/dtd/xhtml-math11-f.dtd" [<!ENTITY mathml 'http://www.w3.org/1998/Math/MathML'>]>
+<html xmlns="http://www.w3.org/1999/xhtml">
+<head><title>complex_bessel_k_integral</title>
+<!-- MathML created with MathCast Equation Editor version 0.92 -->
+</head>
+<body>
+<math xmlns="http://www.w3.org/1998/Math/MathML" display="block">
+  <mrow>
+    <msub>
+      <mi>K</mi>
+      <mi>&#x3b1;</mi>
+    </msub>
+    <mfenced>
+      <mrow>
+        <mi>x</mi>
+      </mrow>
+    </mfenced>
+    <mo>=</mo>
+    <msubsup>
+      <mo>&#x222b;</mo>
+      <mrow>
+        <mi>o</mi>
+      </mrow>
+      <mi>&#x221e;</mi>
+    </msubsup>
+    <msup>
+      <mi>e</mi>
+      <mrow>
+        <mo>&#x2212;</mo>
+        <mi>x</mi>
+        <mtext>cosh</mtext>
+        <mfenced>
+          <mrow>
+            <mi>t</mi>
+          </mrow>
+        </mfenced>
+      </mrow>
+    </msup>
+    <mtext>cosh</mtext>
+    <mfenced>
+      <mrow>
+        <mi>&#x3b1;</mi>
+        <mi>t</mi>
+      </mrow>
+    </mfenced>
+    <mi>d</mi>
+    <mi>t</mi>
+  </mrow>
+</math></body>
+</html>

doc/equations/complex_eta_integral.mml

@@ -0,0 +1,70 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.1 plus MathML 2.0//EN" "http://www.w3.org/TR/MathML2/dtd/xhtml-math11-f.dtd" [<!ENTITY mathml 'http://www.w3.org/1998/Math/MathML'>]>
+<html xmlns="http://www.w3.org/1999/xhtml">
+<head><title>complex_eta_integral</title>
+<!-- MathML created with MathCast Equation Editor version 0.92 -->
+</head>
+<body>
+<math xmlns="http://www.w3.org/1998/Math/MathML" display="block">
+  <mrow>
+    <mi>&#x3b7;</mi>
+    <mfenced>
+      <mrow>
+        <mi>s</mi>
+      </mrow>
+    </mfenced>
+    <mo>=</mo>
+    <munderover>
+      <mo>&#x222b;</mo>
+      <mrow>
+        <mo>&#x2212;</mo>
+        <mi>&#x221e;</mi>
+      </mrow>
+      <mi>&#x221e;</mi>
+    </munderover>
+    <mfrac>
+      <mrow>
+        <mrow>
+          <msup>
+            <mfenced>
+              <mrow>
+                <mfrac>
+                  <mn>1</mn>
+                  <mn>2</mn>
+                </mfrac>
+                <mo>+</mo>
+                <mi>i</mi>
+                <mi>t</mi>
+              </mrow>
+            </mfenced>
+            <mrow>
+              <mo>&#x2212;</mo>
+              <mi>s</mi>
+            </mrow>
+          </msup>
+        </mrow>
+      </mrow>
+      <mrow>
+        <msup>
+          <mi>e</mi>
+          <mrow>
+            <mi>&#x3c0;</mi>
+            <mi>t</mi>
+          </mrow>
+        </msup>
+        <mo>+</mo>
+        <msup>
+          <mi>e</mi>
+          <mrow>
+            <mo>&#x2212;</mo>
+            <mi>&#x3c0;</mi>
+            <mi>t</mi>
+          </mrow>
+        </msup>
+      </mrow>
+    </mfrac>
+    <mi>d</mi>
+    <mi>t</mi>
+  </mrow>
+</math></body>
+</html>

doc/equations/generate.sh

@@ -7,7 +7,7 @@
 #
 # Paths to tools come first, change these to match your system:
 #
-math2svg='m:\download\open\SVGMath-0.3.1\math2svg.py'
+math2svg='d:\download\open\SVGMath-0.3.3\math2svg.py'
 python='/cygdrive/c/Python27/python.exe'
 inkscape='/cygdrive/c/Program Files/Inkscape/inkscape.exe'
 # Image DPI:

doc/equations/sine_integral.mml

@@ -0,0 +1,64 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.1 plus MathML 2.0//EN" "http://www.w3.org/TR/MathML2/dtd/xhtml-math11-f.dtd" [<!ENTITY mathml 'http://www.w3.org/1998/Math/MathML'>]>
+<html xmlns="http://www.w3.org/1999/xhtml">
+<head><title>sine_integral</title>
+<!-- MathML created with MathCast Equation Editor version 0.92 -->
+</head>
+<body>
+<math xmlns="http://www.w3.org/1998/Math/MathML" display="block">
+  <mrow>
+    <mi>S</mi>
+    <mi>i</mi>
+    <mfenced>
+      <mrow>
+        <mi>z</mi>
+      </mrow>
+    </mfenced>
+    <mo>=</mo>
+    <mfrac>
+      <mi>&#x3c0;</mi>
+      <mn>2</mn>
+    </mfrac>
+    <mo>&#x2212;</mo>
+    <msubsup>
+      <mo>&#x222b;</mo>
+      <mrow>
+        <mn>0</mn>
+      </mrow>
+      <mrow>
+        <mfrac>
+          <mi>&#x3c0;</mi>
+          <mn>2</mn>
+        </mfrac>
+      </mrow>
+    </msubsup>
+    <msup>
+      <mi>e</mi>
+      <mrow>
+        <mo>&#x2212;</mo>
+        <mi>z</mi>
+        <mi>cos</mi>
+        <mfenced>
+          <mrow>
+            <mi>t</mi>
+          </mrow>
+        </mfenced>
+      </mrow>
+    </msup>
+    <mi>cos</mi>
+    <mfenced>
+      <mrow>
+        <mi>z</mi>
+        <mi>sin</mi>
+        <mfenced>
+          <mrow>
+            <mi>t</mi>
+          </mrow>
+        </mfenced>
+      </mrow>
+    </mfenced>
+    <mi>d</mi>
+    <mi>t</mi>
+  </mrow>
+</math></body>
+</html>

doc/fp_utilities/float_next.qbk

@@ -20,6 +20,26 @@ __multiprecision __cpp_dec_float and__cpp_bin_float  are fixed at runtime,
 but [*not] a type that extends the representation to provide an exact representation
 for any number, for example [@http://keithbriggs.info/xrc.html XRC eXact Real in C].
 
+The accuracy of mathematical functions can be assessed and displayed in terms of __ULP,
+often as a ulps plot or by binning the differences as a histogram.
+Samples are evaluated using the implementation under test and compared with 'known good'
+representation obtained using a more accurate method.  Other implementations, often using
+arbitrary precision arithmetic, for example __WolframAlpha are one source of references
+values.  The other method, used widely in Boost.Math special functions, it to carry out
+the same algorithm, but using a higher precision type, typically using Boost.Multiprecision
+types like `cpp_bin_float_quad` for 128-bit (about 35 decimal digit precision), or
+`cpp_bin_float_50` (for 50 decimal digit precision).
+
+When converted to a particular machine representation, say `double`, say using a `static_cast`,
+the value is the nearest representation possible for the `double` type.  This value
+cannot be 'wrong' by more than half a __ulp, and can be obtained using the Boost.Math function `ulp`.
+(Unless the algorithm is fundamentally flawed, something that should be revealed by 'sanity'
+checks using some independent sources).
+
+See some discussion and example plots by Cleve Moler of Mathworks
+[@https://blogs.mathworks.com/cleve/2017/01/23/ulps-plots-reveal-math-function-accurary/
+ulps plots reveal math-function accuracy].
+
 [section:nextafter Finding the Next Representable Value in a Specific Direction (nextafter)]
 
 [h4 Synopsis]
@@ -220,13 +240,13 @@ Function `ulp` gives the size of a unit-in-the-last-place for a specified floati
 
 Returns one [@http://en.wikipedia.org/wiki/Unit_in_the_last_place unit in the last place] of ['x].
 
-Corner cases are handled as followes:
+Corner cases are handled as follows:
 
 * If the argument is a NaN, then raises a __domain_error.
 * If the argument is an infinity, then raises an __overflow_error.
 * If the argument is zero then returns the smallest representable value: for example for type
 `double` this would be either `std::numeric_limits<double>::min()` or `std::numeric_limits<double>::denorm_min()`
-depending whether denormals are supported (which have the values 2.`2250738585072014e-308` and `4.9406564584124654e-324` respectively).
+depending whether denormals are supported (which have the values `2.2250738585072014e-308` and `4.9406564584124654e-324` respectively).
 * If the result is too small to represent, then returns the smallest representable value.
 * Always returns a positive value such that `ulp(x) == ulp(-x)`.
 

doc/graphs/diode_iv_plot.svg

@@ -0,0 +1,81 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<svg width="400" height ="300" version="1.1"
+xmlns:svg ="http://www.w3.org/2000/svg"
+xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
+xmlns:cc="http://web.resource.org/cc/"
+xmlns:dc="http://purl.org/dc/elements/1.1/"
+xmlns ="http://www.w3.org/2000/svg"
+>
+<!-- SVG plot written using Boost.Plot program (Creator Jacob Voytko) --> 
+<!-- Use, modification and distribution of Boost.Plot subject to the --> 
+<!-- Boost Software License, Version 1.0.--> 
+<!-- (See accompanying file LICENSE_1_0.txt --> 
+<!-- or copy at http://www.boost.org/LICENSE_1_0.txt) --> 
+
+<!-- SVG Plot Copyright Paul A. Bristow 2016 --> 
+<meta name="copyright" content="Paul A. Bristow" />
+<meta name="date" content="2016" />
+<clipPath id="plot_window"><rect x="85.2" y="54" width="151.1" height="186"/></clipPath>
+<g id="imageBackground" stroke="rgb(255,255,0)" fill="rgb(255,255,255)" stroke-width="2"><rect x="0" y="0" width="400" height="300"/></g>
+<g id="plotBackground" stroke="rgb(119,136,153)" fill="rgb(255,255,255)" stroke-width="2"><rect x="84.2" y="53" width="153.1" height="188"/></g>
+<g id="yMinorGrid" stroke="rgb(200,220,255)" stroke-width="0.5"></g>
+<g id="yMajorGrid" stroke="rgb(200,220,255)" stroke-width="1"><path d="M86.2,84.33 L235.3,84.33 M86.2,53 L235.3,53 M86.2,84.33 L235.3,84.33 M86.2,115.7 L235.3,115.7 M86.2,147 L235.3,147 M86.2,178.3 L235.3,178.3 M86.2,209.7 L235.3,209.7 M86.2,241 L235.3,241 " fill="none"/></g>
+<g id="xMinorGrid" stroke="rgb(200,220,255)" stroke-width="0.5"></g>
+<g id="xMajorGrid" stroke="rgb(200,220,255)" stroke-width="1"><path d="M94.41,53 L94.41,241 M114.8,53 L114.8,241 M135.2,53 L135.2,241 M155.6,53 L155.6,241 M176,53 L176,241 M196.5,53 L196.5,241 M216.9,53 L216.9,241 M237.3,53 L237.3,241 " fill="none"/></g>
+<g id="yAxis" stroke="rgb(0,0,0)" stroke-width="1"><line x1="84.2" y1="53" x2="84.2" y2="241"/></g>
+<g id="xAxis" stroke="rgb(0,0,0)" stroke-width="1"><line x1="79.2" y1="84.33" x2="237.3" y2="84.33"/><line x1="79.2" y1="241" x2="237.3" y2="241"/></g>
+<g id="yMinorTicks" stroke="rgb(0,0,0)" stroke-width="1"><path d="M82.2,78.07 L84.2,78.07 M82.2,71.8 L84.2,71.8 M82.2,65.53 L84.2,65.53 M82.2,59.27 L84.2,59.27 M82.2,84.33 L84.2,84.33 M82.2,90.6 L84.2,90.6 M82.2,96.87 L84.2,96.87 M82.2,103.1 L84.2,103.1 M82.2,109.4 L84.2,109.4 M82.2,115.7 L84.2,115.7 M82.2,121.9 L84.2,121.9 M82.2,128.2 L84.2,128.2 M82.2,134.5 L84.2,134.5 M82.2,140.7 L84.2,140.7 M82.2,147 L84.2,147 M82.2,147 L84.2,147 M82.2,153.3 L84.2,153.3 M82.2,159.5 L84.2,159.5 M82.2,165.8 L84.2,165.8 M82.2,172.1 L84.2,172.1 M82.2,178.3 L84.2,178.3 M82.2,184.6 L84.2,184.6 M82.2,190.9 L84.2,190.9 M82.2,197.1 L84.2,197.1 M82.2,203.4 L84.2,203.4 M82.2,209.7 L84.2,209.7 M82.2,215.9 L84.2,215.9 M82.2,222.2 L84.2,222.2 M82.2,228.5 L84.2,228.5 M82.2,234.7 L84.2,234.7 M82.2,241 L84.2,241 " fill="none"/></g>
+<g id="xMinorTicks" stroke="rgb(0,0,0)" stroke-width="1"><path d="M86.24,241 L86.24,243 M90.32,241 L90.32,243 M98.49,241 L98.49,243 M102.6,241 L102.6,243 M106.7,241 L106.7,243 M110.7,241 L110.7,243 M118.9,241 L118.9,243 M123,241 L123,243 M127.1,241 L127.1,243 M131.1,241 L131.1,243 M139.3,241 L139.3,243 M143.4,241 L143.4,243 M147.5,241 L147.5,243 M151.6,241 L151.6,243 M159.7,241 L159.7,243 M163.8,241 L163.8,243 M167.9,241 L167.9,243 M172,241 L172,243 M180.1,241 L180.1,243 M184.2,241 L184.2,243 M188.3,241 L188.3,243 M192.4,241 L192.4,243 M200.5,241 L200.5,243 M204.6,241 L204.6,243 M208.7,241 L208.7,243 M212.8,241 L212.8,243 M221,241 L221,243 M225,241 L225,243 M229.1,241 L229.1,243 M233.2,241 L233.2,243 " fill="none"/></g>
+<g id="yMajorTicks" stroke="rgb(0,0,0)" stroke-width="2"><path d="M79.2,84.33 L84.2,84.33 M79.2,53 L84.2,53 M79.2,84.33 L84.2,84.33 M79.2,115.7 L84.2,115.7 M79.2,147 L84.2,147 M79.2,178.3 L84.2,178.3 M79.2,209.7 L84.2,209.7 M79.2,241 L84.2,241 " fill="none"/></g>
+<g id="xMajorTicks" stroke="rgb(0,0,0)" stroke-width="2"><path d="M94.41,241 L94.41,246 M114.8,241 L114.8,246 M135.2,241 L135.2,246 M155.6,241 L155.6,246 M176,241 L176,246 M196.5,241 L196.5,246 M216.9,241 L216.9,246 M237.3,241 L237.3,246 " fill="none"/></g>
+<g id="xTicksValues">
+<text x="94.41" y="261.6" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">0.3</text>
+<text x="114.8" y="261.6" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">0.4</text>
+<text x="135.2" y="261.6" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">0.5</text>
+<text x="155.6" y="261.6" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">0.6</text>
+<text x="176" y="261.6" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">0.7</text>
+<text x="196.5" y="261.6" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">0.8</text>
+<text x="216.9" y="261.6" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">0.9</text>
+<text x="237.3" y="261.6" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">1</text></g>
+<g id="yTicksValues">
+<text x="73.2" y="86.73" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">0</text>
+<text x="73.2" y="55.4" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">4</text>
+<text x="73.2" y="86.73" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">0</text>
+<text x="73.2" y="118.1" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">-4</text>
+<text x="73.2" y="149.4" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">-8</text>
+<text x="73.2" y="180.7" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">-12</text>
+<text x="73.2" y="212.1" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">-16</text>
+<text x="73.2" y="243.4" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">-20</text></g>
+<g id="yLabel">
+<text x="42.9" y="147" text-anchor="middle" transform = "rotate(-90 42.9 147 )" font-size="14" font-family="Lucida Sans Unicode">log(current) (A)</text></g>
+<g id="xLabel">
+<text x="160.7" y="276.7" text-anchor="middle" font-size="14" font-family="Lucida Sans Unicode">voltage (V)</text></g>
+<g id="plotLines" stroke-width="2"><g clip-path="url(#plot_window)" stroke="rgb(0,0,0)" stroke-width="0.5"><path d="M94.41,238.3 L94.41,238.3 L114.8,207.7 L135.2,177.2 L155.6,146.8 L176,116.3 L196.5,84.49 L216.9,56.13 " fill="none"/></g>
+<g clip-path="url(#plot_window)" stroke="rgb(0,0,0)" stroke-width="0.5"><path d="M94.41,238.3 L94.41,238.3 L114.8,207.7 L135.2,177.2 L155.6,146.8 L176,117.2 L180.1,111.8 L184.2,106.3 L188.3,102.3 L192.4,100 " fill="none"/></g>
+<g clip-path="url(#plot_window)" stroke="rgb(0,0,255)" stroke-width="1"><path d="M94.41,238.2 L94.41,238.2 L96.45,235.2 L98.49,232.1 L100.5,229.1 L102.6,226.1 L104.6,223 L106.7,220 L108.7,216.9 L110.7,213.9 L112.8,210.8 L114.8,207.8 L116.9,204.7 L118.9,201.7 L120.9,198.6 L123,195.6 L125,192.5 L127.1,189.5 L129.1,186.4 L131.1,183.4 L133.2,180.3 L135.2,177.3 L137.3,174.2 L139.3,171.2 L141.3,168.2 L143.4,165.1 L145.4,162.1 L147.5,159.1 L149.5,156 L151.6,153 L153.6,150 L155.6,147 L157.7,144.1 L159.7,141.2 L161.8,138.4 L163.8,135.6 L165.8,133 L167.9,130.5 L169.9,128.1 L172,125.9 L174,123.9 L176,122 L178.1,120.3 L180.1,118.8 L182.2,117.4 L184.2,116.2 L186.3,115 L188.3,114 L190.3,113 L192.4,112.1 L194.4,111.3 L196.5,110.5 L198.5,109.8 L200.5,109.2 L202.6,108.6 L204.6,108 L206.7,107.4 L208.7,106.9 L210.7,106.4 L212.8,105.9 L214.8,105.5 L216.9,105.1 L218.9,104.7 L221,104.3 L223,103.9 L225,103.5 L227.1,103.2 L229.1,102.9 L231.2,102.6 L233.2,102.2 L235.2,101.9 " fill="none"/></g>
+<g clip-path="url(#plot_window)" stroke="rgb(128,0,128)" stroke-width="1"><path d="M94.41,238.2 L94.41,238.2 L96.45,235.2 L98.49,232.1 L100.5,229.1 L102.6,226.1 L104.6,223 L106.7,220 L108.7,216.9 L110.7,213.9 L112.8,210.8 L114.8,207.8 L116.9,204.7 L118.9,201.7 L120.9,198.6 L123,195.6 L125,192.5 L127.1,189.5 L129.1,186.4 L131.1,183.4 L133.2,180.4 L135.2,177.3 L137.3,174.3 L139.3,171.2 L141.3,168.2 L143.4,165.2 L145.4,162.2 L147.5,159.2 L149.5,156.3 L151.6,153.4 L153.6,150.5 L155.6,147.8 L157.7,145.1 L159.7,142.6 L161.8,140.2 L163.8,138 L165.8,135.9 L167.9,134 L169.9,132.3 L172,130.8 L174,129.4 L176,128.1 L178.1,126.9 L180.1,125.9 L182.2,124.9 L184.2,124 L186.3,123.2 L188.3,122.4 L190.3,121.7 L192.4,121 L194.4,120.4 L196.5,119.8 L198.5,119.2 L200.5,118.7 L202.6,118.2 L204.6,117.7 L206.7,117.3 L208.7,116.9 L210.7,116.5 L212.8,116.1 L214.8,115.7 L216.9,115.3 L218.9,115 L221,114.7 L223,114.3 L225,114 L227.1,113.7 L229.1,113.4 L231.2,113.2 L233.2,112.9 L235.2,112.6 " fill="none"/></g>
+<g clip-path="url(#plot_window)" stroke="rgb(0,128,0)" stroke-width="1"><path d="M94.41,238.2 L94.41,238.2 L96.45,235.2 L98.49,232.1 L100.5,229.1 L102.6,226.1 L104.6,223 L106.7,220 L108.7,216.9 L110.7,213.9 L112.8,210.8 L114.8,207.8 L116.9,204.7 L118.9,201.7 L120.9,198.6 L123,195.6 L125,192.5 L127.1,189.5 L129.1,186.5 L131.1,183.4 L133.2,180.4 L135.2,177.4 L137.3,174.4 L139.3,171.4 L141.3,168.5 L143.4,165.6 L145.4,162.8 L147.5,160 L149.5,157.4 L151.6,154.8 L153.6,152.5 L155.6,150.3 L157.7,148.2 L159.7,146.4 L161.8,144.7 L163.8,143.2 L165.8,141.8 L167.9,140.5 L169.9,139.4 L172,138.3 L174,137.4 L176,136.5 L178.1,135.6 L180.1,134.9 L182.2,134.2 L184.2,133.5 L186.3,132.9 L188.3,132.3 L190.3,131.7 L192.4,131.2 L194.4,130.7 L196.5,130.3 L198.5,129.8 L200.5,129.4 L202.6,129 L204.6,128.6 L206.7,128.2 L208.7,127.9 L210.7,127.5 L212.8,127.2 L214.8,126.9 L216.9,126.6 L218.9,126.3 L221,126 L223,125.7 L225,125.4 L227.1,125.2 L229.1,124.9 L231.2,124.7 L233.2,124.4 L235.2,124.2 " fill="none"/></g>
+<g clip-path="url(#plot_window)" stroke="rgb(255,0,0)" stroke-width="1"><path d="M94.41,238.2 L94.41,238.2 L96.45,235.2 L98.49,232.1 L100.5,229.1 L102.6,226.1 L104.6,223 L106.7,220 L108.7,216.9 L110.7,213.9 L112.8,210.8 L114.8,207.8 L116.9,204.7 L118.9,201.7 L120.9,198.7 L123,195.6 L125,192.6 L127.1,189.6 L129.1,186.6 L131.1,183.6 L133.2,180.7 L135.2,177.8 L137.3,175 L139.3,172.2 L141.3,169.6 L143.4,167.1 L145.4,164.7 L147.5,162.5 L149.5,160.5 L151.6,158.7 L153.6,157 L155.6,155.5 L157.7,154.1 L159.7,152.8 L161.8,151.7 L163.8,150.6 L165.8,149.7 L167.9,148.8 L169.9,148 L172,147.2 L174,146.5 L176,145.9 L178.1,145.2 L180.1,144.6 L182.2,144.1 L184.2,143.6 L186.3,143.1 L188.3,142.6 L190.3,142.2 L192.4,141.8 L194.4,141.3 L196.5,141 L198.5,140.6 L200.5,140.2 L202.6,139.9 L204.6,139.6 L206.7,139.2 L208.7,138.9 L210.7,138.6 L212.8,138.4 L214.8,138.1 L216.9,137.8 L218.9,137.5 L221,137.3 L223,137 L225,136.8 L227.1,136.6 L229.1,136.4 L231.2,136.1 L233.2,135.9 L235.2,135.7 " fill="none"/></g>
+</g>
+<g id="plotPoints" clip-path="url(#plot_window)"><g stroke="rgb(0,0,0)" fill="rgb(211,211,211)"><rect x="92.91" y="236.8" width="3" height="3"/><rect x="113.3" y="206.2" width="3" height="3"/><rect x="133.7" y="175.7" width="3" height="3"/><rect x="154.1" y="145.3" width="3" height="3"/><rect x="174.5" y="115.7" width="3" height="3"/><rect x="178.6" y="110.3" width="3" height="3"/><rect x="182.7" y="104.8" width="3" height="3"/><rect x="186.8" y="100.8" width="3" height="3"/><rect x="190.9" y="98.5" width="3" height="3"/></g>
+</g>
+<g id="legendBackground" stroke="rgb(255,255,0)" fill="rgb(255,255,255)" stroke-width="1"><rect x="251.3" y="53" width="122.7" height="182"/><rect x="251.3" y="53" width="122.7" height="182"/></g>
+<g id="legendPoints"><g stroke="rgb(0,0,0)" fill="rgb(211,211,211)" stroke-width="0.5"><line x1="265.3" y1="81" x2="279.3" y2="81"/></g>
+<g stroke="rgb(0,0,0)" fill="rgb(211,211,211)" stroke-width="0.5"><rect x="263.8" y="107.5" width="3" height="3"/><line x1="286.3" y1="109" x2="300.3" y2="109"/></g>
+<g stroke="rgb(0,0,255)" fill="rgb(255,255,255)" stroke-width="1"><line x1="265.3" y1="137" x2="279.3" y2="137"/></g>
+<g stroke="rgb(128,0,128)" fill="rgb(255,255,255)" stroke-width="1"><line x1="265.3" y1="165" x2="279.3" y2="165"/></g>
+<g stroke="rgb(0,128,0)" fill="rgb(255,255,255)" stroke-width="1"><line x1="265.3" y1="193" x2="279.3" y2="193"/></g>
+<g stroke="rgb(255,0,0)" fill="rgb(255,255,255)" stroke-width="1"><line x1="265.3" y1="221" x2="279.3" y2="221"/></g>
+</g>
+<g id="legendText">
+<text x="286.3" y="81" font-size="14" font-family="Lucida Sans Unicode">I&#x2080;(V)</text>
+<text x="307.3" y="109" font-size="14" font-family="Lucida Sans Unicode">Rs=0 &#x3A9;</text>
+<text x="286.3" y="137" font-size="14" font-family="Lucida Sans Unicode">Rs=2 &#x3A9;</text>
+<text x="286.3" y="165" font-size="14" font-family="Lucida Sans Unicode">Rs=10 &#x3A9;</text>
+<text x="286.3" y="193" font-size="14" font-family="Lucida Sans Unicode">Rs=51 &#x3A9;</text>
+<text x="286.3" y="221" font-size="14" font-family="Lucida Sans Unicode">Rs=249 &#x3A9;</text></g>
+<g id="title">
+<text x="200" y="36" text-anchor="middle" font-size="18" font-family="Lucida Sans Unicode">Diode current versus voltage</text></g>
+<g id="plotXValues"></g>
+<g id="plotYValues"></g>
+</svg>

doc/graphs/dist_graphs.cpp

@@ -10,7 +10,7 @@
     \author John Maddock and Paul A. Bristow
   */
 //  Copyright John Maddock 2008.
-//  Copyright Paul A. Bristow 2008, 2009, 2012
+//  Copyright Paul A. Bristow 2008, 2009, 2012, 2016
 //  Use, modification and distribution are subject to the
 //  Boost Software License, Version 1.0. (See accompanying file
 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
@@ -35,8 +35,9 @@
 
 template <class Dist>
 struct is_discrete_distribution
-   : public boost::mpl::false_{};
+   : public boost::mpl::false_{}; // Default is continuous distribution.
 
+// Some discrete distributions.
 template<class T, class P>
 struct is_discrete_distribution<boost::math::bernoulli_distribution<T,P> >
    : public boost::mpl::true_{};
@@ -68,7 +69,7 @@ struct value_finder
 private:
    Dist m_dist;
    typename Dist::value_type m_value;
-};
+}; // value_finder
 
 template <class Dist>
 class distribution_plotter
@@ -177,7 +178,7 @@ public:
          if(b > m_max_x)
             m_max_x = b;
       }
-   }
+   } // add
 
    void plot(const std::string& title, const std::string& file)
    {
@@ -198,9 +199,12 @@ public:
          m_max_y = 1;
       }
 
+      std::cout << "Plotting " << title << " to " << file << std::endl;
+
       svg_2d_plot plot;
       plot.image_x_size(750);
       plot.image_y_size(400);
+      plot.copyright_holder("John Maddock").copyright_date("2008").boost_license_on(true);
       plot.coord_precision(4); // Avoids any visible steps.
       plot.title_font_size(20);
       plot.legend_title_font_size(15);
@@ -248,19 +252,17 @@ public:
 
       if(!is_discrete_distribution<Dist>::value)
       {
-         //
          // Continuous distribution:
-         //
          for(std::list<std::pair<std::string, Dist> >::const_iterator i = m_distributions.begin();
             i != m_distributions.end(); ++i)
          {
             double x = m_min_x;
-            double interval = (m_max_x - m_min_x) / 200;
+            double continuous_interval = (m_max_x - m_min_x) / 200;
             std::map<double, double> data;
             while(x <= m_max_x)
             {
                data[x] = m_pdf ? pdf(i->second, x) : cdf(i->second, x);
-               x += interval;
+               x += continuous_interval;
             }
             plot.plot(data, i->first)
                .line_on(true)
@@ -275,16 +277,14 @@ public:
       }
       else
       {
-         //
          // Discrete distribution:
-         //
          double x_width = 0.75 / m_distributions.size();
          double x_off = -0.5 * 0.75;
          for(std::list<std::pair<std::string, Dist> >::const_iterator i = m_distributions.begin();
             i != m_distributions.end(); ++i)
          {
             double x = ceil(m_min_x);
-            double interval = 1;
+            double discrete_interval = 1;
             std::map<double, double> data;
             while(x <= m_max_x)
             {
@@ -300,7 +300,7 @@ public:
                data[x + x_off + 0.00001] = p;
                data[x + x_off + x_width] = p;
                data[x + x_off + x_width + 0.00001] = 0;
-               x += interval;
+               x += discrete_interval;
             }
             x_off += x_width;
             svg_2d_plot_series& s = plot.plot(data, i->first);
@@ -312,9 +312,9 @@ public:
             ++color_index;
             color_index = color_index % (sizeof(colors)/sizeof(colors[0]));
          }
-      }
+      } // descrete
       plot.write(file);
-   }
+   } // void plot(const std::string& title, const std::string& file)
 
 private:
    bool m_pdf;
@@ -326,7 +326,7 @@ int main()
 {
   try
   {
-
+   std::cout << "Distribution Graphs" << std::endl;
    distribution_plotter<boost::math::gamma_distribution<> >
       gamma_plotter;
    gamma_plotter.add(boost::math::gamma_distribution<>(0.75), "shape = 0.75");
@@ -718,5 +718,4 @@ int main()
    hyperexponential_plotter3.plot("Hyperexponential Distribution PDF (Different Number of Phases, Same Mean)", "hyperexponential_pdf_samemean.svg");
    */
 
-
 } // int main()

doc/graphs/lambert_w_graph_big_W.svg

@@ -0,0 +1,62 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<svg width="400" height ="300" version="1.1"
+xmlns:svg ="http://www.w3.org/2000/svg"
+xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
+xmlns:cc="http://web.resource.org/cc/"
+xmlns:dc="http://purl.org/dc/elements/1.1/"
+xmlns ="http://www.w3.org/2000/svg"
+>
+<!-- SVG plot written using Boost.Plot program (Creator Jacob Voytko) --> 
+<!-- Use, modification and distribution of Boost.Plot subject to the --> 
+<!-- Boost Software License, Version 1.0.--> 
+<!-- (See accompanying file LICENSE_1_0.txt --> 
+<!-- or copy at http://www.boost.org/LICENSE_1_0.txt) --> 
+
+<!-- SVG Plot Copyright Paul A. Bristow 2018 --> 
+<meta name="copyright" content="Paul A. Bristow" />
+<meta name="date" content="2018" />
+<clipPath id="plot_window"><rect x="76.8" y="54" width="296.2" height="186"/></clipPath>
+<g id="imageBackground" stroke="rgb(255,255,0)" fill="rgb(255,255,255)" stroke-width="2"><rect x="0" y="0" width="400" height="300"/></g>
+<g id="plotBackground" stroke="rgb(119,136,153)" fill="rgb(255,255,255)" stroke-width="2"><rect x="75.8" y="53" width="298.2" height="188"/></g>
+<g id="yMinorGrid" stroke="rgb(200,220,255)" stroke-width="0.5"></g>
+<g id="yMajorGrid" stroke="rgb(200,220,255)" stroke-width="1"><path d="M77.8,220.1 L372,220.1 M77.8,199.2 L372,199.2 M77.8,178.3 L372,178.3 M77.8,157.4 L372,157.4 M77.8,136.6 L372,136.6 M77.8,115.7 L372,115.7 M77.8,94.78 L372,94.78 M77.8,73.89 L372,73.89 M77.8,53 L372,53 M77.8,220.1 L372,220.1 M77.8,241 L372,241 " fill="none"/></g>
+<g id="xMinorGrid" stroke="rgb(200,220,255)" stroke-width="0.5"></g>
+<g id="xMajorGrid" stroke="rgb(200,220,255)" stroke-width="1"><path d="M75.83,53 L75.83,241 M135.5,53 L135.5,241 M195.1,53 L195.1,241 M254.7,53 L254.7,241 M314.4,53 L314.4,241 M374,53 L374,241 M75.83,53 L75.83,241 " fill="none"/></g>
+<g id="yAxis" stroke="rgb(0,0,0)" stroke-width="1"><line x1="75.83" y1="53" x2="75.83" y2="246"/><line x1="75.8" y1="53" x2="75.8" y2="241"/></g>
+<g id="xAxis" stroke="rgb(0,0,0)" stroke-width="1"><line x1="70.8" y1="220.1" x2="374" y2="220.1"/><line x1="70.8" y1="241" x2="374" y2="241"/></g>
+<g id="yMinorTicks" stroke="rgb(0,0,0)" stroke-width="1"><path d="M73.8,215.9 L75.8,215.9 M73.8,211.8 L75.8,211.8 M73.8,207.6 L75.8,207.6 M73.8,203.4 L75.8,203.4 M73.8,195 L75.8,195 M73.8,190.9 L75.8,190.9 M73.8,186.7 L75.8,186.7 M73.8,182.5 L75.8,182.5 M73.8,174.2 L75.8,174.2 M73.8,170 L75.8,170 M73.8,165.8 L75.8,165.8 M73.8,161.6 L75.8,161.6 M73.8,153.3 L75.8,153.3 M73.8,149.1 L75.8,149.1 M73.8,144.9 L75.8,144.9 M73.8,140.7 L75.8,140.7 M73.8,132.4 L75.8,132.4 M73.8,128.2 L75.8,128.2 M73.8,124 L75.8,124 M73.8,119.8 L75.8,119.8 M73.8,111.5 L75.8,111.5 M73.8,107.3 L75.8,107.3 M73.8,103.1 L75.8,103.1 M73.8,98.96 L75.8,98.96 M73.8,90.6 L75.8,90.6 M73.8,86.42 L75.8,86.42 M73.8,82.24 L75.8,82.24 M73.8,78.07 L75.8,78.07 M73.8,69.71 L75.8,69.71 M73.8,65.53 L75.8,65.53 M73.8,61.36 L75.8,61.36 M73.8,57.18 L75.8,57.18 M73.8,220.1 L75.8,220.1 M73.8,224.3 L75.8,224.3 M73.8,228.5 L75.8,228.5 M73.8,232.6 L75.8,232.6 M73.8,236.8 L75.8,236.8 M73.8,241 L75.8,241 " fill="none"/></g>
+<g id="xMinorTicks" stroke="rgb(0,0,0)" stroke-width="1"><path d="M87.76,241 L87.76,243 M99.68,241 L99.68,243 M111.6,241 L111.6,243 M123.5,241 L123.5,243 M147.4,241 L147.4,243 M159.3,241 L159.3,243 M171.2,241 L171.2,243 M183.2,241 L183.2,243 M207,241 L207,243 M219,241 L219,243 M230.9,241 L230.9,243 M242.8,241 L242.8,243 M266.7,241 L266.7,243 M278.6,241 L278.6,243 M290.5,241 L290.5,243 M302.4,241 L302.4,243 M326.3,241 L326.3,243 M338.2,241 L338.2,243 M350.1,241 L350.1,243 M362.1,241 L362.1,243 " fill="none"/></g>
+<g id="yMajorTicks" stroke="rgb(0,0,0)" stroke-width="2"><path d="M70.8,220.1 L75.8,220.1 M70.8,199.2 L75.8,199.2 M70.8,178.3 L75.8,178.3 M70.8,157.4 L75.8,157.4 M70.8,136.6 L75.8,136.6 M70.8,115.7 L75.8,115.7 M70.8,94.78 L75.8,94.78 M70.8,73.89 L75.8,73.89 M70.8,53 L75.8,53 M70.8,220.1 L75.8,220.1 M70.8,241 L75.8,241 " fill="none"/></g>
+<g id="xMajorTicks" stroke="rgb(0,0,0)" stroke-width="2"><path d="M75.83,241 L75.83,246 M135.5,241 L135.5,246 M195.1,241 L195.1,246 M254.7,241 L254.7,246 M314.4,241 L314.4,246 M374,241 L374,246 M75.83,241 L75.83,246 " fill="none"/></g>
+<g id="xTicksValues">
+<text x="75.83" y="261.6" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">0</text>
+<text x="135.5" y="261.6" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">2e3</text>
+<text x="195.1" y="261.6" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">4e3</text>
+<text x="254.7" y="261.6" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">6e3</text>
+<text x="314.4" y="261.6" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">8e3</text>
+<text x="374" y="261.6" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">1e4</text>
+<text x="75.83" y="261.6" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">0</text></g>
+<g id="yTicksValues">
+<text x="64.8" y="222.5" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">0</text>
+<text x="64.8" y="201.6" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">1</text>
+<text x="64.8" y="180.7" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">2</text>
+<text x="64.8" y="159.8" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">3</text>
+<text x="64.8" y="139" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">4</text>
+<text x="64.8" y="118.1" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">5</text>
+<text x="64.8" y="97.18" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">6</text>
+<text x="64.8" y="76.29" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">7</text>
+<text x="64.8" y="55.4" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">8</text>
+<text x="64.8" y="222.5" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">0</text>
+<text x="64.8" y="243.4" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">-1</text></g>
+<g id="yLabel">
+<text x="42.9" y="147" text-anchor="middle" transform = "rotate(-90 42.9 147 )" font-size="14" font-family="Lucida Sans Unicode">W</text></g>
+<g id="xLabel">
+<text x="224.9" y="276.7" text-anchor="middle" font-size="14" font-family="Lucida Sans Unicode">z</text></g>
+<g id="plotLines" stroke-width="2"><g clip-path="url(#plot_window)" stroke="rgb(255,0,0)" stroke-width="1"><path d="M75.82,241 L75.82,241 L77.31,160.5 L78.8,149.4 L80.29,142.8 L81.78,138.1 L83.27,134.3 L84.76,131.2 L86.25,128.6 L87.75,126.3 L89.24,124.3 L90.73,122.5 L92.22,120.9 L93.71,119.4 L95.2,118 L96.69,116.7 L98.18,115.5 L99.67,114.4 L101.2,113.3 L102.7,112.3 L104.1,111.4 L105.6,110.5 L107.1,109.6 L108.6,108.8 L110.1,108 L111.6,107.3 L113.1,106.5 L114.6,105.8 L116.1,105.2 L117.6,104.5 L119.1,103.9 L120.5,103.3 L122,102.7 L123.5,102.2 L125,101.6 L126.5,101.1 L128,100.6 L129.5,100.1 L131,99.58 L132.5,99.11 L134,98.64 L135.5,98.19 L136.9,97.75 L138.4,97.32 L139.9,96.9 L141.4,96.49 L142.9,96.09 L144.4,95.7 L145.9,95.31 L147.4,94.93 L148.9,94.56 L150.4,94.2 L151.9,93.85 L153.3,93.5 L154.8,93.16 L156.3,92.82 L157.8,92.49 L159.3,92.17 L160.8,91.85 L162.3,91.54 L163.8,91.23 L165.3,90.93 L166.8,90.63 L168.3,90.34 L169.7,90.05 L171.2,89.77 L172.7,89.49 L174.2,89.21 L175.7,88.94 L177.2,88.68 L178.7,88.41 L180.2,88.15 L181.7,87.9 L183.2,87.64 L184.7,87.4 L186.1,87.15 L187.6,86.91 L189.1,86.67 L190.6,86.43 L192.1,86.2 L193.6,85.97 L195.1,85.74 L196.6,85.52 L198.1,85.3 L199.6,85.08 L201.1,84.86 L202.5,84.64 L204,84.43 L205.5,84.22 L207,84.02 L208.5,83.81 L210,83.61 L211.5,83.41 L213,83.21 L214.5,83.02 L216,82.82 L217.4,82.63 L218.9,82.44 L220.4,82.25 L221.9,82.07 L223.4,81.88 L224.9,81.7 L226.4,81.52 L227.9,81.34 L229.4,81.16 L230.9,80.99 L232.4,80.81 L233.8,80.64 L235.3,80.47 L236.8,80.3 L238.3,80.13 L239.8,79.97 L241.3,79.8 L242.8,79.64 L244.3,79.48 L245.8,79.32 L247.3,79.16 L248.8,79 L250.2,78.85 L251.7,78.69 L253.2,78.54 L254.7,78.38 L256.2,78.23 L257.7,78.08 L259.2,77.94 L260.7,77.79 L262.2,77.64 L263.7,77.5 L265.2,77.35 L266.6,77.21 L268.1,77.07 L269.6,76.93 L271.1,76.79 L272.6,76.65 L274.1,76.51 L275.6,76.37 L277.1,76.24 L278.6,76.1 L280.1,75.97 L281.6,75.84 L283,75.71 L284.5,75.57 L286,75.44 L287.5,75.32 L289,75.19 L290.5,75.06 L292,74.93 L293.5,74.81 L295,74.68 L296.5,74.56 L298,74.44 L299.4,74.31 L300.9,74.19 L302.4,74.07 L303.9,73.95 L305.4,73.83 L306.9,73.72 L308.4,73.6 L309.9,73.48 L311.4,73.36 L312.9,73.25 L314.4,73.13 L315.8,73.02 L317.3,72.91 L318.8,72.8 L320.3,72.68 L321.8,72.57 L323.3,72.46 L324.8,72.35 L326.3,72.24 L327.8,72.13 L329.3,72.03 L330.8,71.92 L332.2,71.81 L333.7,71.71 L335.2,71.6 L336.7,71.49 L338.2,71.39 L339.7,71.29 L341.2,71.18 L342.7,71.08 L344.2,70.98 L345.7,70.88 L347.2,70.78 L348.6,70.68 L350.1,70.58 L351.6,70.48 L353.1,70.38 L354.6,70.28 L356.1,70.18 L357.6,70.08 L359.1,69.99 L360.6,69.89 L362.1,69.8 L363.6,69.7 L365,69.61 L366.5,69.51 L368,69.42 L369.5,69.32 L371,69.23 L372.5,69.14 L374,69.05 " fill="none"/></g>
+</g>
+<g id="plotPoints" clip-path="url(#plot_window)"></g>
+<g id="title">
+<text x="200" y="36" text-anchor="middle" font-size="18" font-family="Lucida Sans Unicode">Lambert W0 function for larger z.</text></g>
+<g id="plotXValues"></g>
+<g id="plotYValues"></g>
+</svg>

doc/graphs/lambert_wm1_errors_graph.svg

@@ -0,0 +1,55 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<svg width="400" height ="200" version="1.1"
+xmlns:svg ="http://www.w3.org/2000/svg"
+xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
+xmlns:cc="http://web.resource.org/cc/"
+xmlns:dc="http://purl.org/dc/elements/1.1/"
+xmlns ="http://www.w3.org/2000/svg"
+>
+<!-- SVG plot written using Boost.Plot program (Creator Jacob Voytko) --> 
+<!-- Use, modification and distribution of Boost.Plot subject to the --> 
+<!-- Boost Software License, Version 1.0.--> 
+<!-- (See accompanying file LICENSE_1_0.txt --> 
+<!-- or copy at http://www.boost.org/LICENSE_1_0.txt) --> 
+
+<!-- SVG Plot Copyright Paul A. Bristow 2018 --> 
+<meta name="copyright" content="Paul A. Bristow" />
+<meta name="date" content="2018" />
+<clipPath id="plot_window"><rect x="66.8" y="36.5" width="306" height="114"/></clipPath>
+<g id="imageBackground" stroke="rgb(255,255,0)" fill="rgb(255,255,255)" stroke-width="2"><rect x="0" y="0" width="400" height="200"/></g>
+<g id="plotBackground" stroke="rgb(119,136,153)" fill="rgb(255,255,255)" stroke-width="2"><rect x="65.8" y="35.5" width="308" height="116"/></g>
+<g id="yMinorGrid" stroke="rgb(200,220,255)" stroke-width="0.5"></g>
+<g id="yMajorGrid" stroke="rgb(200,220,255)" stroke-width="1"><path d="M67.8,93.3 L372,93.3 M67.8,64.4 L372,64.4 M67.8,35.5 L372,35.5 M67.8,93.3 L372,93.3 M67.8,122 L372,122 M67.8,151 L372,151 " fill="none"/></g>
+<g id="xMinorGrid" stroke="rgb(200,220,255)" stroke-width="0.5"></g>
+<g id="xMajorGrid" stroke="rgb(200,220,255)" stroke-width="1"><path d="M374,35.5 L374,151 M374,35.5 L374,151 M295,35.5 L295,151 M216,35.5 L216,151 M137,35.5 L137,151 " fill="none"/></g>
+<g id="yAxis" stroke="rgb(0,0,0)" stroke-width="1"><line x1="374" y1="35.5" x2="374" y2="156"/><line x1="65.8" y1="35.5" x2="65.8" y2="151"/></g>
+<g id="xAxis" stroke="rgb(0,0,0)" stroke-width="1"><line x1="60.8" y1="93.3" x2="374" y2="93.3"/><line x1="60.8" y1="151" x2="374" y2="151"/></g>
+<g id="yMinorTicks" stroke="rgb(0,0,0)" stroke-width="1"><path d="M63.8,87.5 L65.8,87.5 M63.8,81.7 L65.8,81.7 M63.8,75.9 L65.8,75.9 M63.8,70.2 L65.8,70.2 M63.8,58.6 L65.8,58.6 M63.8,52.8 L65.8,52.8 M63.8,47.1 L65.8,47.1 M63.8,41.3 L65.8,41.3 M63.8,93.3 L65.8,93.3 M63.8,99 L65.8,99 M63.8,105 L65.8,105 M63.8,111 L65.8,111 M63.8,116 L65.8,116 M63.8,122 L65.8,122 M63.8,128 L65.8,128 M63.8,134 L65.8,134 M63.8,139 L65.8,139 M63.8,145 L65.8,145 M63.8,151 L65.8,151 M63.8,151 L65.8,151 " fill="none"/></g>
+<g id="xMinorTicks" stroke="rgb(0,0,0)" stroke-width="1"><path d="M358,151 L358,153 M342,151 L342,153 M327,151 L327,153 M311,151 L311,153 M279,151 L279,153 M263,151 L263,153 M248,151 L248,153 M232,151 L232,153 M200,151 L200,153 M184,151 L184,153 M169,151 L169,153 M153,151 L153,153 M121,151 L121,153 M105,151 L105,153 M89.5,151 L89.5,153 M73.7,151 L73.7,153 " fill="none"/></g>
+<g id="yMajorTicks" stroke="rgb(0,0,0)" stroke-width="2"><path d="M60.8,93.3 L65.8,93.3 M60.8,64.4 L65.8,64.4 M60.8,35.5 L65.8,35.5 M60.8,93.3 L65.8,93.3 M60.8,122 L65.8,122 M60.8,151 L65.8,151 " fill="none"/></g>
+<g id="xMajorTicks" stroke="rgb(0,0,0)" stroke-width="2"><path d="M374,151 L374,156 M374,151 L374,156 M295,151 L295,156 M216,151 L216,156 M137,151 L137,156 " fill="none"/></g>
+<g id="xTicksValues">
+<text x="374" y="172" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">0</text>
+<text x="374" y="172" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">0</text>
+<text x="295" y="172" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">-0.1</text>
+<text x="216" y="172" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">-0.2</text>
+<text x="137" y="172" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">-0.3</text></g>
+<g id="yTicksValues">
+<text x="54.8" y="95.7" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">0</text>
+<text x="54.8" y="66.8" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">2</text>
+<text x="54.8" y="37.9" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">4</text>
+<text x="54.8" y="95.7" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">0</text>
+<text x="54.8" y="125" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">-2</text>
+<text x="54.8" y="153" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">-4</text></g>
+<g id="yLabel">
+<text x="35.9" y="93.3" text-anchor="middle" transform = "rotate(-90 35.9 93.3 )" font-size="9" font-family="Lucida Sans Unicode">W-1 difference (bits)</text></g>
+<g id="xLabel">
+<text x="220" y="180" text-anchor="middle" font-size="9" font-family="Lucida Sans Unicode">z</text></g>
+<g id="plotLines" stroke-width="2"><g clip-path="url(#plot_window)" stroke="rgb(0,0,139)" stroke-width="0.2"><path d="M84.1,78.8 L84.1,78.8 L84.9,64.4 L85.6,78.8 L86.4,49.9 L87.2,64.4 L88,93.3 L88.8,49.9 L89.6,64.4 L90.4,78.8 L91.2,49.9 L92,64.4 L92.8,78.8 L93.5,93.3 L94.3,78.8 L95.1,78.8 L95.9,93.3 L96.7,108 L99.1,64.4 L101,49.9 L101,49.9 L103,93.3 L104,35.5 L105,49.9 L106,35.5 L107,35.5 L109,78.8 L109,64.4 L110,35.5 L111,93.3 L112,93.3 L113,108 L113,108 L114,122 L115,108 L116,35.5 L116,78.8 L118,64.4 L119,49.9 L120,93.3 L120,122 L121,93.3 L122,93.3 L123,78.8 L124,93.3 L124,108 L125,93.3 L126,108 L127,122 L128,49.9 L128,78.8 L129,108 L130,64.4 L131,108 L131,122 L132,64.4 L133,78.8 L134,78.8 L135,78.8 L135,78.8 L136,78.8 L137,64.4 L138,93.3 L139,93.3 L139,93.3 L140,64.4 L142,78.8 L143,49.9 L144,78.8 L145,64.4 L146,49.9 L146,49.9 L147,78.8 L148,78.8 L149,64.4 L150,64.4 L150,93.3 L151,93.3 L152,64.4 L153,78.8 L154,35.5 L154,64.4 L155,49.9 L156,64.4 L157,78.8 L158,49.9 L159,78.8 L160,93.3 L161,93.3 L161,93.3 L162,93.3 L163,93.3 L164,108 L165,108 L165,93.3 L166,108 L167,78.8 L168,108 L169,78.8 L169,93.3 L170,108 L171,108 L172,93.3 L173,93.3 L173,78.8 L174,93.3 L175,108 L176,93.3 L177,93.3 L177,78.8 L178,108 L179,93.3 L180,93.3 L180,93.3 L181,108 L182,93.3 L183,108 L184,93.3 L184,93.3 L185,78.8 L186,108 L187,93.3 L188,93.3 L188,78.8 L189,78.8 L190,93.3 L191,108 L192,78.8 L192,93.3 L193,78.8 L194,108 L195,93.3 L195,93.3 L196,93.3 L197,78.8 L198,93.3 L199,108 L199,93.3 L200,93.3 L201,78.8 L202,108 L203,93.3 L203,108 L204,93.3 L205,93.3 L206,93.3 L207,78.8 L207,93.3 L208,78.8 L209,78.8 L210,93.3 L210,93.3 L211,78.8 L212,93.3 L213,93.3 L214,93.3 L214,78.8 L215,78.8 L216,93.3 L217,78.8 L218,93.3 L218,78.8 L219,78.8 L220,64.4 L221,49.9 L222,93.3 L222,93.3 L223,78.8 L224,93.3 L225,93.3 L225,108 L226,93.3 L227,78.8 L228,93.3 L229,93.3 L229,93.3 L230,93.3 L231,93.3 L232,93.3 L233,78.8 L233,93.3 L234,93.3 L235,108 L236,78.8 L237,78.8 L237,93.3 L238,93.3 L239,93.3 L240,93.3 L240,78.8 L241,78.8 L242,78.8 L243,78.8 L244,78.8 L244,64.4 L245,78.8 L246,78.8 L247,93.3 L248,78.8 L248,93.3 L249,108 L250,93.3 L251,108 L252,93.3 L252,108 L253,108 L254,78.8 L255,78.8 L256,108 L256,78.8 L257,93.3 L258,108 L259,93.3 L259,93.3 L260,93.3 L261,93.3 L262,78.8 L263,78.8 L263,122 L264,93.3 L265,93.3 L266,78.8 L267,108 L267,93.3 L268,122 L269,93.3 L270,78.8 L271,93.3 L271,93.3 L272,93.3 L273,108 L274,122 L274,108 L275,93.3 L276,93.3 L277,93.3 L278,93.3 L278,108 L279,78.8 L280,93.3 L281,93.3 L282,93.3 L282,93.3 L283,93.3 L284,78.8 L285,108 L286,93.3 L286,93.3 L287,78.8 L288,78.8 L289,93.3 L289,108 L290,78.8 L291,93.3 L292,93.3 L293,93.3 L293,93.3 L294,93.3 L295,93.3 L296,93.3 L297,108 L297,122 L298,108 L299,108 L300,93.3 L301,108 L301,78.8 L302,78.8 L303,93.3 L304,93.3 L304,93.3 L305,78.8 L306,93.3 L307,93.3 L308,93.3 L308,93.3 L309,93.3 L310,93.3 L311,93.3 L312,78.8 L312,93.3 L313,93.3 L314,93.3 L315,93.3 L316,78.8 L316,93.3 L317,93.3 L318,93.3 L319,93.3 L320,93.3 L320,93.3 L321,108 L322,93.3 L323,93.3 L323,108 L324,78.8 L325,93.3 L326,93.3 L327,78.8 L327,93.3 L328,108 L329,78.8 L330,93.3 L331,78.8 L331,93.3 L332,108 L333,93.3 L334,93.3 L335,93.3 L335,93.3 L336,108 L337,78.8 L338,93.3 L338,78.8 L339,93.3 L340,93.3 L341,93.3 L342,78.8 L342,78.8 L343,93.3 L344,93.3 L345,93.3 L346,93.3 L346,93.3 L347,93.3 L348,93.3 L349,93.3 L350,93.3 L350,93.3 L351,93.3 L352,93.3 L353,93.3 L353,93.3 L354,93.3 L355,93.3 L356,108 L357,93.3 L357,78.8 L358,93.3 L359,93.3 L360,93.3 L361,93.3 L361,93.3 L362,93.3 L363,93.3 L364,93.3 L365,93.3 L365,93.3 L366,78.8 L367,93.3 L368,93.3 L368,93.3 L369,93.3 L370,93.3 L371,93.3 L372,93.3 L372,108 L373,93.3 " fill="none"/></g>
+</g>
+<g id="plotPoints" clip-path="url(#plot_window)"></g>
+<g id="title">
+<text x="200" y="22" text-anchor="middle" font-size="11" font-family="Lucida Sans Unicode">Lambert W-1 function differences from 'best' for double.</text></g>
+<g id="plotXValues"></g>
+<g id="plotYValues"></g>
+</svg>

doc/graphs/lambert_wm1_prime_graph.svg

@@ -0,0 +1,54 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<svg width="400" height ="300" version="1.1"
+xmlns:svg ="http://www.w3.org/2000/svg"
+xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
+xmlns:cc="http://web.resource.org/cc/"
+xmlns:dc="http://purl.org/dc/elements/1.1/"
+xmlns ="http://www.w3.org/2000/svg"
+>
+<!-- SVG plot written using Boost.Plot program (Creator Jacob Voytko) --> 
+<!-- Use, modification and distribution of Boost.Plot subject to the --> 
+<!-- Boost Software License, Version 1.0.--> 
+<!-- (See accompanying file LICENSE_1_0.txt --> 
+<!-- or copy at http://www.boost.org/LICENSE_1_0.txt) --> 
+
+<!-- SVG Plot Copyright Paul A. Bristow 2018 --> 
+<meta name="copyright" content="Paul A. Bristow" />
+<meta name="date" content="2018" />
+<clipPath id="plot_window"><rect x="85.2" y="54" width="287.8" height="204"/></clipPath>
+<g id="imageBackground" stroke="rgb(255,255,0)" fill="rgb(255,255,255)" stroke-width="2"><rect x="0" y="0" width="400" height="300"/></g>
+<g id="plotBackground" stroke="rgb(119,136,153)" fill="rgb(255,255,255)" stroke-width="2"><rect x="84.2" y="53" width="289.8" height="206"/></g>
+<g id="yMinorGrid" stroke="rgb(200,220,255)" stroke-width="0.5"></g>
+<g id="yMajorGrid" stroke="rgb(200,220,255)" stroke-width="1"><path d="M86.2,53 L372,53 M86.2,121.7 L372,121.7 M86.2,190.3 L372,190.3 M86.2,259 L372,259 " fill="none"/></g>
+<g id="xMinorGrid" stroke="rgb(200,220,255)" stroke-width="0.5"></g>
+<g id="xMajorGrid" stroke="rgb(200,220,255)" stroke-width="1"><path d="M366.9,53 L366.9,259 M366.9,53 L366.9,259 M296.2,53 L296.2,259 M225.6,53 L225.6,259 M154.9,53 L154.9,259 M84.2,53 L84.2,259 " fill="none"/></g>
+<g id="yAxis" stroke="rgb(0,0,0)" stroke-width="1"><line x1="366.9" y1="53" x2="366.9" y2="259"/><line x1="84.2" y1="53" x2="84.2" y2="259"/></g>
+<g id="xAxis" stroke="rgb(0,0,0)" stroke-width="1"><line x1="84.2" y1="53" x2="374" y2="53"/></g>
+<g id="yMinorTicks" stroke="rgb(0,0,0)" stroke-width="1"><path d="M82.2,66.73 L84.2,66.73 M82.2,80.47 L84.2,80.47 M82.2,94.2 L84.2,94.2 M82.2,107.9 L84.2,107.9 M82.2,121.7 L84.2,121.7 M82.2,135.4 L84.2,135.4 M82.2,149.1 L84.2,149.1 M82.2,162.9 L84.2,162.9 M82.2,176.6 L84.2,176.6 M82.2,190.3 L84.2,190.3 M82.2,204.1 L84.2,204.1 M82.2,217.8 L84.2,217.8 M82.2,231.5 L84.2,231.5 M82.2,245.3 L84.2,245.3 M82.2,259 L84.2,259 " fill="none"/></g>
+<g id="xMinorTicks" stroke="rgb(0,0,0)" stroke-width="1"><path d="M352.8,53 L352.8,55 M338.7,53 L338.7,55 M324.5,53 L324.5,55 M310.4,53 L310.4,55 M282.1,53 L282.1,55 M268,53 L268,55 M253.8,53 L253.8,55 M239.7,53 L239.7,55 M211.4,53 L211.4,55 M197.3,53 L197.3,55 M183.2,53 L183.2,55 M169,53 L169,55 M140.7,53 L140.7,55 M126.6,53 L126.6,55 M112.5,53 L112.5,55 M98.34,53 L98.34,55 " fill="none"/></g>
+<g id="yMajorTicks" stroke="rgb(0,0,0)" stroke-width="2"><path d="M79.2,53 L84.2,53 M79.2,121.7 L84.2,121.7 M79.2,190.3 L84.2,190.3 M79.2,259 L84.2,259 " fill="none"/></g>
+<g id="xMajorTicks" stroke="rgb(0,0,0)" stroke-width="2"><path d="M366.9,53 L366.9,58 M366.9,53 L366.9,58 M296.2,53 L296.2,58 M225.6,53 L225.6,58 M154.9,53 L154.9,58 M84.2,53 L84.2,58 " fill="none"/></g>
+<g id="xTicksValues">
+<text x="366.9" y="73.6" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">0</text>
+<text x="366.9" y="73.6" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">0</text>
+<text x="296.2" y="73.6" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">-0.1</text>
+<text x="225.6" y="73.6" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">-0.2</text>
+<text x="154.9" y="73.6" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">-0.3</text>
+<text x="84.2" y="73.6" text-anchor="middle" font-size="12" font-family="Lucida Sans Unicode">-0.4</text></g>
+<g id="yTicksValues">
+<text x="73.2" y="55.4" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">-5</text>
+<text x="73.2" y="124.1" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">-10</text>
+<text x="73.2" y="192.7" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">-15</text>
+<text x="73.2" y="261.4" text-anchor="end" font-size="12" font-family="Lucida Sans Unicode">-20</text></g>
+<g id="yLabel">
+<text x="42.9" y="156" text-anchor="middle" transform = "rotate(-90 42.9 156 )" font-size="14" font-family="Lucida Sans Unicode">W-1'</text></g>
+<g id="xLabel">
+<text x="229.1" y="282.8" text-anchor="middle" font-size="14" font-family="Lucida Sans Unicode">z</text></g>
+<g id="plotLines" stroke-width="2"><g clip-path="url(#plot_window)" stroke="rgb(0,0,255)" stroke-width="1"><path d="M110.5,237.1 L110.5,237.1 L111.2,218.1 L111.9,203.4 L112.6,191.5 L113.3,181.7 L114,173.5 L114.7,166.4 L115.4,160.2 L116.1,154.8 L116.9,150 L117.6,145.7 L118.3,141.9 L119,138.3 L119.7,135.2 L120.4,132.2 L121.1,129.5 L121.8,127.1 L122.5,124.8 L123.2,122.7 L123.9,120.7 L124.6,118.8 L125.3,117.1 L126,115.5 L126.8,113.9 L127.5,112.5 L128.2,111.2 L128.9,109.9 L129.6,108.7 L130.3,107.5 L131,106.5 L131.7,105.4 L132.4,104.5 L133.1,103.5 L133.8,102.7 L134.5,101.8 L135.2,101 L135.9,100.3 L136.6,99.54 L137.4,98.85 L138.1,98.19 L138.8,97.56 L139.5,96.95 L140.2,96.38 L140.9,95.82 L141.6,95.3 L142.3,94.79 L143,94.31 L143.7,93.85 L144.4,93.4 L145.1,92.98 L145.8,92.57 L146.5,92.18 L147.2,91.81 L148,91.46 L148.7,91.12 L149.4,90.79 L150.1,90.48 L150.8,90.18 L151.5,89.9 L152.2,89.63 L152.9,89.37 L153.6,89.12 L154.3,88.88 L155,88.66 L155.7,88.44 L156.4,88.24 L157.1,88.04 L157.9,87.86 L158.6,87.68 L159.3,87.52 L160,87.36 L160.7,87.22 L161.4,87.08 L162.1,86.95 L162.8,86.82 L163.5,86.71 L164.2,86.6 L164.9,86.5 L165.6,86.41 L166.3,86.32 L167,86.25 L167.7,86.17 L168.5,86.11 L169.2,86.05 L169.9,86 L170.6,85.96 L171.3,85.92 L172,85.88 L172.7,85.86 L173.4,85.84 L174.1,85.82 L174.8,85.81 L175.5,85.81 L176.2,85.81 L176.9,85.82 L177.6,85.83 L178.3,85.85 L179.1,85.87 L179.8,85.9 L180.5,85.94 L181.2,85.98 L181.9,86.02 L182.6,86.07 L183.3,86.12 L184,86.18 L184.7,86.24 L185.4,86.31 L186.1,86.38 L186.8,86.46 L187.5,86.54 L188.2,86.63 L189,86.72 L189.7,86.82 L190.4,86.92 L191.1,87.02 L191.8,87.13 L192.5,87.25 L193.2,87.37 L193.9,87.49 L194.6,87.62 L195.3,87.75 L196,87.88 L196.7,88.02 L197.4,88.17 L198.1,88.32 L198.8,88.47 L199.6,88.63 L200.3,88.79 L201,88.96 L201.7,89.13 L202.4,89.3 L203.1,89.48 L203.8,89.66 L204.5,89.85 L205.2,90.04 L205.9,90.24 L206.6,90.44 L207.3,90.65 L208,90.86 L208.7,91.07 L209.5,91.29 L210.2,91.51 L210.9,91.74 L211.6,91.97 L212.3,92.21 L213,92.45 L213.7,92.69 L214.4,92.94 L215.1,93.2 L215.8,93.45 L216.5,93.72 L217.2,93.99 L217.9,94.26 L218.6,94.54 L219.3,94.82 L220.1,95.11 L220.8,95.4 L221.5,95.69 L222.2,96 L222.9,96.3 L223.6,96.61 L224.3,96.93 L225,97.25 L225.7,97.58 L226.4,97.91 L227.1,98.25 L227.8,98.59 L228.5,98.94 L229.2,99.29 L229.9,99.65 L230.7,100 L231.4,100.4 L232.1,100.8 L232.8,101.1 L233.5,101.5 L234.2,101.9 L234.9,102.3 L235.6,102.7 L236.3,103.1 L237,103.5 L237.7,104 L238.4,104.4 L239.1,104.8 L239.8,105.3 L240.6,105.7 L241.3,106.2 L242,106.6 L242.7,107.1 L243.4,107.6 L244.1,108 L244.8,108.5 L245.5,109 L246.2,109.5 L246.9,110 L247.6,110.6 L248.3,111.1 L249,111.6 L249.7,112.2 L250.4,112.7 L251.2,113.3 L251.9,113.8 L252.6,114.4 L253.3,115 L254,115.6 L254.7,116.2 L255.4,116.8 L256.1,117.4 L256.8,118 L257.5,118.6 L258.2,119.3 L258.9,119.9 L259.6,120.6 L260.3,121.3 L261,122 L261.8,122.6 L262.5,123.3 L263.2,124.1 L263.9,124.8 L264.6,125.5 L265.3,126.3 L266,127 L266.7,127.8 L267.4,128.6 L268.1,129.4 L268.8,130.2 L269.5,131 L270.2,131.8 L270.9,132.7 L271.7,133.6 L272.4,134.4 L273.1,135.3 L273.8,136.2 L274.5,137.1 L275.2,138.1 L275.9,139 L276.6,140 L277.3,141 L278,142 L278.7,143 L279.4,144 L280.1,145.1 L280.8,146.1 L281.5,147.2 L282.3,148.3 L283,149.4 L283.7,150.6 L284.4,151.8 L285.1,152.9 L285.8,154.2 L286.5,155.4 L287.2,156.6 L287.9,157.9 L288.6,159.2 L289.3,160.5 L290,161.9 L290.7,163.3 L291.4,164.7 L292.1,166.1 L292.9,167.5 L293.6,169 L294.3,170.5 L295,172.1 L295.7,173.7 L296.4,175.3 L297.1,176.9 L297.8,178.6 L298.5,180.3 L299.2,182.1 L299.9,183.8 L300.6,185.7 L301.3,187.5 L302,189.4 L302.8,191.4 L303.5,193.4 L304.2,195.4 L304.9,197.5 L305.6,199.6 L306.3,201.8 L307,204.1 L307.7,206.3 L308.4,208.7 L309.1,211.1 L309.8,213.5 L310.5,216.1 L311.2,218.7 L311.9,221.3 L312.6,224 L313.4,226.8 L314.1,229.7 L314.8,232.7 L315.5,235.7 L316.2,238.8 L316.9,242 L317.6,245.3 L318.3,248.7 L319,252.2 L319.7,255.8 L320.4,259.5 " fill="none"/></g>
+</g>
+<g id="plotPoints" clip-path="url(#plot_window)"></g>
+<g id="title">
+<text x="200" y="36" text-anchor="middle" font-size="18" font-family="Lucida Sans Unicode">Lambert W-1 prime function.</text></g>
+<g id="plotXValues"></g>
+<g id="plotYValues"></g>
+</svg>

doc/graphs/plot_1d_errors.cpp

@@ -0,0 +1,407 @@
+
+#include <map>
+#include <boost/config.hpp>
+#include <boost/multiprecision/cpp_bin_float.hpp>
+#ifdef BOOST_HAS_FLOAT128
+#include <boost/multiprecision/float128.hpp>
+#endif
+#include <boost/svg_plot/svg_2d_plot.hpp>
+
+template <class Real>
+Real interval_from_range(Real x)
+{
+   BOOST_MATH_STD_USING
+   Real l = floor(log10(x));
+   l = pow(10, l);
+   if (x / l < 2)
+      l /= 10;
+   return l;
+}
+
+
+std::string normalise_filename(std::string name)
+{
+   for(std::string::size_type i = 0; i < name.size(); ++i)
+   {
+      if (!std::isalnum(name[i]))
+         name[i] = '_';
+      else
+         name[i] = std::tolower(name[i]);
+   }
+   return name;
+}
+
+template <class F, class Real>
+void plot_errors_1d(F f, Real start, Real end, unsigned points, const char* function_name, Real max_y_scale = (std::numeric_limits<Real>::max)(), unsigned num_bins = 200)
+{
+   BOOST_MATH_STD_USING
+   std::cout << "Generating points for " << function_name << std::endl;
+   Real pos = start;
+   Real interval = (end - start) / points;
+
+   std::map<Real, Real> points_upper, points_lower;
+
+   Real max_distance(0), min_distance(0), max_error(0), max_error_location(0);
+
+   constexpr unsigned limb_bits = (sizeof(boost::multiprecision::limb_type) * CHAR_BIT);
+   constexpr unsigned mp_digits = (((std::numeric_limits<Real>::digits * 2) / limb_bits + ((std::numeric_limits<Real>::digits * 2) % limb_bits ? 1 : 0))) * limb_bits;
+
+   typedef boost::multiprecision::number<boost::multiprecision::cpp_bin_float<mp_digits, boost::multiprecision::backends::digit_base_2> > mp_type;
+
+   while (pos <= end)
+   {
+      try
+      {
+         Real found_value = f(pos);
+         Real exact_value = static_cast<Real>(f(mp_type(pos)));
+         Real distance = boost::math::sign(found_value - exact_value) * boost::math::epsilon_difference(found_value, exact_value);
+         Real bin = start + ((end - start) / num_bins) * boost::math::itrunc(num_bins * (pos - start) / (end - start));
+         if (points_lower.find(bin) == points_lower.end())
+            points_lower[bin] = 0;
+         if (points_upper.find(bin) == points_upper.end())
+            points_upper[bin] = 0;
+         if (distance > 0)
+         {
+            if (points_upper[bin] < distance)
+               points_upper[bin] = (std::min)(distance, max_y_scale);
+         }
+         else
+         {
+            if (points_lower[bin] > distance)
+               points_lower[bin] = (std::max)(distance, -max_y_scale);
+         }
+         if (max_distance < distance)
+            max_distance = (std::min)(distance, max_y_scale);
+         if (min_distance > distance)
+            min_distance = (std::max)(distance, -max_y_scale);
+         if (fabs(distance) > max_error)
+         {
+            max_error = fabs(distance);
+            max_error_location = pos;
+         }
+         pos += interval;
+      }
+      catch (const std::exception& e)
+      {
+         std::cout << "Found exception at point " << pos << " : " << e.what() << std::endl;
+         pos += interval;
+      }
+   }
+
+   std::cout << "Max error was " << std::setprecision(3) << max_error << " at location " << std::setprecision(std::numeric_limits<Real>::max_digits10) << max_error_location << std::endl;
+
+   boost::svg::svg_2d_plot plot;
+   Real x_start(start), x_end(end);
+   if (end - start > 3)
+   {
+      x_start = floor(start);
+      x_end = ceil(end);
+   }
+   if (min_distance == 0)
+      min_distance = -1;
+   if (max_distance == 0)
+      max_distance = 1;
+
+
+   plot.title(std::string("Errors in ") + function_name).x_range((double)x_start, (double)x_end).image_x_size(700).legend_border_color(boost::svg::lightgray).plot_border_color(boost::svg::lightgray).background_border_color(boost::svg::lightgray)
+      .y_range((int)floor(min_distance), (int)ceil(max_distance)).x_label("x").y_major_interval((double)interval_from_range(max_distance) * 2).x_major_interval((double)interval_from_range(end - start)).legend_on(true).plot_window_on(true).legend_on(false);
+   plot.plot(points_upper).stroke_color(boost::svg::green).fill_color(boost::svg::green).size(1).line_on(true).area_fill(boost::svg::green);
+   plot.plot(points_lower).stroke_color(boost::svg::green).fill_color(boost::svg::green).size(1).line_on(true).area_fill(boost::svg::green);
+
+   plot.write(normalise_filename(function_name) + ".svg");
+
+}
+
+#include <boost/math/special_functions.hpp>
+
+struct digamma_func
+{
+   template <class T>
+   T operator()(T x)
+   {
+      return boost::math::digamma(x);
+   }
+};
+
+struct tgamma_func
+{
+   template <class T>
+   T operator()(T x)
+   {
+      return boost::math::tgamma(x);
+   }
+};
+
+struct lgamma_func
+{
+   template <class T>
+   T operator()(T x)
+   {
+      return boost::math::lgamma(x);
+   }
+};
+
+struct trigamma_func
+{
+   template <class T>
+   T operator()(T x)
+   {
+      return boost::math::tgamma(x);
+   }
+};
+
+struct erf_func
+{
+   template <class T>
+   T operator()(T x)
+   {
+      return boost::math::erf(x);
+   }
+};
+
+struct erfc_func
+{
+   template <class T>
+   T operator()(T x)
+   {
+      return boost::math::erfc(x);
+   }
+};
+
+struct j0_func
+{
+   template <class T>
+   T operator()(T x)
+   {
+      return boost::math::cyl_bessel_j(0, x);
+   }
+};
+
+struct j1_func
+{
+   template <class T>
+   T operator()(T x)
+   {
+      return boost::math::cyl_bessel_j(1, x);
+   }
+};
+
+struct y0_func
+{
+   template <class T>
+   T operator()(T x)
+   {
+      return boost::math::cyl_neumann(0, x);
+   }
+};
+
+struct y1_func
+{
+   template <class T>
+   T operator()(T x)
+   {
+      return boost::math::cyl_neumann(1, x);
+   }
+};
+
+struct i0_func
+{
+   template <class T>
+   T operator()(T x)
+   {
+      return boost::math::cyl_bessel_i(0, x);
+   }
+};
+
+struct i1_func
+{
+   template <class T>
+   T operator()(T x)
+   {
+      return boost::math::cyl_bessel_i(1, x);
+   }
+};
+
+struct k0_func
+{
+   template <class T>
+   T operator()(T x)
+   {
+      return boost::math::cyl_bessel_k(0, x);
+   }
+};
+
+struct k1_func
+{
+   template <class T>
+   T operator()(T x)
+   {
+      return boost::math::cyl_bessel_k(1, x);
+   }
+};
+
+struct ai_func
+{
+   template <class T>
+   T operator()(T x)
+   {
+      return boost::math::airy_ai(x);
+   }
+};
+
+struct aip_func
+{
+   template <class T>
+   T operator()(T x)
+   {
+      return boost::math::airy_ai_prime(x);
+   }
+};
+
+struct bi_func
+{
+   template <class T>
+   T operator()(T x)
+   {
+      return boost::math::airy_bi(x);
+   }
+};
+
+struct bip_func
+{
+   template <class T>
+   T operator()(T x)
+   {
+      return boost::math::airy_bi_prime(x);
+   }
+};
+
+struct ellint_1_func
+{
+   template <class T>
+   T operator()(T x)
+   {
+      return boost::math::ellint_1(x);
+   }
+};
+
+struct ellint_2_func
+{
+   template <class T>
+   T operator()(T x)
+   {
+      return boost::math::ellint_2(x);
+   }
+};
+
+struct ellint_d_func
+{
+   template <class T>
+   T operator()(T x)
+   {
+      return boost::math::ellint_d(x);
+   }
+};
+
+struct zeta_func
+{
+   template <class T>
+   T operator()(T x)
+   {
+      return boost::math::zeta(x);
+   }
+};
+
+struct ei_func
+{
+   template <class T>
+   T operator()(T x)
+   {
+      return boost::math::expint(x);
+   }
+};
+
+int main()
+{
+   plot_errors_1d(digamma_func(), 1e-200, 10.0, 10000, "digamma, double");
+   plot_errors_1d(tgamma_func(), 1e-200, 150.0, 10000, "tgamma, double");
+   plot_errors_1d(lgamma_func(), 1e-200, 1000.0, 10000, "lgamma, double");
+   plot_errors_1d(trigamma_func(), 1e-200, 10.0, 10000, "trigamma, double");
+   plot_errors_1d(erf_func(), -5.0, 5.0, 10000, "erf, double");
+   plot_errors_1d(erfc_func(), -5.0, 30.0, 10000, "erfc, double");
+   plot_errors_1d(j0_func(), 0.0, 50.0, 10000, "j0, double", 50.0);
+   plot_errors_1d(j1_func(), 0.0, 50.0, 10000, "j1, double", 50.0);
+   plot_errors_1d(y0_func(), 1e-100, 50.0, 10000, "y0, double", 50.0);
+   plot_errors_1d(y1_func(), 1e-100, 50.0, 10000, "y1, double", 50.0);
+   plot_errors_1d(i0_func(), 0.0, 50.0, 10000, "i0, double");
+   plot_errors_1d(i1_func(), 0.0, 50.0, 10000, "i1, double");
+   plot_errors_1d(k0_func(), 1e-100, 50.0, 10000, "k0, double");
+   plot_errors_1d(k1_func(), 1e-100, 50.0, 10000, "k1, double");
+   plot_errors_1d(ai_func(), -20.0, 20.0, 10000, "Ai, double", 100.0);
+   plot_errors_1d(bi_func(), -20.0, 20.0, 10000, "Bi, double", 100.0);
+   plot_errors_1d(aip_func(), -20.0, 20.0, 10000, "Ai Prime, double", 100.0);
+   plot_errors_1d(bip_func(), -20.0, 20.0, 10000, "Bi Prime, double", 100.0);
+
+   plot_errors_1d(ellint_1_func(), -1.0, 1.0, 10000, "Elliptic Integral K, double");
+   plot_errors_1d(ellint_2_func(), -1.0, 1.0, 10000, "Elliptic Integral E, double");
+   plot_errors_1d(ellint_d_func(), -1.0, 1.0, 10000, "Elliptic Integral D, double");
+
+   plot_errors_1d(zeta_func(), -20.0, 20.0, 10000, "Zeta, double");
+   plot_errors_1d(ei_func(), -20.0, 20.0, 10000, "Exponential Integral Ei, double");
+
+#if LDBL_MANT_DIG == 64
+   plot_errors_1d(digamma_func(), 1e-200L, 10.0L, 10000, "digamma, 80-bit long double");
+   plot_errors_1d(tgamma_func(), 1e-200L, 150.0L, 10000, "tgamma, 80-bit long double");
+   plot_errors_1d(lgamma_func(), 1e-200L, 1000.0L, 10000, "lgamma, 80-bit long double");
+   plot_errors_1d(trigamma_func(), 1e-200L, 10.0L, 10000, "trigamma, 80-bit long double");
+   plot_errors_1d(erf_func(), -5.0L, 5.0L, 10000, "erf, 80-bit long double");
+   plot_errors_1d(erfc_func(), -5.0L, 120.0L, 10000, "erfc, 80-bit long double");
+   plot_errors_1d(j0_func(), 0.0L, 50.0L, 10000, "j0, 80 bit long double", 50.0L);
+   plot_errors_1d(j1_func(), 0.0L, 50.0L, 10000, "j1, 80 bit long double", 50.0L);
+   plot_errors_1d(y0_func(), 1e-100L, 50.0L, 10000, "y0, 80 bit long double", 50.0L);
+   plot_errors_1d(y1_func(), 1e-100L, 50.0L, 10000, "y1, 80 bit long double", 50.0L);
+   plot_errors_1d(i0_func(), 0.0L, 50.0L, 10000, "i0, 80 bit long double");
+   plot_errors_1d(i1_func(), 0.0L, 50.0L, 10000, "i1, 80 bit long double");
+   plot_errors_1d(k0_func(), 1e-100L, 50.0L, 10000, "k0, 80 bit long double");
+   plot_errors_1d(k1_func(), 1e-100L, 50.0L, 10000, "k1, 80 bit long double");
+   plot_errors_1d(ai_func(), -20.0L, 20.0L, 10000, "Ai, 80 bit long double", 100.0L);
+   plot_errors_1d(bi_func(), -20.0L, 20.0L, 10000, "Bi, 80 bit long double", 100.0L);
+   plot_errors_1d(aip_func(), -20.0L, 20.0L, 10000, "Ai Prime, 80 bit long double", 100.0L);
+   plot_errors_1d(bip_func(), -20.0L, 20.0L, 10000, "Bi Prime, 80 bit long double", 100.0L);
+
+   plot_errors_1d(ellint_1_func(), -1.0L, 1.0L, 10000, "Elliptic Integral K, 80 bit long double");
+   plot_errors_1d(ellint_2_func(), -1.0L, 1.0L, 10000, "Elliptic Integral E, 80 bit long double");
+   plot_errors_1d(ellint_d_func(), -1.0L, 1.0L, 10000, "Elliptic Integral D, 80 bit long double");
+
+   plot_errors_1d(zeta_func(), -20.0L, 20.0L, 10000, "Zeta, 80 bit long double");
+   plot_errors_1d(ei_func(), -20.0L, 20.0L, 10000, "Exponential Integral Ei, 80 bit long double");
+#endif
+#ifdef BOOST_HAS_FLOAT128
+   plot_errors_1d(digamma_func(), boost::multiprecision::float128(1e-200), boost::multiprecision::float128(10.0), 10000, "digamma, __float128");
+   plot_errors_1d(tgamma_func(), boost::multiprecision::float128(1e-200), boost::multiprecision::float128(150.0), 10000, "tgamma, __float128");
+   plot_errors_1d(lgamma_func(), boost::multiprecision::float128(1e-200), boost::multiprecision::float128(1000.0), 10000, "lgamma, __float128");
+   plot_errors_1d(trigamma_func(), boost::multiprecision::float128(1e-200), boost::multiprecision::float128(10.0), 10000, "trigamma, __float128");
+   plot_errors_1d(erf_func(), -boost::multiprecision::float128(5.0), boost::multiprecision::float128(5.0), 10000, "erf, __float128");
+   plot_errors_1d(erfc_func(), -boost::multiprecision::float128(5.0), boost::multiprecision::float128(120.0), 10000, "erfc, __float128");
+   plot_errors_1d(j0_func(), boost::multiprecision::float128(0.0), boost::multiprecision::float128(50.0), 10000, "j0, __float128", boost::multiprecision::float128(50.0));
+   plot_errors_1d(j1_func(), boost::multiprecision::float128(0.0), boost::multiprecision::float128(50.0), 10000, "j1, __float128", boost::multiprecision::float128(50.0));
+   plot_errors_1d(y0_func(), boost::multiprecision::float128(1e-100), boost::multiprecision::float128(50.0), 10000, "y0, __float128", boost::multiprecision::float128(50.0));
+   plot_errors_1d(y1_func(), boost::multiprecision::float128(1e-100), boost::multiprecision::float128(50.0), 10000, "y1, __float128", boost::multiprecision::float128(50.0));
+   plot_errors_1d(i0_func(), boost::multiprecision::float128(0.0), boost::multiprecision::float128(50.0), 10000, "i0, __float128");
+   plot_errors_1d(i1_func(), boost::multiprecision::float128(0.0), boost::multiprecision::float128(50.0), 10000, "i1, __float128");
+   plot_errors_1d(k0_func(), boost::multiprecision::float128(1e-100), boost::multiprecision::float128(50.0), 10000, "k0, __float128");
+   plot_errors_1d(k1_func(), boost::multiprecision::float128(1e-100), boost::multiprecision::float128(50.0), 10000, "k1, __float128");
+   plot_errors_1d(ai_func(), -boost::multiprecision::float128(20.0), boost::multiprecision::float128(20.0), 10000, "Ai, __float128", boost::multiprecision::float128(100.0));
+   plot_errors_1d(bi_func(), -boost::multiprecision::float128(20.0), boost::multiprecision::float128(20.0), 10000, "Bi, __float128", boost::multiprecision::float128(100.0));
+   plot_errors_1d(aip_func(), -boost::multiprecision::float128(20.0), boost::multiprecision::float128(20.0), 10000, "Ai Prime, __float128", boost::multiprecision::float128(100.0));
+   plot_errors_1d(bip_func(), -boost::multiprecision::float128(20.0), boost::multiprecision::float128(20.0), 10000, "Bi Prime, __float128", boost::multiprecision::float128(100.0));
+
+   plot_errors_1d(ellint_1_func(), -boost::multiprecision::float128(1.0), boost::multiprecision::float128(1.0), 10000, "Elliptic Integral K, __float128");
+   plot_errors_1d(ellint_2_func(), -boost::multiprecision::float128(1.0), boost::multiprecision::float128(1.0), 10000, "Elliptic Integral E, __float128");
+   plot_errors_1d(ellint_d_func(), -boost::multiprecision::float128(1.0), boost::multiprecision::float128(1.0), 10000, "Elliptic Integral D, __float128");
+
+   plot_errors_1d(zeta_func(), -boost::multiprecision::float128(20.0), boost::multiprecision::float128(20.0), 10000, "Zeta, __float128");
+   plot_errors_1d(ei_func(), -boost::multiprecision::float128(20.0), boost::multiprecision::float128(20.0), 10000, "Exponential Integral Ei, __float128");
+#endif
+   return 0;
+}