sysdeps/x86_64/fpu/e_powl.S

   1 /* ix87 specific implementation of pow function.
   2    Copyright (C) 1996-1999, 2001, 2004, 2007, 2011-2012
   3    Free Software Foundation, Inc.
   4    This file is part of the GNU C Library.
   5    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996.
   6
   7    The GNU C Library is free software; you can redistribute it and/or
   8    modify it under the terms of the GNU Lesser General Public
   9    License as published by the Free Software Foundation; either
  10    version 2.1 of the License, or (at your option) any later version.
  11
  12    The GNU C Library is distributed in the hope that it will be useful,
  13    but WITHOUT ANY WARRANTY; without even the implied warranty of
  14    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15    Lesser General Public License for more details.
  16
  17    You should have received a copy of the GNU Lesser General Public
  18    License along with the GNU C Library; if not, write to the Free
  19    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
  20    02111-1307 USA.  */
  21
  22 #include <machine/asm.h>
  23
  24         .section .rodata.cst8,"aM",@progbits,8
  25
  26         .p2align 3
  27         ASM_TYPE_DIRECTIVE(one,@object)
  28 one:    .double 1.0
  29         ASM_SIZE_DIRECTIVE(one)
  30         ASM_TYPE_DIRECTIVE(limit,@object)
  31 limit:  .double 0.29
  32         ASM_SIZE_DIRECTIVE(limit)
  33         ASM_TYPE_DIRECTIVE(p63,@object)
  34 p63:    .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43
  35         ASM_SIZE_DIRECTIVE(p63)
  36
  37         .section .rodata.cst16,"aM",@progbits,16
  38
  39         .p2align 3
  40         ASM_TYPE_DIRECTIVE(infinity,@object)
  41 inf_zero:
  42 infinity:
  43         .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
  44         ASM_SIZE_DIRECTIVE(infinity)
  45         ASM_TYPE_DIRECTIVE(zero,@object)
  46 zero:   .double 0.0
  47         ASM_SIZE_DIRECTIVE(zero)
  48         ASM_TYPE_DIRECTIVE(minf_mzero,@object)
  49 minf_mzero:
  50 minfinity:
  51         .byte 0, 0, 0, 0, 0, 0, 0xf0, 0xff
  52 mzero:
  53         .byte 0, 0, 0, 0, 0, 0, 0, 0x80
  54         ASM_SIZE_DIRECTIVE(minf_mzero)
  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         fldl    MO(limit)       // 0.29 : 1.0 : x : y
 150         fld     %st(2)          // x : 0.29 : 1.0 : x : y
 151         fsub    %st(2)          // x-1 : 0.29 : 1.0 : x : y
 152         fabs                    // |x-1| : 0.29 : 1.0 : x : y
 153         fucompp                 // 1.0 : x : y
 154         fnstsw
 155         fxch                    // x : 1.0 : y
 156         test    $0x4500,%eax
 157         jz      7f
 158         fsub    %st(1)          // x-1 : 1.0 : y
 159         fyl2xp1                 // log2(x) : y
 160         jmp     8f
 161
 162 7:      fyl2x                   // log2(x) : y
 163 8:      fmul    %st(1)          // y*log2(x) : y
 164         fxam
 165         fnstsw
 166         andb    $0x45, %ah
 167         cmpb    $0x05, %ah      // is y*log2(x) == ±inf ?
 168         je      28f
 169         fst     %st(1)          // y*log2(x) : y*log2(x)
 170         frndint                 // int(y*log2(x)) : y*log2(x)
 171         fsubr   %st, %st(1)     // int(y*log2(x)) : fract(y*log2(x))
 172         fxch                    // fract(y*log2(x)) : int(y*log2(x))
 173         f2xm1                   // 2^fract(y*log2(x))-1 : int(y*log2(x))
 174         faddl   MO(one)         // 2^fract(y*log2(x)) : int(y*log2(x))
 175         fscale                  // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
 176         fstp    %st(1)          // 2^fract(y*log2(x))*2^int(y*log2(x))
 177         ret
 178
 179 28:     fstp    %st(1)          // y*log2(x)
 180         fldl    MO(one)         // 1 : y*log2(x)
 181         fscale                  // 2^(y*log2(x)) : y*log2(x)
 182         fstp    %st(1)          // 2^(y*log2(x))
 183         ret
 184
 185         // pow(x,±0) = 1
 186         .align ALIGNARG(4)
 187 11:     fstp    %st(0)          // pop y
 188         fldl    MO(one)
 189         ret
 190
 191         // y == ±inf
 192         .align ALIGNARG(4)
 193 12:     fstp    %st(0)          // pop y
 194         fldl    MO(one)         // 1
 195         fldt    8(%rsp)         // x : 1
 196         fabs                    // abs(x) : 1
 197         fucompp                 // < 1, == 1, or > 1
 198         fnstsw
 199         andb    $0x45, %ah
 200         cmpb    $0x45, %ah
 201         je      13f             // jump if x is NaN
 202
 203         cmpb    $0x40, %ah
 204         je      14f             // jump if |x| == 1
 205
 206         shlb    $1, %ah
 207         xorb    %ah, %dl
 208         andl    $2, %edx
 209 #ifdef PIC
 210         lea     inf_zero(%rip),%rcx
 211         fldl    (%rcx, %rdx, 4)
 212 #else
 213         fldl    inf_zero(,%rdx, 4)
 214 #endif
 215         ret
 216
 217         .align ALIGNARG(4)
 218 14:     fldl    MO(one)
 219         ret
 220
 221         .align ALIGNARG(4)
 222 13:     fldt    8(%rsp)         // load x == NaN
 223         ret
 224
 225         .align ALIGNARG(4)
 226         // x is ±inf
 227 15:     fstp    %st(0)          // y
 228         testb   $2, %dh
 229         jz      16f             // jump if x == +inf
 230
 231         // We must find out whether y is an odd integer.
 232         fld     %st             // y : y
 233         fistpll -8(%rsp)        // y
 234         fildll  -8(%rsp)        // int(y) : y
 235         fucomip %st(1),%st
 236         ffreep  %st             // <empty>
 237         jne     17f
 238
 239         // OK, the value is an integer, but is it odd?
 240         mov     -8(%rsp), %eax
 241         mov     -4(%rsp), %edx
 242         andb    $1, %al
 243         jz      18f             // jump if not odd
 244         // It's an odd integer.
 245         shrl    $31, %edx
 246 #ifdef PIC
 247         lea     minf_mzero(%rip),%rcx
 248         fldl    (%rcx, %rdx, 8)
 249 #else
 250         fldl    minf_mzero(,%rdx, 8)
 251 #endif
 252         ret
 253
 254         .align ALIGNARG(4)
 255 16:     fcompl  MO(zero)
 256         fnstsw
 257         shrl    $5, %eax
 258         andl    $8, %eax
 259 #ifdef PIC
 260         lea     inf_zero(%rip),%rcx
 261         fldl    (%rcx, %rax, 1)
 262 #else
 263         fldl    inf_zero(,%rax, 1)
 264 #endif
 265         ret
 266
 267         .align ALIGNARG(4)
 268 17:     shll    $30, %edx       // sign bit for y in right position
 269 18:     shrl    $31, %edx
 270 #ifdef PIC
 271         lea     inf_zero(%rip),%rcx
 272         fldl    (%rcx, %rdx, 8)
 273 #else
 274         fldl    inf_zero(,%rdx, 8)
 275 #endif
 276         ret
 277
 278         .align ALIGNARG(4)
 279         // x is ±0
 280 20:     fstp    %st(0)          // y
 281         testb   $2, %dl
 282         jz      21f             // y > 0
 283
 284         // x is ±0 and y is < 0.  We must find out whether y is an odd integer.
 285         testb   $2, %dh
 286         jz      25f
 287
 288         fld     %st             // y : y
 289         fistpll -8(%rsp)        // y
 290         fildll  -8(%rsp)        // int(y) : y
 291         fucomip %st(1),%st
 292         ffreep  %st             // <empty>
 293         jne     26f
 294
 295         // OK, the value is an integer, but is it odd?
 296         mov     -8(%rsp),%eax
 297         mov     -4(%rsp),%edx
 298         andb    $1, %al
 299         jz      27f             // jump if not odd
 300         // It's an odd integer.
 301         // Raise divide-by-zero exception and get minus infinity value.
 302         fldl    MO(one)
 303         fdivl   MO(zero)
 304         fchs
 305         ret
 306
 307 25:     fstp    %st(0)
 308 26:
 309 27:     // Raise divide-by-zero exception and get infinity value.
 310         fldl    MO(one)
 311         fdivl   MO(zero)
 312         ret
 313
 314         .align ALIGNARG(4)
 315         // x is ±0 and y is > 0.  We must find out whether y is an odd integer.
 316 21:     testb   $2, %dh
 317         jz      22f
 318
 319         fld     %st             // y : y
 320         fistpll -8(%rsp)        // y
 321         fildll  -8(%rsp)        // int(y) : y
 322         fucomip %st(1),%st
 323         ffreep  %st             // <empty>
 324         jne     23f
 325
 326         // OK, the value is an integer, but is it odd?
 327         mov     -8(%rsp),%eax
 328         mov     -4(%rsp),%edx
 329         andb    $1, %al
 330         jz      24f             // jump if not odd
 331         // It's an odd integer.
 332         fldl    MO(mzero)
 333         ret
 334
 335 22:     fstp    %st(0)
 336 23:
 337 24:     fldl    MO(zero)
 338         ret
 339
 340 END(__ieee754_powl)
 341 strong_alias (__ieee754_powl, __powl_finite)