2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2016 Free Software Foundation, Inc.
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; either version 2.1 of the License, or
9 * (at your option) any later version.
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU Lesser General Public License for more details.
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with this program; if not, see <http://www.gnu.org/licenses/>.
19 /************************************************************************/
21 /* MODULE_NAME:halfulp.c */
23 /* FUNCTIONS:halfulp */
24 /* FILES NEEDED: mydefs.h dla.h endian.h */
27 /*Routine halfulp(double x, double y) computes x^y where result does */
28 /*not need rounding. If the result is closer to 0 than can be */
29 /*represented it returns 0. */
30 /* In the following cases the function does not compute anything */
31 /*and returns a negative number: */
32 /*1. if the result needs rounding, */
33 /*2. if y is outside the interval [0, 2^20-1], */
34 /*3. if x can be represented by x=2**n for some integer n. */
35 /************************************************************************/
40 #include <math_private.h>
46 static const int4 tab54
[32] = {
47 262143, 11585, 1782, 511, 210, 107, 63, 42,
48 30, 22, 17, 14, 12, 10, 9, 7,
49 7, 6, 5, 5, 5, 4, 4, 4,
50 3, 3, 3, 3, 3, 3, 3, 3
56 __halfulp (double x
, double y
)
61 double j1
, j2
, j3
, j4
, j5
;
64 if (y
<= 0) /*if power is negative or zero */
67 if (v
.i
[LOW_HALF
] != 0)
70 if (v
.i
[LOW_HALF
] != 0)
72 if ((v
.i
[HIGH_HALF
] & 0x000fffff) != 0)
73 return -10; /* if x =2 ^ n */
74 k
= ((v
.i
[HIGH_HALF
] & 0x7fffffff) >> 20) - 1023; /* find this n */
76 return (z
* y
== -1075.0) ? 0 : -10.0;
80 if (v
.i
[LOW_HALF
] != 0)
84 /* case where x = 2**n for some integer n */
85 if (((v
.i
[HIGH_HALF
] & 0x000fffff) | v
.i
[LOW_HALF
]) == 0)
87 k
= (v
.i
[HIGH_HALF
] >> 20) - 1023;
88 return (((double) k
) * y
== -1075.0) ? 0 : -10.0;
99 n
= (k
& 0x000fffff) | 0x00100000;
100 n
= n
>> (20 - l
); /* n is the odd integer of y */
101 k
= ((k
>> 20) - 1023) - l
; /* y = n*2**k */
116 z
= __ieee754_sqrt (x
);
117 EMULV (z
, z
, u
, uu
, j1
, j2
, j3
, j4
, j5
);
118 if (((u
- x
) + uu
) != 0)
126 /* it is impossible that n == 2, so the mantissa of x must be short */
138 m
= (k
& 0x000fffff) | 0x00100000;
139 m
= m
>> (20 - l
); /* m is the odd integer of x */
141 /* now check whether the length of m**n is at most 54 bits */
143 if (m
> tab54
[n
- 3])
146 /* yes, it is - now compute x**n by simple multiplications */
149 for (k
= 1; k
< n
; k
++)