Skip to content

Commit

Permalink
Created new constexpr DBL variable MIN_ISECT_DEPTH_RETURNED.
Browse files Browse the repository at this point in the history 
Near term need to use a value not MIN_ISECT_DEPTH in blob.cpp to test
Jérôme's github pull request POV-Ray#358. Longer term aim is to drive all returned
intersection depths from shape code to MIN_ISECT_DEPTH_RETURNED. The value
is automatically derived from DBL setting and at double resolves to about
4.44089e-08. On the order of the square root of POV_DBL_EPSILON.

Moving blob.cpp's inside test value INSIDE_TOLERANCE to POV_DBL_EPSILON over
previous recent re-calculation. Over time plan is to move EPSILONs best
nearer a double's step to POV_DBL_EPSILON.

Cleaning up doxygen documentation added during recent solver related
updates.
  • Loading branch information
@wfpokorny
wfpokorny committed Jun 20, 2018
1 parent af78a9b commit 761dd0b
Show file tree
Hide file tree
Showing 4 changed files with 206 additions and 115 deletions.
182 changes: 125 additions & 57 deletions source/base/configbase.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,44 +113,6 @@
#define POV_CPP11_SUPPORTED (__cplusplus >= 201103L)
#endif

/// @}
///
//******************************************************************************
///
/// @name C++11 Constant Expression Support
///
/// The following macros and constant expression (constexpr) functions enable the set up of
/// compile time typed contants.
///
/// References in addtion to the C++11 standard itself:
/// * [sac10-constexpr.pdf](http://www.stroustrup.com/sac10-constexpr.pdf)
/// * [n2235.pdf](http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2007/n2235.pdf)
///
/// @{

/// @def CX_XSTR(A)
/// Macro converts preprocessor macro string to a contant string.
///
/// @note
/// Itself using an macro internale macro called CX_STR(A).
///
#define CX_XSTR(A) CX_STR(A)
#define CX_STR(A) #A

/// @def CX_STRCMP(A,B)
/// constexpr function comparing two strings in compile time constant fashion.
///
/// @param[in] A First string.
/// @param[in] B Second string.
/// @return Integer 0 on matching strings and non zero integer otherwise.
///
constexpr int CX_STRCMP(char const* lhs, char const* rhs)
{
return (('\0' == lhs[0]) && ('\0' == rhs[0])) ? 0
: (lhs[0] != rhs[0]) ? (lhs[0] - rhs[0])
: CX_STRCMP( lhs+1, rhs+1 );
}

/// @}
///
//******************************************************************************
Expand Down Expand Up @@ -357,9 +319,72 @@ constexpr int CX_STRCMP(char const* lhs, char const* rhs)
#ifndef SNGL
#define SNGL float
#endif
/// @}
///
//******************************************************************************
///
/// @name C++11 Constant Expression Support
///
/// The following macros and constant expression (constexpr) functions enable the set up of
/// compile time typed contants.
///
/// References in addtion to the C++11 standard itself:
/// * [sac10-constexpr.pdf](http://www.stroustrup.com/sac10-constexpr.pdf)
/// * [n2235.pdf](http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2007/n2235.pdf)
/// * [github constexpr function library](https://github.com/elbeno/constexpr)
///
/// @attention
/// Do NOT use any CX_ prefixed function at runtime. Use them carefully at
/// compile time. Currently, no error checking is done. Aim is to move burden for
/// complex macro definitions from build systems like autotools, cmake, etc into C++
/// itself where the C++11 standard supports it. If found to fly well in
/// practice, will look to make any constant expression functions more robust.
///
/// @{

/// @def CX_STR(A)
/// Macro used by @ref CX_XSTR(A) in conversion of definition to a string value.
///
/// @def CX_XSTR(A)
/// Macro converts preprocessor macro string to a contant string.
///
/// @note
/// Itself using an macro internal macro called CX_STR(A).
///
#define CX_XSTR(A) CX_STR(A)
#define CX_STR(A) #A

static constexpr int CX_STRCMP(char const* lhs, char const* rhs); // Doxygen requires.
/// constexpr function comparing two strings in compile time constant fashion.
///
/// @param[in] lhs First string.
/// @param[in] rhs Second string.
/// @return Integer 0 on matching strings and non zero integer otherwise.
///
static constexpr int CX_STRCMP(char const* lhs, char const* rhs)
{
return (('\0' == lhs[0]) && ('\0' == rhs[0])) ? 0
: (lhs[0] != rhs[0]) ? (lhs[0] - rhs[0])
: CX_STRCMP( lhs+1, rhs+1 );
}

static constexpr POV_INT64 CX_IPOW(POV_INT64 base, int exp, POV_INT64 result = 1); // Doxygen requires.
/// constexpr function implementing integer pow() in compile time constant fashion.
///
/// @note
/// Optional third argument part of implementation. Call as follows: CX_IPOW(10,3).
///
/// @param[in] base
/// @param[in] exp
/// @param[in] result (do NOT specify on actual first call)
/// @return Integer base raised to exponent.
///
static constexpr POV_INT64 CX_IPOW(POV_INT64 base, int exp, POV_INT64 result) {
return exp < 1 ? result : CX_IPOW(base*base, exp/2, (exp % 2) ? result*base : result);
}

/// @def PRECISE_FLOAT
/// The foating-point data type to use within the polysolve() solver.
/// Optionally more accurate foating-point data type to use within root solver related code.
///
/// Normally PRECISE_FLOAT will be set to @ref DBL. However, 'long double' is a
/// precision guaranteed by the C++ standard to exist and to be at least of
Expand All @@ -370,6 +395,10 @@ constexpr int CX_STRCMP(char const* lhs, char const* rhs)
/// data types. Recent GNU g++ compilers, for example, offer the __float128
/// data type which coorresponds to the IEEE 754-2008 binary128 type.
///
/// The following functions support PRECISE_FLOAT:
/// * polysolve(). Known to users as sturm.
/// * solve_quadratic()
///
/// @note
/// The setting is used to set internal 'constextr int PRECISE_DIG' and
/// 'constexpr DBL PRECISE_EPSILON' typed values.
Expand All @@ -378,14 +407,10 @@ constexpr int CX_STRCMP(char const* lhs, char const* rhs)
/// If @ref PRECISE_FLOAT is set to 'float', 'double' or 'long double', the C++11 standard
/// itself defines via cfloat include the macros: FLT_DIG, DBL_DIG, LDBL_DIG, FLT_EPSILON,
/// DBL_EPSILON, LDBL_EPSILON. These macro values are used to set PRECISE_DIG and
/// PRECISE_EPSILON polysolve() values via the C++ constexpr mechanism.
///
/// @note
/// Users access the polysolve() solver via the 'sturm' keyword though it's always
/// used for polynomials of order 5 or more.
/// PRECISE_EPSILON values via the C++ constexpr mechanism.
///
/// @attention
/// polysolve() math with data types not run in hardware or with awkward bit sizes with
/// Math with data types not run in hardware or with awkward bit sizes with
/// respect to the running computer's data bus size will run substantially slower.
///
#ifndef PRECISE_FLOAT
Expand All @@ -409,6 +434,9 @@ constexpr int CX_STRCMP(char const* lhs, char const* rhs)
///
/// Where @ref PRECISE_FLOAT a C++11 standard floating point type this should be
/// the default sqrt. It would be set to sqrtq if using GNU g++'s quadmath library.
/// If using the inbuilt g++ __float128 capability, set to cast to and from (long double)
/// to avoid the use of quadmath - though this of course effectively limits sqrt()
/// accuracy to 'long double'.
///
#ifndef PRECISE_SQRT
#define PRECISE_SQRT sqrt
Expand All @@ -434,8 +462,17 @@ constexpr int CX_STRCMP(char const* lhs, char const* rhs)
#define PRECISE_EPSLN 1.92592994438724e-34
#endif

/// @def PRECISE_DIG
/// Internal 'constexpr int' value for maximum decimal digits for given @ref PRECISE_FLOAT.
/// @var PRECISE_FLT
/// A 'constexpr int' 0 when @ref PRECISE_FLOAT resolves to float type or !=0 otherwise.
///
/// @var PRECISE_DBL
/// A 'constexpr int' 0 when @ref PRECISE_FLOAT resolves to double type or !=0 otherwise.
///
/// @var PRECISE_LDBL
/// A 'constexpr int' 0 when @ref PRECISE_FLOAT resolves to long double type or !=0 otherwise.
///
/// @var PRECISE_DIG
/// A 'constexpr int' maximum decimal digits given @ref PRECISE_FLOAT.
///
/// Set to C+11 standard value where defined and to @ref PRECISE_DIGITS otherwise.
///
Expand All @@ -450,25 +487,48 @@ constexpr int PRECISE_DIG = PRECISE_FLT ?
: DBL_DIG
: FLT_DIG;

/// @def PRECISE_EPSILON
/// Internal 'constexpr DBL' value*2.0 for minimum epsilon step for given @ref PRECISE_FLOAT.
/// @var PRECISE_EPSILON
/// A 'constexpr DBL' value for minimum epsilon step given @ref PRECISE_FLOAT.
///
/// Set to C+11 standard value where defined and to @ref PRECISE_EPSLN otherwise.
/// Set to C+11 standard value*2.0 where defined and to @ref PRECISE_EPSLN*2.0 otherwise.
///
/// @note
/// Using 2.0 multiplier due maths calculating coefficients for higher order polynomials
/// introducing more than single bit/step error in practice.
///

constexpr DBL PRECISE_EPSILON = PRECISE_FLT ?
PRECISE_DBL ?
PRECISE_LDBL ? PRECISE_EPSLN*2.0
: LDBL_EPSILON*2.0
: DBL_EPSILON*2.0
: FLT_EPSILON*2.0;

/// @def POV_DBL_EPSILON
/// Internal 'constexpr DBL' value for minimum epsilon step for given POV-Ray's @ref DBL.
/// @var POV_DBL_IS_FLT
/// A 'constexpr int' 0 when @ref DBL resolves to float type or !=0 otherwise.
///
/// @var POV_DBL_IS_DBL
/// A 'constexpr int' 0 when @ref DBL resolves to double type or !=0 otherwise.
///
/// @var POV_DBL_IS_LDBL
/// A 'constexpr int' 0 when @ref DBL resolves to long double type or !=0 otherwise.
///
/// @var POV_DBL_DIG
/// A 'constexpr int' value for maximum decimal digits for given @ref DBL.
///
/// Set to C+11 standard value where defined and to @ref PRECISE_DIGITS otherwise.
///
constexpr int POV_DBL_IS_FLT = CX_STRCMP(CX_XSTR(DBL), "float");
constexpr int POV_DBL_IS_DBL = CX_STRCMP(CX_XSTR(DBL), "double");
constexpr int POV_DBL_IS_LDBL = CX_STRCMP(CX_XSTR(DBL), "long double");
constexpr int POV_DBL_DIG = POV_DBL_IS_FLT ?
POV_DBL_IS_DBL ?
POV_DBL_IS_LDBL ? PRECISE_DIGITS
: LDBL_DIG
: DBL_DIG
: FLT_DIG;

/// @var POV_DBL_EPSILON
/// A 'constexpr DBL' value for minimum epsilon step for given POV-Ray's @ref DBL.
///
/// Set to C+11 standard value*2.0 where defined and to @ref PRECISE_EPSLN*2.0 otherwise.
///
Expand All @@ -480,17 +540,25 @@ constexpr DBL PRECISE_EPSILON = PRECISE_FLT ?
/// Using 2.0 multiplier due maths calculating coefficients for higher order polynomials
/// introducing more than single bit/step error in practice.
///
constexpr int POV_DBL_IS_FLT = CX_STRCMP(CX_XSTR(DBL), "float");
constexpr int POV_DBL_IS_DBL = CX_STRCMP(CX_XSTR(DBL), "double");
constexpr int POV_DBL_IS_LDBL = CX_STRCMP(CX_XSTR(DBL), "long double");

constexpr DBL POV_DBL_EPSILON = POV_DBL_IS_FLT ?
POV_DBL_IS_DBL ?
POV_DBL_IS_LDBL ? PRECISE_EPSLN*2.0
: LDBL_EPSILON*2.0
: DBL_EPSILON*2.0
: FLT_EPSILON*2.0;


/// @var MIN_ISECT_DEPTH_RETURNED
/// A 'constexpr DBL' value below which roots from primitive objects are not returned.
///
/// The value will track @ref DBL float type and is very roughly the square root
/// of the determined POV_DBL_EPSILON. The plan is to migrate base shapes to
/// this single value instead of the many different thresholds used today.
/// Aiming for both more accuracy and something which automatically adjust to
/// the DBL type used.
///
constexpr DBL MIN_ISECT_DEPTH_RETURNED = POV_DBL_EPSILON*(DBL)CX_IPOW(10,POV_DBL_DIG/2+1);

/// @}
///
//******************************************************************************
Expand Down
77 changes: 43 additions & 34 deletions source/core/math/polynomialsolver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,57 +50,64 @@ namespace pov
* Local preprocessor defines
******************************************************************************/

/// @def FUDGE_FACTOR2
/// Value defining how close the quartic equation is to being a square
/// of a quadratic in the OLD unused solve_quartic version kept for testing.
/// @var FUDGE_FACTOR2
/// @brief const DBL value defining how close quartic equation is to being a square
/// of a quadratic.
///
/// @note
/// The closer this can get to zero before roots disappear, the less the chance
/// you will get spurious roots.
///
/// @attention
/// Used only in the old unused version of solve_quartic().
// In other words not used.
///
const DBL FUDGE_FACTOR2 = -1.0e-5;

/// @def FUDGE_FACTOR3
/// Similar to @ref FUDGE_FACTOR2 at a later stage of the algebraic solver. Also
/// long used only in the old kept for testing version of solve_quartic().
/// @var FUDGE_FACTOR3
/// @brief const DBL value similar to @ref FUDGE_FACTOR2 at a later stage of the
/// algebraic solver.
///
/// @note
/// @ref FUDGE_FACTOR2 and @ref FUDGE_FACTOR3 have been defined so that quartic
/// equations will properly render on fpu/compiler combinations that only have
/// 64 bit IEEE precision. Do not arbitrarily change any of these values.
///
/// @ref FUDGE_FACTOR2 and @ref FUDGE_FACTOR3 have been defined so that quartic
/// equations will properly render on fpu/compiler combinations that only have
/// 64 bit IEEE precision. Do not arbitrarily change any of these values.
/// If you have a machine with a proper fpu/compiler combo - like a Mac with a
/// 68881, then use the native floating format (96 bits) and tune the values for
/// a little less tolerance: something like: factor2 = -1.0e-7, factor3 =
/// 1.0e-10. Twenty five years later the reality is still double accuracy
/// due use of fastmath (not IEEE compliant) compiling, use of SSE Fused
/// Multiply Add instructions, etc.
///
/// If you have a machine with a proper fpu/compiler combo - like a Mac with a
/// 68881, then use the native floating format (96 bits) and tune the values for
/// a little less tolerance: something like: factor2 = -1.0e-7, factor3 =
/// 1.0e-10. Twenty five years later the reality is still double accuracy
/// due use of fastmath (not IEEE compliant) compiling, use of SSE Fused
/// Multiply Add instructions, etc.
/// @attention
/// Used only in the old unused version of solve_quartic().
// In other words not used.
///
const DBL FUDGE_FACTOR3 = 1.0e-7;

/// @def TWO_M_PI_3
/// Value used in solve_cubic() equal to 2.0 * pi / 3.0.
/// @var TWO_M_PI_3
/// const DBL value used in solve_cubic() equal to 2.0 * pi / 3.0.
///
const DBL TWO_M_PI_3 = 2.0943951023931954923084;

/// @def FOUR_M_PI_3
/// Value used in solve_cubic() equal to 4.0 * pi / 3.0.
/// @var FOUR_M_PI_3
/// const DBL value used in solve_cubic() equal to 4.0 * pi / 3.0.
///
const DBL FOUR_M_PI_3 = 4.1887902047863909846168;

/// @def MAX_ITERATIONS
/// Max number of polysolve sturm chain based bisections.
/// @var MAX_ITERATIONS
/// const int max number of polysolve sturm chain based bisections.
///
/// @note
/// regula_falsa() uses twice this value internally as it can be
/// quite slow to converge in the worst case.
///
const int MAX_ITERATIONS = 65;

/// @def SBISECT_MULT_ROOT_THRESHOLD
/// Value below which multiple roots ignored in sturm chained based bisection
/// and a sigle root at the middle of the current interval is returned.
/// @var SBISECT_MULT_ROOT_THRESHOLD
/// const @ref PRECISE_FLOAT value below which multiple roots ignored in sturm
/// chained based bisection and a single root at the middle of the current
/// interval is returned.
///
/// @note
/// Rays near tangent to surface create extremely close roots and instability
Expand All @@ -110,8 +117,9 @@ const int MAX_ITERATIONS = 65;
///
const PRECISE_FLOAT SBISECT_MULT_ROOT_THRESHOLD = (PRECISE_FLOAT)1e-6;

/// @def REGULA_FALSA_THRESHOLD
/// Threshold below which regula_falsa function is tried when there is a single root.
/// @var REGULA_FALSA_THRESHOLD
/// const @ref PRECISE_FLOAT threshold below which regula_falsa function is tried
/// when there is a single root.
///
/// @note
/// Ray interval max_value - min_value threshold below which regula_falsa
Expand All @@ -127,19 +135,20 @@ const PRECISE_FLOAT SBISECT_MULT_ROOT_THRESHOLD = (PRECISE_FLOAT)1e-6;
///
const PRECISE_FLOAT REGULA_FALSA_THRESHOLD = (PRECISE_FLOAT)1.0;

/// @def RELERROR
/// Smallest relative error along the ray when using the polysolve(), sturm
/// chain bisection / regula-falsi method.
/// @var RELERROR
/// const @ref PRECISE_FLOAT smallest relative error along the ray when using
/// the polysolve(), sturm chain bisection / regula-falsi method.
///
const PRECISE_FLOAT RELERROR = (PRECISE_FLOAT)1.0e-12;

/// @def SMALL_ENOUGH
/// Used to filter determinant value in solve_quadratic() in an unusual way.
/// Used in solve_quartic() to filter values ahead of the trailing quadratics.
/// @var SMALL_ENOUGH
/// const @ref DBL value used to filter determinant value in older
/// solve_quadratic() in an unusual way causing artifacts. Used too in
/// solve_quartic() to filter values ahead of the trailing quadratics.
///
/// @todo
/// Suspect value is larger than it should really be and very likely should
/// not be used at all in solve_quadratic() as it is.
/// not be used at all in solve_quartic() as it is.
///
const DBL SMALL_ENOUGH = 1.0e-10;

Expand Down
Loading

0 comments on commit 761dd0b

Please sign in to comment.