about summary refs log tree commit diff
path: root/sysdeps/libm-i387
diff options
context:
space:
mode:
Diffstat (limited to 'sysdeps/libm-i387')
-rw-r--r--sysdeps/libm-i387/s_cexp.S248
-rw-r--r--sysdeps/libm-i387/s_cexpf.S245
-rw-r--r--sysdeps/libm-i387/s_cexpl.S249
3 files changed, 742 insertions, 0 deletions
diff --git a/sysdeps/libm-i387/s_cexp.S b/sysdeps/libm-i387/s_cexp.S
new file mode 100644
index 0000000000..48e002b2f6
--- /dev/null
+++ b/sysdeps/libm-i387/s_cexp.S
@@ -0,0 +1,248 @@
+/* ix87 specific implementation of complex exponential function for double.
+   Copyright (C) 1997 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 Library General Public License as
+   published by the Free Software Foundation; either version 2 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C Library; see the file COPYING.LIB.  If not,
+   write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+   Boston, MA 02111-1307, USA.  */
+
+#include <sysdep.h>
+
+#ifdef __ELF__
+	.section .rodata
+#else
+	.text
+#endif
+	.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
+	.double	0.0
+	.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
+        call    1f
+1:      popl    %ecx
+        addl    $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx
+#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	2f
+
+	/* 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
+	andb	$0x01, %ah	/* See above why 0x01 is usable here.  */
+	cmpb	$0x01, %ah
+	je	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> */
+	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)
+	fstpl	8(%eax)
+	fstpl	(%eax)
+	ret	$4
+
+	/* The real part is NaN.  */
+	.align ALIGNARG(4)
+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)
diff --git a/sysdeps/libm-i387/s_cexpf.S b/sysdeps/libm-i387/s_cexpf.S
new file mode 100644
index 0000000000..6fd414b045
--- /dev/null
+++ b/sysdeps/libm-i387/s_cexpf.S
@@ -0,0 +1,245 @@
+/* ix87 specific implementation of complex exponential function for double.
+   Copyright (C) 1997 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 Library General Public License as
+   published by the Free Software Foundation; either version 2 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C Library; see the file COPYING.LIB.  If not,
+   write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+   Boston, MA 02111-1307, USA.  */
+
+#include <sysdep.h>
+
+#ifdef __ELF__
+	.section .rodata
+#else
+	.text
+#endif
+	.align ALIGNARG(4)
+	ASM_TYPE_DIRECTIVE(huge_nan_null_null,@object)
+huge_nan_null_null:
+	.byte 0, 0, 0x80, 0x7f
+	.byte 0, 0, 0xc0, 0x7f
+	.float	0.0
+	.float	0.0
+	.byte 0, 0, 0x80, 0x7f
+	.byte 0, 0, 0xc0, 0x7f
+	.float 0.0
+	.byte 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(__cexpf)
+	flds	4(%esp)			/* x */
+	fxam
+	fnstsw
+	flds	8(%esp)			/* y : x */
+#ifdef  PIC
+        call    1f
+1:      popl    %ecx
+        addl    $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx
+#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	2f
+
+	/* 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) */
+	subl	$8, %esp
+	fstps	4(%esp)
+	fstps	(%esp)
+	popl	%eax
+	popl	%edx
+	ret
+
+	/* 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)
+	subl	$8, %esp
+	fstps	4(%esp)
+	fstps	(%esp)
+	popl	%eax
+	popl	%edx
+	ret
+
+	/* The real part is +-inf.  We must make further differences.  */
+	.align ALIGNARG(4)
+1:	fxam			/* y : x */
+	fnstsw
+	movb	%ah, %dl
+	andb	$0x01, %ah	/* See above why 0x01 is usable here.  */
+	cmpb	$0x01, %ah
+	je	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	$6, %edx
+	fstp	%st(0)		/* y */ /* Drop the real part.  */
+	andl	$8, %edx	/* This puts the sign bit of the real part
+				   in bit 3.  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
+	fxch
+	ftst
+	fnstsw
+	fstp	%st(0)
+	shll	$23, %eax
+	andl	$0x80000000, %eax
+	orl	MOX(huge_nan_null_null,%edx,1), %eax
+	movl	MOX(huge_nan_null_null,%edx,1), %ecx
+	movl	%eax, %edx
+	ftst
+	fnstsw
+	fstp	%st(0)
+	shll	$23, %eax
+	andl	$0x80000000, %eax
+	orl	%ecx, %eax
+	ret
+	/* 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
+	fxch
+	ftst
+	fnstsw
+	fstp	%st(0)
+	shll	$23, %eax
+	andl	$0x80000000, %eax
+	orl	MOX(huge_nan_null_null,%edx,1), %eax
+	movl	MOX(huge_nan_null_null,%edx,1), %ecx
+	movl	%eax, %edx
+	ftst
+	fnstsw
+	fstp	%st(0)
+	shll	$23, %eax
+	andl	$0x80000000, %eax
+	orl	%ecx, %eax
+	ret
+
+	/* The real part is +-Inf and the imaginary part is +-0.  So return
+	   +-Inf+-0i.  */
+	.align ALIGNARG(4)
+4:	subl	$4, %esp
+	fstps	(%esp)
+	shrl	$6, %edx
+	fstp	%st(0)
+	andl	$8, %edx
+	movl	MOX(huge_nan_null_null,%edx,1), %eax
+	popl	%edx
+	ret
+
+	/* The real part is +-Inf, the imaginary is also is not finite.  */
+	.align ALIGNARG(4)
+3:	fstp	%st(0)
+	fstp	%st(0)		/* <empty> */
+	movl	%edx, %eax
+	shrl	$6, %edx
+	shll	$3, %eax
+	andl	$8, %edx
+	andl	$16, %eax
+	orl	%eax, %edx
+
+	movl	MOX(huge_nan_null_null,%edx,1), %eax
+	movl	MOX(huge_nan_null_null+4,%edx,1), %edx
+	ret
+
+	/* The real part is NaN.  */
+	.align ALIGNARG(4)
+2:	fstp	%st(0)
+	fstp	%st(0)
+	movl	MO(huge_nan_null_null+4), %eax
+	movl	%eax, %edx
+	ret
+
+END(__cexpf)
+weak_alias (__cexpf, cexpf)
diff --git a/sysdeps/libm-i387/s_cexpl.S b/sysdeps/libm-i387/s_cexpl.S
new file mode 100644
index 0000000000..fa31e74162
--- /dev/null
+++ b/sysdeps/libm-i387/s_cexpl.S
@@ -0,0 +1,249 @@
+/* ix87 specific implementation of complex exponential function for double.
+   Copyright (C) 1997 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 Library General Public License as
+   published by the Free Software Foundation; either version 2 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
+   Library General Public License for more details.
+
+   You should have received a copy of the GNU Library General Public
+   License along with the GNU C Library; see the file COPYING.LIB.  If not,
+   write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+   Boston, MA 02111-1307, USA.  */
+
+#include <sysdep.h>
+
+#ifdef __ELF__
+	.section .rodata
+#else
+	.text
+#endif
+	.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
+	.double	0.0
+	.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(__cexpl)
+	fldt	8(%esp)			/* x */
+	fxam
+	fnstsw
+	fldt	20(%esp)		/* y : x */
+#ifdef  PIC
+        call    1f
+1:      popl    %ecx
+        addl    $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx
+#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	2f
+
+	/* 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.  */
+	fstpt	12(%eax)
+	fstpt	(%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.  */
+	fstpt	12(%eax)
+	fstpt	(%eax)
+	ret	$4
+
+	/* The real part is +-inf.  We must make further differences.  */
+	.align ALIGNARG(4)
+1:	fxam			/* y : x */
+	fnstsw
+	movb	%ah, %dl
+	andb	$0x01, %ah	/* See above why 0x01 is usable here.  */
+	cmpb	$0x01, %ah
+	je	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	$7, %eax
+	andl	$0x8000, %eax
+	orl	%eax, 8(%edx)
+	fstp	%st(0)
+	ftst
+	fnstsw
+	shll	$7, %eax
+	andl	$0x8000, %eax
+	orl	%eax, 20(%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	$7, %eax
+	andl	$0x8000, %eax
+	orl	%eax, 8(%edx)
+	fstp	%st(0)
+	ftst
+	fnstsw
+	shll	$7, %eax
+	andl	$0x8000, %eax
+	orl	%eax, 20(%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.  */
+	fstpt	12(%eax)
+	shrl	$5, %edx
+	fstp	%st(0)
+	andl	$16, %edx
+	fldl	MOX(huge_nan_null_null,%edx,1)
+	fstpt	(%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> */
+	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)
+	fstpt	12(%eax)
+	fstpt	(%eax)
+	ret	$4
+
+	/* The real part is NaN.  */
+	.align ALIGNARG(4)
+2:	fstp	%st(0)
+	fstp	%st(0)
+	movl	4(%esp), %eax		/* Pointer to memory for result.  */
+	fldl	MO(huge_nan_null_null+8)
+	fld	%st(0)
+	fstpt	(%eax)
+	fstpt	12(%eax)
+	ret	$4
+
+END(__cexpl)
+weak_alias (__cexpl, cexpl)