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