1/* 2 * Written by: 3 * J.T. Conklin (jtc@netbsd.org) 4 * Public domain. 5 */ 6 7#include <i387/bsd_asm.h> 8 9/* e^x = 2^(x * log2(e)) */ 10 11ENTRY(exp) 12 /* 13 * If x is +-Inf, then the subtraction would give Inf-Inf = NaN. 14 * Avoid this. Also avoid it if x is NaN for convenience. 15 */ 16 movl 8(%esp),%eax 17 andl $0x7fffffff,%eax 18 cmpl $0x7ff00000,%eax 19 jae x_Inf_or_NaN 20 21 fldl 4(%esp) 22 23 /* 24 * Extended precision is needed to reduce the maximum error from 25 * hundreds of ulps to less than 1 ulp. Switch to it if necessary. 26 * We may as well set the rounding mode to to-nearest and mask traps 27 * if we switch. 28 */ 29 fstcw 4(%esp) 30 movl 4(%esp),%eax 31 andl $0x0300,%eax 32 cmpl $0x0300,%eax /* RC == 0 && PC == 3? */ 33 je 1f /* jump if mode is good */ 34 movl $0x137f,8(%esp) 35 fldcw 8(%esp) 361: 37 fldl2e 38 fmulp /* x * log2(e) */ 39 fst %st(1) 40 frndint /* int(x * log2(e)) */ 41 fst %st(2) 42 fsubrp /* fract(x * log2(e)) */ 43 f2xm1 /* 2^(fract(x * log2(e))) - 1 */ 44 fld1 45 faddp /* 2^(fract(x * log2(e))) */ 46 fscale /* e^x */ 47 fstp %st(1) 48 je 1f 49 fldcw 4(%esp) 501: 51 ret 52 53x_Inf_or_NaN: 54 /* 55 * Return 0 if x is -Inf. Otherwise just return x; when x is Inf 56 * this gives Inf, and when x is a NaN this gives the same result 57 * as (x + x) (x quieted). 58 */ 59 cmpl $0xfff00000,8(%esp) 60 jne x_not_minus_Inf 61 cmpl $0,4(%esp) 62 jne x_not_minus_Inf 63 fldz 64 ret 65 66x_not_minus_Inf: 67 fldl 4(%esp) 68 ret 69END(exp) 70 71 // 72 73/* Enable stack protection */ 74#if defined(__linux__) && defined(__ELF__) 75.section .note.GNU-stack,"",%progbits 76#endif 77