sysdeps/x86_64/fpu/e_powl.S

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