sysdeps/i386/fpu/e_powf.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 p31,@object
  32 p31:    .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x41
  33         ASM_SIZE_DIRECTIVE(p31)
  34
  35         .section .rodata.cst16,"aM",@progbits,16
  36
  37         .p2align 3
  38         .type infinity,@object
  39 inf_zero:
  40 infinity:
  41         .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
  42         ASM_SIZE_DIRECTIVE(infinity)
  43         .type zero,@object
  44 zero:   .double 0.0
  45         ASM_SIZE_DIRECTIVE(zero)
  46         .type minf_mzero,@object
  47 minf_mzero:
  48 minfinity:
  49         .byte 0, 0, 0, 0, 0, 0, 0xf0, 0xff
  50 mzero:
  51         .byte 0, 0, 0, 0, 0, 0, 0, 0x80
  52         ASM_SIZE_DIRECTIVE(minf_mzero)
  53
  54 #ifdef PIC
  55 # define MO(op) op##@GOTOFF(%ecx)
  56 # define MOX(op,x,f) op##@GOTOFF(%ecx,x,f)
  57 #else
  58 # define MO(op) op
  59 # define MOX(op,x,f) op(,x,f)
  60 #endif
  61
  62         .text
  63 ENTRY(__ieee754_powf)
  64         flds    8(%esp) // y
  65         fxam
  66
  67 #ifdef  PIC
  68         LOAD_PIC_REG (cx)
  69 #endif
  70
  71         fnstsw
  72         movb    %ah, %dl
  73         andb    $0x45, %ah
  74         cmpb    $0x40, %ah      // is y == 0 ?
  75         je      11f
  76
  77         cmpb    $0x05, %ah      // is y == ±inf ?
  78         je      12f
  79
  80         cmpb    $0x01, %ah      // is y == NaN ?
  81         je      30f
  82
  83         flds    4(%esp)         // x : y
  84
  85         subl    $4, %esp
  86         cfi_adjust_cfa_offset (4)
  87
  88         fxam
  89         fnstsw
  90         movb    %ah, %dh
  91         andb    $0x45, %ah
  92         cmpb    $0x40, %ah
  93         je      20f             // x is ±0
  94
  95         cmpb    $0x05, %ah
  96         je      15f             // x is ±inf
  97
  98         cmpb    $0x01, %ah
  99         je      32f             // x is NaN
 100
 101         fxch                    // y : x
 102
 103         /* fistpl raises invalid exception for |y| >= 1L<<31.  */
 104         fld     %st             // y : y : x
 105         fabs                    // |y| : y : x
 106         fcompl  MO(p31)         // 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         fistpl  (%esp)          // y : x
 115         fildl   (%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.  */
 122         popl    %edx
 123         cfi_adjust_cfa_offset (-4)
 124         orl     $0, %edx
 125         fstp    %st(0)          // x
 126         jns     4f              // y >= 0, jump
 127         fdivrl  MO(one)         // 1/x          (now referred to as x)
 128         negl    %edx
 129 4:      fldl    MO(one)         // 1 : x
 130         fxch
 131
 132 6:      shrl    $1, %edx
 133         jnc     5f
 134         fxch
 135         fmul    %st(1)          // x : ST*x
 136         fxch
 137 5:      fmul    %st(0), %st     // x*x : ST*x
 138         testl   %edx, %edx
 139         jnz     6b
 140         fstp    %st(0)          // ST*x
 141         ret
 142
 143         /* y is ±NAN */
 144 30:     flds    4(%esp)         // x : y
 145         fldl    MO(one)         // 1.0 : x : y
 146         fucomp  %st(1)          // x : y
 147         fnstsw
 148         sahf
 149         je      31f
 150         fxch                    // y : x
 151 31:     fstp    %st(1)
 152         ret
 153
 154         cfi_adjust_cfa_offset (4)
 155         .align ALIGNARG(4)
 156 2:      /* y is a large integer (so even).  */
 157         fxch                    // x : y
 158         fabs                    // |x| : y
 159         fxch                    // y : x
 160         .align ALIGNARG(4)
 161 3:      /* y is a real number.  */
 162         fxch                    // x : y
 163         fldl    MO(one)         // 1.0 : x : y
 164         fldl    MO(limit)       // 0.29 : 1.0 : x : y
 165         fld     %st(2)          // x : 0.29 : 1.0 : x : y
 166         fsub    %st(2)          // x-1 : 0.29 : 1.0 : x : y
 167         fabs                    // |x-1| : 0.29 : 1.0 : x : y
 168         fucompp                 // 1.0 : x : y
 169         fnstsw
 170         fxch                    // x : 1.0 : y
 171         sahf
 172         ja      7f
 173         fsub    %st(1)          // x-1 : 1.0 : y
 174         fyl2xp1                 // log2(x) : y
 175         jmp     8f
 176
 177 7:      fyl2x                   // log2(x) : y
 178 8:      fmul    %st(1)          // y*log2(x) : y
 179         fst     %st(1)          // y*log2(x) : y*log2(x)
 180         frndint                 // int(y*log2(x)) : y*log2(x)
 181         fsubr   %st, %st(1)     // int(y*log2(x)) : fract(y*log2(x))
 182         fxch                    // fract(y*log2(x)) : int(y*log2(x))
 183         f2xm1                   // 2^fract(y*log2(x))-1 : int(y*log2(x))
 184         faddl   MO(one)         // 2^fract(y*log2(x)) : int(y*log2(x))
 185         fscale                  // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
 186 32:     addl    $4, %esp
 187         cfi_adjust_cfa_offset (-4)
 188         fstp    %st(1)          // 2^fract(y*log2(x))*2^int(y*log2(x))
 189         ret
 190
 191
 192         // pow(x,±0) = 1
 193         .align ALIGNARG(4)
 194 11:     fstp    %st(0)          // pop y
 195         fldl    MO(one)
 196         ret
 197
 198         // y == ±inf
 199         .align ALIGNARG(4)
 200 12:     fstp    %st(0)          // pop y
 201         fldl    MO(one)         // 1
 202         flds    4(%esp)         // x : 1
 203         fabs                    // abs(x) : 1
 204         fucompp                 // < 1, == 1, or > 1
 205         fnstsw
 206         andb    $0x45, %ah
 207         cmpb    $0x45, %ah
 208         je      13f             // jump if x is NaN
 209
 210         cmpb    $0x40, %ah
 211         je      14f             // jump if |x| == 1
 212
 213         shlb    $1, %ah
 214         xorb    %ah, %dl
 215         andl    $2, %edx
 216         fldl    MOX(inf_zero, %edx, 4)
 217         ret
 218
 219         .align ALIGNARG(4)
 220 14:     fldl    MO(one)
 221         ret
 222
 223         .align ALIGNARG(4)
 224 13:     flds    4(%esp)         // load x == NaN
 225         ret
 226
 227         cfi_adjust_cfa_offset (4)
 228         .align ALIGNARG(4)
 229         // x is ±inf
 230 15:     fstp    %st(0)          // y
 231         testb   $2, %dh
 232         jz      16f             // jump if x == +inf
 233
 234         // fistpl raises invalid exception for |y| >= 1L<<31, so test
 235         // that (in which case y is certainly even) before testing
 236         // whether y is odd.
 237         fld     %st             // y : y
 238         fabs                    // |y| : y
 239         fcompl  MO(p31)         // y
 240         fnstsw
 241         sahf
 242         jnc     16f
 243
 244         // We must find out whether y is an odd integer.
 245         fld     %st             // y : y
 246         fistpl  (%esp)          // y
 247         fildl   (%esp)          // int(y) : y
 248         fucompp                 // <empty>
 249         fnstsw
 250         sahf
 251         jne     17f
 252
 253         // OK, the value is an integer.
 254         popl    %edx
 255         cfi_adjust_cfa_offset (-4)
 256         testb   $1, %dl
 257         jz      18f             // jump if not odd
 258         // It's an odd integer.
 259         shrl    $31, %edx
 260         fldl    MOX(minf_mzero, %edx, 8)
 261         ret
 262
 263         cfi_adjust_cfa_offset (4)
 264         .align ALIGNARG(4)
 265 16:     fcompl  MO(zero)
 266         addl    $4, %esp
 267         cfi_adjust_cfa_offset (-4)
 268         fnstsw
 269         shrl    $5, %eax
 270         andl    $8, %eax
 271         fldl    MOX(inf_zero, %eax, 1)
 272         ret
 273
 274         cfi_adjust_cfa_offset (4)
 275         .align ALIGNARG(4)
 276 17:     shll    $30, %edx       // sign bit for y in right position
 277         addl    $4, %esp
 278         cfi_adjust_cfa_offset (-4)
 279 18:     shrl    $31, %edx
 280         fldl    MOX(inf_zero, %edx, 8)
 281         ret
 282
 283         cfi_adjust_cfa_offset (4)
 284         .align ALIGNARG(4)
 285         // x is ±0
 286 20:     fstp    %st(0)          // y
 287         testb   $2, %dl
 288         jz      21f             // y > 0
 289
 290         // x is ±0 and y is < 0.  We must find out whether y is an odd integer.
 291         testb   $2, %dh
 292         jz      25f
 293
 294         // fistpl raises invalid exception for |y| >= 1L<<31, so test
 295         // that (in which case y is certainly even) before testing
 296         // whether y is odd.
 297         fld     %st             // y : y
 298         fabs                    // |y| : y
 299         fcompl  MO(p31)         // y
 300         fnstsw
 301         sahf
 302         jnc     25f
 303
 304         fld     %st             // y : y
 305         fistpl  (%esp)          // y
 306         fildl   (%esp)          // int(y) : y
 307         fucompp                 // <empty>
 308         fnstsw
 309         sahf
 310         jne     26f
 311
 312         // OK, the value is an integer.
 313         popl    %edx
 314         cfi_adjust_cfa_offset (-4)
 315         testb   $1, %dl
 316         jz      27f             // jump if not odd
 317         // It's an odd integer.
 318         // Raise divide-by-zero exception and get minus infinity value.
 319         fldl    MO(one)
 320         fdivl   MO(zero)
 321         fchs
 322         ret
 323
 324         cfi_adjust_cfa_offset (4)
 325 25:     fstp    %st(0)
 326 26:     addl    $4, %esp
 327         cfi_adjust_cfa_offset (-4)
 328 27:     // Raise divide-by-zero exception and get infinity value.
 329         fldl    MO(one)
 330         fdivl   MO(zero)
 331         ret
 332
 333         cfi_adjust_cfa_offset (4)
 334         .align ALIGNARG(4)
 335         // x is ±0 and y is > 0.  We must find out whether y is an odd integer.
 336 21:     testb   $2, %dh
 337         jz      22f
 338
 339         // fistpl raises invalid exception for |y| >= 1L<<31, so test
 340         // that (in which case y is certainly even) before testing
 341         // whether y is odd.
 342         fcoml   MO(p31)         // y
 343         fnstsw
 344         sahf
 345         jnc     22f
 346
 347         fld     %st             // y : y
 348         fistpl  (%esp)          // y
 349         fildl   (%esp)          // int(y) : y
 350         fucompp                 // <empty>
 351         fnstsw
 352         sahf
 353         jne     23f
 354
 355         // OK, the value is an integer.
 356         popl    %edx
 357         cfi_adjust_cfa_offset (-4)
 358         testb   $1, %dl
 359         jz      24f             // jump if not odd
 360         // It's an odd integer.
 361         fldl    MO(mzero)
 362         ret
 363
 364         cfi_adjust_cfa_offset (4)
 365 22:     fstp    %st(0)
 366 23:     addl    $4, %esp        // Don't use pop.
 367         cfi_adjust_cfa_offset (-4)
 368 24:     fldl    MO(zero)
 369         ret
 370
 371 END(__ieee754_powf)
 372 strong_alias (__ieee754_powf, __powf_finite)