u01fixedpt.h (7927B)
1 /* 2 Copyright 2011, D. E. Shaw Research. 3 All rights reserved. 4 5 Redistribution and use in source and binary forms, with or without 6 modification, are permitted provided that the following conditions are 7 met: 8 9 * Redistributions of source code must retain the above copyright 10 notice, this list of conditions, and the following disclaimer. 11 12 * Redistributions in binary form must reproduce the above copyright 13 notice, this list of conditions, and the following disclaimer in the 14 documentation and/or other materials provided with the distribution. 15 16 * Neither the name of D. E. Shaw Research nor the names of its 17 contributors may be used to endorse or promote products derived from 18 this software without specific prior written permission. 19 20 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 21 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 22 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 23 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 24 OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 25 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 26 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 27 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 28 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 29 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 30 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 31 */ 32 #ifndef _random123_ufixed01_dot_h_ 33 #define _random123_ufixed01_dot_h_ 34 35 #include "features/compilerfeatures.h" 36 37 /** @defgroup u01fixedpt The u01fixedpt conversion functions 38 39 These functions convert unsigned W-bit integers to uniformly 40 spaced real values (float or double) between 0.0 and 1.0 with 41 mantissas of M bits. 42 43 PLEASE THINK CAREFULLY BEFORE USING THESE FUNCTIONS. THEY MAY 44 NOT BE WHAT YOU WANT. YOU MAY BE MUCH BETTER SERVED BY THE 45 FUNCTIONS IN ./uniform.hpp. 46 47 These functions produce a finite number *uniformly spaced* values 48 in the range from 0.0 to 1.0 with uniform probability. The price 49 of uniform spacing is that they may not utilize the entire space 50 of possible outputs. E.g., u01fixedpt_closed_open_32_24 will never 51 produce a non-zero value less than 2^-24, even though such values 52 are representable in single-precision floating point. 53 54 There are 12 functions, corresponding to the following choices: 55 56 - W = 32 or 64 57 - M = 24 (float) or 53 (double) 58 - open0 or closed0 : whether the output is open or closed at 0.0 59 - open1 or closed1 : whether the output is open or closed at 1.0 60 61 The W=64 M=24 cases are not implemented. To obtain an M=24 float 62 from a uint64_t, use a cast (possibly with right-shift and bitwise 63 and) to convert some of the bits of the uint64_t to a uint32_t and 64 then use u01fixedpt_x_y_32_float. Note that the 64-bit random integers 65 produced by the Random123 library are random in "all the bits", so 66 with a little extra effort you can obtain two floats this way -- 67 one from the high bits and one from the low bits of the 64-bit 68 value. 69 70 If the output is open at one end, then the extreme 71 value (0.0 or 1.0) will never be returned. Conversely, if the output 72 is closed at one end, then the extreme value is a possible 73 return value. 74 75 The values returned are as follows. All values are returned 76 with equal frequency, except as noted in the closed_closed case: 77 78 closed_open: Let P=min(M,W) 79 there are 2^P possible output values: 80 {0, 1, 2, ..., 2^P-1}/2^P 81 82 open_closed: Let P=min(M,W) 83 there are 2^P possible values: 84 {1, 2, ..., 2^P}/2^P 85 86 open_open: Let P=min(M, W+1) 87 there are 2^(P-1) possible values: 88 {1, 3, 5, ..., 2^P-1}/2^P 89 90 closed_closed: Let P=min(M, W-1) 91 there are 1+2^P possible values: 92 {0, 1, 2, ... 2^P}/2^P 93 The extreme values (0.0 and 1.0) are 94 returned with half the frequency of 95 all others. 96 97 On x86 hardware, especially on 32bit machines, the use of 98 internal 80bit x87-style floating point may result in 99 'bonus' precision, which may cause closed intervals to not 100 be really closed, i.e. the conversions below might not 101 convert UINT{32,64}_MAX to 1.0. This sort of issue is 102 likely to occur when storing the output of a u01fixedpt_*_32_float 103 function in a double, though one can imagine getting extra 104 precision artifacts when going from 64_53 as well. Other 105 artifacts may exist on some GPU hardware. The tests in 106 kat_u01_main.h try to expose such issues, but caveat emptor. 107 108 @cond HIDDEN_FROM_DOXYGEN 109 */ 110 111 /* Hex floats were standardized by C in 1999, but weren't standardized 112 by C++ until 2011. So, we're obliged to write out our constants in 113 decimal, even though they're most naturally expressed in binary. 114 We cross our fingers and hope that the compiler does the compile-time 115 constant arithmetic properly. 116 */ 117 #define R123_0x1p_31f (1.f/(1024.f*1024.f*1024.f*2.f)) 118 #define R123_0x1p_24f (128.f*R123_0x1p_31f) 119 #define R123_0x1p_23f (256.f*R123_0x1p_31f) 120 #define R123_0x1p_32 (1./(1024.*1024.*1024.*4.)) 121 #define R123_0x1p_63 (2.*R123_0x1p_32*R123_0x1p_32) 122 #define R123_0x1p_53 (1024.*R123_0x1p_63) 123 #define R123_0x1p_52 (2048.*R123_0x1p_63) 124 125 /** @endcond */ 126 127 #ifndef R123_USE_U01_DOUBLE 128 #define R123_USE_U01_DOUBLE 1 129 #endif 130 131 #ifdef __cplusplus 132 extern "C"{ 133 #endif 134 135 /* narrowing conversions: uint32_t to float */ 136 R123_CUDA_DEVICE R123_STATIC_INLINE float u01fixedpt_closed_closed_32_float(uint32_t i){ 137 /* N.B. we ignore the high bit, so output is not monotonic */ 138 return ((i&0x7fffffc0) + (i&0x40))*R123_0x1p_31f; /* 0x1.p-31f */ 139 } 140 141 R123_CUDA_DEVICE R123_STATIC_INLINE float u01fixedpt_closed_open_32_float(uint32_t i){ 142 return (i>>8)*R123_0x1p_24f; /* 0x1.0p-24f; */ 143 } 144 145 R123_CUDA_DEVICE R123_STATIC_INLINE float u01fixedpt_open_closed_32_float(uint32_t i){ 146 return (1+(i>>8))*R123_0x1p_24f; /* *0x1.0p-24f; */ 147 } 148 149 R123_CUDA_DEVICE R123_STATIC_INLINE float u01fixedpt_open_open_32_float(uint32_t i){ 150 return (0.5f+(i>>9))*R123_0x1p_23f; /* 0x1.p-23f; */ 151 } 152 153 #if R123_USE_U01_DOUBLE 154 /* narrowing conversions: uint64_t to double */ 155 R123_CUDA_DEVICE R123_STATIC_INLINE double u01fixedpt_closed_closed_64_double(uint64_t i){ 156 /* N.B. we ignore the high bit, so output is not monotonic */ 157 return ((i&R123_64BIT(0x7ffffffffffffe00)) + (i&0x200))*R123_0x1p_63; /* 0x1.p-63; */ 158 } 159 160 R123_CUDA_DEVICE R123_STATIC_INLINE double u01fixedpt_closed_open_64_double(uint64_t i){ 161 return (i>>11)*R123_0x1p_53; /* 0x1.0p-53; */ 162 } 163 164 R123_CUDA_DEVICE R123_STATIC_INLINE double u01fixedpt_open_closed_64_double(uint64_t i){ 165 return (1+(i>>11))*R123_0x1p_53; /* 0x1.0p-53; */ 166 } 167 168 R123_CUDA_DEVICE R123_STATIC_INLINE double u01fixedpt_open_open_64_double(uint64_t i){ 169 return (0.5+(i>>12))*R123_0x1p_52; /* 0x1.0p-52; */ 170 } 171 172 /* widening conversions: u32 to double */ 173 R123_CUDA_DEVICE R123_STATIC_INLINE double u01fixedpt_closed_closed_32_double(uint32_t i){ 174 /* j = i+(i&1) takes on 2^31+1 possible values with a 'trapezoid' distribution: 175 p_j = 1 0 2 0 2 .... 2 0 2 0 1 176 j = 0 1 2 3 4 .... 2^32 177 by converting to double *before* doing the add, we don't wrap the high bit. 178 */ 179 return (((double)(i&1)) + i)*R123_0x1p_32; /* 0x1.p-32; */ 180 } 181 182 R123_CUDA_DEVICE R123_STATIC_INLINE double u01fixedpt_closed_open_32_double(uint32_t i){ 183 return i*R123_0x1p_32; /* 0x1.p-32; */ 184 } 185 186 R123_CUDA_DEVICE R123_STATIC_INLINE double u01fixedpt_open_closed_32_double(uint32_t i){ 187 return (1.+i)*R123_0x1p_32; /* 0x1.p-32; */ 188 } 189 190 R123_CUDA_DEVICE R123_STATIC_INLINE double u01fixedpt_open_open_32_double(uint32_t i){ 191 return (0.5+i)*R123_0x1p_32; /* 0x1.p-32; */ 192 } 193 #endif /* R123_USE_U01_DOUBLE */ 194 195 #ifdef __cplusplus 196 } 197 #endif 198 199 /** @} */ 200 #endif