From Continuous-Time Chaotic Systems to Pseudo Random Number Generators: Analysis and Generalized Methodology

The use of chaotic systems in electronics, such as Pseudo-Random Number Generators (PRNGs), is very appealing. Among them, continuous-time ones are used less because, in addition to having strong temporal correlations, they require further computations to obtain the discrete solutions. Here, the time step and discretization method selection are first studied by conducting a detailed analysis of their effect on the systems’ statistical and chaotic behavior. We employ an approach based on interpreting the time step as a parameter of the new “maps”. From our analysis, it follows that to use them as PRNGs, two actions should be achieved (i) to keep the chaotic oscillation and (ii) to destroy the inner and temporal correlations. We then propose a simple methodology to achieve chaos-based PRNGs with good statistical characteristics and high throughput, which can be applied to any continuous-time chaotic system. We analyze the generated sequences by means of quantifiers based on information theory (permutation entropy, permutation complexity, and causal entropy × complexity plane). We show that the proposed PRNG generates sequences that successfully pass Marsaglia Diehard and NIST (National Institute of Standards and Technology) tests. Finally, we show that its hardware implementation requires very few resources.


Introduction
Many engineering applications require the utilization of random numbers, such as in the area of communication, encryption, codification, and modulation [1][2][3].
The use of chaotic systems, such as Pseudo-Random Number Generators, has recently grown because of the multiple advantages they present over stochastic algorithms [4][5][6][7]. In [8], the authors propose a five-step encryption algorithm. One of these parts is a chaotic systems module, where the system chooses between different number generators. In the conclusions, the authors highlight the importance of its randomness and present digital degradation as a subject to study. It is well known that both chaotic maps and continuoustime chaotic systems have internal correlations (they can be easily seen in 3D plots of their outputs) that prevent them from being used as PRNGs. The picture is even worse in the case of continuous-time chaotic systems that, unlike chaotic maps, present strong temporal correlations. In the literature, there are studies that promise to generate good PRNGs from continuous-time chaotic systems that can be used in cryptography and secure communications. In general, they fail in three main aspects: speed, randomness, and generality. These qualities are essential for any PRNG. A slow PRNG is useless for almost all applications, for example, for real-time encryption. To ensure security, the sequences its chaotic behavior. Therefore, choosing the appropriate step ∆t is not a trivial task. We propose a point of view where the continuous system becomes a discrete map, and the time step used is seen as an extra parameter of this new map. This enables us to characterize the system in terms of ∆t and use statistical quantifiers as well as nonlinear tools to describe the dynamics' evolution of the maps.
Our goal is to propose an extremely simple modification applied to the digitalized continuous-time chaotic system that keeps it oscillating and, at the same time, breaks the internal structures and the time correlation of the outputs, which allows us to apply the discarding method, but discarding a minimum number of bits, so as not to lose speed. Using the standard Marsaglia's Diehard and NIST tests, we show that the resulting map can generate high-quality random numbers to ensure security. Finally, our method is general as it can be applied to any continuous-time chaotic system. We also present the resources needed by a hardware implementation in an FPGA board of the proposed PRNG and compare them with the ones of the original map, showing that the circuit complexity remains almost the same.
The rest of the paper is organized as follows. Section 2 presents the ordinary differential equations that concern us, briefly describes the methods that we consider for the time discretization of those systems and presents the proposed modification over the Euler method. Section 3 gives a short review of the quantifiers employed to characterize the maps' chaoticity and randomness. In Section 4, the obtained results when applying the proposed methodology to the Rössler system are presented. There, we develop new maps that emerge from applying the numerical methods and the proposed modification. Finally, in Section 5, we draw some concluding remarks.

Continuous-Time Chaotic Systems
The general form of a typical continuous-time chaotic system is as follows: where f (t, u) are nonlinear functions of time and the states variables u. Given an initial value u(t 0 ) = u 0 , these systems have a determined evolution. However, there is no general formula to solve this kind of equation; in fact, most of these first-order differential equations cannot be analytically solved [23]. That is why particular time-dependent solutions are most often sought with numerical means. Thus, multiple techniques that approximate the output of the system have emerged. Among the ample range of possibilities for making such a job, the choice of one of them depends on several factors. When an exact reproduction of the continuous system dynamics is required, powerful numerical methods that involve pre-iterations and variable time-steps are mandatory. However, when using these systems, such as PRNGs, the criteria change for choosing which method to use changes. It switches to the ones that allow the output to meet the required properties along with the strong restrictions regarding the hardware implementation, i.e., minimize the required resources and latency, maximize throughput and operation frequency. Considering the above, we focus on fixed integration step methods. In this context, we evaluate the following three well-known numerical methods for the time-digitization of continuous-time systems.

Fourth Order Runge-Kutta Method (RK4)
The main idea of this method is the precalculation of stages at various points using samples of f to obtain the next step [24].
where u t is the discrete-time state variable and ∆t is the time step size.

Heun's Method (HUN)
Heun's method considers the tangent lines to the solution curve at both ends of the interval. This method requires two stages of calculation as follows: Among all numerical procedures for solving ordinary differential equations with a given initial value, the simplest one is Euler's method in which differentials are approximated by a trapezoid with base ∆t, as Equation (4) shows. Euler's method is a one-step algorithm; that is, in order to calculate the variables at the time t + ∆t, it is only necessary to know the values at the previous instant. where:

Modified Euler Proposed Method (EUR_MOD)
Choosing a large ∆t would be believed to help de-correlate the output of the system and thereby improve its randomness. However, the largest ∆t before the system loses its chaotic behavior is not big enough to break its temporal structures. There are numerous proposals to increase the randomness of chaotic systems [25][26][27]. In some of them, postprocessing the outputs is proposed; however, this idea adds hardware and increases latency. Other works propose to disturb the system with external noise, to switch between one or more chaotic maps [28]. Then, complexity is added to the resulting circuit, but the achieved randomness of the final system is just that of the disturbance. As said, the objective is to destroy the temporal correlations while keeping the chaotic oscillation of the system. Furthermore, it is desired to minimize the hardware resources and increase throughput and speed. Our idea is to combine the time-digitalization process with the randomization one. Thus, we have selected Euler's method as is the simplest one and thereby will require the least amount of hardware to be implemented. Based on Equation (5), with the idea of breaking the temporal structure, we apply the following modification: where u t mod 2 returns the remainder of a division after u t is divided by 2. It returns 1 if u t is odd, or 0 if it is even. The parameters p 1 , p 2 and p 3 ∈ [0, 1], so for the particular case where p 1 = p 2 = p 3 = 0, the map converges to ROS EUR map.
The modification consists of incorporating one extra term into each function. This term is simply another state variable multiplied by 0.5. That term will be added or subtracted from the function depending on the parity of the current state variable.

Quantifiers
The time-digitization of continuous systems that turns them into maps generates changes in their dynamics. Chaoticity, stochasticity, and mixing properties change, so the following tools are used to analyze them.

Maximum Lyapunov Exponent
A chaotic orbit (chaotic attractor) is aperiodic, meaning that it never repeats exactly itself, and the oscillation persists for a time tending to infinity. The attractor's movements exhibit sensitive dependence on the initial conditions. This means that two trajectories that start very close, quickly diverge; thus, they will have very different futures. The practical implication of this is that long-term prediction becomes impossible, as small uncertainties are rapidly amplified. The separation δ(t) between two trajectories of the same system that initially differ δ 0 evolves exponentially in the way of Equation (7): Therefore, neighboring trajectories separate exponentially fast. The number λ is called the Lyapunov exponent. When this exponent is positive, it is said that the system has a time horizon beyond which the prediction fails at tolerance a. Actually, λ depends on the trajectory that is being studied. Therefore, it must be averaged over many points of the same trajectory to estimate its true value. In addition, each system has as many Lyapunov exponents as dimensions. The largest of them, known as the maximum Lyapunov exponent (MLE), is of special significance since a positive value indicates the possible existence of chaos [29,30]. Nevertheless, this is a necessary but not sufficient condition of chaoticity since a divergent system can have positive MLE. Therefore, for a system to be chaotic, in addition to having some positive Lyapunov exponent, it must have a bounded nondivergent trajectory in the phase plane.

Bifurcation Diagram
A bifurcation diagram allows studying the changes in the qualitative or topological structure of the trajectories of a dynamical system. It shows the visited values of a system as a function of a certain parameter. It allows differentiating areas of the parameter in which the system behaves like fixed points, periodic orbits, or chaotic attractors [31]. We can say that bifurcation occurs in a dynamical system when a small smooth change of a parameter causes a sudden 'qualitative' or topological change in the dynamical system's behavior.

Probability Density Function (PDF)
The randomness quantifiers used here are functional of the PDF P associated with the data sequence under analysis. The determination of a PDF can be done using several different methods [32], and their applicability depends on particular characteristics of the data, such as stationarity, time series length, parameter variation, and level of noise contamination. The PDF and the sample space are inextricably linked so it is a nontrivial problem to obtain the optimal PDF to extract the desired information. Here, we have employed the Bandt and Pompe approach, and this PDF is able to satisfactorily show the temporal correlations of higher orders [33,34]. The delay method has been used to extract time causal information from the sequences. A delayed reconstruction in D dimensions is formed by the vectors x n given as: The lag or delay time v is the time difference in the number of samples (or in time units τ = v∆t) between adjacent components of the delay vectors. A good estimate of the lag time is very difficult to obtain. If τ is small compared to the internal time scales of the system, successive elements of the delay vectors are strongly correlated, whereas if τ is very large, successive elements are already almost independent. Among the existing proposals, we have adopted the first zero of the autocorrelation function of the signal as the τ value [30]. This algorithm to extract the Bandt-Pompe PDF has been widely addressed and described by previous works [35].

Normalized Shannon Entropy
The well-known normalized Shannon entropy denotes the amount of "disorder" a system presents. It has been shown to be able to successfully characterize determinism and stochasticity of generated sequences [32]. This information theory quantifier is a functional of the probability density function and is defined by the normalized Shannon expression (Equation (9)): where N is the number of elements of the alphabet. We denote permutation entropy (H BP ) as the result of applying the normalized entropy to the PDF proposed by Band and Pompe, which quantifies the causality of the symbolic series discarding amplitude information.

Statistical Complexity Measure
A statistical complexity measure, denoted by C, is a general indicator of structure or correlation. This measure vanishes in the extreme ordered and disordered limits ("boundary conditions"). During the last decade and a half, different measures of statistical complexity have been proposed [36]. Here, we have adopted the functional form introduced by López Ruiz et al. [37] with the modifications advanced by Lamberti et al. [38], given by Equation (10), [39].
where P e is the uniform distribution, and Q j is the so-called "disequilibrium", defined in terms of the extensive Jensen-Shannon divergence, which in turn induces a squared metric, [39] (Equation (11)).
where Q 0 is the normalization constant, Equation (12), and is obtained when the system is deterministic; that is, only one component of P is equal to one, and the remaining components are equal to zero: This quantifier detects internal structures from the symbol source when it is applied to the Bandt and Pompe PDF; thus, we denoted permutation complexity C BP as the resulting quantity.
The juxtaposition in a two-dimensional graph of the quantifiers H BP and C BP has demonstrated to be particularly efficient to reveal properties of the underlying processes from some measurable or observable quantity, called causal Eentropy × complexity plane [40]. High values of C BP correspond to time series with immersed structures, which occurs with chaotic series. On the other hand, the point C BP = 0 and H BP = 1 are that of a sequence with no internal correlations. There are many relevant applications of the H BP × C BP plane; for example, in [34], Rosso et al. use this plane to discriminate between stochastic and chaotic series, in [41], the authors employ it as a tool for distinguishing songs, and Zunino and Ribeiro utilize it to discriminate image textures [42], just to mention a few.

Statistical Randomness Tests
For a sequence to be suitable to be used as PRNG, it is necessary to successfully pass statistical tests. Here, we employed NIST Statistical Test Suite and Marsaglia Diehard tests.

NIST Statistical Test Suite
The NIST SP 800-22 test suite [43] consists of 15 statistical randomness tests that are applied to binary data stream files. It requires the size of each sequence length to be of the order 10 3 to 10 7 . For each test, it yields p-values, and it also checks the proportion of passing sequences and the uniform distribution of the p-values.

Marsaglia Diehard Tests
The 15 statistical tests that make up the Diehard battery should be applied independently over files of several million 32-bit integers. Their output is a statistical p-value. To evidence randomness, each test output should be uniformly distributed between 0 and 1.
The tests should be repeated multiple times with different integer sets to demonstrate the robustness of outcomes.

Results
To show the proposed method, the well-known Rössler system is used here. This continuous-time chaotic system is defined by the following set of coupled ordinary differential Equations [44]: Applying the digitalizing methods mentioned in Section 2, the following maps, which include ∆t as a new parameter, are obtained: We have employed parameters a = 0.2, b = 0.2, and c = 5.7 that assure chaotic behavior of the continuos-time Rössler system. In the case of the ROS EUR_MOD map, the parameters p 1 = p 2 = p 3 = 0.4 were used unless specified otherwise. Since our objective is to utilize the systems as PRNGs, based on subsequent experiments, we have followed three main steps:

1.
First, we analyzed the chaotic behavior when the systems are digitalized in time; focusing on the impact on the dynamic of each discretization method and its dependence on ∆t (Section 4.1). Therefore, we calculate the MLE [29] and bifurcation diagrams of the emerged maps. Note that at this point, we do not consider amplitude discretization of the systems. Therefore, we employ a floating-point arithmetic (IEEE 754 double-precision standard) for the calculations.

2.
The second step deals with the amplitude digitization effect (Section 4.2). Then, we analyze the statistical properties, focusing on achieving the highest randomness.

3.
Finally, we present the hardware implementation of the obtained PRNG that is based on the proposed modification to the system digitalized in time by Euler's method and iterated using signed fixed-point architecture. We also show the resources needed to implement it in an FPGA board (Section 4.3).

Time Digitization Analysis
In all cases, the quantifiers are averaged over 100 surrogates starting at different initial conditions, a transitory of 8 × 10 6 is first deleted, and the maps are then iterated 10 6 times. The minimum ∆t is iterated 10 6 times, and for higher ∆t, the iterations are decreased so as to cover the same attractor window time. In order to understand the maps' behaviors, Figure 1 shows the 3D phase space for some ∆t values of the ROS HUN map. There, it can be seen how attractors change and evolve. It is clearly shown that even though the continuous-time system attractor is blurred by the increase of ∆t, new attractors appear, and these attractors may be even more chaotic than those of the continuous-time systems (depicted by their MLE value). As expected, smaller values of ∆t reproduce an attractor closer to that of the continuous-time system; however, in many cases, they converge to short cycles or fixed points, losing their chaotic behavior.

Topological Analysis
We have used Sprott's method to calculate the MLE of the maps using ∆t as a parameter, [45]. Figure 2a shows how the MLE varies with ∆t in the case of the time-digitalizing Rössler system using the three methods. Each resulting map presents a different behavior regarding the existence of chaoticity with ∆t as a parameter. As it may be supposed, the ROS RK4 map (red line) seems to preserve the chaotic behavior for larger values of ∆t, while the ROS EUR and ROS HUN maps (yellow and blue lines, respectively) behave similarly. In the case of the ROS HUN map, it presents an isolated range of ∆t between ∼0.17 and ∼0.19, where the system behaves chaotically, and it also presents isolated low values of MLE for some time steps indicating low or no chaoticity. The ROS EUR map is the first one that losses its chaotic behavior for higher time steps. It is interesting to note that even though all three maps are derived from the same system for small values of ∆t, the MLE does not show similar values. The reason may be the accumulated arithmetic errors that prevent following the continuous-time attractors. In addition, it can be seen that there are some cases where larger time steps present higher values of MLE.
In Figure 2b, it can be seen that the proposed modification increases the chaotic behavior of the system. The ROS EUR_MOD map presents higher values of MLE than the ROS EUR map. To show the generality of the proposed method, Figure 3 shows the 3D phase of Rössler and Lorenz systems digitalized by Euler's method (ROS EUR and LOR EUR ) in black, and their modified maps (ROS EUR_MOD and LOR EUR_MOD ) in gray. It can be seen how the proposed modification breaks the inner structures and temporal correlations of the sequences and keeps the chaotic behavior. This enables the retainment of more bits for the PRNG output and, in this way, increases the throughput.  Regarding the bifurcation diagram, we have built the diagrams using Poincaré maps, which is the intersection with a certain surface. Then, the bifurcation diagrams show all the visited values by the systems [31]. Figure 4 shows the bifurcation diagram of the ROS HUN map superimposed with the MLE (red line). It can be seen how the MLE is able to effectively predict the chaoticity of the map. Within the chaotic region, some isolated gaps that correspond to low chaoticity can be seen. These gaps match with low values of the MLE. It can be seen that from ∆t ∼ 0.105, it completely loses its chaoticity and it stays on a periodic cycle until ∆t ∼ 0.167. The darker areas of the chaotic region imply that the system, while being in the state of chaos, spends more time there than in the lightly shaded regions. The most interesting places inside that region are the "white spaces", which have an important role in the transition to chaos. The "white regions" and their boundaries also show the instability of the initial conditions, another important aspect of the chaos.

Statistical Analysis
Up to this point, we have only analyzed the chaoticity of the maps; however, we do not have information about the randomness that their outputs present. The applications in which these maps are intended to be used require that their sequences, in addition to being chaotic, have no internal structures and all their possible outputs appear in a balanced way. To evaluate this, we calculate the randomness quantifiers described in Section 3. Each quantifier has been averaged over 100 files. Every surrogate file starts with a different initial condition, a transitory of 8 × 10 6 iterations was first deleted, and the maps were then iterated 10 6 times.
To extract causal informati on by calculating H BP from these observations, we employ here x state variable, v = 1, and D = 6. Figure 5 shows the plane H BP × C BP for the Rössler system time-digitalized with the three mentioned methods using different values of ∆t. The continuous curves correspond to the boundaries of values for the statistical complexity, as a function of the value of the normalized Shannon entropy [38]. It can be seen that in the cases where the system is unmodified (ROS RK4 , ROS HUE , and ROS EUR maps), the quantifiers remain in the same area, that is, strong correlations and poor balance of values. When the proposed modification is applied, the quantifiers move towards the area of chaotic maps. The output of the system slightly improves the balance of its values and also increases its inner correlations. Figure 5. Causal entropy × complexity plane, Rössler system using the three methods of timediscretization for 0.0001 ≤ ∆t to the higher ∆t that could be reached before the maps diverge. Red points are the ROS RK4 map, blue points are the ROS HUE map, and yellow points are the ROS EU map. The ROS EUR and ROS EUR_MOD maps using S (41,38) are the green points and black stars, respectively. Finally, our proposed PRNG (ROS EUR_MOD map S (41,38) considering the 38 least significant bits) are the pink points, and these are the best sequences in terms of randomness (closest to the ideal point H BP = 1, C BP = 0).

Amplitude Digitization Analysis
We iterate the maps using signed fixed-point architecture for analyzing the effect of amplitude digitalization [46]. As said, our goal is to develop a hardware-implemented PRNG, which is why we have selected fixed-point architecture and Euler's discretization method because of the simplicity they mean in terms of hardware design (Equations (14) and (15)).
• ROS EUR map: • ROS EUR_MOD map: We have employed words of wl bits, with fl bits to represent the fractional part in two's complement arithmetic; this architecture is represented by S(wl,fl). Equation (16) outlines the data format used for each state variable.  Figure 6a,b show the 3D phase space of the ROS EUR and ROS EUR_MOD maps, iterated using fixed-point architecture. There, we used ∆t = 0.001. Figure 6a shows the phase space of ROS EUR_MOD map using S (41,38). It can be seen that the phase space does not present significant changes compared to that of the attractor iterated with floating-point arithmetic (Figure 3b). Therefore, the outputs of the proposed PRNG are the k LSB least significant bits of each state variable u t (see Equation (17)). This is a commonly used procedure in PRNGs [10,47].
It is desired to find the minimum wl that keeps oscillating the attractor in a nonperiodic way, and the largest k LSB that produces the output sequences to pass the Marsaglia and NIST tests. In Figure 6b, it can be seen how the obtained sequence does not present any structure and all the space is equally filled.

Statistical Analysis
Returning to Figure 5 where the causal entropy × complexity plane was shown, the red stars correspond to the ROS EUR_MOD map using a signed fixed-point architecture with 41 bits, of which 38 are used to represent the fractional part (S (41,38)). It can be seen that the utilization of fixed-point arithmetic does not influence the statistical properties of the proposed system. Finally, when the most significant bits are discarded (pink point), it can be seen that the sequences reach the ideal point in terms of randomness, H BP = 1 and C BP = 0. Table 1 shows the results obtained when applying the Marsaglia battery of statistical tests to the original system (ROS EUR map) and the modified one (ROS EUR_MOD map). It can be seen that the ROS EUR map needs more bits to pass the tests. There, it is demonstrated that the proposed method keeps the system oscillating and enables it to discard fewer bits of each output. This is due to the fact that the time correlations and internal structures are destroyed. To confirm our proposal's usefulness, we keep k LSB bits of the output of both the original and the modified maps, iterated using S(wl,fl) arithmetic, and test them using Diehard tests. Table 1 shows that the sequences generated by ROS EUR_MOD pass Marsaglia tests using fewer word bits (wl = 41) and allows to keep more bits per iteration for the PRNG (higher k LSB ) than ROS EUR map (wl = 56). Table 2 shows the results of testing the proposed PRNG via the NIST SP 800-22 test suite. In agreement with the values used in the literature [13,[47][48][49], 1000 sequences of length 10 6 bits each have been tested. For a significance level of 0.01 (α = 0.01) and 1000 samples, the minimum pass rate for each statistical test is approximately 980, with the exception of the random excursion (variant) test where it is approximately 597 for a sample size of 611 binary sequences. The proposed PRNG passes all these tests.

Hardware Implementation
Here we show that the proposed modification is extremely simple to implement and assures the required randomness. Figure 7 shows a schematic of a hardware implementation of ROS EUR and ROS EUR_MOD maps. In the latter map, the parameters used were p 1 = p 2 = p 3 = 0.5. Figure 7a shows a schematic of the recursive function for x of a general map obtained by the Euler method applied to a continuous-time chaotic system (recursive functions for y and z are analog, and for simplicity are not shown). There, the f x block receives the three state variables at time t and calculates the next output at time t + 1. In the case of the Rössler system studied here, this term is f x = −y t − z t . Its output is multiplied by ∆t and added to x t in order to generate x t+∆t . This value is then latched by a register at each clock cycle. Figure 7b shows the proposed modified circuit. It can be seen that it consists of just one extra term. This term makes the product of one state variable (in this case z t ) by 0.5. However, for its implementation, no multiplier is required since this term is a right-shifted version of the state variable. Then the least significant bit of the integer part of x t is fed back to select if the new term is added or subtracted. The whole term can be implemented by a positive or negative right-shifted version of z t , as it can be seen in the light blue square of Figure 7b. The figure shows that it requires very few resources. A comparison of the implementation results of the proposed PRNG with three other PRNGs can be seen in Table 3. ROS EUR refers to a PRNG based on the ROS EUR map, which basically consists of implementing the map and then applying the discarding method. The map requires a word length of 56 bits to enable the extraction of 38 bits (of each variable) on each iteration. This word length is the minimum number of bits for ROS EUR to pass the Marsaglia test (see Table 1). We can see that this PRNG requires more resources and that the maximum frequency and throughput are lower than that of the proposed PRNG, which is due to the need to use a larger word size to ensure the randomness of the output. In the third column, the resources used by an implementation of the well-known PRNG Mersenne Twister (MT19937) implemented in an FPGA are shown [50]. We can see that the resource requirement is slightly higher than that of the proposed PRNG. The maximum operating frequency is lower and the achieved throughput is lower as well. The last column corresponds to a continuos-time, chaotic-based PRNG that employs a linear feedback shift register to obtain the transition between Lorenz-like and Chen-like behaviors. Then, authors keep the eight least significant bits and xor them (the discarding method) [13]. Although this generator does not comply with being generic (it only works for the Lü-like chaotic system because it is capable of exhibiting both Lorenz-like and Chen-like chaotic system behaviors for different parameter values), we include it in order to compare its performance.

Conclusions
In this paper, we showed that the digitalization method and the time step have a significant influence on digitalized systems' dynamics and, therefore, on the sequences generated by them. The chaotic behavior and statistical degree of the sequences were analyzed using tools from nonlinear systems analysis and information theory quantifiers. In that way and with the objective of using these systems as PRNGs, we proposed a modification to the map generated by Euler's method that destroys the time correlations of the output and keeps the chaotic oscillation. We have also analyzed the randomness behavior with the amplitude discretization using different precision and data widths. The H BP × C BP plane shows that the three methods of digitalization analyzed produce outputs in almost the same area, poor balance of amplitudes, and strong inner correlations. The proposed method produces both floating and fixed-point architectures moving towards the typical area of the chaotic maps.
Our goal was to demonstrate that our proposed modification to the digitalized Euler system generates the most random output, which is located in the optimum point of H BP × C BP plane (uncorrelated noise). The proposed modification achieves the lower value of wl and higher value of k LSB when passing the Marsaglia Diehard and NIST tests. Our method is general and can be applied to any continuous-time chaotic system.
Regarding portability and reproducibility, which are important PRNG properties, the proposed schematic defines the architecture and precision of the variables and the internal calculations since it is a hardware implementation. Then, identical results will be obtained in any programmable device, and the repeatability of the results will be ensured.
Further study on the basin of attraction of the digitized system would be required to define the set of available seeds for the PRNG.
Finally, we compared the resources needed to implement our and other existing methods to obtain PRNGs. We showed that the proposed PRNG is superior in terms of resources, maximum frequency of operation, and throughput. Data Availability Statement: The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.