sysdeps/i386/fpu/e_powl.S

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