diff options
Diffstat (limited to 'REORG.TODO/sysdeps/ieee754/dbl-64/halfulp.c')
| -rw-r--r-- | REORG.TODO/sysdeps/ieee754/dbl-64/halfulp.c | 152 |
1 files changed, 152 insertions, 0 deletions
diff --git a/REORG.TODO/sysdeps/ieee754/dbl-64/halfulp.c b/REORG.TODO/sysdeps/ieee754/dbl-64/halfulp.c new file mode 100644 index 0000000000..d5f8a010e2 --- /dev/null +++ b/REORG.TODO/sysdeps/ieee754/dbl-64/halfulp.c @@ -0,0 +1,152 @@ +/* + * IBM Accurate Mathematical Library + * written by International Business Machines Corp. + * Copyright (C) 2001-2017 Free Software Foundation, Inc. + * + * 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.1 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 + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program; if not, see <http://www.gnu.org/licenses/>. + */ +/************************************************************************/ +/* */ +/* MODULE_NAME:halfulp.c */ +/* */ +/* FUNCTIONS:halfulp */ +/* FILES NEEDED: mydefs.h dla.h endian.h */ +/* uroot.c */ +/* */ +/*Routine halfulp(double x, double y) computes x^y where result does */ +/*not need rounding. If the result is closer to 0 than can be */ +/*represented it returns 0. */ +/* In the following cases the function does not compute anything */ +/*and returns a negative number: */ +/*1. if the result needs rounding, */ +/*2. if y is outside the interval [0, 2^20-1], */ +/*3. if x can be represented by x=2**n for some integer n. */ +/************************************************************************/ + +#include "endian.h" +#include "mydefs.h" +#include <dla.h> +#include <math_private.h> + +#ifndef SECTION +# define SECTION +#endif + +static const int4 tab54[32] = { + 262143, 11585, 1782, 511, 210, 107, 63, 42, + 30, 22, 17, 14, 12, 10, 9, 7, + 7, 6, 5, 5, 5, 4, 4, 4, + 3, 3, 3, 3, 3, 3, 3, 3 +}; + + +double +SECTION +__halfulp (double x, double y) +{ + mynumber v; + double z, u, uu; +#ifndef DLA_FMS + double j1, j2, j3, j4, j5; +#endif + 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; + 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 */ + k = ((v.i[HIGH_HALF] & 0x7fffffff) >> 20) - 1023; /* find this n */ + z = (double) k; + return (z * y == -1075.0) ? 0 : -10.0; + } + /* 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 */ + 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) + { + m = m << 1; l++; + } + n = (k & 0x000fffff) | 0x00100000; + n = n >> (20 - l); /* n is the odd integer of y */ + k = ((k >> 20) - 1023) - l; /* y = n*2**k */ + if (k > 5) + return -10.0; + if (k > 0) + for (; k > 0; k--) + n *= 2; + if (n > 34) + return -10.0; + k = -k; + if (k > 5) + return -10.0; + + /* now treat x */ + while (k > 0) + { + 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; + + /* 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]; + m = k << 12; + l = 0; + 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; +} |