sysdeps/i386/fpu/e_powf.S

   1 /* ix87 specific implementation of pow function.
   2    Copyright (C) 1996-2017 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, see
  18    <http://www.gnu.org/licenses/>.  */
  19
  20 #include <machine/asm.h>
  21 #include <i386-math-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 p31,@object
  33 p31:    .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x41
  34         ASM_SIZE_DIRECTIVE(p31)
  35
  36         .section .rodata.cst16,"aM",@progbits,16
  37
  38         .p2align 3
  39         .type infinity,@object
  40 inf_zero:
  41 infinity:
  42         .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
  43         ASM_SIZE_DIRECTIVE(infinity)
  44         .type zero,@object
  45 zero:   .double 0.0
  46         ASM_SIZE_DIRECTIVE(zero)
  47         .type minf_mzero,@object
  48 minf_mzero:
  49 minfinity:
  50         .byte 0, 0, 0, 0, 0, 0, 0xf0, 0xff
  51 mzero:
  52         .byte 0, 0, 0, 0, 0, 0, 0, 0x80
  53         ASM_SIZE_DIRECTIVE(minf_mzero)
  54 DEFINE_FLT_MIN
  55
  56 #ifdef PIC
  57 # define MO(op) op##@GOTOFF(%ecx)
  58 # define MOX(op,x,f) op##@GOTOFF(%ecx,x,f)
  59 #else
  60 # define MO(op) op
  61 # define MOX(op,x,f) op(,x,f)
  62 #endif
  63
  64         .text
  65 ENTRY(__ieee754_powf)
  66         flds    8(%esp) // y
  67         fxam
  68
  69 #ifdef  PIC
  70         LOAD_PIC_REG (cx)
  71 #endif
  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         flds    4(%esp)         // x : y
  86
  87         subl    $4, %esp
  88         cfi_adjust_cfa_offset (4)
  89
  90         fxam
  91         fnstsw
  92         movb    %ah, %dh
  93         andb    $0x45, %ah
  94         cmpb    $0x40, %ah
  95         je      20f             // x is ±0
  96
  97         cmpb    $0x05, %ah
  98         je      15f             // x is ±inf
  99
 100         cmpb    $0x01, %ah
 101         je      33f             // x is NaN
 102
 103         fxch                    // y : x
 104
 105         /* fistpl raises invalid exception for |y| >= 1L<<31.  */
 106         fld     %st             // y : y : x
 107         fabs                    // |y| : y : x
 108         fcompl  MO(p31)         // y : x
 109         fnstsw
 110         sahf
 111         jnc     2f
 112
 113         /* First see whether `y' is a natural number.  In this case we
 114            can use a more precise algorithm.  */
 115         fld     %st             // y : y : x
 116         fistpl  (%esp)          // y : x
 117         fildl   (%esp)          // int(y) : y : x
 118         fucomp  %st(1)          // y : x
 119         fnstsw
 120         sahf
 121         jne     3f
 122
 123         /* OK, we have an integer value for y.  */
 124         popl    %edx
 125         cfi_adjust_cfa_offset (-4)
 126         orl     $0, %edx
 127         fstp    %st(0)          // x
 128         jns     4f              // y >= 0, jump
 129         fdivrl  MO(one)         // 1/x          (now referred to as x)
 130         negl    %edx
 131 4:      fldl    MO(one)         // 1 : x
 132         fxch
 133
 134         /* If y is even, take the absolute value of x.  Otherwise,
 135            ensure all intermediate values that might overflow have the
 136            sign of x.  */
 137         testb   $1, %dl
 138         jnz     6f
 139         fabs
 140
 141 6:      shrl    $1, %edx
 142         jnc     5f
 143         fxch
 144         fabs
 145         fmul    %st(1)          // x : ST*x
 146         fxch
 147 5:      fld     %st             // x : x : ST*x
 148         fabs                    // |x| : x : ST*x
 149         fmulp                   // |x|*x : ST*x
 150         testl   %edx, %edx
 151         jnz     6b
 152         fstp    %st(0)          // ST*x
 153         FLT_NARROW_EVAL_UFLOW_NONNAN
 154         ret
 155
 156         /* y is ±NAN */
 157 30:     flds    4(%esp)         // x : y
 158         fldl    MO(one)         // 1.0 : x : y
 159         fucomp  %st(1)          // x : y
 160         fnstsw
 161         sahf
 162         je      31f
 163         fxch                    // y : x
 164 31:     fstp    %st(1)
 165         ret
 166
 167         cfi_adjust_cfa_offset (4)
 168         .align ALIGNARG(4)
 169 2:      /* y is a large integer (so even).  */
 170         fxch                    // x : y
 171         fabs                    // |x| : y
 172         fxch                    // y : x
 173         .align ALIGNARG(4)
 174 3:      /* y is a real number.  */
 175         fxch                    // x : y
 176         fldl    MO(one)         // 1.0 : x : y
 177         fldl    MO(limit)       // 0.29 : 1.0 : x : y
 178         fld     %st(2)          // x : 0.29 : 1.0 : x : y
 179         fsub    %st(2)          // x-1 : 0.29 : 1.0 : x : y
 180         fabs                    // |x-1| : 0.29 : 1.0 : x : y
 181         fucompp                 // 1.0 : x : y
 182         fnstsw
 183         fxch                    // x : 1.0 : y
 184         sahf
 185         ja      7f
 186         fsub    %st(1)          // x-1 : 1.0 : y
 187         fyl2xp1                 // log2(x) : y
 188         jmp     8f
 189
 190 7:      fyl2x                   // log2(x) : y
 191 8:      fmul    %st(1)          // y*log2(x) : y
 192         fst     %st(1)          // y*log2(x) : y*log2(x)
 193         frndint                 // int(y*log2(x)) : y*log2(x)
 194         fsubr   %st, %st(1)     // int(y*log2(x)) : fract(y*log2(x))
 195         fxch                    // fract(y*log2(x)) : int(y*log2(x))
 196         f2xm1                   // 2^fract(y*log2(x))-1 : int(y*log2(x))
 197         faddl   MO(one)         // 2^fract(y*log2(x)) : int(y*log2(x))
 198         fscale                  // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
 199 32:     addl    $4, %esp
 200         cfi_adjust_cfa_offset (-4)
 201         fstp    %st(1)          // 2^fract(y*log2(x))*2^int(y*log2(x))
 202         FLT_NARROW_EVAL_UFLOW_NONNAN
 203         ret
 204
 205         /* x is NaN.  */
 206         cfi_adjust_cfa_offset (4)
 207 33:     addl    $4, %esp
 208         cfi_adjust_cfa_offset (-4)
 209         fstp    %st(1)
 210         ret
 211
 212         // pow(x,±0) = 1
 213         .align ALIGNARG(4)
 214 11:     fstp    %st(0)          // pop y
 215         fldl    MO(one)
 216         ret
 217
 218         // y == ±inf
 219         .align ALIGNARG(4)
 220 12:     fstp    %st(0)          // pop y
 221         fldl    MO(one)         // 1
 222         flds    4(%esp)         // x : 1
 223         fabs                    // abs(x) : 1
 224         fucompp                 // < 1, == 1, or > 1
 225         fnstsw
 226         andb    $0x45, %ah
 227         cmpb    $0x45, %ah
 228         je      13f             // jump if x is NaN
 229
 230         cmpb    $0x40, %ah
 231         je      14f             // jump if |x| == 1
 232
 233         shlb    $1, %ah
 234         xorb    %ah, %dl
 235         andl    $2, %edx
 236         fldl    MOX(inf_zero, %edx, 4)
 237         ret
 238
 239         .align ALIGNARG(4)
 240 14:     fldl    MO(one)
 241         ret
 242
 243         .align ALIGNARG(4)
 244 13:     flds    4(%esp)         // load x == NaN
 245         ret
 246
 247         cfi_adjust_cfa_offset (4)
 248         .align ALIGNARG(4)
 249         // x is ±inf
 250 15:     fstp    %st(0)          // y
 251         testb   $2, %dh
 252         jz      16f             // jump if x == +inf
 253
 254         // fistpl raises invalid exception for |y| >= 1L<<31, so test
 255         // that (in which case y is certainly even) before testing
 256         // whether y is odd.
 257         fld     %st             // y : y
 258         fabs                    // |y| : y
 259         fcompl  MO(p31)         // y
 260         fnstsw
 261         sahf
 262         jnc     16f
 263
 264         // We must find out whether y is an odd integer.
 265         fld     %st             // y : y
 266         fistpl  (%esp)          // y
 267         fildl   (%esp)          // int(y) : y
 268         fucompp                 // <empty>
 269         fnstsw
 270         sahf
 271         jne     17f
 272
 273         // OK, the value is an integer.
 274         popl    %edx
 275         cfi_adjust_cfa_offset (-4)
 276         testb   $1, %dl
 277         jz      18f             // jump if not odd
 278         // It's an odd integer.
 279         shrl    $31, %edx
 280         fldl    MOX(minf_mzero, %edx, 8)
 281         ret
 282
 283         cfi_adjust_cfa_offset (4)
 284         .align ALIGNARG(4)
 285 16:     fcompl  MO(zero)
 286         addl    $4, %esp
 287         cfi_adjust_cfa_offset (-4)
 288         fnstsw
 289         shrl    $5, %eax
 290         andl    $8, %eax
 291         fldl    MOX(inf_zero, %eax, 1)
 292         ret
 293
 294         cfi_adjust_cfa_offset (4)
 295         .align ALIGNARG(4)
 296 17:     shll    $30, %edx       // sign bit for y in right position
 297         addl    $4, %esp
 298         cfi_adjust_cfa_offset (-4)
 299 18:     shrl    $31, %edx
 300         fldl    MOX(inf_zero, %edx, 8)
 301         ret
 302
 303         cfi_adjust_cfa_offset (4)
 304         .align ALIGNARG(4)
 305         // x is ±0
 306 20:     fstp    %st(0)          // y
 307         testb   $2, %dl
 308         jz      21f             // y > 0
 309
 310         // x is ±0 and y is < 0.  We must find out whether y is an odd integer.
 311         testb   $2, %dh
 312         jz      25f
 313
 314         // fistpl raises invalid exception for |y| >= 1L<<31, so test
 315         // that (in which case y is certainly even) before testing
 316         // whether y is odd.
 317         fld     %st             // y : y
 318         fabs                    // |y| : y
 319         fcompl  MO(p31)         // y
 320         fnstsw
 321         sahf
 322         jnc     25f
 323
 324         fld     %st             // y : y
 325         fistpl  (%esp)          // y
 326         fildl   (%esp)          // int(y) : y
 327         fucompp                 // <empty>
 328         fnstsw
 329         sahf
 330         jne     26f
 331
 332         // OK, the value is an integer.
 333         popl    %edx
 334         cfi_adjust_cfa_offset (-4)
 335         testb   $1, %dl
 336         jz      27f             // jump if not odd
 337         // It's an odd integer.
 338         // Raise divide-by-zero exception and get minus infinity value.
 339         fldl    MO(one)
 340         fdivl   MO(zero)
 341         fchs
 342         ret
 343
 344         cfi_adjust_cfa_offset (4)
 345 25:     fstp    %st(0)
 346 26:     addl    $4, %esp
 347         cfi_adjust_cfa_offset (-4)
 348 27:     // Raise divide-by-zero exception and get infinity value.
 349         fldl    MO(one)
 350         fdivl   MO(zero)
 351         ret
 352
 353         cfi_adjust_cfa_offset (4)
 354         .align ALIGNARG(4)
 355         // x is ±0 and y is > 0.  We must find out whether y is an odd integer.
 356 21:     testb   $2, %dh
 357         jz      22f
 358
 359         // fistpl raises invalid exception for |y| >= 1L<<31, so test
 360         // that (in which case y is certainly even) before testing
 361         // whether y is odd.
 362         fcoml   MO(p31)         // y
 363         fnstsw
 364         sahf
 365         jnc     22f
 366
 367         fld     %st             // y : y
 368         fistpl  (%esp)          // y
 369         fildl   (%esp)          // int(y) : y
 370         fucompp                 // <empty>
 371         fnstsw
 372         sahf
 373         jne     23f
 374
 375         // OK, the value is an integer.
 376         popl    %edx
 377         cfi_adjust_cfa_offset (-4)
 378         testb   $1, %dl
 379         jz      24f             // jump if not odd
 380         // It's an odd integer.
 381         fldl    MO(mzero)
 382         ret
 383
 384         cfi_adjust_cfa_offset (4)
 385 22:     fstp    %st(0)
 386 23:     addl    $4, %esp        // Don't use pop.
 387         cfi_adjust_cfa_offset (-4)
 388 24:     fldl    MO(zero)
 389         ret
 390
 391 END(__ieee754_powf)
 392 strong_alias (__ieee754_powf, __powf_finite)