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