uniform.hpp (11442B)
1 /* 2 Copyright 2010-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 33 #ifndef __r123_uniform_dot_hpp 34 #define __r123_uniform_dot_hpp 35 36 /** @defgroup uniform Uniform distribution scalar conversion functions 37 38 This file provides some simple functions that can be used to convert 39 integers of various widths to floats and doubles with various 40 characteristics. It can be used to generate real-valued, uniformly 41 distributed random variables from the random integers produced by 42 the Random123 CBRNGs. 43 44 There are three templated functions: 45 46 - u01: output is as dense as possible in (0,1}, never 0.0. May 47 return 1.0 if and only if the number of output mantissa bits 48 is less than the width of the input. 49 50 - uneg11: output is as dense as possible in {-1,1}, never 0.0. May 51 return 1.0 or -1.0 if and only if the number of output mantissa bits 52 is less than the width of the input. 53 54 - u01fixedpt: output is "fixed point", equispaced, open at both ends, 55 and is never 0.0, 0.5 nor 1.0. 56 57 The behavior of u01 and uneg11 depend on the pre-processor symbol: 58 R123_UNIFORM_FLOAT_STORE. When #defined to a non-zero value, u01 59 and uneg11 declare a volatile intermediate result, with the 60 intention of forcing architectures that have "extra bits" in their 61 floating point registers to more closely conform to IEEE 62 arithmetic. When compiled this way, u01 and uneg11 will be 63 significantly slower, as they will incur a memory write and read on 64 every call. Without it, they may fail the "known answer test" 65 implemented in ut_uniform_IEEEkat.cpp even though they perform 66 perfectly reasonable int to float conversions. We have used 67 this option to get 32-bit x86 to produce the same results as 68 64-bit x86-64 code, but we do not recommend it for normal 69 use. 70 71 Three additional functions are defined when C++11 or newer is in use: 72 73 - u01all 74 - uneg11all 75 - u01fixedptall 76 77 These functions apply the corresponding conversion to every 78 element of their argument, which must be a staticly sized 79 array, e.g., an r123array or a std::array of an integer type. 80 81 This file may not be as portable, and has not been tested as 82 rigorously as other files in the library, e.g., the generators. 83 Nevertheless, we hope it is useful and we encourage developers to 84 copy it and modify it for their own use. We invite comments and 85 improvements. 86 */ 87 88 #include <Random123/features/compilerfeatures.h> 89 #include <limits> 90 #if R123_USE_CXX11_TYPE_TRAITS 91 #include <type_traits> 92 #endif 93 #if __cplusplus >= 201103L 94 #include <array> 95 #endif 96 97 namespace r123{ 98 /** 99 @{ 100 @cond HIDDEN_FROM_DOXYGEN 101 */ 102 103 #if R123_USE_CXX11_TYPE_TRAITS 104 using std::make_signed; 105 using std::make_unsigned; 106 #else 107 // Sigh... We could try to find another <type_traits>, e.g., from 108 // boost or TR1. Or we can do it ourselves in the r123 namespace. 109 // It's not clear which will cause less headache... 110 template <typename T> struct make_signed{}; 111 template <typename T> struct make_unsigned{}; 112 #define R123_MK_SIGNED_UNSIGNED(ST, UT) \ 113 template<> struct make_signed<ST>{ typedef ST type; }; \ 114 template<> struct make_signed<UT>{ typedef ST type; }; \ 115 template<> struct make_unsigned<ST>{ typedef UT type; }; \ 116 template<> struct make_unsigned<UT>{ typedef UT type; } 117 118 R123_MK_SIGNED_UNSIGNED(int8_t, uint8_t); 119 R123_MK_SIGNED_UNSIGNED(int16_t, uint16_t); 120 R123_MK_SIGNED_UNSIGNED(int32_t, uint32_t); 121 R123_MK_SIGNED_UNSIGNED(int64_t, uint64_t); 122 #if R123_USE_GNU_UINT128 123 R123_MK_SIGNED_UNSIGNED(__int128_t, __uint128_t); 124 #endif 125 #undef R123_MK_SIGNED_UNSIGNED 126 #endif 127 128 #if defined(__CUDACC__) || defined(_LIBCPP_HAS_NO_CONSTEXPR) 129 // Amazing! cuda thinks numeric_limits::max() is a __host__ function, so 130 // we can't use it in a device function. 131 // 132 // The LIBCPP_HAS_NO_CONSTEXP test catches situations where the libc++ 133 // library thinks that the compiler doesn't support constexpr, but we 134 // think it does. As a consequence, the library declares 135 // numeric_limits::max without constexpr. This workaround should only 136 // affect a narrow range of compiler/library pairings. 137 // 138 // In both cases, we find max() by computing ~(unsigned)0 right-shifted 139 // by is_signed. 140 template <typename T> 141 R123_CONSTEXPR R123_STATIC_INLINE R123_CUDA_DEVICE T maxTvalue(){ 142 typedef typename make_unsigned<T>::type uT; 143 return (~uT(0)) >> std::numeric_limits<T>::is_signed; 144 } 145 #else 146 template <typename T> 147 R123_CONSTEXPR R123_STATIC_INLINE T maxTvalue(){ 148 return std::numeric_limits<T>::max(); 149 } 150 #endif 151 /** @endcond 152 @} 153 */ 154 155 //! Return a uniform real value in (0, 1] 156 /** 157 @ingroup uniform 158 Input is a W-bit integer (signed or unsigned). It is cast to 159 a W-bit unsigned integer, multiplied by Ftype(2^-W) and added to 160 Ftype(2^(-W-1)). A good compiler should optimize it down to an 161 int-to-float conversion followed by a multiply and an add, which 162 might be fused, depending on the architecture. 163 164 If the input is a uniformly distributed integer, and if Ftype 165 arithmetic follows IEEE754 round-to-nearest rules, then the 166 result is a uniformly distributed floating point number in (0, 1]. 167 168 - The result is never exactly 0.0. 169 - The smallest value returned is 2^-(W-1). 170 - Let M be the number of mantissa bits in Ftype (typically 24 or 53). 171 - If W>M then the largest value retured is 1.0. 172 - If W<=M then the largest value returned is Ftype(1.0 - 2^(-W-1)). 173 */ 174 template <typename Ftype, typename Itype> 175 R123_CUDA_DEVICE R123_STATIC_INLINE Ftype u01(Itype in){ 176 typedef typename make_unsigned<Itype>::type Utype; 177 R123_CONSTEXPR Ftype factor = Ftype(1.)/(Ftype(maxTvalue<Utype>()) + Ftype(1.)); 178 R123_CONSTEXPR Ftype halffactor = Ftype(0.5)*factor; 179 #if R123_UNIFORM_FLOAT_STORE 180 volatile Ftype x = Utype(in)*factor; return x+halffactor; 181 #else 182 return Utype(in)*factor + halffactor; 183 #endif 184 } 185 186 //! Return a signed value in [-1,1] 187 /** 188 @ingroup uniform 189 The argument is converted to a W-bit signed integer, multiplied by Ftype(2^-(W-1)) and 190 then added to Ftype(2^-W). A good compiler should optimize 191 it down to an int-to-float conversion followed by a multiply and 192 an add, which might be fused, depending on the architecture. 193 194 If the input is a uniformly distributed integer, and if Ftype 195 arithmetic follows IEEE754 round-to-nearest rules, then the 196 output is a uniformly distributed floating point number in [-1, 1]. 197 198 - The result is never exactly 0.0. 199 - The smallest absolute value returned is 2^-W 200 - Let M be the number of mantissa bits in Ftype. 201 - If W>M then the largest value retured is 1.0 and the smallest is -1.0. 202 - If W<=M then the largest value returned is the Ftype(1.0 - 2^-W) 203 and the smallest value returned is -Ftype(1.0 - 2^-W). 204 */ 205 template <typename Ftype, typename Itype> 206 R123_CUDA_DEVICE R123_STATIC_INLINE Ftype uneg11(Itype in){ 207 typedef typename make_signed<Itype>::type Stype; 208 R123_CONSTEXPR Ftype factor = Ftype(1.)/(Ftype(maxTvalue<Stype>()) + Ftype(1.)); 209 R123_CONSTEXPR Ftype halffactor = Ftype(0.5)*factor; 210 #if R123_UNIFORM_FLOAT_STORE 211 volatile Ftype x = Stype(in)*factor; return x+halffactor; 212 #else 213 return Stype(in)*factor + halffactor; 214 #endif 215 } 216 217 //! Return a value in (0,1) chosen from a set of equally spaced fixed-point values 218 /** 219 @ingroup uniform 220 Let: 221 - W = width of Itype, e.g., 32 or 64, regardless of signedness. 222 - M = mantissa bits of Ftype, e.g., 24, 53 or 64 223 - B = min(M, W) 224 225 Then the 2^(B-1) possible output values are: 2^-B*{1, 3, 5, ..., 2^B - 1} 226 227 The smallest output is: 2^-B 228 229 The largest output is: 1 - 2^-B 230 231 The output is never exactly 0.0, nor 0.5, nor 1.0. 232 233 The 2^(B-1) possible outputs: 234 - are equally likely, 235 - are uniformly spaced by 2^-(B-1), 236 - are balanced around 0.5 237 */ 238 template <typename Ftype, typename Itype> 239 R123_CUDA_DEVICE R123_STATIC_INLINE Ftype u01fixedpt(Itype in){ 240 typedef typename make_unsigned<Itype>::type Utype; 241 R123_CONSTEXPR int excess = std::numeric_limits<Utype>::digits - std::numeric_limits<Ftype>::digits; 242 if(excess>=0){ 243 R123_CONSTEXPR int ex_nowarn = (excess>=0) ? excess : 0; 244 R123_CONSTEXPR Ftype factor = Ftype(1.)/(Ftype(1.) + Ftype((maxTvalue<Utype>()>>ex_nowarn))); 245 return (1 | (Utype(in)>>ex_nowarn)) * factor; 246 }else 247 return u01<Ftype>(in); 248 } 249 250 #if R123_USE_CXX11_STD_ARRAY 251 252 //! Apply u01 to every item in an r123array, returning a std::array 253 /** @ingroup uniform 254 * Only in C++11 and newer. 255 * The argument type may be any integer collection with a constexpr static_size member, 256 * e.g., an r123array or a std::array of an integer type. 257 */ 258 template <typename Ftype, typename CollType> 259 static inline 260 std::array<Ftype, CollType::static_size> u01all(CollType in) 261 { 262 std::array<Ftype, CollType::static_size> ret; 263 auto p = ret.begin(); 264 for(auto e : in){ 265 *p++ = u01<Ftype>(e); 266 } 267 return ret; 268 } 269 270 //! Apply uneg11 to every item in an r123array, returning a std::array 271 /** @ingroup uniform 272 * Only in C++11 and newer. 273 * The argument type may be any integer collection with a constexpr static_size member, 274 * e.g., an r123array or a std::array of an integer type. 275 */ 276 template <typename Ftype, typename CollType> 277 static inline 278 std::array<Ftype, CollType::static_size> uneg11all(CollType in) 279 { 280 std::array<Ftype, CollType::static_size> ret; 281 auto p = ret.begin(); 282 for(auto e : in){ 283 *p++ = uneg11<Ftype>(e); 284 } 285 return ret; 286 } 287 288 //! Apply u01fixedpt to every item in an r123array, returning a std::array 289 /** @ingroup uniform 290 * Only in C++11 and newer. 291 * The argument type may be any integer collection with a constexpr static_size member, 292 * e.g., an r123array or a std::array of an integer type. 293 */ 294 template <typename Ftype, typename CollType> 295 static inline 296 std::array<Ftype, CollType::static_size> u01fixedptall(CollType in) 297 { 298 std::array<Ftype, CollType::static_size> ret; 299 auto p = ret.begin(); 300 for(auto e : in){ 301 *p++ = u01fixedpt<Ftype>(e); 302 } 303 return ret; 304 } 305 #endif // __cplusplus >= 201103L 306 307 } // namespace r123 308 309 #endif 310