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