Contents

Python Mathematical Computation Performance Benchmark

The cover image was generated by Google Gemini.

Tip
This project was written in October 2024.

While taking the graduate course Simulation Study, the instructor mentioned that the random numbers used in modern computers and programming are mostly not truly random in the strict sense, but rather pseudo-random numbers. Pseudo-random numbers are generated through mathematical algorithms that compute the next value based on the previous state. Therefore, as long as the algorithm, initial value (Seed), and current state are known, the entire random number sequence can theoretically be reproduced. This raises an interesting question: if every pseudo-random number can be calculated from the previous value, where does the very first number come from?

True Random Numbers and Pseudo-random Numbers

To answer this question, we first need to understand that random numbers can generally be divided into two categories: True Random Numbers and Pseudo-random Numbers.

True random numbers are generated from unpredictable physical phenomena in the real world, such as thermal noise produced by hardware circuits, clock jitter, mouse movement patterns, keyboard input intervals, disk I/O operations, and network packet events. Modern operating systems continuously collect these sources of information to build an entropy pool, which serves as the initial source for secure random number generators. Since these sources are difficult to predict, true random numbers provide higher security and are widely used in cryptographic applications, key exchange, HTTPS, SSH, digital signatures, and other security-related fields. However, collecting entropy sources and performing security processing introduce additional computational costs, making true random number generation generally slower.

In contrast, pseudo-random numbers are entirely based on mathematical operations. They are not directly generated from natural phenomena; instead, they use deterministic mathematical formulas to repeatedly calculate the next value from the previous state. As a result, they form sequences with good statistical properties. Although the entire sequence can be reproduced if the initial value and algorithm are known, making them less secure than true random numbers, their advantages include high computational speed, strong reproducibility, and minimal memory requirements for generating extremely long sequences. Therefore, pseudo-random numbers are widely used in Monte Carlo simulation, numerical analysis, statistical simulation, computer games, and general scientific computing.

Algorithms for Pseudo-random Number Generation

Returning to the original question, since pseudo-random numbers depend on the previous state, where does the first number come from? The answer is the initial value (Seed). If users manually specify a Seed, such as random.seed(42) in Python, the program will generate exactly the same random number sequence every time it is executed, which is useful for debugging and reproducing experimental results. If no Seed is specified, Python obtains an initial state with sufficient entropy from the operating system, usually derived from the entropy pool maintained by the operating system, and uses it as the starting point of the pseudo-random number generator.

In other words, modern computer random number generation does not rely entirely on true random numbers. Instead, it uses a small amount of true random information as the initial seed and then applies high-speed mathematical algorithms to efficiently generate a large number of pseudo-random numbers.

One of the most common pseudo-random number generators is the Linear Congruential Generator (LCG). It generates the next value using one multiplication, one addition, and one modulo operation. Its mathematical representation is:

$$ X_{n+1} = (a X_n + c) \bmod (m), $$

where $X_n$ represents the current state and $X_{n+1}$ represents the next state. Due to its simple structure, LCG has very high computational efficiency. However, inappropriate parameter selection may result in short periods or poor statistical properties. Therefore, many improved variants have been developed.

One important variant is the Lehmer Generator, which retains only the multiplication and modulo operations from the original LCG:

$$ X_{n+1} = a X_n \bmod (m). $$

In 1973, P. A. W. Lewis and G. P. Learmonth studied the Lehmer Generator and proposed a set of parameters with good statistical properties. Therefore, some literature refers to implementations using these parameters as the Lewis–Learmonth Generator. Later, in 1988, Park and Miller published the famous paper Random Number Generators: Good Ones Are Hard to Find, recommending the use of $a = 16807$ and $m = 2^{31} - 1 = 2147483647$ as the parameters for the Lehmer Generator. They referred to this implementation as the Minimal Standard Generator. Due to its good statistical properties, simple implementation, and high computational efficiency, it has become one of the most well-known implementations of the Lehmer Generator and has been widely used as a classic example in programming languages, numerical computing libraries, and textbooks.

On the other hand, Python’s standard library random module uses the Mersenne Twister (MT19937) algorithm by default. Compared with traditional LCG methods, Mersenne Twister provides a much longer period of $2^{19937}-1$ and superior statistical properties. Therefore, it can generate higher-quality pseudo-random numbers while maintaining high computational efficiency, making it one of the most widely used general-purpose pseudo-random number generators today.

From the above discussion, we can see that pseudo-random number generation is essentially the result of a sequence of mathematical operations, and the computational methods used by different algorithms directly affect their execution efficiency. In fact, not only random number generators, but also various common mathematical functions in Python, such as addition, subtraction, multiplication, division, exponentiation, logarithms, trigonometric functions, and various probability distributions, require different levels of mathematical computation internally.

Therefore, the performance differences among different operations become an interesting topic worth exploring. In the following sections, this article will use Python’s timeit module to benchmark various common mathematical operations and random number generation methods, and compare the execution time required for each operation in practice.

Python Mathematical Performance Benchmark

To fairly compare the execution efficiency of various operations, this experiment uses Python’s standard library module timeit for performance measurement. Compared with directly using time.time() or time.perf_counter() to record the execution time of a single run, timeit repeatedly executes the specified function a large number of times while minimizing the influence of background processes, system scheduling, and timing errors. Therefore, it is also the method officially recommended by Python for microbenchmarking.

In this experiment, each operation is wrapped as an independent function, and timeit.timeit() is used to execute each function repeatedly for one hundred million times ($10^8$ iterations). The total execution time required to complete all operations is then recorded.

The tested operations include not only basic arithmetic operations such as addition, subtraction, multiplication, division, and modulus, but also common mathematical functions including exponentiation, logarithms, square roots, trigonometric functions, inverse trigonometric functions, and hyperbolic functions.

Furthermore, to compare the efficiency of different pseudo-random number generation methods, several common probability distributions provided by Python’s random module, including Uniform, Normal, Gamma, Beta, and Triangular distributions, are included in the benchmark. In addition, Linear Congruential Generator (LCG) and Lewis–Learmonth Generator are implemented and tested, and their performance is compared with Python’s built-in random.random(), which uses the Mersenne Twister algorithm internally.

Tip
This article compares the overall execution performance of Python programs in practice, rather than the theoretical time complexity of each algorithm. Since every function call, object access, and module-level function wrapper in Python introduces additional computational overhead, the measured execution time reflects not only the computational cost of the underlying algorithms but also the additional overhead introduced by the Python runtime environment. Therefore, the results may vary depending on the hardware configuration of the computer.

The following is the Python code used in this experiment.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
# version 1131002
import timeit
import math
import random
import pandas as pd

# calculate U(0, 1) pseudo-random number by Linear Congruential Generator (LCG)
class LCG:
    def __init__(self, seed):
        if not isinstance(seed, (int, float)):
            raise TypeError('seed must be a number (int or float).')
        self.state = seed

    def random(self):
        self.state = (1103515245 * self.state + 12345) % 2147483648
        return self.state / 2147483648
    
    def randint(self, low, high):
        return low + int(self.random() * (high - low))


# calculate U(0, 1) pseudo-random number by Learmonth-Lewis Generator (LLG)
class LLG:
    def __init__(self, seed):
        if not isinstance(seed, (int, float)):
            raise TypeError('seed must be a number (int or float).')
        self.state = seed
        
    def random(self):
        self.state = (16807 * self.state) % 2147483647
        return self.state / 2147483647
    
    def randint(self, low, high):
        return low + int(self.random() * (high - low + 1))

# function
def test_add(): return a + b
def test_sub(): return a - b
def test_mul(): return a * b
def test_div(): return a / b
def test_mod(): return a % b
def test_pow(): return math.pow(a, b)
def test_abs(): return abs(a)
def test_log(): return math.log(a)
def test_log10(): return math.log10(a)
def test_exp(): return math.exp(a)
def test_sqrt(): return math.sqrt(a)
def test_square(): return a ** 2
def test_sin(): return math.sin(a)
def test_cos(): return math.cos(a)
def test_tan(): return math.tan(a)
def test_asin(): return math.asin(a / b)
def test_acos(): return math.acos(a / b)
def test_atan(): return math.atan(a)
def test_sinh(): return math.sinh(a)
def test_cosh(): return math.cosh(a)
def test_tanh(): return math.tanh(a)
def test_asinh(): return math.asinh(a)
def test_acosh(): return math.acosh(b)
def test_atanh(): return math.atanh(a / b)
def test_uniform(): return random.uniform(0, 1)
def test_normal(): return random.gauss(0, 1)
def test_randint(): return random.randint(0, 100)
def test_exponential(): return random.expovariate(1)
def test_gamma(): return random.gammavariate(2, 2)
def test_beta(): return random.betavariate(2, 5)
def test_triangular(): return random.triangular(0, 10, 5)
def test_mt(): return random.random()
def test_lcg(): return lcg.random()
def test_llg(): return llg.random()

a = 1.234
b = 5.678
lcg = LCG(seed = 42)
llg = LLG(seed = 42)
times = 100000000

operations = [
    ("Addition (+)", test_add),
    ("Subtraction (-)", test_sub),
    ("Multiplication (*)", test_mul),
    ("Division (/)", test_div),
    ("Modulus (%)", test_mod),
    ("Power (pow)", test_pow),
    ("Absolute (abs)", test_abs),
    ("Logarithm (log)", test_log),
    ("Logarithm Base 10 (log10)", test_log10),
    ("Exponential (exp)", test_exp),
    ("Square Root (sqrt)", test_sqrt),
    ("Square (**)", test_square),
    ("Sine (sin)", test_sin),
    ("Cosine (cos)", test_cos),
    ("Tangent (tan)", test_tan),
    ("Arcsine (asin)", test_asin),
    ("Arccosine (acos)", test_acos),
    ("Arctangent (atan)", test_atan),
    ("Hyperbolic Sine (sinh)", test_sinh),
    ("Hyperbolic Cosine (cosh)", test_cosh),
    ("Hyperbolic Tangent (tanh)", test_tanh),
    ("Inverse Hyperbolic Sine (asinh)", test_asinh),
    ("Inverse Hyperbolic Cosine (acosh)", test_acosh),
    ("Inverse Hyperbolic Tangent (atanh)", test_atanh),
    ("Random Uniform", test_uniform),
    ("Random Normal", test_normal),
    ("Random Integer", test_randint),
    ("Exponential Distribution", test_exponential),
    ("Gamma Distribution", test_gamma),
    ("Beta Distribution", test_beta),
    ("Triangular Distribution", test_triangular),
    ("Mersenne Twister (MT)", test_mt),
    ("Linear Congruential Generator (LCG)", test_lcg),
    ("Linear Lagged Generator (LLG)", test_llg),
]

# output
results = []
for name, func in operations:
    time = timeit.timeit(func, number = times)
    results.append((name, time))
    print(f"{name}: {time:.6f} seconds")

min_time = min([time for name, time in results])
output_data = [(name, time, time / min_time) for name, time in results]
df = pd.DataFrame(output_data, columns = ["Operation", "Time (seconds)", "Relative Time (min = 1)"])
df.to_csv("operation_timing_results.csv", index = False)

print("result:\n", df)

After executing the above program, the execution times of various mathematical operations and pseudo-random number generation methods can be obtained. Since each test was repeated $10^8$ times, the times shown in the table represent the accumulated total execution time when each operation is performed a large number of times, rather than the time required for a single operation.

To facilitate comparison of performance differences among different operations, this article also provides the Relative Time. The multiplication operation (Multiplication, *), which has the lowest execution time, is used as the baseline value of 1. The experimental results are shown in the following table:

OperationExecution Time (seconds)Relative Time
Multiplication (*)7.8796761.000
Subtraction (-)7.9226461.005
Addition (+)8.0031401.016
Division (/)8.0130741.017
Modulus (%)8.8231061.120
Absolute (abs)8.6883031.103
Mersenne Twister (MT)10.5629461.341
Square Root (sqrt)12.4810481.584
Exponential (exp)13.0038021.650
Sine (sin)13.1963231.675
Cosine (cos)13.3902401.699
Logarithm Base 10 (log10)13.6604721.734
Arctangent (atan)13.7103291.740
Hyperbolic Sine (sinh)14.4261541.831
Hyperbolic Cosine (cosh)14.4323881.832
Square (**)14.7487661.872
Arccosine (acos)15.0954591.916
Hyperbolic Tangent (tanh)15.1956171.928
Inverse Hyperbolic Cosine (acosh)14.8847971.889
Inverse Hyperbolic Sine (asinh)15.5598971.975
Tangent (tan)15.9032632.018
Power (pow)16.3763592.078
Inverse Hyperbolic Tangent (atanh)16.7084712.120
Logarithm (log)17.3605672.203
Random Uniform28.0619363.561
Exponential Distribution37.6720434.781
Linear Lagged Generator (LLG)36.6612244.653
Linear Congruential Generator (LCG)40.5223695.143
Triangular Distribution54.2954626.891
Random Normal62.2041527.894
Random Integer72.2156409.165
Gamma Distribution141.51939617.960
Beta Distribution281.85035135.769

From the experimental results, it can be observed that the execution times of the four basic arithmetic operations are very similar. Among them, multiplication has the fastest execution speed, requiring approximately 7.88 seconds to complete one hundred million calculations. Subtraction, addition, and division also take around 8 seconds, with only minor differences between them. This indicates that, in modern CPU architectures, basic arithmetic operations can usually be directly executed by hardware, resulting in very low computational costs.

For more complex mathematical functions, the execution times of square roots, exponentiation, logarithms, and trigonometric functions are significantly higher than those of basic arithmetic operations. The natural logarithm function (log) takes approximately 2.20 times longer than multiplication, while the power function (pow) requires approximately 2.08 times the execution time. For trigonometric functions, both sin() and cos() take around 13 seconds, whereas tan() requires approximately 15.9 seconds.

In addition to basic mathematical functions, this experiment also compares the performance of different random number generation methods. The results show that Python’s built-in Mersenne Twister (MT) requires only approximately 10.56 seconds, which is significantly faster than the LCG and LLG implementations written in Python. Although LCG and Lehmer Generator have very simple computational procedures from an algorithmic perspective, additional overhead from Python class method calls, object state updates, and modulo operations causes their actual execution speed to be lower than the highly optimized built-in random module.

For probability distribution functions, the results show that their execution times are generally higher than simply generating uniformly distributed random numbers. For example, random.uniform() requires approximately 28 seconds, while the Normal Distribution requires approximately 62 seconds. This is because many probability distributions do not directly generate random numbers; instead, they first generate basic uniformly distributed random numbers and then apply additional mathematical transformations to obtain the target distribution. As a result, their execution time increases significantly.

Among all tested probability distributions, the Beta Distribution has the highest execution time, requiring approximately 281.85 seconds, which is 35.77 times slower than the basic multiplication operation. This demonstrates that the computational costs of different mathematical functions can vary substantially. Even though these functions ultimately return only a single numerical value, the underlying algorithms used to generate them may have vastly different levels of complexity.

Environment

  • Operating System: Windows 11 25H2
  • Processor: 13th Gen Intel(R) Core(TM) i7-13700 (2.10 GHz)
  • Memory: 32.0 GB
  • Programming Language: Python 3.10.11

References