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