1d25e02daSmrg /* Mersenne Twister pseudo-random number generator functions.
2d25e02daSmrg
3d25e02daSmrg THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY. THEY'RE ALMOST
4d25e02daSmrg CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN
5d25e02daSmrg FUTURE GNU MP RELEASES.
6d25e02daSmrg
7d25e02daSmrg Copyright 2002, 2003, 2006 Free Software Foundation, Inc.
8d25e02daSmrg
9d25e02daSmrg This file is part of the GNU MP Library.
10d25e02daSmrg
11d25e02daSmrg The GNU MP Library is free software; you can redistribute it and/or modify
12f81b1c5bSmrg it under the terms of either:
13f81b1c5bSmrg
14f81b1c5bSmrg * the GNU Lesser General Public License as published by the Free
15f81b1c5bSmrg Software Foundation; either version 3 of the License, or (at your
16d25e02daSmrg option) any later version.
17d25e02daSmrg
18f81b1c5bSmrg or
19f81b1c5bSmrg
20f81b1c5bSmrg * the GNU General Public License as published by the Free Software
21f81b1c5bSmrg Foundation; either version 2 of the License, or (at your option) any
22f81b1c5bSmrg later version.
23f81b1c5bSmrg
24f81b1c5bSmrg or both in parallel, as here.
25f81b1c5bSmrg
26d25e02daSmrg The GNU MP Library is distributed in the hope that it will be useful, but
27d25e02daSmrg WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
28f81b1c5bSmrg or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
29f81b1c5bSmrg for more details.
30d25e02daSmrg
31f81b1c5bSmrg You should have received copies of the GNU General Public License and the
32f81b1c5bSmrg GNU Lesser General Public License along with the GNU MP Library. If not,
33f81b1c5bSmrg see https://www.gnu.org/licenses/. */
34d25e02daSmrg
35d25e02daSmrg #include <stdio.h> /* for NULL */
36d25e02daSmrg
37d25e02daSmrg #include "gmp-impl.h"
38d25e02daSmrg #include "randmt.h"
39d25e02daSmrg
40d25e02daSmrg
41d25e02daSmrg /* This code implements the Mersenne Twister pseudorandom number generator
42d25e02daSmrg by Takuji Nishimura and Makoto Matsumoto. The buffer initialization
43d25e02daSmrg function is different in order to permit seeds greater than 2^32-1.
44d25e02daSmrg
45d25e02daSmrg This file contains a special __gmp_randinit_mt_noseed which excludes the
46d25e02daSmrg seeding function from the gmp_randfnptr_t routines. This is for use by
47d25e02daSmrg mpn_random and mpn_random2 on the global random generator. MT seeding
48d25e02daSmrg uses mpz functions, and we don't want mpn routines dragging mpz functions
49d25e02daSmrg into the link. */
50d25e02daSmrg
51d25e02daSmrg
52d25e02daSmrg /* Default seed to use when the generator is not initialized. */
53d25e02daSmrg #define DEFAULT_SEED 5489 /* was 4357 */
54d25e02daSmrg
55d25e02daSmrg /* Tempering masks. */
56d25e02daSmrg #define MASK_1 0x9D2C5680
57d25e02daSmrg #define MASK_2 0xEFC60000
58d25e02daSmrg
59d25e02daSmrg /* Initial state of buffer when initialized with default seed. */
60d25e02daSmrg static const gmp_uint_least32_t default_state[N] =
61d25e02daSmrg {
62d25e02daSmrg 0xD247B233,0x9E5AA8F1,0x0FFA981B,0x9DCB0980,0x74200F2B,0xA576D044,
63d25e02daSmrg 0xE9F05ADF,0x1538BFF5,0x59818BBF,0xCF9E58D8,0x09FCE032,0x6A1C663F,
64d25e02daSmrg 0x5116E78A,0x69B3E0FA,0x6D92D665,0xD0A8BE98,0xF669B734,0x41AC1B68,
65d25e02daSmrg 0x630423F1,0x4B8D6B8A,0xC2C46DD7,0x5680747D,0x43703E8F,0x3B6103D2,
66d25e02daSmrg 0x49E5EB3F,0xCBDAB4C1,0x9C988E23,0x747BEE0B,0x9111E329,0x9F031B5A,
67d25e02daSmrg 0xECCA71B9,0x2AFE4EF8,0x8421C7ED,0xAC89AFF1,0xAED90DF3,0x2DD74F01,
68d25e02daSmrg 0x14906A13,0x75873FA9,0xFF83F877,0x5028A0C9,0x11B4C41D,0x7CAEDBC4,
69d25e02daSmrg 0x8672D0A7,0x48A7C109,0x8320E59F,0xBC0B3D5F,0x75A30886,0xF9E0D128,
70d25e02daSmrg 0x41AF7580,0x239BB94D,0xC67A3C81,0x74EEBD6E,0xBC02B53C,0x727EA449,
71d25e02daSmrg 0x6B8A2806,0x5853B0DA,0xBDE032F4,0xCE234885,0x320D6145,0x48CC053F,
72d25e02daSmrg 0x00DBC4D2,0xD55A2397,0xE1059B6F,0x1C3E05D1,0x09657C64,0xD07CB661,
73d25e02daSmrg 0x6E982E34,0x6DD1D777,0xEDED1071,0xD79DFD65,0xF816DDCE,0xB6FAF1E4,
74d25e02daSmrg 0x1C771074,0x311835BD,0x18F952F7,0xF8F40350,0x4ECED354,0x7C8AC12B,
75d25e02daSmrg 0x31A9994D,0x4FD47747,0xDC227A23,0x6DFAFDDF,0x6796E748,0x0C6F634F,
76d25e02daSmrg 0xF992FA1D,0x4CF670C9,0x067DFD31,0xA7A3E1A5,0x8CD7D9DF,0x972CCB34,
77d25e02daSmrg 0x67C82156,0xD548F6A8,0x045CEC21,0xF3240BFB,0xDEF656A7,0x43DE08C5,
78d25e02daSmrg 0xDAD1F92F,0x3726C56B,0x1409F19A,0x942FD147,0xB926749C,0xADDC31B8,
79d25e02daSmrg 0x53D0D869,0xD1BA52FE,0x6722DF8C,0x22D95A74,0x7DC1B52A,0x1DEC6FD5,
80d25e02daSmrg 0x7262874D,0x0A725DC9,0xE6A8193D,0xA052835A,0xDC9AD928,0xE59EBB90,
81d25e02daSmrg 0x70DBA9FF,0xD612749D,0x5A5A638C,0x6086EC37,0x2A579709,0x1449EA3A,
82d25e02daSmrg 0xBC8E3C06,0x2F900666,0xFBE74FD1,0x6B35B911,0xF8335008,0xEF1E979D,
83d25e02daSmrg 0x738AB29D,0xA2DC0FDC,0x7696305D,0xF5429DAC,0x8C41813B,0x8073E02E,
84d25e02daSmrg 0xBEF83CCD,0x7B50A95A,0x05EE5862,0x00829ECE,0x8CA1958C,0xBE4EA2E2,
85d25e02daSmrg 0x4293BB73,0x656F7B23,0x417316D8,0x4467D7CF,0x2200E63B,0x109050C8,
86d25e02daSmrg 0x814CBE47,0x36B1D4A8,0x36AF9305,0x308327B3,0xEBCD7344,0xA738DE27,
87d25e02daSmrg 0x5A10C399,0x4142371D,0x64A18528,0x0B31E8B2,0x641057B9,0x6AFC363B,
88d25e02daSmrg 0x108AD953,0x9D4DA234,0x0C2D9159,0x1C8A1A1F,0x310C66BA,0x87AA1070,
89d25e02daSmrg 0xDAC832FF,0x0A433422,0x7AF15812,0x2D8D9BD0,0x995A25E9,0x25326CAC,
90d25e02daSmrg 0xA34384DB,0x4C8421CC,0x4F0315EC,0x29E8649E,0xA7732D6F,0x2E94D3E3,
91d25e02daSmrg 0x7D98A340,0x397C4D74,0x659DB4DE,0x747D4E9A,0xD9DB8435,0x4659DBE9,
92d25e02daSmrg 0x313E6DC5,0x29D104DC,0x9F226CBA,0x452F18B0,0xD0BC5068,0x844CA299,
93d25e02daSmrg 0x782B294E,0x4AE2EB7B,0xA4C475F8,0x70A81311,0x4B3E8BCC,0x7E20D4BA,
94d25e02daSmrg 0xABCA33C9,0x57BE2960,0x44F9B419,0x2E567746,0x72EB757A,0x102CC0E8,
95d25e02daSmrg 0xB07F32B9,0xD0DABD59,0xBA85AD6B,0xF3E20667,0x98D77D81,0x197AFA47,
96d25e02daSmrg 0x518EE9AC,0xE10CE5A2,0x01CF2C2A,0xD3A3AF3D,0x16DDFD65,0x669232F8,
97d25e02daSmrg 0x1C50A301,0xB93D9151,0x9354D3F4,0x847D79D0,0xD5FE2EC6,0x1F7B0610,
98d25e02daSmrg 0xFA6B90A5,0xC5879041,0x2E7DC05E,0x423F1F32,0xEF623DDB,0x49C13280,
99d25e02daSmrg 0x98714E92,0xC7B6E4AD,0xC4318466,0x0737F312,0x4D3C003F,0x9ACC1F1F,
100d25e02daSmrg 0x5F1C926D,0x085FA771,0x185A83A2,0xF9AA159D,0x0B0B0132,0xF98E7A43,
101d25e02daSmrg 0xCD9EBDBE,0x0190CB29,0x10D93FB6,0x3B8A4D97,0x66A65A41,0xE43E766F,
102d25e02daSmrg 0x77BE3C41,0xB9686364,0xCB36994D,0x6846A287,0x567E77F7,0x36178DD8,
103d25e02daSmrg 0xBDE6B1F2,0xB6EFDC64,0x82950324,0x42053F47,0xC09BE51C,0x0942D762,
104d25e02daSmrg 0x35F92C7F,0x367DEC61,0x6EE3D983,0xDBAAF78A,0x265D2C47,0x8EB4BF5C,
105d25e02daSmrg 0x33B232D7,0xB0137E77,0x373C39A7,0x8D2B2E76,0xC7510F01,0x50F9E032,
106d25e02daSmrg 0x7B1FDDDB,0x724C2AAE,0xB10ECB31,0xCCA3D1B8,0x7F0BCF10,0x4254BBBD,
107d25e02daSmrg 0xE3F93B97,0x2305039B,0x53120E22,0x1A2F3B9A,0x0FDDBD97,0x0118561E,
108d25e02daSmrg 0x0A798E13,0x9E0B3ACD,0xDB6C9F15,0xF512D0A2,0x9E8C3A28,0xEE2184AE,
109d25e02daSmrg 0x0051EC2F,0x2432F74F,0xB0AA66EA,0x55128D88,0xF7D83A38,0x4DAE8E82,
110d25e02daSmrg 0x3FDC98D6,0x5F0BD341,0x7244BE1D,0xC7B48E78,0x2D473053,0x43892E20,
111d25e02daSmrg 0xBA0F1F2A,0x524D4895,0x2E10BCB1,0x4C372D81,0x5C3E50CD,0xCF61CC2E,
112d25e02daSmrg 0x931709AB,0x81B3AEFC,0x39E9405E,0x7FFE108C,0x4FBB3FF8,0x06ABE450,
113d25e02daSmrg 0x7F5BF51E,0xA4E3CDFD,0xDB0F6C6F,0x159A1227,0x3B9FED55,0xD20B6F7F,
114d25e02daSmrg 0xFBE9CC83,0x64856619,0xBF52B8AF,0x9D7006B0,0x71165BC6,0xAE324AEE,
115d25e02daSmrg 0x29D27F2C,0x794C2086,0x74445CE2,0x782915CC,0xD4CE6886,0x3289AE7C,
116d25e02daSmrg 0x53DEF297,0x4185F7ED,0x88B72400,0x3C09DC11,0xBCE3AAB6,0x6A75934A,
117d25e02daSmrg 0xB267E399,0x000DF1BF,0x193BA5E2,0xFA3E1977,0x179E14F6,0x1EEDE298,
118d25e02daSmrg 0x691F0B06,0xB84F78AC,0xC1C15316,0xFFFF3AD6,0x0B457383,0x518CD612,
119d25e02daSmrg 0x05A00F3E,0xD5B7D275,0x4C5ECCD7,0xE02CD0BE,0x5558E9F2,0x0C89BBF0,
120d25e02daSmrg 0xA3D96227,0x2832D2B2,0xF667B897,0xD4556554,0xF9D2F01F,0xFA1E3FAE,
121d25e02daSmrg 0x52C2E1EE,0xE5451F31,0x7E849729,0xDABDB67A,0x54BF5E7E,0xF831C271,
122d25e02daSmrg 0x5F1A17E3,0x9D140AFE,0x92741C47,0x48CFABCE,0x9CBBE477,0x9C3EE57F,
123d25e02daSmrg 0xB07D4C39,0xCC21BCE2,0x697708B1,0x58DA2A6B,0x2370DB16,0x6E641948,
124d25e02daSmrg 0xACC5BD52,0x868F24CC,0xCA1DB0F5,0x4CADA492,0x3F443E54,0xC4A4D5E9,
125d25e02daSmrg 0xF00AD670,0xE93C86E0,0xFE90651A,0xDDE532A3,0xA66458DF,0xAB7D7151,
126d25e02daSmrg 0x0E2E775F,0xC9109F99,0x8D96D59F,0x73CEF14C,0xC74E88E9,0x02712DC0,
127d25e02daSmrg 0x04F41735,0x2E5914A2,0x59F4B2FB,0x0287FC83,0x80BC0343,0xF6B32559,
128d25e02daSmrg 0xC74178D4,0xF1D99123,0x383CCC07,0xACC0637D,0x0863A548,0xA6FCAC85,
129d25e02daSmrg 0x2A13EFF0,0xAF2EEDB1,0x41E72750,0xE0C6B342,0x5DA22B46,0x635559E0,
130d25e02daSmrg 0xD2EA40AC,0x10AA98C0,0x19096497,0x112C542B,0x2C85040C,0xA868E7D0,
131d25e02daSmrg 0x6E260188,0xF596D390,0xC3BB5D7A,0x7A2AA937,0xDFD15032,0x6780AE3B,
132d25e02daSmrg 0xDB5F9CD8,0x8BD266B0,0x7744AF12,0xB463B1B0,0x589629C9,0xE30DBC6E,
133d25e02daSmrg 0x880F5569,0x209E6E16,0x9DECA50C,0x02987A57,0xBED3EA57,0xD3A678AA,
134d25e02daSmrg 0x70DD030D,0x0CFD9C5D,0x92A18E99,0xF5740619,0x7F6F0A7D,0x134CAF9A,
135d25e02daSmrg 0x70F5BAE4,0x23DCA7B5,0x4D788FCD,0xC7F07847,0xBCF77DA1,0x9071D568,
136d25e02daSmrg 0xFC627EA1,0xAE004B77,0x66B54BCB,0x7EF2DAAC,0xDCD5AC30,0xB9BDF730,
137d25e02daSmrg 0x505A97A7,0x9D881FD3,0xADB796CC,0x94A1D202,0x97535D7F,0x31EC20C0,
138d25e02daSmrg 0xB1887A98,0xC1475069,0xA6F73AF3,0x71E4E067,0x46A569DE,0xD2ADE430,
139d25e02daSmrg 0x6F0762C7,0xF50876F4,0x53510542,0x03741C3E,0x53502224,0xD8E54D60,
140d25e02daSmrg 0x3C44AB1A,0x34972B46,0x74BFA89D,0xD7D768E0,0x37E605DC,0xE13D1BDF,
141d25e02daSmrg 0x5051C421,0xB9E057BE,0xB717A14C,0xA1730C43,0xB99638BE,0xB5D5F36D,
142d25e02daSmrg 0xE960D9EA,0x6B1388D3,0xECB6D3B6,0xBDBE8B83,0x2E29AFC5,0x764D71EC,
143d25e02daSmrg 0x4B8F4F43,0xC21DDC00,0xA63F657F,0x82678130,0xDBF535AC,0xA594FC58,
144d25e02daSmrg 0x942686BC,0xBD9B657B,0x4A0F9B61,0x44FF184F,0x38E10A2F,0x61910626,
145d25e02daSmrg 0x5E247636,0x7106D137,0xC62802F0,0xBD1D1F00,0x7CC0DCB2,0xED634909,
146d25e02daSmrg 0xDC13B24E,0x9799C499,0xD77E3D6A,0x14773B68,0x967A4FB7,0x35EECFB1,
147d25e02daSmrg 0x2A5110B8,0xE2F0AF94,0x9D09DEA5,0x20255D27,0x5771D34B,0xE1089EE4,
148d25e02daSmrg 0x246F330B,0x8F7CAEE5,0xD3064712,0x75CAFBEE,0xB94F7028,0xED953666,
149d25e02daSmrg 0x5D1975B4,0x5AF81271,0x13BE2025,0x85194659,0x30805331,0xEC9D46C0,
150d25e02daSmrg 0xBC027C36,0x2AF84188,0xC2141B80,0xC02B1E4A,0x04D36177,0xFC50E9D7,
151d25e02daSmrg 0x39CE79DA,0x917E0A00,0xEF7A0BF4,0xA98BD8D1,0x19424DD2,0x9439DF1F,
152d25e02daSmrg 0xC42AF746,0xADDBE83E,0x85221F0D,0x45563E90,0x9095EC52,0x77887B25,
153d25e02daSmrg 0x8AE46064,0xBD43B71A,0xBB541956,0x7366CF9D,0xEE8E1737,0xB5A727C9,
154d25e02daSmrg 0x5076B3E7,0xFC70BACA,0xCE135B75,0xC4E91AA3,0xF0341911,0x53430C3F,
155d25e02daSmrg 0x886B0824,0x6BB5B8B7,0x33E21254,0xF193B456,0x5B09617F,0x215FFF50,
156d25e02daSmrg 0x48D97EF1,0x356479AB,0x6EA9DDC4,0x0D352746,0xA2F5CE43,0xB226A1B3,
157d25e02daSmrg 0x1329EA3C,0x7A337CC2,0xB5CCE13D,0x563E3B5B,0x534E8E8F,0x561399C9,
158d25e02daSmrg 0xE1596392,0xB0F03125,0x4586645B,0x1F371847,0x94EAABD1,0x41F97EDD,
159d25e02daSmrg 0xE3E5A39B,0x71C774E2,0x507296F4,0x5960133B,0x7852C494,0x3F5B2691,
160d25e02daSmrg 0xA3F87774,0x5A7AF89E,0x17DA3F28,0xE9D9516D,0xFCC1C1D5,0xE4618628,
161d25e02daSmrg 0x04081047,0xD8E4DB5F,0xDC380416,0x8C4933E2,0x95074D53,0xB1B0032D,
162d25e02daSmrg 0xCC8102EA,0x71641243,0x98D6EB6A,0x90FEC945,0xA0914345,0x6FAB037D,
163d25e02daSmrg 0x70F49C4D,0x05BF5B0E,0x927AAF7F,0xA1940F61,0xFEE0756F,0xF815369F,
164d25e02daSmrg 0x5C00253B,0xF2B9762F,0x4AEB3CCC,0x1069F386,0xFBA4E7B9,0x70332665,
165d25e02daSmrg 0x6BCA810E,0x85AB8058,0xAE4B2B2F,0x9D120712,0xBEE8EACB,0x776A1112
166d25e02daSmrg };
167d25e02daSmrg
168d25e02daSmrg void
__gmp_mt_recalc_buffer(gmp_uint_least32_t mt[])169d25e02daSmrg __gmp_mt_recalc_buffer (gmp_uint_least32_t mt[])
170d25e02daSmrg {
171d25e02daSmrg gmp_uint_least32_t y;
172d25e02daSmrg int kk;
173d25e02daSmrg
174d25e02daSmrg for (kk = 0; kk < N - M; kk++)
175d25e02daSmrg {
176d25e02daSmrg y = (mt[kk] & 0x80000000) | (mt[kk + 1] & 0x7FFFFFFF);
177d25e02daSmrg mt[kk] = mt[kk + M] ^ (y >> 1) ^ ((y & 0x01) != 0 ? MATRIX_A : 0);
178d25e02daSmrg }
179d25e02daSmrg for (; kk < N - 1; kk++)
180d25e02daSmrg {
181d25e02daSmrg y = (mt[kk] & 0x80000000) | (mt[kk + 1] & 0x7FFFFFFF);
182d25e02daSmrg mt[kk] = mt[kk - (N - M)] ^ (y >> 1) ^ ((y & 0x01) != 0 ? MATRIX_A : 0);
183d25e02daSmrg }
184d25e02daSmrg
185d25e02daSmrg y = (mt[N - 1] & 0x80000000) | (mt[0] & 0x7FFFFFFF);
186d25e02daSmrg mt[N - 1] = mt[M - 1] ^ (y >> 1) ^ ((y & 0x01) != 0 ? MATRIX_A : 0);
187d25e02daSmrg }
188d25e02daSmrg
189d25e02daSmrg
190d25e02daSmrg /* Get nbits bits of output from the generator into dest.
191d25e02daSmrg Note that Mersenne Twister is designed to produce outputs in
192d25e02daSmrg 32-bit words. */
193d25e02daSmrg void
__gmp_randget_mt(gmp_randstate_t rstate,mp_ptr dest,unsigned long int nbits)194d25e02daSmrg __gmp_randget_mt (gmp_randstate_t rstate, mp_ptr dest, unsigned long int nbits)
195d25e02daSmrg {
196d25e02daSmrg gmp_uint_least32_t y;
197d25e02daSmrg int rbits;
198d25e02daSmrg mp_size_t i;
199d25e02daSmrg mp_size_t nlimbs;
200d25e02daSmrg int *pmti;
201d25e02daSmrg gmp_uint_least32_t *mt;
202d25e02daSmrg
203d25e02daSmrg pmti = &((gmp_rand_mt_struct *) RNG_STATE (rstate))->mti;
204d25e02daSmrg mt = ((gmp_rand_mt_struct *) RNG_STATE (rstate))->mt;
205d25e02daSmrg
206d25e02daSmrg nlimbs = nbits / GMP_NUMB_BITS;
207d25e02daSmrg rbits = nbits % GMP_NUMB_BITS;
208d25e02daSmrg
209d25e02daSmrg #define NEXT_RANDOM \
210d25e02daSmrg do \
211d25e02daSmrg { \
212d25e02daSmrg if (*pmti >= N) \
213d25e02daSmrg { \
214d25e02daSmrg __gmp_mt_recalc_buffer (mt); \
215d25e02daSmrg *pmti = 0; \
216d25e02daSmrg } \
217d25e02daSmrg y = mt[(*pmti)++]; \
218d25e02daSmrg y ^= (y >> 11); \
219d25e02daSmrg y ^= (y << 7) & MASK_1; \
220d25e02daSmrg y ^= (y << 15) & MASK_2; \
221d25e02daSmrg y ^= (y >> 18); \
222d25e02daSmrg } \
223d25e02daSmrg while (0)
224d25e02daSmrg
225d25e02daSmrg
226d25e02daSmrg /* Handle the common cases of 32- or 64-bit limbs with fast,
227d25e02daSmrg optimized routines, and the rest of cases with a general
228d25e02daSmrg routine. In all cases, no more than 31 bits are rejected
229d25e02daSmrg for the last limb so that every version of the code is
230d25e02daSmrg consistent with the others. */
231d25e02daSmrg
232d25e02daSmrg #if (GMP_NUMB_BITS == 32)
233d25e02daSmrg
234d25e02daSmrg for (i = 0; i < nlimbs; i++)
235d25e02daSmrg {
236d25e02daSmrg NEXT_RANDOM;
237d25e02daSmrg dest[i] = (mp_limb_t) y;
238d25e02daSmrg }
239d25e02daSmrg if (rbits)
240d25e02daSmrg {
241d25e02daSmrg NEXT_RANDOM;
242d25e02daSmrg dest[nlimbs] = (mp_limb_t) (y & ~(ULONG_MAX << rbits));
243d25e02daSmrg }
244d25e02daSmrg
245d25e02daSmrg #else /* GMP_NUMB_BITS != 32 */
246d25e02daSmrg #if (GMP_NUMB_BITS == 64)
247d25e02daSmrg
248d25e02daSmrg for (i = 0; i < nlimbs; i++)
249d25e02daSmrg {
250d25e02daSmrg NEXT_RANDOM;
251d25e02daSmrg dest[i] = (mp_limb_t) y;
252d25e02daSmrg NEXT_RANDOM;
253d25e02daSmrg dest[i] |= (mp_limb_t) y << 32;
254d25e02daSmrg }
255d25e02daSmrg if (rbits)
256d25e02daSmrg {
257d25e02daSmrg if (rbits < 32)
258d25e02daSmrg {
259d25e02daSmrg NEXT_RANDOM;
260d25e02daSmrg dest[nlimbs] = (mp_limb_t) (y & ~(ULONG_MAX << rbits));
261d25e02daSmrg }
262d25e02daSmrg else
263d25e02daSmrg {
264d25e02daSmrg NEXT_RANDOM;
265d25e02daSmrg dest[nlimbs] = (mp_limb_t) y;
266d25e02daSmrg if (rbits > 32)
267d25e02daSmrg {
268d25e02daSmrg NEXT_RANDOM;
269d25e02daSmrg dest[nlimbs] |=
270d25e02daSmrg ((mp_limb_t) (y & ~(ULONG_MAX << (rbits-32)))) << 32;
271d25e02daSmrg }
272d25e02daSmrg }
273d25e02daSmrg }
274d25e02daSmrg
275d25e02daSmrg #else /* GMP_NUMB_BITS != 64 */
276d25e02daSmrg
277d25e02daSmrg {
278d25e02daSmrg /* Fall back to a general algorithm. This algorithm works by
279d25e02daSmrg keeping a pool of up to 64 bits (2 outputs from MT) acting
280d25e02daSmrg as a shift register from which bits are consumed as needed.
281d25e02daSmrg Bits are consumed using the LSB bits of bitpool_l, and
282d25e02daSmrg inserted via bitpool_h and shifted to the right place. */
283d25e02daSmrg
284d25e02daSmrg gmp_uint_least32_t bitpool_h = 0;
285d25e02daSmrg gmp_uint_least32_t bitpool_l = 0;
286d25e02daSmrg int bits_in_pool = 0; /* Holds number of valid bits in the pool. */
287d25e02daSmrg int bits_to_fill; /* Holds total number of bits to put in
288d25e02daSmrg destination. */
289d25e02daSmrg int bitidx; /* Holds the destination bit position. */
290d25e02daSmrg mp_size_t nlimbs2; /* Number of whole+partial limbs to fill. */
291d25e02daSmrg
292d25e02daSmrg nlimbs2 = nlimbs + (rbits != 0);
293d25e02daSmrg
294d25e02daSmrg for (i = 0; i < nlimbs2; i++)
295d25e02daSmrg {
296d25e02daSmrg bitidx = 0;
297d25e02daSmrg if (i < nlimbs)
298d25e02daSmrg bits_to_fill = GMP_NUMB_BITS;
299d25e02daSmrg else
300d25e02daSmrg bits_to_fill = rbits;
301d25e02daSmrg
302d25e02daSmrg dest[i] = CNST_LIMB (0);
303d25e02daSmrg while (bits_to_fill >= 32) /* Process whole 32-bit blocks first. */
304d25e02daSmrg {
305d25e02daSmrg if (bits_in_pool < 32) /* Need more bits. */
306d25e02daSmrg {
307d25e02daSmrg /* 64-bit right shift. */
308d25e02daSmrg NEXT_RANDOM;
309d25e02daSmrg bitpool_h = y;
310d25e02daSmrg bitpool_l |= (bitpool_h << bits_in_pool) & 0xFFFFFFFF;
311d25e02daSmrg if (bits_in_pool == 0)
312d25e02daSmrg bitpool_h = 0;
313d25e02daSmrg else
314d25e02daSmrg bitpool_h >>= 32 - bits_in_pool;
315d25e02daSmrg bits_in_pool += 32; /* We've got 32 more bits. */
316d25e02daSmrg }
317d25e02daSmrg
318d25e02daSmrg /* Fill a 32-bit chunk. */
319d25e02daSmrg dest[i] |= ((mp_limb_t) bitpool_l) << bitidx;
320d25e02daSmrg bitpool_l = bitpool_h;
321d25e02daSmrg bits_in_pool -= 32;
322d25e02daSmrg bits_to_fill -= 32;
323d25e02daSmrg bitidx += 32;
324d25e02daSmrg }
325d25e02daSmrg
326d25e02daSmrg /* Cover the case where GMP_NUMB_BITS is not a multiple of 32. */
327d25e02daSmrg if (bits_to_fill != 0)
328d25e02daSmrg {
329d25e02daSmrg if (bits_in_pool < bits_to_fill)
330d25e02daSmrg {
331d25e02daSmrg NEXT_RANDOM;
332d25e02daSmrg bitpool_h = y;
333d25e02daSmrg bitpool_l |= (bitpool_h << bits_in_pool) & 0xFFFFFFFF;
334d25e02daSmrg if (bits_in_pool == 0)
335d25e02daSmrg bitpool_h = 0;
336d25e02daSmrg else
337d25e02daSmrg bitpool_h >>= 32 - bits_in_pool;
338d25e02daSmrg bits_in_pool += 32;
339d25e02daSmrg }
340d25e02daSmrg
341d25e02daSmrg dest[i] |= (((mp_limb_t) bitpool_l
342d25e02daSmrg & ~(~CNST_LIMB (0) << bits_to_fill))
343d25e02daSmrg << bitidx);
344d25e02daSmrg bitpool_l = ((bitpool_l >> bits_to_fill)
345d25e02daSmrg | (bitpool_h << (32 - bits_to_fill))) & 0xFFFFFFFF;
346d25e02daSmrg bitpool_h >>= bits_to_fill;
347d25e02daSmrg bits_in_pool -= bits_to_fill;
348d25e02daSmrg }
349d25e02daSmrg }
350d25e02daSmrg }
351d25e02daSmrg
352d25e02daSmrg #endif /* GMP_NUMB_BITS != 64 */
353d25e02daSmrg #endif /* GMP_NUMB_BITS != 32 */
354d25e02daSmrg }
355d25e02daSmrg
356d25e02daSmrg void
__gmp_randclear_mt(gmp_randstate_t rstate)357d25e02daSmrg __gmp_randclear_mt (gmp_randstate_t rstate)
358d25e02daSmrg {
359d25e02daSmrg (*__gmp_free_func) ((void *) RNG_STATE (rstate),
360f81b1c5bSmrg ALLOC (rstate->_mp_seed) * GMP_LIMB_BYTES);
361d25e02daSmrg }
362d25e02daSmrg
363d25e02daSmrg void __gmp_randiset_mt (gmp_randstate_ptr, gmp_randstate_srcptr);
364d25e02daSmrg
365d25e02daSmrg static const gmp_randfnptr_t Mersenne_Twister_Generator_Noseed = {
366d25e02daSmrg NULL,
367d25e02daSmrg __gmp_randget_mt,
368d25e02daSmrg __gmp_randclear_mt,
369d25e02daSmrg __gmp_randiset_mt
370d25e02daSmrg };
371d25e02daSmrg
372d25e02daSmrg void
__gmp_randiset_mt(gmp_randstate_ptr dst,gmp_randstate_srcptr src)373d25e02daSmrg __gmp_randiset_mt (gmp_randstate_ptr dst, gmp_randstate_srcptr src)
374d25e02daSmrg {
375f81b1c5bSmrg const mp_size_t sz = ((sizeof (gmp_rand_mt_struct) - 1) / GMP_LIMB_BYTES) + 1;
376d25e02daSmrg gmp_rand_mt_struct *dstp, *srcp;
377d25e02daSmrg mp_size_t i;
378d25e02daSmrg
379d25e02daSmrg /* Set the generator functions. */
380*671ea119Smrg RNG_FNPTR (dst) = RNG_FNPTR(src);
381d25e02daSmrg
382d25e02daSmrg /* Allocate the MT-specific state. */
383d25e02daSmrg dstp = (gmp_rand_mt_struct *) __GMP_ALLOCATE_FUNC_LIMBS (sz);
384d25e02daSmrg RNG_STATE (dst) = (mp_ptr) dstp;
385d25e02daSmrg ALLOC (dst->_mp_seed) = sz; /* Initialize alloc field to placate Camm. */
386d25e02daSmrg
387d25e02daSmrg /* Copy state. */
388d25e02daSmrg srcp = (gmp_rand_mt_struct *) RNG_STATE (src);
389d25e02daSmrg for (i = 0; i < N; i++)
390d25e02daSmrg dstp->mt[i] = srcp->mt[i];
391d25e02daSmrg
392d25e02daSmrg dstp->mti = srcp->mti;
393d25e02daSmrg }
394d25e02daSmrg
395d25e02daSmrg void
__gmp_randinit_mt_noseed(gmp_randstate_ptr dst)396d25e02daSmrg __gmp_randinit_mt_noseed (gmp_randstate_ptr dst)
397d25e02daSmrg {
398f81b1c5bSmrg const mp_size_t sz = ((sizeof (gmp_rand_mt_struct) - 1) / GMP_LIMB_BYTES) + 1;
399d25e02daSmrg gmp_rand_mt_struct *dstp;
400d25e02daSmrg mp_size_t i;
401d25e02daSmrg
402d25e02daSmrg /* Set the generator functions. */
403d25e02daSmrg RNG_FNPTR (dst) = (void *) &Mersenne_Twister_Generator_Noseed;
404d25e02daSmrg
405d25e02daSmrg /* Allocate the MT-specific state. */
406d25e02daSmrg dstp = (gmp_rand_mt_struct *) __GMP_ALLOCATE_FUNC_LIMBS (sz);
407d25e02daSmrg RNG_STATE (dst) = (mp_ptr) dstp;
408d25e02daSmrg ALLOC (dst->_mp_seed) = sz; /* Initialize alloc field to placate Camm. */
409d25e02daSmrg
410d25e02daSmrg /* Set state for default seed. */
411d25e02daSmrg for (i = 0; i < N; i++)
412d25e02daSmrg dstp->mt[i] = default_state[i];
413d25e02daSmrg
414d25e02daSmrg dstp->mti = WARM_UP % N;
415d25e02daSmrg }
416