sysdeps/libm-i387/s_cexpf.S

   1 /* ix87 specific implementation of complex exponential function for double.
   2    Copyright (C) 1997 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 Library General Public License as
   8    published by the Free Software Foundation; either version 2 of the
   9    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    Library General Public License for more details.
  15
  16    You should have received a copy of the GNU Library General Public
  17    License along with the GNU C Library; see the file COPYING.LIB.  If not,
  18    write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  19    Boston, MA 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         call    1f
  75 1:      popl    %ecx
  76         addl    $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx
  77 #endif
  78         movb    %ah, %dh
  79         andb    $0x45, %ah
  80         cmpb    $0x05, %ah
  81         je      1f                      /* Jump if real part is +-Inf */
  82         cmpb    $0x01, %ah
  83         je      2f                      /* Jump if real part is NaN */
  84
  85         fxam                            /* y : x */
  86         fnstsw
  87         /* If the imaginary part is not finite we return NaN+i NaN, as
  88            for the case when the real part is NaN.  A test for +-Inf and
  89            NaN would be necessary.  But since we know the stack register
  90            we applied `fxam' to is not empty we can simply use one test.
  91            Check your FPU manual for more information.  */
  92         andb    $0x01, %ah
  93         cmpb    $0x01, %ah
  94         je      20f
  95
  96         /* We have finite numbers in the real and imaginary part.  Do
  97            the real work now.  */
  98         fxch                    /* x : y */
  99         fldt    MO(l2e)         /* log2(e) : x : y */
 100         fmulp                   /* x * log2(e) : y */
 101         fld     %st             /* x * log2(e) : x * log2(e) : y */
 102         frndint                 /* int(x * log2(e)) : x * log2(e) : y */
 103         fsubr   %st, %st(1)     /* int(x * log2(e)) : frac(x * log2(e)) : y */
 104         fxch                    /* frac(x * log2(e)) : int(x * log2(e)) : y */
 105         f2xm1                   /* 2^frac(x * log2(e))-1 : int(x * log2(e)) : y */
 106         faddl   MO(one)         /* 2^frac(x * log2(e)) : int(x * log2(e)) : y */
 107         fscale                  /* e^x : int(x * log2(e)) : y */
 108         fst     %st(1)          /* e^x : e^x : y */
 109         fxch    %st(2)          /* y : e^x : e^x */
 110         fsincos                 /* cos(y) : sin(y) : e^x : e^x */
 111         fnstsw
 112         testl   $0x400, %eax
 113         jnz     7f
 114         fmulp   %st, %st(3)     /* sin(y) : e^x : e^x * cos(y) */
 115         fmulp   %st, %st(1)     /* e^x * sin(y) : e^x * cos(y) */
 116         subl    $8, %esp
 117         fstps   4(%esp)
 118         fstps   (%esp)
 119         popl    %eax
 120         popl    %edx
 121         ret
 122
 123         /* We have to reduce the argument to fsincos.  */
 124         .align ALIGNARG(4)
 125 7:      fldt    MO(twopi)       /* 2*pi : y : e^x : e^x */
 126         fxch                    /* y : 2*pi : e^x : e^x */
 127 8:      fprem1                  /* y%(2*pi) : 2*pi : e^x : e^x */
 128         fnstsw
 129         testl   $0x400, %eax
 130         jnz     8b
 131         fstp    %st(1)          /* y%(2*pi) : e^x : e^x */
 132         fsincos                 /* cos(y) : sin(y) : e^x : e^x */
 133         fmulp   %st, %st(3)
 134         fmulp   %st, %st(1)
 135         subl    $8, %esp
 136         fstps   4(%esp)
 137         fstps   (%esp)
 138         popl    %eax
 139         popl    %edx
 140         ret
 141
 142         /* The real part is +-inf.  We must make further differences.  */
 143         .align ALIGNARG(4)
 144 1:      fxam                    /* y : x */
 145         fnstsw
 146         movb    %ah, %dl
 147         testb   $0x01, %ah      /* See above why 0x01 is usable here.  */
 148         jne     3f
 149
 150
 151         /* The real part is +-Inf and the imaginary part is finite.  */
 152         andl    $0x245, %edx
 153         cmpb    $0x40, %dl      /* Imaginary part == 0?  */
 154         je      4f              /* Yes ->  */
 155
 156         fxch                    /* x : y */
 157         shrl    $6, %edx
 158         fstp    %st(0)          /* y */ /* Drop the real part.  */
 159         andl    $8, %edx        /* This puts the sign bit of the real part
 160                                    in bit 3.  So we can use it to index a
 161                                    small array to select 0 or Inf.  */
 162         fsincos                 /* cos(y) : sin(y) */
 163         fnstsw
 164         testl   $0x0400, %eax
 165         jnz     5f
 166         fxch
 167         ftst
 168         fnstsw
 169         fstp    %st(0)
 170         shll    $23, %eax
 171         andl    $0x80000000, %eax
 172         orl     MOX(huge_nan_null_null,%edx,1), %eax
 173         movl    MOX(huge_nan_null_null,%edx,1), %ecx
 174         movl    %eax, %edx
 175         ftst
 176         fnstsw
 177         fstp    %st(0)
 178         shll    $23, %eax
 179         andl    $0x80000000, %eax
 180         orl     %ecx, %eax
 181         ret
 182         /* We must reduce the argument to fsincos.  */
 183         .align ALIGNARG(4)
 184 5:      fldt    MO(twopi)
 185         fxch
 186 6:      fprem1
 187         fnstsw
 188         testl   $0x400, %eax
 189         jnz     6b
 190         fstp    %st(1)
 191         fsincos
 192         fxch
 193         ftst
 194         fnstsw
 195         fstp    %st(0)
 196         shll    $23, %eax
 197         andl    $0x80000000, %eax
 198         orl     MOX(huge_nan_null_null,%edx,1), %eax
 199         movl    MOX(huge_nan_null_null,%edx,1), %ecx
 200         movl    %eax, %edx
 201         ftst
 202         fnstsw
 203         fstp    %st(0)
 204         shll    $23, %eax
 205         andl    $0x80000000, %eax
 206         orl     %ecx, %eax
 207         ret
 208
 209         /* The real part is +-Inf and the imaginary part is +-0.  So return
 210            +-Inf+-0i.  */
 211         .align ALIGNARG(4)
 212 4:      subl    $4, %esp
 213         fstps   (%esp)
 214         shrl    $6, %edx
 215         fstp    %st(0)
 216         andl    $8, %edx
 217         movl    MOX(huge_nan_null_null,%edx,1), %eax
 218         popl    %edx
 219         ret
 220
 221         /* The real part is +-Inf, the imaginary is also is not finite.  */
 222         .align ALIGNARG(4)
 223 3:      fstp    %st(0)
 224         fstp    %st(0)          /* <empty> */
 225         andb    $0x45, %ah
 226         andb    $0x47, %dh
 227         xorb    %dh, %ah
 228         jnz     30f
 229         flds    MO(infinity)    /* Raise invalid exception.  */
 230         fmuls   MO(zero)
 231         fstp    %st(0)
 232 30:     movl    %edx, %eax
 233         shrl    $6, %edx
 234         shll    $3, %eax
 235         andl    $8, %edx
 236         andl    $16, %eax
 237         orl     %eax, %edx
 238
 239         movl    MOX(huge_nan_null_null,%edx,1), %eax
 240         movl    MOX(huge_nan_null_null+4,%edx,1), %edx
 241         ret
 242
 243         /* The real part is NaN.  */
 244         .align ALIGNARG(4)
 245 20:     flds    MO(infinity)            /* Raise invalid exception.  */
 246         fmuls   MO(zero)
 247         fstp    %st(0)
 248 2:      fstp    %st(0)
 249         fstp    %st(0)
 250         movl    MO(huge_nan_null_null+4), %eax
 251         movl    %eax, %edx
 252         ret
 253
 254 END(__cexpf)
 255 weak_alias (__cexpf, cexpf)