]>
Commit | Line | Data |
---|---|---|
e4d82761 UD |
1 | /* |
2 | * IBM Accurate Mathematical Library | |
aeb25823 | 3 | * written by International Business Machines Corp. |
0ac5ae23 | 4 | * Copyright (C) 2001, 2005, 2011 Free Software Foundation |
e4d82761 UD |
5 | * |
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 | |
cc7375ce | 8 | * the Free Software Foundation; either version 2.1 of the License, or |
e4d82761 | 9 | * (at your option) any later version. |
ca58f1db | 10 | * |
e4d82761 UD |
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 | |
c6c6dd48 | 14 | * GNU Lesser General Public License for more details. |
e4d82761 UD |
15 | * |
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 | |
ca58f1db | 18 | * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. |
e4d82761 UD |
19 | */ |
20 | /************************************************************************/ | |
21 | /* */ | |
ca58f1db UD |
22 | /* MODULE_NAME:halfulp.c */ |
23 | /* */ | |
e4d82761 UD |
24 | /* FUNCTIONS:halfulp */ |
25 | /* FILES NEEDED: mydefs.h dla.h endian.h */ | |
26 | /* uroot.c */ | |
27 | /* */ | |
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 | /************************************************************************/ | |
37 | ||
ca58f1db | 38 | #include "endian.h" |
e4d82761 UD |
39 | #include "mydefs.h" |
40 | #include "dla.h" | |
15b3c029 | 41 | #include "math_private.h" |
e4d82761 | 42 | |
403a6325 | 43 | static const int4 tab54[32] = { |
e4d82761 UD |
44 | 262143, 11585, 1782, 511, 210, 107, 63, 42, |
45 | 30, 22, 17, 14, 12, 10, 9, 7, | |
0ac5ae23 UD |
46 | 7, 6, 5, 5, 5, 4, 4, 4, |
47 | 3, 3, 3, 3, 3, 3, 3, 3 }; | |
e4d82761 UD |
48 | |
49 | ||
ca58f1db | 50 | double __halfulp(double x, double y) |
e4d82761 UD |
51 | { |
52 | mynumber v; | |
53 | double z,u,uu,j1,j2,j3,j4,j5; | |
54 | int4 k,l,m,n; | |
55 | if (y <= 0) { /*if power is negative or zero */ | |
56 | v.x = y; | |
ca58f1db | 57 | if (v.i[LOW_HALF] != 0) return -10.0; |
e4d82761 | 58 | v.x = x; |
ca58f1db UD |
59 | if (v.i[LOW_HALF] != 0) return -10.0; |
60 | if ((v.i[HIGH_HALF]&0x000fffff) != 0) return -10; /* if x =2 ^ n */ | |
e4d82761 UD |
61 | k = ((v.i[HIGH_HALF]&0x7fffffff)>>20)-1023; /* find this n */ |
62 | z = (double) k; | |
63 | return (z*y == -1075.0)?0: -10.0; | |
64 | } | |
0ac5ae23 | 65 | /* if y > 0 */ |
e4d82761 UD |
66 | v.x = y; |
67 | if (v.i[LOW_HALF] != 0) return -10.0; | |
ca58f1db | 68 | |
e4d82761 | 69 | v.x=x; |
0ac5ae23 | 70 | /* case where x = 2**n for some integer n */ |
e4d82761 UD |
71 | if (((v.i[HIGH_HALF]&0x000fffff)|v.i[LOW_HALF]) == 0) { |
72 | k=(v.i[HIGH_HALF]>>20)-1023; | |
73 | return (((double) k)*y == -1075.0)?0:-10.0; | |
ca58f1db UD |
74 | } |
75 | ||
e4d82761 UD |
76 | v.x = y; |
77 | k = v.i[HIGH_HALF]; | |
78 | m = k<<12; | |
79 | l = 0; | |
ca58f1db | 80 | while (m) |
e4d82761 UD |
81 | {m = m<<1; l++; } |
82 | n = (k&0x000fffff)|0x00100000; | |
83 | n = n>>(20-l); /* n is the odd integer of y */ | |
84 | k = ((k>>20) -1023)-l; /* y = n*2**k */ | |
85 | if (k>5) return -10.0; | |
86 | if (k>0) for (;k>0;k--) n *= 2; | |
87 | if (n > 34) return -10.0; | |
88 | k = -k; | |
89 | if (k>5) return -10.0; | |
ca58f1db | 90 | |
0ac5ae23 | 91 | /* now treat x */ |
e4d82761 | 92 | while (k>0) { |
ca58f1db | 93 | z = __ieee754_sqrt(x); |
e4d82761 UD |
94 | EMULV(z,z,u,uu,j1,j2,j3,j4,j5); |
95 | if (((u-x)+uu) != 0) break; | |
96 | x = z; | |
97 | k--; | |
98 | } | |
ca58f1db UD |
99 | if (k) return -10.0; |
100 | ||
e4d82761 | 101 | /* it is impossible that n == 2, so the mantissa of x must be short */ |
ca58f1db | 102 | |
e4d82761 UD |
103 | v.x = x; |
104 | if (v.i[LOW_HALF]) return -10.0; | |
105 | k = v.i[HIGH_HALF]; | |
106 | m = k<<12; | |
107 | l = 0; | |
108 | while (m) {m = m<<1; l++; } | |
109 | m = (k&0x000fffff)|0x00100000; | |
110 | m = m>>(20-l); /* m is the odd integer of x */ | |
ca58f1db | 111 | |
0ac5ae23 | 112 | /* now check whether the length of m**n is at most 54 bits */ |
ca58f1db | 113 | |
e4d82761 | 114 | if (m > tab54[n-3]) return -10.0; |
ca58f1db | 115 | |
0ac5ae23 | 116 | /* yes, it is - now compute x**n by simple multiplications */ |
ca58f1db | 117 | |
e4d82761 UD |
118 | u = x; |
119 | for (k=1;k<n;k++) u = u*x; | |
120 | return u; | |
121 | } |