sysdeps/x86_64/fpu/e_powl.S

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