sysdeps/i386/fpu/e_powf.S

   1 /* ix87 specific implementation of pow function.
   2    Copyright (C) 1996-2015 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         /* If y is even, take the absolute value of x.  Otherwise,
 133            ensure all intermediate values that might overflow have the
 134            sign of x.  */
 135         testb   $1, %dl
 136         jnz     6f
 137         fabs
 138
 139 6:      shrl    $1, %edx
 140         jnc     5f
 141         fxch
 142         fabs
 143         fmul    %st(1)          // x : ST*x
 144         fxch
 145 5:      fld     %st             // x : x : ST*x
 146         fabs                    // |x| : x : ST*x
 147         fmulp                   // |x|*x : ST*x
 148         testl   %edx, %edx
 149         jnz     6b
 150         fstp    %st(0)          // ST*x
 151         ret
 152
 153         /* y is ±NAN */
 154 30:     flds    4(%esp)         // x : y
 155         fldl    MO(one)         // 1.0 : x : y
 156         fucomp  %st(1)          // x : y
 157         fnstsw
 158         sahf
 159         je      31f
 160         fxch                    // y : x
 161 31:     fstp    %st(1)
 162         ret
 163
 164         cfi_adjust_cfa_offset (4)
 165         .align ALIGNARG(4)
 166 2:      /* y is a large integer (so even).  */
 167         fxch                    // x : y
 168         fabs                    // |x| : y
 169         fxch                    // y : x
 170         .align ALIGNARG(4)
 171 3:      /* y is a real number.  */
 172         fxch                    // x : y
 173         fldl    MO(one)         // 1.0 : x : y
 174         fldl    MO(limit)       // 0.29 : 1.0 : x : y
 175         fld     %st(2)          // x : 0.29 : 1.0 : x : y
 176         fsub    %st(2)          // x-1 : 0.29 : 1.0 : x : y
 177         fabs                    // |x-1| : 0.29 : 1.0 : x : y
 178         fucompp                 // 1.0 : x : y
 179         fnstsw
 180         fxch                    // x : 1.0 : y
 181         sahf
 182         ja      7f
 183         fsub    %st(1)          // x-1 : 1.0 : y
 184         fyl2xp1                 // log2(x) : y
 185         jmp     8f
 186
 187 7:      fyl2x                   // log2(x) : y
 188 8:      fmul    %st(1)          // y*log2(x) : y
 189         fst     %st(1)          // y*log2(x) : y*log2(x)
 190         frndint                 // int(y*log2(x)) : y*log2(x)
 191         fsubr   %st, %st(1)     // int(y*log2(x)) : fract(y*log2(x))
 192         fxch                    // fract(y*log2(x)) : int(y*log2(x))
 193         f2xm1                   // 2^fract(y*log2(x))-1 : int(y*log2(x))
 194         faddl   MO(one)         // 2^fract(y*log2(x)) : int(y*log2(x))
 195         fscale                  // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
 196 32:     addl    $4, %esp
 197         cfi_adjust_cfa_offset (-4)
 198         fstp    %st(1)          // 2^fract(y*log2(x))*2^int(y*log2(x))
 199         ret
 200
 201
 202         // pow(x,±0) = 1
 203         .align ALIGNARG(4)
 204 11:     fstp    %st(0)          // pop y
 205         fldl    MO(one)
 206         ret
 207
 208         // y == ±inf
 209         .align ALIGNARG(4)
 210 12:     fstp    %st(0)          // pop y
 211         fldl    MO(one)         // 1
 212         flds    4(%esp)         // x : 1
 213         fabs                    // abs(x) : 1
 214         fucompp                 // < 1, == 1, or > 1
 215         fnstsw
 216         andb    $0x45, %ah
 217         cmpb    $0x45, %ah
 218         je      13f             // jump if x is NaN
 219
 220         cmpb    $0x40, %ah
 221         je      14f             // jump if |x| == 1
 222
 223         shlb    $1, %ah
 224         xorb    %ah, %dl
 225         andl    $2, %edx
 226         fldl    MOX(inf_zero, %edx, 4)
 227         ret
 228
 229         .align ALIGNARG(4)
 230 14:     fldl    MO(one)
 231         ret
 232
 233         .align ALIGNARG(4)
 234 13:     flds    4(%esp)         // load x == NaN
 235         ret
 236
 237         cfi_adjust_cfa_offset (4)
 238         .align ALIGNARG(4)
 239         // x is ±inf
 240 15:     fstp    %st(0)          // y
 241         testb   $2, %dh
 242         jz      16f             // jump if x == +inf
 243
 244         // fistpl raises invalid exception for |y| >= 1L<<31, so test
 245         // that (in which case y is certainly even) before testing
 246         // whether y is odd.
 247         fld     %st             // y : y
 248         fabs                    // |y| : y
 249         fcompl  MO(p31)         // y
 250         fnstsw
 251         sahf
 252         jnc     16f
 253
 254         // We must find out whether y is an odd integer.
 255         fld     %st             // y : y
 256         fistpl  (%esp)          // y
 257         fildl   (%esp)          // int(y) : y
 258         fucompp                 // <empty>
 259         fnstsw
 260         sahf
 261         jne     17f
 262
 263         // OK, the value is an integer.
 264         popl    %edx
 265         cfi_adjust_cfa_offset (-4)
 266         testb   $1, %dl
 267         jz      18f             // jump if not odd
 268         // It's an odd integer.
 269         shrl    $31, %edx
 270         fldl    MOX(minf_mzero, %edx, 8)
 271         ret
 272
 273         cfi_adjust_cfa_offset (4)
 274         .align ALIGNARG(4)
 275 16:     fcompl  MO(zero)
 276         addl    $4, %esp
 277         cfi_adjust_cfa_offset (-4)
 278         fnstsw
 279         shrl    $5, %eax
 280         andl    $8, %eax
 281         fldl    MOX(inf_zero, %eax, 1)
 282         ret
 283
 284         cfi_adjust_cfa_offset (4)
 285         .align ALIGNARG(4)
 286 17:     shll    $30, %edx       // sign bit for y in right position
 287         addl    $4, %esp
 288         cfi_adjust_cfa_offset (-4)
 289 18:     shrl    $31, %edx
 290         fldl    MOX(inf_zero, %edx, 8)
 291         ret
 292
 293         cfi_adjust_cfa_offset (4)
 294         .align ALIGNARG(4)
 295         // x is ±0
 296 20:     fstp    %st(0)          // y
 297         testb   $2, %dl
 298         jz      21f             // y > 0
 299
 300         // x is ±0 and y is < 0.  We must find out whether y is an odd integer.
 301         testb   $2, %dh
 302         jz      25f
 303
 304         // fistpl raises invalid exception for |y| >= 1L<<31, so test
 305         // that (in which case y is certainly even) before testing
 306         // whether y is odd.
 307         fld     %st             // y : y
 308         fabs                    // |y| : y
 309         fcompl  MO(p31)         // y
 310         fnstsw
 311         sahf
 312         jnc     25f
 313
 314         fld     %st             // y : y
 315         fistpl  (%esp)          // y
 316         fildl   (%esp)          // int(y) : y
 317         fucompp                 // <empty>
 318         fnstsw
 319         sahf
 320         jne     26f
 321
 322         // OK, the value is an integer.
 323         popl    %edx
 324         cfi_adjust_cfa_offset (-4)
 325         testb   $1, %dl
 326         jz      27f             // jump if not odd
 327         // It's an odd integer.
 328         // Raise divide-by-zero exception and get minus infinity value.
 329         fldl    MO(one)
 330         fdivl   MO(zero)
 331         fchs
 332         ret
 333
 334         cfi_adjust_cfa_offset (4)
 335 25:     fstp    %st(0)
 336 26:     addl    $4, %esp
 337         cfi_adjust_cfa_offset (-4)
 338 27:     // Raise divide-by-zero exception and get infinity value.
 339         fldl    MO(one)
 340         fdivl   MO(zero)
 341         ret
 342
 343         cfi_adjust_cfa_offset (4)
 344         .align ALIGNARG(4)
 345         // x is ±0 and y is > 0.  We must find out whether y is an odd integer.
 346 21:     testb   $2, %dh
 347         jz      22f
 348
 349         // fistpl raises invalid exception for |y| >= 1L<<31, so test
 350         // that (in which case y is certainly even) before testing
 351         // whether y is odd.
 352         fcoml   MO(p31)         // y
 353         fnstsw
 354         sahf
 355         jnc     22f
 356
 357         fld     %st             // y : y
 358         fistpl  (%esp)          // y
 359         fildl   (%esp)          // int(y) : y
 360         fucompp                 // <empty>
 361         fnstsw
 362         sahf
 363         jne     23f
 364
 365         // OK, the value is an integer.
 366         popl    %edx
 367         cfi_adjust_cfa_offset (-4)
 368         testb   $1, %dl
 369         jz      24f             // jump if not odd
 370         // It's an odd integer.
 371         fldl    MO(mzero)
 372         ret
 373
 374         cfi_adjust_cfa_offset (4)
 375 22:     fstp    %st(0)
 376 23:     addl    $4, %esp        // Don't use pop.
 377         cfi_adjust_cfa_offset (-4)
 378 24:     fldl    MO(zero)
 379         ret
 380
 381 END(__ieee754_powf)
 382 strong_alias (__ieee754_powf, __powf_finite)