359 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 void __clc_sincos_piby4(__CLC_DOUBLEN x,10 __CLC_DOUBLEN xx,11 private __CLC_DOUBLEN *sinval,12 private __CLC_DOUBLEN *cosval) {13 // Taylor series for sin(x) is x - x^3/3! + x^5/5! - x^7/7! ...14 // = x * (1 - x^2/3! + x^4/5! - x^6/7! ...15 // = x * f(w)16 // where w = x*x and f(w) = (1 - w/3! + w^2/5! - w^3/7! ...17 // We use a minimax approximation of (f(w) - 1) / w18 // because this produces an expansion in even powers of x.19 // If xx (the tail of x) is non-zero, we add a correction20 // term g(x,xx) = (1-x*x/2)*xx to the result, where g(x,xx)21 // is an approximation to cos(x)*sin(xx) valid because22 // xx is tiny relative to x.23 24 // Taylor series for cos(x) is 1 - x^2/2! + x^4/4! - x^6/6! ...25 // = f(w)26 // where w = x*x and f(w) = (1 - w/2! + w^2/4! - w^3/6! ...27 // We use a minimax approximation of (f(w) - 1 + w/2) / (w*w)28 // because this produces an expansion in even powers of x.29 // If xx (the tail of x) is non-zero, we subtract a correction30 // term g(x,xx) = x*xx to the result, where g(x,xx)31 // is an approximation to sin(x)*sin(xx) valid because32 // xx is tiny relative to x.33 34 const __CLC_DOUBLEN sc1 = -0.166666666666666646259241729;35 const __CLC_DOUBLEN sc2 = 0.833333333333095043065222816e-2;36 const __CLC_DOUBLEN sc3 = -0.19841269836761125688538679e-3;37 const __CLC_DOUBLEN sc4 = 0.275573161037288022676895908448e-5;38 const __CLC_DOUBLEN sc5 = -0.25051132068021699772257377197e-7;39 const __CLC_DOUBLEN sc6 = 0.159181443044859136852668200e-9;40 41 const __CLC_DOUBLEN cc1 = 0.41666666666666665390037e-1;42 const __CLC_DOUBLEN cc2 = -0.13888888888887398280412e-2;43 const __CLC_DOUBLEN cc3 = 0.248015872987670414957399e-4;44 const __CLC_DOUBLEN cc4 = -0.275573172723441909470836e-6;45 const __CLC_DOUBLEN cc5 = 0.208761463822329611076335e-8;46 const __CLC_DOUBLEN cc6 = -0.113826398067944859590880e-10;47 48 __CLC_DOUBLEN x2 = x * x;49 __CLC_DOUBLEN x3 = x2 * x;50 __CLC_DOUBLEN r = (__CLC_DOUBLEN)0.5 * x2;51 __CLC_DOUBLEN t = (__CLC_DOUBLEN)1.0 - r;52 53 __CLC_DOUBLEN sp = __clc_fma(54 __clc_fma(__clc_fma(__clc_fma(sc6, x2, sc5), x2, sc4), x2, sc3), x2, sc2);55 56 __CLC_DOUBLEN cp =57 t +58 __clc_fma(__clc_fma(__clc_fma(__clc_fma(__clc_fma(__clc_fma(cc6, x2, cc5),59 x2, cc4),60 x2, cc3),61 x2, cc2),62 x2, cc1),63 x2 * x2, __clc_fma(x, xx, (1.0 - t) - r));64 65 *sinval =66 x - __clc_fma(-x3, sc1, __clc_fma(__clc_fma(-x3, sp, 0.5 * xx), x2, -xx));67 *cosval = cp;68}69 70_CLC_DEF _CLC_OVERLOAD void __clc_tan_piby4(__CLC_DOUBLEN x, __CLC_DOUBLEN xx,71 private __CLC_DOUBLEN *leadval,72 private __CLC_DOUBLEN *tailval) {73 // 0x3fe921fb54442d1874 const __CLC_DOUBLEN piby4_lead = 7.85398163397448278999e-01;75 // 0x3c81a62633145c0676 const __CLC_DOUBLEN piby4_tail = 3.06161699786838240164e-17;77 78 // In order to maintain relative precision transform using the identity:79 // tan(pi/4-x) = (1-tan(x))/(1+tan(x)) for arguments close to pi/4.80 // Similarly use tan(x-pi/4) = (tan(x)-1)/(tan(x)+1) close to -pi/4.81 82 __CLC_LONGN ca = x > 0.68;83 __CLC_LONGN cb = x < -0.68;84 __CLC_DOUBLEN transform = ca ? 1.0 : 0.0;85 transform = cb ? -1.0 : transform;86 87 __CLC_DOUBLEN tx = __clc_fma(-transform, x, piby4_lead) +88 __clc_fma(-transform, xx, piby4_tail);89 __CLC_LONGN c = ca | cb;90 x = c ? tx : x;91 xx = c ? 0.0 : xx;92 93 // Core Remez [2,3] approximation to tan(x+xx) on the interval [0,0.68].94 __CLC_DOUBLEN t1 = x;95 __CLC_DOUBLEN r = __clc_fma(2.0, x * xx, x * x);96 97 __CLC_DOUBLEN a = __clc_fma(r,98 __clc_fma(r, 0.224044448537022097264602535574e-3,99 -0.229345080057565662883358588111e-1),100 0.372379159759792203640806338901e0);101 102 __CLC_DOUBLEN b =103 __clc_fma(r,104 __clc_fma(r,105 __clc_fma(r, -0.232371494088563558304549252913e-3,106 0.260656620398645407524064091208e-1),107 -0.515658515729031149329237816945e0),108 0.111713747927937668539901657944e1);109 110 __CLC_DOUBLEN t2 = __clc_fma(MATH_DIVIDE(a, b), x * r, xx);111 112 __CLC_DOUBLEN tp = t1 + t2;113 114 // Compute -1.0/(t1 + t2) accurately115 __CLC_DOUBLEN z1 =116 __CLC_AS_GENTYPE(__CLC_AS_ULONGN(tp) & 0xffffffff00000000L);117 __CLC_DOUBLEN z2 = t2 - (z1 - t1);118 __CLC_DOUBLEN trec = -MATH_RECIP(tp);119 __CLC_DOUBLEN trec_top =120 __CLC_AS_GENTYPE(__CLC_AS_ULONGN(trec) & 0xffffffff00000000L);121 122 __CLC_DOUBLEN tpr = __clc_fma(123 __clc_fma(trec_top, z2, __clc_fma(trec_top, z1, 1.0)), trec, trec_top);124 125 __CLC_DOUBLEN tpt = transform * (1.0 - MATH_DIVIDE(2.0 * tp, 1.0 + tp));126 __CLC_DOUBLEN tptr = transform * (MATH_DIVIDE(2.0 * tp, tp - 1.0) - 1.0);127 128 *leadval = c ? tpt : tp;129 *tailval = c ? tptr : tpr;130}131 132// Reduction for medium sized arguments133_CLC_DEF _CLC_OVERLOAD void134__clc_remainder_piby2_medium(__CLC_DOUBLEN x, private __CLC_DOUBLEN *r,135 private __CLC_DOUBLEN *rr,136 private __CLC_INTN *regn) {137 // How many pi/2 is x a multiple of?138 const __CLC_DOUBLEN two_by_pi = 0x1.45f306dc9c883p-1;139 __CLC_DOUBLEN dnpi2 = __clc_trunc(__clc_fma(x, two_by_pi, 0.5));140 141 const __CLC_DOUBLEN piby2_h = -7074237752028440.0 / 0x1.0p+52;142 const __CLC_DOUBLEN piby2_m = -2483878800010755.0 / 0x1.0p+105;143 const __CLC_DOUBLEN piby2_t = -3956492004828932.0 / 0x1.0p+158;144 145 // Compute product of npi2 with 159 bits of 2/pi146 __CLC_DOUBLEN p_hh = piby2_h * dnpi2;147 __CLC_DOUBLEN p_ht = __clc_fma(piby2_h, dnpi2, -p_hh);148 __CLC_DOUBLEN p_mh = piby2_m * dnpi2;149 __CLC_DOUBLEN p_mt = __clc_fma(piby2_m, dnpi2, -p_mh);150 __CLC_DOUBLEN p_th = piby2_t * dnpi2;151 __CLC_DOUBLEN p_tt = __clc_fma(piby2_t, dnpi2, -p_th);152 153 // Reduce to 159 bits154 __CLC_DOUBLEN ph = p_hh;155 __CLC_DOUBLEN pm = p_ht + p_mh;156 __CLC_DOUBLEN t = p_mh - (pm - p_ht);157 __CLC_DOUBLEN pt = p_th + t + p_mt + p_tt;158 t = ph + pm;159 pm = pm - (t - ph);160 ph = t;161 t = pm + pt;162 pt = pt - (t - pm);163 pm = t;164 165 // Subtract from x166 t = x + ph;167 __CLC_DOUBLEN qh = t + pm;168 __CLC_DOUBLEN qt = pm - (qh - t) + pt;169 170 *r = qh;171 *rr = qt;172 *regn = __CLC_CONVERT_INTN(__CLC_CONVERT_LONGN(dnpi2)) & 0x3;173}174 175// Given positive argument x, reduce it to the range [-pi/4,pi/4] using176// extra precision, and return the result in r, rr.177// Return value "regn" tells how many lots of pi/2 were subtracted178// from x to put it in the range [-pi/4,pi/4], mod 4.179_CLC_DEF _CLC_OVERLOAD void180__clc_remainder_piby2_large(__CLC_DOUBLEN x, private __CLC_DOUBLEN *r,181 private __CLC_DOUBLEN *rr,182 private __CLC_INTN *regn) {183 184 __CLC_LONGN ux = __CLC_AS_LONGN(x);185 __CLC_INTN e = __CLC_CONVERT_INTN(ux >> EXPSHIFTBITS_DP64) - EXPBIAS_DP64;186 __CLC_INTN i = __clc_max(23, (e >> 3) + 17);187 __CLC_INTN j = 150 - i;188 __CLC_INTN j16 = j & ~0xf;189 __CLC_DOUBLEN fract_temp;190 191 // The following extracts 192 consecutive bits of 2/pi aligned on an arbitrary192 // byte boundary193 __CLC_ULONGN q0 = __CLC_USE_TABLE(pibits_tbl, j16);194 __CLC_ULONGN q1 = __CLC_USE_TABLE(pibits_tbl, (j16 + 8));195 __CLC_ULONGN q2 = __CLC_USE_TABLE(pibits_tbl, (j16 + 16));196 __CLC_ULONGN q3 = __CLC_USE_TABLE(pibits_tbl, (j16 + 24));197 __CLC_ULONGN q4 = __CLC_USE_TABLE(pibits_tbl, (j16 + 32));198 199 __CLC_UINTN q0s0 = __CLC_CONVERT_UINTN(q0);200 __CLC_UINTN q0s1 = __CLC_CONVERT_UINTN(q0 >> 32);201 202 __CLC_UINTN q1s0 = __CLC_CONVERT_UINTN(q1);203 __CLC_UINTN q1s1 = __CLC_CONVERT_UINTN(q1 >> 32);204 205 __CLC_UINTN q2s0 = __CLC_CONVERT_UINTN(q2);206 __CLC_UINTN q2s1 = __CLC_CONVERT_UINTN(q2 >> 32);207 208 __CLC_UINTN q3s0 = __CLC_CONVERT_UINTN(q3);209 __CLC_UINTN q3s1 = __CLC_CONVERT_UINTN(q3 >> 32);210 211 __CLC_UINTN q4s0 = __CLC_CONVERT_UINTN(q4);212 __CLC_UINTN q4s1 = __CLC_CONVERT_UINTN(q4 >> 32);213 214 __CLC_INTN k = (j >> 2) & 0x3;215 __CLC_INTN c0 = k == 0;216 __CLC_INTN c1 = k == 1;217 __CLC_INTN c2 = k == 2;218 __CLC_INTN c3 = k == 3;219 220 __CLC_UINTN u0, u1, u2, u3, u4, u5, u6;221 222 u0 = c1 ? q0s1 : q0s0;223 u0 = c2 ? q1s0 : u0;224 u0 = c3 ? q1s1 : u0;225 226 u1 = c1 ? q1s0 : q0s1;227 u1 = c2 ? q1s1 : u1;228 u1 = c3 ? q2s0 : u1;229 230 u2 = c1 ? q1s1 : q1s0;231 u2 = c2 ? q2s0 : u2;232 u2 = c3 ? q2s1 : u2;233 234 u3 = c1 ? q2s0 : q1s1;235 u3 = c2 ? q2s1 : u3;236 u3 = c3 ? q3s0 : u3;237 238 u4 = c1 ? q2s1 : q2s0;239 u4 = c2 ? q3s0 : u4;240 u4 = c3 ? q3s1 : u4;241 242 u5 = c1 ? q3s0 : q2s1;243 u5 = c2 ? q3s1 : u5;244 u5 = c3 ? q4s0 : u5;245 246 u6 = c1 ? q3s1 : q3s0;247 u6 = c2 ? q4s0 : u6;248 u6 = c3 ? q4s1 : u6;249 250 __CLC_UINTN v0 = bytealign(u1, u0, j);251 __CLC_UINTN v1 = bytealign(u2, u1, j);252 __CLC_UINTN v2 = bytealign(u3, u2, j);253 __CLC_UINTN v3 = bytealign(u4, u3, j);254 __CLC_UINTN v4 = bytealign(u5, u4, j);255 __CLC_UINTN v5 = bytealign(u6, u5, j);256 257 // Place those 192 bits in 4 48-bit doubles along with correct exponent258 // If i > 1018 we would get subnormals so we scale p up and x down to get the259 // same product260 i = 2 + 8 * i;261 x *= __CLC_CONVERT_BIT_INTN(i > 1018) ? 0x1.0p-136 : 1.0;262 i -= i > 1018 ? 136 : 0;263 264#define doublen_lohi(x, y) \265 __CLC_AS_DOUBLEN(__CLC_CONVERT_ULONGN((x)) & 0xFFFFFFFF | \266 __CLC_CONVERT_ULONGN((y)) << 32)267 268 __CLC_UINTN ua = __CLC_CONVERT_UINTN(EXPBIAS_DP64 + EXPSHIFTBITS_DP64 - i)269 << 20;270 __CLC_DOUBLEN a = doublen_lohi((__CLC_ULONGN)0, ua);271 __CLC_DOUBLEN p0 = doublen_lohi(v0, ua | (v1 & 0xffffU)) - a;272 ua += 0x03000000U;273 a = doublen_lohi((__CLC_ULONGN)0, ua);274 __CLC_DOUBLEN p1 =275 doublen_lohi(((v2 << 16) | (v1 >> 16)), (ua | (v2 >> 16))) - a;276 ua += 0x03000000U;277 a = doublen_lohi((__CLC_ULONGN)0, ua);278 __CLC_DOUBLEN p2 = doublen_lohi(v3, (ua | (v4 & 0xffffU))) - a;279 ua += 0x03000000U;280 a = doublen_lohi((__CLC_ULONGN)0, ua);281 __CLC_DOUBLEN p3 =282 doublen_lohi(((v5 << 16) | (v4 >> 16)), (ua | (v5 >> 16))) - a;283 284#undef doublen_lohi285 286 // Exact multiply287 __CLC_DOUBLEN f0h = p0 * x;288 __CLC_DOUBLEN f0l = __clc_fma(p0, x, -f0h);289 __CLC_DOUBLEN f1h = p1 * x;290 __CLC_DOUBLEN f1l = __clc_fma(p1, x, -f1h);291 __CLC_DOUBLEN f2h = p2 * x;292 __CLC_DOUBLEN f2l = __clc_fma(p2, x, -f2h);293 __CLC_DOUBLEN f3h = p3 * x;294 __CLC_DOUBLEN f3l = __clc_fma(p3, x, -f3h);295 296 // Accumulate product into 4 doubles297 __CLC_DOUBLEN s, t;298 299 __CLC_DOUBLEN f3 = f3h + f2h;300 t = f2h - (f3 - f3h);301 s = f3l + t;302 t = t - (s - f3l);303 304 __CLC_DOUBLEN f2 = s + f1h;305 t = f1h - (f2 - s) + t;306 s = f2l + t;307 t = t - (s - f2l);308 309 __CLC_DOUBLEN f1 = s + f0h;310 t = f0h - (f1 - s) + t;311 s = f1l + t;312 313 __CLC_DOUBLEN f0 = s + f0l;314 315 // Strip off unwanted large integer bits316 f3 = 0x1.0p+10 * __clc_fract(f3 * 0x1.0p-10, &fract_temp);317 f3 += f3 + f2 < 0.0 ? 0x1.0p+10 : 0.0;318 319 // Compute least significant integer bits320 t = f3 + f2;321 __CLC_DOUBLEN di = t - __clc_fract(t, &fract_temp);322 i = __CLC_CONVERT_INTN(__CLC_CONVERT_FLOATN(di));323 324 // Shift out remaining integer part325 f3 -= di;326 s = f3 + f2;327 t = f2 - (s - f3);328 f3 = s;329 f2 = t;330 s = f2 + f1;331 t = f1 - (s - f2);332 f2 = s;333 f1 = t;334 f1 += f0;335 336 // Subtract 1 if fraction is >= 0.5, and update regn337 __CLC_INTN g = __CLC_CONVERT_INTN(f3 >= 0.5 ? 1L : 0L);338 i += g;339 f3 -= __CLC_CONVERT_DOUBLEN(__CLC_CONVERT_FLOATN(g));340 341 // Shift up bits342 s = f3 + f2;343 t = f2 - (s - f3);344 f3 = s;345 f2 = t + f1;346 347 // Multiply precise fraction by pi/2 to get radians348 const __CLC_DOUBLEN p2h = 7074237752028440.0 / 0x1.0p+52;349 const __CLC_DOUBLEN p2t = 4967757600021510.0 / 0x1.0p+106;350 351 __CLC_DOUBLEN rhi = f3 * p2h;352 __CLC_DOUBLEN rlo =353 __clc_fma(f2, p2h, __clc_fma(f3, p2t, __clc_fma(f3, p2h, -rhi)));354 355 *r = rhi + rlo;356 *rr = rlo - (*r - rhi);357 *regn = i & 0x3;358}359