Proposal: Add log1p and expm1 to GHC.Float.Floating

[Henning Thielemann]

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.