sysdeps/i386/fpu/e_powf.S

   1 /* ix87 specific implementation of pow function.
   2    Copyright (C) 1996, 1997, 1999, 2001, 2004 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, write to the Free
  18    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
  19    02111-1307 USA.  */
  20
  21 #include <machine/asm.h>
  22
  23 #ifdef __ELF__
  24         .section .rodata
  25 #else
  26         .text
  27 #endif
  28
  29         .align ALIGNARG(4)
  30         ASM_TYPE_DIRECTIVE(infinity,@object)
  31 inf_zero:
  32 infinity:
  33         .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
  34         ASM_SIZE_DIRECTIVE(infinity)
  35         ASM_TYPE_DIRECTIVE(zero,@object)
  36 zero:   .double 0.0
  37         ASM_SIZE_DIRECTIVE(zero)
  38         ASM_TYPE_DIRECTIVE(minf_mzero,@object)
  39 minf_mzero:
  40 minfinity:
  41         .byte 0, 0, 0, 0, 0, 0, 0xf0, 0xff
  42 mzero:
  43         .byte 0, 0, 0, 0, 0, 0, 0, 0x80
  44         ASM_SIZE_DIRECTIVE(minf_mzero)
  45         ASM_TYPE_DIRECTIVE(one,@object)
  46 one:    .double 1.0
  47         ASM_SIZE_DIRECTIVE(one)
  48         ASM_TYPE_DIRECTIVE(limit,@object)
  49 limit:  .double 0.29
  50         ASM_SIZE_DIRECTIVE(limit)
  51         ASM_TYPE_DIRECTIVE(p31,@object)
  52 p31:    .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x41
  53         ASM_SIZE_DIRECTIVE(p31)
  54
  55 #ifdef PIC
  56 #define MO(op) op##@GOTOFF(%ecx)
  57 #define MOX(op,x,f) op##@GOTOFF(%ecx,x,f)
  58 #else
  59 #define MO(op) op
  60 #define MOX(op,x,f) op(,x,f)
  61 #endif
  62
  63         .text
  64 ENTRY(__ieee754_powf)
  65         flds    8(%esp) // y
  66         fxam
  67
  68 #ifdef  PIC
  69         call    1f
  70 1:      popl    %ecx
  71         addl    $_GLOBAL_OFFSET_TABLE_+[.-1b], %ecx
  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         flds    4(%esp)         // x : y
  87
  88         subl    $4, %esp
  89
  90         fxam
  91         fnstsw
  92         movb    %ah, %dh
  93         andb    $0x45, %ah
  94         cmpb    $0x40, %ah
  95         je      20f             // x is ±0
  96
  97         cmpb    $0x05, %ah
  98         je      15f             // x is ±inf
  99
 100         fxch                    // y : x
 101
 102         /* fistpl raises invalid exception for |y| >= 1L<<31.  */
 103         fld     %st             // y : y : x
 104         fabs                    // |y| : y : x
 105         fcompl  MO(p31)         // y : x
 106         fnstsw
 107         sahf
 108         jnc     2f
 109
 110         /* First see whether `y' is a natural number.  In this case we
 111            can use a more precise algorithm.  */
 112         fld     %st             // y : y : x
 113         fistpl  (%esp)          // y : x
 114         fildl   (%esp)          // int(y) : y : x
 115         fucomp  %st(1)          // y : x
 116         fnstsw
 117         sahf
 118         jne     2f
 119
 120         /* OK, we have an integer value for y.  */
 121         popl    %edx
 122         orl     $0, %edx
 123         fstp    %st(0)          // x
 124         jns     4f              // y >= 0, jump
 125         fdivrl  MO(one)         // 1/x          (now referred to as x)
 126         negl    %edx
 127 4:      fldl    MO(one)         // 1 : x
 128         fxch
 129
 130 6:      shrl    $1, %edx
 131         jnc     5f
 132         fxch
 133         fmul    %st(1)          // x : ST*x
 134         fxch
 135 5:      fmul    %st(0), %st     // x*x : ST*x
 136         testl   %edx, %edx
 137         jnz     6b
 138         fstp    %st(0)          // ST*x
 139         ret
 140
 141         /* y is ±NAN */
 142 30:     flds    4(%esp)         // x : y
 143         fldl    MO(one)         // 1.0 : x : y
 144         fucomp  %st(1)          // x : y
 145         fnstsw
 146         sahf
 147         je      31f
 148         fxch                    // y : x
 149 31:     fstp    %st(1)
 150         ret
 151
 152         .align ALIGNARG(4)
 153 2:      /* y is a real number.  */
 154         fxch                    // x : y
 155         fldl    MO(one)         // 1.0 : x : y
 156         fld     %st(1)          // x : 1.0 : x : y
 157         fsub    %st(1)          // x-1 : 1.0 : x : y
 158         fabs                    // |x-1| : 1.0 : x : y
 159         fcompl  MO(limit)       // 1.0 : x : y
 160         fnstsw
 161         fxch                    // x : 1.0 : y
 162         sahf
 163         ja      7f
 164         fsub    %st(1)          // x-1 : 1.0 : y
 165         fyl2xp1                 // log2(x) : y
 166         jmp     8f
 167
 168 7:      fyl2x                   // log2(x) : y
 169 8:      fmul    %st(1)          // y*log2(x) : y
 170         fst     %st(1)          // y*log2(x) : y*log2(x)
 171         frndint                 // int(y*log2(x)) : y*log2(x)
 172         fsubr   %st, %st(1)     // int(y*log2(x)) : fract(y*log2(x))
 173         fxch                    // fract(y*log2(x)) : int(y*log2(x))
 174         f2xm1                   // 2^fract(y*log2(x))-1 : int(y*log2(x))
 175         faddl   MO(one)         // 2^fract(y*log2(x)) : int(y*log2(x))
 176         fscale                  // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
 177         addl    $4, %esp
 178         fstp    %st(1)          // 2^fract(y*log2(x))*2^int(y*log2(x))
 179         ret
 180
 181
 182         // pow(x,±0) = 1
 183         .align ALIGNARG(4)
 184 11:     fstp    %st(0)          // pop y
 185         fldl    MO(one)
 186         ret
 187
 188         // y == ±inf
 189         .align ALIGNARG(4)
 190 12:     fstp    %st(0)          // pop y
 191         flds    4(%esp)         // x
 192         fabs
 193         fcompl  MO(one)         // < 1, == 1, or > 1
 194         fnstsw
 195         andb    $0x45, %ah
 196         cmpb    $0x45, %ah
 197         je      13f             // jump if x is NaN
 198
 199         cmpb    $0x40, %ah
 200         je      14f             // jump if |x| == 1
 201
 202         shlb    $1, %ah
 203         xorb    %ah, %dl
 204         andl    $2, %edx
 205         fldl    MOX(inf_zero, %edx, 4)
 206         ret
 207
 208         .align ALIGNARG(4)
 209 14:     fldl    MO(one)
 210         ret
 211
 212         .align ALIGNARG(4)
 213 13:     flds    4(%esp)         // load x == NaN
 214         ret
 215
 216         .align ALIGNARG(4)
 217         // x is ±inf
 218 15:     fstp    %st(0)          // y
 219         testb   $2, %dh
 220         jz      16f             // jump if x == +inf
 221
 222         // We must find out whether y is an odd integer.
 223         fld     %st             // y : y
 224         fistpl  (%esp)          // y
 225         fildl   (%esp)          // int(y) : y
 226         fucompp                 // <empty>
 227         fnstsw
 228         sahf
 229         jne     17f
 230
 231         // OK, the value is an integer, but is the number of bits small
 232         // enough so that all are coming from the mantissa?
 233         popl    %edx
 234         testb   $1, %dl
 235         jz      18f             // jump if not odd
 236         movl    %edx, %eax
 237         orl     %edx, %edx
 238         jns     155f
 239         negl    %eax
 240 155:    cmpl    $0x01000000, %eax
 241         ja      18f             // does not fit in mantissa bits
 242         // It's an odd integer.
 243         shrl    $31, %edx
 244         fldl    MOX(minf_mzero, %edx, 8)
 245         ret
 246
 247         .align ALIGNARG(4)
 248 16:     fcompl  MO(zero)
 249         addl    $4, %esp
 250         fnstsw
 251         shrl    $5, %eax
 252         andl    $8, %eax
 253         fldl    MOX(inf_zero, %eax, 1)
 254         ret
 255
 256         .align ALIGNARG(4)
 257 17:     shll    $30, %edx       // sign bit for y in right position
 258         addl    $4, %esp
 259 18:     shrl    $31, %edx
 260         fldl    MOX(inf_zero, %edx, 8)
 261         ret
 262
 263         .align ALIGNARG(4)
 264         // x is ±0
 265 20:     fstp    %st(0)          // y
 266         testb   $2, %dl
 267         jz      21f             // y > 0
 268
 269         // x is ±0 and y is < 0.  We must find out whether y is an odd integer.
 270         testb   $2, %dh
 271         jz      25f
 272
 273         fld     %st             // y : y
 274         fistpl  (%esp)          // y
 275         fildl   (%esp)          // int(y) : y
 276         fucompp                 // <empty>
 277         fnstsw
 278         sahf
 279         jne     26f
 280
 281         // OK, the value is an integer, but is the number of bits small
 282         // enough so that all are coming from the mantissa?
 283         popl    %edx
 284         testb   $1, %dl
 285         jz      27f             // jump if not odd
 286         cmpl    $0xff000000, %edx
 287         jbe     27f             // does not fit in mantissa bits
 288         // It's an odd integer.
 289         // Raise divide-by-zero exception and get minus infinity value.
 290         fldl    MO(one)
 291         fdivl   MO(zero)
 292         fchs
 293         ret
 294
 295 25:     fstp    %st(0)
 296 26:     addl    $4, %esp
 297 27:     // Raise divide-by-zero exception and get infinity value.
 298         fldl    MO(one)
 299         fdivl   MO(zero)
 300         ret
 301
 302         .align ALIGNARG(4)
 303         // x is ±0 and y is > 0.  We must find out whether y is an odd integer.
 304 21:     testb   $2, %dh
 305         jz      22f
 306
 307         fld     %st             // y : y
 308         fistpl  (%esp)          // y
 309         fildl   (%esp)          // int(y) : y
 310         fucompp                 // <empty>
 311         fnstsw
 312         sahf
 313         jne     23f
 314
 315         // OK, the value is an integer, but is the number of bits small
 316         // enough so that all are coming from the mantissa?
 317         popl    %edx
 318         testb   $1, %dl
 319         jz      24f             // jump if not odd
 320         cmpl    $0xff000000, %edx
 321         jae     24f             // does not fit in mantissa bits
 322         // It's an odd integer.
 323         fldl    MO(mzero)
 324         ret
 325
 326 22:     fstp    %st(0)
 327 23:     addl    $4, %esp        // Don't use pop.
 328 24:     fldl    MO(zero)
 329         ret
 330
 331 END(__ieee754_powf)