sysdeps/i386/fpu/s_cexp.S

   1 /* ix87 specific implementation of complex exponential function for double.
   2    Copyright (C) 1997, 2005, 2012 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         .section .rodata
  24
  25         .align ALIGNARG(4)
  26         ASM_TYPE_DIRECTIVE(huge_nan_null_null,@object)
  27 huge_nan_null_null:
  28         .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
  29         .byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f
  30         .double 0.0
  31 zero:   .double 0.0
  32 infinity:
  33         .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
  34         .byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f
  35         .double 0.0
  36         .byte 0, 0, 0, 0, 0, 0, 0, 0x80
  37         ASM_SIZE_DIRECTIVE(huge_nan_null_null)
  38
  39         ASM_TYPE_DIRECTIVE(twopi,@object)
  40 twopi:
  41         .byte 0x35, 0xc2, 0x68, 0x21, 0xa2, 0xda, 0xf, 0xc9, 0x1, 0x40
  42         .byte 0, 0, 0, 0, 0, 0
  43         ASM_SIZE_DIRECTIVE(twopi)
  44
  45         ASM_TYPE_DIRECTIVE(l2e,@object)
  46 l2e:
  47         .byte 0xbc, 0xf0, 0x17, 0x5c, 0x29, 0x3b, 0xaa, 0xb8, 0xff, 0x3f
  48         .byte 0, 0, 0, 0, 0, 0
  49         ASM_SIZE_DIRECTIVE(l2e)
  50
  51         ASM_TYPE_DIRECTIVE(one,@object)
  52 one:    .double 1.0
  53         ASM_SIZE_DIRECTIVE(one)
  54
  55
  56 #ifdef PIC
  57 #define MO(op) op##@GOTOFF(%ecx)
  58 #define MOX(op,x,f) op##@GOTOFF(%ecx,x,f)
  59 #else
  60 #define MO(op) op
  61 #define MOX(op,x,f) op(,x,f)
  62 #endif
  63
  64         .text
  65 ENTRY(__cexp)
  66         fldl    8(%esp)                 /* x */
  67         fxam
  68         fnstsw
  69         fldl    16(%esp)                /* y : x */
  70 #ifdef  PIC
  71         LOAD_PIC_REG (cx)
  72 #endif
  73         movb    %ah, %dh
  74         andb    $0x45, %ah
  75         cmpb    $0x05, %ah
  76         je      1f                      /* Jump if real part is +-Inf */
  77         cmpb    $0x01, %ah
  78         je      2f                      /* Jump if real part is NaN */
  79
  80         fxam                            /* y : x */
  81         fnstsw
  82         /* If the imaginary part is not finite we return NaN+i NaN, as
  83            for the case when the real part is NaN.  A test for +-Inf and
  84            NaN would be necessary.  But since we know the stack register
  85            we applied `fxam' to is not empty we can simply use one test.
  86            Check your FPU manual for more information.  */
  87         andb    $0x01, %ah
  88         cmpb    $0x01, %ah
  89         je      20f
  90
  91         /* We have finite numbers in the real and imaginary part.  Do
  92            the real work now.  */
  93         fxch                    /* x : y */
  94         fldt    MO(l2e)         /* log2(e) : x : y */
  95         fmulp                   /* x * log2(e) : y */
  96         fld     %st             /* x * log2(e) : x * log2(e) : y */
  97         frndint                 /* int(x * log2(e)) : x * log2(e) : y */
  98         fsubr   %st, %st(1)     /* int(x * log2(e)) : frac(x * log2(e)) : y */
  99         fxch                    /* frac(x * log2(e)) : int(x * log2(e)) : y */
 100         f2xm1                   /* 2^frac(x * log2(e))-1 : int(x * log2(e)) : y */
 101         faddl   MO(one)         /* 2^frac(x * log2(e)) : int(x * log2(e)) : y */
 102         fscale                  /* e^x : int(x * log2(e)) : y */
 103         fst     %st(1)          /* e^x : e^x : y */
 104         fxch    %st(2)          /* y : e^x : e^x */
 105         fsincos                 /* cos(y) : sin(y) : e^x : e^x */
 106         fnstsw
 107         testl   $0x400, %eax
 108         jnz     7f
 109         fmulp   %st, %st(3)     /* sin(y) : e^x : e^x * cos(y) */
 110         fmulp   %st, %st(1)     /* e^x * sin(y) : e^x * cos(y) */
 111         movl    4(%esp), %eax           /* Pointer to memory for result.  */
 112         fstpl   8(%eax)
 113         fstpl   (%eax)
 114         ret     $4
 115
 116         /* We have to reduce the argument to fsincos.  */
 117         .align ALIGNARG(4)
 118 7:      fldt    MO(twopi)       /* 2*pi : y : e^x : e^x */
 119         fxch                    /* y : 2*pi : e^x : e^x */
 120 8:      fprem1                  /* y%(2*pi) : 2*pi : e^x : e^x */
 121         fnstsw
 122         testl   $0x400, %eax
 123         jnz     8b
 124         fstp    %st(1)          /* y%(2*pi) : e^x : e^x */
 125         fsincos                 /* cos(y) : sin(y) : e^x : e^x */
 126         fmulp   %st, %st(3)
 127         fmulp   %st, %st(1)
 128         movl    4(%esp), %eax           /* Pointer to memory for result.  */
 129         fstpl   8(%eax)
 130         fstpl   (%eax)
 131         ret     $4
 132
 133         /* The real part is +-inf.  We must make further differences.  */
 134         .align ALIGNARG(4)
 135 1:      fxam                    /* y : x */
 136         fnstsw
 137         movb    %ah, %dl
 138         testb   $0x01, %ah      /* See above why 0x01 is usable here.  */
 139         jne     3f
 140
 141
 142         /* The real part is +-Inf and the imaginary part is finite.  */
 143         andl    $0x245, %edx
 144         cmpb    $0x40, %dl      /* Imaginary part == 0?  */
 145         je      4f              /* Yes ->  */
 146
 147         fxch                    /* x : y */
 148         shrl    $5, %edx
 149         fstp    %st(0)          /* y */ /* Drop the real part.  */
 150         andl    $16, %edx       /* This puts the sign bit of the real part
 151                                    in bit 4.  So we can use it to index a
 152                                    small array to select 0 or Inf.  */
 153         fsincos                 /* cos(y) : sin(y) */
 154         fnstsw
 155         testl   $0x0400, %eax
 156         jnz     5f
 157         fldl    MOX(huge_nan_null_null,%edx,1)
 158         movl    4(%esp), %edx           /* Pointer to memory for result.  */
 159         fstl    8(%edx)
 160         fstpl   (%edx)
 161         ftst
 162         fnstsw
 163         shll    $23, %eax
 164         andl    $0x80000000, %eax
 165         orl     %eax, 4(%edx)
 166         fstp    %st(0)
 167         ftst
 168         fnstsw
 169         shll    $23, %eax
 170         andl    $0x80000000, %eax
 171         orl     %eax, 12(%edx)
 172         fstp    %st(0)
 173         ret     $4
 174         /* We must reduce the argument to fsincos.  */
 175         .align ALIGNARG(4)
 176 5:      fldt    MO(twopi)
 177         fxch
 178 6:      fprem1
 179         fnstsw
 180         testl   $0x400, %eax
 181         jnz     6b
 182         fstp    %st(1)
 183         fsincos
 184         fldl    MOX(huge_nan_null_null,%edx,1)
 185         movl    4(%esp), %edx           /* Pointer to memory for result.  */
 186         fstl    8(%edx)
 187         fstpl   (%edx)
 188         ftst
 189         fnstsw
 190         shll    $23, %eax
 191         andl    $0x80000000, %eax
 192         orl     %eax, 4(%edx)
 193         fstp    %st(0)
 194         ftst
 195         fnstsw
 196         shll    $23, %eax
 197         andl    $0x80000000, %eax
 198         orl     %eax, 12(%edx)
 199         fstp    %st(0)
 200         ret     $4
 201
 202         /* The real part is +-Inf and the imaginary part is +-0.  So return
 203            +-Inf+-0i.  */
 204         .align ALIGNARG(4)
 205 4:      movl    4(%esp), %eax           /* Pointer to memory for result.  */
 206         fstpl   8(%eax)
 207         shrl    $5, %edx
 208         fstp    %st(0)
 209         andl    $16, %edx
 210         fldl    MOX(huge_nan_null_null,%edx,1)
 211         fstpl   (%eax)
 212         ret     $4
 213
 214         /* The real part is +-Inf, the imaginary is also is not finite.  */
 215         .align ALIGNARG(4)
 216 3:      fstp    %st(0)
 217         fstp    %st(0)          /* <empty> */
 218         andb    $0x45, %ah
 219         andb    $0x47, %dh
 220         xorb    %dh, %ah
 221         jnz     30f
 222         fldl    MO(infinity)    /* Raise invalid exception.  */
 223         fmull   MO(zero)
 224         fstp    %st(0)
 225 30:     movl    %edx, %eax
 226         shrl    $5, %edx
 227         shll    $4, %eax
 228         andl    $16, %edx
 229         andl    $32, %eax
 230         orl     %eax, %edx
 231         movl    4(%esp), %eax           /* Pointer to memory for result.  */
 232
 233         fldl    MOX(huge_nan_null_null,%edx,1)
 234         fldl    MOX(huge_nan_null_null+8,%edx,1)
 235         fxch
 236         fstpl   (%eax)
 237         fstpl   8(%eax)
 238         ret     $4
 239
 240         /* The real part is NaN.  */
 241         .align ALIGNARG(4)
 242 20:     fldl    MO(infinity)            /* Raise invalid exception.  */
 243         fmull   MO(zero)
 244         fstp    %st(0)
 245 2:      fstp    %st(0)
 246         fstp    %st(0)
 247         movl    4(%esp), %eax           /* Pointer to memory for result.  */
 248         fldl    MO(huge_nan_null_null+8)
 249         fstl    (%eax)
 250         fstpl   8(%eax)
 251         ret     $4
 252
 253 END(__cexp)
 254 weak_alias (__cexp, cexp)