Replication study challenges and new number formats for chaotic pseudo random number

: Chaotic Pseudo Random Number Generators have been seen as a promising candidate for secure random number generation. Using the logistic map as state transition function, we perform number generation experiments that illustrate the challenges when trying to do a replication study. Those challenges range from uncertainties about the rounding mode in arithmetic hardware over chosen number representations for variables to compiler or programmer decisions on evaluation order for arithmetic expressions. We find that different decisions lead to different streams with different security properties, where we focus on period length. However, descriptions in articles often are not detailed enough to deduce all decisions unambiguously. To address similar problems in other replication studies for security applications, we propose recommendations for descriptions of numerical experiments on security applications to avoid the above challenges. Moreover, we use the results to propose the use of higher-radix and mixed-radix representations to trade storage size for period length, and investigate if exploiting the symmetry of the logistic map function for number representation is advantageous.


Introduction
Chaotic pseudo-random number generators (cPRNGs), in particular the logistic map [1], have been considered as a means for secure PRNGs. Experiments [11,13] revealed that the period lengths are not satisfactory in relation to state size. When studying the descriptions of these experiments, we found that they are missing some implementation details that also influence the period length, even for fixed parameter settings. Hence, replicating those studies experiences some challenges, which we report and from which we derive recommendations for describing experiments including numeric calculations. We also underline experimentally that the implementation details such as chosen number representation, rounding mode and order of evaluation in a formula influence the results.
Beyond the typical number representations integer, i. e. fixed-point, and (binary) floating-point, we next present number representation formats in-between those two that use radix higher than 2 and possibly mixed radix. We demonstrate experimentally that those formats enable to trade accuracy and storage. Thus, our investigations also exhibit relations to the field of approximate computing. Moreover, we investigate if exploiting the symmetry of the logistic map function, i. e. f (x) = f (1 − x) for x ∈ [0; 1] can be mirrored in the number representation in order to achieve longer period lengths.
The remainder of this article is structured as follows. Section 2 gives information on pseudo-random number generation for security applications, especially when using chaotic functions, summarizes basic information on number representations for rational numbers, and discusses related work. Section 3 presents the challenges that we encountered in experiments with chaotic pseudorandom number generation based on logistic map and reports on own experiments to demonstrate those challenges. In Section 4 we present number formats between standard floating-point and fixed-point representations, as well as experimental results that indicate that those formats allow a tradeoff between storage size and achievable period lengths. Section 5 investigates how symmetry in the logistic map function can be exploited in number representations, while Section 6 concludes with recommendations on configuration details to be reported about numerical experiments, and an outlook on future work.

Chaotic pseudo-random number generation
A pseudo-random number generator is a finite state automaton. Upon seeding, it receives input which is transformed into a value for the inner state. When the generator is called from an application, an output is generated by applying an output function on the state, and the state transitions to the follow-up state by applying the state transition function (or update function) on the current state. The follow-up state is unique for each state, as the generator receives no further input until a re-seeding takes place. The pseudo-random number generator should produce outputs that follow the claimed distribution (e. g. if bits are output with uniform distribution, each value should appear 50 % of the time), pass statistical tests of randomness like Diehard [17] (also 0, 1, 0, 1, 0, 1 has uniform distribution but does not look random), and should have a sufficiently long period before output sequences repeat. The latter is unavoidable as the inner state is finite. A good exposition of random number generators can be found in [12].
For using pseudo-random numbers in security applications, it should not be possible to infer the inner state from outputs, or to predict future outputs from current outputs. Finally, it should not be possible to compute from one state previous states or previous outputs [7]. Therefore, while in general the state could be a counter and the complexity rest in the output function, here the state transition function must be one-way. Thus, so-called chaotic functions seemed attractive for use as state transition functions.
The logistic map [1] is a function f : If parameter a is close to 4, then the range [0; a/4] approaches the domain, so the function can be applied iteratively, and it exhibits chaotic behaviour. As f (0) = f (1) = 0, we will identify the two and consider only interval [0; 1) in the sequel. A variant of the logistic map has actually been used in a stream cipher [15].

Finite-length representations for real numbers
For non-integral numbers, the commonly used representation are floating-point numbers [14], which have been standardized by IEEE in 1985, and revised twice since then [8].
Thus, a floating-point representation is comprised of a power of 2 times a mantissa, i. e. (i, M) = 2 −i ⋅ 0.M. In order to avoid redundancies, i. e. to ensure that a representation (i, M) ∈ I i , the mantissa in all intervals except I b−1 must be in the range [0.5; 1), i. e. the first mantissa bit must be 1. This 1 is called implicit as it need not be stored explicitly. Those numbers are called normalized, The numbers in interval I b−1 , which comprises 0, are called de-normalized.
As IEEE standard 754 also covers negative numbers and numbers larger than 1, a floating-point representation with a sign bit s, mantissa bits M of length m and an exponent E of e bits comprises t = s + m + e bits and represents number where b = 2 e−1 − 1 is the bias and ⟨⋅ ⋅ ⋅⟩ denotes binary representation. Exponent value E = 2 e − 1 is used to represent infinity, while exponent value E = 0 is used to represent de-normalized numbers IEEE Standard 854 [9] generalizes floating-point representations to radix or basis B > 2. If the basis is a power of 2, then each mantissa digit in the range {0, 1, . . . , B − 1} can be represented by log 2 B bits. The "implicit 1" cannot be used anymore for higher radix, as the first mantissa digit of a normalized representation will be different from 0, but can have B − 1 possible values, and also the binary representation of this digit has no bit which is 1 for all those values.
Also fixed-point representations [14] can be used to represent interval [0; 1). For other real numbers, a combination of their integral part and their rational part is used for representation. A fixed-point representation M of m ὔ bits for interval [0; 1) represents ⟨M⟩ ⋅ 2 −m ὔ , i. e. the representations are equally spaced over the interval. Addition of fixed-point representations is done just as for integers. Multiplication of fixed-point representations is done by multiplying the mantissas, plus a final right shift of m ὔ bits that realizes the multiplication by 2 −m ὔ inside the brackets of the following formula: Also a rounding can be applied similar to floating-point numbers.
Please note that for m ὔ = t, i. e. similar memory requirements for both representations, the fixed-point representation has denser representations than floating-point representation in the interval [2 −m+1 ; 1), while floatingpoint representations are denser in interval [0; 2 −m ).

Related work
The investigations by Keller and Wiese [11] and Kneusel [13] focus on the influence of parameter a: even for two close values of a ≈ 4, the period length can vary notably. They use floating-point numbers, but do not report different rounding modes or order of evaluation when doing the two multiplications in a ⋅ x ⋅ (1 − x). The thesis by Luber [16] considers different ordering of multiplications for fixed-point representations, but does not consider rounding modes as well. Furthermore, Luber does not investigate period length but statistical properties of output bits derived from state values.
An extension of fixed-point format is described in [4], where two different exponents −m 1 and −m 2 < −m 1 can be used. The choice of exponent is encoded in one bit. This could also be considered as a floating-point format with one exponent bit and mixed radix, and as such be a special case of our proposal in Section 4. In [21], the authors apply signed digit numbers and allow the use of different radix for different parts of the mantissa, yet they do not generalize their approach from signed digit representations to floating-point or fixed-point representations.
[2] point out that higher radix floating-point format is advantageous in reconfigurable hardware, but that normalization (ensuring that not all of the log 2 B uppermost mantissa bits are zero after an operation) often is not applied there. In contrast, our implementation applies normalization.
In [6], the necessary mantissa length (bitwidth) for a given application is analyzed. This approach might be used as a post-processing to our approaches in Sections 3.4 and 4.

Challenges
In this section, we will investigate challenges for replication of execution results, when evaluating the formula a ⋅ x ⋅ (1 − x). Those challenges illustrate that even on such microscopic scale, replication is difficult, and that therefore, replication of experiments on a larger scale must work with the assumption that certain differences will be inevitable.

Order of computations
For real numbers, the order of above multiplications does not matter, as multiplication is commutative and associative. Yet, associativity is not given anymore for finite number representations, because rounding will occur. Hence, depending on the choice, there are three possibilities for choosing the arguments of the first multiplication, when we assume commutativity. Even when the source code is given, we still cannot be sure about the order, as compilers have the possibility to re-order operations, e. g. to improve performance for a given processor architecture. Yet, as long as the compiler product, its version and its compilation flags are specified in detail, this can be avoided. However, when the source code is not given, the implementer of the replication code is free to choose order, and over time several different implementation orders have been seen by the author. While the differences between results from different orders will be tiny, the resulting state graphs will have notable differences, cf. Subsection 3.4.
As an example, we assume a representation with 8 significant binary digits and take a = 3.921875 = 4 − If we multiply x and 1 − x first, this results in 0.21484375 = 0.00110111, which will not be rounded. Further multiplication with a results in 0.84259033203125 = 0.11010111101101 which will be rounded to 0.11011000 = 0.84375.

Rounding modes
Even if the former challenge has been mastered, there is more in the multiplications. So far, we have assumed that rounding goes to the nearest representable number. Yet IEEE 754 standard defines several rounding modes [8], i. e. round towards zero (the rounding applied so far), round towards infinity (always round up to next representation), and round to nearest, i. e. depending on the value of the first digit to be cut away, one rounds up or down.
All the rounding modes are nowadays implemented in hardware by IEEE754 compliant floating-point units in processor cores. The choice of the rounding mode is done either by the BIOS on power up, but might also be changed by the operating system or even by the application program, e. g. via function fesetround() in fenv.h.
As an example, when multiplying x = 0.6875 and 1 − x first and then multiplying with a = 3.921875, the result is 0.11010111101101. When rounding to nearest is used, we get 0.11011001, while the result would change to 0.11011000 when rounding down. Thus, the difference between the two results is 2 −8 , which is a notable 0.33 %. Also here, this difference will lead to different state graph structures, e. g. with respect to cycle lengths.
Unfortunately, rounding modes normally are not mentioned when experimental conditions are described in articles, so everyone seems to assume that round to nearest is used.

Number representation
So far, we have used the straightforward number representation for reals, i. e. floating-point numbers. However, this has some disadvantage. First, two bits of the available format are lost: we are not using negative numbers, as we are focusing on interval [0; 1), and we only use half of the possible exponent values. Second, the representable numbers in IEEE 754 are not evenly distributed over the interval [0; 1) as mentioned in Section 2.
In contrast, a fixed-point number representation comprising 32 bits has 2 32 equally spaced representations in [0; 1). Fixed-point arithmetic is simple when we choose uint32 to implement this scheme. Addition (and subtraction) is done as usual, at least as long as no overflow occurs, which we can exclude here. Multiplication is done by converting the arguments to uint64, doing the multiplication, shifting by 32 bits to the right, and converting back to uint32. Please note that in addition to a round towards zero (by shifting 32 bits to the right), also other rounding modes can be realized by considering the highest bit to be rounded away. Rounding up then is realized by incrementing after the right shift.
Please note that this realization, which has been presented and used in [20], prefers a particular order of multiplications, i. e. first multiply x with 1 − x, to avoid intermediate results larger than 1. However, also other orders are possible with marginal additional effort.

Experiments
As floating-point implementations have been investigated in the past [11,13], we will first focus on fixed-point implementations in our experiments. For 63 equi-spaced values of a, ranging from 4 − 2 −12 = 3.999755859 to 4 − 2 −17 = 3.999992371, we explore the differences between the three different orders of multiplications with respect to the period length of the resulting generator. As there the state graph comprises multiple weakly connected components with differing periods, in each experiment we use 32 equispaced states in the range 2 −6 to 0.5 to start from and report the maximum period. As we expect the largest connected component, which typically also comprises the largest (or second largest) cycle, to comprise about 75 % of all states [5], the chance to miss it with 32 starting points is only 2 −64 which we can tolerate. While the start states are not chosen randomly, they are chosen independently of the particular a, i. e. their choice should be as good as random. The results are depicted in Table 1.
Please note that differences between the rounding modes have been investigated in [10], and that in all except one scenario, round to nearest lead to longest periods. Hence we will use this rounding mode in the rest of this article.
The order where x(1 − x) is computed first leads to shortest cycles on average. Please note that standard deviation is very high, so the average is not too meaningful. However, due to the chaotic nature of the logistic map, this To compare with floating-point implementations, we used an implementation based on IEEE754-compliant floats with round to nearest. For them (cf. Section 5 and Table 3), computation of x(1 − x) first leads to longest cycles, which however range from 91 to 7,802, with an average of 3,172, i. e. are much shorter than cycles from generators with fixed-point implementations, which have the same state size of 32 bits.
This comparison also serves to illustrate a tradeoff that is possible. We might shorten the number of bits used in the fixed-point representation. Table 2 illustrates how average cycle lengths (averaged over 63 values of a as before, with round to nearest an 1 − x multiplied last) are reduced by using shorter number representations. We see that after shortening by 8 bits, i. e. only using 24 bit fixed-point representation, the cycles are about as short as for single precision, i. e. 32 bit, floating-point representation. Thus, if state space is very tight, e. g. in very small devices such as RFID devices, this opens a chance to reduce hardware footprint, i. e. number of gates necessary, and reduce energy consumption for multiplications, targeted towards security needs and energy constraints -similar to other approaches of approximate computing.

Comparing different number formats
Thus, to store a number in this format, we need e = ⌈log 2 d⌉ bits to represent the interval index, and m bits to represent the mantissa (or m − 1 in case of binary floating-point representations, see below). We start with fixed-point representations. They only comprise one interval I 0 , i. e. d = 1, and thus the number of exponent bits is zero. The mantissa bits represent a multiple of 2 −m , i. e. the represented number is r(0, M) = ⟨0.M⟩ = ⟨M⟩ ⋅ 2 −m .
In single precision floating-point representations, d = 127, and interval I 0 contains the denormalized numbers. The exponents form a sequence e i = i − 126, i. e. successive exponents only differ by 1. To ensure that r(i, M) ≥ 2 e i−1 for i ≥ 1, i. e. that r(i, M) ∈ I i , it follows from Eq. (2) that ⟨0.M⟩ ≥ 0.5, i. e. that the uppermost bit of mantissa M must be set (so-called implicit 1). Thus, only m−1 bits must be stored, or put otherwise, if we store m bits, we really work with a m + 1-bit mantissa.
Another example of a possible format is a separation in d = 4 intervals, i. e. using e = 2 exponent bits, with (e i ) i = (−32, −16, −10, 0). Thus, no particular regularity is necessary beyond the restrictions above. The mantissa is shorter than for fixed-point representations, but longer than for floating-point representations, so a possible advantage in cycle length might arise.
If all differences e i − e i−1 are equal, then we can consider the representation to be designed with respect to higher radix 2 e i −e i−1 , if they are not equal, they can be considered multi-radix representations.
Hence, for given total storage size of t = e + m bits, there are different numbers d = 2 e of possible exponents, leading to different mantissa lengths m, and for each d the exponents can have a number of possible distances leading to different number formats, which will in turn lead to different cycle lengths in the logistic map application. In this configuration space, we can search for the representation which is best with respect to a particular application.
Alternatively, one can fix the required performance, i. e. cycle length in our application, and search for the representation with shortest storage size that achieves that performance, taking up the view of approximate computing [19].
Moreover, from all combinations of performance and storage size, the user can draw a plot and identify the Pareto front, i. e. the curve of Pareto-optimal combinations, which is the (typically small) number of combinations that the user can choose from depending on her needs [3].
In general, if e i − e i−1 ≥ m, then all representations of Eq. (2) (if not all mantissa bits are 0) are in I i . Otherwise, the uppermost e i − e i−1 mantissa bits may not all be 0 simultaneously. As there is no similar trick as the implicit 1 if e i − e i−1 ≥ 2, a fraction of 2 e i −e i−1 of the representations does not belong to I i but to I j with j < i, i. e. there is a certain redundancy, i. e. reduction in the total number of distinct representations.
To explore the spectrum of possible number formats, we have implemented a custom library where d, M and the sequence of e i can be specified for a representation type. 1 We investigate total lengths t from 32 down to 24, with mantissa lengths ranging from t (resulting in e = 0) to 24 (resulting in e = t − 24). For each combination of e and m we explore sequences of e i with regular distances e i − e i−1 from the set {1, 2, 4, 6, 8, 10}. So far, only round-to-zero is supported and multiplication with a = 3 + ε is performed last, by first adding x(1 − x) to itself, by computing εx (1 − x) and adding x(1 − x), and finally adding the two intermediate results.
For each number representation, we choose again 63 different values of a as before, and compute the minimum, maximum and average cycle lengths (plus standard deviation).
Like in previous experiments, we see that standard deviation is quite large (ranging 49 % to 69 % of the average) and does not get better if we take larger number of samples. For almost all (31 of 33) combinations of e and m, we find that e i − e i−1 = 1 leads to largest average cycle lengths. In the two exceptions, e i − e i−1 = 2 performs better: in one case (e = 1, m = 29) the difference is only 0.32 %, yet in the other case (e = 1, m 3 0) the difference is 30.9 %, i. e. notably, mainly driven by a really large maximum cycle length achieved. Thus that case can be considered an outlier.
In the following, we further investigate the e, m combinations with distance leading to maximum average cycle length. Figure 1 depicts that the average cycle length is mainly determined by the number of mantissa bits, for each m it is largely constant over the different exponent 1 Available from authors upon request.  lengths e, and the curves are ordered by growing m. Please note that the cycle lengths for e = 0 are taken from Table 2 where different rounding mode and multiplication order are applied, so that the step from e = 0 to e = 1 is not completely fair, which might explain the unusual shape for m = 31. Figure 2 depicts that for a fixed exponent length e, the average cycle length grows with total length t = e + m, i. e. with mantissa length m, and the curves are ordered by shrinking exponent length e. The only exception is a crossing of the curves for e = 1 and e = 2 at the right end of the curves. Figure 3 completes this view by plotting curves for each total length t = e + m. We see that each curve decreases with growing e, i. e. shrinking m, and the curves are ordered by t.
In summary, we find that for our application, distances e i − e i−1 = 1 are best, as in binary floating-point representations, and that fixed-point representation leads indeed to longest average cycle lengths, even when we consider the large space of representations between (binary) floating-point and fixed-point representations. 4 − 2 −12 , the standard deviation does not shrink substantially, e. g. from 1610 to 1523 for symmetric float representations and multiplication with a last. Yet, in some variants standard deviation even increases, e. g. from 1630 to 1688 for conventional floats and multiplication with 1 − x last. The biggest difference is that the average cycle length for ordinary floats where multiplication with a is done last reduces to 2459.7, so that the range of averages is much closer together.

Conclusions
We have presented challenges that occur with respect to replicability in security or cryptographic applications that involve numeric parts originally devised for real numbers. As example, we have chosen the logistic map, which has been one proposal for a chaotic pseudo-random number generator. We have demonstrated by experiments the relevance of considering the differences that possibly result from not arbitrating the ambiguities underlying the above challenges.
As similar problems might occur in other security applications as well, we derive the following recommendations. First, the rounding mode used in the original experiments should be found out by the experimenters and should be reported. Second, the order of computations used in the original experiments should be reported as well, and possibly should be modified by the first experimenters to check if notable differences in result can be the consequence. Finally, different number representations should be considered by experimenters, and results on their influence and reasons for the final choice of representation should be reported.
Moreover, we have investigated higher and mixed radix representations in between fixed-point and binary floating-point representations, plus symmetric representations, and demonstrated how their tradeoff between storage size and performance (in the sense of period length in the application) can be used by users to tailor their application to their needs, thus connecting this investigation to the field of approximate computing.
Future work will comprise a deeper analysis of our experimental results. Also, the scope of numerical challenges will be extended from chaotic random numbers to other branches of IT security. Moreover, the introduction of redundancy by higher radix representations, e. g. by skipping normalization [2], opens the path for information hiding, e. g. storage covert channels [18], in such representations, which we plan to investigate further.
Funding: This work is supported by project SIMARGL (Secure intelligent methods for advanced recognition of malware, stegomalware & information hiding methods, https: //simargl.eu), which has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 833042.