High precision calculations¶
Yade supports high and arbitrary precision Real
type for performing calculations (see [Kozicki2022] for details). All tests and checks pass but still the current support is in testing phase.
The backend library is boost multiprecision
along with corresponding boost math toolkit.
The supported types are following:
type | bits | decimal places [1] | notes |
---|---|---|---|
float |
32 |
6 |
hardware accelerated (not useful, it is only for testing purposes) |
double |
64 |
15 |
hardware accelerated |
long double |
80 |
18 |
hardware accelerated |
boost float128 |
128 |
33 |
depending on processor type it may be hardware accelerated, wrapped by boost |
boost mpfr |
Nbit |
Nbit/(log(2)/log(10)) |
uses external mpfr library, wrapped by boost |
boost cpp_bin_float |
Nbit |
Nbit/(log(2)/log(10)) |
uses boost only, but is slower |
The last two types are arbitrary precision, and their number of bits Nbit
or decimal places is specified as argument during compilation.
Note
See file Real.hpp for details. All Real
types pass the real type concept test from boost concepts. The support for Eigen and CGAL is done with numerical traits.
[1] | The amount of decimal places in this table is the amount of places which are completely determined by the binary represenation. Few additional decimal digits is necessary to fully reconstruct binary representation. A simple python example to demonstrate this fact: for a in range(16): print(1./pow(2.,a)) , shows that every binary digit produces “extra” …25 at the end of decimal representation, but these decimal digits are not completely determined by the binary representation, because for example …37 is impossible to obtain there. More binary bits are necessary to represent …37 , but the …25 was produced by the last available bit. |
Installation¶
The precompiled Yade Daily packages for Ubuntu 22.04 and Debian Bookworm, Trixie are provided for high precision types long double
, float128
and mpfr150
.
To use high precision on other linux distributions Yade has to be compiled and installed from source code by following the
regular installation instructions. With extra following caveats:
Following packages are required to be installed:
python3-mpmath
libmpfr-dev
libmpfrc++-dev
libmpc-dev
(thempfr
andmpc
related packages are necessary only to useboost::multiprecision::mpfr
type). These packages are already listed in the default requirements.A g++ compiler version 9.2.1 or higher is required. It shall be noted that upgrading only the compiler on an existing linux installation (an older one, in which packages for different versions of gcc were not introduced) is difficult and it is not recommended. A simpler solution is to upgrade entire linux installation.
During cmake invocation specify:
- either number of bits as
REAL_PRECISION_BITS=……
, - or number of requested decimal places as
REAL_DECIMAL_PLACES=……
, but not both - optionally to use MPFR specify
ENABLE_MPFR=ON
(isOFF
by default). - optionally decide about using quadruple, octuple or higher precisions with
-DENABLE_MULTI_REAL_HP=ON
(default). This feature is independent of selecting the precision ofReal
type (in point 1. or 2. above) and works even whenReal
is chosen asdouble
(i.e. no special choice is made: the default settings).
The arbitrary precision (
mpfr
orcpp_bin_float
) types are used only when more than 128 bits or more than 39 decimal places are requested. In such case ifENABLE_MPFR=OFF
then the slowercpp_bin_float
type is used. The difference in decimal places between 39 and 33 stems from the fact that 15 bits are used for exponent. Note: a fast quad-double (debian packagelibqd-dev
) implementation with 62 decimal places is in the works with boost multiprecision team.- either number of bits as
Supported modules¶
During compilation several Yade modules can be enabled or disabled by passing an ENABLE_*
command line argument to cmake.
The following table lists which modules are currently working with high precision (those marked with “maybe” were not tested):
ENABLE_* module name |
HP support | cmake default setting | notes |
---|---|---|---|
ENABLE_GUI |
yes | ON |
native support [2] |
ENABLE_CGAL |
yes | ON |
native support [2] |
ENABLE_VTK |
yes | ON |
supported [3] |
ENABLE_OPENMP |
partial | ON |
partial support [4] |
ENABLE_MPI |
maybe | OFF |
not tested [5] |
ENABLE_GTS |
yes | ON |
supported [6] |
ENABLE_GL2PS |
yes | ON |
supported [6] |
ENABLE_LINSOLV |
no | OFF |
not supported [7] |
ENABLE_PARTIALSAT |
no | OFF |
not supported [7] |
ENABLE_PFVFLOW |
no | OFF |
not supported [7] |
ENABLE_TWOPHASEFLOW |
no | OFF |
not supported [7] |
ENABLE_THERMAL |
no | OFF |
not supported [7] |
ENABLE_LBMFLOW |
yes | ON |
supported [6] |
ENABLE_SPH |
maybe | OFF |
not tested [8] |
ENABLE_LIQMIGRATION |
maybe | OFF |
not tested [8] |
ENABLE_MASK_ARBITRARY |
maybe | OFF |
not tested [8] |
ENABLE_PROFILING |
maybe | OFF |
not tested [8] |
ENABLE_POTENTIAL_BLOCKS |
no | OFF |
not supported [9] |
ENABLE_POTENTIAL_PARTICLES |
yes | ON |
supported [10] |
ENABLE_DEFORM |
maybe | OFF |
not tested [8] |
ENABLE_OAR |
maybe | OFF |
not tested [8] |
ENABLE_FEMLIKE |
yes | ON |
supported [6] |
ENABLE_ASAN |
yes | OFF |
supported [6] |
ENABLE_MPFR |
yes | OFF |
native support [2] |
ENABLE_LS_DEM |
no | ON |
not supported [11] |
The unsupported modules are automatically disabled during a high precision cmake
stage.
Footnotes
[2] | (1, 2, 3) This feature is supported natively, which means that specific numerical traits were written for Eigen and for CGAL, as well as GUI and python support was added. |
[3] | VTK is supported via the compatibility layer which converts all numbers down to double type. See below. |
[4] | The OpenMPArrayAccumulator is experimentally supported for long double and float128 . For types mpfr and cpp_bin_float the single-threaded version of accumulator is used. File lib/base/openmp-accu.hpp needs further testing. If in doubt, compile yade with ENABLE_OPENMP=OFF . In all other places OpenMP multithreading should work correctly. |
[5] | MPI support has not been tested and sending data over network hasn’t been tested yet. |
[6] | (1, 2, 3, 4, 5) The module was tested, the yade --test and yade --check pass, as well as most of examples are working. But it hasn’t been tested extensively for all possible use cases. |
[7] | (1, 2, 3, 4, 5) Not supported, the code uses external cholmod library which supports only double type. To make it work a native Eigen solver for linear equations should be used. |
[8] | (1, 2, 3, 4, 5, 6) This feature is OFF by default, the support of this feature has not been tested. |
[9] | Potential blocks use external library coinor for linear programming, this library uses double type only. To make it work a linear programming routine has to be implemented using Eigen or coinor library should start using C++ templates or a converter/wrapper similar to LAPACK library should be used. |
[10] | The module is enabled by default, the yade --test and yade --check pass, as well as most of examples are working. However the calculations are performed at lower double precision. A wrapper/converter layer for LAPACK library has been implemented. To make it work with full precision these routines should be reimplemented using Eigen. |
[11] | Possible future enchancement. See comments there . |
Double, quadruple, octuple and higher precisions¶
Sometimes a critical section of the calculations in C++ would work better if it was performed in the higher precision to guarantee that it will produce the correct result in the default precision. A simple example is solving a system of linear equations (basically inverting a matrix) where some coefficients are very close to zero. Another example of alleviating such problem is the Kahan summation algorithm.
If requirements are satisfied, Yade supports higher precision multipliers in such a way that RealHP<1>
is the Real
type described above, and every higher number is a multiplier of the Real
precision. RealHP<2>
is double precision of RealHP<1>
, RealHP<4>
is quadruple precision and so on. The general formula for amount of decimal places is implemented in RealHP.hpp file and the number of decimal places used is simply a multiple N of decimal places in Real
precision, it is used when native types are not available. The family of available native precision types is listed in the RealHPLadder type list.
All types listed in MathEigenTypes.hpp follow the same naming pattern: Vector3rHP<1>
is the regular Vector3r
and Vector3rHP<N>
for any supported N uses the precision multiplier N. One could then use an Eigen algorithm for solving a system of linear equations with a higher N using MatrixXrHP<N>
to obtain the result with higher precision. Then continuing calculations in default Real
precision, after the critical section is done. The same naming convention is used for CGAL types, e.g. CGAL_AABB_treeHP<N>
which are declared in file AliasCGAL.hpp.
Before we fully move to C++20 standard, one small restriction is in place: the precision multipliers actually supported are determined by these two defines in the RealHPConfig.hpp file:
#define YADE_EIGENCGAL_HP (1)(2)(3)(4)(8)(10)(20)
- the multipliers listed here will work in C++ forRealHP<N>
in CGAL and Eigen. They are cheap in compilation time, but have to be listed here nonetheless. After we move code to C++20 this define will be removed and all multipliers will be supported via single template constraint. This inconvenience arises from the fact that both CGAL and Eigen libraries offer template specializations only for a specific type, not a generalized family of types. Thus this define is used to declare the required template specializations.
Hint
The highest precision available by default N= (20)
corresponds to 300 decimal places when compiling Yade with the default settings, without changing REAL_DECIMAL_PLACES=……
cmake compilation option.
#define YADE_MINIEIGEN_HP (1)(2)
- the precision multipliers listed here are exported to python, they are expensive: each one makes compilation longer by 1 minute. Adding more can be useful only for debugging purposes. The doubleRealHP<2>
type is by default listed here to allow exploring the higher precision types from python. Also please note thatmpmath
supports only one precision at a time. Having differentmpmath
variables with different precision is poorly supported, albeitmpmath
authors promise to improve that in the future. Fortunately this is not a big problem for Yade users because the general goal here is to allow more precise calculations in the critical sections of C++ code, not in python. This problem is partially mitigated by changing mpmath precision each time when aC++
↔python
conversion occurs. So one should keep in mind that the variablempmath.mp.dps
always reflects the precision used by latest conversion performed, even if that conversion took place in GUI (not in the running script). Existingmpmath
variables are not truncated to lower precision, their extra digits are simply ignored untilmpmath.mp.dps
is increased again, however the truncation might occur during assignment.
On some occasions it is useful to have an intuitive up-conversion between C++ types of different precisions, say for example to add RealHP<1>
to RealHP<2>
type. The file UpconversionOfBasicOperatorsHP.hpp serves this purpose. This header is not included by default, because more often than not, adding such two different types will be a mistake (efficiency–wise) and compiler will catch them and complain. After including this header this operation will become possible and the resultant type of such operation will be always the higher precision of the two types used. This file should be included only in .cpp
files. If it was included in any .hpp
file then it could pose problems with C++ type safety and will have unexpected consequences. An example usage of this header is in the following test routine.
Warning
Trying to use N unregistered in YADE_MINIEIGEN_HP
for a Vector3rHP<N>
type inside the YADE_CLASS_BASE_DOC_ATTRS_*
macro to export it to python will not work. Only these N listed in YADE_MINIEIGEN_HP
will work. However it is safe (and intended) to use these from YADE_EIGENCGAL_HP
in the C++ calculations in critical sections of code, without exporting them to python.
Compatibility¶
Python¶
To declare python variables with Real
and RealHP<N>
precision use functions math.Real(…), math.Real1(…), math.Real2(…). Supported are precisions listed in YADE_MINIEIGEN_HP
, but please note the mpmath-conversion-restrictions.
Python has native support for high precision types using mpmath
package. Old Yade scripts that use supported modules can be immediately converted to high precision by switching to yade.minieigenHP
. In order to do so, the following line:
from minieigen import *
has to be replaced with:
from yade.minieigenHP import *
Respectively import minieigen
has to be replaced with import yade.minieigenHP as minieigen
, the old name as minieigen
being used only for the sake of backward compatibility. Then high precision (binary compatible) version of minieigen is used when non double
type is used as Real
.
The RealHP<N>
higher precision vectors and matrices can be accessed in python by using the .HPn
module scope. For example:
import yade.minieigenHP as mne
mne.HP2.Vector3(1,2,3) # produces Vector3 using RealHP<2> precision
mne.Vector3(1,2,3) # without using HPn module scope it defaults to RealHP<1>
The respective math functions such as:
import yade.math as mth
mth.HP2.sqrt(2) # produces square root of 2 using RealHP<2> precision
mth.sqrt(2) # without using HPn module scope it defaults to RealHP<1>
are supported as well and work by using the respective C++ function calls, which is usually faster than the mpmath
functions.
Warning
There may be still some parts of python code that were not migrated to high precision and may not work well with mpmath
module. See debugging section for details.
C++¶
Before introducing high precision it was assumed that Real
is actually a POD double
type. It was possible to use memset(…)
, memcpy(…)
and similar functions on double
. This was not a good approach and even some compiler #pragma
commands were used to silence the compilation warnings. To make Real
work with other types, this assumption had to be removed. A single memcpy(…)
still remains in file openmp-accu.hpp and will have to be removed. In future development such raw memory access functions are to be avoided.
All remaining double
were replaced with Real
and any attempts to use double
type in the code will fail in the gitlab-CI pipeline.
Mathematical functions of all high precision types are wrapped using file MathFunctions.hpp, these are the inline redirections to respective functions of the type that Yade is currently being compiled with. The code will not pass the pipeline checks if std::
is used. All functions that take Real
argument should now call these functions in yade::math::
namespace. Functions which take only Real
arguments may omit math::
specifier and use ADL instead. Examples:
- Call to
std::min(a,b)
is replaced withmath::min(a,b)
, becausea
orb
may beint
(nonReal
) thereforemath::
is necessary. - Call to
std::sqrt(a)
can be replaced with eithersqrt(a)
ormath::sqrt(a)
thanks to ADL, becausea
is alwaysReal
.
If a new mathematical function is needed it has to be added in the following places:
- lib/high-precision/MathFunctions.hpp or lib/high-precision/MathComplexFunctions.hpp or lib/high-precision/MathSpecialFunctions.hpp, depending on function type.
- py/high-precision/_math.cpp, see math module for details.
- py/tests/testMath.py
- py/tests/testMathHelper.py
The tests for a new function are to be added in py/tests/testMath.py in one of these functions: oneArgMathCheck(…):
, twoArgMathCheck(…):
, threeArgMathCheck(…):
. A table of approximate expected error tolerances in self.defaultTolerances
is to be supplemented as well. To determine tolerances with better confidence it is recommended to temporarily increase number of tests in the test loop. To determine tolerances for currently implemented functions a range(1000000)
in the loop was used.
Note
When passing arguments in C++
in function calls it is preferred to use const Real&
rather than to make a copy of the argument as Real
. The reason is following: in non high-precision
regular case both the double
type and the reference have 8 bytes. However float128
is 16 bytes large, while its reference is still only 8 bytes.
So for regular precision, there is no difference. For all higher precision types it is beneficial to use const Real&
as the function argument. Also for const Vector3r&
arguments
the speed gain is larger, even without high precision.
Using higher precisions in C++¶
As mentioned above RealHP<1>
is the Real
type and every higher number is a multiplier of the Real
precision. RealHP<2>
is twice the precision of RealHP<1>
, RealHP<4>
is quadruple precision and so on. In C++ you have access to these higher precision typedefs at all time, so it is possible to write some critical part of an algorithm in higher precision by declaring the respective variables to be of type RealHP<2>
or RealHP<4>
or higher.
String conversions¶
On the python
side it is recommended to use math.Real(…) math.Real1(…), or math.toHP1(…) to declare python
variables and math.radiansHP1(…) to convert angles to radians using full Pi precision.
On the C++
side it is recommended to use yade::math::toString(…) and yade::math::fromStringReal(…) conversion functions instead of boost::lexical_cast<std::string>(…)
. The toString
and its high precision version toStringHP functions (in file RealIO.hpp) guarantee full precision during conversion. It is important to note that std::to_string
does not guarantee this and boost::lexical_cast
does not guarantee this either.
For higher precision types it is possible to control in runtime the precision of C++
↔ python
during the RealHP<N>
string conversion by changing the math.RealHPConfig.extraStringDigits10 static parameter. Each decimal digit needs \(\log_{10}(2)\approx3.3219\) bits. The std::numeric_limits<Real>::digits10
provides information about how many decimal digits are completely determined by binary representation, meaning that these digits are absolutely correct. However to convert back to binary more decimal digits are necessary because \(\log_{2}(10)\approx0.3010299\) decimal digits are used by each bit, and the last digit from std::numeric_limits<Real>::digits10
is not sufficient. In general 3 or more in extraStringDigits10 is enough to have an always working number round tripping. However if one wants to only extract results from python, without feeding them back in to continue calculations then a smaller value of extraStringDigits10 is recommended, like 0 or 1, to avoid a fake sense of having more precision, when it’s not there: these extra decimal digits are not correct in decimal sense. They are only there to have working number round tripping. See also a short discussion about this with boost developers. Also see file RealHPConfig.cpp for more details.
Note
The parameter extraStringDigits10
does not affect double
conversions, because boost::python
uses an internal converter for this particular type. It might be changed in the future if the need arises. E.g. using a class similar to ThinRealWrapper.
It is important to note that creating higher types such as RealHP<2>
from string representation of RealHP<1>
is ambiguous. Consider following example:
import yade.math as mth
mth.HP1.getDecomposedReal(1.23)['bits']
Out[2]: '10011101011100001010001111010111000010100011110101110'
mth.HP2.getDecomposedReal('1.23')['bits'] # passing the same arg in decimal format to HP2 produces nonzero bits after the first 53 bits of HP1
Out[3]: '10011101011100001010001111010111000010100011110101110000101000111101011100001010001111010111000010100011110101110'
mth.HP2.getDecomposedReal(mth.HP1.toHP2(1.23))['bits'] # it is possible to use yade.math.HPn.toHPm(…) conversion, which preserves binary representation
Out[4]: '10011101011100001010001111010111000010100011110101110000000000000000000000000000000000000000000000000000000000000'
Which of these two RealHP<2>
binary representations is more desirable depends on what is needed:
- The best binary approximation of a
1.23
decimal. - Reproducing the 53 binary bits of that number into a higher precision to continue the calculations on the same number which was previously in lower precision.
To achieve 1. simply pass the argument '1.23'
as string. To achieve 2. use math.HPn.toHPm(…) or math.Realn(…) conversion, which maintains binary fidelity using a single static_cast<RealHP<m>>(…). Similar problem is discussed in mpmath and boost documentation.
The difference between toHPn and Realn is following: the functions HPn.toHPm
create a \(m\times n\) matrix converting from RealHP<n>
to RealHP<m>
. When \(n<m\) then extra bits are set to zero (case 2 above, depending on what is required one might say that “precision loss occurs”). The functions math.Real(…), math.Real1(…), math.Real2(…) are aliases to the diagonal of this matrix (case 1 above, depending on what is required one might say that “no conversion loss occurs” when using them).
Hint
All RealHP<N>
function arguments that are of type higher than double
can also accept decimal strings. This allows to preserve precision above python default floating point precision.
Warning
On the contrary all the function arguments that are of type double
can not accept decimal strings. To mitigate that one can use toHPn(…)
converters with string arguments.
Hint
To make debugging of this problem easier the function math.toHP1(…) will raise RuntimeError if the argument is a python float (not a decimal string).
Warning
I cannot stress this problem enough, please try running yade --check
(or yade ./checkGravityRungeKuttaCashKarp54.py
) in precision different than double
after changing this line into g = -9.81
. In this (particular and simple) case the getCurrentPos()
function fails on the python side because low-precision g
is multiplied by high-precision t
.
Complex types¶
Complex numbers are supported as well. All standard C++
functions are available in lib/high-precision/MathComplexFunctions.hpp and also are exported to python
in py/high-precision/_math.cpp. There is a cmake compilation option ENABLE_COMPLEX_MP
which enables
using better complex
types from boost::multiprecision
library for representing ComplexHP<N>
family of types: complex128
, mpc_complex
, cpp_complex
and complex_adaptor
.
It is ON by default whenever possible: for boost version >= 1.71. For older boost the ComplexHP<N>
types are represented by std::complex<RealHP<N>>
instead, which has larger
numerical errors in some mathematical functions.
When using the ENABLE_COMPLEX_MP=ON
(default) the previously mentioned lib/high-precision/UpconversionOfBasicOperatorsHP.hpp is not functional for complex types,
it is a reported problem with the boost library.
When using MPFR type, the libmpc-dev
package has to be installed (mentioned above).
Eigen and CGAL¶
Eigen and CGAL libraries have native high precision support.
- All declarations required by Eigen are provided in files EigenNumTraits.hpp and MathEigenTypes.hpp
- All declarations required by CGAL are provided in files CgalNumTraits.hpp and AliasCGAL.hpp
VTK¶
Since VTK is only used to record results for later viewing in other software, such as paraview, the recording of all decimal places does not seem to be necessary (for now).
Hence all recording commands in C++
convert Real
type down to double
using static_cast<double>
command. This has been implemented via classes vtkPointsReal
, vtkTransformReal
and vtkDoubleArrayFromReal
in file VTKCompatibility.hpp. Maybe VTK in the future will support non double
types. If that will be needed, the interface can be updated there.
LAPACK¶
Lapack is an external library which only supports double
type. Since it is not templatized it is not possible to use it with Real
type. Current solution is to down-convert arguments to double
upon calling linear equation solver (and other functions), then convert them back to Real
. This temporary solution omits all benefits of high precision, so in the future Lapack is to be replaced with Eigen or other templatized libraries which support arbitrary floating point types.
Debugging¶
High precision is still in the experimental stages of implementation. Some errors may occur during use. Not all of these errors are caught by the checks and tests. Following examples may be instructive:
- Trying to use const references to Vector3r members - a type of problem with results in a segmentation fault during runtime.
- A part of python code does not cooperate with mpmath - the checks and tests do not cover all lines of the python code (yet), so more errors like this one are expected. The solution is to put the non compliant python functions into py/high-precision/math.py. Then replace original calls to this function with function in
yade.math
, e.g.numpy.linspace(…)
is replaced withyade.math.linspace(…)
.
The most flexibility in debugging is with the long double
type, because special files ThinRealWrapper.hpp, ThinComplexWrapper.hpp were written for that. They are implemented with boost::operators, using partially ordered field. Note that they do not provide operator++.
A couple of #defines
were introduced in these two files to help debugging more difficult problems:
YADE_IGNORE_IEEE_INFINITY_NAN
- it can be used to detect all occurrences whenNaN
orInf
are used. Also it is recommended to use this define when compiling Yade with-Ofast
flag, without-fno-associative-math -fno-finite-math-only -fsigned-zeros
YADE_WRAPPER_THROW_ON_NAN_INF_REAL
,YADE_WRAPPER_THROW_ON_NAN_INF_COMPLEX
- can be useful for debugging when calculations go all wrong for unknown reason.
Also refer to address sanitizer section, as it is most useful for debugging in many cases.
Hint
If crash is inside a macro, for example YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY
, it is useful to know where inside this macro the problem happens. For this purpose it is possible to use g++
preprocessor to remove the macro and then compile the postprocessed code without the macro. Invoke the preprocessor with some variation of this command:
g++ -E -P core/Body.hpp -I ./ -I /usr/include/eigen3 -I /usr/include/python3.7m > /tmp/Body.hpp
Maybe use clang-format so that this file is more readable:
./scripts/clang-formatter.sh /tmp/Body.hpp
Be careful because such files tend to be large and clang-format is slow. So sometimes it is more useful to only use the last part of the file, where the macro was postprocessed. Then replace the macro in the original file in question, and then continue debugging. But this time it will be revealed where inside a macro the problem occurs.
Note
When asking questions about High Precision it is recommended to start the question title with [RealHP]
.