1 /*
2 * Double-precision SVE pow(x, y) function.
3 *
4 * Copyright (c) 2022-2023, Arm Limited.
5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6 */
7
8 #include "sv_math.h"
9 #include "pl_sig.h"
10 #include "pl_test.h"
11
12 /* This version share a similar algorithm as AOR scalar pow.
13
14 The core computation consists in computing pow(x, y) as
15
16 exp (y * log (x)).
17
18 The algorithms for exp and log are very similar to scalar exp and log.
19 The log relies on table lookup for 3 variables and an order 8 polynomial.
20 It returns a high and a low contribution that are then passed to the exp,
21 to minimise the loss of accuracy in both routines.
22 The exp is based on 8-bit table lookup for scale and order-4 polynomial.
23 The SVE algorithm drops the tail in the exp computation at the price of
24 a lower accuracy, slightly above 1ULP.
25 The SVE algorithm also drops the special treatement of small (< 2^-65) and
26 large (> 2^63) finite values of |y|, as they only affect non-round to nearest
27 modes.
28
29 Maximum measured error is 1.04 ULPs:
30 SV_NAME_D2 (pow) (0x1.3d2d45bc848acp+63, -0x1.a48a38b40cd43p-12)
31 got 0x1.f7116284221fcp-1
32 want 0x1.f7116284221fdp-1. */
33
34 /* Data is defined in v_pow_log_data.c. */
35 #define N_LOG (1 << V_POW_LOG_TABLE_BITS)
36 #define A __v_pow_log_data.poly
37 #define Off 0x3fe6955500000000
38
39 /* Data is defined in v_pow_exp_data.c. */
40 #define N_EXP (1 << V_POW_EXP_TABLE_BITS)
41 #define SignBias (0x800 << V_POW_EXP_TABLE_BITS)
42 #define C __v_pow_exp_data.poly
43 #define SmallExp 0x3c9 /* top12(0x1p-54). */
44 #define BigExp 0x408 /* top12(512.). */
45 #define ThresExp 0x03f /* BigExp - SmallExp. */
46 #define HugeExp 0x409 /* top12(1024.). */
47
48 /* Constants associated with pow. */
49 #define SmallPowX 0x001 /* top12(0x1p-126). */
50 #define BigPowX 0x7ff /* top12(INFINITY). */
51 #define ThresPowX 0x7fe /* BigPowX - SmallPowX. */
52 #define SmallPowY 0x3be /* top12(0x1.e7b6p-65). */
53 #define BigPowY 0x43e /* top12(0x1.749p62). */
54 #define ThresPowY 0x080 /* BigPowY - SmallPowY. */
55
56 /* Check if x is an integer. */
57 static inline svbool_t
sv_isint(svbool_t pg,svfloat64_t x)58 sv_isint (svbool_t pg, svfloat64_t x)
59 {
60 return svcmpeq (pg, svrintz_z (pg, x), x);
61 }
62
63 /* Check if x is real not integer valued. */
64 static inline svbool_t
sv_isnotint(svbool_t pg,svfloat64_t x)65 sv_isnotint (svbool_t pg, svfloat64_t x)
66 {
67 return svcmpne (pg, svrintz_z (pg, x), x);
68 }
69
70 /* Check if x is an odd integer. */
71 static inline svbool_t
sv_isodd(svbool_t pg,svfloat64_t x)72 sv_isodd (svbool_t pg, svfloat64_t x)
73 {
74 svfloat64_t y = svmul_x (pg, x, 0.5);
75 return sv_isnotint (pg, y);
76 }
77
78 /* Returns 0 if not int, 1 if odd int, 2 if even int. The argument is
79 the bit representation of a non-zero finite floating-point value. */
80 static inline int
checkint(uint64_t iy)81 checkint (uint64_t iy)
82 {
83 int e = iy >> 52 & 0x7ff;
84 if (e < 0x3ff)
85 return 0;
86 if (e > 0x3ff + 52)
87 return 2;
88 if (iy & ((1ULL << (0x3ff + 52 - e)) - 1))
89 return 0;
90 if (iy & (1ULL << (0x3ff + 52 - e)))
91 return 1;
92 return 2;
93 }
94
95 /* Top 12 bits (sign and exponent of each double float lane). */
96 static inline svuint64_t
sv_top12(svfloat64_t x)97 sv_top12 (svfloat64_t x)
98 {
99 return svlsr_x (svptrue_b64 (), svreinterpret_u64 (x), 52);
100 }
101
102 /* Returns 1 if input is the bit representation of 0, infinity or nan. */
103 static inline int
zeroinfnan(uint64_t i)104 zeroinfnan (uint64_t i)
105 {
106 return 2 * i - 1 >= 2 * asuint64 (INFINITY) - 1;
107 }
108
109 /* Returns 1 if input is the bit representation of 0, infinity or nan. */
110 static inline svbool_t
sv_zeroinfnan(svbool_t pg,svuint64_t i)111 sv_zeroinfnan (svbool_t pg, svuint64_t i)
112 {
113 return svcmpge (pg, svsub_x (pg, svmul_x (pg, i, 2), 1),
114 2 * asuint64 (INFINITY) - 1);
115 }
116
117 /* Handle cases that may overflow or underflow when computing the result that
118 is scale*(1+TMP) without intermediate rounding. The bit representation of
119 scale is in SBITS, however it has a computed exponent that may have
120 overflown into the sign bit so that needs to be adjusted before using it as
121 a double. (int32_t)KI is the k used in the argument reduction and exponent
122 adjustment of scale, positive k here means the result may overflow and
123 negative k means the result may underflow. */
124 static inline double
specialcase(double tmp,uint64_t sbits,uint64_t ki)125 specialcase (double tmp, uint64_t sbits, uint64_t ki)
126 {
127 double scale;
128 if ((ki & 0x80000000) == 0)
129 {
130 /* k > 0, the exponent of scale might have overflowed by <= 460. */
131 sbits -= 1009ull << 52;
132 scale = asdouble (sbits);
133 return 0x1p1009 * (scale + scale * tmp);
134 }
135 /* k < 0, need special care in the subnormal range. */
136 sbits += 1022ull << 52;
137 /* Note: sbits is signed scale. */
138 scale = asdouble (sbits);
139 double y = scale + scale * tmp;
140 return 0x1p-1022 * y;
141 }
142
143 /* Scalar fallback for special cases of SVE pow's exp. */
144 static inline svfloat64_t
sv_call_specialcase(svfloat64_t x1,svuint64_t u1,svuint64_t u2,svfloat64_t y,svbool_t cmp)145 sv_call_specialcase (svfloat64_t x1, svuint64_t u1, svuint64_t u2,
146 svfloat64_t y, svbool_t cmp)
147 {
148 svbool_t p = svpfirst (cmp, svpfalse ());
149 while (svptest_any (cmp, p))
150 {
151 double sx1 = svclastb (p, 0, x1);
152 uint64_t su1 = svclastb (p, 0, u1);
153 uint64_t su2 = svclastb (p, 0, u2);
154 double elem = specialcase (sx1, su1, su2);
155 svfloat64_t y2 = sv_f64 (elem);
156 y = svsel (p, y2, y);
157 p = svpnext_b64 (cmp, p);
158 }
159 return y;
160 }
161
162 /* Compute y+TAIL = log(x) where the rounded result is y and TAIL has about
163 additional 15 bits precision. IX is the bit representation of x, but
164 normalized in the subnormal range using the sign bit for the exponent. */
165 static inline svfloat64_t
sv_log_inline(svbool_t pg,svuint64_t ix,svfloat64_t * tail)166 sv_log_inline (svbool_t pg, svuint64_t ix, svfloat64_t *tail)
167 {
168 /* x = 2^k z; where z is in range [Off,2*Off) and exact.
169 The range is split into N subintervals.
170 The ith subinterval contains z and c is near its center. */
171 svuint64_t tmp = svsub_x (pg, ix, Off);
172 svuint64_t i = svand_x (pg, svlsr_x (pg, tmp, 52 - V_POW_LOG_TABLE_BITS),
173 sv_u64 (N_LOG - 1));
174 svint64_t k = svasr_x (pg, svreinterpret_s64 (tmp), 52);
175 svuint64_t iz = svsub_x (pg, ix, svand_x (pg, tmp, sv_u64 (0xfffULL << 52)));
176 svfloat64_t z = svreinterpret_f64 (iz);
177 svfloat64_t kd = svcvt_f64_x (pg, k);
178
179 /* log(x) = k*Ln2 + log(c) + log1p(z/c-1). */
180 /* SVE lookup requires 3 separate lookup tables, as opposed to scalar version
181 that uses array of structures. We also do the lookup earlier in the code to
182 make sure it finishes as early as possible. */
183 svfloat64_t invc = svld1_gather_index (pg, __v_pow_log_data.invc, i);
184 svfloat64_t logc = svld1_gather_index (pg, __v_pow_log_data.logc, i);
185 svfloat64_t logctail = svld1_gather_index (pg, __v_pow_log_data.logctail, i);
186
187 /* Note: 1/c is j/N or j/N/2 where j is an integer in [N,2N) and
188 |z/c - 1| < 1/N, so r = z/c - 1 is exactly representible. */
189 svfloat64_t r = svmad_x (pg, z, invc, -1.0);
190 /* k*Ln2 + log(c) + r. */
191 svfloat64_t t1 = svmla_x (pg, logc, kd, __v_pow_log_data.ln2_hi);
192 svfloat64_t t2 = svadd_x (pg, t1, r);
193 svfloat64_t lo1 = svmla_x (pg, logctail, kd, __v_pow_log_data.ln2_lo);
194 svfloat64_t lo2 = svadd_x (pg, svsub_x (pg, t1, t2), r);
195
196 /* Evaluation is optimized assuming superscalar pipelined execution. */
197 svfloat64_t ar = svmul_x (pg, r, -0.5); /* A[0] = -0.5. */
198 svfloat64_t ar2 = svmul_x (pg, r, ar);
199 svfloat64_t ar3 = svmul_x (pg, r, ar2);
200 /* k*Ln2 + log(c) + r + A[0]*r*r. */
201 svfloat64_t hi = svadd_x (pg, t2, ar2);
202 svfloat64_t lo3 = svmla_x (pg, svneg_x (pg, ar2), ar, r);
203 svfloat64_t lo4 = svadd_x (pg, svsub_x (pg, t2, hi), ar2);
204 /* p = log1p(r) - r - A[0]*r*r. */
205 /* p = (ar3 * (A[1] + r * A[2] + ar2 * (A[3] + r * A[4] + ar2 * (A[5] + r *
206 A[6])))). */
207 svfloat64_t a56 = svmla_x (pg, sv_f64 (A[5]), r, A[6]);
208 svfloat64_t a34 = svmla_x (pg, sv_f64 (A[3]), r, A[4]);
209 svfloat64_t a12 = svmla_x (pg, sv_f64 (A[1]), r, A[2]);
210 svfloat64_t p = svmla_x (pg, a34, ar2, a56);
211 p = svmla_x (pg, a12, ar2, p);
212 p = svmul_x (pg, ar3, p);
213 svfloat64_t lo = svadd_x (
214 pg, svadd_x (pg, svadd_x (pg, svadd_x (pg, lo1, lo2), lo3), lo4), p);
215 svfloat64_t y = svadd_x (pg, hi, lo);
216 *tail = svadd_x (pg, svsub_x (pg, hi, y), lo);
217 return y;
218 }
219
220 /* Computes sign*exp(x+xtail) where |xtail| < 2^-8/N and |xtail| <= |x|.
221 The sign_bias argument is SignBias or 0 and sets the sign to -1 or 1. */
222 static inline svfloat64_t
sv_exp_inline(svbool_t pg,svfloat64_t x,svfloat64_t xtail,svuint64_t sign_bias)223 sv_exp_inline (svbool_t pg, svfloat64_t x, svfloat64_t xtail,
224 svuint64_t sign_bias)
225 {
226 /* 3 types of special cases: tiny (uflow and spurious uflow), huge (oflow)
227 and other cases of large values of x (scale * (1 + TMP) oflow). */
228 svuint64_t abstop = svand_x (pg, sv_top12 (x), 0x7ff);
229 /* |x| is large (|x| >= 512) or tiny (|x| <= 0x1p-54). */
230 svbool_t uoflow = svcmpge (pg, svsub_x (pg, abstop, SmallExp), ThresExp);
231
232 /* Conditions special, uflow and oflow are all expressed as uoflow &&
233 something, hence do not bother computing anything if no lane in uoflow is
234 true. */
235 svbool_t special = svpfalse_b ();
236 svbool_t uflow = svpfalse_b ();
237 svbool_t oflow = svpfalse_b ();
238 if (unlikely (svptest_any (pg, uoflow)))
239 {
240 /* |x| is tiny (|x| <= 0x1p-54). */
241 uflow = svcmpge (pg, svsub_x (pg, abstop, SmallExp), 0x80000000);
242 uflow = svand_z (pg, uoflow, uflow);
243 /* |x| is huge (|x| >= 1024). */
244 oflow = svcmpge (pg, abstop, HugeExp);
245 oflow = svand_z (pg, uoflow, svbic_z (pg, oflow, uflow));
246 /* For large |x| values (512 < |x| < 1024) scale * (1 + TMP) can overflow
247 or underflow. */
248 special = svbic_z (pg, uoflow, svorr_z (pg, uflow, oflow));
249 }
250
251 /* exp(x) = 2^(k/N) * exp(r), with exp(r) in [2^(-1/2N),2^(1/2N)]. */
252 /* x = ln2/N*k + r, with int k and r in [-ln2/2N, ln2/2N]. */
253 svfloat64_t z = svmul_x (pg, x, __v_pow_exp_data.n_over_ln2);
254 /* z - kd is in [-1, 1] in non-nearest rounding modes. */
255 svfloat64_t shift = sv_f64 (__v_pow_exp_data.shift);
256 svfloat64_t kd = svadd_x (pg, z, shift);
257 svuint64_t ki = svreinterpret_u64 (kd);
258 kd = svsub_x (pg, kd, shift);
259 svfloat64_t r = x;
260 r = svmls_x (pg, r, kd, __v_pow_exp_data.ln2_over_n_hi);
261 r = svmls_x (pg, r, kd, __v_pow_exp_data.ln2_over_n_lo);
262 /* The code assumes 2^-200 < |xtail| < 2^-8/N. */
263 r = svadd_x (pg, r, xtail);
264 /* 2^(k/N) ~= scale. */
265 svuint64_t idx = svand_x (pg, ki, N_EXP - 1);
266 svuint64_t top
267 = svlsl_x (pg, svadd_x (pg, ki, sign_bias), 52 - V_POW_EXP_TABLE_BITS);
268 /* This is only a valid scale when -1023*N < k < 1024*N. */
269 svuint64_t sbits = svld1_gather_index (pg, __v_pow_exp_data.sbits, idx);
270 sbits = svadd_x (pg, sbits, top);
271 /* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (exp(r) - 1). */
272 svfloat64_t r2 = svmul_x (pg, r, r);
273 svfloat64_t tmp = svmla_x (pg, sv_f64 (C[1]), r, C[2]);
274 tmp = svmla_x (pg, sv_f64 (C[0]), r, tmp);
275 tmp = svmla_x (pg, r, r2, tmp);
276 svfloat64_t scale = svreinterpret_f64 (sbits);
277 /* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there
278 is no spurious underflow here even without fma. */
279 z = svmla_x (pg, scale, scale, tmp);
280
281 /* Update result with special and large cases. */
282 if (unlikely (svptest_any (pg, special)))
283 z = sv_call_specialcase (tmp, sbits, ki, z, special);
284
285 /* Handle underflow and overflow. */
286 svuint64_t sign_bit = svlsr_x (pg, svreinterpret_u64 (x), 63);
287 svbool_t x_is_neg = svcmpne (pg, sign_bit, 0);
288 svuint64_t sign_mask = svlsl_x (pg, sign_bias, 52 - V_POW_EXP_TABLE_BITS);
289 svfloat64_t res_uoflow = svsel (x_is_neg, sv_f64 (0.0), sv_f64 (INFINITY));
290 res_uoflow = svreinterpret_f64 (
291 svorr_x (pg, svreinterpret_u64 (res_uoflow), sign_mask));
292 z = svsel (oflow, res_uoflow, z);
293 /* Avoid spurious underflow for tiny x. */
294 svfloat64_t res_spurious_uflow
295 = svreinterpret_f64 (svorr_x (pg, sign_mask, 0x3ff0000000000000));
296 z = svsel (uflow, res_spurious_uflow, z);
297
298 return z;
299 }
300
301 static inline double
pow_sc(double x,double y)302 pow_sc (double x, double y)
303 {
304 uint64_t ix = asuint64 (x);
305 uint64_t iy = asuint64 (y);
306 /* Special cases: |x| or |y| is 0, inf or nan. */
307 if (unlikely (zeroinfnan (iy)))
308 {
309 if (2 * iy == 0)
310 return issignaling_inline (x) ? x + y : 1.0;
311 if (ix == asuint64 (1.0))
312 return issignaling_inline (y) ? x + y : 1.0;
313 if (2 * ix > 2 * asuint64 (INFINITY) || 2 * iy > 2 * asuint64 (INFINITY))
314 return x + y;
315 if (2 * ix == 2 * asuint64 (1.0))
316 return 1.0;
317 if ((2 * ix < 2 * asuint64 (1.0)) == !(iy >> 63))
318 return 0.0; /* |x|<1 && y==inf or |x|>1 && y==-inf. */
319 return y * y;
320 }
321 if (unlikely (zeroinfnan (ix)))
322 {
323 double_t x2 = x * x;
324 if (ix >> 63 && checkint (iy) == 1)
325 x2 = -x2;
326 /* Without the barrier some versions of clang hoist the 1/x2 and
327 thus division by zero exception can be signaled spuriously. */
328 return (iy >> 63) ? opt_barrier_double (1 / x2) : x2;
329 }
330 return x;
331 }
332
SV_NAME_D2(pow)333 svfloat64_t SV_NAME_D2 (pow) (svfloat64_t x, svfloat64_t y, const svbool_t pg)
334 {
335 /* This preamble handles special case conditions used in the final scalar
336 fallbacks. It also updates ix and sign_bias, that are used in the core
337 computation too, i.e., exp( y * log (x) ). */
338 svuint64_t vix0 = svreinterpret_u64 (x);
339 svuint64_t viy0 = svreinterpret_u64 (y);
340 svuint64_t vtopx0 = svlsr_x (svptrue_b64 (), vix0, 52);
341
342 /* Negative x cases. */
343 svuint64_t sign_bit = svlsr_m (pg, vix0, 63);
344 svbool_t xisneg = svcmpeq (pg, sign_bit, 1);
345
346 /* Set sign_bias and ix depending on sign of x and nature of y. */
347 svbool_t yisnotint_xisneg = svpfalse_b ();
348 svuint64_t sign_bias = sv_u64 (0);
349 svuint64_t vix = vix0;
350 svuint64_t vtopx1 = vtopx0;
351 if (unlikely (svptest_any (pg, xisneg)))
352 {
353 /* Determine nature of y. */
354 yisnotint_xisneg = sv_isnotint (xisneg, y);
355 svbool_t yisint_xisneg = sv_isint (xisneg, y);
356 svbool_t yisodd_xisneg = sv_isodd (xisneg, y);
357 /* ix set to abs(ix) if y is integer. */
358 vix = svand_m (yisint_xisneg, vix0, 0x7fffffffffffffff);
359 vtopx1 = svand_m (yisint_xisneg, vtopx0, 0x7ff);
360 /* Set to SignBias if x is negative and y is odd. */
361 sign_bias = svsel (yisodd_xisneg, sv_u64 (SignBias), sv_u64 (0));
362 }
363
364 /* Special cases of x or y: zero, inf and nan. */
365 svbool_t xspecial = sv_zeroinfnan (pg, vix0);
366 svbool_t yspecial = sv_zeroinfnan (pg, viy0);
367 svbool_t special = svorr_z (pg, xspecial, yspecial);
368
369 /* Small cases of x: |x| < 0x1p-126. */
370 svuint64_t vabstopx0 = svand_x (pg, vtopx0, 0x7ff);
371 svbool_t xsmall = svcmplt (pg, vabstopx0, SmallPowX);
372 if (unlikely (svptest_any (pg, xsmall)))
373 {
374 /* Normalize subnormal x so exponent becomes negative. */
375 svbool_t topx_is_null = svcmpeq (xsmall, vtopx1, 0);
376
377 svuint64_t vix_norm = svreinterpret_u64 (svmul_m (xsmall, x, 0x1p52));
378 vix_norm = svand_m (xsmall, vix_norm, 0x7fffffffffffffff);
379 vix_norm = svsub_m (xsmall, vix_norm, 52ULL << 52);
380 vix = svsel (topx_is_null, vix_norm, vix);
381 }
382
383 /* y_hi = log(ix, &y_lo). */
384 svfloat64_t vlo;
385 svfloat64_t vhi = sv_log_inline (pg, vix, &vlo);
386
387 /* z = exp(y_hi, y_lo, sign_bias). */
388 svfloat64_t vehi = svmul_x (pg, y, vhi);
389 svfloat64_t velo = svmul_x (pg, y, vlo);
390 svfloat64_t vemi = svmls_x (pg, vehi, y, vhi);
391 velo = svsub_x (pg, velo, vemi);
392 svfloat64_t vz = sv_exp_inline (pg, vehi, velo, sign_bias);
393
394 /* Cases of finite y and finite negative x. */
395 vz = svsel (yisnotint_xisneg, sv_f64 (__builtin_nan ("")), vz);
396
397 /* Cases of zero/inf/nan x or y. */
398 if (unlikely (svptest_any (pg, special)))
399 vz = sv_call2_f64 (pow_sc, x, y, vz, special);
400
401 return vz;
402 }
403
404 PL_SIG (SV, D, 2, pow)
405 PL_TEST_ULP (SV_NAME_D2 (pow), 0.55)
406 /* Wide intervals spanning the whole domain but shared between x and y. */
407 #define SV_POW_INTERVAL2(xlo, xhi, ylo, yhi, n) \
408 PL_TEST_INTERVAL2 (SV_NAME_D2 (pow), xlo, xhi, ylo, yhi, n) \
409 PL_TEST_INTERVAL2 (SV_NAME_D2 (pow), xlo, xhi, -ylo, -yhi, n) \
410 PL_TEST_INTERVAL2 (SV_NAME_D2 (pow), -xlo, -xhi, ylo, yhi, n) \
411 PL_TEST_INTERVAL2 (SV_NAME_D2 (pow), -xlo, -xhi, -ylo, -yhi, n)
412 #define EXPAND(str) str##000000000
413 #define SHL52(str) EXPAND (str)
414 SV_POW_INTERVAL2 (0, SHL52 (SmallPowX), 0, inf, 40000)
415 SV_POW_INTERVAL2 (SHL52 (SmallPowX), SHL52 (BigPowX), 0, inf, 40000)
416 SV_POW_INTERVAL2 (SHL52 (BigPowX), inf, 0, inf, 40000)
417 SV_POW_INTERVAL2 (0, inf, 0, SHL52 (SmallPowY), 40000)
418 SV_POW_INTERVAL2 (0, inf, SHL52 (SmallPowY), SHL52 (BigPowY), 40000)
419 SV_POW_INTERVAL2 (0, inf, SHL52 (BigPowY), inf, 40000)
420 SV_POW_INTERVAL2 (0, inf, 0, inf, 1000)
421 /* x~1 or y~1. */
422 SV_POW_INTERVAL2 (0x1p-1, 0x1p1, 0x1p-10, 0x1p10, 10000)
423 SV_POW_INTERVAL2 (0x1.ep-1, 0x1.1p0, 0x1p8, 0x1p16, 10000)
424 SV_POW_INTERVAL2 (0x1p-500, 0x1p500, 0x1p-1, 0x1p1, 10000)
425 /* around estimated argmaxs of ULP error. */
426 SV_POW_INTERVAL2 (0x1p-300, 0x1p-200, 0x1p-20, 0x1p-10, 10000)
427 SV_POW_INTERVAL2 (0x1p50, 0x1p100, 0x1p-20, 0x1p-10, 10000)
428 /* x is negative, y is odd or even integer, or y is real not integer. */
429 PL_TEST_INTERVAL2 (SV_NAME_D2 (pow), -0.0, -10.0, 3.0, 3.0, 10000)
430 PL_TEST_INTERVAL2 (SV_NAME_D2 (pow), -0.0, -10.0, 4.0, 4.0, 10000)
431 PL_TEST_INTERVAL2 (SV_NAME_D2 (pow), -0.0, -10.0, 0.0, 10.0, 10000)
432 PL_TEST_INTERVAL2 (SV_NAME_D2 (pow), 0.0, 10.0, -0.0, -10.0, 10000)
433 /* |x| is inf, y is odd or even integer, or y is real not integer. */
434 SV_POW_INTERVAL2 (inf, inf, 0.5, 0.5, 1)
435 SV_POW_INTERVAL2 (inf, inf, 1.0, 1.0, 1)
436 SV_POW_INTERVAL2 (inf, inf, 2.0, 2.0, 1)
437 SV_POW_INTERVAL2 (inf, inf, 3.0, 3.0, 1)
438 /* 0.0^y. */
439 SV_POW_INTERVAL2 (0.0, 0.0, 0.0, 0x1p120, 1000)
440 /* 1.0^y. */
441 PL_TEST_INTERVAL2 (SV_NAME_D2 (pow), 1.0, 1.0, 0.0, 0x1p-50, 1000)
442 PL_TEST_INTERVAL2 (SV_NAME_D2 (pow), 1.0, 1.0, 0x1p-50, 1.0, 1000)
443 PL_TEST_INTERVAL2 (SV_NAME_D2 (pow), 1.0, 1.0, 1.0, 0x1p100, 1000)
444 PL_TEST_INTERVAL2 (SV_NAME_D2 (pow), 1.0, 1.0, -1.0, -0x1p120, 1000)
445