1/* ix87 specific implementation of pow function.
2   Copyright (C) 1996-2022 Free Software Foundation, Inc.
3   This file is part of the GNU C Library.
4
5   The GNU C Library is free software; you can redistribute it and/or
6   modify it under the terms of the GNU Lesser General Public
7   License as published by the Free Software Foundation; either
8   version 2.1 of the License, or (at your option) any later version.
9
10   The GNU C Library is distributed in the hope that it will be useful,
11   but WITHOUT ANY WARRANTY; without even the implied warranty of
12   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13   Lesser General Public License for more details.
14
15   You should have received a copy of the GNU Lesser General Public
16   License along with the GNU C Library; if not, see
17   <https://www.gnu.org/licenses/>.  */
18
19#include <machine/asm.h>
20#include <i386-math-asm.h>
21#include <libm-alias-finite.h>
22
23	.section .rodata.cst8,"aM",@progbits,8
24
25	.p2align 3
26	.type one,@object
27one:	.double 1.0
28	ASM_SIZE_DIRECTIVE(one)
29	.type limit,@object
30limit:	.double 0.29
31	ASM_SIZE_DIRECTIVE(limit)
32	.type p63,@object
33p63:	.byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43
34	ASM_SIZE_DIRECTIVE(p63)
35	.type p10,@object
36p10:	.byte 0, 0, 0, 0, 0, 0, 0x90, 0x40
37	ASM_SIZE_DIRECTIVE(p10)
38
39	.section .rodata.cst16,"aM",@progbits,16
40
41	.p2align 3
42	.type infinity,@object
43inf_zero:
44infinity:
45	.byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
46	ASM_SIZE_DIRECTIVE(infinity)
47	.type zero,@object
48zero:	.double 0.0
49	ASM_SIZE_DIRECTIVE(zero)
50	.type minf_mzero,@object
51minf_mzero:
52minfinity:
53	.byte 0, 0, 0, 0, 0, 0, 0xf0, 0xff
54mzero:
55	.byte 0, 0, 0, 0, 0, 0, 0, 0x80
56	ASM_SIZE_DIRECTIVE(minf_mzero)
57DEFINE_DBL_MIN
58
59#ifdef PIC
60# define MO(op) op##@GOTOFF(%ecx)
61# define MOX(op,x,f) op##@GOTOFF(%ecx,x,f)
62#else
63# define MO(op) op
64# define MOX(op,x,f) op(,x,f)
65#endif
66
67	.text
68ENTRY(__ieee754_pow)
69	fldl	12(%esp)	// y
70	fxam
71
72#ifdef	PIC
73	LOAD_PIC_REG (cx)
74#endif
75
76	fnstsw
77	movb	%ah, %dl
78	andb	$0x45, %ah
79	cmpb	$0x40, %ah	// is y == 0 ?
80	je	11f
81
82	cmpb	$0x05, %ah	// is y == ±inf ?
83	je	12f
84
85	cmpb	$0x01, %ah	// is y == NaN ?
86	je	30f
87
88	fldl	4(%esp)		// x : y
89
90	subl	$8,%esp
91	cfi_adjust_cfa_offset (8)
92
93	fxam
94	fnstsw
95	movb	%ah, %dh
96	andb	$0x45, %ah
97	cmpb	$0x40, %ah
98	je	20f		// x is ±0
99
100	cmpb	$0x05, %ah
101	je	15f		// x is ±inf
102
103	cmpb	$0x01, %ah
104	je	32f		// x is NaN
105
106	fxch			// y : x
107
108	/* fistpll raises invalid exception for |y| >= 1L<<63.  */
109	fld	%st		// y : y : x
110	fabs			// |y| : y : x
111	fcompl	MO(p63)		// y : x
112	fnstsw
113	sahf
114	jnc	2f
115
116	/* First see whether `y' is a natural number.  In this case we
117	   can use a more precise algorithm.  */
118	fld	%st		// y : y : x
119	fistpll	(%esp)		// y : x
120	fildll	(%esp)		// int(y) : y : x
121	fucomp	%st(1)		// y : x
122	fnstsw
123	sahf
124	jne	3f
125
126	/* OK, we have an integer value for y.  If large enough that
127	   errors may propagate out of the 11 bits excess precision, use
128	   the algorithm for real exponent instead.  */
129	fld	%st		// y : y : x
130	fabs			// |y| : y : x
131	fcompl	MO(p10)		// y : x
132	fnstsw
133	sahf
134	jnc	2f
135	popl	%eax
136	cfi_adjust_cfa_offset (-4)
137	popl	%edx
138	cfi_adjust_cfa_offset (-4)
139	orl	$0, %edx
140	fstp	%st(0)		// x
141	jns	4f		// y >= 0, jump
142	fdivrl	MO(one)		// 1/x		(now referred to as x)
143	negl	%eax
144	adcl	$0, %edx
145	negl	%edx
1464:	fldl	MO(one)		// 1 : x
147	fxch
148
149	/* If y is even, take the absolute value of x.  Otherwise,
150	   ensure all intermediate values that might overflow have the
151	   sign of x.  */
152	testb	$1, %al
153	jnz	6f
154	fabs
155
1566:	shrdl	$1, %edx, %eax
157	jnc	5f
158	fxch
159	fabs
160	fmul	%st(1)		// x : ST*x
161	fxch
1625:	fld	%st		// x : x : ST*x
163	fabs			// |x| : x : ST*x
164	fmulp			// |x|*x : ST*x
165	shrl	$1, %edx
166	movl	%eax, %ecx
167	orl	%edx, %ecx
168	jnz	6b
169	fstp	%st(0)		// ST*x
170#ifdef	PIC
171	LOAD_PIC_REG (cx)
172#endif
173	DBL_NARROW_EVAL_UFLOW_NONNAN
174	ret
175
176	/* y is ±NAN */
17730:	fldl	4(%esp)		// x : y
178	fldl	MO(one)		// 1.0 : x : y
179	fucomp	%st(1)		// x : y
180	fnstsw
181	sahf
182	je	31f
183	fxch			// y : x
18431:	fstp	%st(1)
185	ret
186
187	cfi_adjust_cfa_offset (8)
18832:	addl	$8, %esp
189	cfi_adjust_cfa_offset (-8)
190	fstp	%st(1)
191	ret
192
193	cfi_adjust_cfa_offset (8)
194	.align ALIGNARG(4)
1952:	// y is a large integer (absolute value at least 1L<<10), but
196	// may be odd unless at least 1L<<64.  So it may be necessary
197	// to adjust the sign of a negative result afterwards.
198	fxch			// x : y
199	fabs			// |x| : y
200	fxch			// y : x
201	.align ALIGNARG(4)
2023:	/* y is a real number.  */
203	fxch			// x : y
204	fldl	MO(one)		// 1.0 : x : y
205	fldl	MO(limit)	// 0.29 : 1.0 : x : y
206	fld	%st(2)		// x : 0.29 : 1.0 : x : y
207	fsub	%st(2)		// x-1 : 0.29 : 1.0 : x : y
208	fabs			// |x-1| : 0.29 : 1.0 : x : y
209	fucompp			// 1.0 : x : y
210	fnstsw
211	fxch			// x : 1.0 : y
212	sahf
213	ja	7f
214	fsub	%st(1)		// x-1 : 1.0 : y
215	fyl2xp1			// log2(x) : y
216	jmp	8f
217
2187:	fyl2x			// log2(x) : y
2198:	fmul	%st(1)		// y*log2(x) : y
220	fst	%st(1)		// y*log2(x) : y*log2(x)
221	frndint			// int(y*log2(x)) : y*log2(x)
222	fsubr	%st, %st(1)	// int(y*log2(x)) : fract(y*log2(x))
223	fxch			// fract(y*log2(x)) : int(y*log2(x))
224	f2xm1			// 2^fract(y*log2(x))-1 : int(y*log2(x))
225	faddl	MO(one)		// 2^fract(y*log2(x)) : int(y*log2(x))
226
227	// Before scaling, we must negate if x is negative and y is an
228	// odd integer.
229	testb	$2, %dh
230	jz	291f
231	// x is negative.  If y is an odd integer, negate the result.
232	fldl	20(%esp)	// y : 2^fract(y*log2(x)) : int(y*log2(x))
233	fld	%st		// y : y : 2^fract(y*log2(x)) : int(y*log2(x))
234	fabs			// |y| : y : 2^fract(y*log2(x)) : int(y*log2(x))
235	fcompl	MO(p63)		// y : 2^fract(y*log2(x)) : int(y*log2(x))
236	fnstsw
237	sahf
238	jnc	290f
239
240	// We must find out whether y is an odd integer.
241	fld	%st		// y : y : 2^fract(y*log2(x)) : int(y*log2(x))
242	fistpll	(%esp)		// y : 2^fract(y*log2(x)) : int(y*log2(x))
243	fildll	(%esp)		// int(y) : y : 2^fract(y*log2(x)) : int(y*log2(x))
244	fucompp			// 2^fract(y*log2(x)) : int(y*log2(x))
245	fnstsw
246	sahf
247	jne	291f
248
249	// OK, the value is an integer, but is it odd?
250	popl	%eax
251	cfi_adjust_cfa_offset (-4)
252	popl	%edx
253	cfi_adjust_cfa_offset (-4)
254	andb	$1, %al
255	jz	292f		// jump if not odd
256	// It's an odd integer.
257	fchs
258	jmp	292f
259
260	cfi_adjust_cfa_offset (8)
261290:	fstp	%st(0)		// 2^fract(y*log2(x)) : int(y*log2(x))
262291:	addl	$8, %esp
263	cfi_adjust_cfa_offset (-8)
264292:	fscale			// +/- 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
265	fstp	%st(1)		// +/- 2^fract(y*log2(x))*2^int(y*log2(x))
266	DBL_NARROW_EVAL_UFLOW_NONNAN
267	ret
268
269
270	// pow(x,±0) = 1
271	.align ALIGNARG(4)
27211:	fstp	%st(0)		// pop y
273	fldl	MO(one)
274	ret
275
276	// y == ±inf
277	.align ALIGNARG(4)
27812:	fstp	%st(0)		// pop y
279	fldl	MO(one)		// 1
280	fldl	4(%esp)		// x : 1
281	fabs			// abs(x) : 1
282	fucompp			// < 1, == 1, or > 1
283	fnstsw
284	andb	$0x45, %ah
285	cmpb	$0x45, %ah
286	je	13f		// jump if x is NaN
287
288	cmpb	$0x40, %ah
289	je	14f		// jump if |x| == 1
290
291	shlb	$1, %ah
292	xorb	%ah, %dl
293	andl	$2, %edx
294	fldl	MOX(inf_zero, %edx, 4)
295	ret
296
297	.align ALIGNARG(4)
29814:	fldl	MO(one)
299	ret
300
301	.align ALIGNARG(4)
30213:	fldl	4(%esp)		// load x == NaN
303	ret
304
305	cfi_adjust_cfa_offset (8)
306	.align ALIGNARG(4)
307	// x is ±inf
30815:	fstp	%st(0)		// y
309	testb	$2, %dh
310	jz	16f		// jump if x == +inf
311
312	// fistpll raises invalid exception for |y| >= 1L<<63, so test
313	// that (in which case y is certainly even) before testing
314	// whether y is odd.
315	fld	%st		// y : y
316	fabs			// |y| : y
317	fcompl	MO(p63)		// y
318	fnstsw
319	sahf
320	jnc	16f
321
322	// We must find out whether y is an odd integer.
323	fld	%st		// y : y
324	fistpll	(%esp)		// y
325	fildll	(%esp)		// int(y) : y
326	fucompp			// <empty>
327	fnstsw
328	sahf
329	jne	17f
330
331	// OK, the value is an integer.
332	popl	%eax
333	cfi_adjust_cfa_offset (-4)
334	popl	%edx
335	cfi_adjust_cfa_offset (-4)
336	andb	$1, %al
337	jz	18f		// jump if not odd
338	// It's an odd integer.
339	shrl	$31, %edx
340	fldl	MOX(minf_mzero, %edx, 8)
341	ret
342
343	cfi_adjust_cfa_offset (8)
344	.align ALIGNARG(4)
34516:	fcompl	MO(zero)
346	addl	$8, %esp
347	cfi_adjust_cfa_offset (-8)
348	fnstsw
349	shrl	$5, %eax
350	andl	$8, %eax
351	fldl	MOX(inf_zero, %eax, 1)
352	ret
353
354	cfi_adjust_cfa_offset (8)
355	.align ALIGNARG(4)
35617:	shll	$30, %edx	// sign bit for y in right position
357	addl	$8, %esp
358	cfi_adjust_cfa_offset (-8)
35918:	shrl	$31, %edx
360	fldl	MOX(inf_zero, %edx, 8)
361	ret
362
363	cfi_adjust_cfa_offset (8)
364	.align ALIGNARG(4)
365	// x is ±0
36620:	fstp	%st(0)		// y
367	testb	$2, %dl
368	jz	21f		// y > 0
369
370	// x is ±0 and y is < 0.  We must find out whether y is an odd integer.
371	testb	$2, %dh
372	jz	25f
373
374	// fistpll raises invalid exception for |y| >= 1L<<63, so test
375	// that (in which case y is certainly even) before testing
376	// whether y is odd.
377	fld	%st		// y : y
378	fabs			// |y| : y
379	fcompl	MO(p63)		// y
380	fnstsw
381	sahf
382	jnc	25f
383
384	fld	%st		// y : y
385	fistpll	(%esp)		// y
386	fildll	(%esp)		// int(y) : y
387	fucompp			// <empty>
388	fnstsw
389	sahf
390	jne	26f
391
392	// OK, the value is an integer.
393	popl	%eax
394	cfi_adjust_cfa_offset (-4)
395	popl	%edx
396	cfi_adjust_cfa_offset (-4)
397	andb	$1, %al
398	jz	27f		// jump if not odd
399	// It's an odd integer.
400	// Raise divide-by-zero exception and get minus infinity value.
401	fldl	MO(one)
402	fdivl	MO(zero)
403	fchs
404	ret
405
406	cfi_adjust_cfa_offset (8)
40725:	fstp	%st(0)
40826:	addl	$8, %esp
409	cfi_adjust_cfa_offset (-8)
41027:	// Raise divide-by-zero exception and get infinity value.
411	fldl	MO(one)
412	fdivl	MO(zero)
413	ret
414
415	cfi_adjust_cfa_offset (8)
416	.align ALIGNARG(4)
417	// x is ±0 and y is > 0.  We must find out whether y is an odd integer.
41821:	testb	$2, %dh
419	jz	22f
420
421	// fistpll raises invalid exception for |y| >= 1L<<63, so test
422	// that (in which case y is certainly even) before testing
423	// whether y is odd.
424	fcoml	MO(p63)		// y
425	fnstsw
426	sahf
427	jnc	22f
428
429	fld	%st		// y : y
430	fistpll	(%esp)		// y
431	fildll	(%esp)		// int(y) : y
432	fucompp			// <empty>
433	fnstsw
434	sahf
435	jne	23f
436
437	// OK, the value is an integer.
438	popl	%eax
439	cfi_adjust_cfa_offset (-4)
440	popl	%edx
441	cfi_adjust_cfa_offset (-4)
442	andb	$1, %al
443	jz	24f		// jump if not odd
444	// It's an odd integer.
445	fldl	MO(mzero)
446	ret
447
448	cfi_adjust_cfa_offset (8)
44922:	fstp	%st(0)
45023:	addl	$8, %esp	// Don't use 2 x pop
451	cfi_adjust_cfa_offset (-8)
45224:	fldl	MO(zero)
453	ret
454
455END(__ieee754_pow)
456libm_alias_finite (__ieee754_pow, __pow)
457