1/* ix87 specific implementation of pow function.
2   Copyright (C) 1996, 1997, 1998, 1999, 2001, 2004, 2007
3   Free Software Foundation, Inc.
4   This file is part of the GNU C Library.
5   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996.
6
7   The GNU C Library is free software; you can redistribute it and/or
8   modify it under the terms of the GNU Lesser General Public
9   License as published by the Free Software Foundation; either
10   version 2.1 of the License, or (at your option) any later version.
11
12   The GNU C Library is distributed in the hope that it will be useful,
13   but WITHOUT ANY WARRANTY; without even the implied warranty of
14   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15   Lesser General Public License for more details.
16
17   You should have received a copy of the GNU Lesser General Public
18   License along with the GNU C Library; if not, write to the Free
19   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
20   02111-1307 USA.  */
21
22#include <machine/asm.h>
23
24#ifdef __ELF__
25	.section .rodata
26#else
27	.text
28#endif
29
30	.align ALIGNARG(4)
31	ASM_TYPE_DIRECTIVE(infinity,@object)
32inf_zero:
33infinity:
34	.byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
35	ASM_SIZE_DIRECTIVE(infinity)
36	ASM_TYPE_DIRECTIVE(zero,@object)
37zero:	.double 0.0
38	ASM_SIZE_DIRECTIVE(zero)
39	ASM_TYPE_DIRECTIVE(minf_mzero,@object)
40minf_mzero:
41minfinity:
42	.byte 0, 0, 0, 0, 0, 0, 0xf0, 0xff
43mzero:
44	.byte 0, 0, 0, 0, 0, 0, 0, 0x80
45	ASM_SIZE_DIRECTIVE(minf_mzero)
46	ASM_TYPE_DIRECTIVE(one,@object)
47one:	.double 1.0
48	ASM_SIZE_DIRECTIVE(one)
49	ASM_TYPE_DIRECTIVE(limit,@object)
50limit:	.double 0.29
51	ASM_SIZE_DIRECTIVE(limit)
52	ASM_TYPE_DIRECTIVE(p63,@object)
53p63:
54	.byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43
55	ASM_SIZE_DIRECTIVE(p63)
56
57#ifdef PIC
58#define MO(op) op##(%rip)
59#else
60#define MO(op) op
61#endif
62
63	.text
64ENTRY(__ieee754_powl)
65	fldt	24(%rsp)	// y
66	fxam
67
68
69	fnstsw
70	movb	%ah, %dl
71	andb	$0x45, %ah
72	cmpb	$0x40, %ah	// is y == 0 ?
73	je	11f
74
75	cmpb	$0x05, %ah	// is y == �inf ?
76	je	12f
77
78	cmpb	$0x01, %ah	// is y == NaN ?
79	je	30f
80
81	fldt	8(%rsp)		// x : y
82
83	fxam
84	fnstsw
85	movb	%ah, %dh
86	andb	$0x45, %ah
87	cmpb	$0x40, %ah
88	je	20f		// x is �0
89
90	cmpb	$0x05, %ah
91	je	15f		// x is �inf
92
93	fxch			// y : x
94
95	/* fistpll raises invalid exception for |y| >= 1L<<63.  */
96	fldl	MO(p63)		// 1L<<63 : y : x
97	fld	%st(1)		// y : 1L<<63 : y : x
98	fabs			// |y| : 1L<<63 : y : x
99	fcomip	%st(1), %st	// 1L<<63 : y : x
100	fstp	%st(0)		// y : x
101	jnc	2f
102
103	/* First see whether `y' is a natural number.  In this case we
104	   can use a more precise algorithm.  */
105	fld	%st		// y : y : x
106	fistpll	-8(%rsp)	// y : x
107	fildll	-8(%rsp)	// int(y) : y : x
108	fucomip	%st(1),%st	// y : x
109	jne	2f
110
111	/* OK, we have an integer value for y.  */
112	mov	-8(%rsp),%eax
113	mov	-4(%rsp),%edx
114	orl	$0, %edx
115	fstp	%st(0)		// x
116	jns	4f		// y >= 0, jump
117	fdivrl	MO(one)		// 1/x		(now referred to as x)
118	negl	%eax
119	adcl	$0, %edx
120	negl	%edx
1214:	fldl	MO(one)		// 1 : x
122	fxch
123
1246:	shrdl	$1, %edx, %eax
125	jnc	5f
126	fxch
127	fmul	%st(1)		// x : ST*x
128	fxch
1295:	fmul	%st(0), %st	// x*x : ST*x
130	shrl	$1, %edx
131	movl	%eax, %ecx
132	orl	%edx, %ecx
133	jnz	6b
134	fstp	%st(0)		// ST*x
135	ret
136
137	/* y is �NAN */
13830:	fldt	8(%rsp)		// x : y
139	fldl	MO(one)		// 1.0 : x : y
140	fucomip	%st(1),%st	// x : y
141	je	31f
142	fxch			// y : x
14331:	fstp	%st(1)
144	ret
145
146	.align ALIGNARG(4)
1472:	/* y is a real number.  */
148	fxch			// x : y
149	fldl	MO(one)		// 1.0 : x : y
150	fldl	MO(limit)	// 0.29 : 1.0 : x : y
151	fld	%st(2)		// x : 0.29 : 1.0 : x : y
152	fsub	%st(2)		// x-1 : 0.29 : 1.0 : x : y
153	fabs			// |x-1| : 0.29 : 1.0 : x : y
154	fucompp			// 1.0 : x : y
155	fnstsw
156	fxch			// x : 1.0 : y
157	test	$4500,%eax
158	jz	7f
159	fsub	%st(1)		// x-1 : 1.0 : y
160	fyl2xp1			// log2(x) : y
161	jmp	8f
162
1637:	fyl2x			// log2(x) : y
1648:	fmul	%st(1)		// y*log2(x) : y
165	fxam
166	fnstsw
167	andb	$0x45, %ah
168	cmpb	$0x05, %ah      // is y*log2(x) == �inf ?
169	je	28f
170	fst	%st(1)		// y*log2(x) : y*log2(x)
171	frndint			// int(y*log2(x)) : y*log2(x)
172	fsubr	%st, %st(1)	// int(y*log2(x)) : fract(y*log2(x))
173	fxch			// fract(y*log2(x)) : int(y*log2(x))
174	f2xm1			// 2^fract(y*log2(x))-1 : int(y*log2(x))
175	faddl	MO(one)		// 2^fract(y*log2(x)) : int(y*log2(x))
176	fscale			// 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
177	fstp	%st(1)		// 2^fract(y*log2(x))*2^int(y*log2(x))
178	ret
179
18028:	fstp	%st(1)		// y*log2(x)
181	fldl	MO(one)		// 1 : y*log2(x)
182	fscale			// 2^(y*log2(x)) : y*log2(x)
183	fstp	%st(1)		// 2^(y*log2(x))
184	ret
185
186	// pow(x,�0) = 1
187	.align ALIGNARG(4)
18811:	fstp	%st(0)		// pop y
189	fldl	MO(one)
190	ret
191
192	// y == �inf
193	.align ALIGNARG(4)
19412:	fstp	%st(0)		// pop y
195	fldl	MO(one)		// 1
196	fldt	8(%rsp)		// x : 1
197	fabs			// abs(x) : 1
198	fucompp			// < 1, == 1, or > 1
199	fnstsw
200	andb	$0x45, %ah
201	cmpb	$0x45, %ah
202	je	13f		// jump if x is NaN
203
204	cmpb	$0x40, %ah
205	je	14f		// jump if |x| == 1
206
207	shlb	$1, %ah
208	xorb	%ah, %dl
209	andl	$2, %edx
210#ifdef PIC
211	lea	inf_zero(%rip),%rcx
212	fldl	(%rcx, %rdx, 4)
213#else
214	fldl	inf_zero(,%rdx, 4)
215#endif
216	ret
217
218	.align ALIGNARG(4)
21914:	fldl	MO(one)
220	ret
221
222	.align ALIGNARG(4)
22313:	fldt	8(%rsp)		// load x == NaN
224	ret
225
226	.align ALIGNARG(4)
227	// x is �inf
22815:	fstp	%st(0)		// y
229	testb	$2, %dh
230	jz	16f		// jump if x == +inf
231
232	// We must find out whether y is an odd integer.
233	fld	%st		// y : y
234	fistpll	-8(%rsp)	// y
235	fildll	-8(%rsp)	// int(y) : y
236	fucomip %st(1),%st
237	ffreep	%st		// <empty>
238	jne	17f
239
240	// OK, the value is an integer, but is it odd?
241	mov	-8(%rsp), %eax
242	mov	-4(%rsp), %edx
243	andb	$1, %al
244	jz	18f		// jump if not odd
245	// It's an odd integer.
246	shrl	$31, %edx
247#ifdef PIC
248	lea	minf_mzero(%rip),%rcx
249	fldl	(%rcx, %rdx, 8)
250#else
251	fldl	minf_mzero(,%rdx, 8)
252#endif
253	ret
254
255	.align ALIGNARG(4)
25616:	fcompl	MO(zero)
257	fnstsw
258	shrl	$5, %eax
259	andl	$8, %eax
260#ifdef PIC
261	lea	inf_zero(%rip),%rcx
262	fldl	(%rcx, %rax, 1)
263#else
264	fldl	inf_zero(,%rax, 1)
265#endif
266	ret
267
268	.align ALIGNARG(4)
26917:	shll	$30, %edx	// sign bit for y in right position
27018:	shrl	$31, %edx
271#ifdef PIC
272	lea	inf_zero(%rip),%rcx
273	fldl	(%rcx, %rdx, 8)
274#else
275	fldl	inf_zero(,%rdx, 8)
276#endif
277	ret
278
279	.align ALIGNARG(4)
280	// x is �0
28120:	fstp	%st(0)		// y
282	testb	$2, %dl
283	jz	21f		// y > 0
284
285	// x is �0 and y is < 0.  We must find out whether y is an odd integer.
286	testb	$2, %dh
287	jz	25f
288
289	fld	%st		// y : y
290	fistpll	-8(%rsp)	// y
291	fildll	-8(%rsp)	// int(y) : y
292	fucomip	%st(1),%st
293	ffreep	%st		// <empty>
294	jne	26f
295
296	// OK, the value is an integer, but is it odd?
297	mov	-8(%rsp),%eax
298	mov	-4(%rsp),%edx
299	andb	$1, %al
300	jz	27f		// jump if not odd
301	// It's an odd integer.
302	// Raise divide-by-zero exception and get minus infinity value.
303	fldl	MO(one)
304	fdivl	MO(zero)
305	fchs
306	ret
307
30825:	fstp	%st(0)
30926:
31027:	// Raise divide-by-zero exception and get infinity value.
311	fldl	MO(one)
312	fdivl	MO(zero)
313	ret
314
315	.align ALIGNARG(4)
316	// x is �0 and y is > 0.  We must find out whether y is an odd integer.
31721:	testb	$2, %dh
318	jz	22f
319
320	fld	%st		// y : y
321	fistpll	-8(%rsp)	// y
322	fildll	-8(%rsp)	// int(y) : y
323	fucomip %st(1),%st
324	ffreep	%st		// <empty>
325	jne	23f
326
327	// OK, the value is an integer, but is it odd?
328	mov	-8(%rsp),%eax
329	mov	-4(%rsp),%edx
330	andb	$1, %al
331	jz	24f		// jump if not odd
332	// It's an odd integer.
333	fldl	MO(mzero)
334	ret
335
33622:	fstp	%st(0)
33723:
33824:	fldl	MO(zero)
339	ret
340
341END(__ieee754_powl)
342