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