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

17
编辑:我的主要问题是我想在电脑上复制TI-84 plus RNG算法,这样我就可以用像Javascript或Lua这样的语言编写它,以便更快地测试它。

我尝试使用模拟器,但结果比计算器还慢。

仅针对相关人士:这有另一个问题,但是该问题的答案只是说明如何将已生成的数字传输到计算机上。我不想这样做。我已经尝试过类似的方式,但必须让计算器运行整个周末,但仍未完成。


你能澄清一下你所说的“有人知道TI-84 Plus计算器上的RNG是如何工作的”吗?我不明白你想做什么。你想知道它是如何种子化的吗?你想知道分布有多随机吗?你提到你运行了一个长时间的测试,但没有明确的问题。 - Dan
我编辑了这个问题。 - Potassium Ion
所以你正在尝试“测试它”?到目前为止,您要做什么还不是100%清楚。您只想要像其他问题中那样的大量随机样本吗? - Dan
4个回答

23
所使用的算法来自P. L'Ecuyer的论文“高效便携式组合随机数生成器”。您可以在这里找到该论文,并从这里免费下载。Ti计算器使用的算法位于第747页的RHS侧。我已经包含了一张图片。

L'Ecuyer's Algorithm

我已将此翻译为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)仿真器中获得了输出:

Ti-83 Screenshot

这与我的实现产生的结果相匹配

My computer

我刚刚调整了我的实现输出精度,并得到了以下结果:

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 或类似类型以模拟此行为。再次假设处理器执行类似操作。
编辑
这个网站提供了以下 Ti-Basic 代码实现:This site
: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

这与基于原始论文实现的结果相匹配。


如果您拥有TI-84+ SE,您可以合法使用TI-84+ SE ROM,例如在http://tibasic.com/rom/TI84PlusSE.rom找到的ROM。 - Timtech
您可以从TI网站获取2.55操作系统,但2.43已不再提供。如果您需要下载,请告诉我。我是从TI支持团队那里得到的,但只是因为他们中的某个人仍然在电脑上保留了它。 - Fabian Röling

4
我在Python中实现了rand、randInt、randM和randBin。感谢Richard提供的C代码。所有实现的指令都按预期工作。你也可以在这个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))

3
TI-Basic命令rand所使用的算法是L'Ecuyer's算法,参考TIBasicDev

rand生成一个0到1之间的均匀分布的伪随机数(有时为了简单起见,本页和其他页面会省略伪前缀), rand(n)生成一个由n个0到1之间的均匀分布的伪随机数组成的列表。seed→rand用于初始化内置的伪随机数生成器的种子。出厂默认种子为0。

TI计算器使用L'Ecuyer's算法生成伪随机数。

不幸的是,我没有找到任何由德州仪器公司发布支持这一说法的来源,因此我不能确定这是否是所使用的算法。我也不确定L'Ecuyer's算法确切指的是什么。


P. L’Ecuyer, “Combined Multiple Recursive Random Number Generators”, Operations Research, 44, 5 (1996), 816–822. - Mr. Llama
2
此链接详细描述了该算法:http://tibasicdev.wikidot.com/68k:randseed - John Coleman

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

这并没有回答问题。一旦您拥有足够的声望,您将能够评论任何帖子;相反,提供不需要询问者澄清的答案。- 来自审核 - Brent Worden

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接