Created new constexpr DBL variable MIN_ISECT_DEPTH_RETURNED.

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.

Author: wfpokorny

source/base/configbase.h

@@ -114,44 +114,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 );
-}
-
 /// @}
 ///
 //******************************************************************************
@@ -358,9 +320,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
@@ -371,6 +396,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.
@@ -379,14 +408,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
@@ -410,6 +435,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
@@ -435,8 +463,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.
 ///
@@ -451,25 +488,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.
 ///
@@ -481,17 +541,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);
+
 /// @}
 ///
 //******************************************************************************

source/core/math/polynomialsolver.cpp

@@ -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
@@ -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
@@ -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;