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 } 175