As observed by Berezowski, ⁶ the distribution of prime numbers exhibits chaotic behavior, making their prediction and arrangement highly non-trivial. This was the first thought for selection of a proper seed to generate a sequence of pseudo-random numbers. A seed can be selected as an arbitrary prime number whose magnitude depends on the magnitude of the random numbers that need to be generated. Therefore, we define the seed to be: p 0 = select _ closest _ next _ prime ( random _ integer ( 100 , max ( 1000 , 2 ⌊ m 4 ⌋ ) ) )

Here, m is the number of bits used to represent the integers in the sequence. If m = 32 , then the numbers in the sequence will be bounded by 2 m , meaning they can take values between 0 and 2 32 − 1 (for 32-bit integers). Thus, this function takes an arbitrary integer value between 100 and the maximum value between 2 ⌊ m 4 ⌋ and 1000, and chooses the closest prime number greater than that integer, which gives a wide range of prime numbers for the algorithm to choose from. We use 2 ⌊ m 4 ⌋ as an upper limit to make the initial seed comparatively lower in magnitude to the largest magnitude of random number that the algorithm can generate, as shall be discussed later.

The initial idea for developing this algorithm was to reduce, or, if possible, remove the repeating periodic cycles that appeared in the use of LCG (Linear Congruential Generators), ⁷ one of the most primitive algorithms to generate pseudo-random numbers. This led to the thought of using a non-constant entropy term which would oscillate between increasing and decreasing the current number that is being generated. b = ( − 1 ) i · i · p i − 1

Here, one would notice the presence of the term p i − 1 , which speaks a lot about our algorithm. We do not simply discard the initial seed after generating the first term of the sequence. We keep our notion of the distribution of prime numbers being chaotic in nature, ⁶ and hence use a term p i which is the next prime number after p i − 1 . Thus, p 1 would be the next prime number after p 0 , p 2 would be the next prime number after p 1 and so on.

This, again, is a very simple idea, inspired from LCG ⁷ where the first component of the algorithm is of the form x i = f ( x i − 1 ) , where the previous term is simply multiplied by a constant. We can modify it so that the term has a bit more randomness in the following way: x i = f ( x i − 1 , p i − 1 ) , where the function describes the multiplication of the two terms from its previous iteration.

It is important to mention that every term that the algorithm generates is always limited to m bits. This is achieved by taking a mod 2 m operation on every term that is generated. This gives us the advantage of limiting the number of bits to the necessary amount, i.e., m , as per the use case. In order to add an extra layer of complexity, with the intention to further reduce the chance of periodic cycle formation in the sequence, a bit-wise shifting operation is introduced. The value of this shift itself is not predetermined and depends on a modulus of the iteration counter such that, shift _ amount = imod ⌊ m 2 ⌋

This ensures that the number of bits shifted always falls in the range [ 0 , m 2 − 1 ] . Thus, the ( i − 1 ) th term is right-shifted according to the above idea and can be represented as follows: x i − 1 ≫ ( imod ⌊ m 2 ⌋ ) ≡ ( x i − 1 2 ( imod ⌊ m 2 ⌋ ) )

Since x i − 1 is of m bits, the term produced above is of m-bits. Again, on applying a modulus on the multiplier term, we can also guarantee it to have the same number of bits. With this knowledge, we can easily apply an ‘ExclusiveOR’ or ‘XOR’ operation on these two terms to get a fairly convoluted string of bits. Marsaglia stated that in the context of random number generation, the use of XOR operations has been shown to enhance the efficiency and unpredictability of pseudo-random number generators. As demonstrated by him in 2003, ⁸ XORShift random number generators utilize bit-wise XOR operations to produce sequences of numbers that exhibit desirable statistical properties, making them suitable for a variety of applications. Therefore, these operations can be combined to form the whole m-bit term as follows: x i = ( ( f ( x i − 1 , p i − 1 ) mod 2 m ) ⊕ ( x i − 1 2 ( imod ⌊ m 2 ⌋ ) ) )

The XOR operation is mathematically represented by the ⊕ symbol.

The final output sequence at each step i , denoted by x i , is generated using a combination of Modular Arithmetic, bit-wise manipulation (XORshift), and Entropy Modulation as is a culmination of the ideas discussed above and can be represented as follows: x i = ( ( ( x i − 1 · p i − 1 ) mod 2 m ) ⊕ ( x i − 1 2 ( imod ⌊ m 2 ⌋ ) ) + ( − 1 ) i · i · p i − 1 ) mod 2 m

Where: •

i denotes the iteration counter,

The algorithm for this can be given as follows:

seed ← next prime from random number in range [ 100 , max ( 1000 , 2 m 4 ) ] Initialize sequence as an empty list current ← seed

next_num ← ( current × seed) mod 2 m next_num ← next_num ⊕ ( current ≫ ( imod m 2 ) ) next_num ← ( next_num + ( − 1 ) i × i × seed ) mod 2 m Append next_num to sequence seed ← next prime after seed current ← next_num

return sequence

The algorithm itself is defined in a way to function well with an arbitrary seed. The choice of seed as mentioned in equation ([eq:seed]) shows that this algorithm trivially has the notion of an arbitrary seed selection system, thus, freeing it from the problem of selecting seeds that do not yield good results. It is important to note that the algorithm always selects a prime number as its seed, by definition, even if a non-prime integer is given as an input. In case a non-prime integer is given as an input, the algorithm selects the closest prime number greater than that integer as the seed p 0 . The range of [ 100 , max ( 1000 , 2 ⌊ m 4 ⌋ ) ] has only been introduced as a convention to reduce the computational overhead of just selecting a seed. The only constraint that we’ve faced is when the value of the seed p 0 is very small compared to 2 m . In this case, the first few numbers of the generated sequence are also smaller values and thus shift away from the notion of randomness in the whole distribution. In such a case, the lower bound of the range can be chosen to be a higher number masted on the magnitude of m .

A common problem of PRNG’s are the occurrence of periodic cycles in the sequence. Since the algorithms are deterministic, once a random number is generated that is equal to some previous number in the sequence, the whole algorithm starts repeating itself, or atleast has a chance to do so. The is the most common problem for using LCG’s, ⁷ where one must be really careful in the selection of seeds and parameters for maximal cycle formation. However, in MAXEM, the possibility of a periodic cycle formation has simply been eliminated. The same number may be generated more than once, but the scenario will not lead to the repetition of the whole sequence. We can see that, from equation ([eq:entropy]), the term x i depends directly on the product of the terms i and p i − 1 . It is going to be highly unlikely for bmod 2 m [equation([eq:entropy])] to generate the same term as it had before, because both i and p i are changing. Thus, the chances that the whole sequence repeats is almost equal to null. This has also been verified by simulation after running the algorithm many times, with not a single case of formation of periodic cycles.

As discussed by Ross in Introductory Statistics (Third Edition), ¹¹ non-parametric hypothesis tests, such as the Wald Wolfowitz run test for randomness, are useful when no specific distribution is assumed. The run test evaluates the randomness of a sequence by analyzing the occurrence of ‘runs’—sequences of consecutive identical elements. A significant deviation in the number of runs from what is expected under randomness suggests the sequence may not be random. The Run test has been carried out on the whole sequence with the following hypothesis:

We have failed to reject the Null Hypothesis ( H 0 ) for every single time the algorithm has been run, thus concluding that the sequence thus generated by MAXEM can be considered random.

From the example above ( section 3), here is the sequence of letters that were generated and considered for the run test: [ l , l , l , m , l , m , l , m , m , l , m , l , l , l , m , l , m , m , l , m , l , m , m , m , m , l , l , l , m , m , … ]

Here, l : x i < Median ( { x n } ) and m : x i ≥ Median ( { x n } ) . We can see that the sequence has neither way too many runs, nor way too less of them, suggesting that there is very little chance of patterns forming in the sequence. If u be the total number of runs in a sequence and n l and n m be the number of l ’s and m ’s in the sequence respectively, we can calculate: μ ( u ) = n + 2 n l · n m n σ ( u ) 2 = [ μ ( u ) − 1 ] [ μ ( u ) − 2 ] n − 1 Z = u − μ ( u ) σ ( u ) where μ ( u ) is the expected number of runs and σ ( u ) is the standard error. If the Z-score( Z ) is such that | Z | ≤ 1.96 (taking a 5% significance level of the standard normal variate), then we accept the Null Hypothesis ( H 0 ), else reject it.

Martin-Löf ¹² introduced a formal way to define random sequences using algorithmic complexity. While a sequence might be considered random as a whole, this doesn’t mean that every subsequence will also be random. The randomness of a full sequence doesn’t carry over automatically to its parts, as a subsequence might still follow some patterns depending on how it’s selected.

With this in mind, we tested the sequence generated by MAXEM following a simple algorithm as given below.

Choose an initial value for len Choose an arbitrary number y ∈ [ 0 , n − len ] Select subsequence of length len starting at position y Conduct a run test on each subsequence. Let z = number of times the run test is satisfied Compute success _ percentag e len = z 1000 × 100 Increase len gradually, starting with small intervals and increasing the gap for larger values.

Figure 2 shows the results of local randomness test for len values like 30, 40, 50, 100, 150, 200, 400,…, upto 20000. The high success results prove how the sequence generated by MAXEM is not only random globally, but also locally.

These analyses indicate that MAXEM provides a high level of randomness with good scalability. Benchmarks like TestU01, designed by L’Ecuyer and Simard, ¹³ provide a comprehensive framework to evaluate PRNG performance, and can be used to conduct further analyses on the robustness of the MAXEM algorithm.

The worst-case time complexity of the algorithm is primarily dominated by the modular multiplication and can be given as O ( n · mlogm ) if faster and more optimized multiplication such as Karatsuba’s algorithm ¹⁴ is used. This algorithmic complexity does not take into account the complexity of the separate function that finds the nearest prime number for the seed, as it can be (and should be) done separately or by parallel methods. We can see that this is computationally more expensive that other common PRNGs succh as the Mersenne Twister algorithm. ³ However, the trade-off is justified by increased entropy and unpredictability.

The space complexity of MAXEM remains to be O ( m ) , as is the case with most other common PRNGs.

The major goals of the MAXEM algorithm are to negate the formation of periodic cycles in a sequence of random numbers and also preserve the quality (randomness) of the sequence with the selection of an arbitrary seed. This makes MAXEM stand apart from the other traditional PRNGs like the Linear Congruential Generator. ²

This algorithm, given a high value of m , also has great security. Due to the intrinsic features of this algorithm, reverse engineering this without a good guess on the final seed would be difficult as one would need to guess both the iteration number and the seed, given that the seed is a prime number of x bits, where x is quite large. This task is not impossible but quite arduous. Hence, this algorithm can also be great for use in cryptography ¹⁵ and other security measures.

The ability to adjust the bit range 2 m offers significant flexibility. By selecting appropriate values for m , the generator can be tuned for different use cases, whether it’s cryptographic applications or scientific simulations.

Despite the advantages of the MAXEM algorithm, the inherent high computational complexity can be a limiting factor and make the generation of sequences for large values of m and n very slow to the extent that it becomes unsuitable for real-time applications. However, taking advantage of these modern techniques such as parallel computing can significantly help in the process of improving the efficiency of the system without a sacrifice of quality.

The authors further declare that the data used in the work was not obtained from any external sources. All data supporting the results in this study were generated by the authors using the MAXEM algorithm, using their own code. The full dataset (random sequences, statistical measures, and plots) and the implementation code are publicly available on Zenodo at: https://doi.org/10.5281/zenodo.15631843. ¹⁶

Data are available under the terms of the Creative Commons Attribution 4.0 International license (CC-BY 4.0).

The authors also declare that this study does not involve any human or animal participants, data, or tissue.