sysdeps/i386/fpu/e_pow.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 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 6:      shrdl   $1, %edx, %eax
 148         jnc     5f
 149         fxch
 150         fmul    %st(1)          // x : ST*x
 151         fxch
 152 5:      fmul    %st(0), %st     // x*x : ST*x
 153         shrl    $1, %edx
 154         movl    %eax, %ecx
 155         orl     %edx, %ecx
 156         jnz     6b
 157         fstp    %st(0)          // ST*x
 158         ret
 159
 160         /* y is ±NAN */
 161 30:     fldl    4(%esp)         // x : y
 162         fldl    MO(one)         // 1.0 : x : y
 163         fucomp  %st(1)          // x : y
 164         fnstsw
 165         sahf
 166         je      31f
 167         fxch                    // y : x
 168 31:     fstp    %st(1)
 169         ret
 170
 171         cfi_adjust_cfa_offset (8)
 172 32:     addl    $8, %esp
 173         cfi_adjust_cfa_offset (-8)
 174         fstp    %st(1)
 175         ret
 176
 177         cfi_adjust_cfa_offset (8)
 178         .align ALIGNARG(4)
 179 2:      // y is a large integer (absolute value at least 1L<<10), but
 180         // may be odd unless at least 1L<<64.  So it may be necessary
 181         // to adjust the sign of a negative result afterwards.
 182         fxch                    // x : y
 183         fabs                    // |x| : y
 184         fxch                    // y : x
 185         .align ALIGNARG(4)
 186 3:      /* y is a real number.  */
 187         fxch                    // x : y
 188         fldl    MO(one)         // 1.0 : x : y
 189         fldl    MO(limit)       // 0.29 : 1.0 : x : y
 190         fld     %st(2)          // x : 0.29 : 1.0 : x : y
 191         fsub    %st(2)          // x-1 : 0.29 : 1.0 : x : y
 192         fabs                    // |x-1| : 0.29 : 1.0 : x : y
 193         fucompp                 // 1.0 : x : y
 194         fnstsw
 195         fxch                    // x : 1.0 : y
 196         sahf
 197         ja      7f
 198         fsub    %st(1)          // x-1 : 1.0 : y
 199         fyl2xp1                 // log2(x) : y
 200         jmp     8f
 201
 202 7:      fyl2x                   // log2(x) : y
 203 8:      fmul    %st(1)          // y*log2(x) : y
 204         fst     %st(1)          // y*log2(x) : y*log2(x)
 205         frndint                 // int(y*log2(x)) : y*log2(x)
 206         fsubr   %st, %st(1)     // int(y*log2(x)) : fract(y*log2(x))
 207         fxch                    // fract(y*log2(x)) : int(y*log2(x))
 208         f2xm1                   // 2^fract(y*log2(x))-1 : int(y*log2(x))
 209         faddl   MO(one)         // 2^fract(y*log2(x)) : int(y*log2(x))
 210         fscale                  // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
 211         fstp    %st(1)          // 2^fract(y*log2(x))*2^int(y*log2(x))
 212         testb   $2, %dh
 213         jz      292f
 214         // x is negative.  If y is an odd integer, negate the result.
 215         fldl    20(%esp)        // y : abs(result)
 216         fld     %st             // y : y : abs(result)
 217         fabs                    // |y| : y : abs(result)
 218         fcompl  MO(p63)         // y : abs(result)
 219         fnstsw
 220         sahf
 221         jnc     291f
 222
 223         // We must find out whether y is an odd integer.
 224         fld     %st             // y : y : abs(result)
 225         fistpll (%esp)          // y : abs(result)
 226         fildll  (%esp)          // int(y) : y : abs(result)
 227         fucompp                 // abs(result)
 228         fnstsw
 229         sahf
 230         jne     292f
 231
 232         // OK, the value is an integer, but is it odd?
 233         popl    %eax
 234         cfi_adjust_cfa_offset (-4)
 235         popl    %edx
 236         cfi_adjust_cfa_offset (-4)
 237         andb    $1, %al
 238         jz      290f            // jump if not odd
 239         // It's an odd integer.
 240         fchs
 241 290:    ret
 242         cfi_adjust_cfa_offset (8)
 243 291:    fstp    %st(0)          // abs(result)
 244 292:    addl    $8, %esp
 245         cfi_adjust_cfa_offset (-8)
 246         ret
 247
 248
 249         // pow(x,±0) = 1
 250         .align ALIGNARG(4)
 251 11:     fstp    %st(0)          // pop y
 252         fldl    MO(one)
 253         ret
 254
 255         // y == ±inf
 256         .align ALIGNARG(4)
 257 12:     fstp    %st(0)          // pop y
 258         fldl    MO(one)         // 1
 259         fldl    4(%esp)         // x : 1
 260         fabs                    // abs(x) : 1
 261         fucompp                 // < 1, == 1, or > 1
 262         fnstsw
 263         andb    $0x45, %ah
 264         cmpb    $0x45, %ah
 265         je      13f             // jump if x is NaN
 266
 267         cmpb    $0x40, %ah
 268         je      14f             // jump if |x| == 1
 269
 270         shlb    $1, %ah
 271         xorb    %ah, %dl
 272         andl    $2, %edx
 273         fldl    MOX(inf_zero, %edx, 4)
 274         ret
 275
 276         .align ALIGNARG(4)
 277 14:     fldl    MO(one)
 278         ret
 279
 280         .align ALIGNARG(4)
 281 13:     fldl    4(%esp)         // load x == NaN
 282         ret
 283
 284         cfi_adjust_cfa_offset (8)
 285         .align ALIGNARG(4)
 286         // x is ±inf
 287 15:     fstp    %st(0)          // y
 288         testb   $2, %dh
 289         jz      16f             // jump if x == +inf
 290
 291         // fistpll raises invalid exception for |y| >= 1L<<63, so test
 292         // that (in which case y is certainly even) before testing
 293         // whether y is odd.
 294         fld     %st             // y : y
 295         fabs                    // |y| : y
 296         fcompl  MO(p63)         // y
 297         fnstsw
 298         sahf
 299         jnc     16f
 300
 301         // We must find out whether y is an odd integer.
 302         fld     %st             // y : y
 303         fistpll (%esp)          // y
 304         fildll  (%esp)          // int(y) : y
 305         fucompp                 // <empty>
 306         fnstsw
 307         sahf
 308         jne     17f
 309
 310         // OK, the value is an integer.
 311         popl    %eax
 312         cfi_adjust_cfa_offset (-4)
 313         popl    %edx
 314         cfi_adjust_cfa_offset (-4)
 315         andb    $1, %al
 316         jz      18f             // jump if not odd
 317         // It's an odd integer.
 318         shrl    $31, %edx
 319         fldl    MOX(minf_mzero, %edx, 8)
 320         ret
 321
 322         cfi_adjust_cfa_offset (8)
 323         .align ALIGNARG(4)
 324 16:     fcompl  MO(zero)
 325         addl    $8, %esp
 326         cfi_adjust_cfa_offset (-8)
 327         fnstsw
 328         shrl    $5, %eax
 329         andl    $8, %eax
 330         fldl    MOX(inf_zero, %eax, 1)
 331         ret
 332
 333         cfi_adjust_cfa_offset (8)
 334         .align ALIGNARG(4)
 335 17:     shll    $30, %edx       // sign bit for y in right position
 336         addl    $8, %esp
 337         cfi_adjust_cfa_offset (-8)
 338 18:     shrl    $31, %edx
 339         fldl    MOX(inf_zero, %edx, 8)
 340         ret
 341
 342         cfi_adjust_cfa_offset (8)
 343         .align ALIGNARG(4)
 344         // x is ±0
 345 20:     fstp    %st(0)          // y
 346         testb   $2, %dl
 347         jz      21f             // y > 0
 348
 349         // x is ±0 and y is < 0.  We must find out whether y is an odd integer.
 350         testb   $2, %dh
 351         jz      25f
 352
 353         // fistpll raises invalid exception for |y| >= 1L<<63, so test
 354         // that (in which case y is certainly even) before testing
 355         // whether y is odd.
 356         fld     %st             // y : y
 357         fabs                    // |y| : y
 358         fcompl  MO(p63)         // y
 359         fnstsw
 360         sahf
 361         jnc     25f
 362
 363         fld     %st             // y : y
 364         fistpll (%esp)          // y
 365         fildll  (%esp)          // int(y) : y
 366         fucompp                 // <empty>
 367         fnstsw
 368         sahf
 369         jne     26f
 370
 371         // OK, the value is an integer.
 372         popl    %eax
 373         cfi_adjust_cfa_offset (-4)
 374         popl    %edx
 375         cfi_adjust_cfa_offset (-4)
 376         andb    $1, %al
 377         jz      27f             // jump if not odd
 378         // It's an odd integer.
 379         // Raise divide-by-zero exception and get minus infinity value.
 380         fldl    MO(one)
 381         fdivl   MO(zero)
 382         fchs
 383         ret
 384
 385         cfi_adjust_cfa_offset (8)
 386 25:     fstp    %st(0)
 387 26:     addl    $8, %esp
 388         cfi_adjust_cfa_offset (-8)
 389 27:     // Raise divide-by-zero exception and get infinity value.
 390         fldl    MO(one)
 391         fdivl   MO(zero)
 392         ret
 393
 394         cfi_adjust_cfa_offset (8)
 395         .align ALIGNARG(4)
 396         // x is ±0 and y is > 0.  We must find out whether y is an odd integer.
 397 21:     testb   $2, %dh
 398         jz      22f
 399
 400         // fistpll raises invalid exception for |y| >= 1L<<63, so test
 401         // that (in which case y is certainly even) before testing
 402         // whether y is odd.
 403         fcoml   MO(p63)         // y
 404         fnstsw
 405         sahf
 406         jnc     22f
 407
 408         fld     %st             // y : y
 409         fistpll (%esp)          // y
 410         fildll  (%esp)          // int(y) : y
 411         fucompp                 // <empty>
 412         fnstsw
 413         sahf
 414         jne     23f
 415
 416         // OK, the value is an integer.
 417         popl    %eax
 418         cfi_adjust_cfa_offset (-4)
 419         popl    %edx
 420         cfi_adjust_cfa_offset (-4)
 421         andb    $1, %al
 422         jz      24f             // jump if not odd
 423         // It's an odd integer.
 424         fldl    MO(mzero)
 425         ret
 426
 427         cfi_adjust_cfa_offset (8)
 428 22:     fstp    %st(0)
 429 23:     addl    $8, %esp        // Don't use 2 x pop
 430         cfi_adjust_cfa_offset (-8)
 431 24:     fldl    MO(zero)
 432         ret
 433
 434 END(__ieee754_pow)
 435 strong_alias (__ieee754_pow, __pow_finite)