Generador de números aleatorios de Lehmer

El generador de números aleatorios de Lehmer ^{[ 1 ]} (llamado así por DH Lehmer ), a veces también conocido como generador de números aleatorios de Park-Miller (por Stephen K.  Park y Keith W. Miller ), es un tipo de generador congruencial lineal (LCG) que opera en el grupo multiplicativo de enteros módulo n . La fórmula general es

${\displaystyle X_{k+1}=a\cdot X_{k}{\bmod {m}},}$

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 ₀ es coprima con m .

Otros nombres son generador congruencial lineal multiplicativo (MLCG) ^{[ 2 ]} y generador congruencial multiplicativo (MCG) .

En 1988, Park y Miller ^{[ 3 ]} _{propusieron un generador de números aleatorios} de Lehmer con parámetros particulares m = ^2³¹ − 1 = 2.147.483.647 (un primo de Mersenne M³¹ ) y a = ^7⁵ = 16.807 (una raíz primitiva módulo M³¹ ), 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 ) . Park, Miller y Stockmeyer respondieron a la crítica (1993), ^{[ 6 ]} diciendo:minstd_rand0

Dada la naturaleza dinámica del área, a quienes no son especialistas les resulta difícil decidir qué generador utilizar. «Necesito algo que pueda comprender, implementar y adaptar... no tiene por qué ser de última generación, solo 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 generador de números aleatorios de Lehmer con parámetros m  =  2 ¹⁶  +  1 = 65 537 (un primo de Fermat F ₄ ) y a  =  75 (un módulo de raíz primitiva F ₄ ). ^{[ 7 ] [ 8 ]} El generador de números aleatorios RANF de CRAY es un generador de números aleatorios de Lehmer con el módulo de potencia de dos m = 2 ⁴⁸ y a = 44 485 709 377 909. ^{[ 9 ]} La Biblioteca Científica GNU incluye varios generadores de números aleatorios de la forma Lehmer, incluidos MINSTD, RANF y el infame generador de números aleatorios RANDU de IBM . ^{[ 9 ]}    

Lo más habitual es elegir el módulo como un número primo, lo que simplifica enormemente la elección de una semilla coprima (cualquier valor entre 0 < X ​​y ₀ < m es válido). Esto produce el resultado de mejor calidad, pero introduce cierta complejidad en la implementación, y es poco probable que el rango del resultado se ajuste a la aplicación deseada; para convertirlo al rango deseado se requiere una multiplicación adicional.

El uso de un módulo m que es una potencia de dos facilita 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. Esto se debe a que los k bits más bajos forman por sí solos un generador ^{de k bits} con módulo 2 ; los bits de orden superior nunca afectan _a los bits de orden inferior. ^{[ 10 ]} Los valores Xᵢ 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 usar los bits más significativos; reducir el rango mediante 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 ₀ debe ser impar.

Es posible usar un módulo compuesto, pero el generador debe inicializarse con un valor coprimo a m , o el período se reducirá considerablemente. Por ejemplo, un módulo de F ₅ = 2 ³²  +  1 podría parecer atractivo, ya que las salidas se pueden mapear fácilmente a una palabra de 32 bits 0 ≤ X _i  −  1 < 2 ³² . Sin embargo, una semilla de X ₀  =  6700417 (que divide a 2 ³²  +  1) o cualquier múltiplo daría como resultado una salida con un período de solo 640.

Otro generador con un módulo compuesto es el recomendado por Nakazawa y Nakazawa: ^{[ 13 ]}

• m =134 265 023 ×134 475 827 =18 055 400 005 099 021 ≈ 2 ⁵⁴
• a =7 759 097 958 782 935 (cualquiera de ± a ^±1 (mod m ) también servirá)

Como ambos factores del módulo son menores que ^2³² , es posible mantener el estado módulo cada uno de los factores y construir el valor de salida utilizando el teorema chino del resto , 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 ; combinar (por ejemplo, sumando sus salidas) 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)( m ₂  −  1)···( m _k  −  1)/2 ^{k −1} . ^{[ 2 ] : 744} Un ejemplo de esto es el generador de Wichmann-Hill .

_{Si bien el generador de números aleatorios de Lehmer (RNG de Lehmer )} puede considerarse un caso particular del generador congruencial lineal con c = 0 , se trata de un caso especial que implica ciertas restricciones y propiedades. En particular, para el RNG de Lehmer, la semilla inicial X₀ debe ser coprima con el módulo m , lo cual no se requiere para los generadores congruenciales lineales (LCG) en general. La elección del módulo m y del multiplicador a también es más restrictiva para el RNG de Lehmer. A diferencia de los LCG, el período máximo del RNG de Lehmer es igual a m − 1, y se cumple cuando m es primo y a es una raíz primitiva módulo m .  

Por otro lado, los logaritmos discretos (en base a o cualquier raíz primitiva módulo m ) de X _k en${\displaystyle \mathbb {Z} _ {m}}$representar una secuencia lineal congruencial módulo el totiente de Euler${\displaystyle \varphi (m)}$.

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 ligeramente menor que una potencia de 2 (los primos de Mersenne 2 ³¹  −  1 y 2 ⁶¹  −  1 son populares, al igual que 2 ³²  −  5 y 2 ⁶⁴  −  59), la reducción módulo m = 2 ^e − d puede implementarse de forma más económica que una división general de doble ancho utilizando la identidad 2 ^e ≡ d (mod m ) .

El paso básico de reducción divide el producto en dos partes de e bits, multiplica la parte de mayor valor por d y las suma: ( ax mod 2e ⁾ + d ⌊ ax /2 ^e ⌋ . A continuación, se resta 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 una si d es pequeño y se elige a < m / d . (Esta condición también garantiza que d ⌊ ax /2 ^e ⌋ sea un producto de ancho simple; si no se cumple, se debe calcular un producto de ancho doble).

Cuando el módulo es un primo de Mersenne ( d  =  1), el procedimiento es particularmente sencillo. No solo la multiplicación por d es trivial, sino que la resta condicional puede reemplazarse por un desplazamiento y una suma incondicionales. Para comprobarlo, observe que el algoritmo garantiza que x ≢ 0 (mod m ) , lo que significa que x  =  0 y x  = m son imposibles. Esto evita la necesidad de considerar representaciones equivalentes de e bits del estado; solo los valores cuyos bits superiores no son cero requieren reducción. 

Los bits e bajos del producto ax no pueden representar un valor mayor que m , y los bits altos nunca tendrá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 activado), 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 ³² − 5 radica en asegurar que produzcamos una única representación para valores como 1 ≡ 2 ³² − 4. La solución consiste en sumar temporalmente d , de modo que el rango de valores posibles sea de d a 2 ^e − 1, y reducir los valores mayores que e bits de manera que nunca se generen representaciones menores que d . Finalmente, restando el desplazamiento temporal se obtiene el valor deseado.          

Comencemos asumiendo que tenemos un valor parcialmente reducido y acotado de modo que 0  ≤ y < 2 m = 2 ^{e +1} − 2 d . En este caso, un solo 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:            

0 ≤ y < m = 2 ^e − d

En este caso, y + d < 2 ^e y ′ =  y < m , como se deseaba .   

m ≤ y < 2 m

En este caso, 2 ^e ≤ y + d < 2 ^{e +1} es un número de ( e + 1) bits, y ⌊ ( y + d )/2 ^e ⌋ = 1. Por lo tanto, y ′ = ( y + d ) − 2 ^e + d − d = y − 2 ^e + d = y − m < m , como se deseaba. Debido a que la parte superior multiplicada es d , la suma es al menos d , y restar el desplazamiento nunca causa subdesbordamiento.                              

(En el caso específico de un generador de Lehmer, un estado cero o su imagen y  = m nunca ocurrirá, por lo que un desplazamiento de d − 1 funcionará igual de bien, si resulta más conveniente. Esto reduce el desplazamiento a 0 en el caso de los números primos de Mersenne, cuando d = 1).     

Reducir un producto mayor ax a menos de 2 m  =  2 ^{e +1}  −  2 d se puede hacer mediante uno o más pasos de reducción sin desplazamiento.

Si ad  ≤ m , entonces un paso de reducción adicional es suficiente. 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 la suposición 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 ^e ,        La primera reducción producirá un valor dentro del 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 un coste:

• El módulo debe poder representarse en un entero con signo ; las operaciones aritméticas deben permitir un rango de ± m .
• La elección del multiplicador a está restringida. Requerimos que m mod a ≤ ⌊ m / a ⌋ , lo cual se logra comúnmente eligiendo a  ≤ √ m . 
• Se requiere una división (con resto) por iteración.

Si bien esta técnica es popular para implementaciones portátiles en lenguajes de alto nivel que carecen de operaciones de doble ancho, ^{[ 2 ] : 744} en computadoras modernas la división por una constante generalmente se implementa 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 factorice m = qa + r , es decir, precalcule las constantes auxiliares r = m mod a y q = ⌊ m / a ⌋ = ( m − r )/ a . Luego, en cada iteración, calcule ax ≡ a ( x mod q ) − r ⌊ x / q ⌋ (mod m ) .

Esta igualdad se mantiene porque

${\displaystyle {\begin{aligned}aq&=a\cdot (mr)/a\\&=(a/a)\cdot (mr)\\&=(mr)\\&\equiv -r{\pmod {m}}\\\end{aligned}}}$

Entonces, si factorizamos x = ( x mod q ) + q ⌊ x / q ⌋ , obtenemos:

${\displaystyle {\begin{aligned}ax&=a(x{\bmod {q}})+aq\lfloor x/q\rfloor \\&=a(x{\bmod {q}})+(mr)\lfloor x/q\rfloor \\&\equiv a(x{\bmod {q}})-r\lfloor x/q\rfloor {\pmod {m}}\\\end{alineado}}}$

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 puede calcularse 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 puede reducirse a [0, m −1] con una sola 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 de a mayores 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 con 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 terminarse.    

Utilizando código C , el generador de números aleatorios 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 puede invocarse repetidamente para generar números pseudoaleatorios, siempre que quien la invoque tenga cuidado de inicializar el estado con 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 enteros de 32 bits podría desbordarse.

Para evitar la división de 64 bits, realice la reducción manualmente:

uint32_t lcg_parkmiller ( uint32_t * estado ) { uint64_t producto = ( uint64_t ) * estado * 48271 ; uint32_t x = ( producto & 0x7fffffff ) + ( producto >> 31 );x = ( x & 0x7fffffff ) + ( x >> 31 ); return * estado = x ; }

Para usar únicamente 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 = 0x7fffffff ; const uint32_t A = 48271 ; const uint32_t Q = M / A ; // 44488 const uint32_t R = M % A ; // 3399uint32_t div = * estado / Q ; // máximo: M / Q = A = 48,271 uint32_t rem = * estado % Q ; // máximo: Q - 1 = 44,487int32_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 result = s - t ;si ( resultado < 0 ) resultado += M ;devolver * estado = resultado ; }

o utilice dos multiplicaciones de 16×16 bits:

uint32_t lcg_parkmiller ( uint32_t * estado ) { const uint32_t A = 48271 ;uint32_t low = ( * state & 0x7fff ) * A ; // máx.: 32.767 * 48.271 = 1.581.695.857 = 0x5e46c371 uint32_t high = ( * state >> 15 ) * A ; // máx.: 65.535 * 48.271 = 3.163.439.985 = 0xbc8e4371 uint32_t x = low + (( high & 0xffff ) << 15 ) + ( high >> 16 ); // máx.: 0x5e46c371 + 0x7fff8000 + 0xbc8e = 0xde46ffffx = ( x & 0x7fffffff ) + ( x >> 31 ); return * estado = x ; }

Otro generador de Lehmer popular utiliza el módulo primo 2 ³² −5:

uint32_t lcg_rand ( uint32_t * estado ) { return * estado = ( uint64_t ) * estado * 279470273u % 0xfffffffb ; }

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 & 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 podría ser de 33 bits. producto = ( producto & 0xffffffff ) + 5 * ( uint32_t )( producto >> 32 );producto += 4 ; // Esta suma está garantizada para ser de 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 ¹²⁸ requiere soporte de 128 bits del compilador y utiliza un multiplicador calculado por L'Ecuyer. ^{[ 19 ]} Tiene un período de 2 ¹²⁶ :

estado estático sin signo __int128 ;/* El estado debe inicializarse 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, así que usamos una expresión const unsigned __int128 mult = ( unsigned __int128 ) 0x12e15e35b500f16e << 64 | 0x2e714eb2b37916a5 ; state *= mult ; return state >> 64 ; }

El generador calcula un valor impar de 128 bits y devuelve sus 64 bits superiores.

Este generador supera la prueba BigCrush de TestU01 , pero falla la prueba TMFn de PractRand . Dicha prueba se diseñó para detectar precisamente el defecto de este tipo de generador: dado que el módulo es una potencia de 2, el período del bit menos significativo en la salida es solo 2⁶² ^, en lugar de ^2¹²⁶ . Los generadores congruenciales lineales con un módulo que es una potencia de 2 presentan un comportamiento similar.

La siguiente rutina principal mejora la velocidad del código anterior para cargas de trabajo con números enteros (si el compilador permite optimizar la declaración de la 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, así 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.

• Lehmer, DH (1949). "Métodos matemáticos en unidades de computación a gran escala". Actas del Segundo Simposio sobre Maquinaria de Cálculo Digital a Gran Escala . págs. 141-146 . MR 0044899 .  (Versión publicada en revista: Annals of the Computation Laboratory of Harvard University , vol. 26 (1951)).
• Steve Park, Generadores de números aleatorios

• Los números primos ligeramente menores que una potencia de dos pueden ser útiles para elegir módulos. Parte de Páginas de números primos .