xref: /illumos-gate/usr/src/lib/libm/common/C/pow.c (revision ddc0e0b53c661f6e439e3b7072b3ef353eadb4af)
1 /*
2  * CDDL HEADER START
3  *
4  * The contents of this file are subject to the terms of the
5  * Common Development and Distribution License (the "License").
6  * You may not use this file except in compliance with the License.
7  *
8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9  * or http://www.opensolaris.org/os/licensing.
10  * See the License for the specific language governing permissions
11  * and limitations under the License.
12  *
13  * When distributing Covered Code, include this CDDL HEADER in each
14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15  * If applicable, add the following below this CDDL HEADER, with the
16  * fields enclosed by brackets "[]" replaced with your own identifying
17  * information: Portions Copyright [yyyy] [name of copyright owner]
18  *
19  * CDDL HEADER END
20  */
21 
22 /*
23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
24  */
25 /*
26  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
27  * Use is subject to license terms.
28  */
29 
30 #pragma weak __pow = pow
31 
32 /*
33  * pow(x,y) return x**y
34  *		      n
35  * Method:  Let x =  2   * (1+f)
36  *	1. Compute and return log2(x) in two pieces:
37  *		log2(x) = w1 + w2,
38  *	   where w1 has 24 bits trailing zero.
39  *	2. Perform y*log2(x) by simulating muti-precision arithmetic
40  *	3. Return x**y = exp2(y*log(x))
41  *
42  * Special cases:
43  *	1.  (anything) ** +-0 is 1
44  *	1'. 1 ** (anything)   is 1	(C99; 1 ** +-INF/NAN used to be NAN)
45  *	2.  (anything) ** 1   is itself
46  *	3.  (anything except 1) ** NAN is NAN ("except 1" is C99)
47  *	4.  NAN ** (anything except 0) is NAN
48  *	5.  +-(|x| > 1) **  +INF is +INF
49  *	6.  +-(|x| > 1) **  -INF is +0
50  *	7.  +-(|x| < 1) **  +INF is +0
51  *	8.  +-(|x| < 1) **  -INF is +INF
52  *	9.  -1          ** +-INF is 1	(C99; -1 ** +-INF used to be NAN)
53  *	10. +0 ** (+anything except 0, NAN)               is +0
54  *	11. -0 ** (+anything except 0, NAN, odd integer)  is +0
55  *	12. +0 ** (-anything except 0, NAN)               is +INF
56  *	13. -0 ** (-anything except 0, NAN, odd integer)  is +INF
57  *	14. -0 ** (odd integer) = -( +0 ** (odd integer) )
58  *	15. +INF ** (+anything except 0,NAN) is +INF
59  *	16. +INF ** (-anything except 0,NAN) is +0
60  *	17. -INF ** (anything)  = -0 ** (-anything)
61  *	18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
62  *	19. (-anything except 0 and inf) ** (non-integer) is NAN
63  *
64  * Accuracy:
65  *	pow(x,y) returns x**y nearly rounded. In particular
66  *			pow(integer,integer)
67  *	always returns the correct integer provided it is representable.
68  */
69 
70 #include "libm.h"
71 #include "xpg6.h"	/* __xpg6 */
72 #define	_C99SUSv3_pow	_C99SUSv3_pow_treats_Inf_as_an_even_int
73 
74 static const double zero = 0.0, one = 1.0, two = 2.0;
75 
76 extern const double _TBL_log2_hi[], _TBL_log2_lo[];
77 static const double
78 	two53 = 9007199254740992.0,
79 	A1_hi = 2.8853900432586669921875,
80 	A1_lo = 3.8519259825035041963606002e-8,
81 	A1 = 2.885390081777926817222541963606002026086e+0000,
82 	A2 = 9.617966939207270828380543979852286255862e-0001,
83 	A3 = 5.770807680887875964868853124873696201995e-0001,
84 	B0_hi = 2.8853900432586669921875,
85 	B0_lo = 3.8519259822532793056374320585e-8,
86 	B0 = 2.885390081777926814720293056374320585689e+0000,
87 	B1 = 9.617966939259755138949202350396200257632e-0001,
88 	B2 = 5.770780163585687000782112776448797953382e-0001,
89 	B3 = 4.121985488948771523290174512461778354953e-0001,
90 	B4 = 3.207590534812432970433641789022666850193e-0001;
91 
92 static double
log2_x(double x,double * w)93 log2_x(double x, double *w) {
94 	double f, s, z, qn, h, t;
95 	int *px = (int *) &x;
96 	int *pz = (int *) &z;
97 	int i, j, ix, n;
98 
99 	n = 0;
100 	ix = px[HIWORD];
101 	if (ix >= 0x3fef03f1 && ix < 0x3ff08208) {	/* 65/63 > x > 63/65 */
102 		double f1, v;
103 		f = x - one;
104 		if (((ix - 0x3ff00000) | px[LOWORD]) == 0) {
105 			*w = zero;
106 			return (zero);		/* log2(1)= +0 */
107 		}
108 		qn = one / (two + f);
109 		s = f * qn;				/* |s|<2**-6 */
110 		v = s * s;
111 		h = (double) ((float) s);
112 		f1 = (double) ((float) f);
113 		t = qn * (((f - two * h) - h * f1) - h * (f - f1));
114 								/* s = h+t */
115 		f1 = h * B0_lo + s * (v * (B1 + v * (B2 + v * (B3 + v * B4))));
116 		t = f1 + t * B0;
117 		h *= B0_hi;
118 		s = (double) ((float) (h + t));
119 		*w = t - (s - h);
120 		return (s);
121 	}
122 	if (ix < 0x00100000) {				/* subnormal x */
123 		x *= two53;
124 		n = -53;
125 		ix = px[HIWORD];
126 	}
127 	/* LARGE N */
128 	n += ((ix + 0x1000) >> 20) - 0x3ff;
129 	ix = (ix & 0x000fffff) | 0x3ff00000;		/* scale x to [1,2] */
130 	px[HIWORD] = ix;
131 	i = ix + 0x1000;
132 	pz[HIWORD] = i & 0xffffe000;
133 	pz[LOWORD] = 0;
134 	qn = one / (x + z);
135 	f = x - z;
136 	s = f * qn;
137 	h = (double) ((float) s);
138 	t = qn * ((f - (h + h) * z) - h * f);
139 	j = (i >> 13) & 0x7f;
140 	f = s * s;
141 	t = t * A1 + h * A1_lo;
142 	t += (s * f) * (A2 + f * A3);
143 	qn = h * A1_hi;
144 	s = n + _TBL_log2_hi[j];
145 	h = qn + s;
146 	t += _TBL_log2_lo[j] - ((h - s) - qn);
147 	f = (double) ((float) (h + t));
148 	*w = t - (f - h);
149 	return (f);
150 }
151 
152 extern const double _TBL_exp2_hi[], _TBL_exp2_lo[];
153 static const double		/* poly app of 2^x-1 on [-1e-10,2^-7+1e-10] */
154 	E1 = 6.931471805599453100674958533810346197328e-0001,
155 	E2 = 2.402265069587779347846769151717493815979e-0001,
156 	E3 = 5.550410866475410512631124892773937864699e-0002,
157 	E4 = 9.618143209991026824853712740162451423355e-0003,
158 	E5 = 1.333357676549940345096774122231849082991e-0003;
159 
160 double
pow(double x,double y)161 pow(double x, double y) {
162 	double z, ax;
163 	double y1, y2, w1, w2;
164 	int sbx, sby, j, k, yisint;
165 	int hx, hy, ahx, ahy;
166 	unsigned lx, ly;
167 	int *pz = (int *) &z;
168 
169 	hx = ((int *) &x)[HIWORD];
170 	lx = ((unsigned *) &x)[LOWORD];
171 	hy = ((int *) &y)[HIWORD];
172 	ly = ((unsigned *) &y)[LOWORD];
173 	ahx = hx & ~0x80000000;
174 	ahy = hy & ~0x80000000;
175 	if ((ahy | ly) == 0) {	/* y==zero  */
176 		if ((ahx | lx) == 0)
177 			z = _SVID_libm_err(x, y, 20);	/* +-0**+-0 */
178 		else if ((ahx | (((lx | -lx) >> 31) & 1)) > 0x7ff00000)
179 			z = _SVID_libm_err(x, y, 42);	/* NaN**+-0 */
180 		else
181 			z = one;			/* x**+-0 = 1 */
182 		return (z);
183 	} else if (hx == 0x3ff00000 && lx == 0 &&
184 		(__xpg6 & _C99SUSv3_pow) != 0)
185 		return (one);			/* C99: 1**anything = 1 */
186 	else if (ahx > 0x7ff00000 || (ahx == 0x7ff00000 && lx != 0) ||
187 		ahy > 0x7ff00000 || (ahy == 0x7ff00000 && ly != 0))
188 		return (x * y);	/* +-NaN return x*y; + -> * for Cheetah */
189 				/* includes Sun: 1**NaN = NaN */
190 	sbx = (unsigned) hx >> 31;
191 	sby = (unsigned) hy >> 31;
192 	ax = fabs(x);
193 
194 	/*
195 	 * determine if y is an odd int when x < 0
196 	 * yisint = 0 ... y is not an integer
197 	 * yisint = 1 ... y is an odd int
198 	 * yisint = 2 ... y is an even int
199 	 */
200 	yisint = 0;
201 	if (sbx) {
202 		if (ahy >= 0x43400000)
203 			yisint = 2;		/* even integer y */
204 		else if (ahy >= 0x3ff00000) {
205 			k = (ahy >> 20) - 0x3ff;	/* exponent */
206 			if (k > 20) {
207 				j = ly >> (52 - k);
208 				if ((j << (52 - k)) == ly)
209 					yisint = 2 - (j & 1);
210 			} else if (ly == 0) {
211 				j = ahy >> (20 - k);
212 				if ((j << (20 - k)) == ahy)
213 					yisint = 2 - (j & 1);
214 			}
215 		}
216 	}
217 	/* special value of y */
218 	if (ly == 0) {
219 		if (ahy == 0x7ff00000) {	/* y is +-inf */
220 			if (((ahx - 0x3ff00000) | lx) == 0) {
221 				if ((__xpg6 & _C99SUSv3_pow) != 0)
222 					return (one);
223 						/* C99: (-1)**+-inf = 1 */
224 				else
225 					return (y - y);
226 						/* Sun: (+-1)**+-inf = NaN */
227 			} else if (ahx >= 0x3ff00000)
228 						/* (|x|>1)**+,-inf = inf,0 */
229 				return (sby == 0 ? y : zero);
230 			else			/* (|x|<1)**-,+inf = inf,0 */
231 				return (sby != 0 ? -y : zero);
232 		}
233 		if (ahy == 0x3ff00000) {	/* y is  +-1 */
234 			if (sby != 0) {	/* y is -1 */
235 				if (x == zero)	/* divided by zero */
236 					return (_SVID_libm_err(x, y, 23));
237 				else if (ahx < 0x40000 || ((ahx - 0x40000) |
238 					lx) == 0)	/* overflow */
239 					return (_SVID_libm_err(x, y, 21));
240 				else
241 					return (one / x);
242 			} else
243 				return (x);
244 		}
245 		if (hy == 0x40000000) {		/* y is  2 */
246 			if (ahx >= 0x5ff00000 && ahx < 0x7ff00000)
247 				return (_SVID_libm_err(x, y, 21));
248 							/* x*x overflow */
249 			else if ((ahx < 0x1e56a09e && (ahx | lx) != 0) ||
250 				(ahx == 0x1e56a09e && lx < 0x667f3bcd))
251 				return (_SVID_libm_err(x, y, 22));
252 							/* x*x underflow */
253 			else
254 				return (x * x);
255 		}
256 		if (hy == 0x3fe00000) {
257 			if (!((ahx | lx) == 0 || ((ahx - 0x7ff00000) | lx) ==
258 				0 || sbx == 1))
259 				return (sqrt(x));	/* y is 0.5 and x > 0 */
260 		}
261 	}
262 	/* special value of x */
263 	if (lx == 0) {
264 		if (ahx == 0x7ff00000 || ahx == 0 || ahx == 0x3ff00000) {
265 			/* x is +-0,+-inf,-1 */
266 			z = ax;
267 			if (sby == 1) {
268 				z = one / z;	/* z = |x|**y */
269 				if (ahx == 0)
270 					return (_SVID_libm_err(x, y, 23));
271 			}
272 			if (sbx == 1) {
273 				if (ahx == 0x3ff00000 && yisint == 0)
274 					z = _SVID_libm_err(x, y, 24);
275 					/* neg**non-integral is NaN + invalid */
276 				else if (yisint == 1)
277 					z = -z;	/* (x<0)**odd = -(|x|**odd) */
278 			}
279 			return (z);
280 		}
281 	}
282 	/* (x<0)**(non-int) is NaN */
283 	if (sbx == 1 && yisint == 0)
284 		return (_SVID_libm_err(x, y, 24));
285 	/* Now ax is finite, y is finite */
286 	/* first compute log2(ax) = w1+w2, with 24 bits w1 */
287 	w1 = log2_x(ax, &w2);
288 
289 	/* split up y into y1+y2 and compute (y1+y2)*(w1+w2) */
290 	if (((ly & 0x07ffffff) == 0) || ahy >= 0x47e00000 ||
291 		ahy <= 0x38100000) {
292 		/* no need to split if y is short or too large or too small */
293 		y1 = y * w1;
294 		y2 = y * w2;
295 	} else {
296 		y1 = (double) ((float) y);
297 		y2 = (y - y1) * w1 + y * w2;
298 		y1 *= w1;
299 	}
300 	z = y1 + y2;
301 	j = pz[HIWORD];
302 	if (j >= 0x40900000) {				/* z >= 1024 */
303 		if (!(j == 0x40900000 && pz[LOWORD] == 0))	/* z > 1024 */
304 			return (_SVID_libm_err(x, y, 21));	/* overflow */
305 		else {
306 			w2 = y1 - z;
307 			w2 += y2;
308 							/* rounded to inf */
309 			if (w2 >= -8.008566259537296567160e-17)
310 				return (_SVID_libm_err(x, y, 21));
311 								/* overflow */
312 		}
313 	} else if ((j & ~0x80000000) >= 0x4090cc00) {	/* z <= -1075 */
314 		if (!(j == 0xc090cc00 && pz[LOWORD] == 0))	/* z < -1075 */
315 			return (_SVID_libm_err(x, y, 22));	/* underflow */
316 		else {
317 			w2 = y1 - z;
318 			w2 += y2;
319 			if (w2 <= zero)			/* underflow */
320 				return (_SVID_libm_err(x, y, 22));
321 		}
322 	}
323 	/*
324 	 * compute 2**(k+f[j]+g)
325 	 */
326 	k = (int) (z * 64.0 + (((hy ^ (ahx - 0x3ff00000)) > 0) ? 0.5 : -0.5));
327 	j = k & 63;
328 	w1 = y2 - ((double) k * 0.015625 - y1);
329 	w2 = _TBL_exp2_hi[j];
330 	z = _TBL_exp2_lo[j] + (w2 * w1) * (E1 + w1 * (E2 + w1 * (E3 + w1 *
331 		(E4 + w1 * E5))));
332 	z += w2;
333 	k >>= 6;
334 	if (k < -1021)
335 		z = scalbn(z, k);
336 	else			/* subnormal output */
337 		pz[HIWORD] += k << 20;
338 	if (sbx == 1 && yisint == 1)
339 		z = -z;		/* (-ve)**(odd int) */
340 	return (z);
341 }
342