{"id":1192,"date":"2015-04-17T11:26:15","date_gmt":"2015-04-17T16:26:15","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=1192"},"modified":"2015-05-17T22:54:13","modified_gmt":"2015-05-18T03:54:13","slug":"random-number-generator-mersenne-twister","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2015\/04\/17\/random-number-generator-mersenne-twister\/","title":{"rendered":"Random Number Generators, Mersenne Twister"},"content":{"rendered":"<div class=\"content\"><!--introduction--><p>This is the first of a multi-part series about the MATLAB random number generators.<\/p><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#4d65cb50-0ae9-47cc-91b8-4cc1c360fe0b\">Taking control<\/a><\/li><li><a href=\"#944ea929-6620-4dd3-9b80-1a23b2e2152a\">Mersenne Twister<\/a><\/li><li><a href=\"#e2ee4889-6eb0-4c07-967b-d07b4e281b9d\">Mersenne primes<\/a><\/li><li><a href=\"#3db68c62-bbc1-4707-99fb-46618b0b4cc8\">Why the name?<\/a><\/li><li><a href=\"#4785d541-7c1d-4a09-be7d-067b628acf19\">Algorithm<\/a><\/li><li><a href=\"#f11678ec-bce7-43f5-932e-6ed6f1196407\">Doubles<\/a><\/li><li><a href=\"#67a5ab4f-d5ec-4c8d-848c-236a724bf91d\">Seeds, Streams and State<\/a><\/li><li><a href=\"#17160ee3-6191-497d-9249-56bd8c03aa8c\">Thanks<\/a><\/li><li><a href=\"#c2eb6a1d-9b21-4f5b-bdd8-41004dd5268a\">References<\/a><\/li><\/ul><\/div><h4>Taking control<a name=\"4d65cb50-0ae9-47cc-91b8-4cc1c360fe0b\"><\/a><\/h4><p>If you issue the following commands at any point in any recent version of MATLAB, you will always get this plot.<\/p><pre class=\"codeinput\">   rng <span class=\"string\">default<\/span>\r\n   hist(randn(10000,1),100)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/random_blog_01.png\" alt=\"\"> <p>The <tt>rng<\/tt> command controls the random number generator that is used by the <tt>rand<\/tt>, <tt>randn<\/tt>, and <tt>randi<\/tt> functions.  When called with the <tt>default<\/tt> parameter, <tt>rng<\/tt> resets the generator to the condition that it has when a fresh MATLAB is started.  To see what that condition is, just call <tt>rng<\/tt> by itself.<\/p><pre class=\"codeinput\">   rng\r\n<\/pre><pre class=\"codeoutput\">ans = \r\n     Type: 'twister'\r\n     Seed: 0\r\n    State: [625x1 uint32]\r\n<\/pre><p>We see that the default generator is <tt>'twister'<\/tt>, the default seed is 0, and that the state is a length 625 vector of unsigned 32-bit integers.<\/p><p>If you ask for <tt>help rng<\/tt>, you will get lots of information, including the fact that there are three modern generators.<\/p><pre>Generator              Description\r\n------------------------------------------------------\r\n'twister'              Mersenne Twister\r\n'combRecursive'        Combined Multiple Recursive\r\n'multFibonacci'        Multiplicative Lagged Fibonacci<\/pre><p>The remainder of today's post is about <tt>'twister'<\/tt>.  I will cover the others in future posts.<\/p><h4>Mersenne Twister<a name=\"944ea929-6620-4dd3-9b80-1a23b2e2152a\"><\/a><\/h4><p>Mersenne Twister is, by far, today's most popular pseudorandom number generator.  It is used by every widely distributed mathematical software package.  It has been available as an option in MATLAB since it was invented and has been the default for almost a decade.<\/p><p>Mersenne Twister was developed by professors Makoto Matsumoto and Takuji Nishimura of Hiroshima University almost twenty years ago. Here is their <a href=\"http:\/\/www.math.sci.hiroshima-u.ac.jp\/~m-mat\/MT\/emt.html\">home page<\/a>.  The C source code is <a href=\"http:\/\/www.math.sci.hiroshima-u.ac.jp\/~m-mat\/MT\/MT2002\/CODES\/mt19937ar.c\">available here<\/a>.<\/p><h4>Mersenne primes<a name=\"e2ee4889-6eb0-4c07-967b-d07b4e281b9d\"><\/a><\/h4><p>Mersenne primes are primes of the form <i>2^p - 1<\/i> where <i>p<\/i> itself is prime.  They are named after a French friar who studied them in the early 17th century.  We learn from Wikipedia that the largest known prime number is the Mersenne prime with <i>p<\/i> equal to 57,885,161. The Mersenne Twister has <i>p<\/i> equal to 19937.  This is tiny as far as Mersenne primes go, but huge as far as random number generators are concerned.<\/p><h4>Why the name?<a name=\"3db68c62-bbc1-4707-99fb-46618b0b4cc8\"><\/a><\/h4><p>Matsumoto explains how the name \"Mersenne Twister\" originated in the Mersenne Twister Home Page.<\/p><pre>MT was firstly named \"Primitive Twisted Generalized Feedback Shift\r\n  Register Sequence\" by a historical reason.\r\nMakoto: Prof. Knuth said in his letter \"the name is mouthful.\"\r\nTakuji: ........\r\na few days later\r\nMakoto: Hi, Takkun, How about \"Mersenne Twister?\" Since it uses Mersenne\r\n  primes, and it shows that it has its ancestor Twisted GFSR.\r\nTakuji: Well.\r\nMakoto: It sounds like a jet coaster, so it sounds quite fast, easy to\r\n  remember and easy to pronounce. Moreover, although it is a secret, it\r\n  hides in its name the initials of the inventors.\r\nTakuji: .......\r\nMakoto: Come on, let's go with MT!\r\nTakuji: ....well, affirmative.\r\nLater, we got a letter from Prof. Knuth saying \"it sounds a nice name.\" :-)<\/pre><h4>Algorithm<a name=\"4785d541-7c1d-4a09-be7d-067b628acf19\"><\/a><\/h4><p>The integer portion of the Mersenne twister algorithm does not involve any arithmetic in the sense of addition, subtraction, multiplication or division. All the operations are shifts, and's, or's, and xor's.<\/p><p>All the elements of the state, except the last, are unsigned 32 bit random integers that form a cache which is carefully generated upon startup. This generation is triggered by a <tt>seed<\/tt>, a single integer that initiates the whole process.<\/p><p>The last element of the state is a pointer into the cache.  Each request for a random integer causes an element to be withdrawn from the cache and the pointer incremented.  The element is \"tempered\" with additional logical operations to improve the randomness.  When the pointer reaches the end of the cache, the cache is refilled with another 623 elements.<\/p><p>The algorithm is analyzed by investigating the group theoretic properties of the permutation and tempering operations.  The parameters have been chosen so that the period is the Mersenne prime 2^19937-1.  This period is much longer than any other random number generator proposed before or since and is one of the reasons for MT's popularity.<\/p><p>By design, the results generated satisfy an equidistribution property in a 623-dimensional cube.<\/p><h4>Doubles<a name=\"f11678ec-bce7-43f5-932e-6ed6f1196407\"><\/a><\/h4><p>Here is the function in the Mersenne Twister <a href=\"http:\/\/www.math.sci.hiroshima-u.ac.jp\/~m-mat\/MT\/MT2002\/CODES\/mt19937ar.c\">source code<\/a> that converts a pair of random uint32s into a random double. You can see that it takes the top 27 bits of one int and the top 26 bits of the other, sticks them together, and multiplies by what MATLAB would call <tt>eps\/2<\/tt>.  This is the only place in the code where floating point arithmetic is involved.<\/p><pre>double genrand_res53(void)\r\n{\r\n    unsigned long a=genrand_int32()&gt;&gt;5, b=genrand_int32()&gt;&gt;6;\r\n    return(a*67108864.0+b)*(1.0\/9007199254740992.0);\r\n}<\/pre><p>The result would be zero in the highly unlikely event that both <tt>a<\/tt> and <tt>b<\/tt> are zero.  If that happens, the MATLAB interface rejects the result and calls this function again.  So the smallest double results when <tt>a<\/tt> equals zero and <tt>b<\/tt> equals 1.  The largest double results when both <tt>a<\/tt> and <tt>b<\/tt> are all 1's.  Consequently, the output from <tt>rand<\/tt> is in the closed interval<\/p><p>$$ 2^{-53} \\leq x \\leq 1-2^{-53} $$<\/p><h4>Seeds, Streams and State<a name=\"67a5ab4f-d5ec-4c8d-848c-236a724bf91d\"><\/a><\/h4><p>For more about random number seeds, streams, and state, see Peter Perkins, guest blogger in <a href=\"https:\/\/blogs.mathworks.com\/loren\/2011\/07\/07\/simpler-control-of-random-number-generation-in-matlab\">Loren's blog<\/a>. And, of course, see <a href=\"https:\/\/www.mathworks.com\/help\/releases\/R2011a\/techdoc\/\/math\/example-demos-rngdemo.html\">the documentation<\/a>.<\/p><h4>Thanks<a name=\"17160ee3-6191-497d-9249-56bd8c03aa8c\"><\/a><\/h4><p>Thanks to Peter Perkins for the work he has done on our random number suite over the years, and for enlightening me.<\/p><h4>References<a name=\"c2eb6a1d-9b21-4f5b-bdd8-41004dd5268a\"><\/a><\/h4><p>M. Matsumoto and T. Nishimura. \"Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudorandom Number Generator.\" ACM Transactions on Modeling and Computer Simulation, 8(1):3-30. 1998. Available online at: <a href=\"http:\/\/www.math.sci.hiroshima-u.ac.jp\/~m-mat\/MT\/ARTICLES\/mt.pdf\">&lt;http:\/\/www.math.sci.hiroshima-u.ac.jp\/~m-mat\/MT\/ARTICLES\/mt.pdf<\/a>&gt;.<\/p><script language=\"JavaScript\"> <!-- \r\n    function grabCode_a977e1b84e3742938c54460ea9be8bb6() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='a977e1b84e3742938c54460ea9be8bb6 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' a977e1b84e3742938c54460ea9be8bb6';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        copyright = 'Copyright 2015 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add copyright line at the bottom if specified.\r\n        if (copyright.length > 0) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\\n');\r\n\r\n        d.title = title + ' (MATLAB code)';\r\n        d.close();\r\n    }   \r\n     --> <\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_a977e1b84e3742938c54460ea9be8bb6()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n      the MATLAB code <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; R2014b<br><\/p><\/div><!--\r\na977e1b84e3742938c54460ea9be8bb6 ##### SOURCE BEGIN #####\r\n%% Random Number Generators, Mersenne Twister\r\n% This is the first of a multi-part series about the MATLAB random number\r\n% generators.\r\n\r\n%% Taking control\r\n% If you issue the following commands at any point in any recent version of\r\n% MATLAB, you will always get this plot.\r\n\r\n   rng default\r\n   hist(randn(10000,1),100)\r\n\r\n%%\r\n% The |rng| command controls the random number generator that is used by\r\n% the |rand|, |randn|, and |randi| functions.  When called with the |default|\r\n% parameter, |rng| resets the generator to the condition that it has when\r\n% a fresh MATLAB is started.  To see what that condition is, just call |rng|\r\n% by itself.\r\n\r\n   rng\r\n\r\n%%\r\n% We see that the default generator is |'twister'|, the default seed is 0,\r\n% and that the state is a length 625 vector of unsigned 32-bit integers.\r\n\r\n%%\r\n% If you ask for |help rng|, you will get lots of information, including\r\n% the fact that there are three modern generators.\r\n%\r\n%  Generator              Description\r\n%  REPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASH\r\n%  'twister'              Mersenne Twister\r\n%  'combRecursive'        Combined Multiple Recursive\r\n%  'multFibonacci'        Multiplicative Lagged Fibonacci\r\n\r\n%%\r\n% The remainder of today's post is about |'twister'|.  I will cover the\r\n% others in future posts.\r\n\r\n%% Mersenne Twister\r\n% Mersenne Twister is, by far, today's most popular pseudorandom number\r\n% generator.  It is used by every widely distributed mathematical software\r\n% package.  It has been available as an option in MATLAB since it was\r\n% invented and has been the default for almost a decade.\r\n\r\n%%\r\n% Mersenne Twister was developed by professors Makoto Matsumoto and\r\n% Takuji Nishimura of Hiroshima University almost twenty years ago.\r\n% Here is their <http:\/\/www.math.sci.hiroshima-u.ac.jp\/~m-mat\/MT\/emt.html\r\n% home page>.  The C source code is\r\n% <http:\/\/www.math.sci.hiroshima-u.ac.jp\/~m-mat\/MT\/MT2002\/CODES\/mt19937ar.c\r\n% available here>.\r\n\r\n%% Mersenne primes\r\n% Mersenne primes are primes of the form _2^p - 1_ where _p_ itself is\r\n% prime.  They are named after a French friar who studied them in the\r\n% early 17th century.  We learn from Wikipedia that the largest known\r\n% prime number is the Mersenne prime with _p_ equal to 57,885,161. \r\n% The Mersenne Twister has _p_ equal to 19937.  This is tiny as far as\r\n% Mersenne primes go, but huge as far as random number generators are\r\n% concerned.\r\n\r\n%% Why the name?\r\n% Matsumoto explains how the name \"Mersenne Twister\" originated in the\r\n% Mersenne Twister Home Page.\r\n%\r\n%  MT was firstly named \"Primitive Twisted Generalized Feedback Shift\r\n%    Register Sequence\" by a historical reason.\r\n%  Makoto: Prof. Knuth said in his letter \"the name is mouthful.\" \r\n%  Takuji: ........ \r\n%  a few days later \r\n%  Makoto: Hi, Takkun, How about \"Mersenne Twister?\" Since it uses Mersenne\r\n%    primes, and it shows that it has its ancestor Twisted GFSR. \r\n%  Takuji: Well. \r\n%  Makoto: It sounds like a jet coaster, so it sounds quite fast, easy to\r\n%    remember and easy to pronounce. Moreover, although it is a secret, it\r\n%    hides in its name the initials of the inventors. \r\n%  Takuji: .......\r\n%  Makoto: Come on, let's go with MT! \r\n%  Takuji: ....well, affirmative. \r\n%  Later, we got a letter from Prof. Knuth saying \"it sounds a nice name.\" :-)\r\n\r\n%% Algorithm\r\n% The integer portion of the Mersenne twister algorithm does not involve any\r\n% arithmetic in the sense of addition, subtraction, multiplication or division.\r\n% All the operations are shifts, and's, or's, and xor's.\r\n\r\n%%\r\n% All the elements of the state, except the last, are unsigned 32 bit random\r\n% integers that form a cache which is carefully generated upon startup.\r\n% This generation is triggered by a |seed|, a single integer that initiates\r\n% the whole process.\r\n\r\n%%\r\n% The last element of the state is a pointer into the cache.  Each request\r\n% for a random integer causes an element to be withdrawn from the cache\r\n% and the pointer incremented.  The element is \"tempered\" with additional\r\n% logical operations to improve the randomness.  When the pointer reaches\r\n% the end of the cache, the cache is refilled with another 623 elements.\r\n\r\n%%\r\n% The algorithm is analyzed by investigating the group theoretic properties\r\n% of the permutation and tempering operations.  The parameters have been\r\n% chosen so that the period is the Mersenne prime 2^19937-1.  This period\r\n% is much longer than any other random number generator proposed before or\r\n% since and is one of the reasons for MT's popularity.\r\n\r\n%%\r\n% By design, the results generated satisfy an equidistribution property in a\r\n% 623-dimensional cube.  \r\n\r\n%% Doubles\r\n% Here is the function in the Mersenne Twister \r\n% <http:\/\/www.math.sci.hiroshima-u.ac.jp\/~m-mat\/MT\/MT2002\/CODES\/mt19937ar.c\r\n% source code> that converts a pair of random uint32s into a random double.\r\n% You can see that it takes the top 27 bits of one int and the top 26 bits\r\n% of the other, sticks them together, and multiplies by what MATLAB would\r\n% call |eps\/2|.  This is the only place in the code where floating point\r\n% arithmetic is involved.\r\n%\r\n%  double genrand_res53(void) \r\n%  { \r\n%      unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6; \r\n%      return(a*67108864.0+b)*(1.0\/9007199254740992.0); \r\n%  } \r\n%\r\n% The result would be zero in the highly unlikely event that both |a| and |b|\r\n% are zero.  If that happens, the MATLAB interface rejects the result and calls\r\n% this function again.  So the smallest double results when |a| equals zero\r\n% and |b| equals 1.  The largest double results when both |a| and |b| are\r\n% all 1's.  Consequently, the output from |rand| is in the closed interval\r\n%\r\n% $$ 2^{-53} \\leq x \\leq 1-2^{-53} $$\r\n\r\n\r\n%% Seeds, Streams and State\r\n% For more about random number seeds, streams, and state, see Peter Perkins,\r\n% guest blogger in\r\n% <https:\/\/blogs.mathworks.com\/loren\/2011\/07\/07\/simpler-control-of-random-number-generation-in-matlab\r\n% Loren's blog>.\r\n% And, of course, see\r\n% <https:\/\/www.mathworks.com\/help\/releases\/R2011a\/techdoc\/\/math\/example-demos-rngdemo.html\r\n% the documentation>.\r\n\r\n%% Thanks\r\n% Thanks to Peter Perkins for the work he has done on our random number\r\n% suite over the years, and for enlightening me.\r\n\r\n%% References\r\n% M. Matsumoto and T. Nishimura. \"Mersenne Twister: A 623-Dimensionally\r\n% Equidistributed Uniform Pseudorandom Number Generator.\"\r\n% ACM Transactions on Modeling and Computer Simulation, 8(1):3-30. 1998.\r\n% Available online at:\r\n% <http:\/\/www.math.sci.hiroshima-u.ac.jp\/~m-mat\/MT\/ARTICLES\/mt.pdf \r\n% http:\/\/www.math.sci.hiroshima-u.ac.jp\/~m-mat\/MT\/ARTICLES\/mt.pdf>.\r\n\r\n##### SOURCE END ##### a977e1b84e3742938c54460ea9be8bb6\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img decoding=\"async\"  class=\"img-responsive\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/random_blog_01.png\" onError=\"this.style.display ='none';\" \/><\/div><!--introduction--><p>This is the first of a multi-part series about the MATLAB random number generators.... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2015\/04\/17\/random-number-generator-mersenne-twister\/\">read more >><\/a><\/p>","protected":false},"author":78,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[11,4,8,27],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/1192"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/users\/78"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/comments?post=1192"}],"version-history":[{"count":2,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/1192\/revisions"}],"predecessor-version":[{"id":1194,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/1192\/revisions\/1194"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=1192"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=1192"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=1192"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}