sysdeps/x86_64/fpu/e_powl.S

   1 /* ix87 specific implementation of pow function.
   2    Copyright (C) 1996, 1997, 1998, 1999, 2001, 2004, 2007
   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, write to the Free
  19    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
  20    02111-1307 USA.  */
  21
  22 #include <machine/asm.h>
  23
  24 #ifdef __ELF__
  25         .section .rodata
  26 #else
  27         .text
  28 #endif
  29
  30         .align ALIGNARG(4)
  31         ASM_TYPE_DIRECTIVE(infinity,@object)
  32 inf_zero:
  33 infinity:
  34         .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
  35         ASM_SIZE_DIRECTIVE(infinity)
  36         ASM_TYPE_DIRECTIVE(zero,@object)
  37 zero:   .double 0.0
  38         ASM_SIZE_DIRECTIVE(zero)
  39         ASM_TYPE_DIRECTIVE(minf_mzero,@object)
  40 minf_mzero:
  41 minfinity:
  42         .byte 0, 0, 0, 0, 0, 0, 0xf0, 0xff
  43 mzero:
  44         .byte 0, 0, 0, 0, 0, 0, 0, 0x80
  45         ASM_SIZE_DIRECTIVE(minf_mzero)
  46         ASM_TYPE_DIRECTIVE(one,@object)
  47 one:    .double 1.0
  48         ASM_SIZE_DIRECTIVE(one)
  49         ASM_TYPE_DIRECTIVE(limit,@object)
  50 limit:  .double 0.29
  51         ASM_SIZE_DIRECTIVE(limit)
  52         ASM_TYPE_DIRECTIVE(p63,@object)
  53 p63:
  54         .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43
  55         ASM_SIZE_DIRECTIVE(p63)
  56
  57 #ifdef PIC
  58 #define MO(op) op##(%rip)
  59 #else
  60 #define MO(op) op
  61 #endif
  62
  63         .text
  64 ENTRY(__ieee754_powl)
  65         fldt    24(%rsp)        // y
  66         fxam
  67
  68
  69         fnstsw
  70         movb    %ah, %dl
  71         andb    $0x45, %ah
  72         cmpb    $0x40, %ah      // is y == 0 ?
  73         je      11f
  74
  75         cmpb    $0x05, %ah      // is y == ±inf ?
  76         je      12f
  77
  78         cmpb    $0x01, %ah      // is y == NaN ?
  79         je      30f
  80
  81         fldt    8(%rsp)         // x : y
  82
  83         fxam
  84         fnstsw
  85         movb    %ah, %dh
  86         andb    $0x45, %ah
  87         cmpb    $0x40, %ah
  88         je      20f             // x is ±0
  89
  90         cmpb    $0x05, %ah
  91         je      15f             // x is ±inf
  92
  93         fxch                    // y : x
  94
  95         /* fistpll raises invalid exception for |y| >= 1L<<63.  */
  96         fldl    MO(p63)         // 1L<<63 : y : x
  97         fld     %st(1)          // y : 1L<<63 : y : x
  98         fabs                    // |y| : 1L<<63 : y : x
  99         fcomip  %st(1), %st     // 1L<<63 : y : x
 100         fstp    %st(0)          // y : x
 101         jnc     2f
 102
 103         /* First see whether `y' is a natural number.  In this case we
 104            can use a more precise algorithm.  */
 105         fld     %st             // y : y : x
 106         fistpll -8(%rsp)        // y : x
 107         fildll  -8(%rsp)        // int(y) : y : x
 108         fucomip %st(1),%st      // y : x
 109         jne     2f
 110
 111         /* OK, we have an integer value for y.  */
 112         mov     -8(%rsp),%eax
 113         mov     -4(%rsp),%edx
 114         orl     $0, %edx
 115         fstp    %st(0)          // x
 116         jns     4f              // y >= 0, jump
 117         fdivrl  MO(one)         // 1/x          (now referred to as x)
 118         negl    %eax
 119         adcl    $0, %edx
 120         negl    %edx
 121 4:      fldl    MO(one)         // 1 : x
 122         fxch
 123
 124 6:      shrdl   $1, %edx, %eax
 125         jnc     5f
 126         fxch
 127         fmul    %st(1)          // x : ST*x
 128         fxch
 129 5:      fmul    %st(0), %st     // x*x : ST*x
 130         shrl    $1, %edx
 131         movl    %eax, %ecx
 132         orl     %edx, %ecx
 133         jnz     6b
 134         fstp    %st(0)          // ST*x
 135         ret
 136
 137         /* y is ±NAN */
 138 30:     fldt    8(%rsp)         // x : y
 139         fldl    MO(one)         // 1.0 : x : y
 140         fucomip %st(1),%st      // x : y
 141         je      31f
 142         fxch                    // y : x
 143 31:     fstp    %st(1)
 144         ret
 145
 146         .align ALIGNARG(4)
 147 2:      /* y is a real number.  */
 148         fxch                    // x : y
 149         fldl    MO(one)         // 1.0 : x : y
 150         fldl    MO(limit)       // 0.29 : 1.0 : x : y
 151         fld     %st(2)          // x : 0.29 : 1.0 : x : y
 152         fsub    %st(2)          // x-1 : 0.29 : 1.0 : x : y
 153         fabs                    // |x-1| : 0.29 : 1.0 : x : y
 154         fucompp                 // 1.0 : x : y
 155         fnstsw
 156         fxch                    // x : 1.0 : y
 157         test    $4500,%eax
 158         jz      7f
 159         fsub    %st(1)          // x-1 : 1.0 : y
 160         fyl2xp1                 // log2(x) : y
 161         jmp     8f
 162
 163 7:      fyl2x                   // log2(x) : y
 164 8:      fmul    %st(1)          // y*log2(x) : y
 165         fxam
 166         fnstsw
 167         andb    $0x45, %ah
 168         cmpb    $0x05, %ah      // is y*log2(x) == ±inf ?
 169         je      28f
 170         fst     %st(1)          // y*log2(x) : y*log2(x)
 171         frndint                 // int(y*log2(x)) : y*log2(x)
 172         fsubr   %st, %st(1)     // int(y*log2(x)) : fract(y*log2(x))
 173         fxch                    // fract(y*log2(x)) : int(y*log2(x))
 174         f2xm1                   // 2^fract(y*log2(x))-1 : int(y*log2(x))
 175         faddl   MO(one)         // 2^fract(y*log2(x)) : int(y*log2(x))
 176         fscale                  // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
 177         fstp    %st(1)          // 2^fract(y*log2(x))*2^int(y*log2(x))
 178         ret
 179
 180 28:     fstp    %st(1)          // y*log2(x)
 181         fldl    MO(one)         // 1 : y*log2(x)
 182         fscale                  // 2^(y*log2(x)) : y*log2(x)
 183         fstp    %st(1)          // 2^(y*log2(x))
 184         ret
 185
 186         // pow(x,±0) = 1
 187         .align ALIGNARG(4)
 188 11:     fstp    %st(0)          // pop y
 189         fldl    MO(one)
 190         ret
 191
 192         // y == ±inf
 193         .align ALIGNARG(4)
 194 12:     fstp    %st(0)          // pop y
 195         fldl    MO(one)         // 1
 196         fldt    8(%rsp)         // x : 1
 197         fabs                    // abs(x) : 1
 198         fucompp                 // < 1, == 1, or > 1
 199         fnstsw
 200         andb    $0x45, %ah
 201         cmpb    $0x45, %ah
 202         je      13f             // jump if x is NaN
 203
 204         cmpb    $0x40, %ah
 205         je      14f             // jump if |x| == 1
 206
 207         shlb    $1, %ah
 208         xorb    %ah, %dl
 209         andl    $2, %edx
 210 #ifdef PIC
 211         lea     inf_zero(%rip),%rcx
 212         fldl    (%rcx, %rdx, 4)
 213 #else
 214         fldl    inf_zero(,%rdx, 4)
 215 #endif
 216         ret
 217
 218         .align ALIGNARG(4)
 219 14:     fldl    MO(one)
 220         ret
 221
 222         .align ALIGNARG(4)
 223 13:     fldt    8(%rsp)         // load x == NaN
 224         ret
 225
 226         .align ALIGNARG(4)
 227         // x is ±inf
 228 15:     fstp    %st(0)          // y
 229         testb   $2, %dh
 230         jz      16f             // jump if x == +inf
 231
 232         // We must find out whether y is an odd integer.
 233         fld     %st             // y : y
 234         fistpll -8(%rsp)        // y
 235         fildll  -8(%rsp)        // int(y) : y
 236         fucomip %st(1),%st
 237         ffreep  %st             // <empty>
 238         jne     17f
 239
 240         // OK, the value is an integer, but is it odd?
 241         mov     -8(%rsp), %eax
 242         mov     -4(%rsp), %edx
 243         andb    $1, %al
 244         jz      18f             // jump if not odd
 245         // It's an odd integer.
 246         shrl    $31, %edx
 247 #ifdef PIC
 248         lea     minf_mzero(%rip),%rcx
 249         fldl    (%rcx, %rdx, 8)
 250 #else
 251         fldl    minf_mzero(,%rdx, 8)
 252 #endif
 253         ret
 254
 255         .align ALIGNARG(4)
 256 16:     fcompl  MO(zero)
 257         fnstsw
 258         shrl    $5, %eax
 259         andl    $8, %eax
 260 #ifdef PIC
 261         lea     inf_zero(%rip),%rcx
 262         fldl    (%rcx, %rax, 1)
 263 #else
 264         fldl    inf_zero(,%rax, 1)
 265 #endif
 266         ret
 267
 268         .align ALIGNARG(4)
 269 17:     shll    $30, %edx       // sign bit for y in right position
 270 18:     shrl    $31, %edx
 271 #ifdef PIC
 272         lea     inf_zero(%rip),%rcx
 273         fldl    (%rcx, %rdx, 8)
 274 #else
 275         fldl    inf_zero(,%rdx, 8)
 276 #endif
 277         ret
 278
 279         .align ALIGNARG(4)
 280         // x is ±0
 281 20:     fstp    %st(0)          // y
 282         testb   $2, %dl
 283         jz      21f             // y > 0
 284
 285         // x is ±0 and y is < 0.  We must find out whether y is an odd integer.
 286         testb   $2, %dh
 287         jz      25f
 288
 289         fld     %st             // y : y
 290         fistpll -8(%rsp)        // y
 291         fildll  -8(%rsp)        // int(y) : y
 292         fucomip %st(1),%st
 293         ffreep  %st             // <empty>
 294         jne     26f
 295
 296         // OK, the value is an integer, but is it odd?
 297         mov     -8(%rsp),%eax
 298         mov     -4(%rsp),%edx
 299         andb    $1, %al
 300         jz      27f             // jump if not odd
 301         // It's an odd integer.
 302         // Raise divide-by-zero exception and get minus infinity value.
 303         fldl    MO(one)
 304         fdivl   MO(zero)
 305         fchs
 306         ret
 307
 308 25:     fstp    %st(0)
 309 26:
 310 27:     // Raise divide-by-zero exception and get infinity value.
 311         fldl    MO(one)
 312         fdivl   MO(zero)
 313         ret
 314
 315         .align ALIGNARG(4)
 316         // x is ±0 and y is > 0.  We must find out whether y is an odd integer.
 317 21:     testb   $2, %dh
 318         jz      22f
 319
 320         fld     %st             // y : y
 321         fistpll -8(%rsp)        // y
 322         fildll  -8(%rsp)        // int(y) : y
 323         fucomip %st(1),%st
 324         ffreep  %st             // <empty>
 325         jne     23f
 326
 327         // OK, the value is an integer, but is it odd?
 328         mov     -8(%rsp),%eax
 329         mov     -4(%rsp),%edx
 330         andb    $1, %al
 331         jz      24f             // jump if not odd
 332         // It's an odd integer.
 333         fldl    MO(mzero)
 334         ret
 335
 336 22:     fstp    %st(0)
 337 23:
 338 24:     fldl    MO(zero)
 339         ret
 340
 341 END(__ieee754_powl)