diff --git a/source/base/configbase.h b/source/base/configbase.h index a57d95521..f1ae62b61 100644 --- a/source/base/configbase.h +++ b/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,16 +488,15 @@ 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 @@ -468,8 +504,32 @@ constexpr DBL PRECISE_EPSILON = PRECISE_FLT ? : 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,10 +541,6 @@ 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 @@ -492,6 +548,18 @@ constexpr DBL POV_DBL_EPSILON = POV_DBL_IS_FLT ? : 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); + /// @} /// //****************************************************************************** diff --git a/source/core/math/polynomialsolver.cpp b/source/core/math/polynomialsolver.cpp index 8d30a4728..cd0c9f69d 100644 --- a/source/core/math/polynomialsolver.cpp +++ b/source/core/math/polynomialsolver.cpp @@ -50,47 +50,53 @@ 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 @@ -98,9 +104,10 @@ const DBL FOUR_M_PI_3 = 4.1887902047863909846168; /// 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; diff --git a/source/core/math/polynomialsolver.h b/source/core/math/polynomialsolver.h index 1e299d0c2..3352cd223 100644 --- a/source/core/math/polynomialsolver.h +++ b/source/core/math/polynomialsolver.h @@ -8,7 +8,7 @@ /// @parblock /// /// Persistence of Vision Ray Tracer ('POV-Ray') version 3.8. -/// Copyright 1991-2017 Persistence of Vision Raytracer Pty. Ltd. +/// Copyright 1991-2018 Persistence of Vision Raytracer Pty. Ltd. /// /// POV-Ray is free software: you can redistribute it and/or modify /// it under the terms of the GNU Affero General Public License as @@ -55,6 +55,19 @@ class RenderStatistics; * Global preprocessor defines ******************************************************************************/ +/// @def MAX_ORDER +/// Maximum supported polynomial order. +/// +/// @todo +/// This currently carries a large, fixed per polysolve() call +/// memory allocation on the stack. Size is on the order of +/// (MAX_ORDER+1)*int + PRECISE_FLOAT * (MAX_ORDER+1)^2 +/// which impacts performance and performance stability especially +/// for threads > cores. Allocation based on current equation order +/// would be better. My C++ attempts at this all badly slower. C itself +/// supports a struct "flexible array member" feature, however this +/// not a C++11 standard though some C++ compilers support it. +/// #define MAX_ORDER 35 diff --git a/source/core/shape/blob.cpp b/source/core/shape/blob.cpp index 098cb8801..65facca96 100644 --- a/source/core/shape/blob.cpp +++ b/source/core/shape/blob.cpp @@ -119,34 +119,35 @@ namespace pov * Local preprocessor defines ******************************************************************************/ -/// @def DEPTH_TOLERANCE -/// Minimum returned intersection depth for a valid intersection. +/// @var DEPTH_TOLERANCE +/// const DBL value above which intersections are returned. /// -/// @note -/// Value for @ref DEPTH_TOLERANCE was since v1.0 set to 1.0e-2. Further the -/// value was used for both ling used for the minimum returned and as a mindist -/// in all of the determine_influences()'s intersect_ functions. The -/// latter probably done initially to help the initial solve_quartic() function -/// though using such a large value to determine when to add and remove -/// sub-element influences caused artifacts. In 3.7 when the subsurface (SSLT) -/// feature was added a special conditional existed which set both values to -/// 0.0. Currently a 0.0 value is used for all but the returned intersection -/// where @ref MIN_ISECT_DEPTH is now used. In testing all worked with 0.0, but -/// there is a performance hit to returning all >0.0 intersections. For example, -/// internal refelections at distances which must be ignored in downstream code. +/// Currently set to @ref MIN_ISECT_DEPTH_RETURNED. /// -const DBL DEPTH_TOLERANCE = MIN_ISECT_DEPTH; +/// @ref DEPTH_TOLERANCE since v1.0 was 1.0e-2. Further, the value was used for +/// both minimum returned root depth and as the 'mindist' in all of the +/// determine_influences()'s intersect_<> functions. The latter likely done +/// to help the initial, and long not used, solve_quartic() function. In any +/// case, the 1.0e-2 value causes artifacts. +/// +/// In v3.7 the subsurface (SSLT) feature added a special conditional to set +/// both values to 0.0 which worked without issue. Currently this 0.0 value is +/// used for all but the returned intersection where @ref +/// MIN_ISECT_DEPTH_RETURNED now used. +/// +const DBL DEPTH_TOLERANCE = MIN_ISECT_DEPTH_RETURNED; -/// @def INSIDE_TOLERANCE -/// Density tolerance for inside test. +/// @var INSIDE_TOLERANCE +/// const DBL density tolerance for inside test. +/// +/// Currently set to @ref POV_DBL_EPSILON. /// -/// @note -/// Value for @ref INSIDE_TOLERANCE was since v1.0 set to 1.0e-6 which is a good -/// value if running single floats. With double floats a value of 1e-15 is better. -/// This is an offset to the outside of the 0.0 density surface after the blob -/// threshold is subtracted. Now using C++ standard's DBL_DIG macro to set. +/// @ref INSIDE_TOLERANCE was since v1.0 set to 1.0e-6 which is a good value if +/// running single floats. With double floats much smaller is better as this is +/// an offset to the outside of the 0.0 density surface after the blob threshold +/// is subtracted. /// -const DBL INSIDE_TOLERANCE = ((DBL)1.0/pow((DBL)10.0,(DBL)DBL_DIG)); +const DBL INSIDE_TOLERANCE = POV_DBL_EPSILON; /* Ray enters/exits a component. */ const int ENTERING = 0;