|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r83208 - in sandbox/e_float/doc: bessel_zeros derivative multiprecision_math sine_table
From: e_float_at_[hidden]
Date: 2013-02-28 15:08:37
Author: christopher_kormanyos
Date: 2013-02-28 15:08:35 EST (Thu, 28 Feb 2013)
New Revision: 83208
URL: http://svn.boost.org/trac/boost/changeset/83208
Log:
Added some documentation of Multiprecision examples.
Added:
sandbox/e_float/doc/bessel_zeros/ (props changed)
sandbox/e_float/doc/bessel_zeros/bessel_zeros.pdf (contents, props changed)
sandbox/e_float/doc/bessel_zeros/bessel_zeros.tcp (contents, props changed)
sandbox/e_float/doc/bessel_zeros/bessel_zeros.tex (contents, props changed)
sandbox/e_float/doc/derivative/ (props changed)
sandbox/e_float/doc/derivative/derivative.pdf (contents, props changed)
sandbox/e_float/doc/derivative/derivative.tcp (contents, props changed)
sandbox/e_float/doc/derivative/derivative.tex (contents, props changed)
sandbox/e_float/doc/multiprecision_math/ (props changed)
sandbox/e_float/doc/multiprecision_math/multiprecision_math.pdf (contents, props changed)
sandbox/e_float/doc/multiprecision_math/multiprecision_math.tcp (contents, props changed)
sandbox/e_float/doc/multiprecision_math/multiprecision_math.tex (contents, props changed)
sandbox/e_float/doc/sine_table/ (props changed)
sandbox/e_float/doc/sine_table/sine_table.tcp (contents, props changed)
sandbox/e_float/doc/sine_table/sine_table.tex (contents, props changed)
Added: sandbox/e_float/doc/bessel_zeros/bessel_zeros.pdf
==============================================================================
Binary file. No diff available.
Added: sandbox/e_float/doc/bessel_zeros/bessel_zeros.tcp
==============================================================================
--- (empty file)
+++ sandbox/e_float/doc/bessel_zeros/bessel_zeros.tcp 2013-02-28 15:08:35 EST (Thu, 28 Feb 2013)
@@ -0,0 +1,12 @@
+[FormatInfo]
+Type=TeXnicCenterProjectInformation
+Version=4
+
+[ProjectInfo]
+MainFile=bessel_zeros.tex
+UseBibTeX=1
+UseMakeIndex=0
+ActiveProfile=LaTeX ⇨ PDF
+ProjectLanguage=
+ProjectDialect=
+
Added: sandbox/e_float/doc/bessel_zeros/bessel_zeros.tex
==============================================================================
--- (empty file)
+++ sandbox/e_float/doc/bessel_zeros/bessel_zeros.tex 2013-02-28 15:08:35 EST (Thu, 28 Feb 2013)
@@ -0,0 +1,359 @@
+\documentclass{article}[10pt]
+
+\usepackage{amsmath}
+\usepackage{amsfonts}
+\usepackage{amssymb}
+\usepackage{graphicx}
+\usepackage{listings}
+
+\def\codedefault {\ttfamily}
+
+\lstset
+{
+ language=[ISO]C++,
+ morekeywords={interrupt,static_assert,constexpr,delete,auto,decltype,nullptr,noexcept},
+ framerule=0.40pt,
+ showstringspaces=false,
+ extendedchars=true,
+ basicstyle=\codedefault,
+ commentstyle=\codedefault\itshape,
+ keywordstyle=\codedefault\bfseries,
+ frame=tb,
+ aboveskip={1.1\baselineskip},
+ belowskip={1.1\baselineskip}
+}
+
+\def\cppnineeight {C++$98$}
+\def\cppothree {C++$03$}
+\def\cppox {C++$11$}
+\def\cninenine {C$99$}
+
+\def\trademarksymbolr {$^{\text{\rm{\scriptsize{\textregistered}}}}$}
+\def\trademarksymboltm {$^{\text{\rm{\scriptsize{TM}}}}$}
+
+\def\wolframalpha {{WolframAlpha\trademarksymbolr}}
+
+\begin{document}
+
+\title{Zeros of Cylindrical Bessel and Neumann Functions}
+\maketitle
+
+\noindent
+For every real order~$\nu$, cylindrical Bessel and Neumann
+functions have an infinite number of zeros on the positive
+real axis. The real zeros on the positive real axis
+can be found by solving for the roots of
+
+\begin{eqnarray}
+J_{\nu}
+\left(j_{\nu,\,m}\right)
+&=&
+0
+\nonumber \\
+Y_{\nu}
+\left(y_{\nu,\,m}\right)
+&=&
+0
+\,{\text{.}}
+\end{eqnarray}
+
+\noindent
+Here, $j_{\nu,\,m}$ represents the~$m^{\text{th}}$
+root of the cylindrical Bessel function of order~$\nu$
+and $y_{\nu,\,m}$ represents the~$m^{\text{th}}$
+root of the cylindrical Neumann function of order~$\nu$.
+
+Various methods are used to compute initial estimates
+for~$j_{\nu,\,m}$ and~$y_{\nu,\,m}$, and these will be described
+in detail below. After finding the initial estimate of a given root,
+its precision is subsequently refined to the desired level
+using Newton-Raphson iteration from \lstinline|Boost.Math|'s
+root-finding utilities combined with the functions
+\lstinline|cyl_bessel_j()| and \lstinline|cyl_neumann()|,
+also from \lstinline|Boost.Math|.
+
+Newton iteration requires both $J_{\nu}(x)$ (or~$Y_{\nu}(x)$)
+as well as its derivative. The derivatives of $J_{\nu}(x)$ and $Y_{\nu}(x)$
+with respect to~$x$ are given by
+Eq.~$9$.$1$.$3$ in~\cite{bibitem:abramowitz}.
+In particular,
+
+\begin{eqnarray}
+\dfrac{d}{dx}J_{\nu}(x)
+&=&
+J_{\nu - 1}(x)
+\,-\,
+\frac{\nu}{x}\,J_{\nu}(x)
+\nonumber \\
+\dfrac{d}{dx}Y_{\nu}(x)
+&=&
+Y_{\nu - 1}(x)
+\,-\,
+\frac{\nu}{x}\,Y_{\nu}(x)
+\,{\text{.}}
+\end{eqnarray}
+
+Enumeration of the rank of a root
+(in other words the index of a root)
+begins with one and counts up, in other words
+$m\,=\,1,\,2,\,3,\ldots$,~etc. The value of the
+first root is always greater than zero.
+
+For certain special parameters, cylindrical Bessel functions
+and cylindrical Neumann functions have
+a root at the origin. For example,
+$J_{\nu }(x)$ has a root at the origin for every positive order
+$\nu\,>\,0$ and for every negative integer order
+$\nu\,=\,-n$, with $n\,\in\,\mathbb{N}^{+}$
+and $n\,\ne\,0$.
+In addition, $Y_{\nu }(x)$ has a root at the origin
+for every negative half-integer order
+$\nu\,=\,-n/\,2$, with $n\,\in\,\mathbb{N}^{+}$
+and $n\,\ne\,0$.
+For these special parameter values, the origin with
+a value of~$x\,=\,0$ is provided as the $0^{\text{th}}$
+root generated by {\lstinline|cyl_bessel_j_zero()|}
+and {\lstinline|cyl_neumann_zero()|}.
+
+When calculating initial estimates for the roots
+of Bessel functions, a distinction is made between
+positive order and negative order, and different
+methods are used for these. In addition, different algorithms
+are used for the first root ($m\,=\,1$) and
+for subsequent roots with higher rank ($m\,\ge\,2$).
+Furthermore, estimates of the roots for Bessel functions
+with order above and below a cutoff at~$\nu\,=\,2.2$
+are calculated with different methods.
+
+Calculations of the estimates of $j_{\nu,\,1}$ and $y_{\nu,\,1}$
+with~$0\,\le\,\nu\,<\,2.2$ use empirically tabulated values.
+The coefficients for these have been generated by a
+computer algebra system.
+
+Calculations of the estimates of $j_{\nu,\,1}$ and $y_{\nu,\,1}$
+with~$\nu\,\ge\,2.2$ use
+Eqs.~$9$.$5$.$14$ and~$9$.$5$.$15$
+in~\cite{bibitem:abramowitz}.
+In particular,
+
+\begin{eqnarray}
+j_{\nu,\,1}
+& \sim &
+\nu
+\,+\,
+1.85575\,71\nu^{\frac{1}{3}}
+\,+\,
+1.03315\,0\nu^{-\frac{1}{3}}
+\,-\,
+0.00397\nu^{-1}
+\nonumber \\
+&&
+\,-\,
+0.0908\nu^{-\frac{5}{3}}
+\,+\,
+0.043\nu^{-\frac{7}{3}}
+\,+\,\cdots
+\,{\text{,}}
+\end{eqnarray}
+
+\noindent
+and
+
+\begin{eqnarray}
+y_{\nu,\,1}
+& \sim &
+\nu
+\,+\,
+0.93157\,68\nu^{\frac{1}{3}}
+\,+\,
+0.26035\,1\nu^{-\frac{1}{3}}
+\,+\,
+0.01198\nu^{-1}
+\nonumber \\
+&&
+\,-\,
+0.0060\nu^{-\frac{5}{3}}
+\,-\,
+0.001\nu^{-\frac{7}{3}}
+\,+\,\cdots
+\,{\text{,}}
+\end{eqnarray}
+
+Calculations of the estimates of $j_{\nu,\,m}$ and $y_{\nu,\,m}$
+with rank $m\,\ge\,2$ and~$0\,\le\,\nu\,<\,2.2$ use
+McMahon's approximation, as described in Sect.~$9$.$5$
+of~\cite{bibitem:abramowitz} and Eq.~$9$.$5$.$12$ therein.
+In particular,
+
+\begin{eqnarray}
+j_{\nu,\,m},\,y_{\nu,\,m}
+& \sim &
+\beta
+\,-\,
+\frac{\mu - 1}{8\beta}
+\,-\,
+\frac{4(\mu - 1)(7\mu -31)}{3(8\beta)^3}
+\nonumber \\
+&&
+\,-\,
+\frac{32(\mu - 1)(83\mu^2 - 982\mu +3779)}{15(8\beta)^5}
+\nonumber \\
+&&
+\!\!\!\! - \,
+\frac{64(\mu - 1)(6949\mu^3 - 1\,53855\mu^2 + 15\,85743\mu - 62\,77237)}{105(8\beta)^7}
+\nonumber \\
+&&
+\,-\,\cdots
+\,{\text{,}}
+\label{equation:mcmahon}
+\end{eqnarray}
+
+\noindent
+where $\mu\,=\,4\nu^2$
+and $\beta\,=\,\left(m\,+\,\frac{1}{2}\nu\,-\,\frac{1}{4}\right)\pi$
+for $j_{\nu,\,m}$
+and $\beta\,=\,\left(m\,+\,\frac{1}{2}\nu\,-\,\frac{3}{4}\right)\pi$
+for $y_{\nu,\,m}$.
+
+\vspace{2pt}
+
+Calculations of the estimates of $j_{\nu,\,m}$ and $y_{\nu,\,m}$
+with~$\nu\,\ge\,2.2$ use
+one term in the asymptotic expansion given in
+Eq.~$9$.$5$.$22$ and top line of Eq.~$9$.$5$.$26$
+combined with Eq.~$9$.$3$.$39$, all in~\cite{bibitem:abramowitz}.
+The latter two equations are expressed for
+argument~$x$ greater than one. In summary,
+
+\begin{equation}
+j_{\nu,\,m}
+\,\sim\,
+\nu x(-\zeta)
+\,+\,
+\dfrac{f_{1}(-\zeta)}{\nu}
+\,{\text{,}}
+\end{equation}
+
+\noindent
+where $-\zeta\,=\,\nu^{-2/\,3}a_{m}$ and $a_{m}$ is
+the absolute value of the $m^{\text{th}}$ root
+of $Ai(x)$ on the negative real axis.
+Here, $x\,=\,x(-\zeta)$ is the inverse of the function
+
+\begin{equation}
+\dfrac{2}{3}
+(-\zeta)^{3/\,2}
+\,=\,
+\sqrt{x^{2}\,-\,1}
+\,-\,
+\cos^{-1}\left(\dfrac{1}{x}\right)
+\,{\text{.}}
+\label{equation:uniform}
+\end{equation}
+
+\noindent
+Furthermore,
+
+\begin{equation}
+f_1(-\zeta)
+\,=\,
+\dfrac{1}{2}\,
+x(-\zeta)
+\bigl\{h(-\zeta)\bigr\}^{2}
+b_{0}(-\zeta)\,
+{\text{,}}
+\end{equation}
+
+\noindent
+where
+
+\begin{equation}
+h(-\zeta)
+\,=\,
+\left\{
+\dfrac{4(-\zeta)}{x^2\,-\,1}
+\right\}^{\frac{1}{4}}\,
+{\text{,}}
+\end{equation}
+
+\noindent
+and
+
+\begin{equation}
+b_{0}(-\zeta)
+\,=\,
+-\dfrac{5}{48\zeta^{2}}
+\,+\,
+\dfrac{1}{(-\zeta)^{\frac{1}{2}}}
+\Biggl\{
+\dfrac{5}{24(x^{2}\,-\,1)^{3/\,2}}
+\,+\,
+\dfrac{1}{8(x^{2}\,-\,1)^{1/\,2}}
+\Biggr\}
+{\text{.}}
+\end{equation}
+
+When solving for $x(-\zeta)$ in Eq.~\ref{equation:uniform} above,
+the right-hand-side is expanded to order--$(x^{2})$ in
+a Taylor series for large~$x$. This results in
+
+\begin{equation}
+\dfrac{2}{3}
+(-\zeta)^{3/\,2}
+\,\approx\,
+x
+\,+\,
+\dfrac{1}{2x}
+\,-\,
+\dfrac{\pi}{2}\,
+{\text{.}}
+\end{equation}
+
+\noindent
+The positive root of the resulting quadratic equation
+is used to find an initial estimate~$x(-\zeta)$.
+This initial estimate is subsequently refined with
+several steps of Newton-Raphson iteration
+in Eq.~\ref{equation:uniform}.
+
+Estimates of the roots of cylindrical Bessel functions
+of negative order on the positive real axis are found
+using interlacing relations. For example, the~$m^{\text{th}}$
+root of the cylindrical Bessel function $j_{-\nu,m}$
+is bracketed by the $m^{\text{th}}$ root and the
+$(m\,+\,1)^{\text{th}}$ root of the Bessel function of
+corresponding positive integer order. In other words,
+
+\begin{equation}
+j_{n_{\nu},m}
+\,<\,
+j_{-\nu,m}
+\,<\,
+j_{n_{\nu},m\,+\,1}
+\end{equation}
+
+\noindent
+where $m\,>\,1$ and $n_{\nu}$ represents the integral
+floor of the absolute value of~$\left|-\nu\right|$.
+
+Similar bracketing relations are used to find estimates
+of the roots of Neumann functions of negative order,
+whereby a discontinuity at every negative half-integer
+order needs to be handled.
+
+Bracketing relations do not hold for the first root
+of cylindrical Bessel functions and cylindrical Neumann
+functions with negative order. Therefore, iterative algorithms
+combined with root-finding via bisection are used
+to localize $j_{-\nu,1}$ and $y_{-\nu,1}$.
+
+\begin{thebibliography}{9}
+
+\bibitem{bibitem:abramowitz}
+M.~Abramowitz and I.~A.~Stegun:
+{\textit{Handbook of Mathematical Functions, $9^{th}$ Printing}},
+(Dover Publications, New York, 1972)
+
+\end{thebibliography}
+
+\end{document}
Added: sandbox/e_float/doc/derivative/derivative.pdf
==============================================================================
Binary file. No diff available.
Added: sandbox/e_float/doc/derivative/derivative.tcp
==============================================================================
--- (empty file)
+++ sandbox/e_float/doc/derivative/derivative.tcp 2013-02-28 15:08:35 EST (Thu, 28 Feb 2013)
@@ -0,0 +1,12 @@
+[FormatInfo]
+Type=TeXnicCenterProjectInformation
+Version=4
+
+[ProjectInfo]
+MainFile=derivative.tex
+UseBibTeX=1
+UseMakeIndex=0
+ActiveProfile=LaTeX ⇨ PDF
+ProjectLanguage=
+ProjectDialect=
+
Added: sandbox/e_float/doc/derivative/derivative.tex
==============================================================================
--- (empty file)
+++ sandbox/e_float/doc/derivative/derivative.tex 2013-02-28 15:08:35 EST (Thu, 28 Feb 2013)
@@ -0,0 +1,357 @@
+\documentclass{article}[10pt]
+
+\usepackage{amsmath}
+\usepackage{amsfonts}
+\usepackage{amssymb}
+\usepackage{graphicx}
+\usepackage{listings}
+
+\def\codedefault {\ttfamily}
+
+\lstset
+{
+ language=[ISO]C++,
+ morekeywords={interrupt,static_assert,constexpr,delete,auto,decltype,nullptr,noexcept},
+ framerule=0.40pt,
+ showstringspaces=false,
+ extendedchars=true,
+ basicstyle=\codedefault,
+ commentstyle=\codedefault\itshape,
+ keywordstyle=\codedefault\bfseries,
+ frame=tb,
+ aboveskip={1.1\baselineskip},
+ belowskip={1.1\baselineskip}
+}
+
+\def\cppnineeight {C++$98$}
+\def\cppothree {C++$03$}
+\def\cppox {C++$11$}
+\def\cninenine {C$99$}
+
+\def\trademarksymbolr {$^{\text{\rm{\scriptsize{\textregistered}}}}$}
+\def\trademarksymboltm {$^{\text{\rm{\scriptsize{TM}}}}$}
+
+\def\mathematica {{Mathematica\trademarksymbolr}}
+\def\wolframalpha {{WolframAlpha\trademarksymbolr}}
+
+\newcommand{\deriveval}[2][]{#1|_{#2}}
+
+\begin{document}
+
+\title{Countering Precision Loss with {\codedefault{Boost.Multiprecision}}}
+\maketitle
+
+\noindent
+The following example shows how multiprecision calculations
+can be used to obtain full precision in a numerical derivative
+calculation that suffers from precision loss.
+Consider some well-known central difference rules for
+numerically computing the first derivative of a function
+$f'(x)$ with $x\in\mathbb{R}$.
+
+\begin{eqnarray}
+f'(x)
+&\approx&
+m_{1}\,+\,O(dx^{2}) \nonumber \\
+f'(x)
+&\approx&
+\dfrac{4}{3}m_{1}\,-\,\dfrac{1}{3}m_{2}\,+\,O(dx^{4}) \nonumber \\
+f'(x)
+&\approx&
+\dfrac{3}{2}m_{1}\,-\,\dfrac{3}{5}m_{2}\,+\,\dfrac{1}{10}m_{3}\,+\,O(dx^{6})
+\,{\text{.}}
+\end{eqnarray}
+
+\noindent
+Here, the difference terms $m_{n}$ are given by
+
+\begin{eqnarray}
+m_{1}
+&=&
+\dfrac{f\left(x+dx\right)-f\left(x-dx\right)}{2dx} \nonumber \\
+m_{2}
+&=&
+\dfrac{f\left(x+2dx\right)-f\left(x-2dx\right)}{4dx} \nonumber \\
+m_{3}
+&=&
+\dfrac{f\left(x+3dx\right)-f\left(x-3dx\right)}{6dx}
+\,{\text{,}}
+\end{eqnarray}
+
+\noindent
+where $dx$ is the incremental step-size used when calculating
+the derivative. The third equation represents a three-point
+central difference rule with precision $O(dx^6)$.
+
+One possible implementation of this three-point
+central difference rule is shown below.
+
+\begin{lstlisting}
+template<typename value_type,
+ typename function_type>
+value_type derivative(const value_type x,
+ const value_type dx,
+ function_type function)
+{
+ // Compute the derivative of function using a
+ // three-point central difference rule of O(dx^6).
+
+ const value_type dx2(dx * 2U);
+ const value_type dx3(dx * 3U);
+
+ const value_type m1(( function(x + dx)
+ - function(x - dx)) / 2U);
+ const value_type m2(( function(x + dx2)
+ - function(x - dx2)) / 4U);
+ const value_type m3(( function(x + dx3)
+ - function(x - dx3)) / 6U);
+
+ const value_type fifteen_m1(m1 * 15U);
+ const value_type six_m2 (m2 * 6U);
+ const value_type ten_dx (dx * 10U);
+
+ return ((fifteen_m1 - six_m2) + m3) / ten_dx;
+}
+\end{lstlisting}
+
+Here, the implementation uses a C++ template that can be
+instantiated with various floating-point types
+such as {\lstinline|float|}, {\lstinline|double|},
+{\lstinline|long double|}, or even a user-defined
+floating-point type.
+
+We will now use the {\lstinline|derivative|} template
+with the built-in type {\lstinline|double|} in order to
+numerically compute the derivative of a function.
+Consider the function shown below.
+
+\begin{equation}
+f(x)
+\,=\,
+\sqrt{x^{2}\,-\,1}
+\,-\,
+\cos^{-1}\left(\dfrac{1}{x}\right)
+\end{equation}
+
+We will now take the derivative of this function with
+respect to~$x$ evaluated at $x\,=\,3/\,2$. In other words,
+
+\begin{equation}
+\dfrac{d}{dx}\,
+\Biggl[
+\sqrt{x^{2}\,-\,1}
+\,-\,
+\cos^{-1}\left(\dfrac{1}{x}\right)
+\Biggr]
+\,\deriveval[\Bigg]{x=\frac{3}{2}}
+\,{\text{.}}
+\end{equation}
+
+The expected result is
+
+\begin{equation}
+\dfrac{\sqrt{x^{2}\,-\,1}}{x}
+\,\deriveval[\Bigg]{x=\frac{3}{2}}
+\,=\,
+\dfrac{\sqrt{5}}{3}
+\,\approx\,
+0.74535\,59924\,99929\,89880
+\,{\text{.}}
+\end{equation}
+
+The program below uses the {\lstinline|derivative|} template
+in order to perform the numerical calculation of this derivative.
+The program also compares the numerically-obtained result
+with the expected result and reports the absolute relative error
+scaled to a deviation that can easily be related to the
+number of bits of lost precision.
+
+\begin{lstlisting}
+#include <iostream>
+#include <iomanip>
+#include <limits>
+#include <cmath>
+
+int main(void)
+{
+ const double d =
+ derivative(1.5,
+ std::ldexp(1.0, -9),
+ [](const double& x_) -> double
+ {
+ return std::sqrt((x_ * x_) - 1.0)
+ - std::acos(1.0 / x_);
+ });
+
+ const double rel_error =
+ (d - 0.74535599249992989880)
+ / 0.74535599249992989880;
+
+ const double bit_error =
+ std::abs(rel_error)
+ / std::numeric_limits<double>::epsilon();
+
+ std::cout.precision(
+ std::numeric_limits<double>::digits10);
+
+ std::cout << "derivative: "
+ << d
+ << std::endl;
+
+ std::cout << "expected : "
+ << 0.74535599249992989880
+ << std::endl;
+
+ std::cout << "bit_error : "
+ << unsigned long(bit_error)
+ << std::endl;
+}
+\end{lstlisting}
+
+The result of this program on a system with an eight-byte,
+IEEE-754 conforming floating-point representation for
+{\lstinline|double|} is:
+
+\begin{lstlisting}
+derivative: 0.745355992499951
+expected : 0.74535599249993
+bit_error : 130
+\end{lstlisting}
+
+Here, the {\lstinline|bit_error|} variable
+is calculated from the absolute
+relative error scaled with $\epsilon$, where $\epsilon$
+is {\lstinline|std::numeric_limits<double>::epsilon()|}
+and is equal to the smallest number that differs from
+one that is representable in the floating-point system.
+The derivative has a {\lstinline|bit_error|} of~130,
+corresponding to approximately~7 bits of precision loss.
+
+In this example, we have carefully selected a step-size
+$dx$ of $2^{-9}$ using the library function {\lstinline|ldexp()|}.
+This is done because numerical differentiation
+generally produces the best results when a step-size
+that is exactly representable in the floating-point
+system is used, in other words a pure power of~2
+for an IEEE-754 representation.
+
+If the calculation is repeated using step sizes of
+$2^{-n}$ with $n\,=\,7,\,8,\,9,\,10$, and~$11$,
+the best result is in fact found using
+a step-size of~$2^{-9}$.
+Empirical experiments with other step-sizes in various
+bases on this system have all resulted in worse values.
+So using the step-size of~$2^{-9}$ may, in fact,
+produce the best numerical result that can be obtained
+for the built-in type {\lstinline|double|} when used for
+this numerical derivative on this system.
+
+With $dx\,=\,2^{-9}$, the derivative central difference rule
+produces a very good result. It still, however, suffers
+from approximately~7 bits of precision. A better result
+can be obtained by using a multiprecision type from
+{\lstinline|Boost.Multiprecision|} such as
+{\lstinline|boost::multiprecision::cpp_dec_float|}.
+
+We will now repeat the numerical derivative calculation
+shown above using a multiprecision type. First, we need to
+define a multiprecision type using the facilities in
+{\lstinline|Boost.Multiprecision|}.
+
+\begin{lstlisting}
+#include <boost/multiprecision/cpp_dec_float.hpp>
+
+using boost::multiprecision::number;
+using boost::multiprecision::cpp_dec_float;
+
+#define MP_DIGITS10 \
+ unsigned( \
+ std::numeric_limits<double>::max_digits10 + 5)
+
+typedef cpp_dec_float<MP_DIGITS10> mp_backend;
+
+typedef number<mp_backend> mp_type;
+\end{lstlisting}
+
+Here, the user-defined data type {\lstinline|mp_type|}
+has been defined with~5 decimal digits of precision
+more than the built-in type {\lstinline|double|}.
+
+The following program repeats the derivative calculation.
+This time, however, the multiprecision type {\lstinline|mp_type|}
+is used. In particular,
+
+\begin{lstlisting}
+int main(void)
+{
+ const mp_type mp =
+ derivative(mp_type(mp_type(3) / 2U),
+ mp_type(mp_type(1) / 10000000U),
+ [](const mp_type& x_) -> mp_type
+ {
+ return sqrt((x_ * x_) - 1.0)
+ - acos(1.0 / x_);
+ });
+
+ const double d = mp.convert_to<double>();
+
+ const double rel_error =
+ (d - 0.74535599249992989880)
+ / 0.74535599249992989880;
+
+ const double bit_error =
+ std::abs(rel_error)
+ / std::numeric_limits<double>::epsilon();
+
+ std::cout.precision(
+ std::numeric_limits<double>::digits10);
+
+ std::cout << "derivative: "
+ << d
+ << std::endl;
+
+ std::cout << "expected : "
+ << 0.74535599249992989880
+ << std::endl;
+
+ std::cout << "bit_error : "
+ << unsigned long(bit_error)
+ << std::endl;
+\end{lstlisting}
+
+Note that the derivative calculation is carried out
+with the multiprecision type. The result of the
+multiprecision calculation is subsequently converted to
+the built-in type {\lstinline|double|} using the
+template {\lstinline|convert_to()|} member function.
+
+The result of this program is:
+
+\begin{lstlisting}
+derivative: 0.74535599249993
+expected : 0.74535599249993
+bit_error : 0
+\end{lstlisting}
+
+The resulting bit error is~0. This means that
+the result of the derivative calculation is bit-identical
+with the {\lstinline|double|} representation of the
+expected result, and this is the best result possible
+for the built-in type.
+
+The derivative in this example has a known closed form.
+There are, however, countless situations in numerical
+analysis (and not only for numerical derivatives)
+for which the calculation at hand does not have a known
+closed-form solution or for which the closed-form solution
+is highly inconvenient to use. In such cases, this technique
+may be useful.
+
+This example has shown how multiprecision can be used
+to add extra digits to an ill-conditioned calculation
+that suffers from precision loss. When the result of the
+multiprecision calculation is converted to a built-in type
+such as {\lstinline|double|}, the entire precision of the
+result in {\lstinline|double|} is preserved.
+
+\end{document}
\ No newline at end of file
Added: sandbox/e_float/doc/multiprecision_math/multiprecision_math.pdf
==============================================================================
Binary file. No diff available.
Added: sandbox/e_float/doc/multiprecision_math/multiprecision_math.tcp
==============================================================================
--- (empty file)
+++ sandbox/e_float/doc/multiprecision_math/multiprecision_math.tcp 2013-02-28 15:08:35 EST (Thu, 28 Feb 2013)
@@ -0,0 +1,12 @@
+[FormatInfo]
+Type=TeXnicCenterProjectInformation
+Version=4
+
+[ProjectInfo]
+MainFile=multiprecision_math.tex
+UseBibTeX=1
+UseMakeIndex=0
+ActiveProfile=LaTeX ⇨ PDF
+ProjectLanguage=
+ProjectDialect=
+
Added: sandbox/e_float/doc/multiprecision_math/multiprecision_math.tex
==============================================================================
--- (empty file)
+++ sandbox/e_float/doc/multiprecision_math/multiprecision_math.tex 2013-02-28 15:08:35 EST (Thu, 28 Feb 2013)
@@ -0,0 +1,407 @@
+\documentclass{article}[10pt]
+
+\usepackage{amsmath}
+\usepackage{amsfonts}
+\usepackage{amssymb}
+\usepackage{graphicx}
+\usepackage{listings}
+
+\def\codedefault {\ttfamily}
+
+\lstset
+{
+ language=[ISO]C++,
+ morekeywords={interrupt,static_assert,constexpr,delete,auto,decltype,nullptr,noexcept},
+ framerule=0.40pt,
+ showstringspaces=false,
+ extendedchars=true,
+ basicstyle=\codedefault,
+ commentstyle=\codedefault\itshape,
+ keywordstyle=\codedefault\bfseries,
+ frame=tb,
+ aboveskip={1.1\baselineskip},
+ belowskip={1.1\baselineskip}
+}
+
+\def\cppnineeight {C++$98$}
+\def\cppothree {C++$03$}
+\def\cppox {C++$11$}
+\def\cninenine {C$99$}
+
+\def\trademarksymbolr {$^{\text{\rm{\scriptsize{\textregistered}}}}$}
+\def\trademarksymboltm {$^{\text{\rm{\scriptsize{TM}}}}$}
+
+\def\mathematica {{Mathematica\trademarksymbolr}}
+\def\wolframalpha {{WolframAlpha\trademarksymbolr}}
+
+\begin{document}
+
+\title{Using {\codedefault{Boost.Math}} with {\codedefault{Boost.Multiprecision}}}
+\maketitle
+
+\noindent
+We will now describe a variety of ways to use
+\lstinline|Boost.Math| in combination with
+\lstinline|Boost.Multiprecision|
+in order to perform floating-point numerical calculations
+with high precision.
+
+The \lstinline|Boost.Multiprecision| library can be used
+for computations requiring precision exceeding that
+of standard built-in types such as
+\lstinline|float|, \lstinline|double| and \lstinline|long double|.
+For extended-precision calculations,
+\lstinline|Boost.Multiprecision| supplies a template data
+type called \lstinline|cpp_dec_float|.
+The number of decimal digits of precision is fixed at compile-time
+via template parameter. For example,
+
+\begin{lstlisting}
+#include <boost/multiprecision/cpp_dec_float.hpp>
+
+using boost::multiprecision::cpp_dec_float;
+
+// A floating-point type with 50 digits of precision.
+typedef cpp_dec_float<50> float50;
+\end{lstlisting}
+
+Here, we have defined the local data type
+\lstinline|float50| with with~$50$ decimal digits of
+precision. We can use this data type to print, for example,
+the value of $1/\,7$ to~$50$ digits.
+
+\begin{lstlisting}
+#include <iostream>
+#include <limits>
+
+// 1/7 with 50 digits of precision.
+float50 seventh = float50(1) / 7;
+
+const int d
+ = std::numeric_limits<float50>::digits10;
+
+// Print 1/7 with 50 digits of precision.
+std::cout << std::setprecision(d)
+ << seventh
+ << std::endl;
+\end{lstlisting}
+
+For the sake of convenience, \lstinline|Boost.Multiprecision|
+defines a variety of data types with fixed precision.
+These include, among others, data types with~$50$ and~$100$
+decimal digits of precision, respectively. In particular,
+
+\begin{lstlisting}
+using boost::multiprecision;
+
+cpp_dec_float_50 f50; // Has 50 digits.
+cpp_dec_float_100 f100; // Has 100 digits.
+\end{lstlisting}
+
+For more information on the \lstinline|cpp_dec_float| data type,
+see the tutorial in \lstinline|Boost.Multiprecision|.
+
+Mathematical software usually requires one or more
+known constant values such as~$\pi$ or $\log{2}$,~etc.
+It often makes sense
+to initialize a constant value stored in a built-in type
+from a multiprecision value of higher precision.
+This guarantees that the built-in type will be initialized
+to the the very last bit of its own precision.
+
+We will now obtain the value of~$\pi$ with~$50$ decimal
+digits of precision from \lstinline|Boost.Math|'s
+template \lstinline|pi()| function.
+
+\begin{lstlisting}
+using boost::multiprecision;
+
+cpp_dec_float_50 pi_50;
+ = boost::math::constants::pi<cpp_dec_float_50>();
+\end{lstlisting}
+
+A multiprecision value can be converted to a built-in type
+using the member function \lstinline|convert_to()|.
+In particular, the code below initializes a
+value of type \lstinline|long double| from
+a multiprecision value having~$50$ digits of precision.
+
+\begin{lstlisting}
+const long double pi_ld
+ = pi_50.convert_to<long_double>();
+\end{lstlisting}
+
+One usually needs to compute tables of numbers in
+mathematical software. A fast Fourier transform (FFT),
+for example, may use a table of the values of
+$\sin\left({\pi /\,2^n}\right)$
+in its implementation details. In order to maximize the
+precision in the FFT implementation, the precision of the
+tabulated trigonometric values should exceed that of the
+built-in floating-point type used in the FFT.
+
+The sample below computes a table of the values of
+$\sin\left({\pi /\,2^n}\right)$ in the range
+$1\,\le\,n\le\,31$.
+A precision of~$50$ decimal digits is used.
+
+\begin{lstlisting}
+#include <vector>
+#include <algorithm>
+#include <iostream>
+#include <iomanip>
+#include <iterator>
+#include <limits>
+#include <boost/multiprecision/cpp_dec_float.hpp>
+
+using boost::multiprecision::cpp_dec_float_50;
+using boost::math::constants::pi;
+
+typedef cpp_dec_float_50 float50;
+
+int main(int, char**)
+{
+ std::vector<float50> sin_values(31U);
+
+ unsigned n = 1U;
+
+ // Generate the sine values.
+ std::for_each(
+ sin_values.begin(),
+ sin_values.end(),
+ [&n](float50& y)
+ {
+ y = sin(pi<float50>() / pow(float50(2), n));
+ ++n;
+ });
+
+ // Set the output precision.
+ const int d
+ = std::numeric_limits<float50>::digits10;
+
+ std::cout.precision(d);
+
+ // Print the sine table.
+ std::copy(sin_values.begin(),
+ sin_values.end(),
+ std::ostream_iterator<float50>(std::cout, "\n"));
+}
+\end{lstlisting}
+
+This program makes use of, among other program elements,
+the data type \lstinline|cpp_dec_float_50| from
+\lstinline|Boost.Multiprecision|, the value of $\pi$
+retrieved from \lstinline|Boost.Math|
+and a \cppox\ lambda function combined with
+\lstinline|std::for_each()|.
+
+\pagebreak
+
+The output of this program is shown below.
+
+\begin{lstlisting}
+1
+0.70710678118654752440084436210484903928483593768847
+0.38268343236508977172845998403039886676134456248563
+0.19509032201612826784828486847702224092769161775195
+0.098017140329560601994195563888641845861136673167501
+0.049067674327418014254954976942682658314745363025753
+0.024541228522912288031734529459282925065466119239451
+0.012271538285719926079408261951003212140372319591769
+0.0061358846491544753596402345903725809170578863173913
+0.003067956762965976270145365490919842518944610213452
+0.0015339801862847656123036971502640790799548645752374
+0.00076699031874270452693856835794857664314091945206328
+0.00038349518757139558907246168118138126339502603496474
+0.00019174759731070330743990956198900093346887403385916
+9.5873799095977345870517210976476351187065612851145e-05
+4.7936899603066884549003990494658872746866687685767e-05
+2.3968449808418218729186577165021820094761474895673e-05
+1.1984224905069706421521561596988984804731977538387e-05
+5.9921124526424278428797118088908617299871778780951e-06
+2.9960562263346607504548128083570598118251878683408e-06
+1.4980281131690112288542788461553611206917585861527e-06
+7.4901405658471572113049856673065563715595930217207e-07
+3.7450702829238412390316917908463317739740476297248e-07
+1.8725351414619534486882457659356361712045272098287e-07
+9.3626757073098082799067286680885620193236507169473e-08
+4.681337853654909269511551813854009695950362701667e-08
+2.3406689268274552759505493419034844037886207223779e-08
+1.1703344634137277181246213503238103798093456639976e-08
+5.8516723170686386908097901008341396943900085051757e-09
+2.9258361585343193579282304690689559020175857150074e-09
+1.4629180792671596805295321618659637103742615227834e-09
+\end{lstlisting}
+
+This output can be copied
+as text and readily integrated into a given source code file.
+Alternatively, the output can be written to a text file or
+even be used as part of a self-written automatic code generator.
+
+A computer algebra system can be used to verify
+the results obtained from
+\lstinline|Boost.Math| and \lstinline|Boost.Multiprecision|.
+For example, the \mathematica\
+computer algebra system~\cite{bibitem:mathematica}
+can obtain a similar table with the command:
+
+\begin{verbatim}
+Table[N[Sin[Pi / (2^n)], 50], {n, 1, 31, 1}]
+\end{verbatim}
+
+The \wolframalpha\
+Computational Knowledge Engine~\cite{bibitem:wolframalpha}
+can also be used to generate this table.
+The same command can be used.
+
+In general, calculations of special functions inherently
+lose a few bits of precision (or more)
+due to cancellations in series expansions,
+recurrence relations and the like. It often arises, however,
+that data tables originating from
+special function calculations are used as the basis
+of an even higher-level calculation. Computing data tables
+with extended precision can improve the overall accuracy of
+these kinds of calculations because the data table
+has more precision than necessary.
+
+It may be wise to warm-cache pre-computed high-precision
+table values in a static container such as \lstinline|std::vector|.
+It can also be convenient to store pre-computed high-precision
+data tables in text-based files that can be parsed whenever
+the values need to be initialized and warm-cached.
+
+For our final example, we will use \lstinline|Boost.Math| and
+\lstinline|Boost.Multiprecision| to generate a table of
+high-precision real-valued roots of
+cylindrical Bessel functions $J_{\nu}(x)$ on the
+positive real axis.
+
+The roots of cylindrical Bessel functions can be found
+by solving
+
+\begin{equation}
+J_{\nu}\left(j_{\nu,\,m}\right)
+\,\,=\,\,
+0
+\,{\text{,}}
+\end{equation}
+
+\noindent
+where $j_{\nu,\,m}$ represents the~$m^{\text{th}}$
+root of the cylindrical Bessel function of order~$\nu$.
+
+The calculation uses McMahon's approximation to find
+initial estimates for~$j_{\nu,\,m}$,
+as described in Sect.~$10$.$21$({\ttfamily{vi}})
+of~\cite{bibitem:nisthandbook} and the references therein.
+The precision of a given root is subsequently
+refined to the desired level using Newton iteration.
+
+McMahon's asymptotic approximation for $j_{\nu,\,m}$
+with $\nu\,\ge\,0$ and $m\,\rightarrow\,\infty$
+is given by Eq.~$10$.$21$.$19$ in~\cite{bibitem:nisthandbook}
+
+\begin{eqnarray}
+j_{\nu,\,m}
+& = &
+a
+\,-\,
+\frac{\mu - 1}{8a}
+\,-\,
+\frac{4(\mu - 1)(7\mu -31)}{3(8a)^3}
+\,-\,
+\frac{32(\mu - 1)(83\mu^2 - 982\mu +3779)}{15(8a)^5}
+\nonumber \\
+&&
+\!\!\!\!\!\!\!\! - \,
+\frac{64(\mu - 1)(6949\mu^3 - 1\,53855\mu^2 + 15\,85743\mu - 62\,77237)}{105(8a)^7}
+\,-\,\cdots
+\,{\text{,}}
+\label{equation:mcmahon}
+\end{eqnarray}
+
+\noindent
+where $\mu\,=\,4\nu^2$
+and
+$a\,=\,\left(m\,+\,\frac{1}{2}\nu\,-\,\frac{1}{4}\right)\pi$.
+
+\vspace{2pt}
+
+The specialized expression for $j_{\nu,\,1}$
+is explicitly given by
+Eq.~$10$.$21$.$19$ in~\cite{bibitem:nisthandbook}
+
+\begin{eqnarray}
+j_{\nu,\,1}
+& = &
+\nu
+\,+\,
+1.85575\,71\nu^{\frac{1}{3}}
+\,+\,
+1.03315\,0\nu^{-\frac{1}{3}}
+\,-\,
+0.00397\nu^{-1}
+\nonumber \\
+&&
+\,-\,
+0.0908\nu^{-\frac{5}{3}}
+\,+\,
+0.043\nu^{-\frac{7}{3}}
+\,+\,\cdots
+\,{\text{,}}
+\end{eqnarray}
+
+Newton iteration requires both $J_{\nu}(x)$ as well as
+its derivative. The derivative of $J_{\nu}(x)$
+with respect to~$x$ is given by
+Eq.~$10$.$6$.$2$ in~\cite{bibitem:nisthandbook}
+
+\begin{equation}
+\dfrac{d}{dx}J_{\nu}(x)
+\,=\,
+J_{\nu - 1}(x)
+\,-\,
+\frac{\nu}{x}\,J_{\nu}(x)
+\,{\text{.}}
+\end{equation}
+
+The sample code shown below can be used to find
+the roots of $J_{\nu}(x)$. This program
+can also be found in the link here.
+In this example, a subroutine called
+\lstinline|cyl_bessel_j_zero()| is created
+and used to compute the first~$20$ roots of~$J_\frac{71}{19}(x)$
+for~$x$ on the positive real axis with~$50$ decimal digits
+of precision.
+
+\vspace{2pt}
+
+TBD: Make and show sample code.
+
+The output of the program is shown below.
+
+TBD: Show the program output.
+
+The results of this program can be verified with
+\mathematica\ using the following command:
+
+\begin{verbatim}
+Table[N[BesselJZero[71/19, n], 50], {n, 1, 20, 1}]
+\end{verbatim}
+
+\begin{thebibliography}{9}
+\bibitem{bibitem:mathematica}
+Wolfram: {\textit{Wolfram \mathematica}},\\
+{\ttfamily{http://www.wolfram.com/mathematica}} (2013)
+\bibitem{bibitem:wolframalpha}
+Wolfram: {\textit{\wolframalpha\ Comutational Knowledge Engine}},\\
+{\ttfamily{http://www.wolframalpha.com}} (2013)
+\bibitem{bibitem:nisthandbook}
+F.~W.~J.~Olver, D.~W. Lozier, R.~F. Boisvert and C.~W. Clark:
+{\textit{NIST Handbook of Mathematical Functions}},
+(NIST National Instutute of Standards and Technology
+US Department of Commerce and Cambridge University Press, New York, 2010)
+\end{thebibliography}
+
+\end{document}
Added: sandbox/e_float/doc/sine_table/sine_table.tcp
==============================================================================
--- (empty file)
+++ sandbox/e_float/doc/sine_table/sine_table.tcp 2013-02-28 15:08:35 EST (Thu, 28 Feb 2013)
@@ -0,0 +1,12 @@
+[FormatInfo]
+Type=TeXnicCenterProjectInformation
+Version=4
+
+[ProjectInfo]
+MainFile=sine_table.tex
+UseBibTeX=1
+UseMakeIndex=0
+ActiveProfile=LaTeX ⇨ PDF
+ProjectLanguage=
+ProjectDialect=
+
Added: sandbox/e_float/doc/sine_table/sine_table.tex
==============================================================================
--- (empty file)
+++ sandbox/e_float/doc/sine_table/sine_table.tex 2013-02-28 15:08:35 EST (Thu, 28 Feb 2013)
@@ -0,0 +1,285 @@
+\documentclass{article}[10pt]
+
+\usepackage{amsmath}
+\usepackage{amsfonts}
+\usepackage{amssymb}
+\usepackage{graphicx}
+\usepackage{listings}
+
+\def\codedefault {\ttfamily}
+
+\lstset
+{
+ language=[ISO]C++,
+ morekeywords={interrupt,static_assert,constexpr,delete,auto,decltype,nullptr,noexcept},
+ framerule=0.40pt,
+ showstringspaces=false,
+ extendedchars=true,
+ basicstyle=\codedefault,
+ commentstyle=\codedefault\itshape,
+ keywordstyle=\codedefault\bfseries,
+ frame=tb,
+ aboveskip={1.1\baselineskip},
+ belowskip={1.1\baselineskip}
+}
+
+\def\cppnineeight {C++$98$}
+\def\cppothree {C++$03$}
+\def\cppox {C++$11$}
+\def\cninenine {C$99$}
+
+\def\trademarksymbolr {$^{\text{\rm{\scriptsize{\textregistered}}}}$}
+\def\trademarksymboltm {$^{\text{\rm{\scriptsize{TM}}}}$}
+
+\def\mathematica {{Mathematica\trademarksymbolr}}
+\def\wolframalpha {{WolframAlpha\trademarksymbolr}}
+
+\newcommand{\deriveval}[2][]{#1|_{#2}}
+
+\begin{document}
+
+\title{Making a Sine Table with {\codedefault{Boost.Multiprecision}}}
+\maketitle
+
+\noindent
+We will now describe a variety of ways to use
+\lstinline|Boost.Math| in combination with
+\lstinline|Boost.Multiprecision|
+in order to perform floating-point numerical calculations
+with high precision.
+
+The \lstinline|Boost.Multiprecision| library can be used
+for computations requiring precision exceeding that
+of standard built-in types such as
+\lstinline|float|, \lstinline|double| and \lstinline|long double|.
+For extended-precision calculations,
+\lstinline|Boost.Multiprecision| supplies a template data
+type called \lstinline|cpp_dec_float|.
+The number of decimal digits of precision is fixed at compile-time
+via template parameter. For example,
+
+\begin{lstlisting}
+#include <boost/multiprecision/cpp_dec_float.hpp>
+
+using boost::multiprecision::cpp_dec_float;
+
+// A floating-point type with 50 digits of precision.
+typedef cpp_dec_float<50> float50;
+\end{lstlisting}
+
+Here, we have defined the local data type
+\lstinline|float50| with with~$50$ decimal digits of
+precision. We can use this data type to print, for example,
+the value of $1/\,7$ to~$50$ digits.
+
+\begin{lstlisting}
+#include <iostream>
+#include <limits>
+
+// 1/7 with 50 digits of precision.
+float50 seventh = float50(1) / 7;
+
+const int d
+ = std::numeric_limits<float50>::digits10;
+
+// Print 1/7 with 50 digits of precision.
+std::cout << std::setprecision(d)
+ << seventh
+ << std::endl;
+\end{lstlisting}
+
+For the sake of convenience, \lstinline|Boost.Multiprecision|
+defines a variety of data types with fixed precision.
+These include, among others, data types with~$50$ and~$100$
+decimal digits of precision, respectively. In particular,
+
+\begin{lstlisting}
+using boost::multiprecision;
+
+cpp_dec_float_50 f50; // Has 50 digits.
+cpp_dec_float_100 f100; // Has 100 digits.
+\end{lstlisting}
+
+For more information on the \lstinline|cpp_dec_float| data type,
+see the tutorial in \lstinline|Boost.Multiprecision|.
+
+Mathematical software usually requires one or more
+known constant values such as~$\pi$ or $\log{2}$,~etc.
+It often makes sense
+to initialize a constant value stored in a built-in type
+from a multiprecision value of higher precision.
+This guarantees that the built-in type will be initialized
+to the the very last bit of its own precision.
+
+We will now obtain the value of~$\pi$ with~$50$ decimal
+digits of precision from \lstinline|Boost.Math|'s
+template \lstinline|pi()| function.
+
+\begin{lstlisting}
+using boost::multiprecision;
+
+cpp_dec_float_50 pi_50;
+ = boost::math::constants::pi<cpp_dec_float_50>();
+\end{lstlisting}
+
+A multiprecision value can be converted to a built-in type
+using the member function \lstinline|convert_to()|.
+In particular, the code below initializes a
+value of type \lstinline|long double| from
+a multiprecision value having~$50$ digits of precision.
+
+\begin{lstlisting}
+const long double pi_ld
+ = pi_50.convert_to<long_double>();
+\end{lstlisting}
+
+One usually needs to compute tables of numbers in
+mathematical software. A fast Fourier transform (FFT),
+for example, may use a table of the values of
+$\sin\left({\pi /\,2^n}\right)$
+in its implementation details. In order to maximize the
+precision in the FFT implementation, the precision of the
+tabulated trigonometric values should exceed that of the
+built-in floating-point type used in the FFT.
+
+The sample below computes a table of the values of
+$\sin\left({\pi /\,2^n}\right)$ in the range
+$1\,\le\,n\le\,31$.
+A precision of~$50$ decimal digits is used.
+
+\begin{lstlisting}
+#include <vector>
+#include <algorithm>
+#include <iostream>
+#include <iomanip>
+#include <iterator>
+#include <limits>
+#include <boost/multiprecision/cpp_dec_float.hpp>
+
+using boost::multiprecision::cpp_dec_float_50;
+using boost::math::constants::pi;
+
+typedef cpp_dec_float_50 float50;
+
+int main(int, char**)
+{
+ std::vector<float50> sin_values(31U);
+
+ unsigned n = 1U;
+
+ // Generate the sine values.
+ std::for_each(
+ sin_values.begin(),
+ sin_values.end(),
+ [&n](float50& y)
+ {
+ y = sin(pi<float50>() / pow(float50(2), n));
+ ++n;
+ });
+
+ // Set the output precision.
+ const int d
+ = std::numeric_limits<float50>::digits10;
+
+ std::cout.precision(d);
+
+ // Print the sine table.
+ std::copy(sin_values.begin(),
+ sin_values.end(),
+ std::ostream_iterator<float50>(std::cout, "\n"));
+}
+\end{lstlisting}
+
+This program makes use of, among other program elements,
+the data type \lstinline|cpp_dec_float_50| from
+\lstinline|Boost.Multiprecision|, the value of $\pi$
+retrieved from \lstinline|Boost.Math|
+and a \cppox\ lambda function combined with
+\lstinline|std::for_each()|.
+
+\pagebreak
+
+The output of this program is shown below.
+
+\begin{lstlisting}
+1
+0.70710678118654752440084436210484903928483593768847
+0.38268343236508977172845998403039886676134456248563
+0.19509032201612826784828486847702224092769161775195
+0.098017140329560601994195563888641845861136673167501
+0.049067674327418014254954976942682658314745363025753
+0.024541228522912288031734529459282925065466119239451
+0.012271538285719926079408261951003212140372319591769
+0.0061358846491544753596402345903725809170578863173913
+0.003067956762965976270145365490919842518944610213452
+0.0015339801862847656123036971502640790799548645752374
+0.00076699031874270452693856835794857664314091945206328
+0.00038349518757139558907246168118138126339502603496474
+0.00019174759731070330743990956198900093346887403385916
+9.5873799095977345870517210976476351187065612851145e-05
+4.7936899603066884549003990494658872746866687685767e-05
+2.3968449808418218729186577165021820094761474895673e-05
+1.1984224905069706421521561596988984804731977538387e-05
+5.9921124526424278428797118088908617299871778780951e-06
+2.9960562263346607504548128083570598118251878683408e-06
+1.4980281131690112288542788461553611206917585861527e-06
+7.4901405658471572113049856673065563715595930217207e-07
+3.7450702829238412390316917908463317739740476297248e-07
+1.8725351414619534486882457659356361712045272098287e-07
+9.3626757073098082799067286680885620193236507169473e-08
+4.681337853654909269511551813854009695950362701667e-08
+2.3406689268274552759505493419034844037886207223779e-08
+1.1703344634137277181246213503238103798093456639976e-08
+5.8516723170686386908097901008341396943900085051757e-09
+2.9258361585343193579282304690689559020175857150074e-09
+1.4629180792671596805295321618659637103742615227834e-09
+\end{lstlisting}
+
+This output can be copied
+as text and readily integrated into a given source code file.
+Alternatively, the output can be written to a text file or
+even be used as part of a self-written automatic code generator.
+
+A computer algebra system can be used to verify
+the results obtained from
+\lstinline|Boost.Math| and \lstinline|Boost.Multiprecision|.
+For example, the \mathematica\
+computer algebra system~\cite{bibitem:mathematica}
+can obtain a similar table with the command:
+
+\begin{verbatim}
+Table[N[Sin[Pi / (2^n)], 50], {n, 1, 31, 1}]
+\end{verbatim}
+
+The \wolframalpha\
+Computational Knowledge Engine~\cite{bibitem:wolframalpha}
+can also be used to generate this table.
+The same command can be used.
+
+In general, calculations of special functions inherently
+lose a few bits of precision (or more)
+due to cancellations in series expansions,
+recurrence relations and the like. It often arises, however,
+that data tables originating from
+special function calculations are used as the basis
+of an even higher-level calculation. Computing data tables
+with extended precision can improve the overall accuracy of
+these kinds of calculations because the data table
+has more precision than necessary.
+
+It may be wise to warm-cache pre-computed high-precision
+table values in a static container such as \lstinline|std::vector|.
+It can also be convenient to store pre-computed high-precision
+data tables in text-based files that can be parsed whenever
+the values need to be initialized and warm-cached.
+
+\begin{thebibliography}{9}
+\bibitem{bibitem:mathematica}
+Wolfram: {\textit{Wolfram \mathematica}},\\
+{\ttfamily{http://www.wolfram.com/mathematica}} (2013)
+\bibitem{bibitem:wolframalpha}
+Wolfram: {\textit{\wolframalpha\ Comutational Knowledge Engine}},\\
+{\ttfamily{http://www.wolframalpha.com}} (2013)
+\end{thebibliography}
+
+\end{document}
Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk