1*21508fa7SKévin Petit#include <clc/clc.h> 2*21508fa7SKévin Petit 3*21508fa7SKévin Petit#pragma OPENCL EXTENSION cl_khr_byte_addressable_store : enable 4*21508fa7SKévin Petit 5*21508fa7SKévin Petit#define ROUND_VEC1(out, in, ROUNDF) out = ROUNDF(in); 6*21508fa7SKévin Petit#define ROUND_VEC2(out, in, ROUNDF) \ 7*21508fa7SKévin Petit ROUND_VEC1(out.lo, in.lo, ROUNDF); \ 8*21508fa7SKévin Petit ROUND_VEC1(out.hi, in.hi, ROUNDF); 9*21508fa7SKévin Petit#define ROUND_VEC3(out, in, ROUNDF) \ 10*21508fa7SKévin Petit ROUND_VEC1(out.s0, in.s0, ROUNDF); \ 11*21508fa7SKévin Petit ROUND_VEC1(out.s1, in.s1, ROUNDF); \ 12*21508fa7SKévin Petit ROUND_VEC1(out.s2, in.s2, ROUNDF); 13*21508fa7SKévin Petit#define ROUND_VEC4(out, in, ROUNDF) \ 14*21508fa7SKévin Petit ROUND_VEC2(out.lo, in.lo, ROUNDF); \ 15*21508fa7SKévin Petit ROUND_VEC2(out.hi, in.hi, ROUNDF); 16*21508fa7SKévin Petit#define ROUND_VEC8(out, in, ROUNDF) \ 17*21508fa7SKévin Petit ROUND_VEC4(out.lo, in.lo, ROUNDF); \ 18*21508fa7SKévin Petit ROUND_VEC4(out.hi, in.hi, ROUNDF); 19*21508fa7SKévin Petit#define ROUND_VEC16(out, in, ROUNDF) \ 20*21508fa7SKévin Petit ROUND_VEC8(out.lo, in.lo, ROUNDF); \ 21*21508fa7SKévin Petit ROUND_VEC8(out.hi, in.hi, ROUNDF); 22*21508fa7SKévin Petit 23*21508fa7SKévin Petit#define __FUNC(SUFFIX, VEC_SIZE, TYPE, AS, ROUNDF) \ 24*21508fa7SKévin Petit void _CLC_OVERLOAD vstore_half_##VEC_SIZE(TYPE, size_t, AS half *); \ 25*21508fa7SKévin Petit _CLC_OVERLOAD _CLC_DEF void vstore_half##SUFFIX(TYPE vec, size_t offset, \ 26*21508fa7SKévin Petit AS half *mem) { \ 27*21508fa7SKévin Petit TYPE rounded_vec; \ 28*21508fa7SKévin Petit ROUND_VEC##VEC_SIZE(rounded_vec, vec, ROUNDF); \ 29*21508fa7SKévin Petit vstore_half_##VEC_SIZE(rounded_vec, offset, mem); \ 30*21508fa7SKévin Petit } \ 31*21508fa7SKévin Petit void _CLC_OVERLOAD vstorea_half_##VEC_SIZE(TYPE, size_t, AS half *); \ 32*21508fa7SKévin Petit _CLC_OVERLOAD _CLC_DEF void vstorea_half##SUFFIX(TYPE vec, size_t offset, \ 33*21508fa7SKévin Petit AS half *mem) { \ 34*21508fa7SKévin Petit TYPE rounded_vec; \ 35*21508fa7SKévin Petit ROUND_VEC##VEC_SIZE(rounded_vec, vec, ROUNDF); \ 36*21508fa7SKévin Petit vstorea_half_##VEC_SIZE(rounded_vec, offset, mem); \ 37*21508fa7SKévin Petit } 38*21508fa7SKévin Petit 39*21508fa7SKévin Petit_CLC_DEF _CLC_OVERLOAD float __clc_rtz(float x) { 40*21508fa7SKévin Petit /* Handle nan corner case */ 41*21508fa7SKévin Petit if (isnan(x)) 42*21508fa7SKévin Petit return x; 43*21508fa7SKévin Petit /* RTZ does not produce Inf for large numbers */ 44*21508fa7SKévin Petit if (fabs(x) > 65504.0f && !isinf(x)) 45*21508fa7SKévin Petit return copysign(65504.0f, x); 46*21508fa7SKévin Petit 47*21508fa7SKévin Petit const int exp = (as_uint(x) >> 23 & 0xff) - 127; 48*21508fa7SKévin Petit /* Manage range rounded to +- zero explicitely */ 49*21508fa7SKévin Petit if (exp < -24) 50*21508fa7SKévin Petit return copysign(0.0f, x); 51*21508fa7SKévin Petit 52*21508fa7SKévin Petit /* Remove lower 13 bits to make sure the number is rounded down */ 53*21508fa7SKévin Petit int mask = 0xffffe000; 54*21508fa7SKévin Petit /* Denormals cannot be flushed, and they use different bit for rounding */ 55*21508fa7SKévin Petit if (exp < -14) 56*21508fa7SKévin Petit mask <<= min(-(exp + 14), 10); 57*21508fa7SKévin Petit 58*21508fa7SKévin Petit return as_float(as_uint(x) & mask); 59*21508fa7SKévin Petit} 60*21508fa7SKévin Petit 61*21508fa7SKévin Petit_CLC_DEF _CLC_OVERLOAD float __clc_rti(float x) { 62*21508fa7SKévin Petit /* Handle nan corner case */ 63*21508fa7SKévin Petit if (isnan(x)) 64*21508fa7SKévin Petit return x; 65*21508fa7SKévin Petit 66*21508fa7SKévin Petit const float inf = copysign(INFINITY, x); 67*21508fa7SKévin Petit uint ux = as_uint(x); 68*21508fa7SKévin Petit 69*21508fa7SKévin Petit /* Manage +- infinity explicitely */ 70*21508fa7SKévin Petit if (as_float(ux & 0x7fffffff) > 0x1.ffcp+15f) { 71*21508fa7SKévin Petit return inf; 72*21508fa7SKévin Petit } 73*21508fa7SKévin Petit /* Manage +- zero explicitely */ 74*21508fa7SKévin Petit if ((ux & 0x7fffffff) == 0) { 75*21508fa7SKévin Petit return copysign(0.0f, x); 76*21508fa7SKévin Petit } 77*21508fa7SKévin Petit 78*21508fa7SKévin Petit const int exp = (as_uint(x) >> 23 & 0xff) - 127; 79*21508fa7SKévin Petit /* Manage range rounded to smallest half denormal explicitely */ 80*21508fa7SKévin Petit if (exp < -24) { 81*21508fa7SKévin Petit return copysign(0x1.0p-24f, x); 82*21508fa7SKévin Petit } 83*21508fa7SKévin Petit 84*21508fa7SKévin Petit /* Set lower 13 bits */ 85*21508fa7SKévin Petit int mask = (1 << 13) - 1; 86*21508fa7SKévin Petit /* Denormals cannot be flushed, and they use different bit for rounding */ 87*21508fa7SKévin Petit if (exp < -14) { 88*21508fa7SKévin Petit mask = (1 << (13 + min(-(exp + 14), 10))) - 1; 89*21508fa7SKévin Petit } 90*21508fa7SKévin Petit 91*21508fa7SKévin Petit const float next = nextafter(as_float(ux | mask), inf); 92*21508fa7SKévin Petit return ((ux & mask) == 0) ? as_float(ux) : next; 93*21508fa7SKévin Petit} 94*21508fa7SKévin Petit_CLC_DEF _CLC_OVERLOAD float __clc_rtn(float x) { 95*21508fa7SKévin Petit return ((as_uint(x) & 0x80000000) == 0) ? __clc_rtz(x) : __clc_rti(x); 96*21508fa7SKévin Petit} 97*21508fa7SKévin Petit_CLC_DEF _CLC_OVERLOAD float __clc_rtp(float x) { 98*21508fa7SKévin Petit return ((as_uint(x) & 0x80000000) == 0) ? __clc_rti(x) : __clc_rtz(x); 99*21508fa7SKévin Petit} 100*21508fa7SKévin Petit_CLC_DEF _CLC_OVERLOAD float __clc_rte(float x) { 101*21508fa7SKévin Petit /* Mantisa + implicit bit */ 102*21508fa7SKévin Petit const uint mantissa = (as_uint(x) & 0x7fffff) | (1u << 23); 103*21508fa7SKévin Petit const int exp = (as_uint(x) >> 23 & 0xff) - 127; 104*21508fa7SKévin Petit int shift = 13; 105*21508fa7SKévin Petit if (exp < -14) { 106*21508fa7SKévin Petit /* The default assumes lower 13 bits are rounded, 107*21508fa7SKévin Petit * but it might be more for denormals. 108*21508fa7SKévin Petit * Shifting beyond last == 0b, and qr == 00b is not necessary */ 109*21508fa7SKévin Petit shift += min(-(exp + 14), 15); 110*21508fa7SKévin Petit } 111*21508fa7SKévin Petit int mask = (1 << shift) - 1; 112*21508fa7SKévin Petit const uint grs = mantissa & mask; 113*21508fa7SKévin Petit const uint last = mantissa & (1 << shift); 114*21508fa7SKévin Petit /* IEEE round up rule is: grs > 101b or grs == 100b and last == 1. 115*21508fa7SKévin Petit * exp > 15 should round to inf. */ 116*21508fa7SKévin Petit bool roundup = (grs > (1 << (shift - 1))) || 117*21508fa7SKévin Petit (grs == (1 << (shift - 1)) && last != 0) || (exp > 15); 118*21508fa7SKévin Petit return roundup ? __clc_rti(x) : __clc_rtz(x); 119*21508fa7SKévin Petit} 120*21508fa7SKévin Petit 121*21508fa7SKévin Petit#define __XFUNC(SUFFIX, VEC_SIZE, TYPE, AS) \ 122*21508fa7SKévin Petit __FUNC(SUFFIX, VEC_SIZE, TYPE, AS, __clc_rte) \ 123*21508fa7SKévin Petit __FUNC(SUFFIX##_rtz, VEC_SIZE, TYPE, AS, __clc_rtz) \ 124*21508fa7SKévin Petit __FUNC(SUFFIX##_rtn, VEC_SIZE, TYPE, AS, __clc_rtn) \ 125*21508fa7SKévin Petit __FUNC(SUFFIX##_rtp, VEC_SIZE, TYPE, AS, __clc_rtp) \ 126*21508fa7SKévin Petit __FUNC(SUFFIX##_rte, VEC_SIZE, TYPE, AS, __clc_rte) 127*21508fa7SKévin Petit 128*21508fa7SKévin Petit#define FUNC(SUFFIX, VEC_SIZE, TYPE, AS) __XFUNC(SUFFIX, VEC_SIZE, TYPE, AS) 129*21508fa7SKévin Petit 130*21508fa7SKévin Petit#define __CLC_BODY "vstore_half.inc" 131*21508fa7SKévin Petit#include <clc/math/gentype.inc> 132*21508fa7SKévin Petit#undef __CLC_BODY 133*21508fa7SKévin Petit#undef FUNC 134*21508fa7SKévin Petit#undef __XFUNC 135*21508fa7SKévin Petit#undef __FUNC 136