¿Por qué la nueva biblioteca aleatoria es mejor que std :: rand ()?

Así que vi una charla llamada rand () Considered Nocivo y abogó por el uso del paradigma de distribución del motor de generación de números aleatorios sobre el std::rand()paradigma simple más módulo.

Sin embargo, quería ver las fallas de std::rand()primera mano, así que hice un experimento rápido:

1. Básicamente, escribí 2 funciones getRandNum_Old()y getRandNum_New()eso generó un número aleatorio entre 0 y 5 inclusive usando std::rand()y std::mt19937+ std::uniform_int_distributionrespectivamente.
2. Luego generé 960.000 (divisibles por 6) números aleatorios usando la forma "antigua" y registré las frecuencias de los números 0-5. Luego calculé la desviación estándar de estas frecuencias. Lo que estoy buscando es una desviación estándar lo más baja posible, ya que eso es lo que sucedería si la distribución fuera realmente uniforme.
3. Ejecuté esa simulación 1000 veces y registré la desviación estándar para cada simulación. También registré el tiempo que tomó en milisegundos.
4. Después, hice exactamente lo mismo de nuevo, pero esta vez generando números aleatorios de la "nueva" forma.
5. Finalmente, calculé la desviación estándar y media de la lista de desviaciones estándar tanto para la forma antigua como para la nueva y la desviación media y estándar para la lista de tiempos tomados para la forma antigua y nueva.

Estos fueron los resultados:

[OLD WAY]
Spread
       mean:  346.554406
    std dev:  110.318361
Time Taken (ms)
       mean:  6.662910
    std dev:  0.366301

[NEW WAY]
Spread
       mean:  350.346792
    std dev:  110.449190
Time Taken (ms)
       mean:  28.053907
    std dev:  0.654964

Sorprendentemente, la extensión total de rollos fue la misma para ambos métodos. Es decir, std::mt19937+ std::uniform_int_distributionno era "más uniforme" que simple std::rand()+ %. Otra observación que hice fue que la nueva forma era 4 veces más lenta que la antigua. En general, parecía que estaba pagando un gran costo en velocidad por casi ninguna ganancia en calidad.

¿Mi experimento tiene algún defecto? ¿O std::rand()realmente no es tan malo, y tal vez incluso mejor?

Como referencia, aquí está el código que usé en su totalidad:

#include <cstdio>
#include <random>
#include <algorithm>
#include <chrono>

int getRandNum_Old() {
    static bool init = false;
    if (!init) {
        std::srand(time(nullptr)); // Seed std::rand
        init = true;
    }

    return std::rand() % 6;
}

int getRandNum_New() {
    static bool init = false;
    static std::random_device rd;
    static std::mt19937 eng;
    static std::uniform_int_distribution<int> dist(0,5);
    if (!init) {
        eng.seed(rd()); // Seed random engine
        init = true;
    }

    return dist(eng);
}

template <typename T>
double mean(T* data, int n) {
    double m = 0;
    std::for_each(data, data+n, [&](T x){ m += x; });
    m /= n;
    return m;
}

template <typename T>
double stdDev(T* data, int n) {
    double m = mean(data, n);
    double sd = 0.0;
    std::for_each(data, data+n, [&](T x){ sd += ((x-m) * (x-m)); });
    sd /= n;
    sd = sqrt(sd);
    return sd;
}

int main() {
    const int N = 960000; // Number of trials
    const int M = 1000;   // Number of simulations
    const int D = 6;      // Num sides on die

    /* Do the things the "old" way (blech) */

    int freqList_Old[D];
    double stdDevList_Old[M];
    double timeTakenList_Old[M];

    for (int j = 0; j < M; j++) {
        auto start = std::chrono::high_resolution_clock::now();
        std::fill_n(freqList_Old, D, 0);
        for (int i = 0; i < N; i++) {
            int roll = getRandNum_Old();
            freqList_Old[roll] += 1;
        }
        stdDevList_Old[j] = stdDev(freqList_Old, D);
        auto end = std::chrono::high_resolution_clock::now();
        auto dur = std::chrono::duration_cast<std::chrono::microseconds>(end-start);
        double timeTaken = dur.count() / 1000.0;
        timeTakenList_Old[j] = timeTaken;
    }

    /* Do the things the cool new way! */

    int freqList_New[D];
    double stdDevList_New[M];
    double timeTakenList_New[M];

    for (int j = 0; j < M; j++) {
        auto start = std::chrono::high_resolution_clock::now();
        std::fill_n(freqList_New, D, 0);
        for (int i = 0; i < N; i++) {
            int roll = getRandNum_New();
            freqList_New[roll] += 1;
        }
        stdDevList_New[j] = stdDev(freqList_New, D);
        auto end = std::chrono::high_resolution_clock::now();
        auto dur = std::chrono::duration_cast<std::chrono::microseconds>(end-start);
        double timeTaken = dur.count() / 1000.0;
        timeTakenList_New[j] = timeTaken;
    }

    /* Display Results */

    printf("[OLD WAY]\n");
    printf("Spread\n");
    printf("       mean:  %.6f\n", mean(stdDevList_Old, M));
    printf("    std dev:  %.6f\n", stdDev(stdDevList_Old, M));
    printf("Time Taken (ms)\n");
    printf("       mean:  %.6f\n", mean(timeTakenList_Old, M));
    printf("    std dev:  %.6f\n", stdDev(timeTakenList_Old, M));
    printf("\n");
    printf("[NEW WAY]\n");
    printf("Spread\n");
    printf("       mean:  %.6f\n", mean(stdDevList_New, M));
    printf("    std dev:  %.6f\n", stdDev(stdDevList_New, M));
    printf("Time Taken (ms)\n");
    printf("       mean:  %.6f\n", mean(timeTakenList_New, M));
    printf("    std dev:  %.6f\n", stdDev(timeTakenList_New, M));
}

Prácticamente cualquier implementación de "antiguo" rand()usa un LCG ; Si bien generalmente no son los mejores generadores que existen, generalmente no los verá fallar en una prueba tan básica; la media y la desviación estándar generalmente son correctas incluso con los peores PRNG.

Los errores comunes de las rand()implementaciones "malas", pero lo suficientemente comunes, son:

• baja aleatoriedad de bits de bajo orden;
• período corto;
• baja RAND_MAX;
• cierta correlación entre extracciones sucesivas (en general, los LCG producen números que se encuentran en un número limitado de hiperplanos, aunque esto puede mitigarse de alguna manera).

Aún así, ninguno de estos es específico de la API de rand(). Una implementación particular podría colocar un generador de la familia xorshift detrás de srand/ randy, hablando de manera algorítmica, obtener un PRNG de última generación sin cambios de interfaz, por lo que ninguna prueba como la que hizo mostraría alguna debilidad en la salida.

Editar: @R. señala correctamente que la interfaz rand/ srandestá limitada por el hecho de que srandtoma un unsigned int, por lo que cualquier generador que una implementación pueda poner detrás de ellos está intrínsecamente limitado a UINT_MAXposibles semillas iniciales (y por lo tanto secuencias generadas). De hecho, esto es cierto, aunque la API podría extenderse trivialmente para hacer que srandtome una sobrecarga unsigned long longo agregar una srand(unsigned char *, size_t).

De hecho, el problema real con rand()no es mucho de implementación en principio sino:

• compatibilidad al revés; muchas implementaciones actuales utilizan generadores subóptimos, normalmente con parámetros mal elegidos; un ejemplo notorio es Visual C ++, que tiene un valor RAND_MAXde solo 32767. Sin embargo, esto no se puede cambiar fácilmente, ya que rompería la compatibilidad con el pasado; las personas que usan sranduna semilla fija para simulaciones reproducibles no estarían muy contentas (de hecho, IIRC la implementación mencionada se remonta a las primeras versiones de Microsoft C, o incluso a Lattice C, de mediados de los ochenta);

• interfaz simplista; rand()proporciona un solo generador con el estado global de todo el programa. Si bien esto está perfectamente bien (y en realidad es bastante útil) para muchos casos de uso simples, plantea problemas:

  • con código multiproceso: para solucionarlo, necesita un mutex global, que ralentizaría todo sin ningún motivo y eliminaría cualquier posibilidad de repetibilidad, ya que la secuencia de llamadas se vuelve aleatoria en sí misma, o un estado local de subproceso; este último ha sido adoptado por varias implementaciones (notablemente Visual C ++);
  • si desea una secuencia reproducible "privada" en un módulo específico de su programa que no afecte el estado global.

Finalmente, el randestado de cosas:

• no especifica una implementación real (el estándar C proporciona solo una implementación de muestra), por lo que cualquier programa que esté destinado a producir una salida reproducible (o esperar un PRNG de alguna calidad conocida) en diferentes compiladores debe generar su propio generador;
• no proporciona ningún método multiplataforma para obtener una semilla decente ( time(NULL)no lo es, ya que no es lo suficientemente granular y, a menudo, piense en dispositivos integrados sin RTC, ni siquiera lo suficientemente aleatorio).

De ahí el nuevo <random>encabezado, que intenta arreglar este lío proporcionando algoritmos que son:

• completamente especificado (para que pueda tener una salida reproducible del compilador cruzado y características garantizadas, por ejemplo, rango del generador);
• generalmente de calidad de vanguardia ( desde el momento en que se diseñó la biblioteca ; ver más abajo);
• encapsulado en clases (por lo que no se le impone ningún estado global, lo que evita completamente los problemas de enhebrado y no localidad);

... y un valor predeterminado random_devicetambién para sembrarlos.

Ahora, si me preguntas, me hubiera gustado también una API simple construida encima de esto para los casos "fáciles", "adivina un número" (similar a cómo Python proporciona la API "complicada", pero también la trivial random.randint& Co .usando un PRNG global pre-sembrado para nosotros, personas sencillas a las que no les gustaría ahogarse en dispositivos / motores / adaptadores aleatorios / lo que sea cada vez que queramos extraer un número para los cartones de bingo), pero es cierto que puede hacerlo fácilmente. constrúyalo usted mismo sobre las instalaciones actuales (aunque no sería posible crear la API "completa" sobre una simplista).

Finalmente, para volver a la comparación de rendimiento: como han especificado otros, está comparando un LCG rápido con un Mersenne Twister más lento (pero generalmente considerado de mejor calidad); si está de acuerdo con la calidad de un LCG, puede usar en std::minstd_randlugar de std::mt19937.

De hecho, después de ajustar su función para usar std::minstd_randy evitar variables estáticas inútiles para la inicialización

int getRandNum_New() {
    static std::minstd_rand eng{std::random_device{}()};
    static std::uniform_int_distribution<int> dist{0, 5};
    return dist(eng);
}

Obtengo 9 ms (antiguo) frente a 21 ms (nuevo); finalmente, si me deshago de dist(que, en comparación con el operador de módulo clásico, maneja el sesgo de distribución para el rango de salida no múltiplo del rango de entrada) y vuelvo a lo que está haciendo engetRandNum_Old()

int getRandNum_New() {
    static std::minstd_rand eng{std::random_device{}()};
    return eng() % 6;
}

Lo reduzco a 6 ms (es decir, un 30% más rápido), probablemente porque, a diferencia de la llamada a rand(), std::minstd_randes más fácil de integrar.

Por cierto, hice la misma prueba usando un enrollado a mano (pero bastante conforme a la interfaz de biblioteca estándar) XorShift64*, y es 2,3 veces más rápido que rand()(3,68 ms frente a 8,61 ms); dado que, a diferencia del Mersenne Twister y los diversos LCG proporcionados, pasa los conjuntos de pruebas de aleatoriedad actuales con gran éxito y es increíblemente rápido, hace que te preguntes por qué no está incluido en la biblioteca estándar todavía.

Si repite su experimento con un rango superior a 5, probablemente verá resultados diferentes. Cuando su rango es significativamente más pequeño, no RAND_MAXhay problema para la mayoría de las aplicaciones.

Por ejemplo, si tenemos un RAND_MAXde 25 rand() % 5, producirá números con las siguientes frecuencias:

0: 6
1: 5
2: 5
3: 5
4: 5

Como RAND_MAXse garantiza que es más de 32767 y la diferencia de frecuencias entre lo menos probable y lo más probable es solo 1, para números pequeños la distribución es lo suficientemente aleatoria para la mayoría de los casos de uso.

Primero, sorprendentemente, la respuesta cambia dependiendo de para qué esté usando el número aleatorio. Si va a conducir, digamos, un cambiador de color de fondo aleatorio, usar rand () está perfectamente bien. Si está utilizando un número aleatorio para crear una mano de póquer aleatoria o una clave segura criptográficamente, no está bien.

Previsibilidad: la secuencia 012345012345012345012345 ... proporcionaría una distribución uniforme de cada número en su muestra, pero obviamente no es aleatorio. Para que una secuencia sea aleatoria, el valor de n + 1 no se puede predecir fácilmente por el valor de n (o incluso por los valores de n, n-1, n-2, n-3, etc.) Claramente una secuencia repetida de los mismos dígitos es un caso degenerado, pero una secuencia generada con cualquier generador congruencial lineal puede someterse a análisis; Si usa la configuración predeterminada lista para usar de un LCG común de una biblioteca común, una persona malintencionada podría "romper la secuencia" sin mucho esfuerzo. En el pasado, varios casinos en línea (y algunos tradicionales) sufrieron pérdidas por máquinas que utilizaban generadores de números aleatorios deficientes. Incluso las personas que deberían saberlo mejor se han puesto al día;

Distribución: como se menciona en el video, tomar un módulo de 100 (o cualquier valor que no se pueda dividir uniformemente en la longitud de la secuencia) garantizará que algunos resultados serán al menos un poco más probables que otros resultados. En el universo de 32767 posibles valores iniciales módulo 100, los números del 0 al 66 aparecerán 328/327 (0,3%) con más frecuencia que los valores del 67 al 99; un factor que puede proporcionar una ventaja al atacante.

La respuesta correcta es: depende de lo que quieras decir con "mejor".

Los "nuevos" <random>motores se introdujeron en C ++ hace más de 13 años, por lo que no son realmente nuevos. La biblioteca C rand()se introdujo hace décadas y ha sido muy útil en ese momento para muchas cosas.

La biblioteca estándar de C ++ proporciona tres clases de motores generadores de números aleatorios: Linear Congruential (del cual rand()es un ejemplo), Lagged Fibonacci y Mersenne Twister. Hay compensaciones de cada clase, y cada clase es "mejor" en ciertos aspectos. Por ejemplo, los LCG tienen un estado muy pequeño y, si se eligen los parámetros correctos, son bastante rápidos en los procesadores de escritorio modernos. Los LFG tienen un estado más grande y solo usan recuperaciones de memoria y operaciones de adición, por lo que son muy rápidos en sistemas integrados y microcontroladores que carecen de hardware matemático especializado. El MTG tiene un estado enorme y es lento, pero puede tener una secuencia no repetitiva muy grande con excelentes características espectrales.

Si ninguno de los generadores suministrados es lo suficientemente bueno para su uso específico, la biblioteca estándar de C ++ también proporciona una interfaz para un generador de hardware o su propio motor personalizado. Ninguno de los generadores está destinado a ser utilizado de forma independiente: su uso previsto es a través de un objeto de distribución que proporciona una secuencia aleatoria con una función de distribución de probabilidad particular.

Otra ventaja de <random>over rand()es que rand()usa el estado global, no es reentrante ni seguro para subprocesos, y permite una sola instancia por proceso. Si necesita un control detallado o previsibilidad (es decir, capaz de reproducir un error dado el estado de semilla RNG), entonces rand()es inútil. Los <random>generadores se instancian localmente y tienen un estado serializable (y restaurable).