361 lines · plain
1//===----------------------------------------------------------------------===//2//3// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.4// See https://llvm.org/LICENSE.txt for license information.5// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception6//7//===----------------------------------------------------------------------===//8 9_CLC_DEF _CLC_OVERLOAD __CLC_FLOATN __clc_sinf_piby4(__CLC_FLOATN x,10 __CLC_FLOATN y) {11 // Taylor series for sin(x) is x - x^3/3! + x^5/5! - x^7/7! ...12 // = x * (1 - x^2/3! + x^4/5! - x^6/7! ...13 // = x * f(w)14 // where w = x*x and f(w) = (1 - w/3! + w^2/5! - w^3/7! ...15 // We use a minimax approximation of (f(w) - 1) / w16 // because this produces an expansion in even powers of x.17 18 const __CLC_FLOATN c1 = -0.1666666666e0f;19 const __CLC_FLOATN c2 = 0.8333331876e-2f;20 const __CLC_FLOATN c3 = -0.198400874e-3f;21 const __CLC_FLOATN c4 = 0.272500015e-5f;22 const __CLC_FLOATN c5 = -2.5050759689e-08f; // 0xb2d72f3423 const __CLC_FLOATN c6 = 1.5896910177e-10f; // 0x2f2ec9d324 25 __CLC_FLOATN z = x * x;26 __CLC_FLOATN v = z * x;27 __CLC_FLOATN r = __clc_mad(28 z, __clc_mad(z, __clc_mad(z, __clc_mad(z, c6, c5), c4), c3), c2);29 __CLC_FLOATN ret =30 x - __clc_mad(v, -c1, __clc_mad(z, __clc_mad(y, 0.5f, -v * r), -y));31 32 return ret;33}34 35_CLC_DEF _CLC_OVERLOAD __CLC_FLOATN __clc_cosf_piby4(__CLC_FLOATN x,36 __CLC_FLOATN y) {37 // Taylor series for cos(x) is 1 - x^2/2! + x^4/4! - x^6/6! ...38 // = f(w)39 // where w = x*x and f(w) = (1 - w/2! + w^2/4! - w^3/6! ...40 // We use a minimax approximation of (f(w) - 1 + w/2) / (w*w)41 // because this produces an expansion in even powers of x.42 43 const __CLC_FLOATN c1 = 0.416666666e-1f;44 const __CLC_FLOATN c2 = -0.138888876e-2f;45 const __CLC_FLOATN c3 = 0.248006008e-4f;46 const __CLC_FLOATN c4 = -0.2730101334e-6f;47 const __CLC_FLOATN c5 = 2.0875723372e-09f; // 0x310f74f648 const __CLC_FLOATN c6 = -1.1359647598e-11f; // 0xad47d74e49 50 __CLC_FLOATN z = x * x;51 __CLC_FLOATN r =52 z *53 __clc_mad(54 z,55 __clc_mad(z, __clc_mad(z, __clc_mad(z, __clc_mad(z, c6, c5), c4), c3),56 c2),57 c1);58 59 // if |x| < 0.360 __CLC_FLOATN qx = 0.0f;61 62 __CLC_INTN ix = __CLC_AS_INTN(x) & EXSIGNBIT_SP32;63 64 // 0.78125 > |x| >= 0.365 __CLC_FLOATN xby4 = __CLC_AS_FLOATN(ix - 0x01000000);66 qx = ((ix >= 0x3e99999a) & (ix <= 0x3f480000)) ? xby4 : qx;67 68 // x > 0.7812569 qx = ix > 0x3f480000 ? 0.28125f : qx;70 71 __CLC_FLOATN hz = __clc_mad(z, 0.5f, -qx);72 __CLC_FLOATN a = 1.0f - qx;73 __CLC_FLOATN ret = a - (hz - __clc_mad(z, r, -x * y));74 return ret;75}76 77// Evaluate single precisions sin and cos of value in interval [-pi/4, pi/4]78_CLC_DEF _CLC_OVERLOAD void __clc_sincos_piby4(__CLC_FLOATN x,79 private __CLC_FLOATN *sinval,80 private __CLC_FLOATN *cosval) {81 // Taylor series for sin(x) is x - x^3/3! + x^5/5! - x^7/7! ...82 // = x * (1 - x^2/3! + x^4/5! - x^6/7! ...83 // = x * f(w)84 // where w = x*x and f(w) = (1 - w/3! + w^2/5! - w^3/7! ...85 // We use a minimax approximation of (f(w) - 1) / w86 // because this produces an expansion in even powers of x.87 88 // Taylor series for cos(x) is 1 - x^2/2! + x^4/4! - x^6/6! ...89 // = f(w)90 // where w = x*x and f(w) = (1 - w/2! + w^2/4! - w^3/6! ...91 // We use a minimax approximation of (f(w) - 1 + w/2) / (w*w)92 // because this produces an expansion in even powers of x.93 94 const __CLC_FLOATN sc1 = -0.166666666638608441788607926e0F;95 const __CLC_FLOATN sc2 = 0.833333187633086262120839299e-2F;96 const __CLC_FLOATN sc3 = -0.198400874359527693921333720e-3F;97 const __CLC_FLOATN sc4 = 0.272500015145584081596826911e-5F;98 99 const __CLC_FLOATN cc1 = 0.41666666664325175238031e-1F;100 const __CLC_FLOATN cc2 = -0.13888887673175665567647e-2F;101 const __CLC_FLOATN cc3 = 0.24800600878112441958053e-4F;102 const __CLC_FLOATN cc4 = -0.27301013343179832472841e-6F;103 104 __CLC_FLOATN x2 = x * x;105 106 *sinval = __clc_mad(107 x * x2, __clc_mad(x2, __clc_mad(x2, __clc_mad(x2, sc4, sc3), sc2), sc1),108 x);109 *cosval = __clc_mad(110 x2 * x2, __clc_mad(x2, __clc_mad(x2, __clc_mad(x2, cc4, cc3), cc2), cc1),111 __clc_mad(x2, -0.5f, 1.0f));112}113 114_CLC_DEF _CLC_OVERLOAD __CLC_FLOATN __clc_tanf_piby4(__CLC_FLOATN x,115 __CLC_INTN regn) {116 // Core Remez [1,2] approximation to tan(x) on the interval [0,pi/4].117 __CLC_FLOATN r = x * x;118 119 __CLC_FLOATN a =120 __clc_mad(r, -0.0172032480471481694693109f, 0.385296071263995406715129f);121 122 __CLC_FLOATN b = __clc_mad(123 r,124 __clc_mad(r, 0.01844239256901656082986661f, -0.51396505478854532132342f),125 1.15588821434688393452299f);126 127 __CLC_FLOATN t = __clc_mad(x * r, __clc_native_divide(a, b), x);128 __CLC_FLOATN tr = -MATH_RECIP(t);129 130 return (regn & 1) != 0 ? tr : t;131}132 133_CLC_DEF _CLC_OVERLOAD void __clc_fullMulS(private __CLC_FLOATN *hi,134 private __CLC_FLOATN *lo,135 __CLC_FLOATN a, __CLC_FLOATN b,136 __CLC_FLOATN bh, __CLC_FLOATN bt) {137 if (__CLC_HAVE_HW_FMA32()) {138 __CLC_FLOATN ph = a * b;139 *hi = ph;140 *lo = __clc_fma(a, b, -ph);141 } else {142 __CLC_FLOATN ah = __CLC_AS_FLOATN(__CLC_AS_UINTN(a) & 0xfffff000U);143 __CLC_FLOATN at = a - ah;144 __CLC_FLOATN ph = a * b;145 __CLC_FLOATN pt = __clc_mad(146 at, bt, __clc_mad(at, bh, __clc_mad(ah, bt, __clc_mad(ah, bh, -ph))));147 *hi = ph;148 *lo = pt;149 }150}151 152_CLC_DEF _CLC_OVERLOAD __CLC_FLOATN __clc_removePi2S(private __CLC_FLOATN *hi,153 private __CLC_FLOATN *lo,154 __CLC_FLOATN x) {155 // 72 bits of pi/2156 const __CLC_FLOATN fpiby2_1 = (__CLC_FLOATN)0xC90FDA / 0x1.0p+23f;157 const __CLC_FLOATN fpiby2_1_h = (__CLC_FLOATN)0xC90 / 0x1.0p+11f;158 const __CLC_FLOATN fpiby2_1_t = (__CLC_FLOATN)0xFDA / 0x1.0p+23f;159 160 const __CLC_FLOATN fpiby2_2 = (__CLC_FLOATN)0xA22168 / 0x1.0p+47f;161 const __CLC_FLOATN fpiby2_2_h = (__CLC_FLOATN)0xA22 / 0x1.0p+35f;162 const __CLC_FLOATN fpiby2_2_t = (__CLC_FLOATN)0x168 / 0x1.0p+47f;163 164 const __CLC_FLOATN fpiby2_3 = (__CLC_FLOATN)0xC234C4 / 0x1.0p+71f;165 const __CLC_FLOATN fpiby2_3_h = (__CLC_FLOATN)0xC23 / 0x1.0p+59f;166 const __CLC_FLOATN fpiby2_3_t = (__CLC_FLOATN)0x4C4 / 0x1.0p+71f;167 168 const __CLC_FLOATN twobypi = 0x1.45f306p-1f;169 170 __CLC_FLOATN fnpi2 = __clc_trunc(__clc_mad(x, twobypi, 0.5f));171 172 // subtract n * pi/2 from x173 __CLC_FLOATN rhead, rtail;174 __clc_fullMulS(&rhead, &rtail, fnpi2, fpiby2_1, fpiby2_1_h, fpiby2_1_t);175 __CLC_FLOATN v = x - rhead;176 __CLC_FLOATN rem = v + (((x - v) - rhead) - rtail);177 178 __CLC_FLOATN rhead2, rtail2;179 __clc_fullMulS(&rhead2, &rtail2, fnpi2, fpiby2_2, fpiby2_2_h, fpiby2_2_t);180 v = rem - rhead2;181 rem = v + (((rem - v) - rhead2) - rtail2);182 183 __CLC_FLOATN rhead3, rtail3;184 __clc_fullMulS(&rhead3, &rtail3, fnpi2, fpiby2_3, fpiby2_3_h, fpiby2_3_t);185 v = rem - rhead3;186 187 *hi = v + ((rem - v) - rhead3);188 *lo = -rtail3;189 return fnpi2;190}191 192_CLC_DEF _CLC_OVERLOAD __CLC_INTN __clc_argReductionSmallS(193 private __CLC_FLOATN *r, private __CLC_FLOATN *rr, __CLC_FLOATN x) {194 __CLC_FLOATN fnpi2 = __clc_removePi2S(r, rr, x);195 return __CLC_CONVERT_INTN(fnpi2) & 0x3;196}197 198_CLC_DEF _CLC_OVERLOAD __CLC_INTN __clc_argReductionLargeS(199 private __CLC_FLOATN *r, private __CLC_FLOATN *rr, __CLC_FLOATN x) {200 __CLC_INTN xe = __CLC_AS_INTN((__CLC_AS_UINTN(x) >> 23) - 127);201 __CLC_UINTN xm = 0x00800000U | (__CLC_AS_UINTN(x) & 0x7fffffU);202 203 // 224 bits of 2/PI: . A2F9836E 4E441529 FC2757D1 F534DDC0 DB629599 3C439041204 // FE5163AB205 const __CLC_UINTN b6 = 0xA2F9836EU;206 const __CLC_UINTN b5 = 0x4E441529U;207 const __CLC_UINTN b4 = 0xFC2757D1U;208 const __CLC_UINTN b3 = 0xF534DDC0U;209 const __CLC_UINTN b2 = 0xDB629599U;210 const __CLC_UINTN b1 = 0x3C439041U;211 const __CLC_UINTN b0 = 0xFE5163ABU;212 213 __CLC_UINTN p0, p1, p2, p3, p4, p5, p6, p7, c0, c1;214 215 __CLC_FULL_MUL(xm, b0, c0, p0);216 __CLC_FULL_MAD(xm, b1, c0, c1, p1);217 __CLC_FULL_MAD(xm, b2, c1, c0, p2);218 __CLC_FULL_MAD(xm, b3, c0, c1, p3);219 __CLC_FULL_MAD(xm, b4, c1, c0, p4);220 __CLC_FULL_MAD(xm, b5, c0, c1, p5);221 __CLC_FULL_MAD(xm, b6, c1, p7, p6);222 223 __CLC_UINTN fbits = (__CLC_UINTN)224 + (__CLC_UINTN)23 - __CLC_AS_UINTN(xe);224 225 // shift amount to get 2 lsb of integer part at top 2 bits226 // min: 25 (xe=18) max: 134 (xe=127)227 __CLC_UINTN shift = 256U - 2 - fbits;228 229 // Shift by up to 134/32 = 4 words230 __CLC_INTN c = shift > 31;231 p7 = c ? p6 : p7;232 p6 = c ? p5 : p6;233 p5 = c ? p4 : p5;234 p4 = c ? p3 : p4;235 p3 = c ? p2 : p3;236 p2 = c ? p1 : p2;237 p1 = c ? p0 : p1;238 shift -= (c ? 32U : 0U);239 240 c = shift > 31;241 p7 = c ? p6 : p7;242 p6 = c ? p5 : p6;243 p5 = c ? p4 : p5;244 p4 = c ? p3 : p4;245 p3 = c ? p2 : p3;246 p2 = c ? p1 : p2;247 shift -= (c ? 32U : 0U);248 249 c = shift > 31;250 p7 = c ? p6 : p7;251 p6 = c ? p5 : p6;252 p5 = c ? p4 : p5;253 p4 = c ? p3 : p4;254 p3 = c ? p2 : p3;255 shift -= (c ? 32U : 0U);256 257 c = shift > 31;258 p7 = c ? p6 : p7;259 p6 = c ? p5 : p6;260 p5 = c ? p4 : p5;261 p4 = c ? p3 : p4;262 shift -= (c ? 32U : 0U);263 264 // bitalign cannot handle a shift of 32265 c = shift > 0;266 shift = 32 - shift;267 __CLC_UINTN t7 = bitalign(p7, p6, shift);268 __CLC_UINTN t6 = bitalign(p6, p5, shift);269 __CLC_UINTN t5 = bitalign(p5, p4, shift);270 p7 = c ? t7 : p7;271 p6 = c ? t6 : p6;272 p5 = c ? t5 : p5;273 274 // Get 2 lsb of int part and msb of fraction275 __CLC_INTN i = __CLC_AS_INTN(p7 >> 29U);276 277 // Scoot up 2 more bits so only fraction remains278 p7 = bitalign(p7, p6, 30);279 p6 = bitalign(p6, p5, 30);280 p5 = bitalign(p5, p4, 30);281 282 // Subtract 1 if msb of fraction is 1, i.e. fraction >= 0.5283 __CLC_UINTN flip = (i & 1) != 0 ? 0xFFFFFFFFU : 0U;284 __CLC_UINTN sign = (i & 1) != 0 ? 0x80000000U : 0U;285 p7 = p7 ^ flip;286 p6 = p6 ^ flip;287 p5 = p5 ^ flip;288 289 // Find exponent and shift away leading zeroes and hidden bit290 xe = __CLC_AS_INTN(__clc_clz(p7)) + 1;291 shift = 32 - __CLC_AS_UINTN(xe);292 p7 = bitalign(p7, p6, shift);293 p6 = bitalign(p6, p5, shift);294 295 // Most significant part of fraction296 __CLC_FLOATN q1 =297 __CLC_AS_FLOATN(sign | ((127U - __CLC_AS_UINTN(xe)) << 23U) | p7 >> 9);298 299 // Shift out bits we captured on q1300 p7 = bitalign(p7, p6, 32 - 23);301 302 // Get 24 more bits of fraction in another float, there are not long strings303 // of zeroes here304 __CLC_INTN xxe = __CLC_AS_INTN(__clc_clz(p7)) + 1;305 p7 = bitalign(p7, p6, 32 - xxe);306 __CLC_FLOATN q0 = __CLC_AS_FLOATN(307 sign | ((127U - __CLC_AS_UINTN(xe + 23 + xxe)) << 23U) | p7 >> 9);308 309 // At this point, the fraction q1 + q0 is correct to at least 48 bits310 // Now we need to multiply the fraction by pi/2311 // This loses us about 4 bits312 // pi/2 = C90 FDA A22 168 C23 4C4313 314 const __CLC_FLOATN pio2h = (__CLC_FLOATN)0xc90fda / 0x1.0p+23f;315 const __CLC_FLOATN pio2hh = (__CLC_FLOATN)0xc90 / 0x1.0p+11f;316 const __CLC_FLOATN pio2ht = (__CLC_FLOATN)0xfda / 0x1.0p+23f;317 const __CLC_FLOATN pio2t = (__CLC_FLOATN)0xa22168 / 0x1.0p+47f;318 319 __CLC_FLOATN rh, rt;320 321 if (__CLC_HAVE_HW_FMA32()) {322 rh = q1 * pio2h;323 rt = __clc_fma(q0, pio2h, __clc_fma(q1, pio2t, __clc_fma(q1, pio2h, -rh)));324 } else {325 __CLC_FLOATN q1h = __CLC_AS_FLOATN(__CLC_AS_UINTN(q1) & 0xfffff000);326 __CLC_FLOATN q1t = q1 - q1h;327 rh = q1 * pio2h;328 rt = __clc_mad(329 q1t, pio2ht,330 __clc_mad(q1t, pio2hh,331 __clc_mad(q1h, pio2ht, __clc_mad(q1h, pio2hh, -rh))));332 rt = __clc_mad(q0, pio2h, __clc_mad(q1, pio2t, rt));333 }334 335 __CLC_FLOATN t = rh + rt;336 rt = rt - (t - rh);337 338 *r = t;339 *rr = rt;340 return ((i >> 1) + (i & 1)) & 0x3;341}342 343_CLC_DEF _CLC_OVERLOAD __CLC_INTN __clc_argReductionS(private __CLC_FLOATN *r,344 private __CLC_FLOATN *rr,345 __CLC_FLOATN x) {346 __CLC_INTN is_small = x < (__CLC_FLOATN)0x1.0p+23f;347#ifdef __CLC_SCALAR348 if (is_small)349 return __clc_argReductionSmallS(r, rr, x);350 else351 return __clc_argReductionLargeS(r, rr, x);352#else353 __CLC_FLOATN r1, rr1, r2, rr2;354 __CLC_INTN ret1 = __clc_argReductionSmallS(&r1, &rr1, x);355 __CLC_INTN ret2 = __clc_argReductionLargeS(&r2, &rr2, x);356 *r = is_small ? r1 : r2;357 *rr = is_small ? rr1 : rr2;358 return is_small ? ret1 : ret2;359#endif360}361