sysdeps/i386/fpu/s_cexpf.S

   1 /* ix87 specific implementation of complex exponential function for double.
   2    Copyright (C) 1997, 2005 Free Software Foundation, Inc.
   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
   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.
  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
  14    Lesser General Public License for more details.
  15
  16    You should have received a copy of the GNU Lesser General Public
  17    License along with the GNU C Library; if not, write to the Free
  18    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
  19    02111-1307 USA.  */
  20
  21 #include <sysdep.h>
  22
  23 #ifdef __ELF__
  24         .section .rodata
  25 #else
  26         .text
  27 #endif
  28         .align ALIGNARG(4)
  29         ASM_TYPE_DIRECTIVE(huge_nan_null_null,@object)
  30 huge_nan_null_null:
  31         .byte 0, 0, 0x80, 0x7f
  32         .byte 0, 0, 0xc0, 0x7f
  33         .float  0.0
  34 zero:   .float  0.0
  35 infinity:
  36         .byte 0, 0, 0x80, 0x7f
  37         .byte 0, 0, 0xc0, 0x7f
  38         .float 0.0
  39         .byte 0, 0, 0, 0x80
  40         ASM_SIZE_DIRECTIVE(huge_nan_null_null)
  41
  42         ASM_TYPE_DIRECTIVE(twopi,@object)
  43 twopi:
  44         .byte 0x35, 0xc2, 0x68, 0x21, 0xa2, 0xda, 0xf, 0xc9, 0x1, 0x40
  45         .byte 0, 0, 0, 0, 0, 0
  46         ASM_SIZE_DIRECTIVE(twopi)
  47
  48         ASM_TYPE_DIRECTIVE(l2e,@object)
  49 l2e:
  50         .byte 0xbc, 0xf0, 0x17, 0x5c, 0x29, 0x3b, 0xaa, 0xb8, 0xff, 0x3f
  51         .byte 0, 0, 0, 0, 0, 0
  52         ASM_SIZE_DIRECTIVE(l2e)
  53
  54         ASM_TYPE_DIRECTIVE(one,@object)
  55 one:    .double 1.0
  56         ASM_SIZE_DIRECTIVE(one)
  57
  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
  68 ENTRY(__cexpf)
  69         flds    4(%esp)                 /* x */
  70         fxam
  71         fnstsw
  72         flds    8(%esp)                 /* y : x */
  73 #ifdef  PIC
  74         LOAD_PIC_REG (cx)
  75 #endif
  76         movb    %ah, %dh
  77         andb    $0x45, %ah
  78         cmpb    $0x05, %ah
  79         je      1f                      /* Jump if real part is +-Inf */
  80         cmpb    $0x01, %ah
  81         je      2f                      /* Jump if real part is NaN */
  82
  83         fxam                            /* y : x */
  84         fnstsw
  85         /* If the imaginary part is not finite we return NaN+i NaN, as
  86            for the case when the real part is NaN.  A test for +-Inf and
  87            NaN would be necessary.  But since we know the stack register
  88            we applied `fxam' to is not empty we can simply use one test.
  89            Check your FPU manual for more information.  */
  90         andb    $0x01, %ah
  91         cmpb    $0x01, %ah
  92         je      20f
  93
  94         /* We have finite numbers in the real and imaginary part.  Do
  95            the real work now.  */
  96         fxch                    /* x : y */
  97         fldt    MO(l2e)         /* log2(e) : x : y */
  98         fmulp                   /* x * log2(e) : y */
  99         fld     %st             /* x * log2(e) : x * log2(e) : y */
 100         frndint                 /* int(x * log2(e)) : x * log2(e) : y */
 101         fsubr   %st, %st(1)     /* int(x * log2(e)) : frac(x * log2(e)) : y */
 102         fxch                    /* frac(x * log2(e)) : int(x * log2(e)) : y */
 103         f2xm1                   /* 2^frac(x * log2(e))-1 : int(x * log2(e)) : y */
 104         faddl   MO(one)         /* 2^frac(x * log2(e)) : int(x * log2(e)) : y */
 105         fscale                  /* e^x : int(x * log2(e)) : y */
 106         fst     %st(1)          /* e^x : e^x : y */
 107         fxch    %st(2)          /* y : e^x : e^x */
 108         fsincos                 /* cos(y) : sin(y) : e^x : e^x */
 109         fnstsw
 110         testl   $0x400, %eax
 111         jnz     7f
 112         fmulp   %st, %st(3)     /* sin(y) : e^x : e^x * cos(y) */
 113         fmulp   %st, %st(1)     /* e^x * sin(y) : e^x * cos(y) */
 114         subl    $8, %esp
 115         cfi_adjust_cfa_offset (8)
 116         fstps   4(%esp)
 117         fstps   (%esp)
 118         popl    %eax
 119         cfi_adjust_cfa_offset (-4)
 120         popl    %edx
 121         cfi_adjust_cfa_offset (-4)
 122         ret
 123
 124         /* We have to reduce the argument to fsincos.  */
 125         .align ALIGNARG(4)
 126 7:      fldt    MO(twopi)       /* 2*pi : y : e^x : e^x */
 127         fxch                    /* y : 2*pi : e^x : e^x */
 128 8:      fprem1                  /* y%(2*pi) : 2*pi : e^x : e^x */
 129         fnstsw
 130         testl   $0x400, %eax
 131         jnz     8b
 132         fstp    %st(1)          /* y%(2*pi) : e^x : e^x */
 133         fsincos                 /* cos(y) : sin(y) : e^x : e^x */
 134         fmulp   %st, %st(3)
 135         fmulp   %st, %st(1)
 136         subl    $8, %esp
 137         cfi_adjust_cfa_offset (8)
 138         fstps   4(%esp)
 139         fstps   (%esp)
 140         popl    %eax
 141         cfi_adjust_cfa_offset (-4)
 142         popl    %edx
 143         cfi_adjust_cfa_offset (-4)
 144         ret
 145
 146         /* The real part is +-inf.  We must make further differences.  */
 147         .align ALIGNARG(4)
 148 1:      fxam                    /* y : x */
 149         fnstsw
 150         movb    %ah, %dl
 151         testb   $0x01, %ah      /* See above why 0x01 is usable here.  */
 152         jne     3f
 153
 154
 155         /* The real part is +-Inf and the imaginary part is finite.  */
 156         andl    $0x245, %edx
 157         cmpb    $0x40, %dl      /* Imaginary part == 0?  */
 158         je      4f              /* Yes ->  */
 159
 160         fxch                    /* x : y */
 161         shrl    $6, %edx
 162         fstp    %st(0)          /* y */ /* Drop the real part.  */
 163         andl    $8, %edx        /* This puts the sign bit of the real part
 164                                    in bit 3.  So we can use it to index a
 165                                    small array to select 0 or Inf.  */
 166         fsincos                 /* cos(y) : sin(y) */
 167         fnstsw
 168         testl   $0x0400, %eax
 169         jnz     5f
 170         fxch
 171         ftst
 172         fnstsw
 173         fstp    %st(0)
 174         shll    $23, %eax
 175         andl    $0x80000000, %eax
 176         orl     MOX(huge_nan_null_null,%edx,1), %eax
 177         movl    MOX(huge_nan_null_null,%edx,1), %ecx
 178         movl    %eax, %edx
 179         ftst
 180         fnstsw
 181         fstp    %st(0)
 182         shll    $23, %eax
 183         andl    $0x80000000, %eax
 184         orl     %ecx, %eax
 185         ret
 186         /* We must reduce the argument to fsincos.  */
 187         .align ALIGNARG(4)
 188 5:      fldt    MO(twopi)
 189         fxch
 190 6:      fprem1
 191         fnstsw
 192         testl   $0x400, %eax
 193         jnz     6b
 194         fstp    %st(1)
 195         fsincos
 196         fxch
 197         ftst
 198         fnstsw
 199         fstp    %st(0)
 200         shll    $23, %eax
 201         andl    $0x80000000, %eax
 202         orl     MOX(huge_nan_null_null,%edx,1), %eax
 203         movl    MOX(huge_nan_null_null,%edx,1), %ecx
 204         movl    %eax, %edx
 205         ftst
 206         fnstsw
 207         fstp    %st(0)
 208         shll    $23, %eax
 209         andl    $0x80000000, %eax
 210         orl     %ecx, %eax
 211         ret
 212
 213         /* The real part is +-Inf and the imaginary part is +-0.  So return
 214            +-Inf+-0i.  */
 215         .align ALIGNARG(4)
 216 4:      subl    $4, %esp
 217         cfi_adjust_cfa_offset (4)
 218         fstps   (%esp)
 219         shrl    $6, %edx
 220         fstp    %st(0)
 221         andl    $8, %edx
 222         movl    MOX(huge_nan_null_null,%edx,1), %eax
 223         popl    %edx
 224         cfi_adjust_cfa_offset (-4)
 225         ret
 226
 227         /* The real part is +-Inf, the imaginary is also is not finite.  */
 228         .align ALIGNARG(4)
 229 3:      fstp    %st(0)
 230         fstp    %st(0)          /* <empty> */
 231         andb    $0x45, %ah
 232         andb    $0x47, %dh
 233         xorb    %dh, %ah
 234         jnz     30f
 235         flds    MO(infinity)    /* Raise invalid exception.  */
 236         fmuls   MO(zero)
 237         fstp    %st(0)
 238 30:     movl    %edx, %eax
 239         shrl    $6, %edx
 240         shll    $3, %eax
 241         andl    $8, %edx
 242         andl    $16, %eax
 243         orl     %eax, %edx
 244
 245         movl    MOX(huge_nan_null_null,%edx,1), %eax
 246         movl    MOX(huge_nan_null_null+4,%edx,1), %edx
 247         ret
 248
 249         /* The real part is NaN.  */
 250         .align ALIGNARG(4)
 251 20:     flds    MO(infinity)            /* Raise invalid exception.  */
 252         fmuls   MO(zero)
 253         fstp    %st(0)
 254 2:      fstp    %st(0)
 255         fstp    %st(0)
 256         movl    MO(huge_nan_null_null+4), %eax
 257         movl    %eax, %edx
 258         ret
 259
 260 END(__cexpf)
 261 weak_alias (__cexpf, cexpf)