sysdeps/libm-i387/e_pow.S

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