El generador de números aleatorios Lehmer [1] (llamado así por DH Lehmer ), a veces también conocido como generador de números aleatorios Park-Miller (en honor a Stephen K. Park y Keith W. Miller), es un tipo de generador congruencial lineal (LCG) que opera en un grupo multiplicativo de números enteros módulo n . La fórmula general es
donde el módulo m es un número primo o una potencia de un número primo , el multiplicador a es un elemento de alto orden multiplicativo módulo m (por ejemplo, una raíz primitiva módulo n ), y la semilla X 0 es coprimo a m .
Otros nombres son generador congruencial lineal multiplicativo (MLCG) [2] y generador congruencial multiplicativo (MCG) .
En 1988, Park y Miller [3] sugirieron un generador de números aleatorios de Lehmer con parámetros particulares m = 2 31 − 1 = 2.147.483.647 (un primo de Mersenne M 31 ) y a = 7 5 = 16.807 (una raíz primitiva módulo M 31 ), ahora conocido como MINSTD . Aunque MINSTD fue criticado posteriormente por Marsaglia y Sullivan (1993), [4] [5] todavía se utiliza hoy en día (en particular, en CarbonLib y C++11 's minstd_rand0). Park, Miller y Stockmeyer respondieron a la crítica (1993), [6] diciendo:
Dada la naturaleza dinámica del área, es difícil para los no especialistas tomar decisiones sobre qué generador utilizar. "Denme algo que pueda entender, implementar y trasladar... no tiene que ser de última generación, sólo asegúrense de que sea razonablemente bueno y eficiente". Nuestro artículo y el generador estándar mínimo asociado fueron un intento de responder a esta solicitud. Cinco años después, no vemos la necesidad de modificar nuestra respuesta, salvo sugerir el uso del multiplicador a = 48271 en lugar de 16807.
Esta constante revisada se utiliza en el generador de números aleatorios de C++11minstd_rand .
El Sinclair ZX81 y sus sucesores utilizan el RNG de Lehmer con parámetros m = 2 16 + 1 = 65,537 (un primo de Fermat F 4 ) y a = 75 (una raíz primitiva módulo F 4 ). [7] [8] El generador de números aleatorios CRAY RANF es un RNG de Lehmer con el módulo de potencia de dos m = 2 48 y a = 44,485,709,377,909. [9] La Biblioteca Científica GNU incluye varios generadores de números aleatorios de la forma Lehmer, incluyendo MINSTD, RANF, y el infame generador de números aleatorios de IBM RANDU . [9]
Lo más común es que el módulo se elija como un número primo, lo que hace que la elección de una semilla coprima sea trivial (cualquier 0 < X 0 < m servirá). Esto produce el resultado de mejor calidad, pero introduce cierta complejidad de implementación y es poco probable que el rango del resultado coincida con la aplicación deseada; la conversión al rango deseado requiere una multiplicación adicional.
El uso de un módulo m que es una potencia de dos permite una implementación informática particularmente conveniente, pero tiene un costo: el período es como máximo m /4, y los bits inferiores tienen períodos más cortos que eso. Esto se debe a que los k bits más bajos forman un generador k módulo 2 por sí mismos; los bits de orden superior nunca afectan a los bits de orden inferior. [10] Los valores Xi son siempre impares (el bit 0 nunca cambia), los bits 2 y 1 se alternan ( los 3 bits inferiores se repiten con un período de 2), los 4 bits inferiores se repiten con un período de 4, y así sucesivamente. Por lo tanto, la aplicación que utiliza estos números aleatorios debe utilizar los bits más significativos; reducir a un rango más pequeño utilizando una operación de módulo con un módulo par producirá resultados desastrosos. [11]
Para lograr este período, el multiplicador debe satisfacer a ≡ ±3 (mod 8), [12] y la semilla X 0 debe ser impar.
Es posible utilizar un módulo compuesto, pero el generador debe tener como semilla un valor coprimo con m , o el período se reducirá considerablemente. Por ejemplo, un módulo de F 5 = 2 32 + 1 puede parecer atractivo, ya que las salidas se pueden mapear fácilmente a una palabra de 32 bits 0 ≤ X i − 1 < 2 32 . Sin embargo, una semilla de X 0 = 6700417 (que divide 2 32 + 1) o cualquier múltiplo conduciría a una salida con un período de solo 640.
Otro generador con módulo compuesto es el recomendado por Nakazawa & Nakazawa: [13]
Como ambos factores del módulo son menores que 2 32 , es posible mantener el estado módulo de cada uno de los factores y construir el valor de salida utilizando el teorema del resto chino , utilizando no más de 64 bits de aritmética intermedia. [13] : 70
Una implementación más popular para períodos grandes es un generador congruencial lineal combinado ; la combinación (por ejemplo, sumando sus salidas) de varios generadores es equivalente a la salida de un solo generador cuyo módulo es el producto de los módulos de los generadores componentes. [14] y cuyo período es el mínimo común múltiplo de los períodos componentes. Aunque los períodos compartirán un divisor común de 2, los módulos se pueden elegir de modo que sea el único divisor común y el período resultante sea ( m 1 − 1)( m 2 − 1)···( m k − 1)/2 k −1 . [2] : 744 Un ejemplo de esto es el generador de Wichmann-Hill .
Si bien el generador aleatorio de Lehmer puede considerarse un caso particular del generador congruencial lineal con c = 0 , es un caso especial que implica ciertas restricciones y propiedades. En particular, para el generador aleatorio de Lehmer, la semilla inicial X 0 debe ser coprima con el módulo m , lo que no se requiere para los LCG en general. La elección del módulo m y del multiplicador a también es más restrictiva para el generador aleatorio de Lehmer. A diferencia del LCG, el período máximo del generador aleatorio de Lehmer es igual a m − 1, y es así cuando m es primo y a es una raíz primitiva módulo m .
Por otra parte, los logaritmos discretos (en base a o cualquier raíz primitiva módulo m ) de X k en representan una secuencia congruencial lineal módulo el totiente de Euler .
Un módulo primo requiere el cálculo de un producto de doble ancho y un paso de reducción explícito. Si se utiliza un módulo apenas menor que una potencia de 2 (los primos de Mersenne 2 31 − 1 y 2 61 − 1 son populares, al igual que 2 32 − 5 y 2 64 − 59), la reducción módulo m = 2 e − d se puede implementar de manera más económica que una división de doble ancho general utilizando la identidad 2 e ≡ d (mod m ) .
El paso básico de reducción divide el producto en dos partes de e -bit, multiplica la parte alta por d y las suma: ( ax mod 2 e ) + d ⌊ ax /2 e ⌋ . A esto se le puede seguir restando m hasta que el resultado esté dentro del rango. El número de restas está limitado a ad / m , que se puede limitar fácilmente a uno si d es pequeño y se elige a < m / d . (Esta condición también asegura que d ⌊ ax /2 e ⌋ sea un producto de ancho simple; si se viola, se debe calcular un producto de ancho doble).
Cuando el módulo es un primo de Mersenne ( d = 1), el procedimiento es particularmente simple. No sólo es trivial la multiplicación por d , sino que la resta condicional puede ser reemplazada por un desplazamiento y una adición incondicionales. Para ver esto, note que el algoritmo garantiza que x ≢ 0 (mod m ) , lo que significa que x = 0 y x = m son ambos imposibles. Esto evita la necesidad de considerar representaciones equivalentes de e -bit del estado; sólo los valores donde los bits altos no son cero necesitan reducción.
Los bits e bajos del producto ax no pueden representar un valor mayor que m , y los bits altos nunca contendrán un valor mayor que a − 1 ≤ m − 2. Por lo tanto, el primer paso de reducción produce un valor como máximo m + a − 1 ≤ 2 m − 2 = 2 e +1 − 4. Este es un número de ( e + 1) bits, que puede ser mayor que m (es decir, podría tener el bit e establecido), pero la mitad alta es como máximo 1, y si lo es, los bits e bajos serán estrictamente menores que m . Por lo tanto, ya sea que el bit alto sea 1 o 0, un segundo paso de reducción (suma de las mitades) nunca desbordará los bits e , y la suma será el valor deseado.
Si d > 1, también se puede evitar la resta condicional, pero el procedimiento es más complejo. El desafío fundamental de un módulo como 2 32 − 5 radica en garantizar que produzcamos solo una representación para valores como 1 ≡ 2 32 − 4. La solución es agregar temporalmente d , de modo que el rango de valores posibles sea d hasta 2 e − 1, y reducir los valores mayores que e bits de una manera que nunca genere representaciones menores que d . Finalmente, restar el desplazamiento temporal produce el valor deseado.
Comencemos suponiendo que tenemos un valor parcialmente reducido y acotado de modo que 0 ≤ y < 2 m = 2 e +1 − 2 d . En este caso, un único paso de resta de desplazamiento producirá 0 ≤ y ′ = (( y + d ) mod 2 e ) + d ⌊ ( y + d )/2 e ⌋ − d < m . Para ver esto, consideremos dos casos:
(Para el caso de un generador de Lehmer específicamente, un estado cero o su imagen y = m nunca ocurrirá, por lo que un desplazamiento de d − 1 funcionará de la misma manera, si eso es más conveniente. Esto reduce el desplazamiento a 0 en el caso primo de Mersenne, cuando d = 1).
La reducción de un producto mayor ax a menos de 2 m = 2 e +1 − 2 d se puede realizar mediante uno o más pasos de reducción sin desplazamiento.
Si ad ≤ m , entonces basta con un paso de reducción adicional. Dado que x < m , ax < am ≤ ( a − 1)2 e , y un paso de reducción convierte esto en, como máximo, 2 e − 1 + ( a − 1) d = m + ad − 1. Esto está dentro del límite de 2 m si ad − 1 < m , que es el supuesto inicial.
Si ad > m , entonces es posible que el primer paso de reducción produzca una suma mayor que 2 m = 2 e +1 − 2 d , que es demasiado grande para el paso de reducción final. (También requiere la multiplicación por d para producir un producto mayor que e bits, como se mencionó anteriormente). Sin embargo, siempre que d 2 < 2 e , la primera reducción producirá un valor en el rango requerido para que se aplique el caso anterior de dos pasos de reducción.
Si no se dispone de un producto de doble ancho, se puede utilizar el método de Schrage , [15] [16] también llamado método de factorización aproximada, [17] para calcular ax mod m , pero esto tiene el costo:
Si bien esta técnica es popular para implementaciones portables en lenguajes de alto nivel que carecen de operaciones de doble ancho, [2] : 744 en computadoras modernas la división por una constante se implementa usualmente usando multiplicación de doble ancho, por lo que esta técnica debe evitarse si la eficiencia es una preocupación. Incluso en lenguajes de alto nivel, si el multiplicador a está limitado a √ m , entonces el producto de doble ancho ax puede calcularse usando dos multiplicaciones de ancho simple y reducirse usando las técnicas descritas anteriormente.
Para utilizar el método de Schrage, primero se factoriza m = qa + r , es decir, se precalculan las constantes auxiliares r = m mod a y q = ⌊ m / a ⌋ = ( m − r )/ a . Luego, en cada iteración, se calcula ax ≡ a ( x mod q ) − r ⌊ x / q ⌋ (mod m ) .
Esta igualdad se cumple porque
Entonces, si factorizamos x = ( x mod q ) + q ⌊ x / q ⌋ , obtenemos:
La razón por la que no se desborda es que ambos términos son menores que m . Dado que x mod q < q ≤ m / a , el primer término es estrictamente menor que am / a = m y se puede calcular con un producto de ancho simple.
Si se elige a de modo que r ≤ q (y por lo tanto r / q ≤ 1), entonces el segundo término también es menor que m : r ⌊ x / q ⌋ ≤ rx / q = x ( r / q ) ≤ x (1) = x < m . Por lo tanto, la diferencia se encuentra en el rango [1− m , m −1] y se puede reducir a [0, m −1] con una única suma condicional. [18]
Esta técnica puede extenderse para permitir un r negativo (− q ≤ r < 0), cambiando la reducción final a una resta condicional.
La técnica también puede extenderse para permitir valores mayores de a aplicándola recursivamente. [17] : 102 De los dos términos restados para producir el resultado final, solo el segundo ( r ⌊ x / q ⌋ ) corre el riesgo de desbordarse. Pero esto es en sí mismo una multiplicación modular por una constante de tiempo de compilación r , y puede implementarse mediante la misma técnica. Debido a que cada paso, en promedio, reduce a la mitad el tamaño del multiplicador (0 ≤ r < a , valor promedio ( a −1)/2), esto parecería requerir un paso por bit y ser espectacularmente ineficiente. Sin embargo, cada paso también divide x por un cociente cada vez mayor q = ⌊ m / a ⌋ , y rápidamente se llega a un punto donde el argumento es 0 y la recursión puede terminar.
Usando código C , el RNG de Park-Miller se puede escribir de la siguiente manera:
uint32_t lcg_parkmiller ( uint32_t * estado ) { return * estado = ( uint64_t ) * estado * 48271 % 0x7fffffff ; } Esta función se puede llamar repetidamente para generar números pseudoaleatorios, siempre que el llamador tenga cuidado de inicializar el estado en cualquier número mayor que cero y menor que el módulo. En esta implementación, se requiere aritmética de 64 bits; de lo contrario, el producto de dos números enteros de 32 bits puede desbordarse.
Para evitar la división de 64 bits, haga la reducción a mano:
uint32_t lcg_parkmiller ( uint32_t * estado ) { uint64_t producto = ( uint64_t ) * estado * 48271 ; uint32_t x = ( producto & 0x7ffffffff ) + ( producto >> 31 ); x = ( x & 0x7ffffffff ) + ( x >> 31 ); devolver * estado = x ; } Para utilizar sólo aritmética de 32 bits, utilice el método de Schrage:
uint32_t lcg_parkmiller ( uint32_t * state ) { // Parámetros precalculados para el método de Schrage const uint32_t M = 0x7ffffffff ; const uint32_t A = 48271 ; const uint32_t Q = M / A ; // 44488 const uint32_t R = M % A ; // 3399 uint32_t div = * estado / Q ; // máx: M / Q = A = 48,271 uint32_t rem = * estado % Q ; // máx: Q - 1 = 44,487 int32_t s = rem * A ; // máx.: 44 487 * 48 271 = 2 147 431 977 = 0x7fff3629 int32_t t = div * R ; // máx.: 48 271 * 3 399 = 164 073 129 int32_t resultado = s - t ; si ( resultado < 0 ) resultado += M ; devolver * estado = resultado ; } o utilizar dos multiplicaciones de 16×16 bits:
uint32_t lcg_parkmiller ( uint32_t * estado ) { const uint32_t A = 48271 ; uint32_t bajo = ( * estado y 0x7fff ) * A ; // máx.: 32 767 * 48 271 = 1 581 695 857 = 0x5e46c371 uint32_t alto = ( * estado >> 15 ) * A ; // máx.: 65 535 * 48 271 = 3 163 439 985 = 0xbc8e4371 uint32_t x = bajo + (( alto y 0xffff ) << 15 ) + ( alto >> 16 ); // máx.: 0x5e46c371 + 0x7fff8000 + 0xbc8e = 0xde46ffff x = ( x & 0x7ffffffff ) + ( x >> 31 ); devolver * estado = x ; } Otro generador de Lehmer popular utiliza el módulo primo 2 32 −5:
uint32_t lcg_rand ( uint32_t * estado ) { return * estado = ( uint64_t ) * estado * 279470273u % 0xffffffffb ; } Esto también se puede escribir sin una división de 64 bits:
uint32_t lcg_rand ( uint32_t * estado ) { uint64_t producto = ( uint64_t ) * estado * 279470273u ; uint32_t x ; // No es necesario porque 5 * 279470273 = 0x5349e3c5 cabe en 32 bits. // producto = (producto y 0xffffffff) + 5 * (producto >> 32); // Un multiplicador mayor que 0x33333333 = 858.993.459 lo necesitaría.// El resultado de la multiplicación cabe en 32 bits, pero la suma puede ser de 33 bits. producto = ( producto & 0xffffffff ) + 5 * ( uint32_t )( producto >> 32 ); producto += 4 ; // Se garantiza que esta suma tendrá 32 bits. x = ( uint32_t ) producto + 5 * ( uint32_t )( producto >> 32 ); return * estado = x - 4 ; } Muchos otros generadores de Lehmer tienen buenas propiedades. El siguiente generador de Lehmer módulo 2 128 requiere soporte de 128 bits por parte del compilador y utiliza un multiplicador calculado por L'Ecuyer. [19] Tiene un período de 2 126 :
estado estático sin signo __int128 ; /* El estado debe ser inicializado con un valor impar. */ void seed ( unsigned __int128 seed ) { state = seed << 1 | 1 ; } uint64_t next ( void ) { // GCC no puede escribir literales de 128 bits, por lo que usamos una expresión const unsigned __int128 mult = ( unsigned __int128 ) 0x12e15e35b500f16e << 64 | 0x2e714eb2b37916a5 ; estado *= mult ; devolver estado >> 64 ; } El generador calcula un valor impar de 128 bits y devuelve sus 64 bits superiores.
Este generador pasa la prueba BigCrush de TestU01 , pero no pasa la prueba TMFn de PractRand. Esa prueba se diseñó para detectar exactamente el defecto de este tipo de generador: dado que el módulo es una potencia de 2, el período del bit más bajo en la salida es solo 2 62 , en lugar de 2 126. Los generadores congruenciales lineales con un módulo de potencia de 2 tienen un comportamiento similar .
La siguiente rutina principal mejora la velocidad del código anterior para cargas de trabajo de números enteros (si el compilador permite optimizar la declaración de constante fuera de un bucle de cálculo):
uint64_t next ( void ) { uint64_t result = state >> 64 ; // GCC no puede escribir literales de 128 bits, por lo que usamos una expresión const unsigned __int128 mult = ( unsigned __int128 ) 0x12e15e35b500f16e << 64 | 0x2e714eb2b37916a5 ; state *= mult ; return result ; } Sin embargo, debido a que la multiplicación se difiere, no es adecuada para el hash, ya que la primera llamada simplemente devuelve los 64 bits superiores del estado de la semilla.
El ZX81 utiliza p=65537 y a=75 [...](Tenga en cuenta que el manual del ZX81 afirma incorrectamente que 65537 es un primo de Mersenne que equivale a 2 · 16 − 1. El manual del ZX Spectrum corrigió esto y afirma correctamente que es un primo de Fermat que equivale a 2 · 16 + 1).
El ZX Spectrum utiliza p=65537 y a=75, y almacena algunos bi-1 en la memoria.
El dado se determina mediante aritmética modular, p. ej.,
, ... La función CRAY RANF solo arroja tres de los seis resultados posibles (¡cuyos tres lados dependen de la semilla)!
lrand48() % 6 + 1mz1y mz2calcular cada iteración mz1a, mz2aes más eficiente mantener estas últimas como variables de estado. Además, el módulo de reducción final m (llamado iden el libro) es de un valor menor que 2 m , por lo que puede consistir en una única resta condicional.