sysdeps/i386/fpu/e_pow.S

   1 /* ix87 specific implementation of pow function.
   2    Copyright (C) 1996-1999, 2001, 2004-2005, 2007, 2011-2012
   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, see
  19    <http://www.gnu.org/licenses/>.  */
  20
  21 #include <machine/asm.h>
  22
  23         .section .rodata.cst8,"aM",@progbits,8
  24
  25         .p2align 3
  26         ASM_TYPE_DIRECTIVE(one,@object)
  27 one:    .double 1.0
  28         ASM_SIZE_DIRECTIVE(one)
  29         ASM_TYPE_DIRECTIVE(limit,@object)
  30 limit:  .double 0.29
  31         ASM_SIZE_DIRECTIVE(limit)
  32         ASM_TYPE_DIRECTIVE(p63,@object)
  33 p63:    .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43
  34         ASM_SIZE_DIRECTIVE(p63)
  35
  36         .section .rodata.cst16,"aM",@progbits,16
  37
  38         .p2align 3
  39         ASM_TYPE_DIRECTIVE(infinity,@object)
  40 inf_zero:
  41 infinity:
  42         .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
  43         ASM_SIZE_DIRECTIVE(infinity)
  44         ASM_TYPE_DIRECTIVE(zero,@object)
  45 zero:   .double 0.0
  46         ASM_SIZE_DIRECTIVE(zero)
  47         ASM_TYPE_DIRECTIVE(minf_mzero,@object)
  48 minf_mzero:
  49 minfinity:
  50         .byte 0, 0, 0, 0, 0, 0, 0xf0, 0xff
  51 mzero:
  52         .byte 0, 0, 0, 0, 0, 0, 0, 0x80
  53         ASM_SIZE_DIRECTIVE(minf_mzero)
  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
  68 #ifdef  PIC
  69         LOAD_PIC_REG (cx)
  70 #endif
  71
  72         fnstsw
  73         movb    %ah, %dl
  74         andb    $0x45, %ah
  75         cmpb    $0x40, %ah      // is y == 0 ?
  76         je      11f
  77
  78         cmpb    $0x05, %ah      // is y == ±inf ?
  79         je      12f
  80
  81         cmpb    $0x01, %ah      // is y == NaN ?
  82         je      30f
  83
  84         fldl    4(%esp)         // x : y
  85
  86         subl    $8,%esp
  87         cfi_adjust_cfa_offset (8)
  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         /* fistpll raises invalid exception for |y| >= 1L<<63.  */
 102         fld     %st             // y : y : x
 103         fabs                    // |y| : y : x
 104         fcompl  MO(p63)         // y : x
 105         fnstsw
 106         sahf
 107         jnc     2f
 108
 109         /* First see whether `y' is a natural number.  In this case we
 110            can use a more precise algorithm.  */
 111         fld     %st             // y : y : x
 112         fistpll (%esp)          // y : x
 113         fildll  (%esp)          // int(y) : y : x
 114         fucomp  %st(1)          // y : x
 115         fnstsw
 116         sahf
 117         jne     2f
 118
 119         /* OK, we have an integer value for y.  */
 120         popl    %eax
 121         cfi_adjust_cfa_offset (-4)
 122         popl    %edx
 123         cfi_adjust_cfa_offset (-4)
 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:     fldl    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         cfi_adjust_cfa_offset (8)
 159         .align ALIGNARG(4)
 160 2:      /* y is a real number.  */
 161         fxch                    // x : y
 162         fldl    MO(one)         // 1.0 : x : y
 163         fldl    MO(limit)       // 0.29 : 1.0 : x : y
 164         fld     %st(2)          // x : 0.29 : 1.0 : x : y
 165         fsub    %st(2)          // x-1 : 0.29 : 1.0 : x : y
 166         fabs                    // |x-1| : 0.29 : 1.0 : x : y
 167         fucompp                 // 1.0 : x : y
 168         fnstsw
 169         fxch                    // x : 1.0 : y
 170         sahf
 171         ja      7f
 172         fsub    %st(1)          // x-1 : 1.0 : y
 173         fyl2xp1                 // log2(x) : y
 174         jmp     8f
 175
 176 7:      fyl2x                   // log2(x) : y
 177 8:      fmul    %st(1)          // y*log2(x) : y
 178         fst     %st(1)          // y*log2(x) : y*log2(x)
 179         frndint                 // int(y*log2(x)) : y*log2(x)
 180         fsubr   %st, %st(1)     // int(y*log2(x)) : fract(y*log2(x))
 181         fxch                    // fract(y*log2(x)) : int(y*log2(x))
 182         f2xm1                   // 2^fract(y*log2(x))-1 : int(y*log2(x))
 183         faddl   MO(one)         // 2^fract(y*log2(x)) : int(y*log2(x))
 184         fscale                  // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
 185         addl    $8, %esp
 186         cfi_adjust_cfa_offset (-8)
 187         fstp    %st(1)          // 2^fract(y*log2(x))*2^int(y*log2(x))
 188         ret
 189
 190
 191         // pow(x,±0) = 1
 192         .align ALIGNARG(4)
 193 11:     fstp    %st(0)          // pop y
 194         fldl    MO(one)
 195         ret
 196
 197         // y == ±inf
 198         .align ALIGNARG(4)
 199 12:     fstp    %st(0)          // pop y
 200         fldl    MO(one)         // 1
 201         fldl    4(%esp)         // x : 1
 202         fabs                    // abs(x) : 1
 203         fucompp                 // < 1, == 1, or > 1
 204         fnstsw
 205         andb    $0x45, %ah
 206         cmpb    $0x45, %ah
 207         je      13f             // jump if x is NaN
 208
 209         cmpb    $0x40, %ah
 210         je      14f             // jump if |x| == 1
 211
 212         shlb    $1, %ah
 213         xorb    %ah, %dl
 214         andl    $2, %edx
 215         fldl    MOX(inf_zero, %edx, 4)
 216         ret
 217
 218         .align ALIGNARG(4)
 219 14:     fldl    MO(one)
 220         ret
 221
 222         .align ALIGNARG(4)
 223 13:     fldl    4(%esp)         // load x == NaN
 224         ret
 225
 226         cfi_adjust_cfa_offset (8)
 227         .align ALIGNARG(4)
 228         // x is ±inf
 229 15:     fstp    %st(0)          // y
 230         testb   $2, %dh
 231         jz      16f             // jump if x == +inf
 232
 233         // fistpll raises invalid exception for |y| >= 1L<<63, so test
 234         // that (in which case y is certainly even) before testing
 235         // whether y is odd.
 236         fld     %st             // y : y
 237         fabs                    // |y| : y
 238         fcompl  MO(p63)         // y
 239         fnstsw
 240         sahf
 241         jnc     16f
 242
 243         // We must find out whether y is an odd integer.
 244         fld     %st             // y : y
 245         fistpll (%esp)          // y
 246         fildll  (%esp)          // int(y) : y
 247         fucompp                 // <empty>
 248         fnstsw
 249         sahf
 250         jne     17f
 251
 252         // OK, the value is an integer.
 253         popl    %eax
 254         cfi_adjust_cfa_offset (-4)
 255         popl    %edx
 256         cfi_adjust_cfa_offset (-4)
 257         andb    $1, %al
 258         jz      18f             // jump if not odd
 259         // It's an odd integer.
 260         shrl    $31, %edx
 261         fldl    MOX(minf_mzero, %edx, 8)
 262         ret
 263
 264         cfi_adjust_cfa_offset (8)
 265         .align ALIGNARG(4)
 266 16:     fcompl  MO(zero)
 267         addl    $8, %esp
 268         cfi_adjust_cfa_offset (-8)
 269         fnstsw
 270         shrl    $5, %eax
 271         andl    $8, %eax
 272         fldl    MOX(inf_zero, %eax, 1)
 273         ret
 274
 275         cfi_adjust_cfa_offset (8)
 276         .align ALIGNARG(4)
 277 17:     shll    $30, %edx       // sign bit for y in right position
 278         addl    $8, %esp
 279         cfi_adjust_cfa_offset (-8)
 280 18:     shrl    $31, %edx
 281         fldl    MOX(inf_zero, %edx, 8)
 282         ret
 283
 284         cfi_adjust_cfa_offset (8)
 285         .align ALIGNARG(4)
 286         // x is ±0
 287 20:     fstp    %st(0)          // y
 288         testb   $2, %dl
 289         jz      21f             // y > 0
 290
 291         // x is ±0 and y is < 0.  We must find out whether y is an odd integer.
 292         testb   $2, %dh
 293         jz      25f
 294
 295         // fistpll raises invalid exception for |y| >= 1L<<63, so test
 296         // that (in which case y is certainly even) before testing
 297         // whether y is odd.
 298         fld     %st             // y : y
 299         fabs                    // |y| : y
 300         fcompl  MO(p63)         // y
 301         fnstsw
 302         sahf
 303         jnc     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.
 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         // fistpll raises invalid exception for |y| >= 1L<<63, so test
 343         // that (in which case y is certainly even) before testing
 344         // whether y is odd.
 345         fcoml   MO(p63)         // y
 346         fnstsw
 347         sahf
 348         jnc     22f
 349
 350         fld     %st             // y : y
 351         fistpll (%esp)          // y
 352         fildll  (%esp)          // int(y) : y
 353         fucompp                 // <empty>
 354         fnstsw
 355         sahf
 356         jne     23f
 357
 358         // OK, the value is an integer.
 359         popl    %eax
 360         cfi_adjust_cfa_offset (-4)
 361         popl    %edx
 362         cfi_adjust_cfa_offset (-4)
 363         andb    $1, %al
 364         jz      24f             // jump if not odd
 365         // It's an odd integer.
 366         fldl    MO(mzero)
 367         ret
 368
 369         cfi_adjust_cfa_offset (8)
 370 22:     fstp    %st(0)
 371 23:     addl    $8, %esp        // Don't use 2 x pop
 372         cfi_adjust_cfa_offset (-8)
 373 24:     fldl    MO(zero)
 374         ret
 375
 376 END(__ieee754_pow)
 377 strong_alias (__ieee754_pow, __pow_finite)