Skip to content

Commit 92233d1

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 17bacf0 commit 92233d1

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>
@@ -114,6 +114,44 @@
114114
#define POV_CPP11_SUPPORTED (__cplusplus >= 201103L)
115115
#endif
116116

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

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

0 commit comments

Comments
 (0)