1 //
2 //* \file  rng.cc
3 //* \brief Random number generator.
4 //
5 // This file is part of the libecc package.
6 // Copyright (C) 2002, by
7 //
8 // Carlo Wood, Run on IRC <carlo@alinoe.com>
9 // RSA-1024 0x624ACAD5 1997-01-26                    Sign & Encrypt
10 // Fingerprint16 = 32 EC A7 B6 AC DB 65 A6  F6 F6 55 DD 1C DC FF 61
11 //
12 // This program is free software; you can redistribute it and/or
13 // modify it under the terms of the GNU General Public License
14 // as published by the Free Software Foundation; either version 2
15 // of the License, or (at your option) any later version.
16 //
17 // This program is distributed in the hope that it will be useful,
18 // but WITHOUT ANY WARRANTY; without even the implied warranty of
19 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
20 // GNU General Public License for more details.
21 //
22 // You should have received a copy of the GNU General Public License
23 // along with this program; if not, write to the Free Software
24 // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
25 //
26 
27 #include "sys.h"
28 #include "debug.h"
29 #include <libecc/rng.h>
30 
31 namespace libecc {
32 
33 /**\class rng rng.h libecc/rng.h
34  *
35  * \brief Pseudo Random Number Generator.
36  *
37  * This random number generator was designed from scratch.
38  * It is fast (2.5 million bit/second on a 900 Mhz PC),
39  * has a long period (2^521 - 1) and passes all known statistical
40  * tests for randomness.&nbsp;
41  *
42  * \note
43  * \htmlonly<FONT size=-1>\endhtmlonly
44  * The chosen period of 2^521 - 1 is by far large enough for any purpose;
45  * when 15 billion people each using 10,000 PCs with a clock frequency of 100 Ghz
46  * generate random numbers for the next 10 million years - then the chance that \e any two of
47  * them ever generated a (partly) overlapping sequence of bits is less than
48  * 0.000...(108&nbsp;zeroes)...0001.&nbsp;
49  * In order words, if you chose a good random seed, you \e will be
50  * the only one who ever uses the resulting sequence of output bits.
51  * \htmlonly</FONT>\endhtmlonly
52  *
53  * <H4>Theory</H4>
54  *
55  * The basis of this Random Number Generator (RNG) is a shift register of 521 bit using nine feedback points.&nbsp;
56  * Because such a feedback is linear, it is possible to
57  * write any amount of shift as <TT>&nbsp;NEWSTATE = MATRIX * OLDSTATE</TT>.&nbsp;
58  * The matrix that corresponds with a single shift can easily be determined
59  * as function of the feedback points.&nbsp;
60  * Let's call this matrix <TT>&nbsp;M0</TT>.&nbsp;
61  * The matrix corresponding to a shift of 2 then corresponds to <TT>&nbsp;M1 = M0 * M0</TT>.&nbsp;
62  * A shift of 4 corresponds to <TT>&nbsp;M2 = M1 * M1</TT>, etcetera.
63  *
64  * The period of the RNG will be equal to the number of different internal states that will
65  * be used.&nbsp; In the optimal case all possible internal states, except all zeroes which would
66  * lead to an output of only zeroes, will be used and the period of the RNG
67  * will be equal to 2<SUP>521</SUP> - 1.&nbsp;
68  * This hypothesis can be checked by confirming that <TT>M521 = M0<SUP>2<SUP>521</SUP></SUP> == M0&nbsp;</TT>
69  * and because 2<SUP>521</SUP> - 1 is a prime this value then must be the real period of the RNG.&nbsp;
70  * Of course, one of the feedback points is fixed at bit 520 (feeding to bit 0).&nbsp;
71  * By choosing the other feedback points close to bit 0, a very fast mangling of the bits is achieved.
72  *
73  * \note
74  * \htmlonly<FONT size=-1>\endhtmlonly
75  * Even with a seed of '1' (a single bit) the first resulting 512 bits of output appear to be totally random already.&nbsp;
76  * After 9 calls to generate_512_bits(), each bit in the pool has been set at least once.&nbsp;
77  * \htmlonly</FONT>\endhtmlonly
78  *
79  * <H5>Verifying the period, an example</H5>
80  *
81  * Suppose we have a shift register of 3 bits (2<SUP>3</SUP>-1 = 7 is a prime).&nbsp;
82  * Let's use one feedback point from bit 1 (and from bit 2 of course), then we have:
83  *
84  * <PRE>
85  * bit        0    1    2(output)
86  * state0:    a    b    c
87  * state1:    b    c   a^b
88  * state2:    c   a^b  b^c
89  *   ...
90  * </PRE>
91  *
92  * Which gives, with a starting seed of 001:
93  *
94  * <PRE>
95  * 001
96  * 010
97  * 101
98  * 011
99  * 111
100  * 110
101  * 100
102  * and back to the top
103  * </PRE>
104  *
105  * Or rotated 90 degrees:
106  *
107  * <PRE>
108  * bit
109  *  0  0010111
110  *  1  0101110
111  *  2  1011100
112  * </PRE>
113  *
114  * Each column can be considered as a vector:
115  *
116  * <PRE>
117  * a
118  * b
119  * c
120  * </PRE>
121  *
122  * The matrix representing one shift now is:
123  *
124  * <PRE>
125  * m11 m12 m13  a   a m11 + b m12 + c m13    b
126  * m21 m22 m23  b = a m21 + b m22 + c m23 =  c
127  * m31 m32 m33  c   a m31 + b m32 + c m33    a+b mod 2 = a^b
128  * </PRE>
129  *
130  * From which we can deduce the matrix
131  *
132  * <PRE>
133  *      0 1 0
134  * M0 = 0 0 1
135  *      1 1 0
136  * </PRE>
137  *
138  * and that we have to do our calculation modulo 2.
139  *
140  * Note that
141  *
142  * <PRE>
143  *                0 1 0   0 1 0   0 0 1
144  * M1 = M0 * M0 = 0 0 1   0 0 1 = 1 1 0
145  *                1 1 0   1 1 0   0 1 1
146  *
147  *                0 0 1   0 0 1   0 1 1
148  * M2 = M1 * M1 = 1 1 0   1 1 0 = 1 1 1
149  *                0 1 1   0 1 1   1 0 1
150  * </PRE>
151  *
152  * and
153  *
154  * <PRE>
155  *                0 1 1   0 1 1   0 1 0
156  * M3 = M2 * M2 = 1 1 1   1 1 1 = 0 0 1 = M0
157  *                1 0 1   1 0 1   1 1 0
158  *
159  * </PRE>
160  *
161  * Which proves that this repeats after 2<SUP>3</SUP>-1 = 7 times.
162  *
163  * <H5>Predictability</H5>
164  *
165  * Because the state of the RNG changes in a linear way (we can use matrices),
166  * the output is predictable: the internal state can be calculated from a
167  * small amount of output.
168  *
169  * In the above example assume that the RNG produces the output sequence xyz.&nbsp;
170  * If the orginal state was abc then after three bits of output the state would
171  * be:
172  *
173  * <PRE>
174  *       a   0 1 0   0 0 1   a   1 1 0   a   a+b
175  * M0 M1 b = 0 0 1   1 1 0   b = 0 1 1   b = b+c      (mod 2)
176  *       c   1 1 0   0 1 1   c   1 1 1   c   a+b+c
177  * </PRE>
178  *
179  * while the output would be
180  *
181  * <PRE>
182  * x   0 0 1   <SUP> </SUP> a     0 0 0   <SUP> </SUP> a     0 0 0   <SUP> </SUP> a   1 1 0  a
183  * y = 0 0 0 M0<SUP>1</SUP> b  +  0 0 1 M0<SUP>2</SUP> b  +  0 0 0 M0<SUP>3</SUP> b = 0 1 1  b
184  * z   0 0 0   <SUP> </SUP> c     0 0 0   <SUP> </SUP> c     0 0 1   <SUP> </SUP> c   1 1 1  c
185  * </PRE>
186  *
187  * The inverse of the latter allows one to crack the RNG (calculate the current state from x,y,z).
188  *
189  * <H5>Matrix compression</H5>
190  *
191  * Although it might sound easy to check if the period of a given
192  * RNG is maximal, in practise it still isn't when we want to use
193  * a huge number of bits.&nbsp;
194  * I wrote a program that does 4 Gigabit operations per second,
195  * it was therefore able to do 3 times 127 matrix multiplications
196  * of 127x127 matrices per second.&nbsp;
197  * That means that it is reasonably fast to look for working
198  * feedback points.&nbsp; Let's say, one minute was enough.&nbsp;
199  * But if then we try to use this program to find feedback points
200  * for a shift register of 1279 bits, it suddenly takes (1279/127)<SUP>4</SUP> = 10286 minutes = 1 week!
201  * And 1279 bits is <EM>still</EM> not really big...&nbsp; Moreover, storing large matrices
202  * uses a lot of memory.&nbsp; A matrix of 19937 bits (2<SUP>19937</SUP> - 1 is the period
203  * of the <A HREF="http://www.math.keio.ac.jp/~matumoto/emt.html">Mersenne Twister</A>) would
204  * use 47.4 Mb per matrix already.
205  *
206  * Fortunately, the matrices that we need to prove the period of are not arbitrary.&nbsp;
207  * As can be deduced from the paragraphs above, M0 can be written as
208  *
209  * <PRE>
210  *      0 1 0 0 0 0 0 ... 0 0
211  *      0 0 1 0 0 0 0 ... 0 0
212  *      0 0 0 1 0 0 0 ... 0 0
213  *            .&nbsp;
214  * M0 =       .&nbsp;
215  *            .&nbsp;
216  *      0 0 0 0 0 0 0 ... 1 0
217  *      0 0 0 0 0 0 0 ... 0 1
218  *      1 0 0 ... 1 ... 1 ...
219  * </PRE>
220  *
221  * Where the 1's on the bottom row are at the feedback point places.
222  *
223  * As a result, we can write M0 (for a pxp matrix) as
224  * a vector of rows as follows
225  *
226  * <PRE>
227  *       (0 1 0 0 0 0 0 ... 0)
228  *       (0 1 0 0 0 0 0 ... 0) M0<SUP>1</SUP>
229  *       (0 1 0 0 0 0 0 ... 0) M0<SUP>2</SUP>
230  *       (0 1 0 0 0 0 0 ... 0) M0<SUP>3</SUP>
231  * M0 =       .&nbsp;
232  *            .&nbsp;
233  *            .&nbsp;
234  *       (0 1 0 0 0 0 0 ... 0) M0<SUP>p-1</SUP>
235  * </PRE>
236  *
237  * If we multiply this an arbitrary number of times with M0 (say k times),
238  * then we get:
239  *
240  * <PRE>
241  *      <SUP> </SUP>   (0 1 0 0 0 0 0 ... 0) M0<SUP>k&nbsp;</SUP>           (0 1 0 0 0 0 0 ... 0) M0<SUP>k</SUP>
242  *      <SUP> </SUP>   (0 1 0 0 0 0 0 ... 0) M0<SUP>1</SUP> M0<SUP>k</SUP>        (0 1 0 0 0 0 0 ... 0) M0<SUP>k</SUP> M0<SUP>1</SUP>
243  *      <SUP> </SUP>   (0 1 0 0 0 0 0 ... 0) M0<SUP>2</SUP> M0<SUP>k</SUP>        (0 1 0 0 0 0 0 ... 0) M0<SUP>k</SUP> M0<SUP>2</SUP>
244  *      <SUP> </SUP>   (0 1 0 0 0 0 0 ... 0) M0<SUP>3</SUP> M0<SUP>k</SUP>        (0 1 0 0 0 0 0 ... 0) M0<SUP>k</SUP> M0<SUP>3</SUP>
245  * M0 M0<SUP>k</SUP> =       .&nbsp;                        =          .&nbsp;
246  *      <SUP> </SUP>         .&nbsp;                                   .&nbsp;
247  *      <SUP> </SUP>         .&nbsp;                                   .&nbsp;
248  *      <SUP> </SUP>   (0 1 0 0 0 0 0 ... 0) M0<SUP>p-1</SUP> M0<SUP>k</SUP>       (0 1 0 0 0 0 0 ... 0) M0<SUP>k</SUP> M0<SUP>p-1</SUP>
249  * </PRE>
250  *
251  * which proves that every matrix that is a power of M0 can be expressed in terms of the top row
252  * of the matrix and M0 as follows:
253  *
254  * <PRE>
255  * Row(N+1) = Row(N) M0
256  * </PRE>
257  *
258  * <b>Example</b>
259  *
260  * Consider a RNG of 7 bits with a single feedback points at bit 7 - 3 = 4.&nbsp;
261  * The first 7 powers of M0 for that RNG are
262  *
263  * <PRE>
264  *   <SUP> </SUP>   1 0 0 0 0 0 0
265  *   <SUP> </SUP>   <B>0 1 0 0 0 0 0</B>
266  *   <SUP> </SUP>   0 0 1 0 0 0 0
267  * M0<SUP>0</SUP> = 0 0 0 1 0 0 0
268  *   <SUP> </SUP>   0 0 0 0 1 0 0
269  *   <SUP> </SUP>   0 0 0 0 0 1 0
270  *   <SUP> </SUP>   0 0 0 0 0 0 1
271  *
272  *   <SUP> </SUP>   0 1 0 0 0 0 0       <SUP> </SUP>   0 0 1 0 0 0 0       <SUP> </SUP>   0 0 0 1 0 0 0
273  *   <SUP> </SUP>   <B>0 0 1 0 0 0 0</B>       <SUP> </SUP>   <B>0 0 0 1 0 0 0</B>       <SUP> </SUP>   <B>0 0 0 0 1 0 0</B>
274  *   <SUP> </SUP>   0 0 0 1 0 0 0       <SUP> </SUP>   0 0 0 0 1 0 0       <SUP> </SUP>   0 0 0 0 0 1 0
275  * M0<SUP>1</SUP> = 0 0 0 0 1 0 0  ,  M0<SUP>2</SUP> = 0 0 0 0 0 1 0  ,  M0<SUP>3</SUP> = 0 0 0 0 0 0 1
276  *   <SUP> </SUP>   0 0 0 0 0 1 0       <SUP> </SUP>   0 0 0 0 0 0 1       <SUP> </SUP>   1 0 0 0 1 0 0
277  *   <SUP> </SUP>   0 0 0 0 0 0 1       <SUP> </SUP>   1 0 0 0 1 0 0       <SUP> </SUP>   0 1 0 0 0 1 0
278  *   <SUP> </SUP>   1 0 0 0 1 0 0       <SUP> </SUP>   0 1 0 0 0 1 0       <SUP> </SUP>   0 0 1 0 0 0 1
279  *
280  *   <SUP> </SUP>   0 0 0 0 1 0 0      <SUP> </SUP>   0 0 0 0 0 1 0       <SUP> </SUP>   0 0 0 0 0 0 1
281  *   <SUP> </SUP>   <B>0 0 0 0 0 1 0</B>      <SUP> </SUP>   <B>0 0 0 0 0 0 1</B>       <SUP> </SUP>   <B>1 0 0 0 1 0 0</B>
282  *   <SUP> </SUP>   0 0 0 0 0 0 1      <SUP> </SUP>   1 0 0 0 1 0 0       <SUP> </SUP>   0 1 0 0 0 1 0
283  * M0<SUP>4</SUP> = 1 0 0 0 1 0 0 ,  M0<SUP>5</SUP> = 0 1 0 0 0 1 0  ,  M0<SUP>6</SUP> = 0 0 1 0 0 0 1
284  *   <SUP> </SUP>   0 1 0 0 0 1 0      <SUP> </SUP>   0 0 1 0 0 0 1       <SUP> </SUP>   1 0 0 1 1 0 0
285  *   <SUP> </SUP>   0 0 1 0 0 0 1      <SUP> </SUP>   1 0 0 1 1 0 0       <SUP> </SUP>   0 1 0 0 1 1 0
286  *   <SUP> </SUP>   1 0 0 1 1 0 0      <SUP> </SUP>   0 1 0 0 1 1 0       <SUP> </SUP>   0 0 1 0 0 1 1
287  * </PRE>
288  *
289  * And indeed, each second row of a sequential power of M0 is equivalent to a next row of M0 itself.
290  *
291  * Now let's consider an arbitrary power k of M0 (in the example below k = 48).
292  *
293  * <PRE>
294  *   <SUP> </SUP>   <B>0 0 1 1 0 1 0</B>
295  *   <SUP> </SUP>   0 0 0 1 1 0 1
296  *   <SUP> </SUP>   1 0 0 0 0 1 0
297  * M0<SUP>k</SUP> = 0 1 0 0 0 0 1
298  *   <SUP> </SUP>   1 0 1 0 1 0 0
299  *   <SUP> </SUP>   0 1 0 1 0 1 0
300  *   <SUP> </SUP>   0 0 1 0 1 0 1
301  * </PRE>
302  *
303  * and let's call the first row C<SUB>k</SUB><SUP>T</SUP>
304  *
305  * <PRE>
306  *  <SUB> </SUB>   0
307  *  <SUB> </SUB>   0
308  *  <SUB> </SUB>   1
309  * C<SUB>k</SUB> = 1
310  *  <SUB> </SUB>   0
311  *  <SUB> </SUB>   1
312  *  <SUB> </SUB>   0
313  * </PRE>
314  *
315  * then
316  *
317  * <PRE>
318  * M0<SUP>k</SUP> = C<SUB>k</SUB><SUP>T</SUP> P
319  * </PRE>
320  *
321  * where <TT>P</TT> is the vector of matrices
322  *
323  * <PRE>
324  *     M0<SUP>0</SUP>
325  *     M0<SUP>1</SUP>
326  *     M0<SUP>2</SUP>
327  * P = M0<SUP>3</SUP>
328  *     M0<SUP>4</SUP>
329  *     M0<SUP>5</SUP>
330  *     M0<SUP>6</SUP>
331  * </PRE>
332  *
333  * <H5>Multiplying compressed matrices</H5>
334  *
335  * Let <TT>I<SUB>n</SUB></TT> be the n-th column of <TT>I</TT>
336  * and <TT>P<SUB>n</SUB> = P I<SUB>n</SUB></TT>, a vector of the n-th column of each of the matrices of <TT>P</TT>.
337  * Then the n-th column of M0<SUP>k</SUP> is
338  *
339  * <PRE>
340  * M0<SUP>k</SUP> I<SUB>n</SUB> = C<SUB>k</SUB><SUP>T</SUP> P I<SUB>n</SUB> = C<SUB>k</SUB><SUP>T</SUP> P<SUB>n</SUB>
341  * </PRE>
342  *
343  * or
344  *
345  * <PRE>
346  * I<SUB>n</SUB><SUP>T</SUP> (M0<SUP>k</SUP>)<SUP>T</SUP> = P<SUB>n</SUB><SUP>T</SUP> C<SUB>k</SUB>
347  * </PRE>
348  *
349  * Recall that <TT>C<SUB>k</SUB><SUP>T</SUP></TT> is the top row of <TT>M0<SUP>k</SUP></TT>.&nbsp;
350  * That means we can write
351  *
352  * <PRE>
353  * C<SUB>k</SUB> = (I<SUB>0</SUB><SUP>T</SUP> M0<SUP>k</SUP>)<SUP>T</SUP> = (M0<SUP>k</SUP>)<SUP>T</SUP> I<SUB>0</SUB>
354  * </PRE>
355  *
356  * and can deduce that
357  *
358  * <PRE>
359  * C<SUB>m+k</SUB> = (M0<SUP>m+k</SUP>)<SUP>T</SUP> I<SUB>0</SUB> = (M0<SUP>m</SUP> M0<SUP>k</SUP>)<SUP>T</SUP> I<SUB>0</SUB> = (M0<SUP>k</SUP>)<SUP>T</SUP> (M0<SUP>m</SUP>)<SUP>T</SUP> I<SUB>0</SUB> = (M0<SUP>k</SUP>)<SUP>T</SUP> (I<SUB>0</SUB><SUP>T</SUP> M0<SUP>m</SUP>)<SUP>T</SUP> = (M0<SUP>k</SUP>)<SUP>T</SUP> C<SUB>m</SUB> =
360  * [sum over n] (I<SUB>n</SUB><SUP>T</SUP> (M0<SUP>k</SUP>)<SUP>T</SUP> C<SUB>m</SUB> I<SUB>n</SUB>) =
361  * [sum over n] (P<SUB>n</SUB><SUP>T</SUP> C<SUB>k</SUB> C<SUB>m</SUB> I<SUB>n</SUB>) =
362  * [sum over n] (I<SUB>n</SUB><SUP>T</SUP> P<SUP>T</SUP> C<SUB>k</SUB> C<SUB>m</SUB> I<SUB>n</SUB>) =
363  * P<SUP>T</SUP> C<SUB>k</SUB> C<SUB>m</SUB>
364  * </PRE>
365  *
366  * from which we conclude that bit n of <TT>C<SUB>k+m</SUB></TT> is
367  *
368  * <PRE>
369  * I<SUB>n</SUB><SUP>T</SUP> C<SUB>k+m</SUB> =
370  * I<SUB>n</SUB><SUP>T</SUP> P<SUP>T</SUP> C<SUB>k</SUB> C<SUB>m</SUB> =
371  * P<SUB>n</SUB><SUP>T</SUP> C<SUB>k</SUB> C<SUB>m</SUB>
372  * </PRE>
373  *
374  * Note that <TT>P<SUB>n</SUB><SUP>T</SUP></TT> is a row of rows, hence
375  * <TT>P<SUB>n</SUB><SUP>T</SUP> C<SUB>k</SUB></TT> is a row
376  * which can also be written as <TT>C<SUB>k</SUB><SUP>T</SUP> Q</TT>.&nbsp;
377  * The n-th bit of <TT>C<SUB>k+m</SUB></TT> then becomes
378  *
379  * <PRE>
380  * I<SUB>n</SUB><SUP>T</SUP> C<SUB>k+m</SUB> = C<SUB>k</SUB><SUP>T</SUP> Q<SUB>n</SUB> C<SUB>m</SUB>
381  * </PRE>
382  *
383  * where <TT>Q<SUB>n</SUB></TT> is a matrix consisting of the n-th columns of the elements of <TT>P</TT>,
384  * in other words
385  *
386  * <PRE>
387  * Q<SUB>n</SUB> = Q I<SUB>n</SUB>
388  * </PRE>
389  *
390  * where
391  *
392  * <PRE>
393  * Q = (M0<SUP>0</SUP> M0<SUP>1</SUP> M0<SUP>2</SUP> ... M0<SUP>p-1</SUP>)
394  * </PRE>
395  *
396  * Note that Q<SUB>n</SUB> is symmetric so that <TT>Q<SUB>n</SUB><SUP>T</SUP> = Q<SUB>n</SUB></TT> and therefore
397  * <TT>C<SUB>k</SUB><SUP>T</SUP> Q<SUB>n</SUB> C<SUB>m</SUB> = C<SUB>m</SUB><SUP>T</SUP> Q<SUB>n</SUB> C<SUB>k</SUB></TT>,
398  * as it should be since <TT>M0<SUP>k</SUP> M0<SUP>m</SUP> = M0<SUP>m</SUP> M0<SUP>k</SUP></TT>.
399  *
400  * Using this method it was possible to write a program that checks the period of the RNG
401  * in a time of the order O(p<SUP>2</SUP> * f<SUP>3</SUP>), for pxp matrices with f feedback points,
402  * instead of the order O(p<SUP>4</SUP>).  Note that the following observations have been used
403  * as well; the period of a RNG with feedback points f1, f2, f3 ... is equal to the period
404  * of a RNG with feedback points p - f1, p - f2, p - f3, ...  This allowed us to drastically reduce
405  * the order of the number of feedback points.  Finally note the observation that a RNG only has a
406  * period of 2<SUP>p</SUP> - 1 solution when using an odd number of feedback points.
407  *
408  * Considering that it took me something of the order of an hour to find working feedback points
409  * for a RNG with 521 bits, it would now take about 9 weeks to find feedback points for a RNG
410  * of 19937 bits instead of 262 year as would be the case when using ordinairy matrices.
411  *
412  * <b>Implementation</b>
413  *
414  * The implemented Random Number Generator of libecc uses 521 bits and 9 feedback points
415  * at 2, 3, 7, 13, 31, 61, 131, 151 and 251.&nbsp; These feedback points are chosen to be primes
416  * in order to garantee the least interference.&nbsp; The distance between the feedback points
417  * is every time increased by a factor of two (except for the feedback point at 151 which was
418  * added in order to make the RNG have its maximum period).&nbsp; The reason for this is again
419  * to have the least interference between feedback frequencies, this way the different feedback
420  * points nicely supplement each other in achieving bit mangling over the full range.
421  */
422 
423 static unsigned int const feed1 = rng::S_pool_size - 2;
424 static unsigned int const feed2 = rng::S_pool_size - 3;
425 static unsigned int const feed3 = rng::S_pool_size - 7;
426 static unsigned int const feed4 = rng::S_pool_size - 13;
427 static unsigned int const feed5 = rng::S_pool_size - 31;
428 static unsigned int const feed6 = rng::S_pool_size - 61;
429 static unsigned int const feed7 = rng::S_pool_size - 131;
430 static unsigned int const feed8 = rng::S_pool_size - 151;
431 static unsigned int const feed9 = rng::S_pool_size - 251;
432 
433 /**
434  * \brief Constructor.
435  *
436  * \param seed A bitset of 521 bits, any value is allowed except all zeroes.
437  */
rng(pool_type const & seed)438 rng::rng(pool_type const& seed) : M_out_cnt(0), M_entropy_ptr(M_pool), M_entropy_ptr_end(M_pool + sizeof(M_pool) / 4 - 1),
439     M_head(M_pool, 0), M_fbp1(M_pool, feed1), M_fbp2(M_pool, feed2), M_fbp3(M_pool, feed3), M_fbp4(M_pool, feed4),
440     M_fbp5(M_pool, feed5), M_fbp6(M_pool, feed6), M_fbp7(M_pool, feed7), M_fbp8(M_pool, feed8), M_fbp9(M_pool, feed9)
441 {
442   for (unsigned int d = 0; d < sizeof(M_pool) / 4; ++d)
443     M_pool[d] = seed.digit32(d);
444 }
445 
446 /** \brief Generate 512 bits.
447  *
448  * Fills the internal output buffer with 512 new random bits.
449  * You can retrieve this output with the member function rng::get_512_bits().
450  */
generate_512_bits(void)451 void rng::generate_512_bits(void)
452 {
453   do
454   {
455     bitset_digit_t c = 0;
456     for (bitset_digit_t out_mask = 1; out_mask != 0; out_mask <<= 1)
457     {
458       uint32_t a = M_head.increment_and_test(M_pool);
459       a ^= M_fbp1.increment_and_test(M_pool);
460       a ^= M_fbp2.increment_and_test(M_pool);
461       a ^= M_fbp3.increment_and_test(M_pool);
462       a ^= M_fbp4.increment_and_test(M_pool);
463       a ^= M_fbp5.increment_and_test(M_pool);
464       a ^= M_fbp6.increment_and_test(M_pool);
465       a ^= M_fbp7.increment_and_test(M_pool);
466       a ^= M_fbp8.increment_and_test(M_pool);
467       a ^= M_fbp9.increment_and_test(M_pool);
468       uint32_t b = a >> 16;
469       a ^= b;
470       b = a >> 8;
471       a ^= b;
472       if (libecc::oddnumberofbits[a & 0xff])
473       {
474 	M_head.set();
475 	c |= out_mask;
476       }
477       else
478 	M_head.clear();
479     }
480     M_out.rawdigit(M_out_cnt++) = c;
481     M_out_cnt %= (512 / bitset_digit_bits);
482   }
483   while(M_out_cnt != 0);
484 }
485 
486 /** \brief Add entropy to the random number generator pool.
487  *
488  * By adding entropy to the random number generator, the output
489  * becomes more unpredictable.  Its not really necessary to do
490  * this however, its hard enough to 'guess' the 521 bits of
491  * seed.  It is advisable to use SHA-1 on the output though,
492  * in order to make it harder to reverse engineer the internal
493  * state of the rng from its output.
494  */
add_entropy(uint32_t const * noise,unsigned int number_of_ints)495 void rng::add_entropy(uint32_t const* noise, unsigned int number_of_ints)
496 {
497   for (unsigned int cnt = 0; cnt < number_of_ints; ++cnt)
498   {
499     *M_entropy_ptr ^= noise[cnt];
500     if (++M_entropy_ptr == M_entropy_ptr_end)
501       M_entropy_ptr = M_pool;
502   }
503 }
504 
505 unsigned int const rng::S_pool_size;
506 
507 } // namespace libecc
508