Commit | Line | Data |
---|---|---|
78ed81a3 | 1 | // This may look like C code, but it is really -*- C++ -*- |
2 | /* | |
3 | Copyright (C) 1989 Free Software Foundation | |
4 | ||
5 | This file is part of the GNU C++ Library. This library is free | |
6 | software; you can redistribute it and/or modify it under the terms of | |
7 | the GNU Library General Public License as published by the Free | |
8 | Software Foundation; either version 2 of the License, or (at your | |
9 | option) any later version. This library is distributed in the hope | |
10 | that it will be useful, but WITHOUT ANY WARRANTY; without even the | |
11 | implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR | |
12 | PURPOSE. See the GNU Library General Public License for more details. | |
13 | You should have received a copy of the GNU Library General Public | |
14 | License along with this library; if not, write to the Free Software | |
15 | Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. | |
16 | */ | |
15637ed4 RG |
17 | #ifdef __GNUG__ |
18 | #pragma implementation | |
19 | #endif | |
20 | #include <values.h> | |
21 | #include <assert.h> | |
22 | #include <builtin.h> | |
23 | #include <RNG.h> | |
24 | ||
78ed81a3 | 25 | // These two static fields get initialized by RNG::RNG(). |
26 | PrivateRNGSingleType RNG::singleMantissa; | |
27 | PrivateRNGDoubleType RNG::doubleMantissa; | |
28 | ||
15637ed4 RG |
29 | // |
30 | // The scale constant is 2^-31. It is used to scale a 31 bit | |
31 | // long to a double. | |
32 | // | |
33 | ||
34 | //static const double randomDoubleScaleConstant = 4.656612873077392578125e-10; | |
35 | //static const float randomFloatScaleConstant = 4.656612873077392578125e-10; | |
36 | ||
37 | static char initialized = 0; | |
38 | ||
39 | RNG::RNG() | |
40 | { | |
41 | if (!initialized) | |
42 | { | |
43 | ||
44 | assert (sizeof(double) == 2 * sizeof(unsigned long)); | |
45 | ||
46 | // | |
47 | // The following is a hack that I attribute to | |
48 | // Andres Nowatzyk at CMU. The intent of the loop | |
49 | // is to form the smallest number 0 <= x < 1.0, | |
50 | // which is then used as a mask for two longwords. | |
51 | // this gives us a fast way way to produce double | |
52 | // precision numbers from longwords. | |
53 | // | |
54 | // I know that this works for IEEE and VAX floating | |
55 | // point representations. | |
56 | // | |
57 | // A further complication is that gnu C will blow | |
58 | // the following loop, unless compiled with -ffloat-store, | |
59 | // because it uses extended representations for some of | |
60 | // of the comparisons. Thus, we have the following hack. | |
61 | // If we could specify #pragma optimize, we wouldn't need this. | |
62 | // | |
63 | ||
64 | PrivateRNGDoubleType t; | |
65 | PrivateRNGSingleType s; | |
66 | ||
67 | #if _IEEE == 1 | |
68 | ||
69 | t.d = 1.5; | |
70 | if ( t.u[1] == 0 ) { // sun word order? | |
71 | t.u[0] = 0x3fffffff; | |
72 | t.u[1] = 0xffffffff; | |
73 | } | |
74 | else { | |
75 | t.u[0] = 0xffffffff; // encore word order? | |
76 | t.u[1] = 0x3fffffff; | |
77 | } | |
78 | ||
79 | s.u = 0x3fffffff; | |
80 | #else | |
81 | volatile double x = 1.0; // volatile needed when fp hardware used, | |
82 | // and has greater precision than memory doubles | |
83 | double y = 0.5; | |
84 | do { // find largest fp-number < 2.0 | |
85 | t.d = x; | |
86 | x += y; | |
87 | y *= 0.5; | |
88 | } while (x != t.d && x < 2.0); | |
89 | ||
90 | volatile float xx = 1.0; // volatile needed when fp hardware used, | |
91 | // and has greater precision than memory floats | |
92 | float yy = 0.5; | |
93 | do { // find largest fp-number < 2.0 | |
94 | s.s = xx; | |
95 | xx += yy; | |
96 | yy *= 0.5; | |
97 | } while (xx != s.s && xx < 2.0); | |
98 | #endif | |
99 | // set doubleMantissa to 1 for each doubleMantissa bit | |
100 | doubleMantissa.d = 1.0; | |
101 | doubleMantissa.u[0] ^= t.u[0]; | |
102 | doubleMantissa.u[1] ^= t.u[1]; | |
103 | ||
104 | // set singleMantissa to 1 for each singleMantissa bit | |
105 | singleMantissa.s = 1.0; | |
106 | singleMantissa.u ^= s.u; | |
107 | ||
108 | initialized = 1; | |
109 | } | |
110 | } | |
78ed81a3 | 111 | |
112 | float RNG::asFloat() | |
113 | { | |
114 | PrivateRNGSingleType result; | |
115 | result.s = 1.0; | |
116 | result.u |= (asLong() & singleMantissa.u); | |
117 | result.s -= 1.0; | |
118 | assert( result.s < 1.0 && result.s >= 0); | |
119 | return( result.s ); | |
120 | } | |
121 | ||
122 | double RNG::asDouble() | |
123 | { | |
124 | PrivateRNGDoubleType result; | |
125 | result.d = 1.0; | |
126 | result.u[0] |= (asLong() & doubleMantissa.u[0]); | |
127 | result.u[1] |= (asLong() & doubleMantissa.u[1]); | |
128 | result.d -= 1.0; | |
129 | assert( result.d < 1.0 && result.d >= 0); | |
130 | return( result.d ); | |
131 | } | |
132 |