Last modified: October 02, 2023

This article is written in: 🇵🇱

Liczby losowe i generatory liczb losowych

W języku C++ liczby losowe generuje się za pomocą standardowej biblioteki <random>. Proces losowania zaczyna się od utworzenia generatora liczb pseudolosowych, np. std::mt19937, który bazuje na algorytmie Mersenne Twister. Aby uzyskać bardziej losowe wyniki, generator inicjalizuje się za pomocą unikalnej wartości, zwaną "ziarnem" (ang. seed), co można zrobić np. poprzez std::random_device. Następnie używa się odpowiednich dystrybucji, takich jak std::uniform_int_distribution (dla liczb całkowitych z równomiernym rozkładem) lub std::uniform_real_distribution (dla liczb zmiennoprzecinkowych), aby wygenerować liczby z określonego zakresu. Dzięki tej bibliotece losowanie w C++ jest bardziej elastyczne i daje kontrolę nad różnymi aspektami generowania liczb losowych, w tym nad zakresem i typem wartości.

┌──────────────────────────────────────────┐
│          Generator liczb losowych        │
│  (std::random_device, std::mt19937 itp.) │
└──────────────────────────────────────────┘
                │           ▲
                │ surowe    │ ziarno
                │ bity      │ (opcjonalnie)
                ▼           │
╔═════════════════════════════════════════╗
║   Sekwencja losowych bitów / liczb      ║
╚═════════════════════════════════════════╝
                │
                │ źródło entropii
                ▼
┌─────────────────────────────────────────┐
│             Dystrybucja                 │
│ (std::uniform_int_distribution,         │
│  std::normal_distribution itp.)         │
└─────────────────────────────────────────┘
                │
                │ przekształcenie na
                │ określony rozkład
                ▼
╔═══════════════════════════════════════════════╗
║   Gotowe liczby losowe w zadanym rozkładzie   ║
╚═══════════════════════════════════════════════╝

Liczby losowe

Teraz zapoznamy się z matematycznymi podstawami pojęcia liczby losowej, definiując zmienne losowe oraz ich rozkłady.

W klasycznej teorii prawdopodobieństwa liczbę losową modeluje zmienna losowa

$$ X:(\Omega,\mathcal F,\mathbb P)\longrightarrow (\mathbb R,\mathcal B) $$

gdzie $(\Omega,\mathcal F,\mathbb P)$ to przestrzeń probabilistyczna, a $\mathcal B$ – σ–algebra borelowska.

Dla każdego przedziału $A\subseteq\mathbb R$ prawdopodobieństwo $\mathbb P!\bigl(X\in A\bigr)$ jest ustalone przez rozkład $F_X(x)=\mathbb P(X\le x)$.

Własności oczekiwane

W tej części opisujemy podstawowe miary statystyczne zmiennej losowej, takie jak wartość oczekiwana i wariancja.

Jeśli $X$ ma gęstość $f_X$, to

$$ E[X]=\int_{-\infty}^{\infty}xf_X(x)dx $$

$$ Var[X]=\int_{-\infty}^{\infty}(x- E[X])^{2}f_X(x)dx $$

Prawo wielkich liczb (LLN) gwarantuje zbieżność średniej $\bar X_n$ do $E[X]$ przy $n\to\infty$, a centralne tw. graniczne (CLT) – normalne odchylenie $O(n^{-1/2})$.

Prawdziwa losowość vs. pseudolosowość

Teraz omówimy różnice pomiędzy rzeczywistymi źródłami losowości a deterministycznymi generatorami pseudolosowymi.

Źródła entropii

• Prawdziwe RNG (TRNG). Oparte na zjawiskach fizycznych: szum termiczny, promieniotwórczość, odchylenia czasu taktowania CPU.
• Generator pseudolosowy (PRNG). Deterministyczny algorytm

$$ s_{k+1}=F(s_k)\pmod m,\qquad X_k=g(s_k) $$

który przy zadanym ziarnie $s_0$ tworzy powtarzalną sekwencję. Parametry:

• Okres $p$ – najmniejsze $k>0$ z $s_{n+k}=s_n$.
• Wymiar równomierności – równomierne pokrycie hipersześcianu $[0,1)^d$.
• Test next-bit (dla RNG kryptograficznych) – nieprzewidywalność kolejnego bitu z prawdopodobieństwem istotnie $>\tfrac12$.

Statystyczne testy jakości

Poniższa tabela przedstawia główne klasy testów statystycznych używanych do oceny generatorów pseudolosowych:

Sekwencja testów nie dowodzi losowości, lecz obala ją, gdy statystyki wyjdą poza akceptowalny przedział ufności.

Generowanie liczb losowych

Teraz przedstawimy dostępne w C++ generatory liczb losowych oraz ich główne właściwości i zastosowania.

Przegląd generatorów

Mersenne Twister. Rekurencja:

$$ x_{k+n}=x_{k+m}\oplus\bigl((x_k^{\text{upper}}\|\;x_{k+1}^{\text{lower}})\cdot A\bigr), $$

gdzie $n=624,m=397$. Duży okres i 623-wymiarowa równomierność zapewniają dobry rozkład w praktyce nie-krypto.

Dystrybucje

W tej części omówione zostaną najważniejsze dystrybucje dostępne w bibliotece <random> oraz ich własności.

Dystrybucje ciągłe implementują metodę inwersji

$$ X=F^{-1}(U)\text{ z }U\sim\mathcal U(0,1) $$

lub transformacje specjalne (np. Box-Muller).

Inicjalizacja – właściwości ziarna

Omówimy, jak dobór ziarna wpływa na powtarzalność oraz jakość generowanych sekwencji.

Jeśli dwa generatory otrzymają to samo seed, ich sekwencje będą identyczne:

$$ \forall n\;\;x^{(1)}_n = x^{(2)}_n. $$

W symulacjach replikowalnych podajemy jawnie wartość seeda; w sytuacjach wymagających nieprzewidywalności (kryptografia) używamy std::random_device.

Przykłady użycia (C++ 17)

Poniżej znajdują się przykłady wykorzystania biblioteki <random> w praktycznych scenariuszach, ilustrujące różne zastosowania dystrybucji.

Liczba z przedziału $[a,b]\subset\mathbb Z$

Wygenerujemy liczbę całkowitą z zadanego przedziału.

std::mt19937 gen(std::random_device{}());
std::uniform_int_distribution<int> dist(a,b);
int x = dist(gen);

Własności. $\Pr{X=k}=1/(b-a+1)$. Oczekiwana wartość $E[X]=\tfrac{a+b}{2}$, wariancja $\sigma^2=\tfrac{(b-a+1)^2-1}{12}$.

Rzut monetą – dystrybucja Bernoulliego

Przykład rzutu uczciwą monetą.

std::bernoulli_distribution coin(0.5);
bool isHeads = coin(gen);

Analiza. $X\sim\mathrm{Bern}(p)$: $E[X]=p$, Var[X]=p(1-p). Dla $n$ rzutów błąd względny częstości maleje jak $O!\bigl(n^{-1/2}\bigr)$.

Rzut sześciościenną kostką

Generowanie wyrzutu kostki o sześciu ściankach.

std::uniform_int_distribution<int> d6(1,6);
int result = d6(gen);

Rozkład dyskretny równomierny na ${1,\dots,6}$. $E[X]=3.5,; \sigma^2=\tfrac{35}{12}$.

Generator silnych haseł

Przykład tworzenia losowego hasła o zadanej długości.

std::string alphabet =
    "abcdefghijklmnopqrstuvwxyz"
    "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
    "0123456789!@#$%^&*()-_=+";
std::uniform_int_distribution<size_t> pick(0, alphabet.size()-1);

auto password = [&] (size_t n) {
    std::string s;
    for(size_t i=0;i<n;++i) s+=alphabet[pick(gen)];
    return s;
};

Entropia hasła $H = n\log_2 |\text{alphabet}|$. Dla $n=12$ i $|alphabet| = 74$: $H\approx 71$ bitów.

Metody Monte Carlo

Demonstracja estymacji całki przy pomocy losowania.

std::uniform_real_distribution<double> U(a,b);
double sum = 0;
for(size_t k=0;k<N;++k) sum += g(U(gen));
double I = (b-a)*sum/N;

Symulacyjna estymata całki:

$$ I=\int_{a}^{b}g(x)dx, \qquad
\hat I_N = (b-a)\frac1N\sum_{k=1}^N g(U_k),\quad U_k\sim\mathcal U(a,b). $$

Błąd przeciętny $RMSE=\sigma/\sqrt N$ z $\sigma^2=Var[g(U)]$. Uwaga. Dla funkcji silnie oscylujących warto stosować ważoną próbę (importance sampling) lub stratyfikację.

RNG kryptograficzne (CSPRNG)

Wymagania formalne (Goldwasser–Micali):

1. Jednokierunkowość – brak efektywnego algorytmu do odtworzenia seeda.
2. Test next-bit –

$$ \Pr[G(x)i = 1 \mid G(x)_1, G(x)_2, \dots, G(x){i-1}] = \frac{1}{2} \pm \varepsilon(\lambda). $$

gdzie $\varepsilon$ jest zaniedbywalne względem parametru bezpieczeństwa $\lambda$.

W C++20 brak oficjalnego CSPRNG, ale popularne biblioteki (libsodium, Botan) udostępniają randombytes_buf, crypto_rng. Dla systemowego źródła entropii można użyć:

#include <openssl/rand.h>
unsigned char buf[32];
RAND_bytes(buf, sizeof(buf)); // 256-bit

Zastosowanie różnych dystrybucji

Teraz omówimy zastosowania kilku dystrybucji losowych dostępnych w bibliotece <random>, prezentując ich definicje, własności, implementację w C++ oraz przykłady użycia.

Dystrybucja równomierna

Dystrybucja równomierna generuje wartości o jednakowym prawdopodobieństwie w całym zakresie, co jest przydatne, gdy potrzebujemy równomiernie rozłożyć próbki.

Dyskretna wersja na zbiorze kolejnych liczb całkowitych ${a,\dots,b}$:

$$ \mathbb P(X=k)=\frac1{b-a+1},\qquad k=a,\dots,b $$

Wersja ciągła na przedziale $[a,b]$:

$$ f(x)=\frac1{b-a}\mathbf 1_{[a,b]}(x) $$

$$ F(x)=\frac{x-a}{b-a} $$

W poniższych wzorach przedstawiono momenty i funkcję generującą momenty (MGF), które charakteryzują tę dystrybucję:

$$ E[X]=\frac{a+b}{2} $$

$$ Var[X]=\frac{(b-a)^2}{12} $$

$$ M_X(t)=\frac{e^{tb}-e^{ta}}{t(b-a)}\;(t\ne 0) $$

Implementacja

std::mt19937 gen(std::random_device{}());
std::uniform_int_distribution<int> U(a,b);
int k = U(gen);

Rada. Generator wywołuj raz i przekazuj referencję; kosztowna inicjalizacja std::random_device{} nie powinna znajdować się w pętli.

Przykład: 10 liczb w zakresie $[1,100]$

Poniżej przykład wylosowania dziesięciu liczb z przedziału [1,100]:

std::cout << "Uniform[1,100]: ";
for(int i=0;i<10;++i) std::cout << U(gen) << ' ';

Dystrybucja normalna (Gaussa)

Rozkład normalny, zwany też Gaussa, opisuje zmienną losową o gęstości dzwonowej i jest uzywany w statystyce ze względu na centralne twierdzenie graniczne.

$$ f(x)=\frac1{\sqrt{2\pi}\sigma}\exp!\Bigl[-\frac{(x-\mu)^2}{2\sigma^2}\Bigr], \quad x\in\mathbb R $$

Parametry

Podstawowe parametry wpływające na położenie i rozrzut rozkładu przedstawiono poniżej:

$$ E[X]=\mu $$

$$ Var[X]=\sigma^{2} $$

$$ M_X(t)=\exp!\left(\mu t+\tfrac12\sigma^{2}t^{2}\right) $$

Suma niezależnych $X_i\sim\mathcal N(\mu_i,\sigma_i^2)$ jest również normalna:

$$ \sum_i X_i \sim \mathcal N!\Bigl(\sum \mu_i,\;\sum\sigma_i^2\Bigr) $$

Implementacja

std::normal_distribution<double> N(mu, sigma);
double x = N(gen);

Przykład generowania próbek z rozkładu normalnego o zadanych parametrach:

std::cout << "Normal(50,10): ";
for(int i=0;i<10;++i) std::cout << N(gen) << ' ';

Weryfikacja. Aby sprawdzić poprawność generatora, warto przeprowadzić prostą weryfikację statystyczną:

1. zbuduj histogram,
2. oblicz statystykę K–S $D_n=\sup_x|F_n(x)-\Phi_{\mu,\sigma}(x)|$,
3. odrzucaj generator jeśli $D_n>D_{\alpha,n}$ dla wybranego poziomu $\alpha$.

Dystrybucja Bernoulliego

Dystrybucja Bernoulliego modeluje zdarzenie dychotomiczne, zwracając 1 z prawdopodobieństwem $p$ lub 0 z prawdopodobieństwem $1-p$.

$$ \mathbb P(X=1)=p,\qquad \mathbb P(X=0)=1-p,\quad 0<p<1. $$

Parametry

$$ E[X]=p $$

$$ Var[X]=p(1-p) $$

$$ M_X(t)=1-p+pe^{t} $$

Związek z rozkladem dwumianowym. Suma prób Bernoulliego prowadzi do rozkładu dwumianowego, co ilustruje związek pomiędzy tymi rozkładami:

$$ \sum_{i=1}^{n}X_i \sim \mathrm{Bin}(n,p). $$

Implementacja

double p = 0.3;
std::bernoulli_distribution B(p);
bool success = B(gen);

Przykład: 10 losowań, $p=0,3$

Przykład wykonania dziesięciu niezależnych prób Bernoulliego:

for(int i=0;i<10;++i) std::cout << B(gen) << ' ';

Analiza częstości. Analiza częstościowa pozwala na ocenę zbieżności częstości względnej do $p$ oraz wyznaczenie przedziału ufności 95 %:

$\hat p \pm 1.96\sqrt{p(1-p)/n}$

Tabela poniżej zestawia informacje o omawianych dystrybucjach, umożliwiając szybkie porównanie:

Zalety i wady różnych metod

Teraz porównamy najpopularniejsze metody generowania liczb losowych w C++, omawiając ich zalety oraz ograniczenia w różnych zastosowaniach.

std::random (C++11)

Biblioteka <random> w C++11 dostarcza elastyczny interfejs „generator + dystrybucja” bazujący domyślnie na std::mt19937. Poniżej przyjrzymy się matematycznemu modelowi Mersenne Twistera wraz z jego mocnymi i słabymi stronami.

Model matematyczny

Domyślny wybór – std::mt19937 – jest implementacją algorytmu Mersenne Twistera:

$$ \begin{aligned} x_{k+n} &= x_{k+m}\oplus\Bigl((x_k^{u}\parallel x_{k+1}^{\ell})A\Bigr),\\ n&=624,\;m=397,\\ \text{okres}&=2^{19937}-1\approx 10^{6001},\\ d_{\max}&=623\quad\text{(wymiar równomierności).} \end{aligned} $$

Spektralny test równoległości daje wynik rzędu $2^{-2018}$ – przy wielu zastosowaniach symulacyjnych uznaje się to za „praktycznie idealne”.

Zalety

Wady

1. Trzeba rozumieć dwustopniowy model: generator + dystrybucja.
2. Brak odporności kryptograficznej – Mersenne Twister jest liniowy nad \$\mathbb F_2\$; stan generatora (19937 bitów) można odtworzyć na podstawie 19937 kolejnych wyjść.
3. Koszt inicjalizacji – std::random_device może być wolny (blokuje podczas odczytu z /dev/urandom).
4. Rozmiar stanu – 2,5 KB dla mt19937; przy wielu wątkach oznacza to znaczną ilość pamięci.

rand() – klasyczny LCG

Funkcja rand() to prosty liniowy kongruentny generator (LCG) o minimalnym API, często używany dla szybkiego startu, ale z istotnymi ograniczeniami.

Model matematyczny

Większość implementacji to liniowy generator:

$$ x_{n+1} = (ax_n + c)\;\bmod m, $$

typowe stałe (glibc, MSVC):

$$ a=1103515245,\;c=12345,\;m=2^{31}. $$

Okres ≤ $m=2^{31};(\approx2.1\cdot10^9)$. Najmłodszy bit ma okres 2, co czyni go bezużytecznym w symulacjach binarnych.

Zalety

Wady

1. Krótki okres generacji – przy prędkości 10⁸ liczb/s pełny cykl wyczerpuje się w ok. 21 s.
2. Słabe własności spektralne – pojawiają się korelacje w przestrzeni 3D w symulacjach Monte Carlo.
3. Niejednorodność implementacji – różne kompilatory używają odmiennych stałych.
4. Brak wbudowanych rozkładów – konieczne ręczne przekształcenia (np. metoda Box–Mullera).
5. Łamliwość kryptograficzna – na podstawie zaledwie trzech wyjść można odzyskać parametry generatora.

Zestawienie porównawcze

Poniżej tabela podsumowująca różnice między std::mt19937 a rand():