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

[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)

fldl 4(%esp) // x

fxam  // Is NaN or +-Inf?

fstsw %ax

movb $0x45, %ch

andb %ah, %ch

cmpb $0x40, %ch

je 3f // If +-0, jump.

#ifdef PIC

LOAD_PIC_REG (dx)

#endif

cmpb $0x05, %ch

je 2f // If +-Inf, jump.


 fldt MO(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)

fldl MO(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)

fsubrl MO(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)

ret


2: testl $0x200, %eax // Test sign.

jz 3f // If positive, jump.

fstp %st

fldl MO(minus1) // Set result to -1.0.

3: ret

I don't think we want to get into replicating that logic. =)

-Edward



On Thu, Apr 17, 2014 at 2:48 PM, Henning Thielemann <
schlepptop at henning-thielemann.de> wrote:

> Am 17.04.2014 19:15, schrieb Edward Kmett:
>
>
>  log1p and expm1 are C standard library functions that are important for
>> work with exponentials and logarithms.
>>
>> I propose adding them to the Floating class where it is defined in
>> GHC.Float.
>>
>
> No question, these functions are useful. But I think there should be two
> proposals:
>
> 1) Add log1pFloat, log1pDouble, expm1Float, expm1Double to GHC.Float
> 2) Extend Floating class with log1p and expm1 methods.
>
> I think the first item is unproblematic since it is a simple addition.
> Since FPUs sometimes directly implement log1p and expm1 functions, I wonder
> whether GHC also should support the according machine instructions. E.g.
> x86 has F2XM1 and FYL2XP1 and good old MC68882 had FETOXM1 and FLOGNP1.
>
> The second item means to alter the Floating class which affects all custom
> Floating instances. I think one should add default implementations. They
> don't have an numerical advantage but they save programmers from code
> breakage.
>
>
>
>  We do not have to export these from Prelude. My knee-jerk reaction is
>> that we probably shouldn't.
>>
>
> Not exporting them from Prelude still means to export them from GHC.Float,
> right? I mean, users must be able to implement these methods for custom
> types like extended precision floating point numbers as provided by libqd.
> But there should also be a non-GHC module that exports the full Floating
> class.
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.haskell.org/pipermail/libraries/attachments/20140417/91b37ade/attachment.html>