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     2f
 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 real number.  */
 149         fxch                    // x : y
 150         fldl    MO(one)         // 1.0 : x : y
 151         fldl    MO(limit)       // 0.29 : 1.0 : x : y
 152         fld     %st(2)          // x : 0.29 : 1.0 : x : y
 153         fsub    %st(2)          // x-1 : 0.29 : 1.0 : x : y
 154         fabs                    // |x-1| : 0.29 : 1.0 : x : y
 155         fucompp                 // 1.0 : x : y
 156         fnstsw
 157         fxch                    // x : 1.0 : y
 158         test    $0x4500,%eax
 159         jz      7f
 160         fsub    %st(1)          // x-1 : 1.0 : y
 161         fyl2xp1                 // log2(x) : y
 162         jmp     8f
 163
 164 7:      fyl2x                   // log2(x) : y
 165 8:      fmul    %st(1)          // y*log2(x) : y
 166         fxam
 167         fnstsw
 168         andb    $0x45, %ah
 169         cmpb    $0x05, %ah      // is y*log2(x) == ±inf ?
 170         je      28f
 171         fst     %st(1)          // y*log2(x) : y*log2(x)
 172         frndint                 // int(y*log2(x)) : y*log2(x)
 173         fsubr   %st, %st(1)     // int(y*log2(x)) : fract(y*log2(x))
 174         fxch                    // fract(y*log2(x)) : int(y*log2(x))
 175         f2xm1                   // 2^fract(y*log2(x))-1 : int(y*log2(x))
 176         faddl   MO(one)         // 2^fract(y*log2(x)) : int(y*log2(x))
 177         fscale                  // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
 178         fstp    %st(1)          // 2^fract(y*log2(x))*2^int(y*log2(x))
 179         ret
 180
 181 28:     fstp    %st(1)          // y*log2(x)
 182         fldl    MO(one)         // 1 : y*log2(x)
 183         fscale                  // 2^(y*log2(x)) : y*log2(x)
 184         fstp    %st(1)          // 2^(y*log2(x))
 185         ret
 186
 187         // pow(x,±0) = 1
 188         .align ALIGNARG(4)
 189 11:     fstp    %st(0)          // pop y
 190         fldl    MO(one)
 191         ret
 192
 193         // y == ±inf
 194         .align ALIGNARG(4)
 195 12:     fstp    %st(0)          // pop y
 196         fldl    MO(one)         // 1
 197         fldt    8(%rsp)         // x : 1
 198         fabs                    // abs(x) : 1
 199         fucompp                 // < 1, == 1, or > 1
 200         fnstsw
 201         andb    $0x45, %ah
 202         cmpb    $0x45, %ah
 203         je      13f             // jump if x is NaN
 204
 205         cmpb    $0x40, %ah
 206         je      14f             // jump if |x| == 1
 207
 208         shlb    $1, %ah
 209         xorb    %ah, %dl
 210         andl    $2, %edx
 211 #ifdef PIC
 212         lea     inf_zero(%rip),%rcx
 213         fldl    (%rcx, %rdx, 4)
 214 #else
 215         fldl    inf_zero(,%rdx, 4)
 216 #endif
 217         ret
 218
 219         .align ALIGNARG(4)
 220 14:     fldl    MO(one)
 221         ret
 222
 223         .align ALIGNARG(4)
 224 13:     fldt    8(%rsp)         // load x == NaN
 225         ret
 226
 227         .align ALIGNARG(4)
 228         // x is ±inf
 229 15:     fstp    %st(0)          // y
 230         testb   $2, %dh
 231         jz      16f             // jump if x == +inf
 232
 233         // fistpll raises invalid exception for |y| >= 1L<<63, but y
 234         // may be odd unless we know |y| >= 1L<<64.
 235         fldl    MO(p64)         // 1L<<64 : y
 236         fld     %st(1)          // y : 1L<<64 : y
 237         fabs                    // |y| : 1L<<64 : y
 238         fcomip  %st(1), %st     // 1L<<64 : y
 239         fstp    %st(0)          // y
 240         jnc     16f
 241         fldl    MO(p63)         // p63 : y
 242         fxch                    // y : p63
 243         fprem                   // y%p63 : p63
 244         fstp    %st(1)          // y%p63
 245
 246         // We must find out whether y is an odd integer.
 247         fld     %st             // y : y
 248         fistpll -8(%rsp)        // y
 249         fildll  -8(%rsp)        // int(y) : y
 250         fucomip %st(1),%st
 251         ffreep  %st             // <empty>
 252         jne     17f
 253
 254         // OK, the value is an integer, but is it odd?
 255         mov     -8(%rsp), %eax
 256         mov     -4(%rsp), %edx
 257         andb    $1, %al
 258         jz      18f             // jump if not odd
 259         // It's an odd integer.
 260         shrl    $31, %edx
 261 #ifdef PIC
 262         lea     minf_mzero(%rip),%rcx
 263         fldl    (%rcx, %rdx, 8)
 264 #else
 265         fldl    minf_mzero(,%rdx, 8)
 266 #endif
 267         ret
 268
 269         .align ALIGNARG(4)
 270 16:     fcompl  MO(zero)
 271         fnstsw
 272         shrl    $5, %eax
 273         andl    $8, %eax
 274 #ifdef PIC
 275         lea     inf_zero(%rip),%rcx
 276         fldl    (%rcx, %rax, 1)
 277 #else
 278         fldl    inf_zero(,%rax, 1)
 279 #endif
 280         ret
 281
 282         .align ALIGNARG(4)
 283 17:     shll    $30, %edx       // sign bit for y in right position
 284 18:     shrl    $31, %edx
 285 #ifdef PIC
 286         lea     inf_zero(%rip),%rcx
 287         fldl    (%rcx, %rdx, 8)
 288 #else
 289         fldl    inf_zero(,%rdx, 8)
 290 #endif
 291         ret
 292
 293         .align ALIGNARG(4)
 294         // x is ±0
 295 20:     fstp    %st(0)          // y
 296         testb   $2, %dl
 297         jz      21f             // y > 0
 298
 299         // x is ±0 and y is < 0.  We must find out whether y is an odd integer.
 300         testb   $2, %dh
 301         jz      25f
 302
 303         // fistpll raises invalid exception for |y| >= 1L<<63, but y
 304         // may be odd unless we know |y| >= 1L<<64.
 305         fldl    MO(p64)         // 1L<<64 : y
 306         fld     %st(1)          // y : 1L<<64 : y
 307         fabs                    // |y| : 1L<<64 : y
 308         fcomip  %st(1), %st     // 1L<<64 : y
 309         fstp    %st(0)          // y
 310         jnc     25f
 311         fldl    MO(p63)         // p63 : y
 312         fxch                    // y : p63
 313         fprem                   // y%p63 : p63
 314         fstp    %st(1)          // y%p63
 315
 316         fld     %st             // y : y
 317         fistpll -8(%rsp)        // y
 318         fildll  -8(%rsp)        // int(y) : y
 319         fucomip %st(1),%st
 320         ffreep  %st             // <empty>
 321         jne     26f
 322
 323         // OK, the value is an integer, but is it odd?
 324         mov     -8(%rsp),%eax
 325         mov     -4(%rsp),%edx
 326         andb    $1, %al
 327         jz      27f             // jump if not odd
 328         // It's an odd integer.
 329         // Raise divide-by-zero exception and get minus infinity value.
 330         fldl    MO(one)
 331         fdivl   MO(zero)
 332         fchs
 333         ret
 334
 335 25:     fstp    %st(0)
 336 26:
 337 27:     // Raise divide-by-zero exception and get infinity value.
 338         fldl    MO(one)
 339         fdivl   MO(zero)
 340         ret
 341
 342         .align ALIGNARG(4)
 343         // x is ±0 and y is > 0.  We must find out whether y is an odd integer.
 344 21:     testb   $2, %dh
 345         jz      22f
 346
 347         // fistpll raises invalid exception for |y| >= 1L<<63, but y
 348         // may be odd unless we know |y| >= 1L<<64.
 349         fldl    MO(p64)         // 1L<<64 : y
 350         fxch                    // y : 1L<<64
 351         fcomi   %st(1), %st     // y : 1L<<64
 352         fstp    %st(1)          // y
 353         jnc     22f
 354         fldl    MO(p63)         // p63 : y
 355         fxch                    // y : p63
 356         fprem                   // y%p63 : p63
 357         fstp    %st(1)          // y%p63
 358
 359         fld     %st             // y : y
 360         fistpll -8(%rsp)        // y
 361         fildll  -8(%rsp)        // int(y) : y
 362         fucomip %st(1),%st
 363         ffreep  %st             // <empty>
 364         jne     23f
 365
 366         // OK, the value is an integer, but is it odd?
 367         mov     -8(%rsp),%eax
 368         mov     -4(%rsp),%edx
 369         andb    $1, %al
 370         jz      24f             // jump if not odd
 371         // It's an odd integer.
 372         fldl    MO(mzero)
 373         ret
 374
 375 22:     fstp    %st(0)
 376 23:
 377 24:     fldl    MO(zero)
 378         ret
 379
 380 END(__ieee754_powl)
 381 strong_alias (__ieee754_powl, __powl_finite)