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