TI-84 Plus 随机数生成器算法

TI-84 Plus Random Number Generator Algorithm

编辑: 我的主要问题是我想在我的计算机上复制 TI-84 加 RNG 算法,所以我可以用 [=21= 这样的语言编写它] 或 Lua,以便更快地进行测试。

我试过用模拟器,结果比计算器慢。

仅供有关人士参考:还有一个question是这样的,但是这个问题的答案只是说如何将已经生成的数字传输到计算机上.我不想要这个。我已经尝试过类似的东西,但我不得不整个周末离开计算器 运行,但仍然没有完成。

TI-Basicrand命令使用的算法是根据TIBasicDev的L'Ecuyer算法。

rand generates a uniformly-distributed pseudorandom number (this page and others will sometimes drop the pseudo- prefix for simplicity) between 0 and 1. rand(n) generates a list of n uniformly-distributed pseudorandom numbers between 0 and 1. seed→rand seeds (initializes) the built-in pseudorandom number generator. The factory default seed is 0.

L'Ecuyer's algorithm is used by TI calculators to generate pseudorandom numbers.

遗憾的是,我无法找到德州仪器 (TI) 发布的任何来源来支持这一说法,因此我无法确定这就是所使用的算法。我也不确定 L'Ecuyer 的算法到底指的是什么。

所使用的算法来自 P. L'Ecuyer 的论文 Efficient and portable combined random number generators

你可以找到论文here and download it for free from here

Ti 计算器使用的算法在 p 的 RHS 侧。 747.我附上了一张图片

我已经将其翻译成 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;
}

请注意,初始种子是 s1 = 12345s2 = 67890

并从 Ti-83(抱歉,我找不到 Ti-84 ROM)模拟器获得输出:

这与我的实现产生的结果相符

我刚刚提高了我的实现的输出精度并得到以下结果:

0.9435973904
0.9083188494
0.1466878273
0.5147019439
0.4058096366
0.7338123019
0.04399198693
0.3393625207

请注意,它们在较低有效数字上与 Ti 的结果不同。这可能是两个处理器(Ti 的 Z80 与我的 X86)执行浮点计算的方式不同。如果是这样,将很难克服这个问题。尽管如此,随机数仍将以相同的序列生成(下面有警告),因为该序列仅依赖于精确的整数数学。

我还使用 long 类型来存储中间值。 Ti 实现依赖于整数溢出存在一些风险(我没有仔细阅读 L'Ecuyer 的论文),在这种情况下,您将不得不调整为 int32_t 或类似的类型来模拟此行为。再次假设处理器的性能相似。

编辑

This site提供代码的Ti-Basic实现如下:

: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

我将其翻译成 C++ 进行测试:

#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;
}

结果如下:

0.9435974025
0.908318861
0.1466878292
0.5147019502
0.405809642
0.7338123114
0.04399198747
0.3393625248
0.9954663411
0.2003402617

与基于原论文的实现相匹配。

我在 Python 中实现了 rand、randInt、randM 和 randBin。感谢理查德的 C 代码。所有实施的命令都按预期工作。您也可以在 this 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))
Here is a C++ program that works:    
#include<cmath>
#include<iostream>
#include<iomanip>
using namespace std;
int main()
{

     double seed1 = 12345;
     double seed2 = 67890;
     double mod1 = 2147483563;
     double mod2 = 2147483399;
     double result;
     for(int i=0; i<10; i++)
     {
         seed1 = seed1*40014-mod1*floor((seed1*40014)/mod1);
         seed2 = seed2*40692-mod2*floor((seed2*40692)/mod2);
         result = (seed1 - seed2)/mod1;
         if(result < 0)
            {result = result + 1;}
         cout<<setprecision(10)<<result<<endl;
  }  
return 0;
}