The Middle Squares Method: The First PRNG

And also why you shouldn’t use it

Computers are great at many things, but one thing that they cannot do is generate random numbers on their own. However, computers can fake generating random numbers by using a pseudo-random number generator (PRNG). Unlike true RNGs, PRNGs use algorithms to create seemingly random numbers from an initial starting value, called a seed. PRNGs are deterministic, i.e., repeatable: the same PRNG always outputs the same list of random numbers when starting from the same initial state from the same seed. While using insecure PRNGs in cryptographic applications is a security risk,¹ they are still useful for other areas like procedural generation because of their reproducible outputs. The first ever PRNG, the Middle Squares Method, was developed by John von Neumann in 1946, and later published in 1951.² This PRNG was a crude but important initial step into designing the more complex ones that we use today.

import datetime
import sys

import matplotlib as mpl
import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import pyvis
from pyvis.network import Network

PRNGs as Finite State Machines

We can model a PRNG as a finite state machine (FSM), an abstract computing model that can be in one of a finite number of internal states at a time. To generate pseudorandom numbers, the PRNG performs a transition function \(f\) that changes its current state \(s_i\) to the next state \(s_{i+1}\), then performing an output function \(g\) to convert the new state \(s_{i+1}\) into its corresponding output \(u_{i+1}\). The probability distribution \(\mu\) describes how the initial state \(s_0\) is chosen from the set of all states \(\mathscr{S}\), and the set of all outputs the PRNG can make is \(\mathscr{U}\). Together, these five objects form the mathematical structure \((\mathscr{S}, \mu, f, \mathscr{U}, g)\) that defines a PRNG.³

The middle square function forms the core of the Middle Squares Method; this is the transition function that determines the PRNG’s next state from its current state. First, it takes an input \(a\)-digit number \(s\) in base \(b\) and squares it. It then pads the result \(s^2\) until it contains \(2a\) digits; it chooses the middle \(a\) digits as the next state. For example, if we want to run the Middle Squares Method in decimal (base 10) with \(a = 4\), the next state after 5678 would be 2396; 5678 squared is 32,239,684, which has 2396 as the middle four digits. The value of \(a\) must be an even number for the function to work; if \(a\) was odd, then there would be ambiguity in how to select the middle \(a\) digits.⁴

def middle_square(state: int, digits: int, base: int = 10) -> int:
    """Perform the middle square function on the input state."""
    if digits % 2 != 0:
        msg = "Number of digits must be even."
        raise ValueError(msg)
    if len(np.base_repr(state, base)) > digits:
        msg = f"State {state} is too large for {digits} digits."
        raise ValueError(msg)
    state_squared = np.base_repr(state**2, base).rjust(2 * digits, "0")
    return int(state_squared[digits // 2 : digits // 2 + digits])

The Middle Squares Method’s output for a given state is equal to the state itself, so the output function \(g\) is the identity function. We can determine that \(\mathscr{S}\) and \(\mathscr{U}\) are both equal to the set of all \(a\)-digit numbers in base \(b\). So, to recap:

PRNG Component	Description	Middle Squares Method
\(\mathscr{S}\)	set of all PRNG states	set of all \(a\)-digit numbers in base \(b\)
\(\mu\)	probability distribution on \(\mathscr{S}\) determining initial state \(s_0\)	User randomly selecting a seed
\(f\)	transition function, converts the current state to the next state	The middle square function
\(\mathscr{U}\)	set of all PRNG outputs	set of all \(a\)-digit numbers in base \(b\)
\(g\)	output function, converts a state to its output	identity function

The Middle Squares Method in Action

With these five elements, we can now define the Middle Squares Method.⁵ Let’s generate some numbers to see what its output looks like.

seed = 20230928
ms_outputs = [seed]
for _ in range(30):
    ms_outputs.append(middle_square(ms_outputs[-1], 10))

print(ms_outputs)

[20230928, 4092904477, 8670578466, 9309350629, 91336626, 3423792490, 3550145804, 5352296588, 787659164, 4069586331, 5329054620, 8231429433, 4305104587, 9255050084, 9520573484, 3194642438, 7403066705, 3966386795, 2242075503, 9025611527, 6634363152, 7744326153, 5875640397, 1500748583, 2463093765, 8308951818, 6803138455, 6928378997, 4355260707, 2958259381, 2985652745]

Other than the seed, these look random enough, right? Well, there’s a pretty glaring issue with the Middle Squares Method that we can see by analyzing its structure. Since the Middle Squares Method is a finite state machine, we can use a state diagram to show how the PRNG’s internal states are linked through state changes. We can represent the state diagram as a directed graph (digraph) where each state corresponds to a vertex and each transition with the middle square function corresponds to a directed edge pointing from the starting state to the next state. Let’s look at a simple implementation of the Middle Square Method for \(a = 2\) and \(b = 10\).

DIGITS = 2
BASE = 10

ms_digraph = nx.DiGraph()
ms_digraph.add_nodes_from(range(BASE**DIGITS))
graph_edges = [(s, middle_square(s, DIGITS, BASE)) for s in range(BASE**DIGITS)]
ms_digraph.add_edges_from(graph_edges)

# color nodes by which connected component the corresponding undirected graph
# belongs to
ms_graph = ms_digraph.to_undirected()
for comp_idx, component in enumerate(nx.connected_components(ms_graph)):
    for node in component:
        ms_digraph.nodes[node]["group"] = comp_idx

# Plotting with matplotlib and networkx
fig, ax = plt.subplots(figsize=(16, 9))
pos = nx.nx_agraph.pygraphviz_layout(ms_digraph, prog="dot")
nx.draw_networkx(
    ms_digraph,
    pos,
    with_labels=True,
    node_color=[f"C{ms_digraph.nodes[n]['group']}" for n in ms_digraph.nodes],
)
ax.set_title("Middle Squares State DiGraph")
plt.savefig("middle_squares_net.svg")
plt.close()

# Interactive visualization with pyvis
for i in range(BASE**DIGITS):
    label = str(i).rjust(DIGITS, "0")
    title = f"{label}² = {str(i**2).rjust(2*DIGITS, '0')}"
    ms_digraph.nodes[i]["label"] = label
    ms_digraph.nodes[i]["title"] = title
    ms_digraph.nodes[i]["shape"] = "circle"
ms_net = Network(directed=True, notebook=True, cdn_resources="in_line")
ms_net.from_nx(ms_digraph)

Static Image (matplotlib + networkx)

Middle Squares State DiGraph.

Interactive Graph (pyvis)

ms_net.show("middle_squares_net.html", notebook=True)

middle_squares_net.html

Middle Squares State DiGraph (Interactive Version). Click, drag, and zoom with the mouse to see the relationships between PRNG states.

The Middle Squares Method’s Weakness

If you start at any state and trace the state transitions by following the directions of the arrows, you eventually fall into one of four states that transition to themselves (0, 10, 50, and 60) or into a state that loops back to one you’ve already visited (24 and 57). When either of these situations occur, the Middle Squares Method will become stuck, constantly outputting the sequence of digits forever.⁶

To keep a cryptographic system safe, we have to assume that an attacker outside the system has access to everything except the system’s secrets (like keys) - this includes the internal structure of the RNGs used in the system.⁷ Since the Middle Squares Method always falls into a pattern of constantly repeating the same sequence of digits, that makes it an easy target for an adversary to attack for any cryptographic system it’s a part of. As von Neumann himself put it in the paper describing the Middle Square Method:

“Anyone who considers arithmetical methods of producing random digits is, of course, in a state of sin. For, as has been pointed out several times, there is no such thing as a random number – there are only methods to produce random numbers, and a strict arithmetic procedure of course is not a method.”

–John von Neumann (1951), Various Techniques Used in Connection With Random Digits⁸

Package Versions & Last Run Date

Python 3.13.5 (main, Jun 25 2025, 18:55:22) [GCC 14.2.0]
Matplotlib 3.10.6
NetworkX 3.5
NumPy 2.3.3
PyVis 0.3.2
Last Run 2025-10-04 12:37:01.363688

References

L’Ecuyer, Pierre. “History of Uniform Random Number Generation.” In 2017 Winter Simulation Conference (WSC), 202–30. Las Vegas, NV: IEEE, 2017. https://doi.org/10.1109/WSC.2017.8247790.

Neumann, John von. “Various Techniques Used in Connection With Random Digits.” In Monte Carlo Method, 12:36–38. Los Angeles, California: National Bureau of Standards, 1951. https://catalog.hathitrust.org/Record/003239351.

Stipčević, Mario, and Çetin Kaya Koç. “True Random Number Generators.” In Open Problems in Mathematics and Computational Science, edited by Çetin Kaya Koç, 275–315. Cham: Springer International Publishing, 2014. https://doi.org/10.1007/978-3-319-10683-0_12.

Footnotes

1. Stipčević and Koç, “True Random Number Generators”.

2. L’Ecuyer, “History of Uniform Random Number Generation”.

3. L’Ecuyer.

4. Neumann, “Various Techniques Used in Connection With Random Digits”.

5. L’Ecuyer, “History of Uniform Random Number Generation”.

6. L’Ecuyer.

7. Stipčević and Koç, “True Random Number Generators”.

8. Neumann, “Various Techniques Used in Connection With Random Digits”.