mirror of
https://sourceware.org/git/glibc.git
synced 2024-12-12 22:30:12 +00:00
153 lines
4.5 KiB
C
153 lines
4.5 KiB
C
/*
|
|
* 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;
|
|
}
|