1/*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2017 Free Software Foundation, Inc.
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
8 * the Free Software Foundation; either version 2.1 of the License, or
9 * (at your option) any later version.
10 *
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.
15 *
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with this program; if not, see <http://www.gnu.org/licenses/>.
18 */
19/************************************************************************/
20/* */
21/* MODULE_NAME:halfulp.c */
22/* */
23/* FUNCTIONS:halfulp */
24/* FILES NEEDED: mydefs.h dla.h endian.h */
25/* uroot.c */
26/* */
27/*Routine halfulp(double x, double y) computes x^y where result does */
28/*not need rounding. If the result is closer to 0 than can be */
29/*represented it returns 0. */
30/* In the following cases the function does not compute anything */
31/*and returns a negative number: */
32/*1. if the result needs rounding, */
33/*2. if y is outside the interval [0, 2^20-1], */
34/*3. if x can be represented by x=2**n for some integer n. */
35/************************************************************************/
36
37#include "endian.h"
38#include "mydefs.h"
39#include <dla.h>
40#include <math_private.h>
41
42#ifndef SECTION
43# define SECTION
44#endif
45
46static const int4 tab54[32] = {
47 262143, 11585, 1782, 511, 210, 107, 63, 42,
48 30, 22, 17, 14, 12, 10, 9, 7,
49 7, 6, 5, 5, 5, 4, 4, 4,
50 3, 3, 3, 3, 3, 3, 3, 3
51};
52
53
54double
55SECTION
56__halfulp (double x, double y)
57{
58 mynumber v;
59 double z, u, uu;
60#ifndef DLA_FMS
61 double j1, j2, j3, j4, j5;
62#endif
63 int4 k, l, m, n;
64 if (y <= 0) /*if power is negative or zero */
65 {
66 v.x = y;
67 if (v.i[LOW_HALF] != 0)
68 return -10.0;
69 v.x = x;
70 if (v.i[LOW_HALF] != 0)
71 return -10.0;
72 if ((v.i[HIGH_HALF] & 0x000fffff) != 0)
73 return -10; /* if x =2 ^ n */
74 k = ((v.i[HIGH_HALF] & 0x7fffffff) >> 20) - 1023; /* find this n */
75 z = (double) k;
76 return (z * y == -1075.0) ? 0 : -10.0;
77 }
78 /* if y > 0 */
79 v.x = y;
80 if (v.i[LOW_HALF] != 0)
81 return -10.0;
82
83 v.x = x;
84 /* case where x = 2**n for some integer n */
85 if (((v.i[HIGH_HALF] & 0x000fffff) | v.i[LOW_HALF]) == 0)
86 {
87 k = (v.i[HIGH_HALF] >> 20) - 1023;
88 return (((double) k) * y == -1075.0) ? 0 : -10.0;
89 }
90
91 v.x = y;
92 k = v.i[HIGH_HALF];
93 m = k << 12;
94 l = 0;
95 while (m)
96 {
97 m = m << 1; l++;
98 }
99 n = (k & 0x000fffff) | 0x00100000;
100 n = n >> (20 - l); /* n is the odd integer of y */
101 k = ((k >> 20) - 1023) - l; /* y = n*2**k */
102 if (k > 5)
103 return -10.0;
104 if (k > 0)
105 for (; k > 0; k--)
106 n *= 2;
107 if (n > 34)
108 return -10.0;
109 k = -k;
110 if (k > 5)
111 return -10.0;
112
113 /* now treat x */
114 while (k > 0)
115 {
116 z = __ieee754_sqrt (x);
117 EMULV (z, z, u, uu, j1, j2, j3, j4, j5);
118 if (((u - x) + uu) != 0)
119 break;
120 x = z;
121 k--;
122 }
123 if (k)
124 return -10.0;
125
126 /* it is impossible that n == 2, so the mantissa of x must be short */
127
128 v.x = x;
129 if (v.i[LOW_HALF])
130 return -10.0;
131 k = v.i[HIGH_HALF];
132 m = k << 12;
133 l = 0;
134 while (m)
135 {
136 m = m << 1; l++;
137 }
138 m = (k & 0x000fffff) | 0x00100000;
139 m = m >> (20 - l); /* m is the odd integer of x */
140
141 /* now check whether the length of m**n is at most 54 bits */
142
143 if (m > tab54[n - 3])
144 return -10.0;
145
146 /* yes, it is - now compute x**n by simple multiplications */
147
148 u = x;
149 for (k = 1; k < n; k++)
150 u = u * x;
151 return u;
152}
153