sysdeps/i386/fpu/e_powl.S

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