1 // 2 // Copyright (c) 2017 The Khronos Group Inc. 3 // 4 // Licensed under the Apache License, Version 2.0 (the "License"); 5 // you may not use this file except in compliance with the License. 6 // You may obtain a copy of the License at 7 // 8 // http://www.apache.org/licenses/LICENSE-2.0 9 // 10 // Unless required by applicable law or agreed to in writing, software 11 // distributed under the License is distributed on an "AS IS" BASIS, 12 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13 // See the License for the specific language governing permissions and 14 // limitations under the License. 15 // 16 #ifndef TEST_CONFORMANCE_CLCPP_MATH_FUNCS_REFERENCE_HPP 17 #define TEST_CONFORMANCE_CLCPP_MATH_FUNCS_REFERENCE_HPP 18 19 #include <type_traits> 20 #include <cmath> 21 #include <limits> 22 23 #include "../common.hpp" 24 25 namespace reference 26 { 27 // Reference functions for OpenCL comparison functions that 28 // are not already defined in STL. maxmag(const cl_float & x,const cl_float & y)29 cl_float maxmag(const cl_float& x, const cl_float& y) 30 { 31 if((std::abs)(x) > (std::abs)(y)) 32 { 33 return x; 34 } 35 else if((std::abs)(y) > (std::abs)(x)) 36 { 37 return y; 38 } 39 return (std::fmax)(x, y); 40 } 41 minmag(const cl_float & x,const cl_float & y)42 cl_float minmag(const cl_float& x, const cl_float& y) 43 { 44 if((std::abs)(x) < (std::abs)(y)) 45 { 46 return x; 47 } 48 else if((std::abs)(y) < (std::abs)(x)) 49 { 50 return y; 51 } 52 return (std::fmin)(x, y); 53 } 54 55 // Reference functions for OpenCL exp functions that 56 // are not already defined in STL. exp10(const cl_double & x)57 cl_double exp10(const cl_double& x) 58 { 59 // 10^x = exp2( x * log2(10) ) 60 auto log2_10 = (std::log2)(static_cast<long double>(10.0)); 61 cl_double x_log2_10 = static_cast<cl_double>(x * log2_10); 62 return (std::exp2)(x_log2_10); 63 } 64 65 // Reference functions for OpenCL floating point functions that 66 // are not already defined in STL. fract(cl_double x)67 cl_double2 fract(cl_double x) 68 { 69 // Copied from math_brute_force/reference_math.c 70 cl_double2 r; 71 if((std::isnan)(x)) 72 { 73 r.s[0] = std::numeric_limits<cl_double>::quiet_NaN(); 74 r.s[1] = std::numeric_limits<cl_double>::quiet_NaN(); 75 return r; 76 } 77 78 r.s[0] = (std::modf)(x, &(r.s[1])); 79 if(r.s[0] < 0.0 ) 80 { 81 r.s[0] = 1.0f + r.s[0]; 82 r.s[1] -= 1.0f; 83 if( r.s[0] == 1.0f ) 84 r.s[0] = HEX_FLT(+, 1, fffffe, -, 1); 85 } 86 return r; 87 } 88 remquo(cl_double x,cl_double y)89 cl_double2 remquo(cl_double x, cl_double y) 90 { 91 cl_double2 r; 92 // remquo return the same value that is returned by the 93 // remainder function 94 r.s[0] = (std::remainder)(x,y); 95 // calulcate quo 96 cl_double x_y = (x - r.s[0]) / y; 97 cl_uint quo = (std::abs)(x_y); 98 r.s[1] = quo & 0x0000007fU; 99 if(x_y < 0.0) 100 r.s[1] = -r.s[1]; 101 102 // fix edge cases 103 if(!(std::isnan)(x) && y == 0.0) 104 { 105 r.s[1] = 0; 106 } 107 else if((std::isnan)(x) && (std::isnan)(y)) 108 { 109 r.s[1] = 0; 110 } 111 return r; 112 } 113 114 // Reference functions for OpenCL half_math:: functions that 115 // are not already defined in STL. divide(cl_double x,cl_double y)116 cl_double divide(cl_double x, cl_double y) 117 { 118 return x / y; 119 } 120 recip(cl_double x)121 cl_double recip(cl_double x) 122 { 123 return 1.0 / x; 124 } 125 126 // Reference functions for OpenCL other functions that 127 // are not already defined in STL. mad(cl_double x,cl_double y,cl_double z)128 cl_double mad(cl_double x, cl_double y, cl_double z) 129 { 130 return (x * y) + z; 131 } 132 133 // Reference functions for OpenCL power functions that 134 // are not already defined in STL. rsqrt(const cl_double & x)135 cl_double rsqrt(const cl_double& x) 136 { 137 return cl_double(1.0) / ((std::sqrt)(x)); 138 } 139 powr(const cl_double & x,const cl_double & y)140 cl_double powr(const cl_double& x, const cl_double& y) 141 { 142 //powr(x, y) returns NaN for x < 0. 143 if( x < 0.0 ) 144 return std::numeric_limits<cl_double>::quiet_NaN(); 145 146 //powr ( x, NaN ) returns the NaN for x >= 0. 147 //powr ( NaN, y ) returns the NaN. 148 if((std::isnan)(x) || (std::isnan)(y) ) 149 return std::numeric_limits<cl_double>::quiet_NaN(); 150 151 if( x == 1.0 ) 152 { 153 //powr ( +1, +-inf ) returns NaN. 154 if((std::abs)(y) == INFINITY ) 155 return std::numeric_limits<cl_double>::quiet_NaN(); 156 157 //powr ( +1, y ) is 1 for finite y. (NaN handled above) 158 return 1.0; 159 } 160 161 if( y == 0.0 ) 162 { 163 //powr ( +inf, +-0 ) returns NaN. 164 //powr ( +-0, +-0 ) returns NaN. 165 if( x == 0.0 || x == std::numeric_limits<cl_double>::infinity()) 166 return std::numeric_limits<cl_double>::quiet_NaN(); 167 168 //powr ( x, +-0 ) is 1 for finite x > 0. (x <= 0, NaN, INF already handled above) 169 return 1.0; 170 } 171 172 if( x == 0.0 ) 173 { 174 //powr ( +-0, -inf) is +inf. 175 //powr ( +-0, y ) is +inf for finite y < 0. 176 if( y < 0.0 ) 177 return std::numeric_limits<cl_double>::infinity(); 178 179 //powr ( +-0, y ) is +0 for y > 0. (NaN, y==0 handled above) 180 return 0.0; 181 } 182 183 // x = +inf 184 if( (std::isinf)(x) ) 185 { 186 if( y < 0 ) 187 return 0; 188 return std::numeric_limits<cl_double>::infinity(); 189 } 190 191 double fabsx = (std::abs)(x); 192 double fabsy = (std::abs)(y); 193 194 //y = +-inf cases 195 if( (std::isinf)(fabsy) ) 196 { 197 if( y < 0.0 ) 198 { 199 if( fabsx < 1.0 ) 200 return std::numeric_limits<cl_double>::infinity(); 201 return 0; 202 } 203 if( fabsx < 1.0 ) 204 return 0.0; 205 return std::numeric_limits<cl_double>::infinity(); 206 } 207 return (std::pow)(x, y); 208 } 209 rootn(const cl_double & x,const cl_int n)210 cl_double rootn(const cl_double& x, const cl_int n) 211 { 212 //rootn (x, 0) returns a NaN. 213 if(n == 0) 214 return std::numeric_limits<cl_double>::quiet_NaN(); 215 216 //rootn ( x, n ) returns a NaN for x < 0 and n is even. 217 if(x < 0 && 0 == (n & 1)) 218 return std::numeric_limits<cl_double>::quiet_NaN(); 219 220 if(x == 0.0) 221 { 222 if(n > 0) 223 { 224 //rootn ( +-0, n ) is +0 for even n > 0. 225 if(0 == (n & 1)) 226 { 227 return cl_double(0.0); 228 } 229 //rootn ( +-0, n ) is +-0 for odd n > 0. 230 else 231 { 232 return x; 233 } 234 } 235 else 236 { 237 //rootn ( +-0, n ) is +inf for even n < 0. 238 if(0 == ((-n) & 1)) 239 { 240 return std::numeric_limits<cl_double>::infinity(); 241 } 242 //rootn ( +-0, n ) is +-inf for odd n < 0. 243 else 244 { 245 return (std::copysign)( 246 std::numeric_limits<cl_double>::infinity(), x 247 ); 248 } 249 } 250 } 251 252 cl_double r = (std::abs)(x); 253 r = (std::exp2)((std::log2)(r) / static_cast<cl_double>(n)); 254 return (std::copysign)(r, x); 255 } 256 257 // Reference functions for OpenCL trigonometric functions that 258 // are not already defined in STL. acospi(cl_double x)259 cl_double acospi(cl_double x) 260 { 261 return (std::acos)(x) / CL_M_PI; 262 } 263 asinpi(cl_double x)264 cl_double asinpi(cl_double x) 265 { 266 return (std::asin)(x) / CL_M_PI; 267 } 268 atanpi(cl_double x)269 cl_double atanpi(cl_double x) 270 { 271 return (std::atan)(x) / CL_M_PI; 272 } 273 cospi(cl_double x)274 cl_double cospi(cl_double x) 275 { 276 return (std::cos)(x * CL_M_PI); 277 } 278 sinpi(cl_double x)279 cl_double sinpi(cl_double x) 280 { 281 return (std::sin)(x * CL_M_PI); 282 } 283 tanpi(cl_double x)284 cl_double tanpi(cl_double x) 285 { 286 return (std::tan)(x * CL_M_PI); 287 } 288 atan2(cl_double x,cl_double y)289 cl_double atan2(cl_double x, cl_double y) 290 { 291 #if defined(WIN32) || defined(_WIN32) 292 // Fix edge cases for Windows 293 if ((std::isinf)(x) && (std::isinf)(y)) { 294 cl_double retval = (y > 0) ? CL_M_PI_4 : 3.f * CL_M_PI_4; 295 return (x > 0) ? retval : -retval; 296 } 297 #endif // defined(WIN32) || defined(_WIN32) 298 return (std::atan2)(x, y); 299 } 300 atan2pi(cl_double x,cl_double y)301 cl_double atan2pi(cl_double x, cl_double y) 302 { 303 return ::reference::atan2(x, y) / CL_M_PI; 304 } 305 sincos(cl_double x)306 cl_double2 sincos(cl_double x) 307 { 308 cl_double2 r; 309 r.s[0] = (std::sin)(x); 310 r.s[1] = (std::cos)(x); 311 return r; 312 } 313 } 314 315 #endif // TEST_CONFORMANCE_CLCPP_MATH_FUNCS_REFERENCE_HPP 316