Skip to content

Commit 9fd5f23

Browse files
committed
Implementing suggested PRECISE_FLOAT macro mechanism for polysolve().
Primarily enabling 'long double' support for the polysolve() / 'sturm' solver given 'long double' is part of the C+11 standard and gets most users 80 bits or more as opposed to 64. Further, the additional bits are almost always hardware backed making for reasonable performance. With GNU g++ 64 bit environments '__float128' is a valid setting with no additional library required - though it's very slow. Other extended accuracy floats such as 'double double' enabled, via the PRECISE_ mechanisms, but they involve additional libraries and settings. See code. Note! With this update the polysolve sturm chain is no longer pruned with a default 1e-10 value, but rather a value appropriate for the floating point type. Means better accuracy for some previously bad root results and less accuracy for others. Coming updates will better tune the chain equation pruning.
1 parent 7fc8c63 commit 9fd5f23

File tree

3 files changed

+243
-85
lines changed

3 files changed

+243
-85
lines changed

source/base/configbase.h

Lines changed: 137 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -41,7 +41,7 @@
4141
#include "syspovconfigbase.h"
4242

4343
#include <cstdint>
44-
44+
#include <cfloat>
4545
#include <limits>
4646

4747
#include <boost/version.hpp>
@@ -113,6 +113,44 @@
113113
#define POV_CPP11_SUPPORTED (__cplusplus >= 201103L)
114114
#endif
115115

116+
/// @}
117+
///
118+
//******************************************************************************
119+
///
120+
/// @name C++11 Constant Expression Support
121+
///
122+
/// The following macros and constant expression (constexpr) functions enable the set up of
123+
/// compile time typed contants.
124+
///
125+
/// References in addtion to the C++11 standard itself:
126+
/// * [sac10-constexpr.pdf](http://www.stroustrup.com/sac10-constexpr.pdf)
127+
/// * [n2235.pdf](http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2007/n2235.pdf)
128+
///
129+
/// @{
130+
131+
/// @def CX_XSTR(A)
132+
/// Macro converts preprocessor macro string to a contant string.
133+
///
134+
/// @note
135+
/// Itself using an macro internale macro called CX_STR(A).
136+
///
137+
#define CX_XSTR(A) CX_STR(A)
138+
#define CX_STR(A) #A
139+
140+
/// @def CX_STRCMP(A,B)
141+
/// constexpr function comparing two strings in compile time constant fashion.
142+
///
143+
/// @param[in] A First string.
144+
/// @param[in] B Second string.
145+
/// @return Integer 0 on matching strings and non zero integer otherwise.
146+
///
147+
constexpr int CX_STRCMP(char const* lhs, char const* rhs)
148+
{
149+
return (('\0' == lhs[0]) && ('\0' == rhs[0])) ? 0
150+
: (lhs[0] != rhs[0]) ? (lhs[0] - rhs[0])
151+
: CX_STRCMP( lhs+1, rhs+1 );
152+
}
153+
116154
/// @}
117155
///
118156
//******************************************************************************
@@ -320,6 +358,104 @@
320358
#define SNGL float
321359
#endif
322360

361+
/// @def PRECISE_FLOAT
362+
/// The foating-point data type to use within the polysolve() solver.
363+
///
364+
/// Normally PRECISE_FLOAT will be set to @ref DBL. However, 'long double' is a
365+
/// precision guaranteed by the C++ standard to exist and to be at least of
366+
/// 'double' size. In practice, most environments provide 80bits of extended
367+
/// accuracy in a hardware backed way and some environments provide more bits.
368+
///
369+
/// Further, particular compilers and environments may offer other float
370+
/// data types. Recent GNU g++ compilers, for example, offer the __float128
371+
/// data type which coorresponds to the IEEE 754-2008 binary128 type.
372+
///
373+
/// @note
374+
/// The setting is used to set internal 'constextr int PRECISE_DIG' and
375+
/// 'constexpr DBL PRECISE_EPSILON' typed values.
376+
///
377+
/// @note
378+
/// If @ref PRECISE_FLOAT is set to 'float', 'double' or 'long double', the C++11 standard
379+
/// itself defines via cfloat include the macros: FLT_DIG, DBL_DIG, LDBL_DIG, FLT_EPSILON,
380+
/// DBL_EPSILON, LDBL_EPSILON. These macro values are used to set PRECISE_DIG and
381+
/// PRECISE_EPSILON polysolve() values via the C++ constexpr mechanism.
382+
///
383+
/// @note
384+
/// Users access the polysolve() solver via the 'sturm' keyword though it's always
385+
/// used for polynomials of order 5 or more.
386+
///
387+
/// @attention
388+
/// polysolve() math with data types not run in hardware or with awkward bit sizes with
389+
/// respect to the running computer's data bus size will run substantially slower.
390+
///
391+
#ifndef PRECISE_FLOAT
392+
#define PRECISE_FLOAT DBL
393+
#endif
394+
395+
/// @def PRECISE_FABS
396+
/// The floating point absolute value function name to use with @ref PRECISE_FLOAT.
397+
///
398+
/// Where @ref PRECISE_FLOAT a C++11 standard floating point type this should be
399+
/// the default fabs. When using, for example, __float128 with GNU g++ in a supported 64 bit
400+
/// environment the value would defined to be __builtin_fabsq - as that function is a built-in
401+
/// for 64 bit compiles. It might be set to fabsq if using GNU g++'s quadmath library.
402+
///
403+
/// @note
404+
/// At present PRECISE_FABS, like @ref PRECISE_FLOAT, applies only within polysolve().
405+
///
406+
#ifndef PRECISE_FABS
407+
#define PRECISE_FABS fabs
408+
#endif
409+
410+
/// @def PRECISE_DIGITS
411+
/// Base 10 digits where @ref PRECISE_FLOAT type not a C++11 floating point standard.
412+
///
413+
/// Defaults to an IEEE 754-2008 binary128 value of 33 with the expectation this the most likely
414+
/// non-standard @ref PRECISE_FLOAT usage, but it can be explicitly set as needed.
415+
///
416+
#ifndef PRECISE_DIGITS
417+
#define PRECISE_DIGITS 33
418+
#endif
419+
420+
/// @def PRECISE_EPSLN
421+
/// Best epsilon step where @ref PRECISE_FLOAT type not a C++11 floating point standard.
422+
///
423+
/// Defaults to an IEEE 754-2008 binary128 value of 1.9259e-34 with the expectation this the
424+
/// most likely non-standard @ref PRECISE_FLOAT type, but it can be explicitly set as needed.
425+
///
426+
#ifndef PRECISE_EPSLN
427+
#define PRECISE_EPSLN 1.92592994438724e-34
428+
#endif
429+
430+
/// @def PRECISE_DIG
431+
/// Internal 'constexpr int' value for maximum decimal digits for given @ref PRECISE_FLOAT.
432+
///
433+
/// Set to C+11 standard value where defined and to @ref PRECISE_DIGITS otherwise.
434+
///
435+
constexpr int PRECISE_FLT = CX_STRCMP(CX_XSTR(PRECISE_FLOAT), "float");
436+
constexpr int PRECISE_DBL = CX_STRCMP(CX_XSTR(PRECISE_FLOAT), "double");
437+
constexpr int PRECISE_LDBL = CX_STRCMP(CX_XSTR(PRECISE_FLOAT), "long double");
438+
439+
constexpr int PRECISE_DIG = PRECISE_FLT ?
440+
PRECISE_DBL ?
441+
PRECISE_LDBL ? PRECISE_DIGITS
442+
: LDBL_DIG
443+
: DBL_DIG
444+
: FLT_DIG;
445+
446+
/// @def PRECISE_EPSILON
447+
/// Internal 'constexpr DBL' value for minimum epsilon step for given @ref PRECISE_FLOAT.
448+
///
449+
/// Set to C+11 standard value where defined and to @ref PRECISE_EPSLN otherwise.
450+
///
451+
452+
constexpr DBL PRECISE_EPSILON = PRECISE_FLT ?
453+
PRECISE_DBL ?
454+
PRECISE_LDBL ? PRECISE_EPSLN
455+
: LDBL_EPSILON
456+
: DBL_EPSILON
457+
: FLT_EPSILON;
458+
323459
/// @}
324460
///
325461
//******************************************************************************

0 commit comments

Comments
 (0)