boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Precision updated, quadrature notes added, need prime notes

Browse Source
This commit is contained in:
pabristow
2018-07-19 16:48:04 +01:00
parent aaa38c3fee
commit 8d77dc0516
10 changed files with 205 additions and 55 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/html/index.html
View File

@@ -120,7 +120,7 @@ This manual is also available in <a href="http://sourceforge.net/projects/boost/
</div>
</div>
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
<td align="left"><p><small>Last revised: July 17, 2018 at 16:28:57 GMT</small></p></td>
<td align="left"><p><small>Last revised: July 19, 2018 at 15:38:32 GMT</small></p></td>
<td align="right"><div class="copyright-footer"></div></td>
</tr></table>
<hr>

2
doc/html/indexes/s01.html
View File

@@ -24,7 +24,7 @@
</div>
<div class="section">
<div class="titlepage"><div><div><h2 class="title" style="clear: both">
<a name="id1865620"></a>Function Index</h2></div></div></div>
<a name="id1854868"></a>Function Index</h2></div></div></div>
<p><a class="link" href="s01.html#idx_id_0">2</a> <a class="link" href="s01.html#idx_id_1">4</a> <a class="link" href="s01.html#idx_id_2">A</a> <a class="link" href="s01.html#idx_id_3">B</a> <a class="link" href="s01.html#idx_id_4">C</a> <a class="link" href="s01.html#idx_id_5">D</a> <a class="link" href="s01.html#idx_id_6">E</a> <a class="link" href="s01.html#idx_id_7">F</a> <a class="link" href="s01.html#idx_id_8">G</a> <a class="link" href="s01.html#idx_id_9">H</a> <a class="link" href="s01.html#idx_id_10">I</a> <a class="link" href="s01.html#idx_id_11">J</a> <a class="link" href="s01.html#idx_id_12">K</a> <a class="link" href="s01.html#idx_id_13">L</a> <a class="link" href="s01.html#idx_id_14">M</a> <a class="link" href="s01.html#idx_id_15">N</a> <a class="link" href="s01.html#idx_id_16">O</a> <a class="link" href="s01.html#idx_id_17">P</a> <a class="link" href="s01.html#idx_id_18">Q</a> <a class="link" href="s01.html#idx_id_19">R</a> <a class="link" href="s01.html#idx_id_20">S</a> <a class="link" href="s01.html#idx_id_21">T</a> <a class="link" href="s01.html#idx_id_22">U</a> <a class="link" href="s01.html#idx_id_23">V</a> <a class="link" href="s01.html#idx_id_24">W</a> <a class="link" href="s01.html#idx_id_25">X</a> <a class="link" href="s01.html#idx_id_26">Y</a> <a class="link" href="s01.html#idx_id_27">Z</a></p>
<div class="variablelist"><dl class="variablelist">
<dt>

2
doc/html/indexes/s02.html
View File

@@ -24,7 +24,7 @@
</div>
<div class="section">
<div class="titlepage"><div><div><h2 class="title" style="clear: both">
<a name="id1893898"></a>Class Index</h2></div></div></div>
<a name="id1880806"></a>Class Index</h2></div></div></div>
<p><a class="link" href="s02.html#idx_id_30">A</a> <a class="link" href="s02.html#idx_id_31">B</a> <a class="link" href="s02.html#idx_id_32">C</a> <a class="link" href="s02.html#idx_id_33">D</a> <a class="link" href="s02.html#idx_id_34">E</a> <a class="link" href="s02.html#idx_id_35">F</a> <a class="link" href="s02.html#idx_id_36">G</a> <a class="link" href="s02.html#idx_id_37">H</a> <a class="link" href="s02.html#idx_id_38">I</a> <a class="link" href="s02.html#idx_id_41">L</a> <a class="link" href="s02.html#idx_id_42">M</a> <a class="link" href="s02.html#idx_id_43">N</a> <a class="link" href="s02.html#idx_id_44">O</a> <a class="link" href="s02.html#idx_id_45">P</a> <a class="link" href="s02.html#idx_id_46">Q</a> <a class="link" href="s02.html#idx_id_47">R</a> <a class="link" href="s02.html#idx_id_48">S</a> <a class="link" href="s02.html#idx_id_49">T</a> <a class="link" href="s02.html#idx_id_50">U</a> <a class="link" href="s02.html#idx_id_52">W</a></p>
<div class="variablelist"><dl class="variablelist">
<dt>

2
doc/html/indexes/s03.html
View File

@@ -24,7 +24,7 @@
</div>
<div class="section">
<div class="titlepage"><div><div><h2 class="title" style="clear: both">
<a name="id1892844"></a>Typedef Index</h2></div></div></div>
<a name="id1881997"></a>Typedef Index</h2></div></div></div>
<p><a class="link" href="s03.html#idx_id_58">A</a> <a class="link" href="s03.html#idx_id_59">B</a> <a class="link" href="s03.html#idx_id_60">C</a> <a class="link" href="s03.html#idx_id_61">D</a> <a class="link" href="s03.html#idx_id_62">E</a> <a class="link" href="s03.html#idx_id_63">F</a> <a class="link" href="s03.html#idx_id_64">G</a> <a class="link" href="s03.html#idx_id_65">H</a> <a class="link" href="s03.html#idx_id_66">I</a> <a class="link" href="s03.html#idx_id_69">L</a> <a class="link" href="s03.html#idx_id_71">N</a> <a class="link" href="s03.html#idx_id_72">O</a> <a class="link" href="s03.html#idx_id_73">P</a> <a class="link" href="s03.html#idx_id_75">R</a> <a class="link" href="s03.html#idx_id_76">S</a> <a class="link" href="s03.html#idx_id_77">T</a> <a class="link" href="s03.html#idx_id_78">U</a> <a class="link" href="s03.html#idx_id_79">V</a> <a class="link" href="s03.html#idx_id_80">W</a></p>
<div class="variablelist"><dl class="variablelist">
<dt>

2
doc/html/indexes/s04.html
View File

@@ -24,7 +24,7 @@
</div>
<div class="section">
<div class="titlepage"><div><div><h2 class="title" style="clear: both">
<a name="id1897864"></a>Macro Index</h2></div></div></div>
<a name="id1884456"></a>Macro Index</h2></div></div></div>
<p><a class="link" href="s04.html#idx_id_87">B</a> <a class="link" href="s04.html#idx_id_91">F</a></p>
<div class="variablelist"><dl class="variablelist">
<dt>

2
doc/html/indexes/s05.html
View File

@@ -23,7 +23,7 @@
</div>
<div class="section">
<div class="titlepage"><div><div><h2 class="title" style="clear: both">
<a name="id1900030"></a>Index</h2></div></div></div>
<a name="id1889054"></a>Index</h2></div></div></div>
<p><a class="link" href="s05.html#idx_id_112">2</a> <a class="link" href="s05.html#idx_id_113">4</a> <a class="link" href="s05.html#idx_id_114">A</a> <a class="link" href="s05.html#idx_id_115">B</a> <a class="link" href="s05.html#idx_id_116">C</a> <a class="link" href="s05.html#idx_id_117">D</a> <a class="link" href="s05.html#idx_id_118">E</a> <a class="link" href="s05.html#idx_id_119">F</a> <a class="link" href="s05.html#idx_id_120">G</a> <a class="link" href="s05.html#idx_id_121">H</a> <a class="link" href="s05.html#idx_id_122">I</a> <a class="link" href="s05.html#idx_id_123">J</a> <a class="link" href="s05.html#idx_id_124">K</a> <a class="link" href="s05.html#idx_id_125">L</a> <a class="link" href="s05.html#idx_id_126">M</a> <a class="link" href="s05.html#idx_id_127">N</a> <a class="link" href="s05.html#idx_id_128">O</a> <a class="link" href="s05.html#idx_id_129">P</a> <a class="link" href="s05.html#idx_id_130">Q</a> <a class="link" href="s05.html#idx_id_131">R</a> <a class="link" href="s05.html#idx_id_132">S</a> <a class="link" href="s05.html#idx_id_133">T</a> <a class="link" href="s05.html#idx_id_134">U</a> <a class="link" href="s05.html#idx_id_135">V</a> <a class="link" href="s05.html#idx_id_136">W</a> <a class="link" href="s05.html#idx_id_137">X</a> <a class="link" href="s05.html#idx_id_138">Y</a> <a class="link" href="s05.html#idx_id_139">Z</a></p>
<div class="variablelist"><dl class="variablelist">
<dt>

2
doc/html/math_toolkit/conventions.html
View File

@@ -27,7 +27,7 @@
<a name="math_toolkit.conventions"></a><a class="link" href="conventions.html" title="Document Conventions">Document Conventions</a>
</h2></div></div></div>
<p>
<a class="indexterm" name="id1001609"></a>
<a class="indexterm" name="id985976"></a>
</p>
<p>
This documentation aims to use of the following naming and formatting conventions.

2
doc/html/math_toolkit/navigation.html
View File

@@ -27,7 +27,7 @@
<a name="math_toolkit.navigation"></a><a class="link" href="navigation.html" title="Navigation">Navigation</a>
</h2></div></div></div>
<p>
<a class="indexterm" name="id1001465"></a>
<a class="indexterm" name="id985810"></a>
</p>
<p>
Boost.Math documentation is provided in both HTML and PDF formats.

106
doc/sf/lambert_w.qbk
View File

@@ -37,7 +37,7 @@ It is also called the Omega function, the inverse function of ['f(W) = W e[super
On the z-interval \[0, [infin]\] there is just one real solution.
On the interval (-1/e, 0) there are two real solutions on two branches called variously W0, W+, W-1, Wp, Wm, W-.
Two principal branches functions called `lambert_w0` and `lambert_wm1`
Two principal branches functions called `lambert_w0` and `lambert_wm1`, and their 1st derivative,
are provided by Boost's real-only implementation.
In the graphs below, the red line is the W0 branch, and the blue line the W-1 branch.
@@ -126,7 +126,7 @@ Note the expected zeros for all places up to 50 - and the correct Lambert W resu
(It is just as easy to compute even higher precisions,
at least to thousands of decimal digits, but not shown here for brevity.
See [@../../example/lambert_w_simple_examples.cpp lambert_w_simple_examples.cpp]
).
for comparison of an evaluation at 1000 decimal digit precision with __WolframAlpha).
Policies can be used to control what action to take on errors:
@@ -205,42 +205,82 @@ The symbolic algebra program __Maple also computes Lambert W to an arbitrary pre
[h4:precision Controlling the compromise between Precision and Speed]
This implementation provides good precision and excellent speed.
[h5:small_floats Floating-point types `double` and `float`]
This implementation provides good precision and excellent speed for __fundamental types `float` and `double`.
All the functions usually return values within a few __ulp
for the floating-point type, except for very small arguments very near zero,
and for arguments very close to the singularity at the branch point.
By default, this implementation provides the best possible speed.
Very slightly higher precision and less bias can be obtained by adding a __halley step refinement,
Very slightly average higher precision and less bias might be obtained
by adding a __halley step refinement,
but at the cost of more than doubling the run time.
The 'best' evaluation (the nearest __representable) can be achieved by `static_cast`ing from a
higher precision type, typically a __multiprecision type like `cpp_bin_float_50`,
at the cost of increasing run-time >100-fold;
this has been used here to provide reference values for testing.
[import ../../example/lambert_w_precision_example.cpp]
[lambert_w_precision_1]
The output for ['z = 1.F] is:
[lambert_w_precision_output_1]
showing no improvement from the Halley refinement, but for ['z = 10.F]
we need a Halley refinement (or few) to get the last bit or two correct.
[lambert_w_precision_output_1a]
[h5:big_floats Floating-point types larger than double]
For floating-point types with precision greater than `double` and `float` __fundamental_types,
a `double` evaluation is used as a first approximation followed by Halley refinement,
using a single step where it can be predicted that this will be sufficient,
only using __Halley iteration when necessary.
and only using __halley iteration when necessary.
Higher precision types are always going to be [*very, very much slower].
The 'best' evaluation (the nearest __representable) can be achieved by `static_cast`ing from a
higher precision type, typically a __multiprecision type like `cpp_bin_float_50`,
but at the cost of increasing run-time >100-fold;
this has been used here to provide our reference values for testing.
[import ../../example/lambert_w_precision_example.cpp]
For example, we get a reference value using a high precision type, for example;
`using boost::multiprecision::cpp_bin_float_50;`
that uses Halley iteration to refine until it is as precise as possible for this `cpp_bin_float_50` type.
As a further check we can compare this with a __WolframAlpha computation
using command [^N\[ProductLog\[10.\], 50\]] to get 50 decimal digits
and similarly [^N\[ProductLog\[10.\], 17\]] to get the nearest representable for 64-bit `double` precision.
[lambert_w_precision_reference_w]
giving us the same nearest representable using 64-bit `double` as `1.7455280027406994`.
However, the rational polynomial and Fukushima Schroder approximations are so good for type `float` and `double`
that negligible improvement is gained from a `double` Halley step.
This is shown with [@../../example/lambert_w_precision_example.cpp lambert_w_precision_example.cpp]
for Lambert W0:
[lambert_w_precision_1]
with this output:
[lambert_w_precision_output_1]
and then for W-1:
[lambert_w_precision_2]
with this output:
[lambert_w_precision_output_2]
[h5:differences_distribution Distribution of differences from 'best' `double` evaluations]
The distribution of differences from 'best' are shown in these graphs comparing `double` precision evaluations
with reference 'best' z50 evaluations using `cpp_bin_float_50` type reduced to `double` with `static_cast<double(z50)` :
[graph lambert_w_errors_graph]
[graph lambert_wm1_errors_graph]
As noted in the implementation section, the distribution of these differences is somewhat biased
for Lambert W-1 and this might be reduced using a `double` Halley step at small runtime cost.
But if you are serisouly concerned to get really precise computations,
the only way is using a higher precision type and then reduce to the desired type.
Fortunately, __multiprecision makes this very easy to program, if much slower.
[h4:edge_cases Edge and Corner cases]
[h5:w0_edges w0 Branch]
@@ -430,19 +470,16 @@ For example, table below a test of Lambert W0 10000 values of argument z coverin
[[after Halley step] [ 9710 ] [ 288 ] [ 2 ] [0] [ 292 ] [22]]
] [/table:lambert_w0_Fukushima W0]
Lambert W0 values computed using the Fukushima method with Schroeder refinement gave about 1/6 lambert_w0 values that are
one bit different from the 'best', and < 1% that are a few bits 'wrong'.
If a Halley refinement step is added, only 1 in 30 are even one bit different, and only 2 two bits 'wrong'.
[table:lambert_w0_plus_halley Rational polynomial Lambert W0 and typical improvement from a single Halley step.
[[Method] [Exact] [One_bit] [Two_bits] [Few_bits] [inexact] [bias]]
[[rational/polynomial] [7135] [ 2863 ] [ 2 ] [0] [ 2867 ] [ -59] ]
[[after Halley step ] [ 9724 ] [273] [ 3 ] [0] [ 279 ] [5] ]
] [/table:lambert_w0_plus_halley]
With the rational polynomial approximation method, there are a third one bit from the best and none more than two bits.
Adding a Halley step (or iteration) reduces the number that are one bit different from about a third down to one in 30;
this is unavoidable 'computational noise'.
@@ -696,7 +733,22 @@ precision can be much lower, as might be expected.
[@../../example/lambert_w_simple_examples.cpp lambert_w_simple_examples.cpp]
[@../../test/test_lambert_w.cpp test_lambert_w.cpp] contains routine tests using __boost_test.
Two macros provide optional tests if #defined before BOOST_TEST_MAIN.
[h5:quadrature_testing Testing with quadrature]
A further method of testing over a wide range of argument z values was devised by Nick Thompson
(partly to test the recently added quadrature routines).
These are definite integral formulas involving the W function that are exactly known constants,
for example, LambertW0(1/(z[sup2]) == [sqrt](2[pi]),
see [@https://en.wikipedia.org/wiki/Lambert_W_function#Definite_integrals Definite Integrals].
Some care was needed to avoid overflow and underflow as the function must evaluate to a finite result over the entire range.
These are not routinely evaluated by the Boost testers,
but can be run by defining a macro BOOST_MATH_LAMBERT_W_INTEGRALS.
Three macros provide optional tests if #defined before BOOST_TEST_MAIN.
* BOOST_MATH_TEST_MULTIPRECISION Adds test for __multiprecision types.

138
example/lambert_w_precision_example.cpp
View File

@@ -14,6 +14,8 @@
#include <boost/exception/exception.hpp> // boost::exception
#include <boost/math/constants/constants.hpp> // For exp_minus_one == 3.67879441171442321595523770161460867e-01.
#include <boost/math/policies/policy.hpp>
#include <boost/math/special_functions/next.hpp> // for float_distance.
#include <boost/math/special_functions/relative_difference.hpp> // for relative and epsilon difference.
// Built-in/fundamental GCC float128 or Intel Quad 128-bit type, if available.
#ifdef __GNUC__
@@ -64,47 +66,144 @@ int main()
using boost::math::policies::domain_error;
using boost::math::policies::throw_on_error;
//[lambert_w_precision_reference_w
using boost::multiprecision::cpp_bin_float_50;
using boost::math::float_distance;
cpp_bin_float_50 z("10."); // Note use a decimal digit string, not a double 10.
cpp_bin_float_50 r;
std::cout.precision(std::numeric_limits<cpp_bin_float_50>::digits10);
r = lambert_w0(z); // Default policy.
std::cout << "lambert_w0(z) cpp_bin_float_50 = " << r << std::endl;
//lambert_w0(z) cpp_bin_float_50 = 1.7455280027406993830743012648753899115352881290809
// [N[productlog[10], 50]] == 1.7455280027406993830743012648753899115352881290809
std::cout.precision(std::numeric_limits<double>::max_digits10);
std::cout << "lambert_w0(z) static_cast from cpp_bin_float_50 = "
<< static_cast<double>(r) << std::endl;
// double lambert_w0(z) static_cast from cpp_bin_float_50 = 1.7455280027406994
// [N[productlog[10], 17]] == 1.7455280027406994
std::cout << "bits different from Wolfram = "
<< static_cast<int>(float_distance(static_cast<double>(r), 1.7455280027406994))
<< std::endl; // 0
//] [/lambert_w_precision_reference_w]
//[lambert_w_precision_0
std::cout.precision(std::numeric_limits<float>::max_digits10); // Show all potentially significant decimal digits,
std::cout << std::showpoint << std::endl; // and show any significant trailing zeros too.
double x = 10.;
float x = 10.;
std::cout << "Lambert W (" << x << ") = " << lambert_w0(x) << std::endl;
//
//] [/lambert_w_precision_0]
/*
//[lambert_w_precision_output_0
Lambert W (10.0000000) = 1.74552800
//] [/lambert_w_precision_output_0]
*/
{ // Lambert W0 Halley step example
//[lambert_w_precision_1
using boost::math::lambert_w_detail::lambert_w_halley_step;
using boost::math::epsilon_difference;
using boost::math::relative_difference;
std::cout << std::showpoint << std::endl; // and show any significant trailing zeros too.
std::cout.precision(std::numeric_limits<double>::max_digits10); // 17 decimal digits for double.
cpp_bin_float_50 z50("1.23"); // Note: use a decimal digit string, not a double 1.23!
double z = static_cast<double>(z50);
cpp_bin_float_50 w50;
w50 = lambert_w0(z50);
std::cout.precision(std::numeric_limits<cpp_bin_float_50>::max_digits10); // 50 decimal digits.
std::cout << "Reference Lambert W (" << z << ") =\n "
<< w50 << std::endl;
std::cout.precision(std::numeric_limits<double>::max_digits10); // 17 decimal digits for double.
double wr = static_cast<double>(w50);
std::cout << "Reference Lambert W (" << z << ") = " << wr << std::endl;
double z = 10.;
double w = lambert_w0(z);
std::cout << "Lambert W (" << x << ") = " << lambert_w0(z) << std::endl;
std::cout << "Rat/poly Lambert W (" << z << ") = " << lambert_w0(z) << std::endl;
// Add a Halley step to the value obtained from rational polynomial approximation.
double ww = lambert_w_halley_step(lambert_w0(z), z);
std::cout << "Lambert W (" << x << ") = " << lambert_w_halley_step(lambert_w0(z), z) << std::endl;
std::cout << "Halley Step Lambert W (" << z << ") = " << lambert_w_halley_step(lambert_w0(z), z) << std::endl;
std::cout << "difference from Halley step " << w - ww << std::endl;
std::cout << "absolute difference from Halley step = " << w - ww << std::endl;
std::cout << "relative difference from Halley step = " << relative_difference(w, ww) << std::endl;
std::cout << "epsilon difference from Halley step = " << epsilon_difference(w, ww) << std::endl;
std::cout << "epsilon for float = " << std::numeric_limits<double>::epsilon() << std::endl;
std::cout << "bits different from Halley step = " << static_cast<int>(float_distance(w, ww)) << std::endl;
//] [/lambert_w_precision_1]
/*
/*
//[lambert_w_precision_output_1
Reference Lambert W (1.2299999999999999822364316059974953532218933105468750) =
0.64520356959320237759035605255334853830173300262666480
Reference Lambert W (1.2300000000000000) = 0.64520356959320235
Rat/poly Lambert W (1.2300000000000000) = 0.64520356959320224
Halley Step Lambert W (1.2300000000000000) = 0.64520356959320235
absolute difference from Halley step = -1.1102230246251565e-16
relative difference from Halley step = 1.7207329236029286e-16
epsilon difference from Halley step = 0.77494921535422934
epsilon for float = 2.2204460492503131e-16
bits different from Halley step = 1
//] [/lambert_w_precision_output_1]
*/
} // Lambert W0 Halley step example
{ // Lambert W-1 Halley step example
//[lambert_w_precision_2
using boost::math::lambert_w_detail::lambert_w_halley_step;
using boost::math::epsilon_difference;
using boost::math::relative_difference;
std::cout << std::showpoint << std::endl; // and show any significant trailing zeros too.
std::cout.precision(std::numeric_limits<double>::max_digits10); // 17 decimal digits for double.
cpp_bin_float_50 z50("-0.123"); // Note: use a decimal digit string, not a double -1.234!
double z = static_cast<double>(z50);
cpp_bin_float_50 wm1_50;
wm1_50 = lambert_wm1(z50);
std::cout.precision(std::numeric_limits<cpp_bin_float_50>::max_digits10); // 50 decimal digits.
std::cout << "Reference Lambert W-1 (" << z << ") =\n "
<< wm1_50 << std::endl;
std::cout.precision(std::numeric_limits<double>::max_digits10); // 17 decimal digits for double.
double wr = static_cast<double>(wm1_50);
std::cout << "Reference Lambert W-1 (" << z << ") = " << wr << std::endl;
double w = lambert_wm1(z);
std::cout << "Rat/poly Lambert W-1 (" << z << ") = " << lambert_wm1(z) << std::endl;
// Add a Halley step to the value obtained from rational polynomial approximation.
double ww = lambert_w_halley_step(lambert_wm1(z), z);
std::cout << "Halley Step Lambert W (" << z << ") = " << lambert_w_halley_step(lambert_wm1(z), z) << std::endl;
std::cout << "absolute difference from Halley step = " << w - ww << std::endl;
std::cout << "relative difference from Halley step = " << relative_difference(w, ww) << std::endl;
std::cout << "epsilon difference from Halley step = " << epsilon_difference(w, ww) << std::endl;
std::cout << "epsilon for float = " << std::numeric_limits<double>::epsilon() << std::endl;
std::cout << "bits different from Halley step = " << static_cast<int>(float_distance(w, ww)) << std::endl;
//] [/lambert_w_precision_2]
}
/*
//[lambert_w_precision_output_2
Reference Lambert W-1 (-0.12299999999999999822364316059974953532218933105468750) =
-3.2849102557740360179084675531714935199110302996513384
Reference Lambert W-1 (-0.12300000000000000) = -3.2849102557740362
Rat/poly Lambert W-1 (-0.12300000000000000) = -3.2849102557740357
Halley Step Lambert W (-0.12300000000000000) = -3.2849102557740362
absolute difference from Halley step = 4.4408920985006262e-16
relative difference from Halley step = 1.3519066740696092e-16
epsilon difference from Halley step = 0.60884463935795785
epsilon for float = 2.2204460492503131e-16
bits different from Halley step = -1
//] [/lambert_w_precision_output_2]
*/
// z = 10.F needs Halley to get the last bit or two correct.
//[lambert_w_precision_output_1a
//] [/lambert_w_precision_output_1a]
*/
// Similar example using cpp_bin_float_quad (128-bit floating-point types).
@@ -126,6 +225,7 @@ int main()
cpp_dec_float_50 z("10");
cpp_dec_float_50 r;
std::cout.precision(std::numeric_limits<cpp_dec_float_50>::digits10);
r = lambert_w0(z); // Default policy.
std::cout << "lambert_w0(z) cpp_dec_float_50 = " << r << std::endl; // 0.56714329040978387299996866221035554975381578718651
std::cout.precision(std::numeric_limits<cpp_bin_float_quad>::max_digits10);
@@ -152,8 +252,6 @@ int main()
// lambert_w0(z) cpp_dec_float_50 cast to quad = 1.745528002740699383074301264875389837
// lambert_w0(z) cpp_dec_float_50 cast to quad = 1.74552800274069938307430126487539
}
// std::cout << "Lambert W (" << xd << ") = " << lambert_w0(xd) << std::endl; //
}
catch (std::exception& ex)
{