2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001, 2005, 2011 Free Software Foundation
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 };
55 __halfulp(double x
, double y
)
60 double j1
,j2
,j3
,j4
,j5
;
63 if (y
<= 0) { /*if power is negative or zero */
65 if (v
.i
[LOW_HALF
] != 0) return -10.0;
67 if (v
.i
[LOW_HALF
] != 0) return -10.0;
68 if ((v
.i
[HIGH_HALF
]&0x000fffff) != 0) return -10; /* if x =2 ^ n */
69 k
= ((v
.i
[HIGH_HALF
]&0x7fffffff)>>20)-1023; /* find this n */
71 return (z
*y
== -1075.0)?0: -10.0;
75 if (v
.i
[LOW_HALF
] != 0) return -10.0;
78 /* case where x = 2**n for some integer n */
79 if (((v
.i
[HIGH_HALF
]&0x000fffff)|v
.i
[LOW_HALF
]) == 0) {
80 k
=(v
.i
[HIGH_HALF
]>>20)-1023;
81 return (((double) k
)*y
== -1075.0)?0:-10.0;
90 n
= (k
&0x000fffff)|0x00100000;
91 n
= n
>>(20-l
); /* n is the odd integer of y */
92 k
= ((k
>>20) -1023)-l
; /* y = n*2**k */
93 if (k
>5) return -10.0;
94 if (k
>0) for (;k
>0;k
--) n
*= 2;
95 if (n
> 34) return -10.0;
97 if (k
>5) return -10.0;
101 z
= __ieee754_sqrt(x
);
102 EMULV(z
,z
,u
,uu
,j1
,j2
,j3
,j4
,j5
);
103 if (((u
-x
)+uu
) != 0) break;
109 /* it is impossible that n == 2, so the mantissa of x must be short */
112 if (v
.i
[LOW_HALF
]) return -10.0;
116 while (m
) {m
= m
<<1; l
++; }
117 m
= (k
&0x000fffff)|0x00100000;
118 m
= m
>>(20-l
); /* m is the odd integer of x */
120 /* now check whether the length of m**n is at most 54 bits */
122 if (m
> tab54
[n
-3]) return -10.0;
124 /* yes, it is - now compute x**n by simple multiplications */
127 for (k
=1;k
<n
;k
++) u
= u
*x
;