ix87 specific implementation for complex exponential function for float.
authordrepper <drepper>
Thu, 27 Mar 1997 19:10:16 +0000 (19:10 +0000)
committerdrepper <drepper>
Thu, 27 Mar 1997 19:10:16 +0000 (19:10 +0000)
sysdeps/libm-i387/s_cexpf.S [new file with mode: 0644]

diff --git a/sysdeps/libm-i387/s_cexpf.S b/sysdeps/libm-i387/s_cexpf.S
new file mode 100644 (file)
index 0000000..6fd414b
--- /dev/null
@@ -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)