summaryrefslogtreecommitdiff
path: root/src/regress/lib/libc/cephes/drand.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/regress/lib/libc/cephes/drand.c')
-rw-r--r--src/regress/lib/libc/cephes/drand.c174
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
69static int sx = 1;
70static int sy = 10000;
71static int sz = 3000;
72
73static union {
74 double d;
75 unsigned short s[4];
76} unkans;
77
78/* This function implements the three
79 * congruential generators.
80 */
81
82static int ranwh()
83{
84int r, s;
85
86/* sx = sx * 171 mod 30269 */
87r = sx/177;
88s = sx - 177 * r;
89sx = 171 * s - 2 * r;
90if( sx < 0 )
91 sx += 30269;
92
93
94/* sy = sy * 172 mod 30307 */
95r = sy/176;
96s = sy - 176 * r;
97sy = 172 * s - 35 * r;
98if( sy < 0 )
99 sy += 30307;
100
101/* sz = 170 * sz mod 30323 */
102r = sz/178;
103s = sz - 178 * r;
104sz = 170 * s - 63 * r;
105if( sz < 0 )
106 sz += 30323;
107/* The results are in static sx, sy, sz. */
108return 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
119int drand( a )
120double *a;
121{
122unsigned short r;
123#ifdef DEC
124unsigned short s, t;
125#endif
126
127/* This algorithm of Wichmann and Hill computes a floating point
128 * result:
129 */
130ranwh();
131unkans.d = sx/30269.0 + sy/30307.0 + sz/30323.0;
132r = unkans.d;
133unkans.d -= r;
134unkans.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
144ranwh();
145r = 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 */
155s = unkans.s[2];
156t = s & 07; /* save these bits to put in at the bottom */
157s &= 0177770;
158s |= (r >> 13) & 07;
159unkans.s[2] = s;
160t |= r << 3;
161unkans.s[3] = t;
162#endif
163
164#ifdef IBMPC
165unkans.s[0] = r;
166#endif
167
168#ifdef MIEEE
169unkans.s[3] = r;
170#endif
171
172*a = unkans.d;
173return 0;
174}