sysdeps/i386/fpu/s_cexpl.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(__cexpl)
  66         fldt    8(%esp)                 /* x */
  67         fxam
  68         fnstsw
  69         fldt    20(%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         fstpt   12(%eax)
 113         fstpt   (%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         fstpt   12(%eax)
 130         fstpt   (%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         fld     %st
 160         fstpt   12(%edx)
 161         fstpt   (%edx)
 162         ftst
 163         fnstsw
 164         shll    $7, %eax
 165         andl    $0x8000, %eax
 166         orl     %eax, 8(%edx)
 167         fstp    %st(0)
 168         ftst
 169         fnstsw
 170         shll    $7, %eax
 171         andl    $0x8000, %eax
 172         orl     %eax, 20(%edx)
 173         fstp    %st(0)
 174         ret     $4
 175         /* We must reduce the argument to fsincos.  */
 176         .align ALIGNARG(4)
 177 5:      fldt    MO(twopi)
 178         fxch
 179 6:      fprem1
 180         fnstsw
 181         testl   $0x400, %eax
 182         jnz     6b
 183         fstp    %st(1)
 184         fsincos
 185         fldl    MOX(huge_nan_null_null,%edx,1)
 186         movl    4(%esp), %edx           /* Pointer to memory for result.  */
 187         fld     %st
 188         fstpt   12(%edx)
 189         fstpt   (%edx)
 190         ftst
 191         fnstsw
 192         shll    $7, %eax
 193         andl    $0x8000, %eax
 194         orl     %eax, 8(%edx)
 195         fstp    %st(0)
 196         ftst
 197         fnstsw
 198         shll    $7, %eax
 199         andl    $0x8000, %eax
 200         orl     %eax, 20(%edx)
 201         fstp    %st(0)
 202         ret     $4
 203
 204         /* The real part is +-Inf and the imaginary part is +-0.  So return
 205            +-Inf+-0i.  */
 206         .align ALIGNARG(4)
 207 4:      movl    4(%esp), %eax           /* Pointer to memory for result.  */
 208         fstpt   12(%eax)
 209         shrl    $5, %edx
 210         fstp    %st(0)
 211         andl    $16, %edx
 212         fldl    MOX(huge_nan_null_null,%edx,1)
 213         fstpt   (%eax)
 214         ret     $4
 215
 216         /* The real part is +-Inf, the imaginary is also is not finite.  */
 217         .align ALIGNARG(4)
 218 3:      fstp    %st(0)
 219         fstp    %st(0)          /* <empty> */
 220         andb    $0x45, %ah
 221         andb    $0x47, %dh
 222         xorb    %dh, %ah
 223         jnz     30f
 224         fldl    MO(infinity)    /* Raise invalid exception.  */
 225         fmull   MO(zero)
 226         fstp    %st(0)
 227 30:     movl    %edx, %eax
 228         shrl    $5, %edx
 229         shll    $4, %eax
 230         andl    $16, %edx
 231         andl    $32, %eax
 232         orl     %eax, %edx
 233         movl    4(%esp), %eax           /* Pointer to memory for result.  */
 234
 235         fldl    MOX(huge_nan_null_null,%edx,1)
 236         fldl    MOX(huge_nan_null_null+8,%edx,1)
 237         fxch
 238         fstpt   (%eax)
 239         fstpt   12(%eax)
 240         ret     $4
 241
 242         /* The real part is NaN.  */
 243         .align ALIGNARG(4)
 244 20:     fldl    MO(infinity)            /* Raise invalid exception.  */
 245         fmull   MO(zero)
 246         fstp    %st(0)
 247 2:      fstp    %st(0)
 248         fstp    %st(0)
 249         movl    4(%esp), %eax           /* Pointer to memory for result.  */
 250         fldl    MO(huge_nan_null_null+8)
 251         fld     %st(0)
 252         fstpt   (%eax)
 253         fstpt   12(%eax)
 254         ret     $4
 255
 256 END(__cexpl)
 257 weak_alias (__cexpl, cexpl)