fix exp asm

exponents (base 2) near 16383 were broken due to (1) wrong cutoff, and
(2) inability to fit the necessary range of scalings into a long
double value.

as a solution, we fall back to using frndint/fscale for insanely large
exponents, and also have to special-case infinities here to avoid
inf-inf generating nan.

thankfully the costly code never runs in normal usage cases.

src/math/i386/exp.s

@@ -68,21 +68,19 @@ exp:
 .type exp2,@function
 exp2:
 	fldl 4(%esp)
-1:	mov $0x47000000,%eax
+1:	pushl $0x467ff000
-	push %eax
 	flds (%esp)
-	shl $7,%eax
+	xorl %eax,%eax
-	push %eax
+	pushl $0x80000000
-	add %eax,%eax
 	push %eax
 	fld %st(1)
 	fabs
 	fucom %st(1)
 	fnstsw
+	fstp %st(0)
+	fstp %st(0)
 	sahf
 	ja 2f
-	fstp %st(0)
-	fstp %st(0)
 	fld %st(0)
 	fistpl 8(%esp)
 	fildl 8(%esp)
@@ -99,22 +97,23 @@ exp2:
 	add $12,%esp
 	ret
 
-2:	fstp %st(0)
+2:	fld %st(0)
+	fstpt (%esp)
+	mov 9(%esp),%ah
+	and $0x7f,%ah
+	cmp $0x7f,%ah
+	jne 1f
+	decb 9(%esp)
 	fstp %st(0)
-	fsts 8(%esp)
-	mov 8(%esp),%eax
-	lea (%eax,%eax),%ecx
-	cmp $0xff000000,%ecx
-	ja 2f
-	fstp %st(0)
-	xor %ecx,%ecx
-	inc %ecx
-	add %eax,%eax
-	jc 1f
-	mov $0x7ffe,%ecx
-1:	mov %ecx,8(%esp)
 	fldt (%esp)
-	fld %st(0)
+1:	fld %st(0)
-	fmulp
+	frndint
-2:	add $12,%esp
+	fxch %st(1)
+	fsub %st(1)
+	f2xm1
+	fld1
+	faddp
+	fscale
+	fstp %st(1)
+	add $12,%esp
 	ret