xref: /illumos-gate/usr/src/lib/libm/common/Q/powl.c (revision 25c28e83)
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 #if defined(ELFOBJ)
31 #pragma weak powl = __powl
32 #endif
33 
34 #include "libm.h"
35 #include "xpg6.h"	/* __xpg6 */
36 #define	_C99SUSv3_pow	_C99SUSv3_pow_treats_Inf_as_an_even_int
37 
38 #if defined(__sparc)
39 #define	i0	0
40 #define	i1	1
41 #define	i2	2
42 #define	i3	3
43 
44 static const long double zero = 0.0L, one = 1.0L, two = 2.0L;
45 
46 extern const long double _TBL_logl_hi[], _TBL_logl_lo[];
47 
48 static const long double
49 	two113 = 10384593717069655257060992658440192.0L,
50 	ln2hi = 6.931471805599453094172319547495844850203e-0001L,
51 	ln2lo = 1.667085920830552208890449330400379754169e-0025L,
52 	A2 = 6.666666666666666666666666666666091393804e-0001L,
53 	A3 = 4.000000000000000000000000407167070220671e-0001L,
54 	A4 = 2.857142857142857142730077490612903681164e-0001L,
55 	A5 = 2.222222222222242577702836920812882605099e-0001L,
56 	A6 = 1.818181816435493395985912667105885828356e-0001L,
57 	A7 = 1.538537835211839751112067512805496931725e-0001L,
58 	B1 = 6.666666666666666666666666666666666667787e-0001L,
59 	B2 = 3.999999999999999999999999999999848524411e-0001L,
60 	B3 = 2.857142857142857142857142865084581075070e-0001L,
61 	B4 = 2.222222222222222222222010781800643808497e-0001L,
62 	B5 = 1.818181818181818185051442171337036403674e-0001L,
63 	B6 = 1.538461538461508363540720286292008207673e-0001L,
64 	B7 = 1.333333333506731842033180638329317108428e-0001L,
65 	B8 = 1.176469984587418890634302788283946761670e-0001L,
66 	B9 = 1.053794891561452331722969901564862497132e-0001L;
67 
68 static long double
69 logl_x(long double x, long double *w) {
70 	long double f, f1, v, s, z, qn, h, t;
71 	int *px = (int *) &x;
72 	int *pz = (int *) &z;
73 	int i, j, ix, n;
74 
75 	n = 0;
76 	ix = px[i0];
77 	if (ix > 0x3ffef03f && ix < 0x3fff0820) {	/* 65/63 > x > 63/65 */
78 		f = x - one;
79 		z = f * f;
80 		if (((ix - 0x3fff0000) | px[i1] | px[i2] | px[i3]) == 0) {
81 			*w = zero;
82 			return (zero);	/* log(1)= +0 */
83 		}
84 		qn = one / (two + f);
85 		s = f * qn;	/* |s|<2**-6 */
86 		v = s * s;
87 		h = (long double) (2.0 * (double) s);
88 		f1 = (long double) ((double) f);
89 		t = ((two * (f - h) - h * f1) - h * (f - f1)) * qn +
90 			s * (v * (B1 + v * (B2 + v * (B3 + v * (B4 +
91 			v * (B5 + v * (B6 + v * (B7 + v * (B8 + v * B9)))))))));
92 		s = (long double) ((double) (h + t));
93 		*w = t - (s - h);
94 		return (s);
95 	}
96 	if (ix < 0x00010000) {	/* subnormal x */
97 		x *= two113;
98 		n = -113;
99 		ix = px[i0];
100 	}
101 	/* LARGE_N */
102 	n += ((ix + 0x200) >> 16) - 0x3fff;
103 	ix = (ix & 0x0000ffff) | 0x3fff0000;	/* scale x to [1,2] */
104 	px[i0] = ix;
105 	i = ix + 0x200;
106 	pz[i0] = i & 0xfffffc00;
107 	pz[i1] = pz[i2] = pz[i3] = 0;
108 	qn = one / (x + z);
109 	f = x - z;
110 	s = f * qn;
111 	f1 = (long double) ((double) f);
112 	h = (long double) (2.0 * (double) s);
113 	t = qn * ((two * (f - z * h) - h * f1) - h * (f - f1));
114 	j = (i >> 10) & 0x3f;
115 	v = s * s;
116 	qn = (long double) n;
117 	t += qn * ln2lo + _TBL_logl_lo[j];
118 	t += s * (v * (A2 + v * (A3 + v * (A4 + v * (A5 + v * (A6 +
119 		v * A7))))));
120 	v = qn * ln2hi + _TBL_logl_hi[j];
121 	s = h + v;
122 	t += (h - (s - v));
123 	z = (long double) ((double) (s + t));
124 	*w = t - (z - s);
125 	return (z);
126 }
127 
128 extern const long double _TBL_expl_hi[], _TBL_expl_lo[];
129 static const long double
130 	invln2_32 = 4.616624130844682903551758979206054839765e+1L,
131 	ln2_32hi = 2.166084939249829091928849858592451515688e-2L,
132 	ln2_32lo = 5.209643502595475652782654157501186731779e-27L,
133 	ln2_64 = 1.083042469624914545964425189778400898568e-2L;
134 
135 long double
136 powl(long double x, long double y) {
137 	long double z, ax;
138 	long double y1, y2, w1, w2;
139 	int sbx, sby, j, k, yisint, m;
140 	int hx, lx, hy, ly, ahx, ahy;
141 	int *pz = (int *) &z;
142 	int *px = (int *) &x;
143 	int *py = (int *) &y;
144 
145 	hx = px[i0];
146 	lx = px[i1] | px[i2] | px[i3];
147 	hy = py[i0];
148 	ly = py[i1] | py[i2] | py[i3];
149 	ahx = hx & ~0x80000000;
150 	ahy = hy & ~0x80000000;
151 
152 	if ((ahy | ly) == 0)
153 		return (one);		/* x**+-0 = 1 */
154 	else if (hx == 0x3fff0000 && lx == 0 &&
155 		(__xpg6 & _C99SUSv3_pow) != 0)
156 		return (one);		/* C99: 1**anything = 1 */
157 	else if (ahx > 0x7fff0000 || (ahx == 0x7fff0000 && lx != 0) ||
158 		ahy > 0x7fff0000 || (ahy == 0x7fff0000 && ly != 0))
159 		return (x + y);		/* +-NaN return x+y */
160 					/* includes Sun: 1**NaN = NaN */
161 	sbx = (unsigned) hx >> 31;
162 	sby = (unsigned) hy >> 31;
163 	ax = fabsl(x);
164 	/*
165 	 * determine if y is an odd int when x < 0
166 	 * yisint = 0 ... y is not an integer
167 	 * yisint = 1 ... y is an odd int
168 	 * yisint = 2 ... y is an even int
169 	 */
170 	yisint = 0;
171 	if (sbx) {
172 		if (ahy >= 0x40700000)	/* if |y|>=2**113 */
173 			yisint = 2;	/* even integer y */
174 		else if (ahy >= 0x3fff0000) {
175 			k = (ahy >> 16) - 0x3fff;	/* exponent */
176 			if (k > 80) {
177 				j = ((unsigned) py[i3]) >> (112 - k);
178 				if ((j << (112 - k)) == py[i3])
179 					yisint = 2 - (j & 1);
180 			} else if (k > 48) {
181 				j = ((unsigned) py[i2]) >> (80 - k);
182 				if ((j << (80 - k)) == py[i2])
183 					yisint = 2 - (j & 1);
184 			} else if (k > 16) {
185 				j = ((unsigned) py[i1]) >> (48 - k);
186 				if ((j << (48 - k)) == py[i1])
187 					yisint = 2 - (j & 1);
188 			} else if (ly == 0) {
189 				j = ahy >> (16 - k);
190 				if ((j << (16 - k)) == ahy)
191 					yisint = 2 - (j & 1);
192 			}
193 		}
194 	}
195 
196 	/* special value of y */
197 	if (ly == 0) {
198 		if (ahy == 0x7fff0000) {	/* y is +-inf */
199 			if (((ahx - 0x3fff0000) | lx) == 0) {
200 				if ((__xpg6 & _C99SUSv3_pow) != 0)
201 					return (one);
202 						/* C99: (-1)**+-inf = 1 */
203 				else
204 					return (y - y);
205 						/* Sun: (+-1)**+-inf = NaN */
206 			} else if (ahx >= 0x3fff0000)
207 						/* (|x|>1)**+,-inf = inf,0 */
208 				return (sby == 0 ? y : zero);
209 			else			/* (|x|<1)**-,+inf = inf,0 */
210 				return (sby != 0 ? -y : zero);
211 		} else if (ahy == 0x3fff0000) {	/* y is +-1 */
212 			if (sby != 0)
213 				return (one / x);
214 			else
215 				return (x);
216 		} else if (hy == 0x40000000)	/* y is 2 */
217 			return (x * x);
218 		else if (hy == 0x3ffe0000) {	/* y is 0.5 */
219 			if (!((ahx | lx) == 0 || ((ahx - 0x7fff0000) | lx) ==
220 				0))
221 				return (sqrtl(x));
222 		}
223 	}
224 
225 	/* special value of x */
226 	if (lx == 0) {
227 		if (ahx == 0x7fff0000 || ahx == 0 || ahx == 0x3fff0000) {
228 							/* x is +-0,+-inf,+-1 */
229 			z = ax;
230 			if (sby == 1)
231 				z = one / z;	/* z = 1/|x| if y is negative */
232 			if (sbx == 1) {
233 				if (ahx == 0x3fff0000 && yisint == 0)
234 					z = zero / zero;
235 						/* (-1)**non-int is NaN */
236 				else if (yisint == 1)
237 					z = -z;	/* (x<0)**odd = -(|x|**odd) */
238 			}
239 			return (z);
240 		}
241 	}
242 
243 	/* (x<0)**(non-int) is NaN */
244 	if (sbx == 1 && yisint == 0)
245 		return (zero / zero);	/* should be volatile */
246 
247 	/* Now ax is finite, y is finite */
248 	/* first compute log(ax) = w1+w2, with 53 bits w1 */
249 	w1 = logl_x(ax, &w2);
250 
251 	/* split up y into y1+y2 and compute (y1+y2)*(w1+w2) */
252 	if (ly == 0 || ahy >= 0x43fe0000) {
253 		y1 = y * w1;
254 		y2 = y * w2;
255 	} else {
256 		y1 = (long double) ((double) y);
257 		y2 = (y - y1) * w1 + y * w2;
258 		y1 *= w1;
259 	}
260 	z = y1 + y2;
261 	j = pz[i0];
262 	if ((unsigned) j >= 0xffff0000) {		/* NaN or -inf */
263 		if (sbx == 1 && yisint == 1)
264 			return (one / z);
265 		else
266 			return (-one / z);
267 	} else if ((j & ~0x80000000) < 0x3fc30000) {	/* |x|<2^-60 */
268 		if (sbx == 1 && yisint == 1)
269 			return (-one - z);
270 		else
271 			return (one + z);
272 	} else if (j > 0) {
273 		if (j > 0x400d0000) {
274 			if (sbx == 1 && yisint == 1)
275 				return (scalbnl(-one, 20000));
276 			else
277 				return (scalbnl(one, 20000));
278 		}
279 		k = (int) (invln2_32 * (z + ln2_64));
280 	} else {
281 		if ((unsigned) j > 0xc00d0000) {
282 			if (sbx == 1 && yisint == 1)
283 				return (scalbnl(-one, -20000));
284 			else
285 				return (scalbnl(one, -20000));
286 		}
287 		k = (int) (invln2_32 * (z - ln2_64));
288 	}
289 	j = k & 0x1f;
290 	m = k >> 5;
291 	{
292 		/* rational approximation coeffs for [-(ln2)/64,(ln2)/64] */
293 		long double
294 			t1 = 1.666666666666666666666666666660876387437e-1L,
295 			t2 = -2.777777777777777777777707812093173478756e-3L,
296 			t3 = 6.613756613756613482074280932874221202424e-5L,
297 			t4 = -1.653439153392139954169609822742235851120e-6L,
298 			t5 = 4.175314851769539751387852116610973796053e-8L;
299 		long double t = (long double) k;
300 
301 		w1 = (y2 - (t * ln2_32hi - y1)) - t * ln2_32lo;
302 		t = w1 * w1;
303 		w2 = (w1 - t * (t1 + t * (t2 + t * (t3 + t * (t4 + t * t5))))) -
304 			two;
305 		z = _TBL_expl_hi[j] - ((_TBL_expl_hi[j] * (w1 + w1)) / w2 -
306 			_TBL_expl_lo[j]);
307 	}
308 	j = m + (pz[i0] >> 16);
309 	if (j && (unsigned) j < 0x7fff)
310 		pz[i0] += m << 16;
311 	else
312 		z = scalbnl(z, m);
313 
314 	if (sbx == 1 && yisint == 1)
315 		z = -z;		/* (-ve)**(odd int) */
316 	return (z);
317 }
318 #else
319 #error Unsupported Architecture
320 #endif	/* defined(__sparc) */
321