Am 17.04.2014 23:02, schrieb Edward Kmett:
I'm not advocating that we use the opcodes themselves.
That is pretty much tantamount to numerical precision suicide as most of them have pretty arcane restrictions on their valid input ranges that vary by platform and expect you to play games with scaling or newton-raphson as well. You need a fair bit of wrapping around the available floating point opcodes.
e.g. in glibc on x86 expm1 looks like:
.text
ENTRY(__expm1)
fldl4(%esp)// x
fxam// Is NaN or +-Inf?
fstsw%ax
movb$0x45, %ch
andb%ah, %ch
cmpb$0x40, %ch
je3f// If +-0, jump.
#ifdefPIC
LOAD_PIC_REG (dx)
#endif
cmpb$0x05, %ch
je2f// If +-Inf, jump.
fldtMO(l2e)// log2(e) : x
fmulp// log2(e)*x
fld%st// log2(e)*x : log2(e)*x
frndint// int(log2(e)*x) : log2(e)*x
fsubr%st, %st(1)// int(log2(e)*x) : fract(log2(e)*x)
fxch// fract(log2(e)*x) : int(log2(e)*x)
f2xm1// 2^fract(log2(e)*x)-1 : int(log2(e)*x)
fscale// 2^(log2(e)*x)-2^int(log2(e)*x) : int(log2(e)*x)
fxch// int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
fldlMO(one)// 1 : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
fscale// 2^int(log2(e)*x) : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
fsubrlMO(one)// 1-2^int(log2(e)*x) : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
fstp%st(1)// 1-2^int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
fsubrp%st, %st(1)// 2^(log2(e)*x)
I guess this is 2^(log2(e)*x)-1, otherwise the effort was unnecessary. :-) But I see that there is enough code to justify a sub-routine.