sysdeps/x86_64/fpu/e_powl.S

   1 /* ix87 specific implementation of pow function.
   2    Copyright (C) 1996, 1997, 1998, 1999, 2001, 2004 Free Software Foundation, Inc.
   3    This file is part of the GNU C Library.
   4    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996.
   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 <machine/asm.h>
  22
  23 #ifdef __ELF__
  24         .section .rodata
  25 #else
  26         .text
  27 #endif
  28
  29         .align ALIGNARG(4)
  30         ASM_TYPE_DIRECTIVE(infinity,@object)
  31 inf_zero:
  32 infinity:
  33         .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
  34         ASM_SIZE_DIRECTIVE(infinity)
  35         ASM_TYPE_DIRECTIVE(zero,@object)
  36 zero:   .double 0.0
  37         ASM_SIZE_DIRECTIVE(zero)
  38         ASM_TYPE_DIRECTIVE(minf_mzero,@object)
  39 minf_mzero:
  40 minfinity:
  41         .byte 0, 0, 0, 0, 0, 0, 0xf0, 0xff
  42 mzero:
  43         .byte 0, 0, 0, 0, 0, 0, 0, 0x80
  44         ASM_SIZE_DIRECTIVE(minf_mzero)
  45         ASM_TYPE_DIRECTIVE(one,@object)
  46 one:    .double 1.0
  47         ASM_SIZE_DIRECTIVE(one)
  48         ASM_TYPE_DIRECTIVE(limit,@object)
  49 limit:  .double 0.29
  50         ASM_SIZE_DIRECTIVE(limit)
  51         ASM_TYPE_DIRECTIVE(p63,@object)
  52 p63:
  53         .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43
  54         ASM_SIZE_DIRECTIVE(p63)
  55
  56 #ifdef PIC
  57 #define MO(op) op##(%rip)
  58 #else
  59 #define MO(op) op
  60 #endif
  61
  62         .text
  63 ENTRY(__ieee754_powl)
  64         fldt    24(%rsp)        // y
  65         fxam
  66
  67
  68         fnstsw
  69         movb    %ah, %dl
  70         andb    $0x45, %ah
  71         cmpb    $0x40, %ah      // is y == 0 ?
  72         je      11f
  73
  74         cmpb    $0x05, %ah      // is y == ±inf ?
  75         je      12f
  76
  77         cmpb    $0x01, %ah      // is y == NaN ?
  78         je      30f
  79
  80         fldt    8(%rsp)         // x : y
  81
  82         fxam
  83         fnstsw
  84         movb    %ah, %dh
  85         andb    $0x45, %ah
  86         cmpb    $0x40, %ah
  87         je      20f             // x is ±0
  88
  89         cmpb    $0x05, %ah
  90         je      15f             // x is ±inf
  91
  92         fxch                    // y : x
  93
  94         /* fistpll raises invalid exception for |y| >= 1L<<63.  */
  95         fldl    MO(p63)         // 1L<<63 : y : x
  96         fld     %st(1)          // y : 1L<<63 : y : x
  97         fabs                    // |y| : 1L<<63 : y : x
  98         fcomip  %st(1), %st     // 1L<<63 : y : x
  99         fstp    %st(0)          // y : x
 100         jnc     2f
 101
 102         /* First see whether `y' is a natural number.  In this case we
 103            can use a more precise algorithm.  */
 104         fld     %st             // y : y : x
 105         fistpll -8(%rsp)        // y : x
 106         fildll  -8(%rsp)        // int(y) : y : x
 107         fucomip %st(1),%st      // y : x
 108         jne     2f
 109
 110         /* OK, we have an integer value for y.  */
 111         mov     -8(%rsp),%eax
 112         mov     -4(%rsp),%edx
 113         orl     $0, %edx
 114         fstp    %st(0)          // x
 115         jns     4f              // y >= 0, jump
 116         fdivrl  MO(one)         // 1/x          (now referred to as x)
 117         negl    %eax
 118         adcl    $0, %edx
 119         negl    %edx
 120 4:      fldl    MO(one)         // 1 : x
 121         fxch
 122
 123 6:      shrdl   $1, %edx, %eax
 124         jnc     5f
 125         fxch
 126         fmul    %st(1)          // x : ST*x
 127         fxch
 128 5:      fmul    %st(0), %st     // x*x : ST*x
 129         shrl    $1, %edx
 130         movl    %eax, %ecx
 131         orl     %edx, %ecx
 132         jnz     6b
 133         fstp    %st(0)          // ST*x
 134         ret
 135
 136         /* y is ±NAN */
 137 30:     fldt    8(%rsp)         // x : y
 138         fldl    MO(one)         // 1.0 : x : y
 139         fucomip %st(1),%st      // x : y
 140         je      31f
 141         fxch                    // y : x
 142 31:     fstp    %st(1)
 143         ret
 144
 145         .align ALIGNARG(4)
 146 2:      /* y is a real number.  */
 147         fxch                    // x : y
 148         fldl    MO(one)         // 1.0 : x : y
 149         fld     %st(1)          // x : 1.0 : x : y
 150         fsub    %st(1)          // x-1 : 1.0 : x : y
 151         fabs                    // |x-1| : 1.0 : x : y
 152         fcompl  MO(limit)       // 1.0 : x : y
 153         fnstsw
 154         fxch                    // x : 1.0 : y
 155         test    $4500,%eax
 156         jz      7f
 157         fsub    %st(1)          // x-1 : 1.0 : y
 158         fyl2xp1                 // log2(x) : y
 159         jmp     8f
 160
 161 7:      fyl2x                   // log2(x) : y
 162 8:      fmul    %st(1)          // y*log2(x) : y
 163         fxam
 164         fnstsw
 165         andb    $0x45, %ah
 166         cmpb    $0x05, %ah      // is y*log2(x) == ±inf ?
 167         je      28f
 168         fst     %st(1)          // y*log2(x) : y*log2(x)
 169         frndint                 // int(y*log2(x)) : y*log2(x)
 170         fsubr   %st, %st(1)     // int(y*log2(x)) : fract(y*log2(x))
 171         fxch                    // fract(y*log2(x)) : int(y*log2(x))
 172         f2xm1                   // 2^fract(y*log2(x))-1 : int(y*log2(x))
 173         faddl   MO(one)         // 2^fract(y*log2(x)) : int(y*log2(x))
 174         fscale                  // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
 175         fstp    %st(1)          // 2^fract(y*log2(x))*2^int(y*log2(x))
 176         ret
 177
 178 28:     fstp    %st(1)          // y*log2(x)
 179         fldl    MO(one)         // 1 : y*log2(x)
 180         fscale                  // 2^(y*log2(x)) : y*log2(x)
 181         fstp    %st(1)          // 2^(y*log2(x))
 182         ret
 183
 184         // pow(x,±0) = 1
 185         .align ALIGNARG(4)
 186 11:     fstp    %st(0)          // pop y
 187         fldl    MO(one)
 188         ret
 189
 190         // y == ±inf
 191         .align ALIGNARG(4)
 192 12:     fstp    %st(0)          // pop y
 193         fldt    8(%rsp)         // x
 194         fabs
 195         fcompl  MO(one)         // < 1, == 1, or > 1
 196         fnstsw
 197         andb    $0x45, %ah
 198         cmpb    $0x45, %ah
 199         je      13f             // jump if x is NaN
 200
 201         cmpb    $0x40, %ah
 202         je      14f             // jump if |x| == 1
 203
 204         shlb    $1, %ah
 205         xorb    %ah, %dl
 206         andl    $2, %edx
 207 #ifdef PIC
 208         lea     inf_zero(%rip),%rcx
 209         fldl    (%rcx, %rdx, 4)
 210 #else
 211         fldl    inf_zero(,%rdx, 4)
 212 #endif
 213         ret
 214
 215         .align ALIGNARG(4)
 216 14:     fldl    MO(one)
 217         ret
 218
 219         .align ALIGNARG(4)
 220 13:     fldt    8(%rsp)         // load x == NaN
 221         ret
 222
 223         .align ALIGNARG(4)
 224         // x is ±inf
 225 15:     fstp    %st(0)          // y
 226         testb   $2, %dh
 227         jz      16f             // jump if x == +inf
 228
 229         // We must find out whether y is an odd integer.
 230         fld     %st             // y : y
 231         fistpll -8(%rsp)        // y
 232         fildll  -8(%rsp)        // int(y) : y
 233         fucomip %st(1),%st
 234         ffreep  %st             // <empty>
 235         jne     17f
 236
 237         // OK, the value is an integer, but is it odd?
 238         mov     -8(%rsp), %eax
 239         mov     -4(%rsp), %edx
 240         andb    $1, %al
 241         jz      18f             // jump if not odd
 242         // It's an odd integer.
 243         shrl    $31, %edx
 244 #ifdef PIC
 245         lea     minf_mzero(%rip),%rcx
 246         fldl    (%rcx, %rdx, 8)
 247 #else
 248         fldl    minf_mzero(,%rdx, 8)
 249 #endif
 250         ret
 251
 252         .align ALIGNARG(4)
 253 16:     fcompl  MO(zero)
 254         fnstsw
 255         shrl    $5, %eax
 256         andl    $8, %eax
 257 #ifdef PIC
 258         lea     inf_zero(%rip),%rcx
 259         fldl    (%rcx, %rax, 1)
 260 #else
 261         fldl    inf_zero(,%rax, 1)
 262 #endif
 263         ret
 264
 265         .align ALIGNARG(4)
 266 17:     shll    $30, %edx       // sign bit for y in right position
 267 18:     shrl    $31, %edx
 268 #ifdef PIC
 269         lea     inf_zero(%rip),%rcx
 270         fldl    (%rcx, %rdx, 8)
 271 #else
 272         fldl    inf_zero(,%rdx, 8)
 273 #endif
 274         ret
 275
 276         .align ALIGNARG(4)
 277         // x is ±0
 278 20:     fstp    %st(0)          // y
 279         testb   $2, %dl
 280         jz      21f             // y > 0
 281
 282         // x is ±0 and y is < 0.  We must find out whether y is an odd integer.
 283         testb   $2, %dh
 284         jz      25f
 285
 286         fld     %st             // y : y
 287         fistpll -8(%rsp)        // y
 288         fildll  -8(%rsp)        // int(y) : y
 289         fucomip %st(1),%st
 290         ffreep  %st             // <empty>
 291         jne     26f
 292
 293         // OK, the value is an integer, but is it odd?
 294         mov     -8(%rsp),%eax
 295         mov     -4(%rsp),%edx
 296         andb    $1, %al
 297         jz      27f             // jump if not odd
 298         // It's an odd integer.
 299         // Raise divide-by-zero exception and get minus infinity value.
 300         fldl    MO(one)
 301         fdivl   MO(zero)
 302         fchs
 303         ret
 304
 305 25:     fstp    %st(0)
 306 26:
 307 27:     // Raise divide-by-zero exception and get infinity value.
 308         fldl    MO(one)
 309         fdivl   MO(zero)
 310         ret
 311
 312         .align ALIGNARG(4)
 313         // x is ±0 and y is > 0.  We must find out whether y is an odd integer.
 314 21:     testb   $2, %dh
 315         jz      22f
 316
 317         fld     %st             // y : y
 318         fistpll -8(%rsp)        // y
 319         fildll  -8(%rsp)        // int(y) : y
 320         fucomip %st(1),%st
 321         ffreep  %st             // <empty>
 322         jne     23f
 323
 324         // OK, the value is an integer, but is it odd?
 325         mov     -8(%rsp),%eax
 326         mov     -4(%rsp),%edx
 327         andb    $1, %al
 328         jz      24f             // jump if not odd
 329         // It's an odd integer.
 330         fldl    MO(mzero)
 331         ret
 332
 333 22:     fstp    %st(0)
 334 23:
 335 24:     fldl    MO(zero)
 336         ret
 337
 338 END(__ieee754_powl)