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