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 "ep_log.h" 24#include <clc/clc.h> 25#include <clc/clcmacro.h> 26#include <clc/math/math.h> 27 28_CLC_OVERLOAD _CLC_DEF float acosh(float x) { 29 uint ux = as_uint(x); 30 31 // Arguments greater than 1/sqrt(epsilon) in magnitude are 32 // approximated by acosh(x) = ln(2) + ln(x) 33 // For 2.0 <= x <= 1/sqrt(epsilon) the approximation is 34 // acosh(x) = ln(x + sqrt(x*x-1)) */ 35 int high = ux > 0x46000000U; 36 int med = ux > 0x40000000U; 37 38 float w = x - 1.0f; 39 float s = w*w + 2.0f*w; 40 float t = x*x - 1.0f; 41 float r = sqrt(med ? t : s) + (med ? x : w); 42 float v = (high ? x : r) - (med ? 1.0f : 0.0f); 43 float z = log1p(v) + (high ? 0x1.62e430p-1f : 0.0f); 44 45 z = ux >= PINFBITPATT_SP32 ? x : z; 46 z = x < 1.0f ? as_float(QNANBITPATT_SP32) : z; 47 48 return z; 49} 50 51_CLC_UNARY_VECTORIZE(_CLC_OVERLOAD _CLC_DEF, float, acosh, float) 52 53#ifdef cl_khr_fp64 54#pragma OPENCL EXTENSION cl_khr_fp64 : enable 55 56_CLC_OVERLOAD _CLC_DEF double acosh(double x) { 57 const double recrteps = 0x1.6a09e667f3bcdp+26; // 1/sqrt(eps) = 9.49062656242515593767e+07 58 //log2_lead and log2_tail sum to an extra-precise version of log(2) 59 const double log2_lead = 0x1.62e42ep-1; 60 const double log2_tail = 0x1.efa39ef35793cp-25; 61 62 // Handle x >= 128 here 63 int xlarge = x > recrteps; 64 double r = x + sqrt(fma(x, x, -1.0)); 65 r = xlarge ? x : r; 66 67 int xexp; 68 double r1, r2; 69 __clc_ep_log(r, &xexp, &r1, &r2); 70 71 double dxexp = xexp + xlarge; 72 r1 = fma(dxexp, log2_lead, r1); 73 r2 = fma(dxexp, log2_tail, r2); 74 75 double ret1 = r1 + r2; 76 77 // Handle 1 < x < 128 here 78 // We compute the value 79 // t = x - 1.0 + sqrt(2.0*(x - 1.0) + (x - 1.0)*(x - 1.0)) 80 // using simulated quad precision. 81 double t = x - 1.0; 82 double u1 = t * 2.0; 83 84 // (t,0) * (t,0) -> (v1, v2) 85 double v1 = t * t; 86 double v2 = fma(t, t, -v1); 87 88 // (u1,0) + (v1,v2) -> (w1,w2) 89 r = u1 + v1; 90 double s = (((u1 - r) + v1) + v2); 91 double w1 = r + s; 92 double w2 = (r - w1) + s; 93 94 // sqrt(w1,w2) -> (u1,u2) 95 double p1 = sqrt(w1); 96 double a1 = p1*p1; 97 double a2 = fma(p1, p1, -a1); 98 double temp = (((w1 - a1) - a2) + w2); 99 double p2 = MATH_DIVIDE(temp * 0.5, p1); 100 u1 = p1 + p2; 101 double u2 = (p1 - u1) + p2; 102 103 // (u1,u2) + (t,0) -> (r1,r2) 104 r = u1 + t; 105 s = ((u1 - r) + t) + u2; 106 // r1 = r + s; 107 // r2 = (r - r1) + s; 108 // t = r1 + r2; 109 t = r + s; 110 111 // For arguments 1.13 <= x <= 1.5 the log1p function is good enough 112 double ret2 = log1p(t); 113 114 ulong ux = as_ulong(x); 115 double ret = x >= 128.0 ? ret1 : ret2; 116 117 ret = ux >= 0x7FF0000000000000 ? x : ret; 118 ret = x == 1.0 ? 0.0 : ret; 119 ret = ((ux & SIGNBIT_DP64) != 0UL | x < 1.0) ? as_double(QNANBITPATT_DP64) : ret; 120 121 return ret; 122} 123 124_CLC_UNARY_VECTORIZE(_CLC_OVERLOAD _CLC_DEF, double, acosh, double) 125 126#endif 127 128#ifdef cl_khr_fp16 129 130#pragma OPENCL EXTENSION cl_khr_fp16 : enable 131 132_CLC_DEFINE_UNARY_BUILTIN_FP16(acosh) 133 134#endif 135