1/* 2 * Copyright (c) 2014 Advanced Micro Devices, Inc. 3 * 4 * Permission is hereby granted, free of charge, to any person obtaining a copy 5 * of this software and associated documentation files (the "Software"), to deal 6 * in the Software without restriction, including without limitation the rights 7 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 8 * copies of the Software, and to permit persons to whom the Software is 9 * furnished to do so, subject to the following conditions: 10 * 11 * The above copyright notice and this permission notice shall be included in 12 * all copies or substantial portions of the Software. 13 * 14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 17 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 18 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 19 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN 20 * THE SOFTWARE. 21 */ 22 23#include <clc/clc.h> 24 25#include "math.h" 26#include "tables.h" 27#include "sincos_helpers.h" 28 29#define bitalign(hi, lo, shift) \ 30 ((hi) << (32 - (shift))) | ((lo) >> (shift)); 31 32#define bytealign(src0, src1, src2) \ 33 ((uint) (((((long)(src0)) << 32) | (long)(src1)) >> (((src2) & 3)*8))) 34 35_CLC_DEF float __clc_sinf_piby4(float x, float y) { 36 // Taylor series for sin(x) is x - x^3/3! + x^5/5! - x^7/7! ... 37 // = x * (1 - x^2/3! + x^4/5! - x^6/7! ... 38 // = x * f(w) 39 // where w = x*x and f(w) = (1 - w/3! + w^2/5! - w^3/7! ... 40 // We use a minimax approximation of (f(w) - 1) / w 41 // because this produces an expansion in even powers of x. 42 43 const float c1 = -0.1666666666e0f; 44 const float c2 = 0.8333331876e-2f; 45 const float c3 = -0.198400874e-3f; 46 const float c4 = 0.272500015e-5f; 47 const float c5 = -2.5050759689e-08f; // 0xb2d72f34 48 const float c6 = 1.5896910177e-10f; // 0x2f2ec9d3 49 50 float z = x * x; 51 float v = z * x; 52 float r = mad(z, mad(z, mad(z, mad(z, c6, c5), c4), c3), c2); 53 float ret = x - mad(v, -c1, mad(z, mad(y, 0.5f, -v*r), -y)); 54 55 return ret; 56} 57 58_CLC_DEF float __clc_cosf_piby4(float x, float y) { 59 // Taylor series for cos(x) is 1 - x^2/2! + x^4/4! - x^6/6! ... 60 // = f(w) 61 // where w = x*x and f(w) = (1 - w/2! + w^2/4! - w^3/6! ... 62 // We use a minimax approximation of (f(w) - 1 + w/2) / (w*w) 63 // because this produces an expansion in even powers of x. 64 65 const float c1 = 0.416666666e-1f; 66 const float c2 = -0.138888876e-2f; 67 const float c3 = 0.248006008e-4f; 68 const float c4 = -0.2730101334e-6f; 69 const float c5 = 2.0875723372e-09f; // 0x310f74f6 70 const float c6 = -1.1359647598e-11f; // 0xad47d74e 71 72 float z = x * x; 73 float r = z * mad(z, mad(z, mad(z, mad(z, mad(z, c6, c5), c4), c3), c2), c1); 74 75 // if |x| < 0.3 76 float qx = 0.0f; 77 78 int ix = as_int(x) & EXSIGNBIT_SP32; 79 80 // 0.78125 > |x| >= 0.3 81 float xby4 = as_float(ix - 0x01000000); 82 qx = (ix >= 0x3e99999a) & (ix <= 0x3f480000) ? xby4 : qx; 83 84 // x > 0.78125 85 qx = ix > 0x3f480000 ? 0.28125f : qx; 86 87 float hz = mad(z, 0.5f, -qx); 88 float a = 1.0f - qx; 89 float ret = a - (hz - mad(z, r, -x*y)); 90 return ret; 91} 92 93_CLC_DEF float __clc_tanf_piby4(float x, int regn) 94{ 95 // Core Remez [1,2] approximation to tan(x) on the interval [0,pi/4]. 96 float r = x * x; 97 98 float a = mad(r, -0.0172032480471481694693109f, 0.385296071263995406715129f); 99 100 float b = mad(r, 101 mad(r, 0.01844239256901656082986661f, -0.51396505478854532132342f), 102 1.15588821434688393452299f); 103 104 float t = mad(x*r, native_divide(a, b), x); 105 float tr = -MATH_RECIP(t); 106 107 return regn & 1 ? tr : t; 108} 109 110_CLC_DEF void __clc_fullMulS(float *hi, float *lo, float a, float b, float bh, float bt) 111{ 112 if (HAVE_HW_FMA32()) { 113 float ph = a * b; 114 *hi = ph; 115 *lo = fma(a, b, -ph); 116 } else { 117 float ah = as_float(as_uint(a) & 0xfffff000U); 118 float at = a - ah; 119 float ph = a * b; 120 float pt = mad(at, bt, mad(at, bh, mad(ah, bt, mad(ah, bh, -ph)))); 121 *hi = ph; 122 *lo = pt; 123 } 124} 125 126_CLC_DEF float __clc_removePi2S(float *hi, float *lo, float x) 127{ 128 // 72 bits of pi/2 129 const float fpiby2_1 = (float) 0xC90FDA / 0x1.0p+23f; 130 const float fpiby2_1_h = (float) 0xC90 / 0x1.0p+11f; 131 const float fpiby2_1_t = (float) 0xFDA / 0x1.0p+23f; 132 133 const float fpiby2_2 = (float) 0xA22168 / 0x1.0p+47f; 134 const float fpiby2_2_h = (float) 0xA22 / 0x1.0p+35f; 135 const float fpiby2_2_t = (float) 0x168 / 0x1.0p+47f; 136 137 const float fpiby2_3 = (float) 0xC234C4 / 0x1.0p+71f; 138 const float fpiby2_3_h = (float) 0xC23 / 0x1.0p+59f; 139 const float fpiby2_3_t = (float) 0x4C4 / 0x1.0p+71f; 140 141 const float twobypi = 0x1.45f306p-1f; 142 143 float fnpi2 = trunc(mad(x, twobypi, 0.5f)); 144 145 // subtract n * pi/2 from x 146 float rhead, rtail; 147 __clc_fullMulS(&rhead, &rtail, fnpi2, fpiby2_1, fpiby2_1_h, fpiby2_1_t); 148 float v = x - rhead; 149 float rem = v + (((x - v) - rhead) - rtail); 150 151 float rhead2, rtail2; 152 __clc_fullMulS(&rhead2, &rtail2, fnpi2, fpiby2_2, fpiby2_2_h, fpiby2_2_t); 153 v = rem - rhead2; 154 rem = v + (((rem - v) - rhead2) - rtail2); 155 156 float rhead3, rtail3; 157 __clc_fullMulS(&rhead3, &rtail3, fnpi2, fpiby2_3, fpiby2_3_h, fpiby2_3_t); 158 v = rem - rhead3; 159 160 *hi = v + ((rem - v) - rhead3); 161 *lo = -rtail3; 162 return fnpi2; 163} 164 165_CLC_DEF int __clc_argReductionSmallS(float *r, float *rr, float x) 166{ 167 float fnpi2 = __clc_removePi2S(r, rr, x); 168 return (int)fnpi2 & 0x3; 169} 170 171#define FULL_MUL(A, B, HI, LO) \ 172 LO = A * B; \ 173 HI = mul_hi(A, B) 174 175#define FULL_MAD(A, B, C, HI, LO) \ 176 LO = ((A) * (B) + (C)); \ 177 HI = mul_hi(A, B); \ 178 HI += LO < C 179 180_CLC_DEF int __clc_argReductionLargeS(float *r, float *rr, float x) 181{ 182 int xe = (int)(as_uint(x) >> 23) - 127; 183 uint xm = 0x00800000U | (as_uint(x) & 0x7fffffU); 184 185 // 224 bits of 2/PI: . A2F9836E 4E441529 FC2757D1 F534DDC0 DB629599 3C439041 FE5163AB 186 const uint b6 = 0xA2F9836EU; 187 const uint b5 = 0x4E441529U; 188 const uint b4 = 0xFC2757D1U; 189 const uint b3 = 0xF534DDC0U; 190 const uint b2 = 0xDB629599U; 191 const uint b1 = 0x3C439041U; 192 const uint b0 = 0xFE5163ABU; 193 194 uint p0, p1, p2, p3, p4, p5, p6, p7, c0, c1; 195 196 FULL_MUL(xm, b0, c0, p0); 197 FULL_MAD(xm, b1, c0, c1, p1); 198 FULL_MAD(xm, b2, c1, c0, p2); 199 FULL_MAD(xm, b3, c0, c1, p3); 200 FULL_MAD(xm, b4, c1, c0, p4); 201 FULL_MAD(xm, b5, c0, c1, p5); 202 FULL_MAD(xm, b6, c1, p7, p6); 203 204 uint fbits = 224 + 23 - xe; 205 206 // shift amount to get 2 lsb of integer part at top 2 bits 207 // min: 25 (xe=18) max: 134 (xe=127) 208 uint shift = 256U - 2 - fbits; 209 210 // Shift by up to 134/32 = 4 words 211 int c = shift > 31; 212 p7 = c ? p6 : p7; 213 p6 = c ? p5 : p6; 214 p5 = c ? p4 : p5; 215 p4 = c ? p3 : p4; 216 p3 = c ? p2 : p3; 217 p2 = c ? p1 : p2; 218 p1 = c ? p0 : p1; 219 shift -= (-c) & 32; 220 221 c = shift > 31; 222 p7 = c ? p6 : p7; 223 p6 = c ? p5 : p6; 224 p5 = c ? p4 : p5; 225 p4 = c ? p3 : p4; 226 p3 = c ? p2 : p3; 227 p2 = c ? p1 : p2; 228 shift -= (-c) & 32; 229 230 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 shift -= (-c) & 32; 237 238 c = shift > 31; 239 p7 = c ? p6 : p7; 240 p6 = c ? p5 : p6; 241 p5 = c ? p4 : p5; 242 p4 = c ? p3 : p4; 243 shift -= (-c) & 32; 244 245 // bitalign cannot handle a shift of 32 246 c = shift > 0; 247 shift = 32 - shift; 248 uint t7 = bitalign(p7, p6, shift); 249 uint t6 = bitalign(p6, p5, shift); 250 uint t5 = bitalign(p5, p4, shift); 251 p7 = c ? t7 : p7; 252 p6 = c ? t6 : p6; 253 p5 = c ? t5 : p5; 254 255 // Get 2 lsb of int part and msb of fraction 256 int i = p7 >> 29; 257 258 // Scoot up 2 more bits so only fraction remains 259 p7 = bitalign(p7, p6, 30); 260 p6 = bitalign(p6, p5, 30); 261 p5 = bitalign(p5, p4, 30); 262 263 // Subtract 1 if msb of fraction is 1, i.e. fraction >= 0.5 264 uint flip = i & 1 ? 0xffffffffU : 0U; 265 uint sign = i & 1 ? 0x80000000U : 0U; 266 p7 = p7 ^ flip; 267 p6 = p6 ^ flip; 268 p5 = p5 ^ flip; 269 270 // Find exponent and shift away leading zeroes and hidden bit 271 xe = clz(p7) + 1; 272 shift = 32 - xe; 273 p7 = bitalign(p7, p6, shift); 274 p6 = bitalign(p6, p5, shift); 275 276 // Most significant part of fraction 277 float q1 = as_float(sign | ((127 - xe) << 23) | (p7 >> 9)); 278 279 // Shift out bits we captured on q1 280 p7 = bitalign(p7, p6, 32-23); 281 282 // Get 24 more bits of fraction in another float, there are not long strings of zeroes here 283 int xxe = clz(p7) + 1; 284 p7 = bitalign(p7, p6, 32-xxe); 285 float q0 = as_float(sign | ((127 - (xe + 23 + xxe)) << 23) | (p7 >> 9)); 286 287 // At this point, the fraction q1 + q0 is correct to at least 48 bits 288 // Now we need to multiply the fraction by pi/2 289 // This loses us about 4 bits 290 // pi/2 = C90 FDA A22 168 C23 4C4 291 292 const float pio2h = (float)0xc90fda / 0x1.0p+23f; 293 const float pio2hh = (float)0xc90 / 0x1.0p+11f; 294 const float pio2ht = (float)0xfda / 0x1.0p+23f; 295 const float pio2t = (float)0xa22168 / 0x1.0p+47f; 296 297 float rh, rt; 298 299 if (HAVE_HW_FMA32()) { 300 rh = q1 * pio2h; 301 rt = fma(q0, pio2h, fma(q1, pio2t, fma(q1, pio2h, -rh))); 302 } else { 303 float q1h = as_float(as_uint(q1) & 0xfffff000); 304 float q1t = q1 - q1h; 305 rh = q1 * pio2h; 306 rt = mad(q1t, pio2ht, mad(q1t, pio2hh, mad(q1h, pio2ht, mad(q1h, pio2hh, -rh)))); 307 rt = mad(q0, pio2h, mad(q1, pio2t, rt)); 308 } 309 310 float t = rh + rt; 311 rt = rt - (t - rh); 312 313 *r = t; 314 *rr = rt; 315 return ((i >> 1) + (i & 1)) & 0x3; 316} 317 318_CLC_DEF int __clc_argReductionS(float *r, float *rr, float x) 319{ 320 if (x < 0x1.0p+23f) 321 return __clc_argReductionSmallS(r, rr, x); 322 else 323 return __clc_argReductionLargeS(r, rr, x); 324} 325 326#ifdef cl_khr_fp64 327 328#pragma OPENCL EXTENSION cl_khr_fp64 : enable 329 330// Reduction for medium sized arguments 331_CLC_DEF void __clc_remainder_piby2_medium(double x, double *r, double *rr, int *regn) { 332 // How many pi/2 is x a multiple of? 333 const double two_by_pi = 0x1.45f306dc9c883p-1; 334 double dnpi2 = trunc(fma(x, two_by_pi, 0.5)); 335 336 const double piby2_h = -7074237752028440.0 / 0x1.0p+52; 337 const double piby2_m = -2483878800010755.0 / 0x1.0p+105; 338 const double piby2_t = -3956492004828932.0 / 0x1.0p+158; 339 340 // Compute product of npi2 with 159 bits of 2/pi 341 double p_hh = piby2_h * dnpi2; 342 double p_ht = fma(piby2_h, dnpi2, -p_hh); 343 double p_mh = piby2_m * dnpi2; 344 double p_mt = fma(piby2_m, dnpi2, -p_mh); 345 double p_th = piby2_t * dnpi2; 346 double p_tt = fma(piby2_t, dnpi2, -p_th); 347 348 // Reduce to 159 bits 349 double ph = p_hh; 350 double pm = p_ht + p_mh; 351 double t = p_mh - (pm - p_ht); 352 double pt = p_th + t + p_mt + p_tt; 353 t = ph + pm; pm = pm - (t - ph); ph = t; 354 t = pm + pt; pt = pt - (t - pm); pm = t; 355 356 // Subtract from x 357 t = x + ph; 358 double qh = t + pm; 359 double qt = pm - (qh - t) + pt; 360 361 *r = qh; 362 *rr = qt; 363 *regn = (int)(long)dnpi2 & 0x3; 364} 365 366// Given positive argument x, reduce it to the range [-pi/4,pi/4] using 367// extra precision, and return the result in r, rr. 368// Return value "regn" tells how many lots of pi/2 were subtracted 369// from x to put it in the range [-pi/4,pi/4], mod 4. 370 371_CLC_DEF void __clc_remainder_piby2_large(double x, double *r, double *rr, int *regn) { 372 373 long ux = as_long(x); 374 int e = (int)(ux >> 52) - 1023; 375 int i = max(23, (e >> 3) + 17); 376 int j = 150 - i; 377 int j16 = j & ~0xf; 378 double fract_temp; 379 380 // The following extracts 192 consecutive bits of 2/pi aligned on an arbitrary byte boundary 381 uint4 q0 = USE_TABLE(pibits_tbl, j16); 382 uint4 q1 = USE_TABLE(pibits_tbl, (j16 + 16)); 383 uint4 q2 = USE_TABLE(pibits_tbl, (j16 + 32)); 384 385 int k = (j >> 2) & 0x3; 386 int4 c = (int4)k == (int4)(0, 1, 2, 3); 387 388 uint u0, u1, u2, u3, u4, u5, u6; 389 390 u0 = c.s1 ? q0.s1 : q0.s0; 391 u0 = c.s2 ? q0.s2 : u0; 392 u0 = c.s3 ? q0.s3 : u0; 393 394 u1 = c.s1 ? q0.s2 : q0.s1; 395 u1 = c.s2 ? q0.s3 : u1; 396 u1 = c.s3 ? q1.s0 : u1; 397 398 u2 = c.s1 ? q0.s3 : q0.s2; 399 u2 = c.s2 ? q1.s0 : u2; 400 u2 = c.s3 ? q1.s1 : u2; 401 402 u3 = c.s1 ? q1.s0 : q0.s3; 403 u3 = c.s2 ? q1.s1 : u3; 404 u3 = c.s3 ? q1.s2 : u3; 405 406 u4 = c.s1 ? q1.s1 : q1.s0; 407 u4 = c.s2 ? q1.s2 : u4; 408 u4 = c.s3 ? q1.s3 : u4; 409 410 u5 = c.s1 ? q1.s2 : q1.s1; 411 u5 = c.s2 ? q1.s3 : u5; 412 u5 = c.s3 ? q2.s0 : u5; 413 414 u6 = c.s1 ? q1.s3 : q1.s2; 415 u6 = c.s2 ? q2.s0 : u6; 416 u6 = c.s3 ? q2.s1 : u6; 417 418 uint v0 = bytealign(u1, u0, j); 419 uint v1 = bytealign(u2, u1, j); 420 uint v2 = bytealign(u3, u2, j); 421 uint v3 = bytealign(u4, u3, j); 422 uint v4 = bytealign(u5, u4, j); 423 uint v5 = bytealign(u6, u5, j); 424 425 // Place those 192 bits in 4 48-bit doubles along with correct exponent 426 // If i > 1018 we would get subnormals so we scale p up and x down to get the same product 427 i = 2 + 8*i; 428 x *= i > 1018 ? 0x1.0p-136 : 1.0; 429 i -= i > 1018 ? 136 : 0; 430 431 uint ua = (uint)(1023 + 52 - i) << 20; 432 double a = as_double((uint2)(0, ua)); 433 double p0 = as_double((uint2)(v0, ua | (v1 & 0xffffU))) - a; 434 ua += 0x03000000U; 435 a = as_double((uint2)(0, ua)); 436 double p1 = as_double((uint2)((v2 << 16) | (v1 >> 16), ua | (v2 >> 16))) - a; 437 ua += 0x03000000U; 438 a = as_double((uint2)(0, ua)); 439 double p2 = as_double((uint2)(v3, ua | (v4 & 0xffffU))) - a; 440 ua += 0x03000000U; 441 a = as_double((uint2)(0, ua)); 442 double p3 = as_double((uint2)((v5 << 16) | (v4 >> 16), ua | (v5 >> 16))) - a; 443 444 // Exact multiply 445 double f0h = p0 * x; 446 double f0l = fma(p0, x, -f0h); 447 double f1h = p1 * x; 448 double f1l = fma(p1, x, -f1h); 449 double f2h = p2 * x; 450 double f2l = fma(p2, x, -f2h); 451 double f3h = p3 * x; 452 double f3l = fma(p3, x, -f3h); 453 454 // Accumulate product into 4 doubles 455 double s, t; 456 457 double f3 = f3h + f2h; 458 t = f2h - (f3 - f3h); 459 s = f3l + t; 460 t = t - (s - f3l); 461 462 double f2 = s + f1h; 463 t = f1h - (f2 - s) + t; 464 s = f2l + t; 465 t = t - (s - f2l); 466 467 double f1 = s + f0h; 468 t = f0h - (f1 - s) + t; 469 s = f1l + t; 470 471 double f0 = s + f0l; 472 473 // Strip off unwanted large integer bits 474 f3 = 0x1.0p+10 * fract(f3 * 0x1.0p-10, &fract_temp); 475 f3 += f3 + f2 < 0.0 ? 0x1.0p+10 : 0.0; 476 477 // Compute least significant integer bits 478 t = f3 + f2; 479 double di = t - fract(t, &fract_temp); 480 i = (float)di; 481 482 // Shift out remaining integer part 483 f3 -= di; 484 s = f3 + f2; t = f2 - (s - f3); f3 = s; f2 = t; 485 s = f2 + f1; t = f1 - (s - f2); f2 = s; f1 = t; 486 f1 += f0; 487 488 // Subtract 1 if fraction is >= 0.5, and update regn 489 int g = f3 >= 0.5; 490 i += g; 491 f3 -= (float)g; 492 493 // Shift up bits 494 s = f3 + f2; t = f2 -(s - f3); f3 = s; f2 = t + f1; 495 496 // Multiply precise fraction by pi/2 to get radians 497 const double p2h = 7074237752028440.0 / 0x1.0p+52; 498 const double p2t = 4967757600021510.0 / 0x1.0p+106; 499 500 double rhi = f3 * p2h; 501 double rlo = fma(f2, p2h, fma(f3, p2t, fma(f3, p2h, -rhi))); 502 503 *r = rhi + rlo; 504 *rr = rlo - (*r - rhi); 505 *regn = i & 0x3; 506} 507 508 509_CLC_DEF double2 __clc_sincos_piby4(double x, double xx) { 510 // Taylor series for sin(x) is x - x^3/3! + x^5/5! - x^7/7! ... 511 // = x * (1 - x^2/3! + x^4/5! - x^6/7! ... 512 // = x * f(w) 513 // where w = x*x and f(w) = (1 - w/3! + w^2/5! - w^3/7! ... 514 // We use a minimax approximation of (f(w) - 1) / w 515 // because this produces an expansion in even powers of x. 516 // If xx (the tail of x) is non-zero, we add a correction 517 // term g(x,xx) = (1-x*x/2)*xx to the result, where g(x,xx) 518 // is an approximation to cos(x)*sin(xx) valid because 519 // xx is tiny relative to x. 520 521 // Taylor series for cos(x) is 1 - x^2/2! + x^4/4! - x^6/6! ... 522 // = f(w) 523 // where w = x*x and f(w) = (1 - w/2! + w^2/4! - w^3/6! ... 524 // We use a minimax approximation of (f(w) - 1 + w/2) / (w*w) 525 // because this produces an expansion in even powers of x. 526 // If xx (the tail of x) is non-zero, we subtract a correction 527 // term g(x,xx) = x*xx to the result, where g(x,xx) 528 // is an approximation to sin(x)*sin(xx) valid because 529 // xx is tiny relative to x. 530 531 const double sc1 = -0.166666666666666646259241729; 532 const double sc2 = 0.833333333333095043065222816e-2; 533 const double sc3 = -0.19841269836761125688538679e-3; 534 const double sc4 = 0.275573161037288022676895908448e-5; 535 const double sc5 = -0.25051132068021699772257377197e-7; 536 const double sc6 = 0.159181443044859136852668200e-9; 537 538 const double cc1 = 0.41666666666666665390037e-1; 539 const double cc2 = -0.13888888888887398280412e-2; 540 const double cc3 = 0.248015872987670414957399e-4; 541 const double cc4 = -0.275573172723441909470836e-6; 542 const double cc5 = 0.208761463822329611076335e-8; 543 const double cc6 = -0.113826398067944859590880e-10; 544 545 double x2 = x * x; 546 double x3 = x2 * x; 547 double r = 0.5 * x2; 548 double t = 1.0 - r; 549 550 double sp = fma(fma(fma(fma(sc6, x2, sc5), x2, sc4), x2, sc3), x2, sc2); 551 552 double cp = t + fma(fma(fma(fma(fma(fma(cc6, x2, cc5), x2, cc4), x2, cc3), x2, cc2), x2, cc1), 553 x2*x2, fma(x, xx, (1.0 - t) - r)); 554 555 double2 ret; 556 ret.lo = x - fma(-x3, sc1, fma(fma(-x3, sp, 0.5*xx), x2, -xx)); 557 ret.hi = cp; 558 559 return ret; 560} 561 562#endif 563