sysdeps/x86_64/fpu/e_powl.S

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