sysdeps/libm-i387/s_cexp.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, 0, 0, 0, 0, 0xf0, 0x7f
  32         .byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f
  33         .double 0.0
  34         .double 0.0
  35         .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
  36         .byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f
  37         .double 0.0
  38         .byte 0, 0, 0, 0, 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(__cexp)
  68         fldl    8(%esp)                 /* x */
  69         fxam
  70         fnstsw
  71         fldl    16(%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         movl    4(%esp), %eax           /* Pointer to memory for result.  */
 116         fstpl   8(%eax)
 117         fstpl   (%eax)
 118         ret     $4
 119
 120         /* We have to reduce the argument to fsincos.  */
 121         .align ALIGNARG(4)
 122 7:      fldt    MO(twopi)       /* 2*pi : y : e^x : e^x */
 123         fxch                    /* y : 2*pi : e^x : e^x */
 124 8:      fprem1                  /* y%(2*pi) : 2*pi : e^x : e^x */
 125         fnstsw
 126         testl   $0x400, %eax
 127         jnz     8b
 128         fstp    %st(1)          /* y%(2*pi) : e^x : e^x */
 129         fsincos                 /* cos(y) : sin(y) : e^x : e^x */
 130         fmulp   %st, %st(3)
 131         fmulp   %st, %st(1)
 132         movl    4(%esp), %eax           /* Pointer to memory for result.  */
 133         fstpl   8(%eax)
 134         fstpl   (%eax)
 135         ret     $4
 136
 137         /* The real part is +-inf.  We must make further differences.  */
 138         .align ALIGNARG(4)
 139 1:      fxam                    /* y : x */
 140         fnstsw
 141         movb    %ah, %dl
 142         andb    $0x01, %ah      /* See above why 0x01 is usable here.  */
 143         cmpb    $0x01, %ah
 144         je      3f
 145
 146
 147         /* The real part is +-Inf and the imaginary part is finite.  */
 148         andl    $0x245, %edx
 149         cmpb    $0x40, %dl      /* Imaginary part == 0?  */
 150         je      4f              /* Yes ->  */
 151
 152         fxch                    /* x : y */
 153         shrl    $5, %edx
 154         fstp    %st(0)          /* y */ /* Drop the real part.  */
 155         andl    $16, %edx       /* This puts the sign bit of the real part
 156                                    in bit 4.  So we can use it to index a
 157                                    small array to select 0 or Inf.  */
 158         fsincos                 /* cos(y) : sin(y) */
 159         fnstsw
 160         testl   $0x0400, %eax
 161         jnz     5f
 162         fldl    MOX(huge_nan_null_null,%edx,1)
 163         movl    4(%esp), %edx           /* Pointer to memory for result.  */
 164         fstl    8(%edx)
 165         fstpl   (%edx)
 166         ftst
 167         fnstsw
 168         shll    $23, %eax
 169         andl    $0x80000000, %eax
 170         orl     %eax, 4(%edx)
 171         fstp    %st(0)
 172         ftst
 173         fnstsw
 174         shll    $23, %eax
 175         andl    $0x80000000, %eax
 176         orl     %eax, 12(%edx)
 177         fstp    %st(0)
 178         ret     $4
 179         /* We must reduce the argument to fsincos.  */
 180         .align ALIGNARG(4)
 181 5:      fldt    MO(twopi)
 182         fxch
 183 6:      fprem1
 184         fnstsw
 185         testl   $0x400, %eax
 186         jnz     6b
 187         fstp    %st(1)
 188         fsincos
 189         fldl    MOX(huge_nan_null_null,%edx,1)
 190         movl    4(%esp), %edx           /* Pointer to memory for result.  */
 191         fstl    8(%edx)
 192         fstpl   (%edx)
 193         ftst
 194         fnstsw
 195         shll    $23, %eax
 196         andl    $0x80000000, %eax
 197         orl     %eax, 4(%edx)
 198         fstp    %st(0)
 199         ftst
 200         fnstsw
 201         shll    $23, %eax
 202         andl    $0x80000000, %eax
 203         orl     %eax, 12(%edx)
 204         fstp    %st(0)
 205         ret     $4
 206
 207         /* The real part is +-Inf and the imaginary part is +-0.  So return
 208            +-Inf+-0i.  */
 209         .align ALIGNARG(4)
 210 4:      movl    4(%esp), %eax           /* Pointer to memory for result.  */
 211         fstpl   8(%eax)
 212         shrl    $5, %edx
 213         fstp    %st(0)
 214         andl    $16, %edx
 215         fldl    MOX(huge_nan_null_null,%edx,1)
 216         fstpl   (%eax)
 217         ret     $4
 218
 219         /* The real part is +-Inf, the imaginary is also is not finite.  */
 220         .align ALIGNARG(4)
 221 3:      fstp    %st(0)
 222         fstp    %st(0)          /* <empty> */
 223         movl    %edx, %eax
 224         shrl    $5, %edx
 225         shll    $4, %eax
 226         andl    $16, %edx
 227         andl    $32, %eax
 228         orl     %eax, %edx
 229         movl    4(%esp), %eax           /* Pointer to memory for result.  */
 230
 231         fldl    MOX(huge_nan_null_null,%edx,1)
 232         fldl    MOX(huge_nan_null_null+8,%edx,1)
 233         fstpl   8(%eax)
 234         fstpl   (%eax)
 235         ret     $4
 236
 237         /* The real part is NaN.  */
 238         .align ALIGNARG(4)
 239 2:      fstp    %st(0)
 240         fstp    %st(0)
 241         movl    4(%esp), %eax           /* Pointer to memory for result.  */
 242         fldl    MO(huge_nan_null_null+8)
 243         fstl    (%eax)
 244         fstpl   8(%eax)
 245         ret     $4
 246
 247 END(__cexp)
 248 weak_alias (__cexp, cexp)