diff options
Diffstat (limited to 'sysdeps/ieee754/dbl-64/halfulp.c')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/halfulp.c | 52 |
1 files changed, 25 insertions, 27 deletions
diff --git a/sysdeps/ieee754/dbl-64/halfulp.c b/sysdeps/ieee754/dbl-64/halfulp.c index 7901ec4e36..929ca91c50 100644 --- a/sysdeps/ieee754/dbl-64/halfulp.c +++ b/sysdeps/ieee754/dbl-64/halfulp.c @@ -5,9 +5,9 @@ * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by - * the Free Software Foundation; either version 2 of the License, or + * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. - * + * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the @@ -15,12 +15,12 @@ * * You should have received a copy of the GNU Lesser General Public License * along with this program; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ /************************************************************************/ /* */ -/* MODULE_NAME:halfulp.c */ -/* */ +/* MODULE_NAME:halfulp.c */ +/* */ /* FUNCTIONS:halfulp */ /* FILES NEEDED: mydefs.h dla.h endian.h */ /* uroot.c */ @@ -35,11 +35,11 @@ /*3. if x can be represented by x=2**n for some integer n. */ /************************************************************************/ -#include "endian.h" +#include "endian.h" #include "mydefs.h" #include "dla.h" -double usqrt(double x); +double __ieee754_sqrt(double x); int4 tab54[32] = { 262143, 11585, 1782, 511, 210, 107, 63, 42, @@ -48,17 +48,17 @@ int4 tab54[32] = { 3, 3, 3, 3, 3, 3, 3, 3 }; -double halfulp(double x, double y) +double __halfulp(double x, double y) { mynumber v; double z,u,uu,j1,j2,j3,j4,j5; int4 k,l,m,n; if (y <= 0) { /*if power is negative or zero */ v.x = y; - if (v.i[LOW_HALF] != 0) return -10.0; + if (v.i[LOW_HALF] != 0) return -10.0; v.x = x; - if (v.i[LOW_HALF] != 0) return -10.0; - if ((v.i[HIGH_HALF]&0x000fffff) != 0) return -10; /* if x =2 ^ n */ + if (v.i[LOW_HALF] != 0) return -10.0; + if ((v.i[HIGH_HALF]&0x000fffff) != 0) return -10; /* if x =2 ^ n */ k = ((v.i[HIGH_HALF]&0x7fffffff)>>20)-1023; /* find this n */ z = (double) k; return (z*y == -1075.0)?0: -10.0; @@ -66,19 +66,19 @@ double halfulp(double x, double y) /* if y > 0 */ v.x = y; if (v.i[LOW_HALF] != 0) return -10.0; - + v.x=x; - /* case where x = 2**n for some integer n */ + /* case where x = 2**n for some integer n */ if (((v.i[HIGH_HALF]&0x000fffff)|v.i[LOW_HALF]) == 0) { k=(v.i[HIGH_HALF]>>20)-1023; return (((double) k)*y == -1075.0)?0:-10.0; - } - + } + v.x = y; k = v.i[HIGH_HALF]; m = k<<12; l = 0; - while (m) + while (m) {m = m<<1; l++; } n = (k&0x000fffff)|0x00100000; n = n>>(20-l); /* n is the odd integer of y */ @@ -88,19 +88,19 @@ double halfulp(double x, double y) if (n > 34) return -10.0; k = -k; if (k>5) return -10.0; - + /* now treat x */ while (k>0) { - z = usqrt(x); + z = __ieee754_sqrt(x); EMULV(z,z,u,uu,j1,j2,j3,j4,j5); if (((u-x)+uu) != 0) break; x = z; k--; } - if (k) return -10.0; - + if (k) return -10.0; + /* it is impossible that n == 2, so the mantissa of x must be short */ - + v.x = x; if (v.i[LOW_HALF]) return -10.0; k = v.i[HIGH_HALF]; @@ -109,16 +109,14 @@ double halfulp(double x, double y) while (m) {m = m<<1; l++; } m = (k&0x000fffff)|0x00100000; m = m>>(20-l); /* m is the odd integer of x */ - + /* now check whether the length of m**n is at most 54 bits */ - + if (m > tab54[n-3]) return -10.0; - + /* yes, it is - now compute x**n by simple multiplications */ - + u = x; for (k=1;k<n;k++) u = u*x; return u; } - - |