sysdeps/i386/fpu/e_pow.S

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