2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001, 2005 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, write to the Free Software
18 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
20 /************************************************************************/
22 /* MODULE_NAME:halfulp.c */
24 /* FUNCTIONS:halfulp */
25 /* FILES NEEDED: mydefs.h dla.h endian.h */
28 /*Routine halfulp(double x, double y) computes x^y where result does */
29 /*not need rounding. If the result is closer to 0 than can be */
30 /*represented it returns 0. */
31 /* In the following cases the function does not compute anything */
32 /*and returns a negative number: */
33 /*1. if the result needs rounding, */
34 /*2. if y is outside the interval [0, 2^20-1], */
35 /*3. if x can be represented by x=2**n for some integer n. */
36 /************************************************************************/
41 #include "math_private.h"
43 double __ieee754_sqrt(double x
);
45 static const int4 tab54
[32] = {
46 262143, 11585, 1782, 511, 210, 107, 63, 42,
47 30, 22, 17, 14, 12, 10, 9, 7,
48 7, 6, 5, 5, 5, 4, 4, 4,
49 3, 3, 3, 3, 3, 3, 3, 3 };
52 double __halfulp(double x
, double y
)
55 double z
,u
,uu
,j1
,j2
,j3
,j4
,j5
;
57 if (y
<= 0) { /*if power is negative or zero */
59 if (v
.i
[LOW_HALF
] != 0) return -10.0;
61 if (v
.i
[LOW_HALF
] != 0) return -10.0;
62 if ((v
.i
[HIGH_HALF
]&0x000fffff) != 0) return -10; /* if x =2 ^ n */
63 k
= ((v
.i
[HIGH_HALF
]&0x7fffffff)>>20)-1023; /* find this n */
65 return (z
*y
== -1075.0)?0: -10.0;
69 if (v
.i
[LOW_HALF
] != 0) return -10.0;
72 /* case where x = 2**n for some integer n */
73 if (((v
.i
[HIGH_HALF
]&0x000fffff)|v
.i
[LOW_HALF
]) == 0) {
74 k
=(v
.i
[HIGH_HALF
]>>20)-1023;
75 return (((double) k
)*y
== -1075.0)?0:-10.0;
84 n
= (k
&0x000fffff)|0x00100000;
85 n
= n
>>(20-l
); /* n is the odd integer of y */
86 k
= ((k
>>20) -1023)-l
; /* y = n*2**k */
87 if (k
>5) return -10.0;
88 if (k
>0) for (;k
>0;k
--) n
*= 2;
89 if (n
> 34) return -10.0;
91 if (k
>5) return -10.0;
95 z
= __ieee754_sqrt(x
);
96 EMULV(z
,z
,u
,uu
,j1
,j2
,j3
,j4
,j5
);
97 if (((u
-x
)+uu
) != 0) break;
103 /* it is impossible that n == 2, so the mantissa of x must be short */
106 if (v
.i
[LOW_HALF
]) return -10.0;
110 while (m
) {m
= m
<<1; l
++; }
111 m
= (k
&0x000fffff)|0x00100000;
112 m
= m
>>(20-l
); /* m is the odd integer of x */
114 /* now check whether the length of m**n is at most 54 bits */
116 if (m
> tab54
[n
-3]) return -10.0;
118 /* yes, it is - now compute x**n by simple multiplications */
121 for (k
=1;k
<n
;k
++) u
= u
*x
;