sysdeps/i386/fpu/e_powl.S

   1 /* ix87 specific implementation of pow function.
   2    Copyright (C) 1996, 1997, 1998, 1999, 2001, 2004, 2005, 2007
   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         LOAD_PIC_REG (cx)
  71 #endif
  72
  73         fnstsw
  74         movb    %ah, %dl
  75         andb    $0x45, %ah
  76         cmpb    $0x40, %ah      // is y == 0 ?
  77         je      11f
  78
  79         cmpb    $0x05, %ah      // is y == ±inf ?
  80         je      12f
  81
  82         cmpb    $0x01, %ah      // is y == NaN ?
  83         je      30f
  84
  85         fldt    4(%esp)         // x : y
  86
  87         subl    $8,%esp
  88         cfi_adjust_cfa_offset (8)
  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         /* fistpll raises invalid exception for |y| >= 1L<<63.  */
 103         fld     %st             // y : y : x
 104         fabs                    // |y| : y : x
 105         fcompl  MO(p63)         // 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         fistpll (%esp)          // y : x
 114         fildll  (%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    %eax
 122         cfi_adjust_cfa_offset (-4)
 123         popl    %edx
 124         cfi_adjust_cfa_offset (-4)
 125         orl     $0, %edx
 126         fstp    %st(0)          // x
 127         jns     4f              // y >= 0, jump
 128         fdivrl  MO(one)         // 1/x          (now referred to as x)
 129         negl    %eax
 130         adcl    $0, %edx
 131         negl    %edx
 132 4:      fldl    MO(one)         // 1 : x
 133         fxch
 134
 135 6:      shrdl   $1, %edx, %eax
 136         jnc     5f
 137         fxch
 138         fmul    %st(1)          // x : ST*x
 139         fxch
 140 5:      fmul    %st(0), %st     // x*x : ST*x
 141         shrl    $1, %edx
 142         movl    %eax, %ecx
 143         orl     %edx, %ecx
 144         jnz     6b
 145         fstp    %st(0)          // ST*x
 146         ret
 147
 148         /* y is ±NAN */
 149 30:     fldt    4(%esp)         // x : y
 150         fldl    MO(one)         // 1.0 : x : y
 151         fucomp  %st(1)          // x : y
 152         fnstsw
 153         sahf
 154         je      31f
 155         fxch                    // y : x
 156 31:     fstp    %st(1)
 157         ret
 158
 159         cfi_adjust_cfa_offset (8)
 160         .align ALIGNARG(4)
 161 2:      /* 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         fxam
 180         fnstsw
 181         andb    $0x45, %ah
 182         cmpb    $0x05, %ah      // is y*log2(x) == ±inf ?
 183         je      28f
 184         fst     %st(1)          // y*log2(x) : y*log2(x)
 185         frndint                 // int(y*log2(x)) : y*log2(x)
 186         fsubr   %st, %st(1)     // int(y*log2(x)) : fract(y*log2(x))
 187         fxch                    // fract(y*log2(x)) : int(y*log2(x))
 188         f2xm1                   // 2^fract(y*log2(x))-1 : int(y*log2(x))
 189         faddl   MO(one)         // 2^fract(y*log2(x)) : int(y*log2(x))
 190         fscale                  // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
 191         addl    $8, %esp
 192         cfi_adjust_cfa_offset (-8)
 193         fstp    %st(1)          // 2^fract(y*log2(x))*2^int(y*log2(x))
 194         ret
 195
 196         cfi_adjust_cfa_offset (8)
 197 28:     fstp    %st(1)          // y*log2(x)
 198         fldl    MO(one)         // 1 : y*log2(x)
 199         fscale                  // 2^(y*log2(x)) : y*log2(x)
 200         addl    $8, %esp
 201         cfi_adjust_cfa_offset (-8)
 202         fstp    %st(1)          // 2^(y*log2(x))
 203         ret
 204
 205         // pow(x,±0) = 1
 206         .align ALIGNARG(4)
 207 11:     fstp    %st(0)          // pop y
 208         fldl    MO(one)
 209         ret
 210
 211         // y == ±inf
 212         .align ALIGNARG(4)
 213 12:     fstp    %st(0)          // pop y
 214         fldl    MO(one)         // 1
 215         fldt    4(%esp)         // x : 1
 216         fabs                    // abs(x) : 1
 217         fucompp                 // < 1, == 1, or > 1
 218         fnstsw
 219         andb    $0x45, %ah
 220         cmpb    $0x45, %ah
 221         je      13f             // jump if x is NaN
 222
 223         cmpb    $0x40, %ah
 224         je      14f             // jump if |x| == 1
 225
 226         shlb    $1, %ah
 227         xorb    %ah, %dl
 228         andl    $2, %edx
 229         fldl    MOX(inf_zero, %edx, 4)
 230         ret
 231
 232         .align ALIGNARG(4)
 233 14:     fldl    MO(one)
 234         ret
 235
 236         .align ALIGNARG(4)
 237 13:     fldt    4(%esp)         // load x == NaN
 238         ret
 239
 240         cfi_adjust_cfa_offset (8)
 241         .align ALIGNARG(4)
 242         // x is ±inf
 243 15:     fstp    %st(0)          // y
 244         testb   $2, %dh
 245         jz      16f             // jump if x == +inf
 246
 247         // We must find out whether y is an odd integer.
 248         fld     %st             // y : y
 249         fistpll (%esp)          // y
 250         fildll  (%esp)          // int(y) : y
 251         fucompp                 // <empty>
 252         fnstsw
 253         sahf
 254         jne     17f
 255
 256         // OK, the value is an integer, but is it odd?
 257         popl    %eax
 258         cfi_adjust_cfa_offset (-4)
 259         popl    %edx
 260         cfi_adjust_cfa_offset (-4)
 261         andb    $1, %al
 262         jz      18f             // jump if not odd
 263         // It's an odd integer.
 264         shrl    $31, %edx
 265         fldl    MOX(minf_mzero, %edx, 8)
 266         ret
 267
 268         cfi_adjust_cfa_offset (8)
 269         .align ALIGNARG(4)
 270 16:     fcompl  MO(zero)
 271         addl    $8, %esp
 272         cfi_adjust_cfa_offset (-8)
 273         fnstsw
 274         shrl    $5, %eax
 275         andl    $8, %eax
 276         fldl    MOX(inf_zero, %eax, 1)
 277         ret
 278
 279         cfi_adjust_cfa_offset (8)
 280         .align ALIGNARG(4)
 281 17:     shll    $30, %edx       // sign bit for y in right position
 282         addl    $8, %esp
 283         cfi_adjust_cfa_offset (-8)
 284 18:     shrl    $31, %edx
 285         fldl    MOX(inf_zero, %edx, 8)
 286         ret
 287
 288         cfi_adjust_cfa_offset (8)
 289         .align ALIGNARG(4)
 290         // x is ±0
 291 20:     fstp    %st(0)          // y
 292         testb   $2, %dl
 293         jz      21f             // y > 0
 294
 295         // x is ±0 and y is < 0.  We must find out whether y is an odd integer.
 296         testb   $2, %dh
 297         jz      25f
 298
 299         fld     %st             // y : y
 300         fistpll (%esp)          // y
 301         fildll  (%esp)          // int(y) : y
 302         fucompp                 // <empty>
 303         fnstsw
 304         sahf
 305         jne     26f
 306
 307         // OK, the value is an integer, but is it odd?
 308         popl    %eax
 309         cfi_adjust_cfa_offset (-4)
 310         popl    %edx
 311         cfi_adjust_cfa_offset (-4)
 312         andb    $1, %al
 313         jz      27f             // jump if not odd
 314         // It's an odd integer.
 315         // Raise divide-by-zero exception and get minus infinity value.
 316         fldl    MO(one)
 317         fdivl   MO(zero)
 318         fchs
 319         ret
 320
 321         cfi_adjust_cfa_offset (8)
 322 25:     fstp    %st(0)
 323 26:     addl    $8, %esp
 324         cfi_adjust_cfa_offset (-8)
 325 27:     // Raise divide-by-zero exception and get infinity value.
 326         fldl    MO(one)
 327         fdivl   MO(zero)
 328         ret
 329
 330         cfi_adjust_cfa_offset (8)
 331         .align ALIGNARG(4)
 332         // x is ±0 and y is > 0.  We must find out whether y is an odd integer.
 333 21:     testb   $2, %dh
 334         jz      22f
 335
 336         fld     %st             // y : y
 337         fistpll (%esp)          // y
 338         fildll  (%esp)          // int(y) : y
 339         fucompp                 // <empty>
 340         fnstsw
 341         sahf
 342         jne     23f
 343
 344         // OK, the value is an integer, but is it odd?
 345         popl    %eax
 346         cfi_adjust_cfa_offset (-4)
 347         popl    %edx
 348         cfi_adjust_cfa_offset (-4)
 349         andb    $1, %al
 350         jz      24f             // jump if not odd
 351         // It's an odd integer.
 352         fldl    MO(mzero)
 353         ret
 354
 355         cfi_adjust_cfa_offset (8)
 356 22:     fstp    %st(0)
 357 23:     addl    $8, %esp        // Don't use 2 x pop
 358         cfi_adjust_cfa_offset (-8)
 359 24:     fldl    MO(zero)
 360         ret
 361
 362 END(__ieee754_powl)