sysdeps/i386/fpu/e_powl.S

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