diff options
Diffstat (limited to 'src/regress/lib/libc/cephes/drand.c')
-rw-r--r-- | src/regress/lib/libc/cephes/drand.c | 174 |
1 files changed, 0 insertions, 174 deletions
diff --git a/src/regress/lib/libc/cephes/drand.c b/src/regress/lib/libc/cephes/drand.c deleted file mode 100644 index 7f7000b4e8..0000000000 --- a/src/regress/lib/libc/cephes/drand.c +++ /dev/null | |||
@@ -1,174 +0,0 @@ | |||
1 | /* $OpenBSD: drand.c,v 1.1 2011/07/02 18:11:01 martynas Exp $ */ | ||
2 | |||
3 | /* | ||
4 | * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> | ||
5 | * | ||
6 | * Permission to use, copy, modify, and distribute this software for any | ||
7 | * purpose with or without fee is hereby granted, provided that the above | ||
8 | * copyright notice and this permission notice appear in all copies. | ||
9 | * | ||
10 | * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES | ||
11 | * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF | ||
12 | * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR | ||
13 | * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES | ||
14 | * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN | ||
15 | * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF | ||
16 | * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. | ||
17 | */ | ||
18 | |||
19 | /* drand.c | ||
20 | * | ||
21 | * Pseudorandom number generator | ||
22 | * | ||
23 | * | ||
24 | * | ||
25 | * SYNOPSIS: | ||
26 | * | ||
27 | * double y, drand(); | ||
28 | * | ||
29 | * drand( &y ); | ||
30 | * | ||
31 | * | ||
32 | * | ||
33 | * DESCRIPTION: | ||
34 | * | ||
35 | * Yields a random number 1.0 <= y < 2.0. | ||
36 | * | ||
37 | * The three-generator congruential algorithm by Brian | ||
38 | * Wichmann and David Hill (BYTE magazine, March, 1987, | ||
39 | * pp 127-8) is used. The period, given by them, is | ||
40 | * 6953607871644. | ||
41 | * | ||
42 | * Versions invoked by the different arithmetic compile | ||
43 | * time options DEC, IBMPC, and MIEEE, produce | ||
44 | * approximately the same sequences, differing only in the | ||
45 | * least significant bits of the numbers. The UNK option | ||
46 | * implements the algorithm as recommended in the BYTE | ||
47 | * article. It may be used on all computers. However, | ||
48 | * the low order bits of a double precision number may | ||
49 | * not be adequately random, and may vary due to arithmetic | ||
50 | * implementation details on different computers. | ||
51 | * | ||
52 | * The other compile options generate an additional random | ||
53 | * integer that overwrites the low order bits of the double | ||
54 | * precision number. This reduces the period by a factor of | ||
55 | * two but tends to overcome the problems mentioned. | ||
56 | * | ||
57 | */ | ||
58 | |||
59 | #include "mconf.h" | ||
60 | |||
61 | |||
62 | /* Three-generator random number algorithm | ||
63 | * of Brian Wichmann and David Hill | ||
64 | * BYTE magazine, March, 1987 pp 127-8 | ||
65 | * | ||
66 | * The period, given by them, is (p-1)(q-1)(r-1)/4 = 6.95e12. | ||
67 | */ | ||
68 | |||
69 | static int sx = 1; | ||
70 | static int sy = 10000; | ||
71 | static int sz = 3000; | ||
72 | |||
73 | static union { | ||
74 | double d; | ||
75 | unsigned short s[4]; | ||
76 | } unkans; | ||
77 | |||
78 | /* This function implements the three | ||
79 | * congruential generators. | ||
80 | */ | ||
81 | |||
82 | static int ranwh() | ||
83 | { | ||
84 | int r, s; | ||
85 | |||
86 | /* sx = sx * 171 mod 30269 */ | ||
87 | r = sx/177; | ||
88 | s = sx - 177 * r; | ||
89 | sx = 171 * s - 2 * r; | ||
90 | if( sx < 0 ) | ||
91 | sx += 30269; | ||
92 | |||
93 | |||
94 | /* sy = sy * 172 mod 30307 */ | ||
95 | r = sy/176; | ||
96 | s = sy - 176 * r; | ||
97 | sy = 172 * s - 35 * r; | ||
98 | if( sy < 0 ) | ||
99 | sy += 30307; | ||
100 | |||
101 | /* sz = 170 * sz mod 30323 */ | ||
102 | r = sz/178; | ||
103 | s = sz - 178 * r; | ||
104 | sz = 170 * s - 63 * r; | ||
105 | if( sz < 0 ) | ||
106 | sz += 30323; | ||
107 | /* The results are in static sx, sy, sz. */ | ||
108 | return 0; | ||
109 | } | ||
110 | |||
111 | /* drand.c | ||
112 | * | ||
113 | * Random double precision floating point number between 1 and 2. | ||
114 | * | ||
115 | * C callable: | ||
116 | * drand( &x ); | ||
117 | */ | ||
118 | |||
119 | int drand( a ) | ||
120 | double *a; | ||
121 | { | ||
122 | unsigned short r; | ||
123 | #ifdef DEC | ||
124 | unsigned short s, t; | ||
125 | #endif | ||
126 | |||
127 | /* This algorithm of Wichmann and Hill computes a floating point | ||
128 | * result: | ||
129 | */ | ||
130 | ranwh(); | ||
131 | unkans.d = sx/30269.0 + sy/30307.0 + sz/30323.0; | ||
132 | r = unkans.d; | ||
133 | unkans.d -= r; | ||
134 | unkans.d += 1.0; | ||
135 | |||
136 | /* if UNK option, do nothing further. | ||
137 | * Otherwise, make a random 16 bit integer | ||
138 | * to overwrite the least significant word | ||
139 | * of unkans. | ||
140 | */ | ||
141 | #ifdef UNK | ||
142 | /* do nothing */ | ||
143 | #else | ||
144 | ranwh(); | ||
145 | r = sx * sy + sz; | ||
146 | #endif | ||
147 | |||
148 | #ifdef DEC | ||
149 | /* To make the numbers as similar as possible | ||
150 | * in all arithmetics, the random integer has | ||
151 | * to be inserted 3 bits higher up in a DEC number. | ||
152 | * An alternative would be put it 3 bits lower down | ||
153 | * in all the other number types. | ||
154 | */ | ||
155 | s = unkans.s[2]; | ||
156 | t = s & 07; /* save these bits to put in at the bottom */ | ||
157 | s &= 0177770; | ||
158 | s |= (r >> 13) & 07; | ||
159 | unkans.s[2] = s; | ||
160 | t |= r << 3; | ||
161 | unkans.s[3] = t; | ||
162 | #endif | ||
163 | |||
164 | #ifdef IBMPC | ||
165 | unkans.s[0] = r; | ||
166 | #endif | ||
167 | |||
168 | #ifdef MIEEE | ||
169 | unkans.s[3] = r; | ||
170 | #endif | ||
171 | |||
172 | *a = unkans.d; | ||
173 | return 0; | ||
174 | } | ||