Generating Even More Chaotic Instances in Hardware

Abstract

It is well known that multiplication inside a computer does not follow the associative property because of roundoff effects. It is possible to use this fact to generate other different chaotic instances of chaotic maps or oscillators when a multiplication of three terms appears. Chaos is very sensitive to small changes in the initial conditions and amplifies these small rounding effects. We use this condition to build different chaotic instances, which give different results, of the Lü oscillator and the 2D map, and we show one application to create new instances of a pseudo random number generator using the 2D map. Both chaotic systems are simulated in software and in hardware within an FPGA where another 144 different 2D map instances and 81 different Lü oscillators can be created. To best of our knowledge, it is the first paper that analyze the construction of new chaotic entities by using the roundoff effects.

1. Introduction

A chaotic map [1,2] is a discrete system which generally has a very simple mathematical description. A chaotic oscillator [3] is a continuous time system which is described by three or more differential equations. Those two systems must have a large positive Lyapunov exponent and high entropy to produce chaos [2].

The fixed and floating point arithmetics can be performed in current microprocessors and microcontrollers. It will be shown that the multiplication in both arithmetic systems does not follow the associative property, this is $\left(ab\right)c\ne a\left(bc\right)\ne \left(ac\right)b$ for different numbers a, b, and c. These details must be known for any person that studies Computer Science, and in general these details are given as understood in any application. This characteristic that multiplication does not follow the associative property will be used to build different random number generators from a same chaotic map, and also to build different chaotic oscillators.

In [4,5] were already observed that in the implementation of chaotic maps, the order of the multiplications must be specified if one requires to reproduce the results shown in the literature. Authors in [6] take advantage of these roundoff errors to propose a simple method to determine a lower bound error and the critical time for the pseudo-orbits in chaotic oscillators, and also its connection to the value of the maximum Lyapunov exponent. In [7] the degradation produced by the fixed point implementation of chaotic maps is studied and compared with their floating point implementation.

The motivation of this work was to observe that the implementation in hardware of a chaotic entity must follow carefully the simulation in software, or vice versa, if one requires to obtain the same results. But also we observe that a modified implementation follows the same dynamical behavior. Then, we try in this work to analyze how new entities, of chaotic oscillators or chaotic maps, new in the way that produce different results, can be constructed by modifying the order of the multiplications in their fixed point implementations.

2. Multiplication Inside a Computer Is Not Associative

In algebra, three numbers as $abc$ can be multiplied as $\left(ab\right)c$ or $a\left(bc\right)$ or $\left(ac\right)b$ given all the same result, where the three numbers a, b, and c can be integers of real numbers. This is the associative property which is fulfilled for these numbers. Because of the associativity, grouping parenthesis can be omitted without ambiguity. Inside a computer $\left(ab\right)c\ne a\left(bc\right)\ne \left(ac\right)b$. We are going to explain why the multiplication does not follow the associative property, both in floating point arithmetic or even in fixed point arithmetic.

A floating point (FP) number is composed as $s·2^e$, where s is the significand and e is the exponent. Both s and e are integers numbers. If p bits are used for the significand, it could take values from 0 to $2^p-1$. A FP number is normalized if it has one single number different than zero in the integer part, for example, the number 1234 in decimal is normalized as $1.234\times {10}^3$, and here four decimal numbers are used in the significand. The multiplication of two FP numbers $s_1·2^{e_1}$, and $s_2·2^{e_2}$ is performed multiplying the significands as $s=s_1·s_2$ and the exponents are summed ($e=e_1+e_2$), and also both results are rounded and the final result is normalized. Much more details are involved in FP arithmetic [8]. The multiplication of the two significands is the source than numbers calculated inside the computer are not associative. The multiplication result must be maintained in the same size as the significands, thus that result must be truncated or rounded to the nearest number. We will show this effect with one example. Given these three FP numbers: $1.0789\times 10$, $1.9999\times {10}^{-1}$, and $5.1110\times {10}^0$, multiplying the significands of the first two numbers gives:

  (1.0789 $\times 10$) × (1.9999 $\times {10}^{-1}$) = 2.15769211 $\times {10}^0$ $=2.1577\times {10}^0$, and multiplying this result by the third number results.

(2.1577 $\times {10}^0$) × (5.1110 $\times {10}^0$) = 11.02800470 $\times {10}^0$ $=\mathbf{1}.\mathbf{1028}\times \mathbf{10}$. And now, multiplying the second and the third number results.

(1.9999 $\times {10}^{-1}$) × (5.1110 $\times {10}^0$) = 10.22148890 $\times {10}^{-1}$ $=1.0221\times {10}^0$, and multiplying this last result by the first number gives:

(1.0221 $\times {10}^0$) × (1.0789 $\times 10$) = 1.10274369 $\times 10$ $=\mathbf{1}.\mathbf{1027}\times \mathbf{10}$.  

Both results $11.028$ and $11.027$ differ in 0.001. This small difference will always occur, regardless of whether the multiplication result is truncated or rounded to the nearest number.

Now we are going to analyze the fixed point arithmetic. The notation $A(a.b)$ will be used here to represent a set of numbers that uses a bits in the integer part, and b bits in the fractional part. Each number is of size $a+b+1$ bits (plus the sign bit) and indeed this number is stored as a signed integer of size $a+b+1$ bits.

For a number $x\in A(a.b)$, if the negative numbers of size $a+b+1$ are represented in complement 2 notation [9], the range of numbers that can be represented is

$$-2^a\le x\le 2^a-2^{-b}.$$

(1)

Summing up two numbers $A(a.b)$ results in a number $A\left(\right(a+1).b)$ [9]. The multiplication of two numbers $A(a.b)$ results in a number $A\left(\right(2a+1).2b)$ [9]. It is possible to verify these results by applying the respective operation to the two extreme numbers in Equation (1).

A microprocessor offers the sum and multiplication of two integer numbers and the result is stored in a number of the same size that the operands. In a hardware design for a given application, one must use a big enough number to store the sum of two $A(a.b)$ numbers, and the result to multiply two $a.b$ numbers must be returned to a $A(a.b)$ number. The easiest way to perform this is by truncating the result: the resulted $A(2a.2b)$ number is shifted b bits to the right, and again the number of bits a must be big enough to store the resulted number.

The binary point in the representation $A(a.b)$ is virtual. With three bits for $a+b$ it is possible to represent numbers $A\left(2.1\right)$, $A\left(1.2\right)$ or $A\left(0.3\right)$, which need an adder of 4 bits and a multiplier of two numbers of 4 bits.

The small differences among the multiplications of three terms can be exploited by chaotic systems to generate different behaviors. The chaotic systems of a map and an oscillator will be explained in the next section.

3. Two Chaotic Systems: Maps and Oscillators

The 2D map was proposed by Sprott in 1993 [10] and is defined by

$$\begin{array}{cc}\hfill x_{n+1}& =a_1+a_2{x}_n+a_3{x}_n^2+a_4{x}_n{y}_n+a_5{y}_n+a_6{y}_n^2,\hfill \\ \hfill y_{n+1}& =a_7+a_8{x}_n+a_9{x}_n^2+a_{10}{x}_n{y}_n+a_{11}{y}_n+a_{12}{y}_n^2.\hfill \end{array}$$

(2)

It is possible to observe in Equation (2) that this simple 2D chaotic map uses only sums and multiplications. In [10] the values of the coefficients $a_i$, for $i=\{1,2\cdots 12\}$ in Equation (2) that produce chaos were searched by a computer program. Twelve different sets of values for these coefficients are given in [10]. In this work we use the values $\{-0.6$, $-0.1$, $1.1$, $0.2$, $-0.8$, $0.6$, $-0.7$, $0.7$, $0.7$, $0.3$, $0.6$, $0.9\}$.

There are six quadratic terms in Equation (2), three terms in the expression for $x_{n+1}$: $a_3{x}_n^2$, $a_4{x}_n{y}_n$, and  $a_6{y}_n^2$; and another three quadratics terms for $y_{n+1}$: $a_9{x}_n^2$, $a_{10}{x}_n{y}_n$, and $a_{12}{y}_n^2$. The terms with $x_n^2$ and $y_n^2$ can be calculated in two forms, i.e., $\left(a_3{x}_n\right)x_n$ and $a_3\left(x_n^2\right)$. The other terms can be arranged in three forms, i.e., $a_4{x}_n{y}_n$ can be calculated as $a_4\left(x_n{y}_n\right)$, or  $\left(a_4{x}_n\right)y_n$, or  $\left(a_4{y}_n\right)x_n$. In this manner there exist $2·3·2=12$ different forms to associate the terms to calculate the products for $x_{n+1}$, and the same number of forms to calculate $y_{n+1}$. Therefore, there exist $12·12=144$ different forms to calculate the 2D map. The small differences in the multiplications results are amplified by the chaos, indeed this is a definition for chaos [1,3,11]. 144 different pseudo random number generators as in [11] can be implemented using this 2D map. This result will be shown in Section 4.

The Lü oscillator Chapters 4 and 6 in [9] is defined as

$$\begin{array}{cc}\hfill \dot{x}& =y,\hfill \\ \hfill \dot{y}& =z,\hfill \\ \hfill \dot{z}& =-ax-by-cz+df\left(x\right),\hfill \end{array}$$

(3)

where a, b, c, d are fixed coefficient values, and f function will be defined later. If Equation (3) is integrated using the Euler method, it results in

$$\begin{array}{cc}\hfill x_{n+1}& =x_n+h{y}_n,\hfill \\ \hfill y_{n+1}& =y_n+h{z}_n,\hfill \\ \hfill z_{n+1}& =z_n-ha{x}_n-hb{y}_n-hc{z}_n+hdf\left(x_n\right),\hfill \end{array}$$

(4)

where h is the integration step value. Expanding the parameter h in Equation (4), four terms with two multiplication each one (and different associations on them) can be obtained. An initial state $(x_0,y_0,z_0)$ must be given to apply Equation (4).

Observe in Equation (4) that for the calculation of $z_{n+1}$ there are four terms with two multiplications each, thus $3^4=81$ different oscillators can be built by changing the association in which these multiplications are calculated.

Function f in Equations (3) and (4) is a piece-wise linear function. For an oscillator with p even number of scrolls, $p\ge 2$, the form of f has p plateaus, and $p-1$ transitions between each pair of plateaus. If the width between each pair of transitions is w, the transition has a slope m, and the length of each transition is $2\alpha$, function f has the form in Equation (5).

$$f\left(x\right)=\left\{\begin{array}{cc}m(x-sw)+2sk,\hfill & \mathrm{if}|x-sw|<\alpha ,\hfill \\ (2s+1)k,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(5)

where s is calculated as $s=\lfloor x/w\rfloor$ and $s=\lfloor (p-1)/2\rfloor$ if $s>\lfloor (p-1)/2\rfloor$. In Equation (5), k is a scale factor for f in each plateau. Function f is even, thus when x is negative, it is calculated as $-f(-x)$.

It is interesting to say here that different oscillators can be also obtained in double precision and single precision floating point numbers. But, because these floating point numbers have a much more higher resolution than the used $A\left(4.20\right)$ numbers for the fixed point arithmetic, the difference in the oscillations are obtained at much more simulation time than with fixed point. Figure 1 shows the graphs of two oscillators implemented with double precision FP, single precision FP, and fixed point arithmetic with $A\left(4.20\right)$ numbers, when these two oscillators start to diverge; oscillator 1 is obtained with the term $h{a}_2{x}_n$ (see Equation (4)) calculated as $\left(h{a}_2\right)x_n$, and the oscillator 2 is obtained with $h\left(a_2{x}_n\right)$. The other terms are multiplied following the usual left-to-right precedence for the multiplication.

The four bits size for the integer part in the $A\left(4.20\right)$ numbers were obtained from the simulation with real numbers. The value of 20 for the number of bits for the fractional part was fixed arbitrarily when the oscillator started to work. Perhaps in a given application, this number of bits in the fractional part could be justified or optimized according to the needs of such application.

To measure where the oscillators implemented with fixed point arithmetic start to diverge, the following experiment was performed. Thirty oscillators calculated with code 0000 and other thirty with the code 1000; $y_0$ is fixed equal to 0.1, $x_0$ is randomly obtained within the interval $[0,1]$. Figure 2 shows the mean and ± the standard deviation of the correlation calculated using a window of size 21 samples between the signals calculated for each pair of oscillators. Observing this figure, the signal is uncorrelated from sample 575.

4. Pseudo Random Number Generator

The 2D map described in Section 3 is used in [11] to build a Pseudo Random Number Generator (PRNG) as follows. Fixed point arithmetic with numbers $A\left(3.60\right)$ is used. All coefficient in Equation (2) are lesser than 1, and there are six sumatory terms in Equation (2), thus at most the results could be equal to 6, and three bits are used for the integer part. The number of bits in the fractional part is more arbitrary, with 60 bits can be used numbers of 64 bits in total (A long in C with a 64 bits architecture).

A sequence of 16 bits at each time $n+1$ can be generated as

$$s_{n+1}=\{x_{n+1}\mathrm{mod}256,y_{n+1}\mathrm{mod}256\},$$

(6)

where the initial values $x_0$ and $y_0$ must be given within the interval $[-0.5,0.5]$. With this interval and 60 bits in the fractional part $2^{60}$ different initial values could be given for $x_0$ or $y_0$, thus a key of 120 bits can be used for this PRNG configuration. The size of this key can be increased by increasing the number of bits in the fractional part.

We write a C program to implement the PRNG. The associativity of the terms in the 2D map is coded with the string $abaaba$ where $a\in \{0,1\}$ and $b\in \{0,1,2\}$. This code is passed as a argument to the C program (it is publicly available at https://cs.cinvestav.mx/~fraga/Instances.tar.gz (accessed on 7 November 2022)). The letter a code the terms in Equation (2) which have two forms to realize the multiplications, i.e. the terms with $x_n^2$ and $y_n^2$. The term $a_i{x}_n^3$ be realized as $\left(a_i{x}_n\right)x_n$ or $a_i\left(x_n{x}_n\right)$. The letter b code the multiplications of terms $a_i{x}_n{y}_n$ which can be performed as $\left(a_i{x}_n\right)y_n$, $a_i\left(x_n{y}_n\right)$, or $\left(y_n{a}_i\right)x_n$. The first three letters in the code are to express the associativities to obtain the result of $x_{n+1}$, and the second three letters are used to express the forms to multiply and to obtain the term $y_{n+1}$. Code equal to 0 express the usual left-associativity used in multiplications in the C language.

The first 32 bits generated with four examples of the PRNG are shown in

Table 1. Four output examples of the first 32 bytes generated with the PRNG with the codes shown in the first column.

Code	First 32 Bytes Generated (in Hexadecimal)
000000	3910fd015f5d7e 6213a4eef9f04b74
	a070911d681e2a 7f3a421278a9c41e
100000	3910ed013f55ba 660fa2eff4ed4463
	6de28ed09a4a5ab 5e00900a7d67dff
010000	3910fd015f5c7e 6213a4eef8f04970
	afd1008e89d3a4 782c551274b4b315
020000	3910fd014f5b7e 6312a3eef6ee4769
	8ef1afcf93c1a7 6c1a6e0d76c4990c

. The initial condition used is $x_0=0.1$ and $y_0=0.2$. As we can see in this Table 1, the first 4 bytes (73910fd0) are all equal, after theses 4 bytes the generated bytes start to be different.

The National Institute of Standards and Technology (NIST) of the United States of America created a free software for testing pseudo and random number generators [(P)RNGs]. Specifically, NIST said that the test is for (P)RNGs for cryptographic applications, where randomness is a crucial characteristic [12]. Thus, we test if four generated sequences with the conditions in Table 1 are random by applying the NIST suite of tests [12] to 100 sequences of ${10}^6$ bits. Note that these sequences only vary in the order of the multiplications because the four sequences are generated with the same initial conditions. The four sequences pass all the NIST tests as it it shown in

Table 2. Results of applying NIST tests to the four sequences of Table 1.

	Code	000000		100000		010000		020000
	Test Name	p-Value	Prop.	p-Value	Prop.	p-Value	Prop.	p-Value	Prop.
1	Frequency	0.04281	0.99	0.89776	0.99	0.67869	0.99	0.99425	1.00
2	Block Frequency	0.55442	0.99	0.18156	0.98	0.11539	0.98	0.53415	1.00
3	Cumulative Sums	0.08367	0.99	0.16457	0.99	0.42656	1.00	0.94130	0.99
4	Runs	0.83431	0.98	0.55442	1.00	0.53415	0.99	0.07572	1.00
5	Longest Run	0.75976	0.99	0.81654	0.99	0.94631	0.99	0.65793	0.99
6	Rank	0.59555	0.98	0.69931	0.98	0.63712	1.00	0.43727	0.99
7	FFT	0.92408	0.99	0.00888	0.99	0.91141	0.97	0.43727	1.00
8	NonOverlapping Template	0.47649	0.99	0.486662	0.99	0.50935	0.99	0.48339	0.99
9	Overlapping Template	0.41902	0.98	0.554420	0.99	0.33454	1.00	0.15376	0.97
10	Universal	0.65793	0.98	0.38383	0.99	0.00163	1.00	0.38383	0.99
11	Approximate Entropy	0.08052	0.99	0.03757	0.98	0.00095	0.98	0.15376	1.00
12	Random Excursions	0.19609	0.99	0.80218	0.99	0.43332	0.99	0.70287	0.99
13	Random ExcursionsVariant	0.31675	0.99	0.66026	1.00	0.72320	0.99	0.97041	0.98
14	Linear Complexity	0.35049	0.99	0.30413	1.00	0.35049	1.00	0.36692	0.99
15	Serial	0.68652	0.97	0.84311	1.00	0.13677	1.00	0.16766	1.00

(the p-value must be greater than 0.01, and the proportion value must be greater than 0.96).

To measure where the generated sequences start to be uncorrelated another experiment was performed. With the set of five values $\{-0.5,-0.4,0.0,0.1,0.2\}$, 10 combinations of two values $(x_0,y_0)$ can be calculated; 20 more pairs of values are generated as $(x_0+ϵ,y_0)$ and $(x_0,y_0+ϵ)$ with $ϵ=2^{-60}$. The mean and standard deviation of the 30 correlations between sequences generated by the PRNG with code 000000 and 100000 using the 30 explained initial conditions are shown in Figure 3. A window of size 19 bits was used to calculate the correlations. From this figure, one can observe that sequences are uncorrelated after bit 250, or around 32 bytes. Therefore, the first generated 32 bytes by the PRNG can be discharged to ensure uncorrelated sequences.

5. Hardware Designs

The FPGA implementations of the Lü oscillator and 2D map were carried out from Equations (2) and (4). We can observe that there are three types of multiplications, $ab$, $a{b}^2$, and $abc$. Figure 4 and Figure 5 present the design for two composite multipliers, $abc$ and $a{b}^2$, respectively; note that both multipliers are made up of two $ab$ simple multipliers, represented by a black-filled multiplication symbol. Furthermore, to shorten the critical path, a register was placed at the output of each simple multiplier, this means, the multiplier result will be valid after two clock cycles. The $abc$ multiplier has three multiplexers that allow selecting some of the input combinations $a\left(bc\right)$, or $\left(ab\right)c$, or $b\left(ac\right)$, the control signal for these muxes is two bits. The multiplier $a{b}^2$, has two muxes that are controlled by a one-bit signal and allows to select between the combinations $\left(ab\right)b$ or $a\left(bb\right)$.

In Figure 6, we can observe the block diagram of the 2D map described from Equation (2). The circuit has five 64-bits multipliers: two simple, two $a{b}^2$, and one $abc$. The point-fixed setting is the sign bit, three bits for the integer part, and 60 bits for the fractional part. The code is selected by means of a three-bit signal, where the most significand bit corresponds to a, the control signal for the $a{b}^2$ multiplier; the two least significand bits correspond to b, the control signal for the $abc$ multiplier.

The design of the Lü chaotic oscillator is presented in Figure 7 and Figure 8. The design shown in Figure 8 uses 11 multipliers, from them 7 are simple multipliers and 4 are $abc$ multipliers. The configuration of the simple multiplier is 24 bits, where one bit is used for the sign, three for the integer part and 20 bits in the fractional part. Figure 7 presents the block diagram of the function f in Equation (5), this function is implemented only with simple multipliers. Due to the data dependence, it is necessary to wait for two clock cycles to obtain a result.

Particularly, for Figure 8, $x_{n+1}$ and $y_{n+1}$ are computed with simple multipliers, while $z_{n+1}$ uses four $abc$-type multipliers and other simple arithmetic operations. Each $abc$ multiplier is connected to an independent 2-bit signal to be able to select among the 81 different oscillators. The circuit delivers new output values every five clock cycles.

The 2D map and the Lü oscillator were implemented on the Zybo Z7-20 development board with Zynq-7000 ARM/FPGA SoC.

Table 3. Area utilization for the implemented designs. Mark ‘*’ corresponds to those designs with a selector of the multiplication associations.

Design	LUTs	FFs	DPS
OscLu	464	322	20
OscLu*	706	322	24
2DMap	2629	1149	216
2DMap*	9304	1654	219

shows the area resources utilization for the implementations of two versions, with (OscLU* and 2DMap*) and without the logic for selecting the combination of multipliers (OscLU and 2DMap). There is an increment of LUTs, and this is related to the deployment of the multiplexers needed. The OscLU* area increases about 1.52 times the design OscLu; meanwhile, the 2DMap area increases about 3.53 times the design OscLu. However, if it is considered that now it is possible to obtain 144 different 2D map instances and 81 different Lü oscillators can be created, the use of additions LUTs is justified. The development board and the oscillator signal on the oscilloscope are shown in Figure 9.

Two custom Vivado’s Intellectual Properties (IPs) were created with AXI bus interface, the controller registers functioned to send and receive data to the IP. The mux select signals were declared as external outputs and were directly mapped to the board switches. The write functions send the initial conditions, while the read functions retrieve the results from the IP. The results were sent over the serial port, and it was verified that the software simulation obtains the same results as the multipliers implemented in hardware directly. Finally, for OscLU*, the values of x and y were passed through a DAC in order to observe the real operation of the oscillator as it can be seen in Figure 9. During the experiment, the board switches changed the association of the multiplications and the change of the phase diagram on the oscilloscope can be visualized.

The VHDL code of the hardware designs is publicly available at https://cs.cinvestav.mx/~fraga/Instances.tar.gz (accessed on 7 November 2022).

6. Discussion

In the implementation of a chaotic entity with m equations, $n_{mi}$ number of terms with i variables, $i\ge 3$, in the m equation, and p total number of terms starting with k number of variables, $k\ge 3$, the upper bond in the number of different instances that can be obtained is equal to

$$\prod _{j=1}^m\prod _{i=k}^p{o}_i^{n_{ji}},$$

(7)

where $o_i$ is the number of different associations in the multiplications that can be obtained with one term with i variables. $o_i$ can be obtained recursively as

$$o_3=C(3,2)=3,o_{i+1}=C(i+1,2)o_i,$$

where $C(k,2)$ is the number of combinations in a set of k elements taken 2 elements at a time without repetition.

For example, for the Lü oscillator in Equation (4) there is one equation with 4 terms ($n_{13}=4$) with 3 variables each one, thus we have $\left(1\right)3^4=81$ different implementations of this oscillator.

It could be very interesting to design chaotic entities with terms with more than three variables because the number of instances grows exponentially with the number of such terms. A chaotic entity with two terms with four variables each one could produce ${\left[\left(6\right)\left(3\right)\right]}^2=324$ different instances.

We believe that in applications where several PRNG are needed, our design could be useful. At least it is not necessary to built several PRNGs based in different equations; if the chaotic entity has triplets terms, then several instances of it can be built. Still we are searching for specific applications of the proposed designs.

FPGA has been seen as a intermediate design between the software and a dedicate design on an specific ASIC chip. In a microcontroler or in a Raspberry PI the experiment with the 2D map can be realized using the version in software which is publicly provided. In an architecture of 32 bits, the multiplication of two numbers of 64 bits must be programmed [11].

7. Conclusions

We created new instances, which give different results with the same initial conditions, of both the 2D map and Lü chaotic oscillators. The new instances are generated by different association of the terms when a three-terms multiplication appears. We build the instances both in software and in hardware, within an FPGA, where another 144 different 2D map instances and 81 different Lü oscillators can be created. We show one application to generate different Pseudo Random Number Generators using the 2D map. We show that effectively new instances of chaotic objects can be created in software or in hardware by using different association in their multiplications.