sysdeps/i386/fpu/s_cexpl.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(__cexpl)
  69         fldt    8(%esp)                 /* x */
  70         fxam
  71         fnstsw
  72         fldt    20(%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         fstpt   12(%eax)
 116         fstpt   (%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         fstpt   12(%eax)
 133         fstpt   (%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         fld     %st
 163         fstpt   12(%edx)
 164         fstpt   (%edx)
 165         ftst
 166         fnstsw
 167         shll    $7, %eax
 168         andl    $0x8000, %eax
 169         orl     %eax, 8(%edx)
 170         fstp    %st(0)
 171         ftst
 172         fnstsw
 173         shll    $7, %eax
 174         andl    $0x8000, %eax
 175         orl     %eax, 20(%edx)
 176         fstp    %st(0)
 177         ret     $4
 178         /* We must reduce the argument to fsincos.  */
 179         .align ALIGNARG(4)
 180 5:      fldt    MO(twopi)
 181         fxch
 182 6:      fprem1
 183         fnstsw
 184         testl   $0x400, %eax
 185         jnz     6b
 186         fstp    %st(1)
 187         fsincos
 188         fldl    MOX(huge_nan_null_null,%edx,1)
 189         movl    4(%esp), %edx           /* Pointer to memory for result.  */
 190         fld     %st
 191         fstpt   12(%edx)
 192         fstpt   (%edx)
 193         ftst
 194         fnstsw
 195         shll    $7, %eax
 196         andl    $0x8000, %eax
 197         orl     %eax, 8(%edx)
 198         fstp    %st(0)
 199         ftst
 200         fnstsw
 201         shll    $7, %eax
 202         andl    $0x8000, %eax
 203         orl     %eax, 20(%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         fstpt   12(%eax)
 212         shrl    $5, %edx
 213         fstp    %st(0)
 214         andl    $16, %edx
 215         fldl    MOX(huge_nan_null_null,%edx,1)
 216         fstpt   (%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         fstpt   (%eax)
 242         fstpt   12(%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         fld     %st(0)
 255         fstpt   (%eax)
 256         fstpt   12(%eax)
 257         ret     $4
 258
 259 END(__cexpl)
 260 weak_alias (__cexpl, cexpl)