sysdeps/x86_64/fpu/e_powl.S

   1 /* ix87 specific implementation of pow function.
   2    Copyright (C) 1996-2016 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, see
  18    <http://www.gnu.org/licenses/>.  */
  19
  20 #include <machine/asm.h>
  21 #include <x86_64-math-asm.h>
  22
  23         .section .rodata.cst8,"aM",@progbits,8
  24
  25         .p2align 3
  26         .type one,@object
  27 one:    .double 1.0
  28         ASM_SIZE_DIRECTIVE(one)
  29         .type p2,@object
  30 p2:     .byte 0, 0, 0, 0, 0, 0, 0x10, 0x40
  31         ASM_SIZE_DIRECTIVE(p2)
  32         .type p63,@object
  33 p63:    .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43
  34         ASM_SIZE_DIRECTIVE(p63)
  35         .type p64,@object
  36 p64:    .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x43
  37         ASM_SIZE_DIRECTIVE(p64)
  38         .type p78,@object
  39 p78:    .byte 0, 0, 0, 0, 0, 0, 0xd0, 0x44
  40         ASM_SIZE_DIRECTIVE(p78)
  41         .type pm79,@object
  42 pm79:   .byte 0, 0, 0, 0, 0, 0, 0, 0x3b
  43         ASM_SIZE_DIRECTIVE(pm79)
  44
  45         .section .rodata.cst16,"aM",@progbits,16
  46
  47         .p2align 3
  48         .type infinity,@object
  49 inf_zero:
  50 infinity:
  51         .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
  52         ASM_SIZE_DIRECTIVE(infinity)
  53         .type zero,@object
  54 zero:   .double 0.0
  55         ASM_SIZE_DIRECTIVE(zero)
  56         .type minf_mzero,@object
  57 minf_mzero:
  58 minfinity:
  59         .byte 0, 0, 0, 0, 0, 0, 0xf0, 0xff
  60 mzero:
  61         .byte 0, 0, 0, 0, 0, 0, 0, 0x80
  62         ASM_SIZE_DIRECTIVE(minf_mzero)
  63 DEFINE_LDBL_MIN
  64
  65 #ifdef PIC
  66 # define MO(op) op##(%rip)
  67 #else
  68 # define MO(op) op
  69 #endif
  70
  71         .text
  72 ENTRY(__ieee754_powl)
  73         fldt    24(%rsp)        // y
  74         fxam
  75
  76
  77         fnstsw
  78         movb    %ah, %dl
  79         andb    $0x45, %ah
  80         cmpb    $0x40, %ah      // is y == 0 ?
  81         je      11f
  82
  83         cmpb    $0x05, %ah      // is y == ±inf ?
  84         je      12f
  85
  86         cmpb    $0x01, %ah      // is y == NaN ?
  87         je      30f
  88
  89         fldt    8(%rsp)         // x : y
  90
  91         fxam
  92         fnstsw
  93         movb    %ah, %dh
  94         andb    $0x45, %ah
  95         cmpb    $0x40, %ah
  96         je      20f             // x is ±0
  97
  98         cmpb    $0x05, %ah
  99         je      15f             // x is ±inf
 100
 101         cmpb    $0x01, %ah
 102         je      31f             // x is NaN
 103
 104         fxch                    // y : x
 105
 106         /* fistpll raises invalid exception for |y| >= 1L<<63.  */
 107         fldl    MO(p63)         // 1L<<63 : y : x
 108         fld     %st(1)          // y : 1L<<63 : y : x
 109         fabs                    // |y| : 1L<<63 : y : x
 110         fcomip  %st(1), %st     // 1L<<63 : y : x
 111         fstp    %st(0)          // y : x
 112         jnc     2f
 113
 114         /* First see whether `y' is a natural number.  In this case we
 115            can use a more precise algorithm.  */
 116         fld     %st             // y : y : x
 117         fistpll -8(%rsp)        // y : x
 118         fildll  -8(%rsp)        // int(y) : y : x
 119         fucomip %st(1),%st      // y : x
 120         je      9f
 121
 122         // If y has absolute value at most 0x1p-79, then any finite
 123         // nonzero x will result in 1.  Saturate y to those bounds to
 124         // avoid underflow in the calculation of y*log2(x).
 125         fldl    MO(pm79)        // 0x1p-79 : y : x
 126         fld     %st(1)          // y : 0x1p-79 : y : x
 127         fabs                    // |y| : 0x1p-79 : y : x
 128         fcomip  %st(1), %st     // 0x1p-79 : y : x
 129         fstp    %st(0)          // y : x
 130         jnc     3f
 131         fstp    %st(0)          // pop y
 132         fldl    MO(pm79)        // 0x1p-79 : x
 133         testb   $2, %dl
 134         jnz     3f              // y > 0
 135         fchs                    // -0x1p-79 : x
 136         jmp     3f
 137
 138 9:      /* OK, we have an integer value for y.  Unless very small
 139            (we use < 4), use the algorithm for real exponent to avoid
 140            accumulation of errors.  */
 141         fldl    MO(p2)          // 4 : y : x
 142         fld     %st(1)          // y : 4 : y : x
 143         fabs                    // |y| : 4 : y : x
 144         fcomip  %st(1), %st     // 4 : y : x
 145         fstp    %st(0)          // y : x
 146         jnc     3f
 147         mov     -8(%rsp),%eax
 148         mov     -4(%rsp),%edx
 149         orl     $0, %edx
 150         fstp    %st(0)          // x
 151         jns     4f              // y >= 0, jump
 152         fdivrl  MO(one)         // 1/x          (now referred to as x)
 153         negl    %eax
 154         adcl    $0, %edx
 155         negl    %edx
 156 4:      fldl    MO(one)         // 1 : x
 157         fxch
 158
 159         /* If y is even, take the absolute value of x.  Otherwise,
 160            ensure all intermediate values that might overflow have the
 161            sign of x.  */
 162         testb   $1, %al
 163         jnz     6f
 164         fabs
 165
 166 6:      shrdl   $1, %edx, %eax
 167         jnc     5f
 168         fxch
 169         fabs
 170         fmul    %st(1)          // x : ST*x
 171         fxch
 172 5:      fld     %st             // x : x : ST*x
 173         fabs                    // |x| : x : ST*x
 174         fmulp                   // |x|*x : ST*x
 175         shrl    $1, %edx
 176         movl    %eax, %ecx
 177         orl     %edx, %ecx
 178         jnz     6b
 179         fstp    %st(0)          // ST*x
 180         LDBL_CHECK_FORCE_UFLOW_NONNAN
 181         ret
 182
 183         /* y is ±NAN */
 184 30:     fldt    8(%rsp)         // x : y
 185         fldl    MO(one)         // 1.0 : x : y
 186         fucomip %st(1),%st      // x : y
 187         je      31f
 188         fxch                    // y : x
 189 31:     fstp    %st(1)
 190         ret
 191
 192         .align ALIGNARG(4)
 193 2:      // y is a large integer (absolute value at least 1L<<63).
 194         // If y has absolute value at least 1L<<78, then any finite
 195         // nonzero x will result in 0 (underflow), 1 or infinity (overflow).
 196         // Saturate y to those bounds to avoid overflow in the calculation
 197         // of y*log2(x).
 198         fldl    MO(p78)         // 1L<<78 : y : x
 199         fld     %st(1)          // y : 1L<<78 : y : x
 200         fabs                    // |y| : 1L<<78 : y : x
 201         fcomip  %st(1), %st     // 1L<<78 : y : x
 202         fstp    %st(0)          // y : x
 203         jc      3f
 204         fstp    %st(0)          // pop y
 205         fldl    MO(p78)         // 1L<<78 : x
 206         testb   $2, %dl
 207         jz      3f              // y > 0
 208         fchs                    // -(1L<<78) : x
 209         .align ALIGNARG(4)
 210 3:      /* y is a real number.  */
 211         subq    $40, %rsp
 212         cfi_adjust_cfa_offset (40)
 213         fstpt   16(%rsp)        // x
 214         fstpt   (%rsp)          // <empty>
 215         call    HIDDEN_JUMPTARGET (__powl_helper)       // <result>
 216         addq    $40, %rsp
 217         cfi_adjust_cfa_offset (-40)
 218         ret
 219
 220         // pow(x,±0) = 1
 221         .align ALIGNARG(4)
 222 11:     fstp    %st(0)          // pop y
 223         fldl    MO(one)
 224         ret
 225
 226         // y == ±inf
 227         .align ALIGNARG(4)
 228 12:     fstp    %st(0)          // pop y
 229         fldl    MO(one)         // 1
 230         fldt    8(%rsp)         // x : 1
 231         fabs                    // abs(x) : 1
 232         fucompp                 // < 1, == 1, or > 1
 233         fnstsw
 234         andb    $0x45, %ah
 235         cmpb    $0x45, %ah
 236         je      13f             // jump if x is NaN
 237
 238         cmpb    $0x40, %ah
 239         je      14f             // jump if |x| == 1
 240
 241         shlb    $1, %ah
 242         xorb    %ah, %dl
 243         andl    $2, %edx
 244 #ifdef PIC
 245         lea     inf_zero(%rip),%rcx
 246         fldl    (%rcx, %rdx, 4)
 247 #else
 248         fldl    inf_zero(,%rdx, 4)
 249 #endif
 250         ret
 251
 252         .align ALIGNARG(4)
 253 14:     fldl    MO(one)
 254         ret
 255
 256         .align ALIGNARG(4)
 257 13:     fldt    8(%rsp)         // load x == NaN
 258         ret
 259
 260         .align ALIGNARG(4)
 261         // x is ±inf
 262 15:     fstp    %st(0)          // y
 263         testb   $2, %dh
 264         jz      16f             // jump if x == +inf
 265
 266         // fistpll raises invalid exception for |y| >= 1L<<63, but y
 267         // may be odd unless we know |y| >= 1L<<64.
 268         fldl    MO(p64)         // 1L<<64 : y
 269         fld     %st(1)          // y : 1L<<64 : y
 270         fabs                    // |y| : 1L<<64 : y
 271         fcomip  %st(1), %st     // 1L<<64 : y
 272         fstp    %st(0)          // y
 273         jnc     16f
 274         fldl    MO(p63)         // p63 : y
 275         fxch                    // y : p63
 276         fprem                   // y%p63 : p63
 277         fstp    %st(1)          // y%p63
 278
 279         // We must find out whether y is an odd integer.
 280         fld     %st             // y : y
 281         fistpll -8(%rsp)        // y
 282         fildll  -8(%rsp)        // int(y) : y
 283         fucomip %st(1),%st
 284         ffreep  %st             // <empty>
 285         jne     17f
 286
 287         // OK, the value is an integer, but is it odd?
 288         mov     -8(%rsp), %eax
 289         mov     -4(%rsp), %edx
 290         andb    $1, %al
 291         jz      18f             // jump if not odd
 292         // It's an odd integer.
 293         shrl    $31, %edx
 294 #ifdef PIC
 295         lea     minf_mzero(%rip),%rcx
 296         fldl    (%rcx, %rdx, 8)
 297 #else
 298         fldl    minf_mzero(,%rdx, 8)
 299 #endif
 300         ret
 301
 302         .align ALIGNARG(4)
 303 16:     fcompl  MO(zero)
 304         fnstsw
 305         shrl    $5, %eax
 306         andl    $8, %eax
 307 #ifdef PIC
 308         lea     inf_zero(%rip),%rcx
 309         fldl    (%rcx, %rax, 1)
 310 #else
 311         fldl    inf_zero(,%rax, 1)
 312 #endif
 313         ret
 314
 315         .align ALIGNARG(4)
 316 17:     shll    $30, %edx       // sign bit for y in right position
 317 18:     shrl    $31, %edx
 318 #ifdef PIC
 319         lea     inf_zero(%rip),%rcx
 320         fldl    (%rcx, %rdx, 8)
 321 #else
 322         fldl    inf_zero(,%rdx, 8)
 323 #endif
 324         ret
 325
 326         .align ALIGNARG(4)
 327         // x is ±0
 328 20:     fstp    %st(0)          // y
 329         testb   $2, %dl
 330         jz      21f             // y > 0
 331
 332         // x is ±0 and y is < 0.  We must find out whether y is an odd integer.
 333         testb   $2, %dh
 334         jz      25f
 335
 336         // fistpll raises invalid exception for |y| >= 1L<<63, but y
 337         // may be odd unless we know |y| >= 1L<<64.
 338         fldl    MO(p64)         // 1L<<64 : y
 339         fld     %st(1)          // y : 1L<<64 : y
 340         fabs                    // |y| : 1L<<64 : y
 341         fcomip  %st(1), %st     // 1L<<64 : y
 342         fstp    %st(0)          // y
 343         jnc     25f
 344         fldl    MO(p63)         // p63 : y
 345         fxch                    // y : p63
 346         fprem                   // y%p63 : p63
 347         fstp    %st(1)          // y%p63
 348
 349         fld     %st             // y : y
 350         fistpll -8(%rsp)        // y
 351         fildll  -8(%rsp)        // int(y) : y
 352         fucomip %st(1),%st
 353         ffreep  %st             // <empty>
 354         jne     26f
 355
 356         // OK, the value is an integer, but is it odd?
 357         mov     -8(%rsp),%eax
 358         mov     -4(%rsp),%edx
 359         andb    $1, %al
 360         jz      27f             // jump if not odd
 361         // It's an odd integer.
 362         // Raise divide-by-zero exception and get minus infinity value.
 363         fldl    MO(one)
 364         fdivl   MO(zero)
 365         fchs
 366         ret
 367
 368 25:     fstp    %st(0)
 369 26:
 370 27:     // Raise divide-by-zero exception and get infinity value.
 371         fldl    MO(one)
 372         fdivl   MO(zero)
 373         ret
 374
 375         .align ALIGNARG(4)
 376         // x is ±0 and y is > 0.  We must find out whether y is an odd integer.
 377 21:     testb   $2, %dh
 378         jz      22f
 379
 380         // fistpll raises invalid exception for |y| >= 1L<<63, but y
 381         // may be odd unless we know |y| >= 1L<<64.
 382         fldl    MO(p64)         // 1L<<64 : y
 383         fxch                    // y : 1L<<64
 384         fcomi   %st(1), %st     // y : 1L<<64
 385         fstp    %st(1)          // y
 386         jnc     22f
 387         fldl    MO(p63)         // p63 : y
 388         fxch                    // y : p63
 389         fprem                   // y%p63 : p63
 390         fstp    %st(1)          // y%p63
 391
 392         fld     %st             // y : y
 393         fistpll -8(%rsp)        // y
 394         fildll  -8(%rsp)        // int(y) : y
 395         fucomip %st(1),%st
 396         ffreep  %st             // <empty>
 397         jne     23f
 398
 399         // OK, the value is an integer, but is it odd?
 400         mov     -8(%rsp),%eax
 401         mov     -4(%rsp),%edx
 402         andb    $1, %al
 403         jz      24f             // jump if not odd
 404         // It's an odd integer.
 405         fldl    MO(mzero)
 406         ret
 407
 408 22:     fstp    %st(0)
 409 23:
 410 24:     fldl    MO(zero)
 411         ret
 412
 413 END(__ieee754_powl)
 414 strong_alias (__ieee754_powl, __powl_finite)