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