[thirdparty/glibc.git] / sysdeps / i386 / fpu / s_cexp.S

/* ix87 specific implementation of complex exponential function for double.
   Copyright (C) 1997, 2005, 2012 Free Software Foundation, Inc.
   This file is part of the GNU C Library.
   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.

   The GNU C Library is free software; you can redistribute it and/or
   modify it under the terms of the GNU Lesser General Public
   License as published by the Free Software Foundation; either
   version 2.1 of the License, or (at your option) any later version.

   The GNU C Library is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   Lesser General Public License for more details.

   You should have received a copy of the GNU Lesser General Public
   License along with the GNU C Library; if not, see
   <http://www.gnu.org/licenses/>.  */

#include <sysdep.h>

	.section .rodata

	.align ALIGNARG(4)
	ASM_TYPE_DIRECTIVE(huge_nan_null_null,@object)
huge_nan_null_null:
	.byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
	.byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f
	.double	0.0
zero:	.double	0.0
infinity:
	.byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
	.byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f
	.double 0.0
	.byte 0, 0, 0, 0, 0, 0, 0, 0x80
	ASM_SIZE_DIRECTIVE(huge_nan_null_null)

	ASM_TYPE_DIRECTIVE(twopi,@object)
twopi:
	.byte 0x35, 0xc2, 0x68, 0x21, 0xa2, 0xda, 0xf, 0xc9, 0x1, 0x40
	.byte 0, 0, 0, 0, 0, 0
	ASM_SIZE_DIRECTIVE(twopi)

	ASM_TYPE_DIRECTIVE(l2e,@object)
l2e:
	.byte 0xbc, 0xf0, 0x17, 0x5c, 0x29, 0x3b, 0xaa, 0xb8, 0xff, 0x3f
	.byte 0, 0, 0, 0, 0, 0
	ASM_SIZE_DIRECTIVE(l2e)

	ASM_TYPE_DIRECTIVE(one,@object)
one:	.double 1.0
	ASM_SIZE_DIRECTIVE(one)


#ifdef PIC
#define MO(op) op##@GOTOFF(%ecx)
#define MOX(op,x,f) op##@GOTOFF(%ecx,x,f)
#else
#define MO(op) op
#define MOX(op,x,f) op(,x,f)
#endif

	.text
ENTRY(__cexp)
	fldl	8(%esp)			/* x */
	fxam
	fnstsw
	fldl	16(%esp)		/* y : x */
#ifdef  PIC
	LOAD_PIC_REG (cx)
#endif
	movb	%ah, %dh
	andb	$0x45, %ah
	cmpb	$0x05, %ah
	je	1f			/* Jump if real part is +-Inf */
	cmpb	$0x01, %ah
	je	2f			/* Jump if real part is NaN */

	fxam				/* y : x */
	fnstsw
	/* If the imaginary part is not finite we return NaN+i NaN, as
	   for the case when the real part is NaN.  A test for +-Inf and
	   NaN would be necessary.  But since we know the stack register
	   we applied `fxam' to is not empty we can simply use one test.
	   Check your FPU manual for more information.  */
	andb	$0x01, %ah
	cmpb	$0x01, %ah
	je	20f

	/* We have finite numbers in the real and imaginary part.  Do
	   the real work now.  */
	fxch			/* x : y */
	fldt	MO(l2e)		/* log2(e) : x : y */
	fmulp			/* x * log2(e) : y */
	fld	%st		/* x * log2(e) : x * log2(e) : y */
	frndint			/* int(x * log2(e)) : x * log2(e) : y */
	fsubr	%st, %st(1)	/* int(x * log2(e)) : frac(x * log2(e)) : y */
	fxch			/* frac(x * log2(e)) : int(x * log2(e)) : y */
	f2xm1			/* 2^frac(x * log2(e))-1 : int(x * log2(e)) : y */
	faddl	MO(one)		/* 2^frac(x * log2(e)) : int(x * log2(e)) : y */
	fscale			/* e^x : int(x * log2(e)) : y */
	fst	%st(1)		/* e^x : e^x : y */
	fxch	%st(2)		/* y : e^x : e^x */
	fsincos			/* cos(y) : sin(y) : e^x : e^x */
	fnstsw
	testl	$0x400, %eax
	jnz	7f
	fmulp	%st, %st(3)	/* sin(y) : e^x : e^x * cos(y) */
	fmulp	%st, %st(1)	/* e^x * sin(y) : e^x * cos(y) */
	movl	4(%esp), %eax		/* Pointer to memory for result.  */
	fstpl	8(%eax)
	fstpl	(%eax)
	ret	$4

	/* We have to reduce the argument to fsincos.  */
	.align ALIGNARG(4)
7:	fldt	MO(twopi)	/* 2*pi : y : e^x : e^x */
	fxch			/* y : 2*pi : e^x : e^x */
8:	fprem1			/* y%(2*pi) : 2*pi : e^x : e^x */
	fnstsw
	testl	$0x400, %eax
	jnz	8b
	fstp	%st(1)		/* y%(2*pi) : e^x : e^x */
	fsincos			/* cos(y) : sin(y) : e^x : e^x */
	fmulp	%st, %st(3)
	fmulp	%st, %st(1)
	movl	4(%esp), %eax		/* Pointer to memory for result.  */
	fstpl	8(%eax)
	fstpl	(%eax)
	ret	$4

	/* The real part is +-inf.  We must make further differences.  */
	.align ALIGNARG(4)
1:	fxam			/* y : x */
	fnstsw
	movb	%ah, %dl
	testb	$0x01, %ah	/* See above why 0x01 is usable here.  */
	jne	3f


	/* The real part is +-Inf and the imaginary part is finite.  */
	andl	$0x245, %edx
	cmpb	$0x40, %dl	/* Imaginary part == 0?  */
	je	4f		/* Yes ->  */

	fxch			/* x : y */
	shrl	$5, %edx
	fstp	%st(0)		/* y */ /* Drop the real part.  */
	andl	$16, %edx	/* This puts the sign bit of the real part
				   in bit 4.  So we can use it to index a
				   small array to select 0 or Inf.  */
	fsincos			/* cos(y) : sin(y) */
	fnstsw
	testl	$0x0400, %eax
	jnz	5f
	fldl	MOX(huge_nan_null_null,%edx,1)
	movl	4(%esp), %edx		/* Pointer to memory for result.  */
	fstl	8(%edx)
	fstpl	(%edx)
	ftst
	fnstsw
	shll	$23, %eax
	andl	$0x80000000, %eax
	orl	%eax, 4(%edx)
	fstp	%st(0)
	ftst
	fnstsw
	shll	$23, %eax
	andl	$0x80000000, %eax
	orl	%eax, 12(%edx)
	fstp	%st(0)
	ret	$4
	/* We must reduce the argument to fsincos.  */
	.align ALIGNARG(4)
5:	fldt	MO(twopi)
	fxch
6:	fprem1
	fnstsw
	testl	$0x400, %eax
	jnz	6b
	fstp	%st(1)
	fsincos
	fldl	MOX(huge_nan_null_null,%edx,1)
	movl	4(%esp), %edx		/* Pointer to memory for result.  */
	fstl	8(%edx)
	fstpl	(%edx)
	ftst
	fnstsw
	shll	$23, %eax
	andl	$0x80000000, %eax
	orl	%eax, 4(%edx)
	fstp	%st(0)
	ftst
	fnstsw
	shll	$23, %eax
	andl	$0x80000000, %eax
	orl	%eax, 12(%edx)
	fstp	%st(0)
	ret	$4

	/* The real part is +-Inf and the imaginary part is +-0.  So return
	   +-Inf+-0i.  */
	.align ALIGNARG(4)
4:	movl	4(%esp), %eax		/* Pointer to memory for result.  */
	fstpl	8(%eax)
	shrl	$5, %edx
	fstp	%st(0)
	andl	$16, %edx
	fldl	MOX(huge_nan_null_null,%edx,1)
	fstpl	(%eax)
	ret	$4

	/* The real part is +-Inf, the imaginary is also is not finite.  */
	.align ALIGNARG(4)
3:	fstp	%st(0)
	fstp	%st(0)		/* <empty> */
	andb	$0x45, %ah
	andb	$0x47, %dh
	xorb	%dh, %ah
	jnz	30f
	fldl	MO(infinity)	/* Raise invalid exception.  */
	fmull	MO(zero)
	fstp	%st(0)
30:	movl	%edx, %eax
	shrl	$5, %edx
	shll	$4, %eax
	andl	$16, %edx
	andl	$32, %eax
	orl	%eax, %edx
	movl	4(%esp), %eax		/* Pointer to memory for result.  */

	fldl	MOX(huge_nan_null_null,%edx,1)
	fldl	MOX(huge_nan_null_null+8,%edx,1)
	fxch
	fstpl	(%eax)
	fstpl	8(%eax)
	ret	$4

	/* The real part is NaN.  */
	.align ALIGNARG(4)
20:	fldl	MO(infinity)		/* Raise invalid exception.  */
	fmull	MO(zero)
	fstp	%st(0)
2:	fstp	%st(0)
	fstp	%st(0)
	movl	4(%esp), %eax		/* Pointer to memory for result.  */
	fldl	MO(huge_nan_null_null+8)
	fstl	(%eax)
	fstpl	8(%eax)
	ret	$4

END(__cexp)
weak_alias (__cexp, cexp)
Commit	Line	Data
993b3242	1	/* ix87 specific implementation of complex exponential function for double.
622c86f4	2	Copyright (C) 1997, 2005, 2012 Free Software Foundation, Inc.
993b3242 UD	3	This file is part of the GNU C Library.
	4	Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
	5
	6	The GNU C Library is free software; you can redistribute it and/or
41bdb6e2 AJ	7	modify it under the terms of the GNU Lesser General Public
	8	License as published by the Free Software Foundation; either
	9	version 2.1 of the License, or (at your option) any later version.
993b3242 UD	10
	11	The GNU C Library is distributed in the hope that it will be useful,
	12	but WITHOUT ANY WARRANTY; without even the implied warranty of
	13	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
41bdb6e2	14	Lesser General Public License for more details.
993b3242	15
41bdb6e2	16	You should have received a copy of the GNU Lesser General Public
59ba27a6 PE	17	License along with the GNU C Library; if not, see
59ba27a6 PE	18	<http://www.gnu.org/licenses/>. */
993b3242 UD	19
	20	#include <sysdep.h>
	21
993b3242	22	.section .rodata
622c86f4	23
993b3242 UD	24	.align ALIGNARG(4)
	25	ASM_TYPE_DIRECTIVE(huge_nan_null_null,@object)
	26	huge_nan_null_null:
	27	.byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
	28	.byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f
	29	.double 0.0
779ae82e UD	30	zero: .double 0.0
779ae82e UD	31	infinity:
993b3242 UD	32	.byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
	33	.byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f
	34	.double 0.0
	35	.byte 0, 0, 0, 0, 0, 0, 0, 0x80
	36	ASM_SIZE_DIRECTIVE(huge_nan_null_null)
	37
	38	ASM_TYPE_DIRECTIVE(twopi,@object)
	39	twopi:
	40	.byte 0x35, 0xc2, 0x68, 0x21, 0xa2, 0xda, 0xf, 0xc9, 0x1, 0x40
	41	.byte 0, 0, 0, 0, 0, 0
	42	ASM_SIZE_DIRECTIVE(twopi)
	43
	44	ASM_TYPE_DIRECTIVE(l2e,@object)
	45	l2e:
	46	.byte 0xbc, 0xf0, 0x17, 0x5c, 0x29, 0x3b, 0xaa, 0xb8, 0xff, 0x3f
	47	.byte 0, 0, 0, 0, 0, 0
	48	ASM_SIZE_DIRECTIVE(l2e)
	49
	50	ASM_TYPE_DIRECTIVE(one,@object)
	51	one: .double 1.0
	52	ASM_SIZE_DIRECTIVE(one)
	53
	54
	55	#ifdef PIC
	56	#define MO(op) op##@GOTOFF(%ecx)
	57	#define MOX(op,x,f) op##@GOTOFF(%ecx,x,f)
	58	#else
	59	#define MO(op) op
	60	#define MOX(op,x,f) op(,x,f)
	61	#endif
	62
	63	.text
	64	ENTRY(__cexp)
	65	fldl 8(%esp) /* x */
	66	fxam
	67	fnstsw
	68	fldl 16(%esp) /* y : x */
	69	#ifdef PIC
fee732e5	70	LOAD_PIC_REG (cx)
993b3242 UD	71	#endif
	72	movb %ah, %dh
	73	andb $0x45, %ah
	74	cmpb $0x05, %ah
	75	je 1f /* Jump if real part is +-Inf */
	76	cmpb $0x01, %ah
	77	je 2f /* Jump if real part is NaN */
	78
	79	fxam /* y : x */
	80	fnstsw
	81	/* If the imaginary part is not finite we return NaN+i NaN, as
	82	for the case when the real part is NaN. A test for +-Inf and
	83	NaN would be necessary. But since we know the stack register
	84	we applied `fxam' to is not empty we can simply use one test.
	85	Check your FPU manual for more information. */
	86	andb $0x01, %ah
	87	cmpb $0x01, %ah
779ae82e	88	je 20f
993b3242 UD	89
	90	/* We have finite numbers in the real and imaginary part. Do
	91	the real work now. */
	92	fxch /* x : y */
	93	fldt MO(l2e) /* log2(e) : x : y */
	94	fmulp /* x * log2(e) : y */
	95	fld %st /* x * log2(e) : x * log2(e) : y */
	96	frndint /* int(x * log2(e)) : x * log2(e) : y */
	97	fsubr %st, %st(1) /* int(x * log2(e)) : frac(x * log2(e)) : y */
	98	fxch /* frac(x * log2(e)) : int(x * log2(e)) : y */
	99	f2xm1 /* 2^frac(x * log2(e))-1 : int(x * log2(e)) : y */
	100	faddl MO(one) /* 2^frac(x * log2(e)) : int(x * log2(e)) : y */
	101	fscale /* e^x : int(x * log2(e)) : y */
	102	fst %st(1) /* e^x : e^x : y */
	103	fxch %st(2) /* y : e^x : e^x */
	104	fsincos /* cos(y) : sin(y) : e^x : e^x */
	105	fnstsw
	106	testl $0x400, %eax
	107	jnz 7f
	108	fmulp %st, %st(3) /* sin(y) : e^x : e^x * cos(y) */
	109	fmulp %st, %st(1) /* e^x * sin(y) : e^x * cos(y) */
	110	movl 4(%esp), %eax /* Pointer to memory for result. */
	111	fstpl 8(%eax)
	112	fstpl (%eax)
	113	ret $4
	114
	115	/* We have to reduce the argument to fsincos. */
	116	.align ALIGNARG(4)
	117	7: fldt MO(twopi) /* 2pi : y : e^x : e^x /
	118	fxch /* y : 2pi : e^x : e^x /
	119	8: fprem1 /* y%(2pi) : 2pi : e^x : e^x */
	120	fnstsw
	121	testl $0x400, %eax
	122	jnz 8b
	123	fstp %st(1) /* y%(2pi) : e^x : e^x /
	124	fsincos /* cos(y) : sin(y) : e^x : e^x */
	125	fmulp %st, %st(3)
	126	fmulp %st, %st(1)
	127	movl 4(%esp), %eax /* Pointer to memory for result. */
	128	fstpl 8(%eax)
	129	fstpl (%eax)
	130	ret $4
	131
	132	/* The real part is +-inf. We must make further differences. */
	133	.align ALIGNARG(4)
	134	1: fxam /* y : x */
	135	fnstsw
	136	movb %ah, %dl
779ae82e UD	137	testb $0x01, %ah /* See above why 0x01 is usable here. */
779ae82e UD	138	jne 3f
993b3242 UD	139
	140
	141	/* The real part is +-Inf and the imaginary part is finite. */
	142	andl $0x245, %edx
	143	cmpb $0x40, %dl /* Imaginary part == 0? */
	144	je 4f /* Yes -> */
	145
	146	fxch /* x : y */
	147	shrl $5, %edx
	148	fstp %st(0) /* y / / Drop the real part. */
	149	andl $16, %edx /* This puts the sign bit of the real part
	150	in bit 4. So we can use it to index a
	151	small array to select 0 or Inf. */
	152	fsincos /* cos(y) : sin(y) */
	153	fnstsw
	154	testl $0x0400, %eax
	155	jnz 5f
	156	fldl MOX(huge_nan_null_null,%edx,1)
	157	movl 4(%esp), %edx /* Pointer to memory for result. */
	158	fstl 8(%edx)
	159	fstpl (%edx)
	160	ftst
	161	fnstsw
	162	shll $23, %eax
	163	andl $0x80000000, %eax
	164	orl %eax, 4(%edx)
	165	fstp %st(0)
	166	ftst
	167	fnstsw
	168	shll $23, %eax
	169	andl $0x80000000, %eax
	170	orl %eax, 12(%edx)
	171	fstp %st(0)
	172	ret $4
	173	/* We must reduce the argument to fsincos. */
	174	.align ALIGNARG(4)
	175	5: fldt MO(twopi)
	176	fxch
	177	6: fprem1
	178	fnstsw
	179	testl $0x400, %eax
	180	jnz 6b
	181	fstp %st(1)
	182	fsincos
	183	fldl MOX(huge_nan_null_null,%edx,1)
	184	movl 4(%esp), %edx /* Pointer to memory for result. */
	185	fstl 8(%edx)
	186	fstpl (%edx)
	187	ftst
	188	fnstsw
	189	shll $23, %eax
	190	andl $0x80000000, %eax
	191	orl %eax, 4(%edx)
	192	fstp %st(0)
	193	ftst
	194	fnstsw
	195	shll $23, %eax
	196	andl $0x80000000, %eax
	197	orl %eax, 12(%edx)
	198	fstp %st(0)
	199	ret $4
	200
	201	/* The real part is +-Inf and the imaginary part is +-0. So return
	202	+-Inf+-0i. */
203	.align ALIGNARG(4)
204	4: movl 4(%esp), %eax /* Pointer to memory for result. */
205	fstpl 8(%eax)
206	shrl $5, %edx
207	fstp %st(0)
208	andl $16, %edx
209	fldl MOX(huge_nan_null_null,%edx,1)
210	fstpl (%eax)
211	ret $4
212
213	/* The real part is +-Inf, the imaginary is also is not finite. */
214	.align ALIGNARG(4)
215	3: fstp %st(0)
216	fstp %st(0) /* <empty> */
779ae82e UD	217	andb $0x45, %ah
	218	andb $0x47, %dh
	219	xorb %dh, %ah
	220	jnz 30f
	221	fldl MO(infinity) /* Raise invalid exception. */
	222	fmull MO(zero)
	223	fstp %st(0)
	224	30: movl %edx, %eax
993b3242 UD	225	shrl $5, %edx
	226	shll $4, %eax
	227	andl $16, %edx
	228	andl $32, %eax
	229	orl %eax, %edx
	230	movl 4(%esp), %eax /* Pointer to memory for result. */
	231
	232	fldl MOX(huge_nan_null_null,%edx,1)
	233	fldl MOX(huge_nan_null_null+8,%edx,1)
40a55d20	234	fxch
993b3242	235	fstpl (%eax)
40a55d20	236	fstpl 8(%eax)
993b3242 UD	237	ret $4
	238
	239	/* The real part is NaN. */
	240	.align ALIGNARG(4)
779ae82e UD	241	20: fldl MO(infinity) /* Raise invalid exception. */
	242	fmull MO(zero)
	243	fstp %st(0)
993b3242 UD	244	2: fstp %st(0)
	245	fstp %st(0)
	246	movl 4(%esp), %eax /* Pointer to memory for result. */
	247	fldl MO(huge_nan_null_null+8)
	248	fstl (%eax)
	249	fstpl 8(%eax)
	250	ret $4
	251
	252	END(__cexp)
	253	weak_alias (__cexp, cexp)