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