algorithm - simulacion - numeros aleatorios en excel
Algoritmo generador de números aleatorios TI-84 Plus (3)
El algoritmo que se utiliza proviene del papel Generadores de números aleatorios combinados eficientes y portátiles de P. L''Ecuyer.
Puede encontrar el documento aquí y descargarlo gratis desde aquí .
El algoritmo utilizado por las calculadoras Ti está en el lado RHS de p. 747. He incluido una imagen.
Lo traduje a un programa C ++
#include <iostream>
#include <iomanip>
using namespace std;
long s1,s2;
double Uniform(){
long Z,k;
k = s1 / 53668;
s1 = 40014*(s1-k*53668)-k*12211;
if(s1<0)
s1 = s1+2147483563;
k = s2/52774;
s2 = 40692*(s2-k*52774)-k*3791;
if(s2<0)
s2 = s2+2147483399;
Z=s1-s2;
if(Z<1)
Z = Z+2147483562;
return Z*(4.656613e-10);
}
int main(){
s1 = 12345; //Gotta love these seed values!
s2 = 67890;
for(int i=0;i<10;i++)
cout<<std::setprecision(10)<<Uniform()<<endl;
}
Tenga en cuenta que las semillas iniciales son s1 = 12345 y s2 = 67890 .
Y obtuve un resultado de un emulador Ti-83 (lo siento, no pude encontrar un ROM Ti-84):
Esto coincide con lo que produce mi implementación
Acabo de generar la precisión de salida en mi implementación y obtener los siguientes resultados:
0.9435973904
0.9083188494
0.1466878273
0.5147019439
0.4058096366
0.7338123019
0.04399198693
0.3393625207
Tenga en cuenta que divergen de los resultados de Ti en los dígitos menos significativos. Esto puede ser una diferencia en la forma en que los dos procesadores (Ti''s Z80 versus mi X86) realizan cálculos de punto flotante. Si es así, será difícil superar este problema. No obstante, los números aleatorios aún se generarán en la misma secuencia (con la advertencia a continuación) ya que la secuencia se basa únicamente en las matemáticas enteras, que son exactas.
También utilicé el tipo long para almacenar valores intermedios. Existe el riesgo de que la implementación de Ti dependa del desbordamiento de enteros (no leí demasiado el documento de L''Ecuyer), en cuyo caso tendrías que ajustar a int32_t o un tipo similar para emular este comportamiento. Asumiendo, una vez más, que los procesadores funcionan de manera similar.
Editar
Este sitio proporciona una implementación Ti-Basic del código de la siguiente manera:
:2147483563→mod1
:2147483399→mod2
:40014→mult1
:40692→mult2
#The RandSeed Algorithm
:abs(int(n))→n
:If n=0 Then
: 12345→seed1
: 67890→seed2
:Else
: mod(mult1*n,mod1)→seed1
: mod(n,mod2)→seed2
:EndIf
#The rand() Algorithm
:Local result
:mod(seed1*mult1,mod1)→seed1
:mod(seed2*mult2,mod2)→seed2
:(seed1-seed2)/mod1→result
:If result<0
: result+1→result
:Return result
Lo traduje a C ++ para probarlo:
#include <iostream>
#include <iomanip>
using namespace std;
long mod1 = 2147483563;
long mod2 = 2147483399;
long mult1 = 40014;
long mult2 = 40692;
long seed1,seed2;
void Seed(int n){
if(n<0) //Perform an abs
n = -n;
if(n==0){
seed1 = 12345; //Gotta love these seed values!
seed2 = 67890;
} else {
seed1 = (mult1*n)%mod1;
seed2 = n%mod2;
}
}
double Generate(){
double result;
seed1 = (seed1*mult1)%mod1;
seed2 = (seed2*mult2)%mod2;
result = (double)(seed1-seed2)/(double)mod1;
if(result<0)
result = result+1;
return result;
}
int main(){
Seed(0);
for(int i=0;i<10;i++)
cout<<setprecision(10)<<Generate()<<endl;
}
Esto dio los siguientes resultados:
0.9435974025
0.908318861
0.1466878292
0.5147019502
0.405809642
0.7338123114
0.04399198747
0.3393625248
0.9954663411
0.2003402617
que coinciden con los logrados con la implementación basada en el documento original.
Editar: mi pregunta principal es que quiero replicar el algoritmo TI-84 plus RNG en mi computadora, así que puedo escribirlo en un lenguaje como JavaScript o Lua, para probarlo más rápido.
Intenté usar un emulador, pero resultó ser más lento que la calculadora.
Solo para las personas involucradas: Hay otra pregunta como esta, pero la respuesta a esa pregunta solo dice cómo transferir los números ya generados a la computadora. No quiero esto Ya probé algo así, pero tuve que dejar la calculadora funcionando todo el fin de semana, y todavía no estaba hecho.
El algoritmo utilizado por el comando TI-Basic rand es el algoritmo de L''Ecuyer según TIBasicDev .
rand genera un número pseudoaleatorio distribuido uniformemente (esta página y otras a veces soltarán el pseudo prefijo por simplicidad) entre 0 y 1. rand (n) genera una lista de n números pseudoaleatorios distribuidos uniformemente entre 0 y 1. seed → rand semillas (inicializa) el generador de números pseudoaleatorio incorporado. La semilla predeterminada de fábrica es 0.
El algoritmo de L''Ecuyer es utilizado por las calculadoras TI para generar números pseudoaleatorios.
Lamentablemente, no he podido encontrar ninguna fuente publicada por Texas Instruments que respalde esta afirmación, por lo que no puedo asegurar que este sea el algorthm utilizado. Tampoco estoy seguro de a qué se refiere exactamente el algoritmo de L''Ecuyer.
Implementé rand, randInt, randM y randBin en Python. Gracias Richard por el código C. Todos los comandos implementados funcionan como se esperaba. También puedes encontrarlo en este Gist .
import math
class TIprng(object):
def __init__(self):
self.mod1 = 2147483563
self.mod2 = 2147483399
self.mult1 = 40014
self.mult2 = 40692
self.seed1 = 12345
self.seed2 = 67890
def seed(self, n):
n = math.fabs(math.floor(n))
if (n == 0):
self.seed1 = 12345
self.seed2 = 67890
else:
self.seed1 = (self.mult1 * n) % self.mod1
self.seed2 = (n)% self.mod2
def rand(self, times = 0):
# like TI, this will return a list (array in python) if times == 1,
# or an integer if times isn''t specified
if not(times):
self.seed1 = (self.seed1 * self.mult1) % self.mod1
self.seed2 = (self.seed2 * self.mult2)% self.mod2
result = (self.seed1 - self.seed2)/self.mod1
if(result<0):
result = result+1
return result
else:
return [self.rand() for _ in range(times)]
def randInt(self, minimum, maximum, times = 0):
# like TI, this will return a list (array in python) if times == 1,
# or an integer if times isn''t specified
if not(times):
if (minimum < maximum):
return (minimum + math.floor((maximum- minimum + 1) * self.rand()))
else:
return (maximum + math.floor((minimum - maximum + 1) * self.rand()))
else:
return [self.randInt(minimum, maximum) for _ in range(times)]
def randBin(self, numtrials, prob, times = 0):
if not(times):
return sum([(self.rand() < prob) for _ in range(numtrials)])
else:
return [self.randBin(numtrials, prob) for _ in range(times)]
def randM(self, rows, columns):
# this will return an array of arrays
matrixArr = [[0 for x in range(columns)] for x in range(rows)]
# we go from bottom to top, from right to left
for row in reversed(range(rows)):
for column in reversed(range(columns)):
matrixArr[row][column] = self.randInt(-9, 9)
return matrixArr
testPRNG = TIprng()
testPRNG.seed(0)
print(testPRNG.randInt(0,100))
testPRNG.seed(0)
print(testPRNG.randM(3,4))