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 zero:   .double 0.0
  35 infinity:
  36         .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
  37         .byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f
  38         .double 0.0
  39         .byte 0, 0, 0, 0, 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(__cexp)
  69         fldl    8(%esp)                 /* x */
  70         fxam
  71         fnstsw
  72         fldl    16(%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         movl    4(%esp), %eax           /* Pointer to memory for result.  */
 117         fstpl   8(%eax)
 118         fstpl   (%eax)
 119         ret     $4
 120
 121         /* We have to reduce the argument to fsincos.  */
 122         .align ALIGNARG(4)
 123 7:      fldt    MO(twopi)       /* 2*pi : y : e^x : e^x */
 124         fxch                    /* y : 2*pi : e^x : e^x */
 125 8:      fprem1                  /* y%(2*pi) : 2*pi : e^x : e^x */
 126         fnstsw
 127         testl   $0x400, %eax
 128         jnz     8b
 129         fstp    %st(1)          /* y%(2*pi) : e^x : e^x */
 130         fsincos                 /* cos(y) : sin(y) : e^x : e^x */
 131         fmulp   %st, %st(3)
 132         fmulp   %st, %st(1)
 133         movl    4(%esp), %eax           /* Pointer to memory for result.  */
 134         fstpl   8(%eax)
 135         fstpl   (%eax)
 136         ret     $4
 137
 138         /* The real part is +-inf.  We must make further differences.  */
 139         .align ALIGNARG(4)
 140 1:      fxam                    /* y : x */
 141         fnstsw
 142         movb    %ah, %dl
 143         testb   $0x01, %ah      /* See above why 0x01 is usable here.  */
 144         jne     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         andb    $0x45, %ah
 224         andb    $0x47, %dh
 225         xorb    %dh, %ah
 226         jnz     30f
 227         fldl    MO(infinity)    /* Raise invalid exception.  */
 228         fmull   MO(zero)
 229         fstp    %st(0)
 230 30:     movl    %edx, %eax
 231         shrl    $5, %edx
 232         shll    $4, %eax
 233         andl    $16, %edx
 234         andl    $32, %eax
 235         orl     %eax, %edx
 236         movl    4(%esp), %eax           /* Pointer to memory for result.  */
 237
 238         fldl    MOX(huge_nan_null_null,%edx,1)
 239         fldl    MOX(huge_nan_null_null+8,%edx,1)
 240         fxch
 241         fstpl   (%eax)
 242         fstpl   8(%eax)
 243         ret     $4
 244
 245         /* The real part is NaN.  */
 246         .align ALIGNARG(4)
 247 20:     fldl    MO(infinity)            /* Raise invalid exception.  */
 248         fmull   MO(zero)
 249         fstp    %st(0)
 250 2:      fstp    %st(0)
 251         fstp    %st(0)
 252         movl    4(%esp), %eax           /* Pointer to memory for result.  */
 253         fldl    MO(huge_nan_null_null+8)
 254         fstl    (%eax)
 255         fstpl   8(%eax)
 256         ret     $4
 257
 258 END(__cexp)
 259 weak_alias (__cexp, cexp)