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         .float  0.0
  35         .byte 0, 0, 0x80, 0x7f
  36         .byte 0, 0, 0xc0, 0x7f
  37         .float 0.0
  38         .byte 0, 0, 0, 0x80
  39         ASM_SIZE_DIRECTIVE(huge_nan_null_null)
  40
  41         ASM_TYPE_DIRECTIVE(twopi,@object)
  42 twopi:
  43         .byte 0x35, 0xc2, 0x68, 0x21, 0xa2, 0xda, 0xf, 0xc9, 0x1, 0x40
  44         .byte 0, 0, 0, 0, 0, 0
  45         ASM_SIZE_DIRECTIVE(twopi)
  46
  47         ASM_TYPE_DIRECTIVE(l2e,@object)
  48 l2e:
  49         .byte 0xbc, 0xf0, 0x17, 0x5c, 0x29, 0x3b, 0xaa, 0xb8, 0xff, 0x3f
  50         .byte 0, 0, 0, 0, 0, 0
  51         ASM_SIZE_DIRECTIVE(l2e)
  52
  53         ASM_TYPE_DIRECTIVE(one,@object)
  54 one:    .double 1.0
  55         ASM_SIZE_DIRECTIVE(one)
  56
  57
  58 #ifdef PIC
  59 #define MO(op) op##@GOTOFF(%ecx)
  60 #define MOX(op,x,f) op##@GOTOFF(%ecx,x,f)
  61 #else
  62 #define MO(op) op
  63 #define MOX(op,x,f) op(,x,f)
  64 #endif
  65
  66         .text
  67 ENTRY(__cexpf)
  68         flds    4(%esp)                 /* x */
  69         fxam
  70         fnstsw
  71         flds    8(%esp)                 /* y : x */
  72 #ifdef  PIC
  73         call    1f
  74 1:      popl    %ecx
  75         addl    $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx
  76 #endif
  77         movb    %ah, %dh
  78         andb    $0x45, %ah
  79         cmpb    $0x05, %ah
  80         je      1f                      /* Jump if real part is +-Inf */
  81         cmpb    $0x01, %ah
  82         je      2f                      /* Jump if real part is NaN */
  83
  84         fxam                            /* y : x */
  85         fnstsw
  86         /* If the imaginary part is not finite we return NaN+i NaN, as
  87            for the case when the real part is NaN.  A test for +-Inf and
  88            NaN would be necessary.  But since we know the stack register
  89            we applied `fxam' to is not empty we can simply use one test.
  90            Check your FPU manual for more information.  */
  91         andb    $0x01, %ah
  92         cmpb    $0x01, %ah
  93         je      2f
  94
  95         /* We have finite numbers in the real and imaginary part.  Do
  96            the real work now.  */
  97         fxch                    /* x : y */
  98         fldt    MO(l2e)         /* log2(e) : x : y */
  99         fmulp                   /* x * log2(e) : y */
 100         fld     %st             /* x * log2(e) : x * log2(e) : y */
 101         frndint                 /* int(x * log2(e)) : x * log2(e) : y */
 102         fsubr   %st, %st(1)     /* int(x * log2(e)) : frac(x * log2(e)) : y */
 103         fxch                    /* frac(x * log2(e)) : int(x * log2(e)) : y */
 104         f2xm1                   /* 2^frac(x * log2(e))-1 : int(x * log2(e)) : y */
 105         faddl   MO(one)         /* 2^frac(x * log2(e)) : int(x * log2(e)) : y */
 106         fscale                  /* e^x : int(x * log2(e)) : y */
 107         fst     %st(1)          /* e^x : e^x : y */
 108         fxch    %st(2)          /* y : e^x : e^x */
 109         fsincos                 /* cos(y) : sin(y) : e^x : e^x */
 110         fnstsw
 111         testl   $0x400, %eax
 112         jnz     7f
 113         fmulp   %st, %st(3)     /* sin(y) : e^x : e^x * cos(y) */
 114         fmulp   %st, %st(1)     /* e^x * sin(y) : e^x * cos(y) */
 115         subl    $8, %esp
 116         fstps   4(%esp)
 117         fstps   (%esp)
 118         popl    %eax
 119         popl    %edx
 120         ret
 121
 122         /* We have to reduce the argument to fsincos.  */
 123         .align ALIGNARG(4)
 124 7:      fldt    MO(twopi)       /* 2*pi : y : e^x : e^x */
 125         fxch                    /* y : 2*pi : e^x : e^x */
 126 8:      fprem1                  /* y%(2*pi) : 2*pi : e^x : e^x */
 127         fnstsw
 128         testl   $0x400, %eax
 129         jnz     8b
 130         fstp    %st(1)          /* y%(2*pi) : e^x : e^x */
 131         fsincos                 /* cos(y) : sin(y) : e^x : e^x */
 132         fmulp   %st, %st(3)
 133         fmulp   %st, %st(1)
 134         subl    $8, %esp
 135         fstps   4(%esp)
 136         fstps   (%esp)
 137         popl    %eax
 138         popl    %edx
 139         ret
 140
 141         /* The real part is +-inf.  We must make further differences.  */
 142         .align ALIGNARG(4)
 143 1:      fxam                    /* y : x */
 144         fnstsw
 145         movb    %ah, %dl
 146         andb    $0x01, %ah      /* See above why 0x01 is usable here.  */
 147         cmpb    $0x01, %ah
 148         je      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         movl    %edx, %eax
 226         shrl    $6, %edx
 227         shll    $3, %eax
 228         andl    $8, %edx
 229         andl    $16, %eax
 230         orl     %eax, %edx
 231
 232         movl    MOX(huge_nan_null_null,%edx,1), %eax
 233         movl    MOX(huge_nan_null_null+4,%edx,1), %edx
 234         ret
 235
 236         /* The real part is NaN.  */
 237         .align ALIGNARG(4)
 238 2:      fstp    %st(0)
 239         fstp    %st(0)
 240         movl    MO(huge_nan_null_null+4), %eax
 241         movl    %eax, %edx
 242         ret
 243
 244 END(__cexpf)
 245 weak_alias (__cexpf, cexpf)