sysdeps/i386/fpu/e_powl.S

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