xref: /illumos-gate/usr/src/lib/libm/i386/src/powf.S (revision 5d9d9091) 
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 * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
23 */
24/*
25 * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
26 * Use is subject to license terms.
27 */
28
29        .file "powf.s"
30
31/ Note: 0^SNaN should not signal "invalid" but this implementation
32/ does because y is placed on the NPX stack.
33
34/ Special cases:
35/
36/ x ** 0 is 1
37/ 1 ** y is 1				(C99)
38/ x ** NaN is NaN
39/ NaN ** y (except 0) is NaN
40/ x ** 1 is x
41/ +-(|x| > 1) **  +inf is +inf
42/ +-(|x| > 1) **  -inf is +0
43/ +-(|x| < 1) **  +inf is +0
44/ +-(|x| < 1) **  -inf is +inf
45/ (-1) ** +-inf is +1			(C99)
46/ +0 ** +y (except 0, NaN)              is +0
47/ -0 ** +y (except 0, NaN, odd int)     is +0
48/ +0 ** -y (except 0, NaN)              is +inf (z flag)
49/ -0 ** -y (except 0, NaN, odd int)     is +inf (z flag)
50/ -0 ** y (odd int)			is - (+0 ** x)
51/ +inf ** +y (except 0, NaN)    	is +inf
52/ +inf ** -y (except 0, NaN)    	is +0
53/ -inf ** +-y (except 0, NaN)   	is -0 ** -+y (NO z flag)
54/ x ** -1 is 1/x
55/ x ** 2 is x*x
56/ -x ** y (an integer) is (-1)**(y) * (+x)**(y)
57/ x ** y (x negative & y not integer) is NaN (i flag)
58
59#include "libm.h"
60LIBM_ANSI_PRAGMA_WEAK(powf,function)
61#include "libm_protos.h"
62#include "xpg6.h"
63
64	.data
65	.align	4
66negzero:
67	.float	-0.0
68half:
69	.float	0.5
70one:
71	.float	1.0
72negone:
73	.float	-1.0
74two:
75	.float	2.0
76Snan:
77	.long	0x7f800001
78pinfinity:
79	.long	0x7f800000
80ninfinity:
81	.long	0xff800000
82
83
84	ENTRY(powf)
85	pushl	%ebp
86	movl	%esp,%ebp
87	PIC_SETUP(1)
88
89	flds	8(%ebp)			/ x
90	fxam				/ determine class of x
91	fnstsw	%ax			/ store status in %ax
92	movb	%ah,%dh			/ %dh <- condition code of x
93
94	flds	12(%ebp)		/ y , x
95	fxam				/ determine class of y
96	fnstsw	%ax			/ store status in %ax
97	movb	%ah,%dl			/ %dl <- condition code of y
98
99	call	.pow_main		/// LOCAL
100	PIC_WRAPUP
101	leave
102	ret
103
104.pow_main:
105	/ x ** 0 is 1
106	movb	%dl,%cl
107	andb	$0x45,%cl
108	cmpb	$0x40,%cl		/ C3=1 C2=0 C1=? C0=0 when +-0
109	jne	1f
110	fstp	%st(0)			/ x
111	fstp	%st(0)			/ stack empty
112	fld1				/ 1
113	ret
114
1151:	/ y is not zero
116	PIC_G_LOAD(movzwl,__xpg6,eax)
117	andl	$_C99SUSv3_pow_treats_Inf_as_an_even_int,%eax
118	cmpl	$0,%eax
119	je	1f
120
121	/ C99: 1 ** anything is 1
122	fld1				/ 1, y, x
123	fucomp	%st(2)			/ y, x
124	fnstsw	%ax			/ store status in %ax
125	sahf				/ 80387 flags in %ax to 80386 flags
126	jp	1f			/ so that pow(NaN1,NaN2) returns NaN2
127	jne	1f
128	fstp	%st(0)			/ x
129	ret
130
1311:
132	/ x ** NaN is NaN
133	movb	%dl,%cl
134	andb	$0x45,%cl
135	cmpb	$0x01,%cl		/ C3=0 C2=0 C1=? C0=1 when +-NaN
136	jne	1f
137	fstp	%st(1)			/ y
138	ret
139
1401:	/ y is not NaN
141	/ NaN ** y (except 0) is NaN
142	movb	%dh,%cl
143	andb	$0x45,%cl
144	cmpb	$0x01,%cl		/ C3=0 C2=0 C1=? C0=1 when +-NaN
145	jne	1f
146	fstp	%st(0)			/ x
147	ret
148
1491:	/ x is not NaN
150	/ x ** 1 is x
151	fcoms	PIC_L(one)		/ y , x
152	fnstsw	%ax			/ store status in %ax
153	sahf				/ 80387 flags in %ax to 80386 flags
154	jne	1f
155	fstp	%st(0)			/ x
156	ret
157
1581:	/ y is not 1
159	/ +-(|x| > 1) **  +inf is +inf
160	/ +-(|x| > 1) **  -inf is +0
161	/ +-(|x| < 1) **  +inf is +0
162	/ +-(|x| < 1) **  -inf is +inf
163	/ +-(|x| = 1) ** +-inf is NaN
164	movb	%dl,%cl
165	andb	$0x47,%cl
166	cmpb	$0x05,%cl		/ C3=0 C2=1 C1=0 C0=1 when +inf
167	je	.yispinf
168	cmpb	$0x07,%cl		/ C3=0 C2=1 C1=1 C0=1 when -inf
169	je	.yisninf
170
171	/ +0 ** +y (except 0, NaN)		is +0
172	/ -0 ** +y (except 0, NaN, odd int)	is +0
173	/ +0 ** -y (except 0, NaN)		is +inf (z flag)
174	/ -0 ** -y (except 0, NaN, odd int)	is +inf (z flag)
175	/ -0 ** y (odd int)			is - (+0 ** x)
176	movb	%dh,%cl
177	andb	$0x47,%cl
178	cmpb	$0x40,%cl		/ C3=1 C2=0 C1=0 C0=0 when +0
179	je	.xispzero
180	cmpb	$0x42,%cl		/ C3=1 C2=0 C1=1 C0=0 when -0
181	je	.xisnzero
182
183	/ +inf ** +y (except 0, NaN)	is +inf
184	/ +inf ** -y (except 0, NaN)	is +0
185	/ -inf ** +-y (except 0, NaN)	is -0 ** -+y (NO z flag)
186	movb	%dh,%cl
187	andb	$0x47,%cl
188	cmpb	$0x05,%cl		/ C3=0 C2=1 C1=0 C0=1 when +inf
189	je	.xispinf
190	cmpb	$0x07,%cl		/ C3=0 C2=1 C1=1 C0=1 when -inf
191	je	.xisninf
192
193	/ x ** -1 is 1/x
194	fcoms	PIC_L(negone)		/ y , x
195	fnstsw	%ax			/ store status in %ax
196	sahf				/ 80387 flags in %ax to 80386 flags
197	jne	1f
198	fld	%st(1)			/ x , y , x
199	fdivrs	PIC_L(one)		/ 1/x , y , x
200	jmp	.signok			/ check for over/underflow
201
2021:	/ y is not -1
203	/ x ** 2 is square(x)
204	fcoms	PIC_L(two)		/ y , x
205	fnstsw	%ax			/ store status in %ax
206	sahf				/ 80387 flags in %ax to 80386 flags
207	jne	1f
208	fld	%st(1)			/ x , y , x
209	fld	%st(0)			/ x , x , y , x
210	fmulp				/ x^2 , y , x
211	jmp	.signok			/ check for over/underflow
212
2131:	/ y is not 2
214	/ x ** 1/2 is sqrt(x)
215	fcoms	PIC_L(half)		/ y , x
216	fnstsw	%ax			/ store status in %ax
217	sahf				/ 80387 flags in %ax to 80386 flags
218	jne	1f
219	fld	%st(1)			/ x , y , x
220	fsqrt				/ sqrt(x) , y , x
221	jmp	.signok			/ check for over/underflow
222
2231:	/ y is not 2
224	/ make copies of x & y
225	fld	%st(1)			/ x , y , x
226	fld	%st(1)			/ y , x , y , x
227
228	/ -x ** y (an integer) is (-1)**(y) * (+x)**(y)
229	/ x ** y (x negative & y not integer) is  NaN
230	movl	$0,%ecx			/ track whether to flip sign of result
231	fld	%st(1)			/ x , y , x , y , x
232	ftst				/ compare %st(0) with 0
233	fnstsw	%ax			/ store status in %ax
234	sahf				/ 80387 flags in %ax to 80386 flags
235	fstp	%st(0)			/ y , x , y , x
236	ja	.merge			/ x > 0
237	/ x < 0
238	call	.y_is_int
239	cmpl	$0,%ecx
240	jne	1f
241	/ x < 0 & y != int so x**y = NaN (i flag)
242	fstp	%st(0)			/ x , y , x
243	fstp	%st(0)			/ y , x
244	fstp	%st(0)			/ y , x
245	fstp	%st(0)			/ y , x
246	fldz
247	fdiv	%st,%st(0)		/ 0/0
248	ret
249
2501:	/ x < 0 & y = int
251	fxch				/ x , y , y , x
252	fchs				/ px = -x , y , y , x
253	fxch				/ y , px , y , x
254.merge:
255	/ px > 0
256	fxch				/ px , y , y , x
257
258	/ x**y   =   exp(y*ln(x))
259	fyl2x				/ t=y*log2(px) , y , x
260	fld	%st(0)			/ t , t , y , x
261	frndint				/ [t] , t , y , x
262	fxch				/ t , [t] , y , x
263	fucom
264	fnstsw	%ax			/ store status in %ax
265	sahf				/ 80387 flags in %ax to 80386 flags
266	je	1f			/ px = int
267	fsub    %st(1),%st		/ t-[t] , [t] , y , x
268	f2xm1				/ 2**(t-[t])-1 , [t] , y , x
269	fadds	PIC_L(one)		/ 2**(t-[t]) , [t] , y , x
270	fscale				/ 2**t = px**y , [t] , y , x
271	jmp	2f
2721:
273	fstp    %st(0)                  / t=[t] , y , x
274	fld1                            / 1 , t , y , x
275	fscale                          / 1*2**t = x**y , t , y , x
2762:
277	fstp	%st(1)			/ x**y , y , x
278	cmpl	$1,%ecx
279	jne	.signok
280	fchs				/ change sign since x<0 & y=-int
281.signok:
282	subl	$4,%esp
283	fstps	(%esp)			/ round to single precision
284	flds	(%esp)			/ place result on NPX stack
285	addl	$4,%esp
286	fstp	%st(2)			/ y , x**y
287	fstp	%st(0)			/ x**y
288	ret
289
290/ ------------------------------------------------------------------------
291
292.xispinf:
293	ftst				/ compare %st(0) with 0
294	fnstsw	%ax			/ store status in %ax
295	sahf				/ 80387 flags in %ax to 80386 flags
296	ja	.retpinf		/ y > 0
297	jmp	.retpzero		/ y < 0
298
299.xisninf:
300	/ -inf ** +-y is -0 ** -+y
301	fchs				/ -y , x
302	flds	PIC_L(negzero)		/ -0 , -y , x
303	fstp	%st(2)			/ -y , -0
304	jmp	.xisnzero
305
306.yispinf:
307	fld	%st(1)			/ x , y , x
308	fabs				/ |x| , y , x
309	fcomps	PIC_L(one)		/ y , x
310	fnstsw	%ax			/ store status in %ax
311	sahf				/ 80387 flags in %ax to 80386 flags
312	je	.retponeorinvalid	/ x == -1	C99
313	ja	.retpinf		/ |x| > 1
314	jmp	.retpzero		/ |x| < 1
315
316.yisninf:
317	fld	%st(1)			/ x , y , x
318	fabs				/ |x| , y , x
319	fcomps	PIC_L(one)		/ y , x
320	fnstsw	%ax			/ store status in %ax
321	sahf				/ 80387 flags in %ax to 80386 flags
322	je	.retponeorinvalid	/ x == -1	C99
323	ja	.retpzero		/ |x| > 1
324	jmp	.retpinf		/ |x| < 1
325
326.xispzero:
327	/ y cannot be 0 or NaN ; stack has	y , x
328	ftst				/ compare %st(0) with 0
329	fnstsw	%ax			/ store status in %ax
330	sahf				/ 80387 flags in %ax to 80386 flags
331	ja	.retpzero		/ y > 0
332	/ x = +0 & y < 0 so x**y = +inf
333	jmp	.retpinfzflag		/ ret +inf & z flag
334
335.xisnzero:
336	/ y cannot be 0 or NaN ; stack has	y , x
337	call	.y_is_int
338	cmpl	$1,%ecx
339	jne	1f			/ y is not an odd integer
340	/ y is an odd integer
341	ftst				/ compare %st(0) with 0
342	fnstsw	%ax			/ store status in %ax
343	sahf				/ 80387 flags in %ax to 80386 flags
344	ja	.retnzero		/ y > 0
345	/ x = -0 & y < 0 (odd int)	return -inf (z flag)
346	/ x = -inf & y != 0 or NaN	return -inf (NO z flag)
347	movb	%dh,%cl
348	andb	$0x45,%cl
349	cmpb	$0x05,%cl		/ C3=0 C2=1 C1=? C0=1 when +-inf
350	je	2f
351	fdiv	%st,%st(1)		/ y / x, x (raise z flag)
3522:
353	fstp	%st(0)			/ x
354	fstp	%st(0)			/ stack empty
355	flds	PIC_L(ninfinity)	/ -inf
356	ret
357
3581:	/ y is not an odd integer
359	ftst				/ compare %st(0) with 0
360	fnstsw	%ax			/ store status in %ax
361	sahf				/ 80387 flags in %ax to 80386 flags
362	ja	.retpzero		/ y > 0
363	/ x = -0 & y < 0 (not odd int)	return +inf (z flag)
364	/ x = -inf & y not 0 or NaN 	return +inf (NO z flag)
365	movb	%dh,%cl
366	andb	$0x45,%cl
367	cmpb	$0x05,%cl		/ C3=0 C2=1 C1=? C0=1 when +-inf
368	jne	.retpinfzflag		/ ret +inf & divide-by-0 flag
369	jmp	.retpinf		/ return +inf (NO z flag)
370
371.retpzero:
372	fstp	%st(0)			/ x
373	fstp	%st(0)			/ stack empty
374	fldz				/ +0
375	ret
376
377.retnzero:
378	fstp	%st(0)			/ x
379	fstp	%st(0)			/ stack empty
380	flds	PIC_L(negzero)		/ -0
381	ret
382
383.retponeorinvalid:
384	PIC_G_LOAD(movzwl,__xpg6,eax)
385	andl	$_C99SUSv3_pow_treats_Inf_as_an_even_int,%eax
386	cmpl	$0,%eax
387	je	1f
388	fstp	%st(0)			/ x
389	fstp	%st(0)			/ stack empty
390	fld1				/ 1
391	ret
392
3931:
394	fstp	%st(0)			/ x
395	fstp	%st(0)			/ stack empty
396	flds	PIC_L(Snan)		/ Q NaN (i flag)
397	fwait
398	ret
399
400.retpinf:
401	fstp	%st(0)			/ x
402	fstp	%st(0)			/ stack empty
403	flds	PIC_L(pinfinity)	/ +inf
404	ret
405
406.retpinfzflag:
407	fstp	%st(0)			/ x
408	fstp	%st(0)			/ stack empty
409	fldz
410	fdivrs	PIC_L(one)		/ 1/0
411	ret
412
413/ Set %ecx to 2 if y is an even integer, 1 if y is an odd integer,
414/ 0 otherwise.  Assume y is not zero.  Do not raise inexact or modify
415/ %edx.
416.y_is_int:
417	movl	12(%ebp),%eax
418	andl	$0x7fffffff,%eax	/ |y|
419	cmpl	$0x4b800000,%eax
420	jae	1f			/ |y| >= 2^24, an even int
421	cmpl	$0x3f800000,%eax
422	jb	2f			/ |y| < 1, can't be an int
423	movl	%eax,%ecx
424	sarl	$23,%ecx
425	subl	$150,%ecx
426	negl	%ecx			/ 23 - unbiased exponent of y
427	bsfl	%eax,%eax		/ index of least sig. 1 bit
428	cmpl	%ecx,%eax
429	jb	2f
430	ja	1f
431	movl	$1,%ecx
432	ret
4331:
434	movl	$2,%ecx
435	ret
4362:
437	xorl	%ecx,%ecx
438	ret
439	.align	4
440	SET_SIZE(powf)
441