These files are a collection of the most basic mathematical constants for
C/C++ programs.
These are like to be useful even by programs just computing areas of circles!
Mathematical constants (like pi) are most unlikely to change, whereas
physical ‘constants‘ (like h) are certain to vary.
Physical constants usually have a known uncertainty, usually greater than even
the accuracy of float representations, whereas mathematical constants can be
calculated to an arbitrarily high precision.
Therefore these two types will never be stored in the same file to avoid
confusion.
A special feature of this collection is that they are a little more accurate than can be represented by current and foreseeable floating-point hardware. The collection is convenient in that it does not require a search of reference sources and is useful because it ensures that the values are as accurate as possible and as portable as possible between different systems. Small differences in constants can lead to puzzling differences between computations. The values are slightly, but significantly, more accurate than is usually achieved by using the compiler to calculate constants.
The constants have all been calculated using software working with precision unlimited (except by compute time). The precision selected (40 decimal digits) exceeds the accuracy of foreseeable floating-point hardware.
The objective is to achieve the full accuracy possible for all real-life computations. This has no extra cost to the user, but reduces irritating, and often confusing and hard to trace effects, caused by the inevitable limited precision of floating-point calculations. At least these show as a spurious least-significant digit; at worst slightly inaccurate constants cause algorithms to fail because internal comparisons just fail.
Different applications need constants in different formats. Five types of file are provided:
1. Files which contain C Macros which #define and #undef constants.
2. Files that contain const float constants. Useful if only float precision is required.
3. Files that contain const double constants. Useful for most purposes.
4. Files that contain const long double constants. Useful where long double precision is needed.
5. Files that contain class constant template function constants. Useful where more than one precision is needed and to facilitate optimisation by compilers.
All files contain the same constants and identical numeric values.
To avoid including an excessive number of unwanted constants, the constants are subdivided into several files, the exact distribution and content yet to be agreed.
The file define_constants.h consists of a list of C style macros that define a large collection of mathematical constants very accurately, for example:
#define BOOST_PI 3.14159265358979323846264338327950288L /* pi */
Individual MACRO value(s) can be used directly, for example:
const double PI = BOOST_PI;
or with casting to other types, for example:
const long double PI = (long double)BOOST_PI;
Unlike
most header files, the file can be used in many ways in C and C++, and indeed
many other languages.
(Note that is it deliberately a .h file and not the Boost convention of a .hpp
file. This is to allow it to be used in
C programs as well as C++ programs).
Corresponding files that contain matching #undef statements called undefine*.h are also provided to reduce namespace pollution and the risk of name clashes or confusion.
Files containing const float, const double and const long double values can be #included or
Individual literal value(s) strings can be copied and pasted into programs, for example:
const float pi = 3.14159265358979323846264338327950288F;
(Note that the suffix is F to avoid the need to cast long double to float).
const long double pi = 3.14159265358979323846264338327950288L;
(Note the suffix is L because the literal should be a long double).
(Note that to avoid any loss of precision for long double definitions, the MACRO constants are defined as long double by the suffix L, but this means that casting is needed or desirable to avoid errors or warnings for less precise types double or float.)
Const double pi = (double)BOOST_PI; /* C style cast */
Or
C++ static cast down to double and/or float precisions.
const long double pi = static_cast<double>(BOOST_PI); //
C++ static cast
const long float piF =
static_cast<float>(BOOST_PI); // C++ static cast
To document helpfully, it is sensible to static_cast all definitions, for example:
const long double piF = static_cast<long double>(BOOST_PI); // C++ static cast
(Note one needs two different names, pi and piF, if different precisions must coexist (unless one uses constant class template functions described below).
3 The header file constants.h can be included in a separate translation unit:
#include “constants.h”extern const long double pi = static_cast<double>(BOOST_PI); // C++ static cast
extern const long double e = static_cast<double>(BOOST_E); // C++ static cast
will make only pi and e visible to other translation units through the linker. The names can be chosen to avoid name clashes.
4 The simpler scheme defining constants within a translation unit
const long double pi = static_cast<double>(BOOST_PI); // C++ static cast
may cause some name clashes, so a matching file undef_constants.h is provided to undefine ALL the macros defined by define_constants.h. This undefined file should be included immediately after the last reference to any BOOST macro constant.
5 When multiple precision constants are needed in the same program, (and namespace pollution is to be avoided) class constant template functions in files like function_constants.hpp can be used:
// Select an entire set of constants with a single using directive:
using namespace boost::math::double_constants;
// Or float constants...
// using namespace boost::math::float_constants;
// So naive users (mea culpa) can just say "pi",
// and get boost::math::double_constants::pi.
cout << "pi = " << pi << "\n";
double r = 2.0;
cout << "area = " << pi * r * r << "\n";
// Naughty users CAN'T write pi = 3.; !!! // but they get a cryptic error message:
// C2678: binary '=' : no operator defined which takes a left-hand operand of type
// 'const struct boost::math::constant<struct boost::math::pi_tag,double>'
// Or they can explicitly access fully specified constants at any time:
cout.precision(std::numeric_limits<float>::digits10 + 2);
// Ensure all significant digits and two noisy digits are displayed.
cout << "float pi = " << boost::math::float_constants::pi << "\n";
cout.precision(std::numeric_limits<long double>::digits10 + 2);
cout << "long double pi = " << boost::math::long_double_constants::pi << "\n";
using namespace boost::math;
cout.precision(std::numeric_limits<float>::digits10 + 2);
cout << "float pi = " << float_constants::pi << "\n";
cout.precision(std::numeric_limits<double>::digits10 + 2);
cout << "double pi = " << double_constants::pi << "\n";
// But CANNOT switch to float constants.
// using namespace boost::math::float_constants;
// cout << "pi = " << pi << "\n"; // error C2872: 'pi' : ambiguous symbol
// This is probably a useful safeguard.
// For types that are not specialised, one could write:
const short piI = (short)boost::math::float_constants::pi;
cout << "short pi = " << piI << '\n'; // 3!
// Access via the templated function call may allow better optimisation?
float pifloat = boost::math::constant<pi_tag, float >();
// If a template is not specialised, then fails at link time.
This scheme was devised by Michael Kenniston (after much discussion on Boost.org) for the following reasons:
· Only
need to write one interface (not counting macros-for-C), which is really a
separate library from the C++ stuff)
· No
macros: all names are in namespaces.
· Naive
users can just say "pi" instead of "pi()".
· Average
users can switch "pi" between float and double by changing one
"using" directive.
· No
"double rounding" (i.e. no multiple conversions).
· Platform
implementers can use whatever evil incantations they desire for maximum
accuracy.
· Everything
is implemented with inline functions to ease constant folding.
· It
is (relatively) easy to extend with new constants and representation types.
· Purists
can access all constants via templated functions.
· No
"static-initialization-order" issues, because the static variables contain
no data to initialise.
· No
partial template specialization, so it works on MSVC.
· No
non-type template parameters, so it should work on Borland.
Same basic technique will probably work for converting read access to any statically
allocated class object to an inline function call, including in particular
physical constants with units/dimensions.
6 C++ allows other methods of avoiding name collisions, for example using named and unnamed namespaces.
The advantages of pre-computed constants are in accuracy and precision, reliability, and most important portability.
It may be worth using single constants rather than using the preprocessor and/or compiler to calculate a derived value. For example, it is a little better to use twoPi rather than pi + pi or 2. * p. This is because the single constants are produced with a higher precision than the floating-point hardware and/or software that the compiler is also certain to use for its computation. (The differences are especially significant for trig and pow function!)
The constants are entirely fixed at compile time, so there is no code to call functions, for example, log, nor the need to include header files which define those functions, for example <cmath>. This may speed compilation and reduce file accesses and dependency.
There is no evaluation of the function at program start time, nor wasted evaluation if the constant is never used, so programs will be a bit quicker.
Portability is improved because there is no compile or run-time evaluation of constants, there can be no difference caused by different versions of functions, for example in different version of C math libraries, nor in different version of hardware floating point evaluation. Although in reality, this is rather unlikely to be a problem, the effort and time wasted proving that it is not a problem is likely to be significant.
The most potentially serious loss of accuracy occurs with transcendental functions, most serious of all, the exponentiation function x raised to power y, pow(x, y) for x**y or x^y where y may not be an integer.
If the C program float.h is available, the macro FLT_MANT_DIG may provide the number of digits (usually radix 2) in the mantissa (significand) that limits precision (but not the maximum or minimum values that are limited mainly by maximum exponent).
/* The number of base FLT_RADIX digits in the mantissa. */#define _DBL_RADIX 2 /* exponent radix */
#define DBL_MANT_DIG 53 /* # of bits in mantissa */
#define LDBL_MANT_DIG 64 /* # of bits in mantissa */
The same values are available using Std C++ <limits>, for example as:
int double_radix = std::numeric_limits<long double>::radix; // usually 2 for IEEE FP
int double_mantissa_digits = std::numeric_limits<long double>::digits; // Usually 53 for 64-bit double.
int double_digits10 = std::numeric_limits<long double>::digits10; // 15 for 64-bit double.
(But the 128-bit representation used by the AIX OS warns that that this value, defined as 106-bit for AIX long double, may be not tell the whole story – the actual number of bits of precision can vary.)
To ensure no loss of accuracy when reading in, all floating point constants should be specified to TWO decimal digits MORE than the value of ISO Standard 14882:1998 C++ Library limits for std::numeric_limits<type>:: . For example std::numeric_limits<double>::digits10 (equivalent to Standard C Library float.h DBL_DIG) is 15 for double using IEEE 754/IEC559 standard 64-bit real, 11-bit exponent, 52+1 = 53 bit significand (mantissa) using Motorola 680X0 or Intel 8087 floating point coprocessor as incorporated in all Pentium chips. So double constants must provide at least 17 decimal digits.
All the constants here should meet this criterion, and all also meet C++ long double extended precision if the compiler uses IEEE 80-bit real, 15-bit exponent, 64-bit significand (or mantissa) format for which constants will need 21 decimal digits, and all will even meet the requirements for 128-bit real arithmetic providing significand (mantissa) accuracy of 106-bit used by the Sun Sharc processor, the AIX OS, and which can be emulated in software using, for example C++ NTL quad_float class by Victor Shoup (derived from the doubledouble class of Keith Briggs).
However, if the Intel Itanium IA-64 128-bit quad_precision (112-bit significand, 15 exponent, and sign) is (as planned) implemented in hardware, slightly more precise than 36 decimal digit constants might be required. The number of decimal digits chosen has been restricted to avoid overflow of non-compliant compilers (which should accept, but ignore, unlimited numbers of digits). (See www.intel.com , ADAG.pdf for IA-64 Application Developers Architecture Guide).
Floating-point hardware is unlikely to provide higher precision for some time to come. (For applications that require a higher precision, software solutions are available - see references below).
X86 floating-point hardware includes instructions for providing floating-point constants pi and log base 2 (10). These are accurate to 64, including guard and sticky bits, which is higher than for any constant read into a register or memory by the compiler, even from a more precise constant decimal digit string.
Internally, the Intel X86 floating-point hardware holds all numbers in extended precision real format (64 significand bits compared to external double of 53 significand bits). Extra bits are mainly to support rounding mode correctly. Calculations performed using floating point register to register operations will be performed with this 64-bit accuracy. But operations that use external memory will have lower accuracy. The compiler makes the choice between register and memory, so expecting completely portable accuracy is probably too optimistic; however using accurate literal constants probably gives the best chance of portability (but not necessarily the best possible accuracy).
Alas the complexity of floating point calculations is such that calculations cannot be absolutely guaranteed to give the same result, even if both use hardware that conforms to the IEEE 754/IEC559 floating-point specification. This is most likely to be a significant problem when iterative calculations at the limit of accuracy are involved, perhaps leading to failure to converge, or finding a different solution. These differences are most troublesome when trying to prove the equivalence of difference systems for computer validation. Because one can predict that differences may arise, one should not pessimistically assume that different output means that either systems is faulty, nor assume that the systems are not equivalent.
One method of ensuring that the IEEE extended double precision floating-point (typically the Intel X86/Pentium) is not used is to alter the appropriate bit in the processor “control word”. This can only be done in assembler. Extended precision significand is 64-bit, double is only 53-bit, and memory values are always 53-bit (plus sign and exponent to total 64 bits). (See www.intel.com 27064004.pdf for 80C187 80-bit Math Coprocessor datasheet).
The default value of the code word bit is NOT extended for Microsoft Windows, but IS extended for gcc/Linux.
The third way to fix the problem is to 'force' all intermediate floating point results into memory. This is not an 'ideal' fix, since it is not fully equivalent to 53-bit precision (because of double rounding), but it is reported to work (although a full proof of correctness in this case is not available). See X86ControlWord.cpp for an example.
NTL's quad_float code (used to generate these constants) ensures use of 53-bit significands on all platforms by storing intermediate results in local variables declared to be 'volatile'. This is the solution to the problem that NTL uses if it detects the problem and can't fix it using the GNU 'asm' hack mentioned above. This solution should work on any platform that faithfully implements 'volatile' according to the ANSI C standard. See www.shoup.net.
Note that MSVC++ compiler up to version 6 SP 5 implements long double exactly as double, IEEE 64-bit floating-point format. (The compiler does NOT use IEEE 80-bit extended precision available within the Intel Pentium or X86 floating point co-processor processor family instruction set).
Pre-processor should calculate any long integers using 36 bits (P J Plauger, The Standard C Library, p 78) but none of the constants in constants.h are currently of integer type. (Some are floating point versions of integer values, for example: unity, zero, two, three …).
For most values, Victor Shoup's NTL has been used to compute (and/or cross-check) the values in this file with a significand of 150 bits (about 50 decimal digits) and value are output as 40 decimal digit strings. There is no typing or editing to introduce errors.
The
numeric string values in the file constants.h are checked either one or both of
two ways:
1 Checks that the values are consistent with values calculated by the compiler
preprocessor.
2 Comparison with values in references tables: for example in D. E. Knuth,
Appendix A (scanned in from the printed book and converted to a text file), or
M. Abramowitz & Irene A Stegun.
Some values are be tested with straightforward checks against compiler preprocessor calculated values, for example comparing:
if (twoPi != 2. * pi) // may have a problem?
or C macro assert (which causes program termination if it fails - a more noticeable conclusion, but one that does not allow one to continue testing other constants).
assert(twoPi == 2. * pi)
but many naively formulated tests appear to fail to the extent that the two values differ by the machine least significant significand bit, so instead it is better to test the absolute relative difference between values thus:
if (fabs(twoPi - 2.L * pi)/ twoPi >=
epsilon) // Do have a problem with
accuracy!
(Epsilon is the smallest floating-point representation of the least-significant significand bit. See below.)
The equivalent using the C macro assert is:
assert(fabs(sqrtTwo * sqrtTwo - 2.)/2. <=
eps);
It is best to choose the constant to divide by to make the test relative to the constant value, rather than recomputing the check value.
Typical epsilon for type double for Intel X86 IEEE754/IEC559 hardware is 2.22e-016.
Minimum
values specified by ISO 14882:1998 section 18.2.1 are defined as equivalent to
the C Standard constants:
#define FLT_EPSILON <float rvalue <=
10^(-5)>
#define DBL_EPSILON <double
rvalue <= 10^(-9)>
#define LDBL_EPSILON <long double
rvalue <= 10^(-9)>
The macros yield the smallest X of type float, double or long
double such that 1.0 + X != 1.0.
C++
versions of these are provided by the ISO 14882:1998 Standard numeric limits
template, for example:
std::numeric_limits<float>::epsilon();
std::numeric_limits<double>::epsilon()
std::numeric_limits<long
double>::epsilon();
This is shown in the header file comments, and also includes a timestamp and the version of NTL used.
Paul
A Bristow, pbristow@hetp.u-net.com,
Oct 2001
Floating-point – real type stored as sign, exponent and significand. Unlike integer types, most floating-point values cannot be stored exactly.
Exponent –power of 2 by which to multiply the significand (implicitly binary). Always signed, using one bit, and often offset or biased for storage convenience.
Mantissa – see significand.
Significand – (previously called mantissa) binary fraction where most significant bit represents 0.5, next 0.25, then 0.125 … The most significant bit (0.5) is usually not stored and is assumed to be set, so the significand always represents a fraction greater than 0.5, (but just less than 1).
Acknowledgements:
The
following contributed to discussions on the BOOST.org forum:
K Hagan
Petr Kocmid
P J Plauger
Herve Bronnimann
John Maddock
Peter Schmittteckert
Jes Maurer
Peter Dimov
Beman Dawes
David Abrahams
Ed Bray
Matt Austern
Gabriel Dos Reis
Eric Ford
Michael Kenniston
Author: Paul A. Bristow, hetp Chromatography Quantitation,
Copyright 2001 Paul A Bristow and hetp, 16 Oct 01 15:26