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
  8expl:
  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
 261:	faddp
 27	ret
 28
 29		# should be 0x1.71547652b82fe178p0L == 0x3fff b8aa3b29 5c17f0bc
 30		# it will be wrong on non-nearest rounding mode
 312:	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
1001:	add $48, %esp
101	ret