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 (después de Stephen K. Park y Keith W. Miller), es un tipo de generador congruencial lineal (LCG) que opera en grupo multiplicativo de 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 una es un elemento de alta multiplicativo orden modulo m (por ejemplo, una primitiva raíz módulo n ), y la semilla X 0 es primos entre sí a m .
Otros nombres son generador congruencial lineal multiplicativo (MLCG) [2] y generador congruencial multiplicativo (MCG) .
Parámetros de uso común
En 1988, Park y Miller [3] sugirieron un RNG 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 (un módulo de raíz primitivo M 31 ), ahora conocido como MINSTD . Aunque MINSTD más tarde fue criticada por Marsaglia y Sullivan (1993), [4] [5] que todavía está en uso hoy en día (en particular, en CarbonLib y C ++ 11 's minstd_rand0
). Park, Miller y Stockmeyer respondieron a las críticas (1993), [6] diciendo:
Dada la naturaleza dinámica del área, es difícil para los no especialistas tomar decisiones sobre qué generador utilizar. "Dame algo que pueda entender, implementar y adaptar ... no es necesario que sea de vanguardia, solo asegúrate 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 más que sugerir el uso del multiplicador a = 48271 en lugar de 16807.
Esta constante revisado se utiliza en C ++ 11 's minstd_rand
generador de números aleatorios.
El Sinclair ZX81 y sus sucesores utilizan el Lehmer RNG con los parámetros m = 2 16 +1 = 65,537 (un Fermat primo F 4 ) y a = 75 (un módulo raíz primitivo F 4 ). [7] 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. [8] 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 de IBM RANDU . [8]
Elección de módulo
Más comúnmente, el módulo se elige 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 en la 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 lo convierte en una implementación de computadora particularmente conveniente, pero tiene un costo: el período es como máximo m / 4, y los bits bajos tienen períodos más cortos que eso. Esto se debe a que los k bits bajos forman un generador de módulo 2 k por sí mismos; los bits de orden superior nunca afectan a los bits de orden inferior. [9] Los valores X i son siempre impares (el bit 0 nunca cambia), los bits 2 y 1 se alternan (los 3 bits bajos se repiten con un período de 2), los 4 bits bajos se repiten con un período de 4, y así sucesivamente. Por lo tanto, la aplicación que usa estos números aleatorios debe usar los bits más significativos; reducir a un rango más pequeño usando una operación de módulo con un módulo uniforme producirá resultados desastrosos. [10]
Para lograr este período, el multiplicador debe satisfacer un ≡ ± 3 (mod 8) [11] y la semilla X 0 debe ser impar.
El uso de un módulo compuesto es posible, pero el generador debe ser sembrado con un coprimero valor a 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 asignar 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.
Una implementación más popular para períodos extensos 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 de componentes. [12] 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 2 ··· ( m k - 1) / 2 k −1 . [2] : 744 Un ejemplo de esto es el generador de Wichmann – Hill .
Relación con LCG
Si bien el Lehmer RNG puede verse como un caso particular del generador congruencial lineal con c = 0 , es un caso especial que implica ciertas restricciones y propiedades. En particular, para el Lehmer RNG, la semilla inicial X 0 debe ser coprime al módulo m que no se requiere para los LCG en general. La elección del módulo my el multiplicador a también es más restrictiva para el Lehmer RNG. En contraste con LCG, el período máximo del Lehmer RNG es igual a m −1 y es tal cuando m es primo y a es una raíz primitiva módulo m .
Por otro lado, los logaritmos discretos (para basar a o cualquier raíz primitiva módulo m ) de X k enrepresentar una secuencia congruencial lineal módulo el coeficiente de Euler .
Implementación
Un módulo primo requiere el cálculo de un producto de doble ancho y un paso de reducción explícito. Si se usa un módulo un poco 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), el módulo de reducción 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 de reducción básico divide el producto en dos partes e -bit, multiplica la parte alta por d y las suma: ( ax mod 2 e ) + d ⌊ ax / 2 e ⌋ . A esto le puede seguir restando m hasta que el resultado esté dentro del rango. El número de sustracciones se limita a ad / m , que puede ser fácilmente limitado a uno si d es pequeño y un < m / d se elige. (Esta condición también asegura que d ⌊ ax / 2 e ⌋ es 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 solo la multiplicación por d es trivial, sino que la resta condicional se puede reemplazar por un cambio y una suma incondicionales. Para ver esto, tenga en cuenta 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 -bit del estado; sólo los valores en los que los bits altos son distintos de cero necesitan reducción.
Los bajos e bits del producto hacha no pueden representar un valor mayor que m , y los bits altos nunca se llevará a cabo una mayor valor que un -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 bit 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 producimos solo una representación para valores como 1 ≡ 2 32 −4. La solución es sumar temporalmente d para que el rango de valores posibles sea d a 2 e −1, y reducir valores mayores que e bits de una manera que nunca genere representaciones menores que d . Finalmente, restando el desplazamiento temporal se obtiene el valor deseado.
Comience asumiendo que tenemos un valor y parcialmente reducido que está 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, considere dos casos:
- 0 ≤ y < m = 2 e - d
- En este caso, y + d <2 e y y ′ = y < m , como se desee.
- m ≤ y <2 m
- En este caso, 2 e ≤ y + d <2 e +1 es un número de e + 1 bit, 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 desee. Debido a que la parte alta multiplicada es d , la suma es al menos dy restar el desplazamiento nunca causa subdesbordamiento.
(Para el caso de un generador de Lehmer específicamente, nunca se producirá un estado cero o su imagen y = m , por lo que un desplazamiento de d −1 funcionará de la misma manera, si es más conveniente. Esto reduce el desplazamiento a 0 en el Caso principal de Mersenne, cuando d = 1.)
La reducción de un mayor producto hacha a menos de 2 m = 2 e 1 - 2 d puede hacerse mediante una o más etapas de reducción sin un desplazamiento.
Si ad ≤ m , entonces es suficiente 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 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 <2 e , la primera reducción producirá un valor en el rango requerido para el caso anterior de dos pasos de reducción para aplicar.
El método de Schrage
Si un producto de doble ancho no está disponible, el método de Schrage , [13] [14] también llamado método de factorización aproximada, [15] puede usarse para calcular ax mod m , pero esto tiene el costo:
- El módulo debe ser representable 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 ⌋ , comúnmente se logra 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 se reduce la eficiencia. una preocupación. Incluso en lenguajes de alto nivel, si el multiplicador a está limitado a √ m , entonces el producto ax de doble ancho puede calcularse usando dos multiplicaciones de ancho simple y reducirse usando las técnicas descritas anteriormente.
Para utilizar el método de Schrage, primer factor m = qa + r , es decir, calcular previamente las constantes auxiliares r = m mod una y q = ⌊ m / una ⌋ = ( m - r ) / una . Luego, en cada iteración, calcule ax ≡ a ( x mod q ) - r ⌊ x / q ⌋ (mode m ) .
Esta igualdad se mantiene porque
así que 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 = my puede calcularse con un producto de ancho simple.
Si una se elige de modo que r ≤ q (y por tanto r / q ≤ 1), entonces el segundo término es también menos de 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 sola suma condicional. [dieciséis]
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 se puede ampliar para permitir una mayor aplicándola de forma recursiva. [15] : 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ía espectacularmente ineficiente. Sin embargo, cada paso también divide x por un cociente q = ⌊ m / a ⌋ cada vez mayor , y rápidamente se alcanza un punto en el que el argumento es 0 y la recursividad puede terminar.
Ejemplo de código C99
Usando código C , el Park-Miller RNG se puede escribir de la siguiente manera:
uint32_t lcg_parkmiller ( uint32_t * estado ) { retorno * estado = ( uint64_t ) * estado * 48271 % 0x7fffffff ; }
Esta función se puede llamar repetidamente para generar números pseudoaleatorios, siempre que la persona que llama 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 enteros de 32 bits puede desbordarse.
Para evitar la división de 64 bits, realice la reducción a mano:
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 ); retorno * estado = x ; }
Para usar solo aritmética de 32 bits, use 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 ; // max: M / Q = A = 48,271 uint32_t rem = * estado % Q ; // max: Q - 1 = 44,487int32_t s = rem * A ; // máx: 44,487 * 48,271 = 2,147,431,977 = 0x7fff3629 int32_t t = div * R ; // max: 48,271 * 3,399 = 164,073,129 int32_t resultado = s - t ;si ( resultado < 0 ) resultado + = M ;return * estado = resultado ; }
o use dos multiplicaciones de 16 × 16 bits:
uint32_t lcg_parkmiller ( uint32_t * estado ) { const uint32_t A = 48271 ;uint32_t bajo = ( * estado & 0x7fff ) * A ; // max: 32,767 * 48,271 = 1,581,695,857 = 0x5e46c371 uint32_t high = ( * estado >> 15 ) * A ; // max: 65,535 * 48,271 = 3,163,439,985 = 0xbc8e4371 uint32_t x = bajo + (( alto & 0xffff ) << 15 ) + ( alto >> 16 ); // máx: 0x5e46c371 + 0x7fff8000 + 0xbc8e = 0xde46ffffx = ( x & 0x7fffffff ) + ( x >> 31 ); retorno * estado = x ; }
Otro generador de Lehmer popular usa el módulo primo 2 32 −5:
uint32_t lcg_rand ( uint32_t * estado ) { retorno * 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 puede ser de 33 bits. producto = ( producto & 0xffffffff ) + 5 * ( uint32_t ) ( producto >> 32 );producto + = 4 ; // Esta suma está garantizada en 32 bits. x = ( uint32_t ) producto + 5 * ( uint32_t ) ( producto >> 32 ); retorno * estado = x - 4 ; }
Muchos otros generadores Lehmer tienen buenas propiedades. El siguiente generador de Lehmer módulo 2 128 requiere compatibilidad con 128 bits del compilador y utiliza un multiplicador calculado por L'Ecuyer. [17] Tiene un período de 2 126 :
estado __int128 sin signo estático ; / * El estado se debe sembrar con un valor impar. * / Vacío de semillas ( sin firmar __int128 semilla ) { estado = semilla << 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 ; estado de retorno >> 64 ; }
El generador calcula un valor impar de 128 bits y devuelve sus 64 bits superiores.
Este generador pasa BigCrush de TestU01 , pero falla la prueba TMFn de PractRand . Esa prueba ha sido diseñada 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 congruentes 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 enteras (si el compilador permite optimizar la declaración constante fuera de un ciclo de cálculo):
uint64_t siguiente ( vacío ) { uint64_t resultado = estado >> 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 ; estado * = mult ; devolver resultado ; }
Sin embargo, debido a que la multiplicación se aplaza, no es adecuado para el hash, ya que la primera llamada simplemente devuelve los 64 bits superiores del estado de inicialización.
Referencias
- ^ WH Payne; JR Rabung; TP Bogyo (1969). "Codificación del generador de números pseudoaleatorios de Lehmer" (PDF) . Comunicaciones de la ACM . 12 (2): 85–86. doi : 10.1145 / 362848.362860 .
- ^ a b c L'Ecuyer, Pierre (junio de 1988). "Generadores de números aleatorios combinados eficientes y portátiles" (PDF) . [Comunicaciones de la ACM] . 31 (6): 742–774. doi : 10.1145 / 62959.62969 .
- ^ Park, Stephen K .; Miller, Keith W. (1988). "Generadores de números aleatorios: los buenos son difíciles de encontrar" (PDF) . Comunicaciones de la ACM . 31 (10): 1192–1201. doi : 10.1145 / 63039.63042 .
- ^ Marsaglia, George (1993). "Correspondencia técnica: observaciones sobre la elección e implementación de generadores de números aleatorios" (PDF) . Comunicaciones de la ACM . 36 (7): 105–108. doi : 10.1145 / 159544.376068 .
- ^ Sullivan, Stephen (1993). "Correspondencia técnica: otra prueba de aleatoriedad" (PDF) . Comunicaciones de la ACM . 36 (7): 108. doi : 10.1145 / 159544.376068 .
- ^ Park, Stephen K .; Miller, Keith W .; Stockmeyer, Paul K. (1988). "Correspondencia técnica: respuesta" (PDF) . Comunicaciones de la ACM . 36 (7): 108-110. doi : 10.1145 / 159544.376068 .
- ^ Vickers, Steve (1981). "Capítulo 5 - Funciones" . Programación básica ZX81 (2ª ed.). Sinclair Research Ltd.
El ZX81 usa p = 65537 & a = 75 [...]
(Tenga en cuenta que el manual de ZX81 indica incorrectamente que 65537 es un número primo de Mersenne que equivale a 2 16 -1. El manual de ZX Spectrum lo solucionó y afirma correctamente que es un número primo de Fermat que equivale a 2 16 +1). - ^ a b Biblioteca científica GNU: otros generadores de números aleatorios
- ^ Knuth, Donald (1981). Algoritmos seminuméricos . El arte de la programación informática . 2 (2ª ed.). Reading, MA: Addison-Wesley Professional. págs. 12-14.
- ^ Bique, Stephen; Rosenberg, Robert (mayo de 2009). Generación rápida de permutaciones y números pseudoaleatorios de alta calidad usando MPI y OpenMP en el Cray XD1 . Cray User Group 2009.
El dado se determina usando aritmética modular, por ejemplo,
lrand48() % 6 + 1,
... ¡La función CRAY RANF solo arroja tres de los seis resultados posibles (cuyos tres lados dependen de la semilla)! - ^ Greenberger, Martin (abril de 1961). "Notas sobre un nuevo generador de números pseudoaleatorios". Revista de la ACM . 8 (2): 163-167. doi : 10.1145 / 321062.321065 .
- ^ L'Ecuyer, Pierre; Tezuka, Shu (octubre de 1991). "Propiedades estructurales para dos clases de generadores de números aleatorios combinados". Matemáticas de la Computación . 57 (196): 735–746. doi : 10.2307 / 2938714 . JSTOR 2938714 .
- ^ Schrage, Linus (junio de 1979). "Un generador de números aleatorios Fortran más portátil" (PDF) . Transacciones ACM en software matemático (TOMS) . 5 (2): 132-138. CiteSeerX 10.1.1.470.6958 . doi : 10.1145 / 355826.355828 .
- ^ Jain, Raj (9 de julio de 2010). "Análisis de rendimiento de sistemas informáticos Capítulo 26: Generación de números aleatorios" (PDF) . págs. 19-22 . Consultado el 31 de octubre de 2017 .
- ^ a b L'Ecuyer, Pierre; Côté, Serge (marzo de 1991). "Implementación de un paquete de números aleatorios con funciones de división". Transacciones ACM en software matemático . 17 (1): 98-111. doi : 10.1145 / 103147.103158 . Esto explora varias implementaciones diferentes de multiplicación modular por una constante.
- ^ Fenerty, Paul (11 de septiembre de 2006). "Método de Schrage" . Consultado el 31 de octubre de 2017 .
- ^ L'Ecuyer, Pierre (enero de 1999). "Tablas de generadores congruenciales lineales de diferentes tamaños y buena estructura de celosía" (PDF) . Matemáticas de la Computación . 68 (225): 249–260. CiteSeerX 10.1.1.34.1024 . doi : 10.1090 / s0025-5718-99-00996-5 .
- Lehmer, DH (1949). "Métodos matemáticos en unidades informáticas a gran escala". Actas de un segundo simposio sobre maquinaria de cálculo digital a gran escala . pp. 141 -146. Señor 0044899 .(versión de la revista: Anales del Laboratorio de Computación de la Universidad de Harvard , Vol. 26 (1951)).
- Steve Park, generadores de números aleatorios
enlaces externos
- Los primos de menos de una potencia de dos pueden ser útiles para elegir módulos. Parte de Prime Pages .