sysdeps/i386/fpu/e_pow.S

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