src/math/x32/expl.s

   1 # exp(x) = 2^hi + 2^hi (2^lo - 1)
   2 # where hi+lo = log2e*x with 128bit precision
   3 # exact log2e*x calculation depends on nearest rounding mode
   4 # using the exact multiplication method of Dekker and Veltkamp
   5
   6 .global expl
   7 .type expl,@function
   8 expl:
   9         fldt 8(%esp)
  10
  11                 # interesting case: 0x1p-32 <= |x| < 16384
  12                 # check if (exponent|0x8000) is in [0xbfff-32, 0xbfff+13]
  13         mov 16(%esp), %ax
  14         or $0x8000, %ax
  15         sub $0xbfdf, %ax
  16         cmp $45, %ax
  17         jbe 2f
  18         test %ax, %ax
  19         fld1
  20         js 1f
  21                 # if |x|>=0x1p14 or nan return 2^trunc(x)
  22         fscale
  23         fstp %st(1)
  24         ret
  25                 # if |x|<0x1p-32 return 1+x
  26 1:      faddp
  27         ret
  28
  29                 # should be 0x1.71547652b82fe178p0L == 0x3fff b8aa3b29 5c17f0bc
  30                 # it will be wrong on non-nearest rounding mode
  31 2:      fldl2e
  32         sub $48, %esp
  33                 # hi = log2e_hi*x
  34                 # 2^hi = exp2l(hi)
  35         fmul %st(1),%st
  36         fld %st(0)
  37         fstpt (%esp)
  38         fstpt 16(%esp)
  39         fstpt 32(%esp)
  40         call exp2l@PLT
  41                 # if 2^hi == inf return 2^hi
  42         fld %st(0)
  43         fstpt (%esp)
  44         cmpw $0x7fff, 8(%esp)
  45         je 1f
  46         fldt 32(%esp)
  47         fldt 16(%esp)
  48                 # fpu stack: 2^hi x hi
  49                 # exact mult: x*log2e
  50         fld %st(1)
  51                 # c = 0x1p32+1
  52         movq $0x41f0000000100000,%rax
  53         pushq %rax
  54         fldl (%esp)
  55                 # xh = x - c*x + c*x
  56                 # xl = x - xh
  57         fmulp
  58         fld %st(2)
  59         fsub %st(1), %st
  60         faddp
  61         fld %st(2)
  62         fsub %st(1), %st
  63                 # yh = log2e_hi - c*log2e_hi + c*log2e_hi
  64         movq $0x3ff7154765200000,%rax
  65         pushq %rax
  66         fldl (%esp)
  67                 # fpu stack: 2^hi x hi xh xl yh
  68                 # lo = hi - xh*yh + xl*yh
  69         fld %st(2)
  70         fmul %st(1), %st
  71         fsubp %st, %st(4)
  72         fmul %st(1), %st
  73         faddp %st, %st(3)
  74                 # yl = log2e_hi - yh
  75         movq $0x3de705fc2f000000,%rax
  76         pushq %rax
  77         fldl (%esp)
  78                 # fpu stack: 2^hi x lo xh xl yl
  79                 # lo += xh*yl + xl*yl
  80         fmul %st, %st(2)
  81         fmulp %st, %st(1)
  82         fxch %st(2)
  83         faddp
  84         faddp
  85                 # log2e_lo
  86         movq $0xbfbe,%rax
  87         pushq %rax
  88         movq $0x82f0025f2dc582ee,%rax
  89         pushq %rax
  90         fldt (%esp)
  91         add $40,%esp
  92                 # fpu stack: 2^hi x lo log2e_lo
  93                 # lo += log2e_lo*x
  94                 # return 2^hi + 2^hi (2^lo - 1)
  95         fmulp %st, %st(2)
  96         faddp
  97         f2xm1
  98         fmul %st(1), %st
  99         faddp
 100 1:      add $48, %esp
 101         ret