sysdeps/i386/fpu/e_powl.S

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