Dynamics of the Hénon Map in the Digital Domain

Nonlinear maps are usually implemented using finite-precision floating-point formats both in practical applications and in theoretical investigations. In the digital domain, the size of the state space is finite and every trajectory after a finite number of iterations reaches a cycle. It is therefore important to study the influence of rounding errors and the finiteness of the state space on properties of nonlinear maps such as the number of cycles, their periods, sizes of basins of attraction of cycles, and average convergence times. In this work, a thorough analysis of the dynamics of finite-precision implementations of the Hénon map is carried out. Six computational formulas and three popular finite-precision floating-point formats are considered. An exhaustive search is performed to find all cycles existing for single-precision floating-point implementations. Interval methods are used to reduce the number of initial conditions that must be considered. An efficient graph-based algorithm is designed to find basins of attraction. For the double-precision and extended-precision implementations, statistical methods are utilized to find cycles and to estimate sizes of their basins of attraction. Properties of observed cycles and corresponding dynamical phenomena are thoroughly discussed.


I. INTRODUCTION
N UMERICAL simulations are one of the main tools used to study nonlinear maps. Nonlinear maps implemented in the digital domain are used in various applications, including pseudo-random number generators [1], [2], [3], [4], [5], secure communication [6], image encryption [7], [8], and chaotic cryptography [9], [10], [11]. Finite-precision implementations change dynamical properties of nonlinear maps and may promote various phenomena including the existence of short periodic attractors which are undesirable in practical applications. Hence, it is important to investigate the influence of computational accuracy and the computational formula used on properties of nonlinear maps.
Various methods to improve properties of chaotic maps implemented in the digital domain have been proposed. This includes methods to increase the length of generated trajectories before repeating itself and methods to improve randomness of generated trajectories. A varying parameter compensation method is introduced in [1]. A reseeding-mixing method to extend periods of observed cycles and to improve Manuscript  the statistical properties of generated trajectories is carried out in [12]. Dynamic parameter-control chaotic system framework in which output of one chaotic map is used to control the parameter of another chaotic map, is introduced in [3]. Cascade chaotic systems are used to improve the randomness and security level of pseudo-random number generators in [13]. A congruential generator based method is proposed to counteract the performance degradation of digital chaos in [14]. When a nonlinear map is implemented in the digital domain the state space is finite. As a consequence, each trajectory, after a finite number of steps, must visit a point in the state space that was visited before. From then on, the trajectory remains periodic. The part of a trajectory before reaching the steady-state periodic solution is called a transient. Assuming that the size of the state space is M, the maximum length of a periodic trajectory that may be observed is M. In practice, observed trajectories have much shorter periods. The eventual periodicity of all trajectories has important consequences from both the theoretical and practical points of view. The problem is especially important for non-linear maps that exhibit chaotic behavior. Under certain assumptions, a chaotic attractor contains an infinite number of unstable periodic orbits. The measure of the set of unstable periodic orbits is zero, which means that there is a zero probability of selecting an initial point with periodic behavior. An arbitrarily small deviation from the position of an unstable periodic orbit results in the divergence of the trajectory from this orbit. However, when calculations are carried out in finite precision, it may happen that some of these unstable periodic orbits could actually be observed in simulations. This is caused by a nonzero probability that a trajectory initiated close to the position of an unstable periodp orbit returns to the initial point after a number of iterations, which is equal to p or its low multiple.
The key properties of chaotic map implementations are the lengths of observed cycles and the statistical properties of generated trajectories. The degradation of dynamical properties of nonlinear maps caused by finite-precision implementations has been studied by many researchers. The influence of rounding errors on dynamical behaviors of the logistic map is studied in [15]. The authors study the length of the longest limit cycle and the average transient length as a function of the size of the state space. It is shown that for the logistic map, in simulations, one observes cycles of length of the order √ M, where M is the size of the state space. It is argued that although periodic trajectories are relatively short, different types of behavior are correctly recognized in numerical simulations, and that certain statistical properties of observed trajectories, such as Lyapunov exponents, are preserved. The relation between the expected period of observed orbits and the maximum rounding error for the Ikeda map is studied in [16]. Dynamical degradation of This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ piecewise linear chaotic maps caused by finite-precision implementations is studied in [17]. Statistical study of the influence of double-precision errors on the behavior of the logistic map is carried out in [18]. The influence of the rounding method on properties of trajectories generated by the logistic map is studied in [19]. Five rounding modes are considered, and it is shown that periods of cycles and transient lengths depend on the rounding mode. In [9], fixed-point arithmetic is used to study cycle lengths of finite-precision implementation of the logistic map and the influence of correlation properties on the performance of chaos-based cryptosystems is discussed.
Properties of finite-precision implementations of nonlinear maps may be studied using the idea of state transition diagrams [15]. A similar approach under the name of state-mapping networks is used in [20] to analyze dynamical behaviors of the logistic map and the tent map implemented using fixed-point arithmetic. Analysis of dynamical properties of the logistic map implemented with up to 16 bits of precision is carried out in [14]. Dynamical properties of the logistic map implemented using single-precision and double-precision computations are analyzed in [21]. The graph-based approach is used in [22] to analyze the dynamics of the logistic map.
In principle, to find all cycles existing in a finite-precision implementation of a nonlinear map one should compute trajectories initiated in all points belonging to the state space. When the size of the state space is small and all states can be stored in a computer memory, then a complete characterization of cycles and their basins of attraction can be carried out using the graph-based approach [14], [15], [20], [21]. In this approach, one selects a point in the state space and computes its trajectory until a cycle is found or a point visited before is reached. After all points are visited, the graph structure is known and cycles and their basins of attraction can be easily identified. This approach has been used to study the logistic map and the tent map implemented using fixed-point arithmetic and low precision floating-point arithmetic. Similar calculations have been carried out for the single-precision floating-point implementation of the logistic map in [21]. To find all cycles in larger state spaces one may use a trajectory-based approach, where the exhaustive search is carried out in a reduced search space. This method has been applied to the logistic map implemented in doubleprecision (compare [21]). Most of the results available in the literature involve nonlinear maps of dimension one and lowprecision accuracy. For higher dimensional maps, a complete characterization of the dynamics of the map implemented in the digital domain is a challenging research problem even for the single-precision accuracy.
In this work, we develop methods which may be used to study higher-dimensional maps implemented using floatingpoint formats. As an example, we study dynamical properties of Hénon map [23]-a well-known example of a chaotic twodimensional map-implemented using three popular binary floating-point formats: the single-precision (binary32) format, the double-precision (binary64) format, and the extendedprecision format. Two selections of parameter values are considered: the classical parameter values for which the famous chaotic attractor is observed in simulations and parameter values for which a stable period-28 orbit exists. Due to the size of the state space, one cannot use methods developed for one-dimensional systems to analyse the Hénon map for any of the considered floating point formats. For the single-precision format, a new method is proposed to find all cycles existing in the system and to construct basins of attraction in case of small basins. For the double-precision and the extended-precision format, statistical approach is used to find cycles and sizes of their basins of attraction. Properties of observed trajectories are discussed and it is verified, whether observations on the number of cycles and their average length reported in the literature [15], [16], [20] are valid for the Hénon map.
The main contribution and the novelty of this work is the development of methods to study the dynamics of nonlinear maps in the digital domain and carrying out analysis of the Hénon map implemented using different computational formulas in most popular floating point formats.
All computer programs are developed in C++ language and compiled using the g++ compiler, version 9.4.0. For the interval arithmetic support the CAPD library [24] is used. Computation times are reported for a single core 3.4 GHz processor. Parallel computations are used to reduce the wall computation time. The C++ code implementing algorithms presented in this work and the data necessary to carry out the computations are available at http://www.zet.agh. edu.pl/henon/tcasi2023.
The layout of the remaining part of the paper is as follows. In Section II, the definition of the Hénon map is recalled and its finite precision implementations are described. Computational tools for the analysis of finite-precision implementations of nonlinear maps are presented in Section III. Analysis of finite-precision implementations of the Hénon map for two sets of parameter values is carried out in Section IV. Conclusions are presented in Section V.

II. THE HÉNON MAP IN FINITE PRECISION
The Hénon map [23] is a two-parameter map of a plane defined by The linear change of coordinates (x, y) → (ȳ, bx) converts (1) to Digital implementations of maps (1) and (2) are equivalent, provided that the order of elementary arithmetic operations is the same. In the following, we use the latter version because it permits a more efficient search algorithm for cycles. We consider six computational formulas for equation (2): h(x, y) = (y, 1 + bx − ayy), h(x, y) = (y, bx − ayy + 1), h(x, y) = (y, bx − yya + 1), h(x, y) = (y, 1 − ayy + bx).
(3f) In the following, we consider three finite-precision floatingpoint formats: the single-precision binary floating-point format (binary32, float in C++), the double-precision binary floating-point format (binary64, double in C++), and the extended-precision format (long double in C++). The first two types are defined in the IEEE 754 floating-point computation standard [25] and are commonly used in microprocessors. Some important properties of these representations are collected in Table I: m 10 is the number of decimal digits of precision for a given data type, m 2 specifies the number of base 2 digits in the mantissa part, and ε is the smallest positive number such that the arithmetic operation 1 + ε carried out using a given data type returns a value larger than 1.
In the last row of Table I computation speeds are compared. For each data type, the time t needed to compute 10 8 iterations of the Hénon map is reported. Computations carried out using the first two data types take the same time. The singleprecision format is used primarily to save memory space and not to gain speed. Note that careful implementation of computational formulas in C++ is needed when using the float type. In C++ the constants are of double type by default and therefore the constant 1.0F should be used instead of 1.0 in the formulas (3). If the constant 1.0 is used, then all variables are cast to the double type, which takes extra time, and the computations are slower. Computations using the long double format are more than two times slower than those for the first two formats.

III. ANALYSIS OF NONLINEAR MAPS IN THE DIGITAL DOMAIN
In this section, a set of tools for studying nonlinear maps implemented using finite-precision arithmetic is presented. We describe methods to estimate the size of the state space, compute the inverse of a map, find all cycles, their basins of attraction, and average convergence times (transient lengths). The methods are illustrated using the Hénon map h : R 2 → R 2 as an example. After minor adjustments, the methods can be applied to other nonlinear maps.
In the digital domain, the number of representable points in the state space is finite. For an n-dimensional map a representable point is a point in R n which can be exactly represented using a selected finite-precision arithmetic. The dynamics of a finite precision implementation of a map can be represented as a graph with the vertices being representable points and the edges being transitions between these points as defined by the action of the map. The graph has a structure of c connected components, where c is the number of cycles. Each connected component is the basin of attraction of a single cycle. In the characterization of finite-precision implementations of maps, the most important problems are the number of cycles and their properties such as the period, the size of the basin of attraction, and the average convergence time.
The basin of attraction of a given cycle is the set of initial points such that trajectories based at these points after a finite number of iterations reach this cycle. The probability c prob of reaching a given cycle starting from a random initial condition belonging to the state space can be computed as the ratio of the basin size s and the size of the state space M.
The convergence time for a given representable point is defined as the number of iterations needed to reach the steady state. For a point that belongs to a cycle, the convergence time is zero. The average convergence time t aver and the maximum convergence time t max for a given cycle are defined as the average and maximum convergence times calculated over the whole basin of attraction of this cycle. The average (maximum) length of a trajectory before repeating itself can be computed as the sum of the period p of the steady-state cycle and the average (maximum) convergence time t aver (t max ).

A. Size of the State Space
where k x and k y are integers and δ > 0 is a real number. For a fixed δ each box v k , is uniquely defined by integers k x and k y . Next, for each box v k an enclosure of its image is found using interval arithmetic computations. The result obtained is an interval vector containing images of all points belonging to the box v k . Boxes v j which have non-empty intersection with the enclosure are identified. This process creates a graph representation of possible transitions between boxes. In this representation, boxes are graph vertices and non-forbidden transitions are graph edges. The reduction of the size of the state space is obtained by removing boxes containing transient dynamics only, which is equivalent to finding an enclosure of the invariant part of the map. In this procedure, all boxes which are not a beginning of any edge are removed. Similarly, we remove all boxes which are not an end of any edge. A detailed description of this procedure is given in [26]. Decreasing the box size δ increases the number of boxes in the covering of the invariant part but decreases its area and hence decreases the number of representable points belonging to these boxes. The value of δ should be selected as small as possible under the condition that the covering can be stored in the computer memory. Using this procedure with boxes of size δ = 2 −17 yields the covering E = {v k } N k=1 of the invariant part composed of N = 25736445 boxes (compare also [27]) with the area below 0.0014 which is significantly smaller than the area of the rectangle [−1.5, 1.5] × [−1.5, 1.5]. Let us denote by the set of representable points belonging to the boxes in the covering E. In the following, the set is considered as the state space of finite-precision implementations of the Hénon map for all computational formulas (3). The size of is M f ≈ 2.6 · 10 13 and M d ≈ 5.9 · 10 31 for single-and double-precision representations, respectively. For a single box the number of representable points belonging to this box is calculated as a product of representable numbers in the x and y coordinates. M f and M d are computed as the sum of representable points in each box over all boxes belonging to E. Note that the number of representable points in a given box is not constant despite the fact that the boxes are of the same size. For example, for the single-precision representation, the number of representable points belonging to a box in the covering E varies from 4225 to 1.2 · 10 11 . This huge difference is a consequence of using subnormal numbers defined in the IEEE 754 standard. With subnormal numbers, the density of representable numbers in a neighborhood of zero is significantly increased.

B. Finding All Cycles
As mentioned before, in order to find all cycles, one should, in general, consider all points belonging to the state space and for each initial point find the cycle to which its trajectory converges. In the case of the Hénon map, the size of the state space is very large even for the single-precision floating-point format (M f ≈ 2.6 · 10 13 ). The number of initial points to be considered can be reduced by skipping certain boxes. The idea is to introduce an ordering of boxes and to skip a box if all trajectories passing through this box has to pass through boxes larger according to the selected ordering. An ordering should permit removing boxes containing many representable points (i.e., those that lie close to the axes x = 0 and y = 0). The following ordering is used during the computations: v k v l if and only if k x > l x or k x = l x and k y > l y , where (k x , k y ) and (l x , l y ) are pairs of integers defining boxes v k and v l , respectively. The box v k is skipped if its image or preimage under the mth iteration of the map is enclosed in the set of boxes E k = {v ∈ E : v v k }. This elimination procedure does not require recalculation of images or preimages of boxes. It can be carried out based on the information on possible transitions between boxes which is stored in the graph representation of the dynamics of the map. Applying this procedure to the set of boxes E with m ∈ {1, 2, . . . , 10} leads to the set E composed of N = 1183485 boxes. Using m > 10 increases the computation time without significant improvement of the results. The number of representable points in a single box belonging to E varies from 4225 to 33153. The number of initial points that must be considered c ← the cycle for the initial condition (x, y) 19: end for 21: end for 22: return C 23: end function to find all cycles is reduced from M f ≈ 2.6·10 13 representable points belonging to boxes in the covering E to M f ≈ 6.5 · 10 9 representable points in E . In the last step for each initial point in E the corresponding steady state cycle is computed. The complete procedure to find all cycles is presented as the Algorithm 1.
For the double-precision and extended-precision computations the presented method to find all cycles does not work due to very large sizes of the state space. Cycles can be found by computing steady state trajectories starting from random initial conditions. This method does not guarantee finding all cycles and will be referred to as the statistical approach. Cycles with a large basin of attraction (which is equivalent to a high convergence probability), are more likely to be found. The probability that a cycle with the convergence probability c prob is detected starting from N random initial conditions is 1 − (1 − c prob ) N . To detect such a cycle with the probability larger than 0.95 one should use N such that 1 − (1 − c prob ) N > 0.95. It follows that with N random initial conditions a cycle with the convergence probability c prob ≥ c N = 1− N √ 0.05 is detected with the probability larger than 0.95. For example, for N = 100 and N = 10000 we have c N ≈ 0.02951 and c N ≈ 0.0002995, respectively.

C. Computing Preimages in the Digital Domain
The (infinite precision) Hénon map (2)  return (y,1+b*x-y*y*a) 3: end function 4: function h −1 (x, y) evaluate the inverse 5: return ((y-1-a*x*x)/b,x) 6: end function 7: function FINDPREIMAGE(x , y) 8: S ← ∅ initialize the preimage 9: 10: repeat 11: if (x , y ) = (x, y) then 13: S ← S ∪ {(x , y )} update the preimage 14: end if 15: x ← the smallest representable number above x 16: until y ≤ y 17: repeat 19: if (x , y ) = (x, y) then 21: S ← S ∪ {(x , y )} update the preimage 22: end if 23: x ← the largest representable number below x 24: until y ≥ y is found using the selected computational formula. If z = z, then z is added to the preimage h −1 (z). Next, x is increased to be the next representable number, and the point z = (x , y ) is added to h −1 (z) if z = (x , y ) = h(z ) = z. Increasing x and adding new points to h −1 (z) is continued as long as y ≤ y (the stopping condition is y > y). A similar procedure is carried out with decreasing x starting from x = (y −1−ax 2 )/b and using the stopping condition y < y. A pseudo code for the algorithm to compute the preimage of the version (3a) is presented as the Algorithm 2. Other versions are obtained by changing the line 2 in the algorithm.
The time needed to find the preimage using the Algorithm 2 depends linearly on the preimage size. For the Hénon map, the size of the preimage in case of single-precision computations varies from zero (empty set) to approximately 1.77 · 10 9 . For example, the preimage of z = (0, 1.25) contains 7 points and the preimage of (1.25, −1.1875) = h(z) = h(0, 1.25) contains 1760908629 points. Extremely large preimage sets exist when the preimage of z contains a point (x , y ) with the first coordinate being zero or close to zero. In this case, many points of the form (x , y ) with x being a subnormal number produce the image z.

D. Constructing Basins of Attraction
The basin of attraction of a cycle may be constructed in a recursive manner by finding preimages of points belonging to the cycle, their preimages, and so on until there are no preimages in the state space. The procedure to compute the basin of attraction is as follows. First, we compute the preimage P of the cycle, that is, the set of points outside of the cycle whose images belong to the cycle. For each point z ∈ P its basin of attraction (the set of points whose trajectories reach z) has a directed rooted tree structure with z being the root. The basin of attraction of z ∈ P can be found using the Depth First Search (DFS) algorithm [28]. The DFS algorithm is continued until a leaf (a point without preimages) is reached or a preimage is outside the state space . The advantage of this procedure is that we do not need to store the basin of attraction in a computer memory. As a consequence, the procedure can be applied to relatively large basins. In Section IV, it is shown that the proposed approach can successfully handle basins of size up to 2 · 10 11 .
The sizes of larger basins can be estimated using a statistical method. In this approach, a large number of initial points in the state space is selected and for each initial point the steady state solution is identified. The probability c prob of convergence to a given cycle is calculated as the ratio of the number of initial points whose trajectories converge to this cycle and the total number of initial points. The size of the basin can be estimated as s = c prob · M where M is the number of representable points in .
In infinite precision, to distinguish whether a cycle represents a stable orbit, one may compute the Lyapunov spectrum of a trajectory along this cycle. For a two-dimensional map, the Lyapunov exponents λ 1 , λ 2 computed along a cycle characterize the rate of attracting or repelling nearby trajectories in different directions of the state space. We assume that the Lyapunov exponents are sorted in descending order, i.e., λ 1 ≥ λ 2 . Note that the Lyapunov spectrum computed for different cycles is usually different. If all Lyapunov exponents are negative, then the cycle is stable. This indicates that the cycle corresponds to a periodic attractor of the system. If the largest Lyapunov exponent is positive, then the orbit is unstable. In this case, we may expect that the cycle corresponds to a chaotic trajectory and the observed periodic trajectory is an artifact caused by rounding errors. For the calculation of Lyapunov exponents, the QR decomposition-based method is used. During the computations a trajectory (z 1 , z 2 , . . . , z n ) of the length n is considered. The matrix Q 0 is initialized to be the identity matrix. At the kth step the Jacobian matrix J k = h (z k ) is computed and the matrix J k Q k−1 is decomposed as a product Q k R k of an orthogonal matrix Q k and an upper triangular matrix R k . The i th Lyapunov exponent is computed as the sum of logarithms of the i th diagonal entries of matrices R k divided by n (for details see [29]). When Lyapunov exponents are calculated along a cycle, it is sufficient to consider a single pass of a trajectory along this cycle. To obtain accurate results, the QR decomposition algorithm should be applied for several iterations before starting calculations of the Lyapunov exponents. This ensures that the matrices Q k are properly aligned for all points along the considered trajectory. From the Lyapunov spectrum (λ 1 , λ 2 ) one may compute the Lyapunov dimension [30] where k = max j : For two-dimensional systems, if λ 1 < 0 then the Lyapunov dimension is zero. Otherwise, the Lyapunov dimension is not less than one.

IV. DYNAMICS OF FINITE-PRECISION IMPLEMENTATIONS
OF THE HÉNON MAP: CASE STUDIES In this section, we study the effects of finite-precision implementations of the Hénon map on the dynamics of the system. We consider two sets of parameter values for which in simulations one observes a chaotic trajectory and a stable periodic trajectory, respectively.

A. Chaotic Case
First, consider the case with classical parameter values a = 1.4 and b = 0.3. for which the famous Hénon attractor is observed in simulations (see Fig. 1). The question of whether the Hénon attractor is chaotic remains an open problem [31].
Parameter Let us start with the single-precision format. Using the method presented in Section III all cycles existing for the computational formulas (3) have been found. The results obtained for different versions v of the computational formula are presented in Table II. For each cycle, we report its period p and the position of the periodic point (x, y) with the smallest x coordinate to permit verification of the results. We also report Lyapunov exponents λ 1,2 computed along the orbit and the Lyapunov dimension d. To evaluate the performance of the statistical approach N = 500000 initial points belonging to the covering E are selected randomly and for each initial point its steady state is found. The average computation time is below 0.01 second per initial point. For each cycle found using the statistical approach, we report the probability c prob of convergence for this cycle. The size s of the basin of attraction is reported for small basins for which the graph-based approach with the DFS algorithm is successful. For large basins the basin size can be approximately calculated as s ≈ c prob · M where M is the size of the state space. t aver and t max are the average and maximum convergence times (i.e., the average and the maximum number of iterations needed to reach a given cycle). Their values are exact in cases where the graph-based approach is successful in finding the basin of attraction. For larger basins, the values of t aver and t max are estimated using the statistical approach.
First, let us note that each computational formula leads to different results. It follows that the choice of the computational formula has a significant influence on the dynamics of finite-precision implementations. The number of cycles varies from 5 to 15. The shortest period is between 1 and 32. The longest cycle has a period between 14119 and 55754. Note that several low-period orbits are not found using the statistical approach (the field c prob is empty).
The graph-based approach is successful in finding the basin of attraction if the basin size is below 2.1 · 10 11 , which corresponds to the probability of convergence c prob ≤ 0.007. The size of the state space can be estimated by dividing the basin size s by the convergence probability for the cycles with both numbers reported. Computing this estimate for cycles with basin size greater than 10 9 gives results in the range [10 13 , 5 · 10 13 ], which is consistent with the estimate M f ≈ 2.6 · 10 13 obtained in Section III.
The longest cycle is predicted to be reached from the majority of initial conditions (compare [15], [16]). For example, in [15] it is observed that, typically, 90% of the states evolve to the longest cycle and that the number of distinct cycles is less than 5. The results presented in Table II only partially confirm these hypotheses. In two out of the six cases (versions (3c) and (3e)) the longest cycle does not have the largest basin. For the version (3c) the relatively short period-619 orbit has the largest basin and the convergence probability above 82%. In three cases, the probability of convergence to the cycle with the largest basin is below 80%. Moreover, the number of cycles for the version (3a) is 15 including 11 cycles found using the statistical approach. Another interesting case is observed for the version (3f), where the very short period-13 orbit has a relatively large basin of attraction with the size above 2·10 8 . This orbit has been found using the statistical approach. It follows that with a non-negligible probability, one may encounter short cycles starting from random initial conditions.
Let us now discuss the hypothesis that the length p of the longest cycle is of the order of √ M (compare [15]). This is equivalent to stating that M is of the order of p 2 . With the length p ≈ 5·10 4 we obtain M ≈ p 2 = 2.5·10 9 . This is much The discrepancy is likely to be caused by the fact that a large number of points in the state space (for example, points with a coordinate being a subnormal number) are not reachable (i.e., have an empty preimage). In this context, a more precise estimate of the size of the state space is based on the set E which contains approximately 6.5 · 10 9 representable points. For all cycles, the largest Lyapunov exponent and the Lyapunov dimension are positive. In most cases, the largest Lyapunov exponent is close to the value of the largest Lyapunov exponent λ 1 ≈ 0.419 observed for higher precision computations. The largest deviations are seen for short periodic orbits: λ 1 ≈ 1.182 for the period-1 cycles existing for versions (3a) and (3e), λ 1 ≈ 0.654 for the period-1 cycles existing for versions (3a) and (3b), λ 1 ≈ 0.5391 for the period-4 cycle (version (3f)), and λ 1 ≈ 0.1988 for the period-13 cycle (version (3f)). These cycles correspond to short unstable periodic orbits of the infinite-precision Hénon map. Their presence in finite-precision implementations is a result of rounding errors. The small maximum Lyapunov exponent λ 1 ≈ 0.1988 observed for the period-13 cycle (version (3f)) means that this cycles is weakly repelling. This explains the large basin of attraction of this cycle when compared to other low period cycles. In all cases, the average length of a trajectory before it repeats itself belongs to the interval of [5 · 10 4 , 8 · 10 4 ]. It follows that for the float type it does not make sense to consider trajectories with more than 10 5 iterations.
Let us now consider the double-precision floating-point format. In this case, the graph-based approach does not work due to extremely long computation time needed to consider all initial conditions. Therefore, we are limited to the statistical approach. For the sake of brevity, only the results for versions (3a) and (3f) are presented. For each version N = 3000 initial conditions are considered and steady-state behavior is found. The average computation time is 4 hours per initial point. Computation times are much longer than for the single-precision case due to much longer transients. The results are presented in Table III.
The number of cycles found is 6 and 5 for versions (3a) and (3f), respectively. The periods of the longest cycles are p = 2875179482 and p = 7516037528, which corresponds to estimates of the size of the state space p 2 ≈ 8.3 · 10 18 and p 2 ≈ 5.6 · 10 19 . These estimates are significantly smaller than the number M d ≈ 5.9 · 10 31 of representable points in the covering E. It is expected that the difference is related to the usage of subnormal numbers, similarly as for the singleprecision floating-point format. For both versions, the longest cycle has the largest basin of attraction. The probability of convergence to the longest cycle is c prob = 0.864 and c prob = 0.988 for versions (3a) and (3f), respectively. The probability of convergence to the second longest cycle for the version (3a) is c prob ≈ 9%. For the version (3f) the convergence probability for each cycle apart from the longest one is below 1%. In general, the probability of convergence is lower for shorter orbits. There are, however, some exceptions. For example, for the version (3a) the probability of convergence to the orbit #4 is more than twice that of orbit #3, despite the fact that the period is almost two times shorter.
Values of Lyapunov exponents are very close to each other for all cycles. This means that each cycle represents a chaotic trajectory and can be used to characterize behavior of the infinite-precision system in terms of Lyapunov exponents and the Lyapunov dimension. The average length of a trajectory before repeating itself is above 10 10 .
One should keep in mind that although short cycles have not been found, it is very likely that there exist short cycles with periods below 100 as in the single-precision case. Such orbits are difficult to find using the statistical approach due to small basins of attraction.
The results obtained for the extended-precision (long double) format are shown in Table IV. The number of initial conditions considered for each case is N = 250. The average computation time is 4 and 2 days per initial point for versions (3a) and (3f), respectively. The longer average computation time for version (3a) is the effect of longer convergence times. The lower number of cycles found (2 and 3 for versions (3a) and (3f), respectively) compared to the double-precision case is related to the lower number of initial conditions considered. The results regarding Lyapunov exponents and the Lyapunov dimension are practically the same for all orbits. The average convergence time is above 10 12 and is 100 times longer than for the double-precision case. It follows that the extended-precision format is sufficient for most applications since longer trajectories are usually not needed.

B. Periodic Case
Let us now consider the parameter values a = 1.399999999 99999991118, b = 0.29999997749050000273 for which a stable period-28 orbit exists (compare [32]). Fig. 2 shows a short part of a chaotic transient trajectory (blue dots) and the periodic attractor (red circles). The coordinates of the point   The cycles existing for single-precision implementation found using the graph-based approach are reported in Table V. For the version (3a) the number of cycles is 10 and the longest cycle has period p = 26581. The convergence probability for this cycle is relatively low c prob ≈ 55.2%. Seven cycles are found using the statistical approach with N = 500000 initial points. There are two fixed points (period-one orbits) located very close to the position of the true fixed point Their basins of attraction are very small. For the version (3f) the number of cycles is 8 including 6 cycles found using the statistical approach. The longest cycle with the period p = 184481 attracts more than 95% of initial conditions. The shortest period-2 orbit is located very close to the position of the fixed point x = y ≈ 0.631354477.
Note that period-28 orbit existing for the infinite-precision implementation is not detected. None of the orbits found is located close to the position of the stable period-28 orbit. The reason is that the single-precision implementation does not offer sufficient accuracy to observe this orbit. The point (a, b) = (1.399999976, 0.299999982) in the parameter space used in the computations does not belong to the periodic window for which the stable period-28 orbit exists. If follows that single-precision computations cannot be used to study properties of the considered period-28 stable orbit.
The results obtained for the double-precision implementation using the statistical approach with N = 1000 initial points are presented in Table VI. The average computation time is 40 minutes per initial point. For the version (3a) the majority of initial conditions are attracted to the period-532 orbit. Note that p = 532 is a multiple of 28 and the position of the cycle is close to the position of the stable period-28 orbit. The Lyapunov exponents computed along this orbit are both negative. Almost 20% of the initial conditions are attracted to long cycles with a positive Lyapunov exponent. Even more interesting phenomena are observed for the version (3f). In this case, more than 56% of the initial points are attracted to long cycles, which behave in a chaotic manner. The largest Lyapunov exponent calculated along these cycles is positive, indicating that the behavior of the map is chaotic. Approximately 43% of initial points are attracted to period-224 and period-196 orbits with negative Lyapunov exponents. The periods p = 224 and p = 196 are multiples of 28, and the orbit positions are very close to the position of the stable period-28 orbit. In all cases, the average convergence time is above 10 9 , which means that we need to wait long time before we can decide what type of steady-state behavior exists.
From the results presented in Table VI it follows that the observed chaotic trajectory does not necessarily indicate that the underlying dynamics of the infinite precision system is indeed chaotic. The reason is that the transient times needed to converge to a stable periodic orbit may be very long. For long transients there is a non-negligible probability that before convergence a transient trajectory hits a point in the state space which was visited before, thus forming a long periodic trajectory not related to the existing periodic attractor.
The results obtained for the extended-precision implementation based on the calculations of trajectories started in N = 1000 randomly selected initial points are reported in Table VII. The average computation time is 5 hours per initial point and is significantly smaller than for the chaotic case, which permits testing more initial points. All tested trajectories are attracted to a cycle located very close to the true period-28 orbit. All periods are equal to 28 or are low multiples of 28. For the versions (3a) and (3f) there are two and three distinct period-28 orbits, respectively. Lyapunov exponents computed along the cycles are all negative and are equal within the reported accuracy (four significant decimal digits). Average convergence times are below 10 10 for all cycles, which is approximately 100 times shorter than for the chaotic case.
We may conclude that the extended-precision implementation is better suited to study dynamical behaviors of the Hénon map for parameter values a = 1.4, b = 0.2999999774905. Due to the larger size of the state space, the probability that a transient trajectory hits a cycle before converging to the stable periodic orbit is lower than for the double-precision case.
On the other hand, one may say that the use of a lower precision introduces more chaos in the system (rounding errors are larger). For the double-precision implementation some of the observed cycles resemble chaotic trajectories with positive Lyapunov exponents, while for the extended-precision all steady-state trajectories are low-period cycles with negative Lyapunov exponents.

V. CONCLUSION
A set of tools was presented for the analysis of finite-precision implementations of nonlinear maps. A method to reduce the search space to find all cycles was proposed. An efficient algorithm to compute preimages of a given representable point was described. A DFS based algorithm was designed to find basins of attraction of cycles.
The proposed methods were used to analyze finite-precision implementations of the Hénon map for chaotic and periodic cases. It was shown that different computational formulas lead to different dynamical properties of finite-precision implementations in spite of the fact that these formulas are equivalent from the mathematical point of view. It follows that when results of computer simulations of nonlinear maps are reported, it is necessary to precisely describe the way the map is implemented to permit reproducibility of the results. Different computational formulas and precisions should be tested to make sure that the observed phenomena are not artifacts caused by rounding errors. For the single-precision floatingpoint implementation, all cycles were found. The structure of the state space was investigated. Basins of attraction of cycles were calculated using the graph-based approach for basins with sizes below 2 · 10 11 . Sizes of larger basins were estimated using the statistical approach. For double-precision and extended-precision implementations, cycles and their basins of attraction were studied using the statistical approach. Relations between cycle lengths, precision used, and the size of the state space were studied. It was shown that finite-precision implementations of nonlinear maps may support cycles with very low periods. This may degrade properties of non-linear maps in practical applications.
The results show that dynamical properties of finiteprecision implementations, such as the number of cycles, their lengths, and the probability of convergence to a given cycle, depend on the order of elementary operations used to evaluate the map. It follows that the selection of the computational formula used in finite-precision computations is important both for applications of nonlinear maps and for the results of theoretical studies. In applications of nonlinear maps one should test all possible evaluation methods and select the one with best dynamical properties.
From the results obtained for the double-precision implementations it follows that rounding errors may strongly disturb the simulation results. Chaotic-like trajectories observed in simulations do not necessarily imply that the dynamics of the infinite-precision system is chaotic. The results presented should alert researchers studying nonlinear maps to dangers of making conclusions about dynamical properties based purely on computer simulations.
The methods presented and the results obtained for the Hénon map can be used by other researchers to help design simulation studies and applications of nonlinear maps. The methods after minor modifications may be used to study other nonlinear maps implemented in the digital domain.