Lines Matching +full:range +full:- +full:double
2 * Double-precision x^y function.
4 * Copyright (c) 2018-2024, Arm Limited.
5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
20 #define SmallExp 0x3c9 /* top12(0x1p-54). */
22 #define ThresExp 0x03f /* BigExp - SmallExp. */
30 #define SmallPowX 0x001 /* top12(0x1p-126). */
32 #define ThresPowX 0x7fe /* BigPowX - SmallPowX. */
33 #define SmallPowY 0x3be /* top12(0x1.e7b6p-65). */
35 #define ThresPowY 0x080 /* BigPowY - SmallPowY. */
37 /* Top 12 bits of a double (sign and exponent bits). */
39 top12 (double x) in top12()
46 normalized in the subnormal range using the sign bit for the exponent. */
47 static inline double
48 log_inline (uint64_t ix, double *tail) in log_inline()
50 /* x = 2^k z; where z is in range [Off,2*Off) and exact. in log_inline()
51 The range is split into N subintervals. in log_inline()
53 uint64_t tmp = ix - Off; in log_inline()
54 int i = (tmp >> (52 - V_POW_LOG_TABLE_BITS)) & (N_LOG - 1); in log_inline()
56 uint64_t iz = ix - (tmp & 0xfffULL << 52); in log_inline()
57 double z = asdouble (iz); in log_inline()
58 double kd = (double) k; in log_inline()
60 /* log(x) = k*Ln2 + log(c) + log1p(z/c-1). */ in log_inline()
61 double invc = __v_pow_log_data.invc[i]; in log_inline()
62 double logc = __v_pow_log_data.logc[i]; in log_inline()
63 double logctail = __v_pow_log_data.logctail[i]; in log_inline()
66 |z/c - 1| < 1/N, so r = z/c - 1 is exactly representible. */ in log_inline()
67 double r = fma (z, invc, -1.0); in log_inline()
70 double t1 = kd * __v_pow_log_data.ln2_hi + logc; in log_inline()
71 double t2 = t1 + r; in log_inline()
72 double lo1 = kd * __v_pow_log_data.ln2_lo + logctail; in log_inline()
73 double lo2 = t1 - t2 + r; in log_inline()
76 double ar = As[0] * r; in log_inline()
77 double ar2 = r * ar; in log_inline()
78 double ar3 = r * ar2; in log_inline()
80 double hi = t2 + ar2; in log_inline()
81 double lo3 = fma (ar, r, -ar2); in log_inline()
82 double lo4 = t2 - hi + ar2; in log_inline()
83 /* p = log1p(r) - r - A[0]*r*r. */ in log_inline()
84 double p = (ar3 in log_inline()
87 double lo = lo1 + lo2 + lo3 + lo4 + p; in log_inline()
88 double y = hi + lo; in log_inline()
89 *tail = hi - y + lo; in log_inline()
97 a double. (int32_t)KI is the k used in the argument reduction and exponent
100 static inline double
101 special_case (double tmp, uint64_t sbits, uint64_t ki) in special_case()
103 double scale, y; in special_case()
108 sbits -= 1009ull << 52; in special_case()
113 /* k < 0, need special care in the subnormal range. */ in special_case()
122 range to avoid double rounding that can cause 0.5+E/2 ulp error where in special_case()
123 E is the worst-case ulp error outside the subnormal range. So this in special_case()
124 is only useful if the goal is better than 1 ulp worst-case error. */ in special_case()
125 double hi, lo, one = 1.0; in special_case()
127 one = -1.0; in special_case()
128 lo = scale - y + scale * tmp; in special_case()
130 lo = one - hi + y + lo; in special_case()
131 y = (hi + lo) - one; in special_case()
136 force_eval_double (opt_barrier_double (0x1p-1022) * 0x1p-1022); in special_case()
139 y = 0x1p-1022 * y; in special_case()
143 /* Computes sign*exp(x+xtail) where |xtail| < 2^-8/N and |xtail| <= |x|.
144 The sign_bias argument is SignBias or 0 and sets the sign to -1 or 1. */
145 static inline double
146 exp_inline (double x, double xtail, uint32_t sign_bias) in exp_inline()
149 if (unlikely (abstop - SmallExp >= ThresExp)) in exp_inline()
151 if (abstop - SmallExp >= 0x80000000) in exp_inline()
155 return sign_bias ? -1.0 : 1.0; in exp_inline()
165 double res_uoflow = asuint64 (x) >> 63 ? 0.0 : INFINITY; in exp_inline()
166 return sign_bias ? -res_uoflow : res_uoflow; in exp_inline()
173 /* exp(x) = 2^(k/N) * exp(r), with exp(r) in [2^(-1/2N),2^(1/2N)]. */ in exp_inline()
174 /* x = ln2/N*k + r, with int k and r in [-ln2/2N, ln2/2N]. */ in exp_inline()
175 double z = InvLn2N * x; in exp_inline()
176 double kd = round (z); in exp_inline()
178 double r = x - kd * Ln2HiN - kd * Ln2LoN; in exp_inline()
179 /* The code assumes 2^-200 < |xtail| < 2^-8/N. */ in exp_inline()
182 uint64_t idx = ki & (N_EXP - 1); in exp_inline()
183 uint64_t top = (ki + sign_bias) << (52 - V_POW_EXP_TABLE_BITS); in exp_inline()
184 /* This is only a valid scale when -1023*N < k < 1024*N. */ in exp_inline()
186 /* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (exp(r) - 1). */ in exp_inline()
188 double r2 = r * r; in exp_inline()
189 double tmp = r + r2 * Cs[0] + r * r2 * (Cs[1] + r * Cs[2]); in exp_inline()
192 double scale = asdouble (sbits); in exp_inline()
193 /* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there in exp_inline()
198 /* Computes exp(x+xtail) where |xtail| < 2^-8/N and |xtail| <= |x|.
201 static double NOINLINE
202 exp_nosignbias (double x, double xtail) in exp_nosignbias()
205 if (unlikely (abstop - SmallExp >= ThresExp)) in exp_nosignbias()
208 if (abstop - SmallExp >= 0x80000000) in exp_nosignbias()
221 /* exp(x) = 2^(k/N) * exp(r), with exp(r) in [2^(-1/2N),2^(1/2N)]. */ in exp_nosignbias()
222 /* x = ln2/N*k + r, with k integer and r in [-ln2/2N, ln2/2N]. */ in exp_nosignbias()
223 double z = InvLn2N * x; in exp_nosignbias()
224 double kd = round (z); in exp_nosignbias()
226 double r = x - kd * Ln2HiN - kd * Ln2LoN; in exp_nosignbias()
227 /* The code assumes 2^-200 < |xtail| < 2^-8/N. */ in exp_nosignbias()
230 uint64_t idx = ki & (N_EXP - 1); in exp_nosignbias()
231 uint64_t top = ki << (52 - V_POW_EXP_TABLE_BITS); in exp_nosignbias()
232 /* This is only a valid scale when -1023*N < k < 1024*N. */ in exp_nosignbias()
234 /* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (tail + exp(r) - 1). */ in exp_nosignbias()
235 double r2 = r * r; in exp_nosignbias()
236 double tmp = r + r2 * Cs[0] + r * r2 * (Cs[1] + r * Cs[2]); in exp_nosignbias()
239 double scale = asdouble (sbits); in exp_nosignbias()
240 /* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there in exp_nosignbias()
246 the bit representation of a non-zero finite floating-point value. */
255 if (iy & ((1ULL << (0x3ff + 52 - e)) - 1)) in checkint()
257 if (iy & (1ULL << (0x3ff + 52 - e))) in checkint()
266 return 2 * i - 1 >= 2 * asuint64 (INFINITY) - 1; in zeroinfnan()
269 static double NOINLINE
270 pow_scalar_special_case (double x, double y) in pow_scalar_special_case()
280 if (unlikely (topx - SmallPowX >= ThresPowX in pow_scalar_special_case()
281 || (topy & 0x7ff) - SmallPowY >= ThresPowY)) in pow_scalar_special_case()
284 and if |y| < 2^-54 / 1075 ~= 0x1.e7b6p-65 then pow(x,y) = +-1. */ in pow_scalar_special_case()
285 /* Special cases: (x < 0x1p-126 or inf or nan) or in pow_scalar_special_case()
286 (|y| < 0x1p-65 or |y| >= 0x1p63 or nan). */ in pow_scalar_special_case()
299 return 0.0; /* |x|<1 && y==inf or |x|>1 && y==-inf. */ in pow_scalar_special_case()
304 double x2 = x * x; in pow_scalar_special_case()
307 x2 = -x2; in pow_scalar_special_case()
316 /* Here x and y are non-zero finite. */ in pow_scalar_special_case()
332 if ((topy & 0x7ff) - SmallPowY >= ThresPowY) in pow_scalar_special_case()
337 /* |y| < 2^-65, x^y ~= 1 + y*log(x). */ in pow_scalar_special_case()
352 ix -= 52ULL << 52; in pow_scalar_special_case()
356 double lo; in pow_scalar_special_case()
357 double hi = log_inline (ix, &lo); in pow_scalar_special_case()
358 double ehi = y * hi; in pow_scalar_special_case()
359 double elo = y * lo + fma (y, hi, -ehi); in pow_scalar_special_case()