Re: [lmi] use boost::math:expm1() and log1p()

[Greg Chicares]

On 2009-04-05 19:20Z, Václav Slavík wrote:
> 
> On Apr 4, 2009, at 20:26, Greg Chicares wrote:
>> Before we consider switching to a different boost version or using
>> boost's build system, let's measure the speed and accuracy of these
>> functions. Since you've already built this, could I ask you to run
>> the math_functors_test$(EXEEXT) unit test and post the results?

Accuracy is fine, and speed seems to be unaffected as long as the
transcendental functions provided by the system are used.

> Before:
> 
>    Speed test: pow
>      1.763e-06 s =       1763 ns, mean of 5672 iterations
>    Speed test: expm1l
>      1.587e-06 s =       1587 ns, mean of 6300 iterations

Could I ask what compiler (and version) you're using? I get:

  Speed test: pow
    4.920e-006 s =       4920 ns, mean of 2033 iterations
  Speed test: expm1l
    1.078e-006 s =       1078 ns, mean of 9273 iterations

indicating that the pow() family is really slow for me. I think I may
be using an old libmingwex version whose pow() is slower but more
accurate than msvc's. But that doesn't affect lmi's speed as long as
we instead use log1p() and expm1() in compound-interest calculations,
as we generally (but, defectively, not everywhere) do, because our
other uses of pow() are rare and trivial.

> After:
> 
>    Speed test: pow
>      1.289e-06 s =       1289 ns, mean of 7758 iterations
>    Speed test: expm1
>       1.553e-06 s =       1553 ns, mean of 6439 iterations

Might your compiler treat pow() as an intrinsic? I kind of wonder how
anything could be faster than the 'expm1.c' implementation, because
AFAIK powers require evaluating both an exponential and a logarithmic
function on the x86. This comment in 'include/math.h':

  /* 7.12.6.3 The expm1 functions */
  /* TODO: These could be inlined */

may be part of the explanation. But I'd expect a pow() implementation
to contain at least as many special cases as 'expm1.c' so that it can
handle EDOM and ERANGE.

This is interesting, too:

>      long double precision, expm1 and log1p
>        -0.0039841072324612947578
[...]
>      double precision, expm1 and log1p
>        -0.0039841072324612947578

In this "After" case, switching the hardware to use a 53-bit mantissa
doesn't affect the result, though it did in the "Before" case:

       00000000001111111111
       12345678901234567890
                     ^ 15 == DBL_DIG
  -0.0039841072324612938904 Before
  -0.0039841072324612947578 After

so maybe an intermediate store-as-double is elided. But of course I
don't object to any increase in accuracy.