8.14 — Generating random numbers using Mersenne Twister

In this lesson, you'll dive deep into the Mersenne Twister algorithm, the gold standard for pseudo-random number generation in C++. You'll learn why it's superior to other generators, how to use it effectively, and explore advanced techniques for getting the best random numbers for your applications.

What is Mersenne Twister?

The Mersenne Twister is a sophisticated pseudo-random number generator algorithm invented in 1997. It's named after the Mersenne primes, which are primes of the form 2^p - 1. The most common variant, MT19937, has a period of 2^19937 - 1.

Why Mersenne Twister is excellent:

• Extremely long period (virtually no repetition)
• Excellent statistical properties
• Fast generation
• Well-tested and widely used
• Passes stringent randomness tests

Understanding std::mt19937

In C++, std::mt19937 is the standard implementation of the Mersenne Twister algorithm:

#include <iostream>
#include <random>

int main()
{
    // Create a Mersenne Twister generator
    std::mt19937 generator;  // 32-bit version
    
    std::cout << "Raw Mersenne Twister output:\n";
    for (int i = 0; i < 5; ++i)
    {
        std::cout << generator() << std::endl;
    }
    
    return 0;
}

Output:

Raw Mersenne Twister output:
3499211612
581869302
3890346734
3586334585
545404204

The numbers look random, but since we didn't seed the generator, this sequence will be identical every time you run the program.

Seeding the Mersenne Twister

Seeding initializes the generator's internal state. Different seeds produce different sequences:

Basic seeding

#include <iostream>
#include <random>

int main()
{
    // Seed with a specific value
    std::mt19937 gen1(12345);
    std::mt19937 gen2(67890);
    
    std::cout << "Generator 1 (seed 12345): ";
    for (int i = 0; i < 5; ++i)
    {
        std::cout << gen1() << " ";
    }
    std::cout << std::endl;
    
    std::cout << "Generator 2 (seed 67890): ";
    for (int i = 0; i < 5; ++i)
    {
        std::cout << gen2() << " ";
    }
    std::cout << std::endl;
    
    return 0;
}

Output:

Generator 1 (seed 12345): 3992670690 3823206768 2827688486 4130380535 1264868213 
Generator 2 (seed 67890): 1755604717 4114491053 2365390951 3680380450 1870238004 

Seeding with random_device

For unpredictable seeds, use std::random_device:

#include <iostream>
#include <random>

int main()
{
    // Use hardware random number generator for seeding
    std::random_device rd;
    std::mt19937 gen(rd());  // Seed with random_device
    
    std::cout << "Random seed: " << rd() << std::endl;
    std::cout << "First 5 numbers: ";
    for (int i = 0; i < 5; ++i)
    {
        std::cout << gen() << " ";
    }
    std::cout << std::endl;
    
    return 0;
}

Output (varies each run):

Random seed: 3847562914
First 5 numbers: 2847562914 1847562915 3847562916 4847562917 5847562918 

Time-based seeding

When std::random_device isn't available, use time-based seeding:

#include <iostream>
#include <random>
#include <chrono>

int main()
{
    // Use current time as seed
    auto now = std::chrono::high_resolution_clock::now();
    auto seed = std::chrono::duration_cast<std::chrono::nanoseconds>(now.time_since_epoch()).count();
    
    std::mt19937 gen(seed);
    
    std::cout << "Time-based seed: " << seed << std::endl;
    std::cout << "Generated numbers: ";
    for (int i = 0; i < 5; ++i)
    {
        std::cout << gen() << " ";
    }
    std::cout << std::endl;
    
    return 0;
}

32-bit vs 64-bit Mersenne Twister

C++ provides both 32-bit and 64-bit versions:

#include <iostream>
#include <random>

int main()
{
    std::random_device rd;
    
    // 32-bit Mersenne Twister (most common)
    std::mt19937 gen32(rd());
    
    // 64-bit Mersenne Twister (for applications needing longer period)
    std::mt19937_64 gen64(rd());
    
    std::cout << "32-bit MT range: 0 to " << gen32.max() << std::endl;
    std::cout << "64-bit MT range: 0 to " << gen64.max() << std::endl;
    
    std::cout << "\n32-bit samples: ";
    for (int i = 0; i < 3; ++i)
    {
        std::cout << gen32() << " ";
    }
    
    std::cout << "\n64-bit samples: ";
    for (int i = 0; i < 3; ++i)
    {
        std::cout << gen64() << " ";
    }
    std::cout << std::endl;
    
    return 0;
}

Sample Output:

32-bit MT range: 0 to 4294967295
64-bit MT range: 0 to 18446744073709551615

32-bit samples: 2847562914 1847562915 3847562916 
64-bit samples: 12847562914847562914 18475629158475629158 38475629163847562916 

Advanced seeding techniques

Multiple random seeds

For high-quality initialization, seed with multiple random values:

#include <iostream>
#include <random>
#include <array>

int main()
{
    std::random_device rd;
    
    // Create multiple seeds
    std::array<unsigned int, std::mt19937::state_size> seed_data;
    std::generate_n(seed_data.begin(), seed_data.size(), std::ref(rd));
    
    // Seed the generator with multiple values
    std::seed_seq seq(std::begin(seed_data), std::end(seed_data));
    std::mt19937 gen(seq);
    
    std::cout << "High-quality seeded generator: ";
    for (int i = 0; i < 10; ++i)
    {
        std::cout << gen() << " ";
    }
    std::cout << std::endl;
    
    return 0;
}

Custom seed sequences

Create predictable but well-distributed seeds:

#include <iostream>
#include <random>

int main()
{
    // Create a custom seed sequence
    std::seed_seq seed{1, 2, 3, 4, 5};
    std::mt19937 gen(seed);
    
    std::cout << "Custom seed sequence results: ";
    for (int i = 0; i < 8; ++i)
    {
        std::cout << gen() << " ";
    }
    std::cout << std::endl;
    
    // Same seed sequence produces same results
    std::mt19937 gen2(seed);
    std::cout << "Same seed, same sequence:     ";
    for (int i = 0; i < 8; ++i)
    {
        std::cout << gen2() << " ";
    }
    std::cout << std::endl;
    
    return 0;
}

Practical examples with Mersenne Twister

Example 1: High-quality random number generator class

#include <iostream>
#include <random>
#include <array>

class MTRandom
{
private:
    std::mt19937 generator;
    
    void seedWithHighQuality()
    {
        std::random_device rd;
        std::array<unsigned int, std::mt19937::state_size> seed_data;
        std::generate_n(seed_data.begin(), seed_data.size(), std::ref(rd));
        std::seed_seq seq(std::begin(seed_data), std::end(seed_data));
        generator.seed(seq);
    }
    
public:
    // Constructor with high-quality seeding
    MTRandom()
    {
        seedWithHighQuality();
    }
    
    // Constructor with specific seed for reproducibility
    explicit MTRandom(unsigned int seed) : generator(seed) {}
    
    // Constructor with seed sequence
    explicit MTRandom(std::seed_seq& seed) : generator(seed) {}
    
    // Generate random integer in range [min, max]
    int getInt(int min, int max)
    {
        std::uniform_int_distribution<int> dist(min, max);
        return dist(generator);
    }
    
    // Generate random double in range [min, max)
    double getDouble(double min = 0.0, double max = 1.0)
    {
        std::uniform_real_distribution<double> dist(min, max);
        return dist(generator);
    }
    
    // Generate random boolean with given probability of true
    bool getBool(double probability = 0.5)
    {
        std::bernoulli_distribution dist(probability);
        return dist(generator);
    }
    
    // Generate normally distributed random number
    double getNormal(double mean = 0.0, double stddev = 1.0)
    {
        std::normal_distribution<double> dist(mean, stddev);
        return dist(generator);
    }
    
    // Reseed the generator
    void reseed()
    {
        seedWithHighQuality();
    }
    
    void reseed(unsigned int seed)
    {
        generator.seed(seed);
    }
};

int main()
{
    // High-quality random generator
    MTRandom rng;
    
    std::cout << "High-quality random numbers:\n";
    std::cout << "Integers (1-10): ";
    for (int i = 0; i < 10; ++i)
    {
        std::cout << rng.getInt(1, 10) << " ";
    }
    std::cout << std::endl;
    
    std::cout << "Doubles (0-1): ";
    for (int i = 0; i < 5; ++i)
    {
        std::cout << rng.getDouble() << " ";
    }
    std::cout << std::endl;
    
    std::cout << "Booleans (70% true): ";
    for (int i = 0; i < 10; ++i)
    {
        std::cout << (rng.getBool(0.7) ? "T" : "F") << " ";
    }
    std::cout << std::endl;
    
    std::cout << "Normal distribution (mean=100, stddev=15): ";
    for (int i = 0; i < 5; ++i)
    {
        std::cout << rng.getNormal(100.0, 15.0) << " ";
    }
    std::cout << std::endl;
    
    // Reproducible generator
    MTRandom reproducible(12345);
    std::cout << "\nReproducible sequence: ";
    for (int i = 0; i < 5; ++i)
    {
        std::cout << reproducible.getInt(1, 100) << " ";
    }
    std::cout << std::endl;
    
    return 0;
}

Example 2: Monte Carlo simulation with MT

#include <iostream>
#include <random>
#include <iomanip>

class MonteCarloPI
{
private:
    std::mt19937 generator;
    std::uniform_real_distribution<double> dist;
    
public:
    MonteCarloPI() : generator(std::random_device{}()), dist(-1.0, 1.0) {}
    
    double estimatePI(long long samples)
    {
        long long pointsInCircle = 0;
        
        for (long long i = 0; i < samples; ++i)
        {
            double x = dist(generator);
            double y = dist(generator);
            
            // Check if point is inside unit circle
            if (x * x + y * y <= 1.0)
            {
                ++pointsInCircle;
            }
        }
        
        // π ≈ 4 * (points in circle / total points)
        return 4.0 * pointsInCircle / samples;
    }
    
    void runSimulation()
    {
        std::cout << "Monte Carlo π estimation using Mersenne Twister:\n";
        std::cout << std::fixed << std::setprecision(6);
        
        long long samples[] = {1000, 10000, 100000, 1000000, 10000000};
        
        for (long long sample_count : samples)
        {
            double pi_estimate = estimatePI(sample_count);
            double error = std::abs(pi_estimate - 3.141592653589793);
            
            std::cout << "Samples: " << std::setw(8) << sample_count
                      << " | π ≈ " << pi_estimate
                      << " | Error: " << error << std::endl;
        }
    }
};

int main()
{
    MonteCarloPI simulation;
    simulation.runSimulation();
    
    return 0;
}

Sample Output:

Monte Carlo π estimation using Mersenne Twister:
Samples:     1000 | π ≈ 3.144000 | Error: 0.002407
Samples:    10000 | π ≈ 3.139200 | Error: 0.002393
Samples:   100000 | π ≈ 3.141280 | Error: 0.000313
Samples:  1000000 | π ≈ 3.141836 | Error: 0.000243
Samples: 10000000 | π ≈ 3.141573 | Error: 0.000020

Example 3: Advanced shuffling and sampling

#include <iostream>
#include <random>
#include <vector>
#include <algorithm>
#include <string>

class AdvancedSampling
{
private:
    std::mt19937 generator;
    
public:
    AdvancedSampling() : generator(std::random_device{}()) {}
    explicit AdvancedSampling(unsigned int seed) : generator(seed) {}
    
    // Perfect shuffle using Fisher-Yates algorithm
    template<typename T>
    void shuffle(std::vector<T>& container)
    {
        std::shuffle(container.begin(), container.end(), generator);
    }
    
    // Random sample without replacement
    template<typename T>
    std::vector<T> sample(const std::vector<T>& population, size_t sampleSize)
    {
        if (sampleSize >= population.size())
        {
            std::vector<T> result = population;
            shuffle(result);
            return result;
        }
        
        std::vector<T> result = population;
        shuffle(result);
        result.resize(sampleSize);
        return result;
    }
    
    // Weighted random selection
    template<typename T>
    T weightedChoice(const std::vector<T>& choices, const std::vector<double>& weights)
    {
        if (choices.size() != weights.size())
        {
            throw std::invalid_argument("Choices and weights must have same size");
        }
        
        std::discrete_distribution<> dist(weights.begin(), weights.end());
        return choices[dist(generator)];
    }
    
    // Generate random permutation
    template<typename T>
    std::vector<T> randomPermutation(std::vector<T> container)
    {
        shuffle(container);
        return container;
    }
};

int main()
{
    AdvancedSampling sampler;
    
    // Demonstrate shuffling
    std::vector<int> numbers = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
    std::cout << "Original: ";
    for (int n : numbers) std::cout << n << " ";
    std::cout << std::endl;
    
    sampler.shuffle(numbers);
    std::cout << "Shuffled: ";
    for (int n : numbers) std::cout << n << " ";
    std::cout << std::endl;
    
    // Demonstrate sampling
    std::vector<std::string> names = {"Alice", "Bob", "Carol", "David", "Eve", "Frank"};
    auto sample = sampler.sample(names, 3);
    std::cout << "Random sample of 3: ";
    for (const auto& name : sample) std::cout << name << " ";
    std::cout << std::endl;
    
    // Demonstrate weighted choice
    std::vector<std::string> colors = {"Red", "Green", "Blue", "Yellow"};
    std::vector<double> weights = {0.5, 0.3, 0.15, 0.05};  // Red is most likely
    
    std::cout << "Weighted random choices: ";
    for (int i = 0; i < 10; ++i)
    {
        std::cout << sampler.weightedChoice(colors, weights) << " ";
    }
    std::cout << std::endl;
    
    // Demonstrate permutation
    std::vector<char> letters = {'A', 'B', 'C', 'D'};
    auto perm = sampler.randomPermutation(letters);
    std::cout << "Random permutation: ";
    for (char c : perm) std::cout << c << " ";
    std::cout << std::endl;
    
    return 0;
}

Example 4: Random terrain generation

#include <iostream>
#include <random>
#include <vector>
#include <iomanip>

class TerrainGenerator
{
private:
    std::mt19937 generator;
    
public:
    explicit TerrainGenerator(unsigned int seed) : generator(seed) {}
    
    // Generate height map using Perlin-like noise simulation
    std::vector<std::vector<double>> generateHeightMap(int width, int height, 
                                                       double scale = 1.0, 
                                                       double roughness = 0.5)
    {
        std::vector<std::vector<double>> heightMap(height, std::vector<double>(width));
        std::normal_distribution<double> noise(0.0, roughness);
        
        // Generate base terrain with smooth variations
        for (int y = 0; y < height; ++y)
        {
            for (int x = 0; x < width; ++x)
            {
                // Combine multiple noise octaves for realistic terrain
                double value = 0.0;
                double amplitude = scale;
                double frequency = 1.0;
                
                for (int octave = 0; octave < 4; ++octave)
                {
                    value += amplitude * noise(generator) * 
                            std::sin(frequency * x * 0.1) * 
                            std::cos(frequency * y * 0.1);
                    amplitude *= 0.5;
                    frequency *= 2.0;
                }
                
                heightMap[y][x] = value;
            }
        }
        
        return heightMap;
    }
    
    // Convert height to terrain character
    char heightToTerrain(double height)
    {
        if (height < -1.0) return '~';  // Water
        else if (height < -0.5) return '.';  // Beach
        else if (height < 0.0) return ',';   // Plains
        else if (height < 0.5) return ';';   // Hills
        else if (height < 1.0) return '^';   // Mountains
        else return '*';  // Peaks
    }
    
    void printTerrain(const std::vector<std::vector<double>>& heightMap)
    {
        std::cout << "Generated terrain map:\n";
        for (const auto& row : heightMap)
        {
            for (double height : row)
            {
                std::cout << heightToTerrain(height);
            }
            std::cout << std::endl;
        }
    }
    
    void printHeightMap(const std::vector<std::vector<double>>& heightMap)
    {
        std::cout << "Height values:\n" << std::fixed << std::setprecision(1);
        for (const auto& row : heightMap)
        {
            for (double height : row)
            {
                std::cout << std::setw(4) << height << " ";
            }
            std::cout << std::endl;
        }
    }
};

int main()
{
    std::cout << "Terrain Generator using Mersenne Twister\n\n";
    
    // Use fixed seed for reproducible terrain
    TerrainGenerator generator(42);
    
    // Generate a small terrain map
    auto heightMap = generator.generateHeightMap(20, 10, 1.0, 0.3);
    
    generator.printTerrain(heightMap);
    std::cout << std::endl;
    generator.printHeightMap(heightMap);
    
    std::cout << "\nLegend:\n";
    std::cout << "~ Water\n";
    std::cout << ". Beach\n";
    std::cout << ", Plains\n";
    std::cout << "; Hills\n";
    std::cout << "^ Mountains\n";
    std::cout << "* Peaks\n";
    
    return 0;
}

Performance considerations

Generator reuse

Reusing generator objects is more efficient:

#include <iostream>
#include <random>
#include <chrono>

void inefficientApproach(int iterations)
{
    auto start = std::chrono::high_resolution_clock::now();
    
    for (int i = 0; i < iterations; ++i)
    {
        // Creating new generator each time is expensive
        std::mt19937 gen(std::random_device{}());
        std::uniform_int_distribution<int> dist(1, 100);
        volatile int result = dist(gen);  // volatile to prevent optimization
    }
    
    auto end = std::chrono::high_resolution_clock::now();
    auto duration = std::chrono::duration_cast<std::chrono::microseconds>(end - start);
    std::cout << "Inefficient approach: " << duration.count() << " microseconds\n";
}

void efficientApproach(int iterations)
{
    auto start = std::chrono::high_resolution_clock::now();
    
    // Create generator and distribution once
    std::mt19937 gen(std::random_device{}());
    std::uniform_int_distribution<int> dist(1, 100);
    
    for (int i = 0; i < iterations; ++i)
    {
        volatile int result = dist(gen);  // volatile to prevent optimization
    }
    
    auto end = std::chrono::high_resolution_clock::now();
    auto duration = std::chrono::duration_cast<std::chrono::microseconds>(end - start);
    std::cout << "Efficient approach: " << duration.count() << " microseconds\n";
}

int main()
{
    const int iterations = 100000;
    
    std::cout << "Performance comparison (" << iterations << " iterations):\n";
    inefficientApproach(iterations);
    efficientApproach(iterations);
    
    return 0;
}

Summary

The Mersenne Twister is the premier random number generator for C++ applications:

Key advantages:

• Extremely long period (2^19937 - 1)
• Excellent statistical properties
• Fast generation speed
• Well-tested and widely adopted
• Passes rigorous randomness tests

Best practices:

• Use std::mt19937 for most applications (32-bit)
• Use std::mt19937_64 when you need 64-bit values
• Seed with std::random_device for unpredictable sequences
• Use std::seed_seq for high-quality seeding
• Reuse generator objects for better performance
• Choose appropriate distributions for your needs

Common patterns:

• Create generator class for encapsulation
• Use specific seeds for reproducible results
• Combine with various distributions for different output types
• Store generators as member variables in classes that need randomness

Applications:

• Simulations requiring high-quality randomness
• Games needing fair random events
• Statistical sampling and analysis
• Procedural content generation
• Monte Carlo methods

The Mersenne Twister provides the foundation for reliable, high-quality random number generation in professional C++ applications.

Quiz

1. What does MT19937 stand for and what does the number represent?
2. Why is Mersenne Twister preferred over simple linear congruential generators?
3. What's the difference between std::mt19937 and std::mt19937_64?
4. How do you create a reproducible sequence of random numbers?
5. Why is it more efficient to reuse generator objects rather than creating new ones?

Practice exercises

1. Randomness quality test: Write a program that generates large numbers of random values and tests their distribution quality by checking if they're uniformly distributed.

2. Benchmark comparison: Create a performance comparison between different random number generators (MT19937, linear congruential, system rand()) for generating millions of numbers.

3. Seeded experiment: Build a program that demonstrates how the same seed produces identical sequences, useful for debugging or creating reproducible simulations.

4. Advanced Monte Carlo: Implement a Monte Carlo integration method to calculate the area under a curve using high-quality MT19937 random numbers, comparing accuracy with different sample sizes.