ut_uniform_IEEEkat.cpp (13238B)
1 /* 2 Copyright 2013, 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 /* ut_uniform_IEEEkat.cpp - a "Known Answer Test" for uniform.hpp 34 35 This code tests that the compilation environment reproduces 36 exactly the behavior of u01, uneg11 and u01fixedpt on an x86-64 37 system with strict IEEE arithmetic. It is likely to fail on 38 systems that use 80-bit internal registers (e.g., 32-bit x86), and 39 systems that are smart enough to fuse floating point multiply and 40 add into a single, rounded-only-once instruction (e.g., PowerPC, 41 Fermi, newer ARMs, Haswell, Itanium, etc.). Failures in these 42 cases are *not* necessarily problematic. */ 43 44 #include <Random123/uniform.hpp> 45 #include <stdio.h> 46 #include <stdlib.h> 47 #include <map> 48 #include <string> 49 50 using namespace r123; 51 52 std::map<std::string, long double> katmap; 53 54 // Don't inline this. It's called thousands of times in fill_katmap 55 // and blows out the optimizer in older versions of clang and open64 56 // if it's inlined. 57 void insert(const char *s, long double v){ 58 katmap[std::string(s)] = v; 59 } 60 61 void fill_katmap(){ 62 #include "ut_uniform_IEEEkatvectors.hpp" 63 #if 0 // helpful for debug, 64 for(std::map<std::string, long double>::iterator p=katmap.begin(); p!=katmap.end(); ++p){ 65 fprintf(stderr, "%s -> %La\n", p->first.c_str(), p->second); 66 } 67 #endif 68 } 69 70 template <typename T> 71 typename r123::make_unsigned<T>::type U(T x){ return x; } 72 73 template <typename T> 74 typename r123::make_signed<T>::type S(T x){ return x; } 75 76 bool checking = true; 77 int nfail = 0; 78 int nuntested = 0; 79 int notfound = 0; 80 81 #define DO1(T, expr, astr) DoOne<T>(#expr, astr, expr) 82 83 #define DO(i, astr) do{ \ 84 ChkSignInvariance(#i, i); \ 85 if(std::numeric_limits<float>::digits == 24){ \ 86 DO1(float, u01<float>(i), astr); \ 87 DO1(float, uneg11<float>(i), astr); \ 88 DO1(float, u01fixedpt<float>(i), astr); \ 89 }else{ \ 90 printf("UNTESTED: %s: float does not have a 24 bit mantissa\n", #i); \ 91 nuntested++; \ 92 } \ 93 if(std::numeric_limits<double>::digits == 53){ \ 94 DO1(double, u01<double>(i), astr); \ 95 DO1(double, uneg11<double>(i), astr); \ 96 DO1(double, u01fixedpt<double>(i), astr); \ 97 }else{ \ 98 printf("UNTESTED: %s: double does not have a 53 bit mantissa\n", #i); \ 99 nuntested++; \ 100 } \ 101 if(std::numeric_limits<long double>::digits == 64){ \ 102 DO1(long double, u01<long double>(i), astr); \ 103 DO1(long double, uneg11<long double>(i), astr); \ 104 DO1(long double, u01fixedpt<long double>(i), astr); \ 105 }else{ \ 106 printf("UNTESTED: %s: long double does not have a 64 bit mantissa\n", #i); \ 107 nuntested++; \ 108 } \ 109 } while(0) 110 111 // u01, uneg11 and u01fixedpt should all depend on the bits, but not the 112 // signedness of their argument. The templated functions S(i) and U(i) 113 // return their argument cast to an approprite signed and unsigned type. 114 // ChkSignInvariance verifies that. 115 template <typename IType> 116 void ChkSignInvariance(const std::string& s, IType i){ 117 if( u01<float>(S(i)) != u01<float>(U(i)) ){ 118 printf("INVARIANT FAILURE: u01<float>(Signed(x)) != u01<float>(Unsigned(x)) x=%s\n", s.c_str()); 119 nfail++; 120 } 121 if( uneg11<float>(S(i)) != uneg11<float>(U(i)) ){ 122 printf("INVARIANT FAILURE: uneg11<float>(Signed(x)) != uneg11<float>(Unsigned(x)) x=%s\n", s.c_str()); 123 nfail++; 124 } 125 if( u01fixedpt<float>(S(i)) != u01fixedpt<float>(U(i)) ){ 126 printf("INVARIANT FAILURE: u01<float>(Signed(x)) != u01<float>(Unsigned(x)) x=%s\n", s.c_str()); 127 nfail++; 128 } 129 130 if( u01<double>(S(i)) != u01<double>(U(i)) ){ 131 printf("INVARIANT FAILURE: u01<double>(Signed(x)) != u01<double>(Unsigned(x)) x=%s\n", s.c_str()); 132 nfail++; 133 } 134 if( uneg11<double>(S(i)) != uneg11<double>(U(i)) ){ 135 printf("INVARIANT FAILURE: uneg11<double>(Signed(x)) != uneg11<double>(Unsigned(x)) x=%s\n", s.c_str()); 136 nfail++; 137 } 138 if( u01fixedpt<double>(S(i)) != u01fixedpt<double>(U(i)) ){ 139 printf("INVARIANT FAILURE: u01<double>(Signed(x)) != u01<double>(Unsigned(x)) x=%s\n", s.c_str()); 140 nfail++; 141 } 142 143 if( u01<long double>(S(i)) != u01<long double>(U(i)) ){ 144 printf("INVARIANT FAILURE: u01<long double>(Signed(x)) != u01<long double>(Unsigned(x)) x=%s\n", s.c_str()); 145 nfail++; 146 } 147 if( uneg11<long double>(S(i)) != uneg11<long double>(U(i)) ){ 148 printf("INVARIANT FAILURE: uneg11<long double>(Signed(x)) != uneg11<long double>(Unsigned(x)) x=%s\n", s.c_str()); 149 nfail++; 150 } 151 if( u01fixedpt<long double>(S(i)) != u01fixedpt<long double>(U(i)) ){ 152 printf("INVARIANT FAILURE: u01<long double>(Signed(x)) != u01<long double>(Unsigned(x)) x=%s\n", s.c_str()); 153 nfail++; 154 } 155 } 156 157 template <typename T> 158 void DoOne(const std::string s, const char* astr, volatile T x){ 159 std::string ss = s + " a=" + astr; 160 volatile long double ldx = x; 161 if(checking){ 162 if( katmap.find(ss) == katmap.end() ){ 163 printf("NOT FOUND: katmap[%s]\n", ss.c_str()); 164 notfound++; 165 }else{ 166 if(ldx!=katmap[ss]){ 167 printf("MISMATCH: %s: computed=%.21Lg reference=%.21Lg\n", ss.c_str(), ldx, katmap[ss]); 168 nfail++; 169 } 170 } 171 }else{ 172 printf("insert(\"%s\", %#.21LgL);\n", ss.c_str(), ldx); 173 } 174 } 175 176 void DO3264(int a){ 177 const uint32_t maxu32 = std::numeric_limits<uint32_t>::max(); 178 const uint64_t maxu64 = std::numeric_limits<uint64_t>::max(); 179 const uint32_t minu32 = std::numeric_limits<uint32_t>::min(); 180 const uint64_t minu64 = std::numeric_limits<uint64_t>::min(); 181 182 const int32_t maxi32 = std::numeric_limits<int32_t>::max(); 183 const int64_t maxi64 = std::numeric_limits<int64_t>::max(); 184 const int32_t mini32 = std::numeric_limits<int32_t>::min(); 185 const int64_t mini64 = std::numeric_limits<int64_t>::min(); 186 187 char astr[32]; 188 sprintf(astr, "%d", a); 189 190 DO( minu32 + uint32_t(a), astr ); 191 DO( minu64 + uint64_t(a), astr ); 192 DO( mini32 + int32_t(a), astr ); 193 DO( mini64 + int64_t(a), astr ); 194 195 DO( maxu32 - uint32_t(a), astr ); 196 DO( maxu64 - uint64_t(a), astr ); 197 DO( maxi32 - int32_t(a), astr ); 198 DO( maxi64 - int64_t(a), astr ); 199 } 200 201 int main(int argc, char **argv){ 202 if(argc>1){ 203 checking = false; 204 printf("/* This file was created by '%s %s' on a reference\n" 205 " platform, and is #included in the recompilation of %s\n" 206 " on a target platform. When %s is run with no arguments\n" 207 " on the target platform, it asserts that the values computed\n" 208 " on the target platform match the reference values recorded here.\n" 209 " These reference values were computed on an x86_64 using 32-bit,\n" 210 " 64-bit and 80-bit IEEE arithmetic for float, double and long double\n" 211 " respectively. Other platforms with different representations of\n" 212 " floating point values or different conventions for how intermediates\n" 213 " are stored and rounded will almost certainly fail these tests\n" 214 " even though their results might be perfectly valid.\n" 215 "*/\n", argv[0], argv[1], argv[0], argv[0]); 216 } 217 fill_katmap(); 218 219 220 DO3264(0); 221 DO3264(1); 222 DO3264(2); 223 DO3264(3); 224 DO3264(4); 225 DO3264(5); 226 227 DO3264(63); 228 DO3264(64); 229 DO3264(65); 230 231 DO3264(127); 232 DO3264(128); 233 DO3264(129); 234 235 DO3264(191); 236 DO3264(192); 237 DO3264(193); 238 239 DO3264(255); 240 DO3264(256); 241 DO3264(257); 242 243 DO3264(319); 244 DO3264(320); 245 DO3264(321); 246 247 DO3264(382); 248 DO3264(383); 249 DO3264(384); 250 251 DO3264(639); 252 DO3264(640); 253 DO3264(641); 254 255 DO3264(1023); 256 DO3264(1024); 257 DO3264(1025); 258 259 DO3264(3070); 260 DO3264(3071); 261 DO3264(3072); 262 263 DO3264(5119); 264 DO3264(5120); 265 DO3264(5121); 266 267 if(notfound){ 268 printf("// %s: WARNING: %d tests were not checked because reference values were not compiled in\n", 269 argv[0], notfound); 270 } 271 if(nuntested){ 272 printf("// %s: WARNING: %d tests were not performed because the floating point rep does not match the IEEE format used to compute reference values\n", 273 argv[0], nuntested); 274 } 275 if(nfail){ 276 printf("// %s: FAILED %d Known-Answer-Tests failed\n", argv[0], nfail); 277 printf("Such failures may be due to non-IEEE arithmetic on your platform. In some \n" 278 "cases, you may be able to recover IEEE arithmetic by pre-defining the\n" 279 "pp-symbol R123_UNIFORM_FLOAT_STORE to a non-zero value, e.g., adding\n" 280 "-DR123_UNIFORM_FLOAT_STORE=1 to the compile command line. On some\n" 281 "systems (notably, 32-bit x86 architectures) this will prevent use of\n" 282 "extra-wide internal floating point registers and will recover IEEE\n" 283 "arithmetic. Unfortunately, this will make u01 and uneg11 significantly\n" 284 "slower, so you may not wish to define it in production code. As far\n" 285 "as we know, the floating point values returned with the symbol unset\n" 286 "are perfectly reasonable. They simply don't perfectly match the\n" 287 "values computed on our reference x86-64 platform with IEEE arithmetic\n"); 288 }else{ 289 printf("// %s: SUCCESS\n", argv[0]); 290 } 291 292 return !!nfail; 293 }