RandomMT19937.st
changeset 3340 8d3e3b5152c5
child 3348 d8e779ae0ab1
equal deleted inserted replaced
3339:72f6d17dc008 3340:8d3e3b5152c5
       
     1 "{ Package: 'stx:libbasic2' }"
       
     2 
       
     3 Object subclass:#RandomMT19937
       
     4 	instanceVariableNames:'n m matrixA upperMask lowerMask mt mti'
       
     5 	classVariableNames:''
       
     6 	poolDictionaries:''
       
     7 	category:'Magnitude-Numbers'
       
     8 !
       
     9 
       
    10 !RandomMT19937 class methodsFor:'documentation'!
       
    11 
       
    12 documentation
       
    13 "
       
    14     NO WARRANTY
       
    15 
       
    16     A Mersienne Twister based on code by Takuji Nishimura and Makoto Matsumoto.
       
    17 
       
    18     Before using, initialize the state by using initGen(seed)  
       
    19     or initByArray(init_key, key_length).
       
    20 
       
    21     RandomMT19937 new nextInteger
       
    22     (RandomMT19937 new:5489) nextInteger
       
    23     (RandomMT19937 newByArray:{16r123 . 16r234 .16r345 . 16r456}) nextInteger
       
    24 
       
    25     Please read:
       
    26         Wikipedia article on Mersienne Twister end esp. MT19937
       
    27         MT home page http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
       
    28 
       
    29     [see also:]
       
    30         RandomGenerator - the default; uses the machines /dev/random if available
       
    31         Random  - fast, but generates less quality random numbers
       
    32         RandomTT800 - another random generator
       
    33         RandomParkMiller - another random generator
       
    34 
       
    35     [author:]
       
    36         Original algorithm by Takuji Nishimura and Makoto Matsumoto
       
    37         Ported to Smalltalk by Claus Gittinger.
       
    38 "
       
    39 !
       
    40 
       
    41 test
       
    42 "
       
    43     compare generated numbers a gainst reference implementation
       
    44     (see http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.out)
       
    45                                                                                 [exBegin]
       
    46     |l rnd|
       
    47 
       
    48     rnd := (RandomMT19937 newByArray:{16r123 . 16r234 .16r345 . 16r456}).
       
    49 
       
    50     #(
       
    51         1067595299  955945823  477289528 4107218783 4228976476 
       
    52         3344332714 3355579695  227628506  810200273 2591290167 
       
    53         2560260675 3242736208  646746669 1479517882 4245472273 
       
    54         1143372638 3863670494 3221021970 1773610557 1138697238 
       
    55         1421897700 1269916527 2859934041 1764463362 3874892047 
       
    56         3965319921   72549643 2383988930 2600218693 3237492380 
       
    57         2792901476  725331109  605841842  271258942  715137098 
       
    58         3297999536 1322965544 4229579109 1395091102 3735697720 
       
    59         2101727825 3730287744 2950434330 1661921839 2895579582 
       
    60         2370511479 1004092106 2247096681 2111242379 3237345263 
       
    61         4082424759  219785033 2454039889 3709582971  835606218 
       
    62         2411949883 2735205030  756421180 2175209704 1873865952 
       
    63         2762534237 4161807854 3351099340  181129879 3269891896 
       
    64          776029799 2218161979 3001745796 1866825872 2133627728 
       
    65           34862734 1191934573 3102311354 2916517763 1012402762 
       
    66         2184831317 4257399449 2899497138 3818095062 3030756734 
       
    67         1282161629  420003642 2326421477 2741455717 1278020671 
       
    68         3744179621  271777016 2626330018 2560563991 3055977700 
       
    69         4233527566 1228397661 3595579322 1077915006 2395931898 
       
    70         1851927286 3013683506 1999971931 3006888962 1049781534 
       
    71         1488758959 3491776230  104418065 2448267297 3075614115 
       
    72         3872332600  891912190 3936547759 2269180963 2633455084 
       
    73         1047636807 2604612377 2709305729 1952216715  207593580 
       
    74         2849898034  670771757 2210471108  467711165  263046873 
       
    75         3569667915 1042291111 3863517079 1464270005 2758321352 
       
    76         3790799816 2301278724 3106281430    7974801 2792461636 
       
    77          555991332  621766759 1322453093  853629228  686962251 
       
    78         1455120532  957753161 1802033300 1021534190 3486047311 
       
    79         1902128914 3701138056 4176424663 1795608698  560858864 
       
    80         3737752754 3141170998 1553553385 3367807274  711546358 
       
    81         2475125503  262969859  251416325 2980076994 1806565895 
       
    82          969527843 3529327173 2736343040 2987196734 1649016367 
       
    83         2206175811 3048174801 3662503553 3138851612 2660143804 
       
    84         1663017612 1816683231  411916003 3887461314 2347044079 
       
    85         1015311755 1203592432 2170947766 2569420716  813872093 
       
    86         1105387678 1431142475  220570551 4243632715 4179591855 
       
    87         2607469131 3090613241  282341803 1734241730 1391822177 
       
    88         1001254810  827927915 1886687171 3935097347 2631788714 
       
    89         3905163266  110554195 2447955646 3717202975 3304793075 
       
    90         3739614479 3059127468  953919171 2590123714 1132511021 
       
    91         3795593679 2788030429  982155079 3472349556  859942552 
       
    92         2681007391 2299624053  647443547  233600422  608168955 
       
    93         3689327453 1849778220 1608438222 3968158357 2692977776 
       
    94         2851872572  246750393 3582818628 3329652309 4036366910 
       
    95         1012970930  950780808 3959768744 2538550045  191422718 
       
    96         2658142375 3276369011 2927737484 1234200027 1920815603 
       
    97         3536074689 1535612501 2184142071 3276955054  428488088 
       
    98         2378411984 4059769550 3913744741 2732139246   64369859 
       
    99         3755670074  842839565 2819894466 2414718973 1010060670 
       
   100         1839715346 2410311136  152774329 3485009480 4102101512 
       
   101         2852724304  879944024 1785007662 2748284463 1354768064 
       
   102         3267784736 2269127717 3001240761 3179796763  895723219 
       
   103          865924942 4291570937   89355264 1471026971 4114180745 
       
   104         3201939751 2867476999 2460866060 3603874571 2238880432 
       
   105         3308416168 2072246611 2755653839 3773737248 1709066580 
       
   106         4282731467 2746170170 2832568330  433439009 3175778732 
       
   107           26248366 2551382801  183214346 3893339516 1928168445 
       
   108         1337157619 3429096554 3275170900 1782047316 4264403756 
       
   109         1876594403 4289659572 3223834894 1728705513 4068244734 
       
   110         2867840287 1147798696  302879820 1730407747 1923824407 
       
   111         1180597908 1569786639  198796327  560793173 2107345620 
       
   112         2705990316 3448772106 3678374155  758635715  884524671 
       
   113          486356516 1774865603 3881226226 2635213607 1181121587 
       
   114         1508809820 3178988241 1594193633 1235154121  326117244 
       
   115         2304031425  937054774 2687415945 3192389340 2003740439 
       
   116         1823766188 2759543402   10067710 1533252662 4132494984 
       
   117           82378136  420615890 3467563163  541562091 3535949864 
       
   118         2277319197 3330822853 3215654174 4113831979 4204996991 
       
   119         2162248333 3255093522 2219088909 2978279037  255818579 
       
   120         2859348628 3097280311 2569721123 1861951120 2907080079 
       
   121         2719467166  998319094 2521935127 2404125338  259456032 
       
   122         2086860995 1839848496 1893547357 2527997525 1489393124 
       
   123         2860855349   76448234 2264934035  744914583 2586791259 
       
   124         1385380501   66529922 1819103258 1899300332 2098173828 
       
   125         1793831094  276463159  360132945 4178212058  595015228 
       
   126          177071838 2800080290 1573557746 1548998935  378454223 
       
   127         1460534296 1116274283 3112385063 3709761796  827999348 
       
   128         3580042847 1913901014  614021289 4278528023 1905177404 
       
   129           45407939 3298183234 1184848810 3644926330 3923635459 
       
   130         1627046213 3677876759  969772772 1160524753 1522441192 
       
   131          452369933 1527502551  832490847 1003299676 1071381111 
       
   132         2891255476  973747308 4086897108 1847554542 3895651598 
       
   133         2227820339 1621250941 2881344691 3583565821 3510404498 
       
   134          849362119  862871471  797858058 2867774932 2821282612 
       
   135         3272403146 3997979905  209178708 1805135652    6783381 
       
   136         2823361423  792580494 4263749770  776439581 3798193823 
       
   137         2853444094 2729507474 1071873341 1329010206 1289336450 
       
   138         3327680758 2011491779   80157208  922428856 1158943220 
       
   139         1667230961 2461022820 2608845159  387516115 3345351910 
       
   140         1495629111 4098154157 3156649613 3525698599 4134908037 
       
   141          446713264 2137537399 3617403512  813966752 1157943946 
       
   142         3734692965 1680301658 3180398473 3509854711 2228114612 
       
   143         1008102291  486805123  863791847 3189125290 1050308116 
       
   144         3777341526 4291726501  844061465 1347461791 2826481581 
       
   145          745465012 2055805750 4260209475 2386693097 2980646741 
       
   146          447229436 2077782664 1232942813 4023002732 1399011509 
       
   147         3140569849 2579909222 3794857471  900758066 2887199683 
       
   148         1720257997 3367494931 2668921229  955539029 3818726432 
       
   149         1105704962 3889207255 2277369307 2746484505 1761846513 
       
   150         2413916784 2685127085 4240257943 1166726899 4215215715 
       
   151         3082092067 3960461946 1663304043 2087473241 4162589986 
       
   152         2507310778 1579665506  767234210  970676017  492207530 
       
   153         1441679602 1314785090 3262202570 3417091742 1561989210 
       
   154         3011406780 1146609202 3262321040 1374872171 1634688712 
       
   155         1280458888 2230023982  419323804 3262899800   39783310 
       
   156         1641619040 1700368658 2207946628 2571300939 2424079766 
       
   157          780290914 2715195096 3390957695  163151474 2309534542 
       
   158         1860018424  555755123  280320104 1604831083 2713022383 
       
   159         1728987441 3639955502  623065489 3828630947 4275479050 
       
   160         3516347383 2343951195 2430677756  635534992 3868699749 
       
   161          808442435 3070644069 4282166003 2093181383 2023555632 
       
   162         1568662086 3422372620 4134522350 3016979543 3259320234 
       
   163         2888030729 3185253876 4258779643 1267304371 1022517473 
       
   164          815943045  929020012 2995251018 3371283296 3608029049 
       
   165         2018485115  122123397 2810669150 1411365618 1238391329 
       
   166         1186786476 3155969091 2242941310 1765554882  279121160 
       
   167         4279838515 1641578514 3796324015   13351065  103516986 
       
   168         1609694427  551411743 2493771609 1316337047 3932650856 
       
   169         4189700203  463397996 2937735066 1855616529 2626847990 
       
   170           55091862 3823351211  753448970 4045045500 1274127772 
       
   171         1124182256   92039808 2126345552  425973257  386287896 
       
   172         2589870191 1987762798 4084826973 2172456685 3366583455 
       
   173         3602966653 2378803535 2901764433 3716929006 3710159000 
       
   174         2653449155 3469742630 3096444476 3932564653 2595257433 
       
   175          318974657 3146202484  853571438  144400272 3768408841 
       
   176          782634401 2161109003  570039522 1886241521   14249488 
       
   177         2230804228 1604941699 3928713335 3921942509 2155806892 
       
   178          134366254  430507376 1924011722  276713377  196481886 
       
   179         3614810992 1610021185 1785757066  851346168 3761148643 
       
   180         2918835642 3364422385 3012284466 3735958851 2643153892 
       
   181         3778608231 1164289832  205853021 2876112231 3503398282 
       
   182         3078397001 3472037921 1748894853 2740861475  316056182 
       
   183         1660426908  168885906  956005527 3984354789  566521563 
       
   184         1001109523 1216710575 2952284757 3834433081 3842608301 
       
   185         2467352408 3974441264 3256601745 1409353924 1329904859 
       
   186         2307560293 3125217879 3622920184 3832785684 3882365951 
       
   187         2308537115 2659155028 1450441945 3532257603 3186324194 
       
   188         1225603425 1124246549  175808705 3009142319 2796710159 
       
   189         3651990107  160762750 1902254979 1698648476 1134980669 
       
   190          497144426 3302689335 4057485630 3603530763 4087252587 
       
   191          427812652  286876201  823134128 1627554964 3745564327 
       
   192         2589226092 4202024494   62878473 3275585894 3987124064 
       
   193         2791777159 1916869511 2585861905 1375038919 1403421920 
       
   194           60249114 3811870450 3021498009 2612993202  528933105 
       
   195         2757361321 3341402964 2621861700  273128190 4015252178 
       
   196         3094781002 1621621288 2337611177 1796718448 1258965619 
       
   197         4241913140 2138560392 3022190223 4174180924  450094611 
       
   198         3274724580  617150026 2704660665 1469700689 1341616587 
       
   199          356715071 1188789960 2278869135 1766569160 2795896635 
       
   200           57824704 2893496380 1235723989 1630694347 3927960522 
       
   201          428891364 1814070806 2287999787 4125941184 3968103889 
       
   202         3548724050 1025597707 1404281500 2002212197   92429143 
       
   203         2313943944 2403086080 3006180634 3561981764 1671860914 
       
   204         1768520622 1803542985  844848113 3006139921 1410888995 
       
   205         1157749833 2125704913 1789979528 1799263423  741157179 
       
   206         2405862309  767040434 2655241390 3663420179 2172009096 
       
   207         2511931187 1680542666  231857466 1154981000  157168255 
       
   208         1454112128 3505872099 1929775046 2309422350 2143329496 
       
   209         2960716902  407610648 2938108129 2581749599  538837155 
       
   210         2342628867  430543915  740188568 1937713272 3315215132 
       
   211         2085587024 4030765687  766054429 3517641839  689721775 
       
   212         1294158986 1753287754 4202601348 1974852792   33459103 
       
   213         3568087535 3144677435 1686130825 4134943013 3005738435 
       
   214         3599293386  426570142  754104406 3660892564 1964545167 
       
   215          829466833  821587464 1746693036 1006492428 1595312919 
       
   216         1256599985 1024482560 1897312280 2902903201  691790057 
       
   217         1037515867 3176831208 1968401055 2173506824 1089055278 
       
   218         1748401123 2941380082  968412354 1818753861 2973200866 
       
   219         3875951774 1119354008 3988604139 1647155589 2232450826 
       
   220         3486058011 3655784043 3759258462  847163678 1082052057 
       
   221          989516446 2871541755 3196311070 3929963078  658187585 
       
   222         3664944641 2175149170 2203709147 2756014689 2456473919 
       
   223         3890267390 1293787864 2830347984 3059280931 4158802520 
       
   224         1561677400 2586570938  783570352 1355506163   31495586 
       
   225         3789437343 3340549429 2092501630  896419368  671715824 
       
   226         3530450081 3603554138 1055991716 3442308219 1499434728 
       
   227         3130288473 3639507000   17769680 2259741420  487032199 
       
   228         4227143402 3693771256 1880482820 3924810796  381462353 
       
   229         4017855991 2452034943 2736680833 2209866385 2128986379 
       
   230          437874044  595759426  641721026 1636065708 3899136933 
       
   231          629879088 3591174506  351984326 2638783544 2348444281 
       
   232         2341604660 2123933692  143443325 1525942256  364660499 
       
   233          599149312  939093251 1523003209  106601097  376589484 
       
   234         1346282236 1297387043  764598052 3741218111  933457002 
       
   235         1886424424 3219631016  525405256 3014235619  323149677 
       
   236         2038881721 4100129043 2851715101 2984028078 1888574695 
       
   237         2014194741 3515193880 4180573530 3461824363 2641995497 
       
   238         3179230245 2902294983 2217320456 4040852155 1784656905 
       
   239         3311906931   87498458 2752971818 2635474297 2831215366 
       
   240         3682231106 2920043893 3772929704 2816374944  309949752 
       
   241         2383758854  154870719  385111597 1191604312 1840700563 
       
   242          872191186 2925548701 1310412747 2102066999 1504727249 
       
   243         3574298750 1191230036 3330575266 3180292097 3539347721 
       
   244          681369118 3305125752 3648233597  950049240 4173257693 
       
   245         1760124957  512151405  681175196  580563018 1169662867 
       
   246         4015033554 2687781101  699691603 2673494188 1137221356 
       
   247          123599888  472658308 1053598179 1012713758 3481064843 
       
   248         3759461013 3981457956 3830587662 1877191791 3650996736 
       
   249          988064871 3515461600 4089077232 2225147448 1249609188 
       
   250         2643151863 3896204135 2416995901 1397735321 3460025646 
       
   251     ) do:[:expected |
       
   252         |got|
       
   253 
       
   254         got := rnd nextInteger.
       
   255         self assert:(got = expected)
       
   256     ].
       
   257                                                                                 [exEnd]
       
   258     "
       
   259 ! !
       
   260 
       
   261 !RandomMT19937 class methodsFor:'instance creation'!
       
   262 
       
   263 new
       
   264     ^ self basicNew 
       
   265         initialize;
       
   266         initGenRand:(Time millisecondClockValue)
       
   267 !
       
   268 
       
   269 new:seed
       
   270     ^ self basicNew 
       
   271         initialize;
       
   272         initGenRand:seed
       
   273 !
       
   274 
       
   275 newByArray:keyArray
       
   276     ^ self basicNew 
       
   277         initialize;
       
   278         initByArray:keyArray
       
   279 ! !
       
   280 
       
   281 !RandomMT19937 methodsFor:'initialization'!
       
   282 
       
   283 initByArray:keyArray
       
   284     |i j k kLen|
       
   285 
       
   286     self initGenRand:19650218.
       
   287 
       
   288     kLen := keyArray size.
       
   289     i := 1. j := 0.
       
   290     k := (n max: kLen).
       
   291 
       
   292     [k > 0] whileTrue:[
       
   293         |t|
       
   294 
       
   295         t := (((mt at:i+1) bitXor: (((mt at:i) bitXor: ((mt at:i) >> 30)) * 1664525))
       
   296                         + (keyArray at:j+1) + j). 
       
   297         mt at:i+1 put:(t bitAnd:16rFFFFFFFF).
       
   298         i := i+1. j := j+1.
       
   299         (i >= n) ifTrue:[ mt at:1 put:(mt at:n). i := 1. ].
       
   300         (j >= kLen) ifTrue:[ j := 0].
       
   301         k := k - 1.
       
   302     ].
       
   303 
       
   304     k := n-1.
       
   305     [k > 0] whileTrue:[
       
   306         |t|
       
   307 
       
   308         t := ((mt at:i+1) bitXor:(((mt at:i) bitXor: ((mt at:i) >> 30)) * 1566083941)) - i.
       
   309         mt at:i+1 put:(t bitAnd:16rFFFFFFFF).
       
   310         i := i+1.
       
   311         (i >= n) ifTrue:[ mt at:1 put:(mt at:n). i := 1. ].
       
   312         k := k - 1.
       
   313     ].
       
   314 
       
   315     mt at:1 put:16r80000000. "/ MSB is 1; assuring non-zero initial array 
       
   316 !
       
   317 
       
   318 initGenRand:seedUsed
       
   319     |s|
       
   320 
       
   321     s := seedUsed.
       
   322 
       
   323     mt at:1 put:(s bitAnd:16rFFFFFFFF). 
       
   324     mti := 1.
       
   325     [mti < n] whileTrue:[
       
   326         |t|
       
   327 
       
   328         t := (1812433253
       
   329                 * ((mt at:(mti)) bitXor:((mt at:(mti)) >> 30)))
       
   330              + mti.
       
   331 
       
   332         mt at:(mti+1) put:(t bitAnd:16rFFFFFFFF).
       
   333         mti := mti + 1.
       
   334     ].
       
   335 
       
   336     "
       
   337      self new nextInteger
       
   338     "
       
   339 !
       
   340 
       
   341 initialize
       
   342     n := 624.
       
   343     m := 397.
       
   344     matrixA := 16r9908b0df.      "/ constant vector a
       
   345     upperMask := 16r80000000.   "/ most significant w-r bits 
       
   346     lowerMask := 16r7fffffff.   "/ least significant r bits 
       
   347 
       
   348     mt := Array new:n.          "/ the array for the state vector 
       
   349     mti := n + 1.               "/ mti==N+1 means mt[N] is not initialized */
       
   350 ! !
       
   351 
       
   352 !RandomMT19937 methodsFor:'random numbers'!
       
   353 
       
   354 nextBoolean
       
   355     "generates the next integer in 0..FFFFFFFF"
       
   356 
       
   357     ^ self nextInteger > 16r7FFFFFFF
       
   358 !
       
   359 
       
   360 nextInteger
       
   361     "generates the next integer in 0..FFFFFFFF"
       
   362 
       
   363     | y mag01 |
       
   364 
       
   365     mag01 := { 0 . matrixA }.
       
   366     "/ mag01[x] = x * MATRIX_A  for x=0,1 */
       
   367 
       
   368     (mti >= n) ifTrue:[ "/ generate N words at one time
       
   369         |kk|
       
   370 
       
   371         (mti == (n+1)) ifTrue:[     "/ if init_genrand() has not been called,
       
   372             self initGenRand:5489.      "/ a default initial seed is used
       
   373         ].
       
   374 
       
   375         kk := 0.
       
   376         [kk < (n-m)] whileTrue:[
       
   377             y := ((mt at:kk+1) bitAnd:upperMask)
       
   378                     bitOr:((mt at:(kk+1+1)) bitAnd:lowerMask).
       
   379             mt at:(kk+1) put:(((mt at:(kk+m+1)) bitXor:(y>>1)) bitXor: (mag01 at:(y bitAnd:1)+1)).
       
   380             kk := kk + 1.
       
   381         ].
       
   382         [kk < (n-1)] whileTrue:[
       
   383             y := ((mt at:(kk+1)) bitAnd:upperMask)
       
   384                     bitOr:((mt at:(kk+1+1)) bitAnd:lowerMask).
       
   385             mt at:(kk+1) put:(((mt at:(kk+(m-n)+1)) bitXor:(y>>1)) bitXor:(mag01 at:(y bitAnd:1)+1)).
       
   386             kk := kk + 1.
       
   387         ].
       
   388 
       
   389         y := ((mt at:(n-1+1)) bitAnd:upperMask) bitOr:((mt at:1) bitAnd:lowerMask).
       
   390         mt at:n-1+1 put:(((mt at:m) bitXor:(y>>1)) bitXor: (mag01 at:(y bitAnd:1)+1)).
       
   391 
       
   392         mti := 0.
       
   393     ].
       
   394 
       
   395     y := mt at:(mti+1). mti := mti + 1.
       
   396 
       
   397     "/ Tempering 
       
   398     y := y bitXor:(y >> 11).
       
   399     y := y bitXor:((y << 7) bitAnd: 16r9d2c5680).
       
   400     y := y bitXor:((y << 15) bitAnd: 16refc60000).
       
   401     y := y bitXor:(y >> 18).
       
   402 
       
   403     ^ y
       
   404 ! !
       
   405 
       
   406 !RandomMT19937 class methodsFor:'documentation'!
       
   407 
       
   408 version
       
   409     ^ '$Header: /cvs/stx/stx/libbasic2/RandomMT19937.st,v 1.1 2014-10-01 12:43:46 cg Exp $'
       
   410 !
       
   411 
       
   412 version_CVS
       
   413     ^ '$Header: /cvs/stx/stx/libbasic2/RandomMT19937.st,v 1.1 2014-10-01 12:43:46 cg Exp $'
       
   414 ! !
       
   415