Implementation of high-precision computation capabilities into the open-source dynamic simulation framework YADE

This paper deals with the implementation of arbitrary precision calculations into the open-source discrete element framework YADE published under the GPL-2+ free software license. This new capability paves the way for the simulation framework to be used in many new fields such as quantum mechanics. The implementation details and associated gains in the accuracy of the results are discussed. Besides the"standard"double (64 bits) type, support for the following high-precision types is added: long double (80 bits), float128 (128 bits), mpfr_float_backend (arbitrary precision) and cpp_bin_float (arbitrary precision). Benchmarks are performed to quantify the additional computational cost involved with the new supported precisions. Finally, a simple calculation of a chaotic triple pendulum is performed to demonstrate the new capabilities and the effect of different precisions on the simulation result.


Introduction
The advent of a new era of scientific computing has been predicted in the literature, one in which the numerical precision required for computation is as important as the algorithms and data structures in the program [6,50,51,64]. An appealing example of a simple computation gone wrong was presented in the talk "Are we just getting wrong answers faster?" of Stadtherr in 1998 [85]. An exhaustive list of such computations along with a very detailed analysis can be found in [51].
Many examples exist where low-precision calculations resulted in disasters. The military identified an accumulated error in multiplication by a constant factor of 0.1, which has no exact binary representation, as the cause for a Patriot missile failure on 25 February 1991, which resulted in several fatalities [12]. If more bits were used to represent a number, the explosion of an Ariane 5 rocket launched by the European Space Agency on 4 June 1996 could have been prevented [58,56,49,11] as it was a result of an inappropriate conversion from a 64 bit floating point number into a 16 bit signed integer. Indeed, the 64 bit floating point number was too big to be represented as a 16 bit signed integer. On 14 May 1992 the rendezvous between the shuttle Endeavour and the Intelsat 603 spacecraft nearly failed. The problem was traced back to a mismatch in precision [74,43]. More catastrophic failures related to the lack of precision are discussed in [74,43]. In 2012 it was predicted that most future technical computing will be performed by people with only basic training in numerical analysis or none at all [59,36,6,51]. High-precision computation is an

cmake flag description
High-precision support in present YADE version 37 .

ENABLE LBMFLOW
Fluid-solid interaction in granular media with coupled Lattice Boltzmann/Discrete Element Method [60]. ENABLE POTENTIAL PARTICLES Arbitrarily shaped convex particle described as a 2nd degree polynomial potential function [14].
Selected YADE features with high-precision support.

ENABLE ASAN
AddressSanitizer allows detection of memory errors, memory leaks, heap corruption errors and out-of-bounds accesses [83].

ENABLE OPENMP
OpenMP threads parallelization, full support for double, long double, float128 types 8 .
Modules under development for high-precision support.

ENABLE MPI
MPI environment for massively parallel computation [98].

ENABLE NRQM
Quantum dynamics simulations of diatomic molecules including photoinduced transitions between the coupled states [46,45].
using Real was to allow replacing its definition with other possible precisions 3 . Hence, the same strategy was followed for other types used in the calculations, such as vectors and matrices. Per definition the last letter in the type name indicates its underlying type, e.g. 'Vector3r v;' is a 3D vector v ∈ Q 3 ∈ R 3 , and Vector2i is a 2D vector of integers (where Q is a subset of rational numbers Q, which are representable by the currently used precision: Q ∈ Q ∈ R; the name Real is used instead of Rational or FloatingPoint for the sake of brevity).
In the presented work, the goal to use high precision is achieved by using the C++ operator overloading functionality and the boost::multiprecision library. A simplified dependency diagram of YADE is shown in Fig. 1. The layered structure of YADE remains nearly the same as in the original paper by Kozicki and Donzé [53]. It is built on top of several well established libraries (marked with orange in Fig. 1) as discussed in Section 2.3. Some changes were necessary in the structure of the framework (marked with green in Fig. 1) as highlighted in Section 2.5. The top row in Fig. 1 indicates selected YADE modules with respective citations listed in Tab. 1. It should be noted that YADE relies on many external libraries to expand its functionality which can result in a demanding server setup.
The Boost library [24] provides convenient wrappers for other high-precision types with the perspective of adding more such types in the future 4  The process of adding high-precision support to YADE was divided into several stages which are described in the subsections below 6 .

Preparations
To fully take advantage of the C++ Argument Dependent Lookup (ADL), the entire YADE codebase was moved into namespace yade 7 , thus using the C++ standard capabilities to modularize the namespaces for each software package. Similarly, the libraries used by YADE such as Boost [24], CGAL [93] and 6 also see the consolidated merge request: !383 7 see: !284 4 EIGEN [37] reside in their respective boost, CGAL and Eigen namespaces. After this change, all potential naming conflicts between math functions or types in YADE and these libraries were eliminated.
Before introducing high precision into YADE it was assumed that Real is actually a Plain Old Data (POD) double type. It was possible to use the old C-style memset, memcpy, memcmp and memmove functions which used raw-memory access. However, by doing so the modern C++ structure used by other highprecision types was completely ignored. For example, the MPFR type may reserve memory and inside its structure store a pointer to it. Trying to set its value to zero by invoking memset (which sets that pointer to nullptr) leads to a memory leak and a subsequent program failure. In order to make Real work with other types, this assumption had to be removed. Hence, memset calls were replaced with std::fill calls, which when invoked with a POD type reduce to a (possibly faster) version of memset optimized for a particular type in terms of chunk size used for writing to the memory. In addition, C++ template specialization mechanisms allow for invoking with a non-POD type which then utilizes the functionality provided by this specific type, such as calling the specific constructors. All places in the code which used these four rawmemory access functions were improved to work with the non-POD Real type 8 . For similar reasons one should not rely on storing an address of the n th component of a Vector3r or Vector2r 9 .
Next, all remaining occurrences of double were replaced with Real 10 and the high-precision compilation and testing was added to the gitlab Continuous Integration (CI) testing pipeline, which guarantees that any future attempts to use double type in the code will fail before merging such changes into the main branch. Next the Real type was moved from global namespace into yade namespace 11 to eliminate any potential problems with namespace pollution 12 .

Library compatibility
In order to be able to properly interface YADE with all other libraries it was important to make sure that mathematical functions (see Tab. 4) are called for the appropriate type. For example, the EIGEN library would have to call the high-precision sqrt function when invoking a normalize function on a Vector3r in order to properly calculate vector length. Several steps were necessary to achieve this. First, an inline redirection 13 to these functions was implemented in namespace yade::math in the file MathFunctions.hpp. Next, all invocations in YADE to math functions in the std namespace were replaced with calls to these functions in the yade::math namespace 14 . Functions which take only Real arguments may omit math namespace specifier and use ADL instead. Also some fixes were done in EIGEN and CGAL 15 , although they did not affect YADE directly since it was possible to workaround them.
The C++ type traits is a template metaprogramming technique which allows one to customize program behavior (also called polymorphism) depending on the type used [2,67,89,96]. This decision is done by the compiler (conditional compilation) due to inspecting the types in the compilation stage (this is called static polymorphism). Advanced C++ libraries provide hooks (numerical traits) to allow library users to inform the library about the used precision type. The numerical traits were implemented in YADE for the libraries EIGEN and CGAL 16 as these were the only libraries supporting such a solution at the time of 8 see: !381, with one exception which is yet to be evaluated and hence mpfr and cpp bin float types work with OpenMP, but do not take full advantage of CPU cache size in class OpenMPArrayAccumulator 9 see: !406 10 these changes were divided into several smaller merge requests: !326, !376, !394; there were also a couple of changes such as MatrixXd → MatrixXr and Vector3d → Vector3r. 11 see: !364 12 usually such errors manifest themselves as very unrelated problems, which are notoriously difficult to debug, e.g. due to the fact that an incorrect type (with the same name) is used; see: #57 and https://bugs.launchpad.net/yade/+bug/528509. 13 The recommended practice in such cases is to use the Argument Dependent Lookup (ADL) which lets the compiler pick the best match from all the available candidate functions [2,67,89,96]. No ambiguity is possible, because such situations would always result in a compiler error. This was done by employing the C++ directives using std::function ; and using boost::multiprecision::function ; for the respective function and then calling the function unqualified (without namespace qualifier) in the MathFunctions.hpp file. 14  The LAPACK compatibility layer was provided as well, this time highlighting the problems of interfacing with languages which do not support static polymorphism. The routines in LAPACK are written in a mix of Fortran and C, and have no capability to use high-precision numerical traits like EIGEN and CGAL. The only way to do this (apart from switching to another library) was to down-convert the arguments to double upon calling LAPACK routines (e.g. a routine to solve a linear system) then up-converting the results to Real. This was the first step to phase out YADE's dependency on LAPACK. With this approach the legacy code works even when high precision is enabled although the obtained results are low-precision 19 . Additionally, this allows one to test the new high-precision code against the low-precision version when replacing these function calls with appropriate function calls from another library such as EIGEN in the future. Fortunately, only two YADE modules depend on LAPACK: potential particles and the flow engine [21]. The latter also depends on CHOLMOD, which also supports the double type only, hence it is not shown in Tab. 1. Nevertheless, a similar solution as currently implemented for LAPACK can be used in the future to remove the current dependency on CHOLMOD.

Double, quadruple and higher precisions
Sometimes a critical section of the computations in C++ would work better if performed in a higher precision 1 . This would also guarantee that the overall results in the default precision are correct. The RealHP<N> types serve this purpose. In analogy to float and double types used on older systems, the types RealHP<2>, RealHP<4> and RealHP<N> correspond to double, quadruple and higher multipliers of the Real precision selected during compilation, e.g. with REAL DECIMAL PLACES 5 , respectively. A simple example where this can be useful is solving a system of linear equations where some coefficients are almost zero. The old rule of thumb to ,,perform all computation in arithmetic with somewhat more than twice as many significant digits as are deemed significant in the data and are desired in the final results" works well in many cases [50]. Nevertheless, maintaining a high quality scientific software package without being able to use, when necessary, arithmetic precision twice as wide can badly inflate costs of development and maintenance [50]. On the one hand, there might be additional costs for the theoretical formulation of such tricky single-precision problems. On the other hand, the cost of extra demand for processor cycles and memory when using RealHP<N> types is picayune when compared with the cost of a numerically adept mathematician's time [51]. Hence, the new RealHP<N> makes high and multiple-precision simulations more accessible to the researcher community.
The support for higher precision multipliers was added in YADE 20 in such a way that RealHP<1> is the Real type from Tab. 2 and every higher number N is a multiplier of the Real precision. All other types follow the same naming pattern: Vector3rHP<1> is the regular Vector3r and Vector3rHP<N> uses the precision multiplier N. A similar concept is used for CGAL types (e.g. CGALtriangleHP<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, after the critical code section, one could potentially continue the calculations in the default Real precision. On the Python side the mathematical functions for the higher precision types are accessible via yade.math.HP2.*. By default only the RealHP<2> is exported 17 see: !412 and file OpenGLWrapper.hpp. If the need for drawing on screen with precision higher than double arises (e.g. at high zoom levels) it will be rectified in the future.
18 see: !400 and VTKCompatibility.hpp. If the VTK display software will start supporting high precision, this solution can be readily improved.
19 see: !379 and LapackCompatibility.cpp. 20 see: !496 6 to Python. One can export to Python all the higher types for debugging purposes by adjusting #define YADE MINIEIGEN HP in the file RealHPConfig.hpp. 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>. The file UpconversionOfBasicOperatorsHP.hpp serves this purpose. After including this header, operations using two different precision types are possible and the resultant type of such operation will always be the higher precision of the two types. This header should be used with caution (and only in .cpp files) in order to still be able to take advantage of the C++ static type checking mechanisms. As mentioned in the introduction, this type checking whether a number is being converted to a fewer digits representation can prevent mistakes such as the explosion of the rocket Ariane 5 [58,56,49,11].

Backward compatibility with older YADE scripts
In the present work, preserving the backward compatibility with existing older YADE Python scripts was of prime importance. To obtain this the MiniEigen Python library had to be incorporated into YADE's codebase. The reason for this was the following: python3-minieigen was a binary package, precompiled using double. Thus any attempt of importing MiniEigen into a YADE Python script (i.e. using from minieigen import *) when YADE was using a non-double type resulted in failure. This, combined with the new capability in YADE to use any of the current and future supported types (see Tab. 2) would place a requirement on python3-minieigen that it either becomes a header-only library or is precompiled with all possible high-precision types. It was concluded that integrating its source directly into YADE is the most reasonable solution. Hence, old YADE scripts that use supported modules 37 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 for the sake of backward compatibility with the rest of the script.
Python has native support 21 for high-precision types using the mpmath Python package. However, it shall be noted that although the coverage of YADE's basic testing and checking (i.e. yade --test and yade --check) is fairly large, there may still be some parts of Python code that were not yet migrated to high precision and may not work well with the mpmath module. If such problems occur in the future, the solution is to put the non compliant Python function into the py/high-precision/math.py file 22 .
A typical way of ensuring correct treatment of Real in Python scripts is to initialize Python variables using yade.math.Real(arg). If the initial argument is not an integer and not an mpmath type then it has to be passed as a string (e.g. yade.math.Real('9.81')) to prevent Python from converting it to double. Without this special initialization step a mistake can appear in the Python script where the default Python floating-point type double is for example multiplied or added to the Real type resulting in a loss of precision 23 .

Testing
It should be noted that it is near to impossible to be absolutely certain about the lack of an error in a code [51]. Therefore, to briefly test the implementation of all mathematical functions available in C++ in all precisions, the following test was implemented in RealHPDiagnostics.cpp and run using the Python script testMath.py. Each available function was evaluated 2 × 10 7 times with on average evenly spaced  pseudo-random argument in the range (−100, 100). These 2 × 10 7 evaluations were divided into two sets. The first 10 7 evaluations were performed with uniformly distributed pseudo-random numbers in the range (−100, 100). The second 10 7 evaluations were done by randomly displacing a set of 10 7 equidistant points in the range (−100, 100), where each point was randomly shifted by less than ±0.5 of the distance between the points. This random shift was to lower the chances of duplicate calculations on the same argument after adjusting to the function domain. The 2 × 10 7 arguments were subsequently modified to match the domain argument range of each function using simple operations, such as abs(•) or fmod(abs(•),2)-1. The obtained result for each evaluation was then compared against its respective RealHP<4> type with four times higher precision. Care was taken to exactly use the same argument for a higher precision function call. The arguments were randomized and adjusted to the function domain range in the lower precision RealHP<1>, then the argument was converted to the higher precision by using static cast, thereby ensuring that all the extra bits in the higher precision are set to zero.
The difference expressed in terms of Units in the Last Place (ULP) [51] was calculated. The obtained errors are listed in Tab. 4. During the tests a bug in the implementation of the tgamma function for boost::multiprecision::float128 was discovered but it was immediately fixed by the Boost developers 24 . Some other bug reports 25 instigated a discussion about possible ways to fix the few problems found with the cpp bin float type which can be seen in the last column of Tab. 4. A smaller version of this test, with only 2 × 10 4 pseudo-random evaluations 26 was then added to the standard yade --test invocation.
Finally, an AddressSanitizer [83,10] was employed to additionally check the correctness of the implementation in the code and to quickly locate memory access bugs. Several critical errors were fixed due to the reports of this sanity checker. This tool is now integrated into the Continuous Integration (CI) pipeline for the whole YADE project to prevent introduction of such errors in the future [10] (make asan HP job in the GitLab CI pipeline).

Benchmark
A benchmark 27 yade --stdperformance -j16 (16 OpenMP threads [23,71]) on a PC with two Intel E5-2687W v2 @ 3.40GHz processors (each of the two having 8 cores resulting in a total of 16 cores or 32 threads if hyperthreading is enabled) was performed to assess performance of higher precision types. The benchmark consists of a simple gravity deposition of spherical particles into a box, a typical simulation performed in YADE. A spherical packing with 10,000 spheres is released under gravity within a rectangular box. The spheres are allowed to settle in the box (Fig. 2). The simulation runs for 7,000 iterations and the performance is reported in terms of iterations per wallclock seconds. This standardized test (hence --stdperformance in the name) is constructed in such a way that almost all the computation happens on the C++ side, only the calculation of the wallclock time is done in Python. Obviously doing more calculations in Python will make any script slower. Hence, any calculation in Python should be kept to a minimum.
Since the benchmark results strongly depend on other processes running on the system, the test was performed at highest process priority after first making sure that all unrelated processes are stopped (via pkill -SIGSTOP command). The benchmark was repeated at least 50 times for each precision type and compiler settings. The average calculation speedx (in iterations per seconds) was determined for each precision type and data points not meeting the criterion 2σ < x i −x < 2σ (σ is the standard deviation) were considered outliers. On average 4% of data points per bar was removed, the largest amount removed was 10% (5 data points) on two occasions. Hence, each bar in Fig. 3   A summary of all the benchmark results is shown in Fig. 3 along with the relevant standard deviations. The performance is indicated in terms of iterations per seconds. The tests were performed 28,30 for the seven different precision types (Fig. 3A-G) listed in Tab. 3 for three different optimization settings: default cmake settings (Fig. 3a), with SSE vectorization enabled (Fig. 3b), and with maximum optimizations offered by the compiler but without vectorization 29 (Fig. 3c). The lack of significant improvement in the third case (Fig. 3c) shows that the code is already well optimized and the compiler cannot optimize it any further, except for long double and float128 types and the gcc compiler where a 2% speed gain can be observed ( Fig. 3Cc and Dc versus Ca and Da). It is interesting to note that the clang compiler systematically produced a code that runs about 4 to 9% faster than the gcc or icpc compilers. Intel compiler users should be careful, because the -fast switch might result in a performance loss of around 2 to 10%, depending on particular settings (Fig. 3c). Code vectorization (using the SSE assembly instruction set, an experimental feature, Fig. 3b) provides about 1 to 3% speed gain, however this effect is often smaller than the σ error bars. Enabling intel hyperthreading (HT) did not affect the results more than the standard deviation error of the benchmark. The float128 results for the intel compiler stand out with a 5% speed gain (Fig. 3D). However, not all mathematical functions are currently available for this precision in icpc and to get this test to work a crippled branch 30 was prepared for the tests with some of the mathematical functions disabled. The missing mathematical functions 31 were not required for these particular calculations to work. The clang compiler does not support 32 float128 type yet. The average speed difference between each precision is listed in Tab. 3. The run time increase with precision in the MPFR library is roughly O(N log(N )) (where N is the number of digits used) but it is application specific and strongly depends on the type of simulation performed [44,31].
It shall be noted that currently YADE does not fully take advantage of the SSE assembly instructions (cmake -DVECTORIZE=1) because Vector3r is a three component type, while a four component class Eigen::AlignedVector3 33 is suggested in the EIGEN library but it is not completely functional yet. In the future, this class can be improved in EIGEN and then used in YADE. 28 also see https://gitlab.com/yade-dev/trunk/-/tree/benchmarkGcc 29 the test with SSE and maximum optimizations was also performed but the results were simply additive, thus they were not included here. Also sometimes they produced the following error due to memory alignment problems: http://eigen.tuxfamily. org/dox-devel/group__TopicUnalignedArrayAssert.html, because the operands of an SSE assembly SIMD instruction set must have their addresses to be a multiple of 32 or 64 bytes, and the compiler could not always guarantee this. 30     (c) with maximum optimizations offered by the compiler (gcc, clang: -Ofast -fno-associative-math -fno-finite-math-only -fsigned-zeros § and additionally for native: -march=native -mtune=native; intel icpc -fast § ) but without vectorization. † on a PC with two Intel E5-2687W v2 @ 3.40GHz processors with 16 cores and 32 threads. ‡ via command echo off > /sys/devices/system/cpu/smt/control or by a BIOS setting (HT is also known as intel SMT). § the extra three flags are used by CGAL; in intel compilers -fast enforces all processor model native optimizations.

Problem description
A simple simulation of a triple elastic pendulum system was performed to check the effect of high precision in practice. Triple pendulums are considered highly chaotic as they provide an irregular and complex system response [3]. Numerical modeling of such systems can show the benefits of using high precision. A single thread was used (i.e. yade -j1) during the simulations to avoid numerical artifacts arising from different ordering of arithmetic operations performed by multiple threads. This allowed focusing on high precision and eliminating the non-deterministic effect of parallel calculations (e.g. arithmetic operations performed in a different order result in a different ULP error in the last bits [36,76,41]) and to have a completely reproducible simulation.
The numerical setup of the model represents the chain consisting of three identical elastic pendulums (see attached Listing 1 and on gitlab). The pendulums are represented by a massless elastic rod (a long-range normal interaction) and mass points. The latter are modeled using spheres with radius r = 0.001 m and density ρ = 1 kg/m 3 , noting that the masses of the spheres are lumped into a point. The rods are modeled by a normal interaction using cohesive interaction physics. The length of the chain is L = 0.1 m. Each rod is 1/30 m long and the normal stiffness of the interaction is k = 100 N/m. The strength (i.e. cohesion) is set to an artificial high value (10 7 N/m 2 ) so that the chain cannot break. Hence, the behavior of the rods can be assumed purely elastic. The initial position of the chain is α = −20 • relative to the horizontal plane (see Fig. 4a, t = 0 s). Gravity g = 9.81 m/s 2 is acting on the chain elements as they are moving. The process was simulated with time steps ∆t equal to 10 −5 s, 10 −6 s, 10 −7 s, and 10 −8 s. The results obtained using different precision are discussed in the following subsections. First, the evolution of the angles in the pendulum movement are discussed with ∆t = 10 −5 s. Second, the effect of damping is shown with ∆t = 10 −5 s. Then the effect of using various time steps ∆t is discussed. Finally, the total energy conservation is examined for various time steps.

Pendulum movement
Numerical damping was not used in this simulation series to avoid any energy loss and for the purity of the numerical results. The simulations were carried out with the different precisions listed in Tab. 5 and a time step of ∆t = 10 −5 s. Angles between the two rods were constantly monitored and saved with a period of 10 −4 s for further analysis and comparison. After the simulation was performed and the data was gathered, the Pearson correlation coefficient [88] was calculated for all data sets. The simulation with 150 decimal places (type mpfr 150) was used as the reference solution as it has the highest number of decimal places. The product of the two angles between the three rods was used as an input parameter for the calculation of the correlation coefficient. The data for the correlation was placed in chunks with each having 500 elements (10 −4 s ×500 = 0.005 s of the simulation). Then the scipy.stats.pearsonr function from SciPy [97] was   employed for the calculation of the Pearson correlation coefficient p. For further reference, the point in time when the correlation p between two simulations falls below p < 0.9 is marked as t s . The latter corresponds to the time from the start of the simulation until the correlation is lost and in the following it is denoted correlation duration. Fig. 4 shows snapshots of the evolution of the movement for three different precisions. At the beginning of the simulation the angles between the rods are the same. However, from a certain point in time onward the angles are starting to differ. Indeed, at t = 2.5 s (Fig. 4d), the green line (15 decimal places, type double) has clearly another state compared to the other two with 33 and 150 decimal places (types float128 and mpfr 150 respectively). The snapshot at t = 5.25 s (Fig. 4g) demonstrates the beginning of the deviation of the simulation with 33 decimal places (type float128). Thereafter all the pendulums are moving very differently and no correlation is observed. It can be concluded that using higher precision increases the time when accurate calculation results are obtained which is also reflected by the correlation duration t s listed in Tab. 5. Fig. 5 presents the correlation coefficient as a function of time. The graphs provide a more accurate representation of the point in time when the correlation disappears. One can see that there is a positive linear correlation with p = 1 at the beginning of the simulations. This means initially the rods are moving identically. The lowest precision curve float is starting to jump between correlation values of p = 1 and p = −1 at around t = 1.1 s (Fig. 5a), which means that no visible correlation is observed any more. For all the higher precision simulations the drop off happens later and progressively with increasing precision as summarized in Tab. 5.
Type double and long double have a correlation of p < 0.9 after 2.5 s (Fig. 5b, comparing 15 with 150 decimal places) and 3.1 s (Fig. 5c, 18 vs. 150 decimal places) respectively. The same tendency is seen from Boost float128 type (Fig. 5d, 33 vs. 150 decimal places) which deviates at approximately 5.1 s. Both simulations with 62 decimal places start to deviate at around 9.9 s. This clearly demonstrates that the level of precision, i.e. the number of decimal places, has an influence on the accuracy of the simulation results. Sometimes the decrease of the correlation happens suddenly and sometimes it starts to decrease slowly before it decreases rapidly. It can clearly be seen that the simulations with higher precision are showing better results and are closer to the reference solution calculated with 150 decimal places. Nevertheless, higher precision requires much more time for the simulation and more computing resources.

Effect of damping
To study the importance of precision in combination with other parameters, the same simulations as in Section 5.2 were carried out with a numerical damping coefficient equal to 5 · 10 −3 . Numerical damping is generally applied to dissipate energy. In this particular test case, numerical damping can be interpreted as the slowing down of the pendulum oscillation. The global damping mechanism was used, as described in the original DEM publication of Cundall and Strack [22]. Global damping acts on the absolute velocities of the simulation bodies and is implemented in the NewtonIntegrator class in the source code of YADE. Global damping slows down all affected bodies based on their current velocities.
The damping coefficient was chosen so that an effect of different precisions can be seen on the whole system. If the damping coefficient is too high (> 10 −1 ), the system loses its whole energy very quickly and no visible differences are seen. Too small damping coefficients (< 10 −3 ) lead to very slow energy dissipation. Fig. 6 shows the development of the angle between the first and the second rod of the pendulum during the simulation with a global damping coefficient equal to 5 · 10 −3 . One can clearly see the differences in simulation results based on different precisions. The float precision simulation (blue curve, Fig. 6a) indicates the largest deviation from all other simulations. Higher precisions gradually provide results that are closer to the simulation with the highest precision (boost mpfr 150 decimal places).
Since the angle between rods is oscillating rapidly, as can be seen on the inset in Fig. 6, the raw data was smoothed using the Savitzky-Golay filter [79] for better visibility. The filter removes most of the noise, and there are no visible differences between cpp bin float, mpfr 62 and mpfr 150. The three curves are basically overlapping each other.

Effect of time step ∆t
The time step ∆t is one of the most important parameters influencing the simulation results. Hence, simulations with time step values of 10 −5 s, 10 −6 s, 10 −7 s, and 10 −8 s were performed with different precisions. The values for the correlation duration were recorded and plotted in Fig. 7 for all precisions and time steps considered. It can be clearly seen that the correlation duration increases with increasing precision. This highlights once more that higher precision is required for higher confidence. It can also be seen that generally the correlation duration increases with decreasing time step. This is due the fact that a smaller time step results in a smaller integration error (e.g. the leapfrog integration scheme used in YADE has the error proportional to ∆t 2 per iteration). This tendency is also marked on the figure with a black arrow pointing in the direction of the smaller time step ∆t. The opposite result is observed for float. This is because 6 decimal places are not enough to work with ∆t = 10 −8 s. It shall be noted that the implementation of a more precise time integration scheme is out of the scope of the current paper. In the current symplectic leapfrog integration scheme, as implemented in the present version of NewtonIntegrator, the positions and velocities are leapfrogging over each other. There is no jump-start option which would allow to start the simulation with both position and velocity defined at t = 0. This means that the initial velocities are declared at t = − ∆t 2 and initial positions at t = 0 or initial velocities are declared at t = 0 and initial positions at t = ∆t 2 . Which of these two it is, is only a formal choice 34 , since both interpretations are valid. Therefore a more detailed comparison between the time steps is not carried out as the initial conditions for each simulation with a different ∆t differ slightly, i.e. the starting velocity declared in the script is interpreted as being defined at t = − 10 −5 2 or at t = − 10 −8 2 , thus resulting in slightly different simulations 36 .

Energy conservation
In the following, the total energy in the system is analyzed for different precisions and time step values of 10 −5 s, 10 −6 s, 10 −7 s, and 10 −8 s. No numerical damping is considered. The total energy in the system is calculated as the sum of the elastic energy in the interactions and the kinetic and potential energy of the mass points (i.e. spheres). As already pointed out in the previous section, the symplectic leapfrog integration scheme is used. This means that velocities and positions are not known at the same time. Hence, the velocities needed for an accurate calculation of the kinetic energy are taken as an average of the velocities from the current and the next iteration. All numerical results are compared to the reference solution which was calculated using 150 decimal places, similar as in the previous sections. Fig. 8 shows two typical results obtained from the study using a time step of ∆t = 10 −6 s. The results obtained with the other time steps have a similar trend and are not included for brevity. The top graphs show the evolution of the energy balance where each energy component is divided by the total reference energy to give an energy ratio. The total reference energy is calculated using 150 decimal places. The total energy ratio should be equal to 1 throughout the simulation. The bottom graphs show the absolute error in total energy calculated as the absolute difference between the total energy given by a specific precision and the constant total energy calculated using 150 decimal places. The results obtained using float are depicted in Fig. 8a. It can be seen that the absolute error in total energy starts to increase drastically after about 4 s. This is much later than the correlation duration. From the energy plot it can also be seen that energy is continuously added to the system from this time onward. This makes the simulation not only incorrect but also unstable. A different observation can be drawn from Fig. 8b where the results of a typical simulation with float128 are shown. It can be seen that the energy balance is stable and the absolute error is several order of magnitudes smaller. In addition, the absolute error does not have an increasing trend. Instead, it bounces around, i.e. it increases initially and decreases thereafter, and never goes above a certain threshold value. This is also a reflection of the symplectic leapfrog integration scheme. It should be noted that the results for double, long double, mpfr 62 and cpp bin float 62 are very similar to Fig. 8b and, hence, not shown for brevity. Fig. 9 summarizes the results for all precisions and all time steps. As noted previously, a detailed comparison between different time steps does not make sense because of the different initial velocities. Nevertheless, a qualitative comparison is valid. It can be seen that the maximum absolute error in energy balance for float is many orders of magnitudes larger than for the other precisions (note that the vertical axis uses logarithmic scale). The data also indicates that the error is almost constant for all other precisions. This clearly highlights the effectiveness and reliability of the symplectic leapfrog integration scheme implemented in YADE. However the error is many orders of magnitude larger than the ULP error of the higher precision types. For example to achieve maximum absolute error of 10 −30 [J] for float128, further decreasing the time step is not practical. Different approaches, such as higher order symplectic methods, have to be employed [68,69] 36 . Like in previous section, a smaller time step results in a smaller absolute error (this tendency is indicated by the black arrow), except for float where 6 decimal places are not enough to work with the smaller time steps.

Conclusions and future perspectives
The obtained results show that using high precision has a pronounced influence on the simulation results and the calculation speed. Higher precisions provide more accurate results and reduce numerical errors which in some fields can be beneficial. It can also be concluded that high precision is essential for research of highly chaotic systems (Fig. 5). Nevertheless, increasing the number of decimal places in the code leads to a higher CPU load and raises the calculation times (Fig. 3, Tab. 3).
Updating an existing software with a large codebase to bring the flexibility of arbitrary precision can be a challenging and error-prone process which might require drastic refactoring. A good test coverage of the code (unit and integration tests: the Continuous Integration pipeline in YADE) is highly recommended before beginning of such refactoring to ensure code integrity [10]. Also the employment of AddressSanitizer [83] is highly desirable to prevent heavy memory errors such as heap corruptions, memory leaks, and out-of-bounds accesses.
Simulations with the triple pendulum show that the results are starting to be different after a few seconds of the simulation time because of different precisions. The higher the precision the longer the results remain true to the highest precision tested. The same effect, within each precision ( Fig. 7 and Fig. 9), occurs when using smaller time steps ∆t, because it has a smaller time integration error per iteration. Applying damping can significantly smooth this effect (Fig. 6).
The new high-precision functionality added to YADE does not negatively affect the existing computational performance (i.e. simulations with double precision), because the choice of precision is done at compilation time and is dispatched during compilation via the C++ static polymorphism template mechanisms [89,96].
Concluding this work, the main modules of YADE now fully support two aspects of arbitrary precision (see Tab. 1 and Fig. 1): 1. Selecting the base precision of Real from Tab. 2 (which is an alias for RealHP<1>). 2. Using RealHP<N> in the critical C++ sections of the numerical algorithms (Section 2.4) 1 .
These new arbitrary precision capabilities can be used in several different ways in the management of numerical error [5]: 3. To periodically test YADE computation algorithms to check if some of them are becoming numerically sensitive. 4. To determine how many digits in the obtained intermediate and final results are reliable. 5. To debug the code in order to find the lines of code which produce numerical errors, using the method described in details in chapter 14 of [51]. 6. To fix numerical errors that were found by changing the critical part of the computation to use a higher precision type like RealHP<2> or RealHP<4> as suggested in [50].
The current research focus is to: 7. Add quantum dynamics calculations to YADE using the time integration algorithm which can have the error smaller than the numerical ULP error of any of the high-precision types: the Kosloff method [91,52,80,92] based on the rapidly converging Chebychev polynomial expansion of an exponential propagator 35 . 8. Add unit systems support, because there are also software errors related to unit systems, for example in 10 November 1999 the NASA's Mars Climate Orbiter was lost in space because of mixing SI and imperial units [87,74,43].
Possible future research avenues, opened by the present work, include: 9. Add more precise time integration algorithms. The problems mentioned in Section 5.4 are well known. Albeit symplectic integrators are particularly good [75,90], a better time integration method with smaller ∆t will be able to fully take advantage of the new high-precision capabilities (Section 5.5 and Fig. 9). There are three possible research directions: (a) Use the Boost Odeint library [24] with higher order methods 1,36 .
(b) Use the work on time integrators by Omelyan et al. [68,69]. These approaches suggests that it is potentially possible to decrease the run time more than 50-fold with the same computation effort upon switching to long double, float128 or higher types. The algorithms focus on reducing the truncation errors and eliminating errors introduced by the computation of forces. Such a smaller error allows to use a larger time step which will more than compensate the speed loss due to high-precision calculations. Of course, the standard considerations for the time step [17,70] would have to be re-derived. (c) Investigate whether the exponential propagator approach presented in [80] or in [72,82] could be used in YADE as a general solution for ODEs, similarly to [69,42,19,75,65,30], regardless if that is a classical or a quantum dynamics system. 10. Enhance all auxiliary modules of YADE for the use of high precision 37 . 11. Use the interval computation approach to reduce problems with numerical reproducibility in parallel computations by using boost::multiprecision::mpfi float as the backend for RealHP<N> type [76,66]. 12. Use different rounding modes to run the same computation for more detailed testing of numerical algorithms [51].
Overall, based on the presented work, the architecture of YADE now offers an opportunity to adjust its precision according to the needs of its user. A wide operating system support and simple installation procedure enable forming multidisciplinary teams for computational physics simulations in the Unified Science Environment (USE) [64]. This will expand the spectrum of tasks that can be solved, improve the results and reduce numerical errors. Of course, this option not only complicates the architecture and the source code, but also imposes a restriction on the choice of a programming language. Table 4: Maximum error after 2 × 10 7 function evaluations expressed in terms of Units in the Last Place (ULP) † , calculated by the absolute value of boost::math::float distance ‡ between functions from C++ standard library or boost::multiprecision when compared with its respective RealHP<4> (having four times higher precision) § .