Un algoritmo de generación de números pseudoaleatorios

En la nota a pie del artículo La privacidad no existe citaba un artículo de la Wikipedia acerca de cómo generar números aleatorios mediante software. En esta entrada voy a poner un ejemplo sencillo desarrollado en Python. Obviamente Python tiene sus propias funciones para ello, pero nada nos impide desarrollar nuestra propia función.

En primer lugar veamos la teoría que hay detrás. Estos algoritmos lo que realmente generan son números pseudoaleatorios porque usan procedimientos sistemáticos. El método más empleado es el de congruencia lineal. Se eligen cuatro números enteros: un módulo m, un multiplicador a, el incremento c y la semilla x0 Estos números deben ser números enteros positivos mayores de 0. El artículo de la Wikipedia también impone que  a, c y x0 deben ser menores que m pero no encuentro ninguna razón para ello. De hecho, al final el artículo entra en contradicción al afirmar que el estándar POSIX para C define a = 1103515245 y m = 32768 Con lo que a > m

Una semilla válida la podemos obtener del reloj interno de la CPU: el RTC. En el ejemplo usaré el método time() de la librería time. Éste nos da el tiempo que ha pasado desde el tiempo Unix: el 1 de Enero de 1970. Como se puede ver, devuelve un número en coma flotante:

>>> time.time()
1487416084.126488

Como necesitamos un entero lo convertimos, el número de segundos será suficiente para el ejemplo.

Siguiendo con la teoría, la semilla nos dará el primer número aleatorio y este mismo número será la semilla para el siguiente. Partiendo de una semilla inicial xn sería:

xn+1 = (axn + c) mod m

Como a, c y m usaré los valores del estándar Posix para C. Este programa genera 10 números aleatorios:

import time

xn = int(time.time()) #Semilla

for i in range(10):
    xn1 = (1103515245 * xn + 12345) % 32768
    print(xn1)
    xn = xn1

Si lo ejecuto este es el resultado:

13566
9311
7852
10101
8970
10107
31128
17905
26070
16471

Lo ejecuto otra vez lo más rápidamente posible:

875
16072
23393
9862
27463
32628
19613
9490
29923
7904

Como vemos hay bastante aleatoriedad. Suficiente para videojuegos y modelos de simulación. Aunque entre ambas ejecuciones sólo pase un segundo los resultados no tienen nada que ver. Si necesitamos diferentes ejecuciones en menos de un segundo podríamos crear un entero fruto de la concatenación de la parte entera y la parte decimal del número que devuelve time()

Llegará un momento en que el resultado obtenido coincida con alguno de los resultados ya obtenidos. Cuando eso pase sabremos que la serie ha acabado y a continuación se repetirán los mismos números. Veámoslo usando unos valores más pequeños: m = 7, a = 5, c = 3 y x0 = 2:

x1 = 5x0 + 3 = 5*2 + 3 = 13 mod 7 = 6
x2 = 5x1 + 3 = 5*6 + 3 = 33 mod 7 = 5
x3 = 5x2 + 3 = 5*5 + 3 = 28 mod 7 = 0
x4 = 5x3 + 3 = 5*0 + 3 = 3 mod 7 = 3
x5 = 5x4 + 3 = 5*3 + 3 = 18 mod 7 = 4
x6 = 5x5 + 3 = 5*4 + 3 = 23 mod 7 = 2
x7 = 5x6 + 3 = 5*2 + 3 = 13 mod 7 = 6

Por lo tanto la serie que se repetirá ad infinitum es: 6, 5, 0, 3, 4, 2, …

De hecho, podríamos decir que la serie empieza por el 2, pues ya x6 = x0. Ahí tenemos la primera coincidencia, antes de llegar a x7 En el algoritmo en Python no se muestra x0 para no tener el primer número pseudoaleatorio de la serie demasiado parecido en diferentes ejecuciones.

En el mejor de los casos posibles el número de elementos de la serie será igual a m. Una mala elección de los valores iniciales m, c y x0 conduciría a una serie muy corta. Según el teorema de Hull-Dobell un generador congruencial tiene una secuencia completa (el número de elementos es igual a m) si y solo si:

  1. m y c son primos entre si.
  2. a -1 es divisible por todos los factores primos de m.
  3. Si 4 divide a m, entonces 4 divide a a – 1.

Además, los algoritmos más eficientes escogerán como m una potencia de 2, pues la operación módulo podrá ser calculada simplemente truncando los 32 o 64 bits más a la derecha.

El estándar POSIX para C ha escogido cuidadosamente los valores cumpliendo con las condiciones del teorema de Hull-Dobell. Veamos:

  1. 1103515245 y 32768 son primos entre si.
  2. El único factor primo de 32768 es 2: 32768 = 215 y 1103515244 | 2 Además es eficiente al ser potencia de 2.
  3. Se cumple tanto 32768 | 4 como 1103515244 | 4

Por lo tanto el algoritmo anterior debería dar una serie 32768 números únicos. Veamos si es cierto con una sencilla modificación del programa anterior que irá añadiendo números a la tupla hasta que se produzca la primera repetición:

import time

xn = int(time.time())
total = []
repeated = False

while (not repeated):
    xn1 = (1103515245 * xn + 12345) % 32768
    xn = xn1
    repeated = xn1 in total
    if (not repeated):
        total.append(xn1)
    
print(total)
print(len(total))

No voy a mostrar la tupla resultante por ser demasiado larga, será bastante con ver el número de elementos que la componen (segundo print):

32768

El número de elementos de la tupla resultante siempre será el mismo por más veces que lo ejecutemos.

Si necesitamos números aleatorios entre 0 y 1 haríamos xn / m :

xn = int(time.time())

for i in range(10):
    xn1 = ((1103515245 * xn + 12345) % 32768) / 32768.00
    print(xn1)
    xn = xn1

Los dos decimales .00 añadidos a c son un “truco” para que Python imprima el resultado xn con decimales: le estamos diciendo al intérprete que queremos un float como resultado de la división. Este es uno de los resultados que da su ejecución:

0.654296875
0.87919062376
0.538866591363
0.578234577039
0.358640995601
0.190647192963
0.728419539036
0.0794397516329
0.638488339469
0.501262168

Deja un comentario

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *

Este sitio usa Akismet para reducir el spam. Aprende cómo se procesan los datos de tus comentarios.