Menu
▾
▴
Re: [Mingw-users] expm1()
|
From: Greg C. <chi...@co...> - 2005-01-20 19:24:17
|
On 2005-01-19 08:27 AM, Greg Chicares wrote:
> On 2005-01-19 04:01 AM, Danny Smith wrote:
>
>> Greg Chicares wrote:
>>
>>> log1p() would be much more useful if its companion expm1()
>>> [C99 7.12.6.3] were provided as well. [...]
>>
>> Off top of head, it would be quite easy to do
>>
>> if (x < small_value)
>> {
>> /* Taylor series expansion
>> return y such that additional term is smaller than machine eps */
>> }
>> else
>> return exp(x) - 1;
>>
>> Not very fast, but I'll bet I can find a good old Fortran IV routine
>> that has
>> been tested as accurate.
>
> The hardware already provides F2XM1, which is both fast and accurate.
>
> I'd provide a patch, except that it would take many months to get
> assignment papers signed by my employer. But maybe I can do other
> work that makes it trivial for someone else who doesn't have that
> problem.
[...]
> Maybe I could provide a narrative description of changes to
> mingw/runtime/mingwex/math/expl.c
> that would be easy to translate. The intel instruction set actually
> treats expm1() as the general case, and you have to undo the 'minus
> one' part to implement plain exp(). It might be as simple as [...]
...yes, it is that simple, and here's a clean-room 'not-a-patch'
that should make it easy to implement without copyright or
license issues:
// Make a copy of mingwex/math/expl.c .
//
// Change the function name (I happened to use 'NEW_FUNCTION').
//
// The original function uses an asm instruction that already
// calculates expm1(), which is what we want; then it adds one,
// which is not what we want, so get rid of the code that adds one.
//
// The inline-asm function should be used only within its domain of
// (-.5*M_LN2, +.5*M_LN2)
// A cover function ought to test whether its argument is outside
// that domain and, if it is, simply call exp() and return its
// result minus one.
And here's an optional test suite that doesn't have to enter any
MinGW sources, in case you can't rely on my copyright disclaimer.
// Test suite for a copyright-free expm1() derived from mingwex/math/expl.c .
// Written by Gregory W. Chicares 2005-01-20 and placed in the public domain.
// Compiled and tested with '-ggdb -W -Wall -std=gnu99 -pedantic'
// using MinGW gcc-3.2.3 and -3.4.2 .
// [define your NEW_FUNCTION here for testing]
#include <float.h>
#include <math.h>
#include <stdio.h>
// Daily-interest tests assume 365 days in a year, which is a
// typical practice in the financial industry.
int const n = 365;
void test0(double i)
{
double powi_daily = pow(1.0 + i, 1.0 / n);
double i_ = pow(powi_daily, n) - 1.0;
printf("Daily interest compounding with pow():\n");
printf("%40.20f %s\n", i , " = annual effective rate");
printf("%40.20f %s\n", powi_daily, " = daily effective rate");
printf("%40.20f %s\n", i_, " = daily rate converted back to annual");
printf("%40.20f %s\n", i_ - i, " = absolute error");
printf("%40.0e %s\n", (i_ - i) / i, " = relative error");
printf("\n");
}
void test1(double i)
{
double logi_daily = log(1.0 + i) / n;
double i_ = exp(n * logi_daily) - 1.0;
printf("Daily interest compounding with log() and exp():\n");
printf("%40.20f %s\n", i , " = annual effective rate");
printf("%40.20f %s\n", logi_daily, " = daily effective rate");
printf("%40.20f %s\n", i_, " = daily rate converted back to annual");
printf("%40.20f %s\n", i_ - i, " = absolute error");
printf("%40.0e %s\n", (i_ - i) / i, " = relative error");
printf("\n");
}
void test2(long double i)
{
double log1p_i_daily = log1pl(i) / n;
double i_ = NEW_FUNCTION(n * log1p_i_daily);
printf("Daily interest compounding with log1p() and expm1():\n");
printf("%40.20f %s\n", (double)i , " = annual effective rate");
printf("%40.20f %s\n", log1p_i_daily , " = daily effective rate");
printf("%40.20f %s\n", i_, " = daily rate converted back to annual");
printf("%40.20f %s\n", i_ - (double)i, " = absolute error");
printf("%40.0e %s\n", (i_ - (double)i) / (double)i, " = relative error");
printf("\n");
}
void test3(long double d)
{
long double d0 = NEW_FUNCTION(d);
long double d1 = expl(d) - 1.0;
printf("%25.20g %25.20g %20.0g\n", (double)d, (double)d0, (double)(d1 - d0));
}
int main()
{
// Test three methods of daily interest compounding to demonstrate
// why expm1() is important.
test0(.0001);
test1(.0001);
test2(.0001);
test0(.0000001);
test1(.0000001);
test2(.0000001);
// Test domain. An appropriate cover function would call the inline
// asm routine only if the argument is within its domain; otherwise,
// calling exp() and subtracting one is good enough.
//
// The modified mingwex function is accurate within (-.5*ln2, +.5*ln2).
// That's the domain where accuracy is most greatly improved compared
// to 'exp() - 1.0'. Probably that domain could be doubled by setting
// the hardware rounding flags before each FRNDINT, but that would
// make the function slower while achieving little.
// MinGW's include/math.h defines
// #define M_LN2 0.69314718055994530942
// Calculate half of that: 0.34657359027997265471
// and actually use 0.346573590279972
// (manually truncated *down* to IEC 60559 DBL_DIG digits) for safety.
// The actual result of exp() or expm1() for arguments of +/- DBL_MIN
// may not be very useful, but it's important to demonstrate that
// no bizarre error occurs at those extrema.
printf("Test domain of inline asm function.\n\n");
printf("%25s %25s %20s\n", "d", "expm1(d)", "exp(d) - expm1(d)");
test3( 1.0);
test3( 0.346573590279973); // Outside of domain.
test3( 0.346573590279972); // Within domain.
test3( 0.3);
test3( 0.2);
test3( 0.1);
test3( 0.0001);
test3( 0.0000001);
test3( 0.0000000001);
test3( DBL_MIN);
test3( 0.0);
test3(-DBL_MIN);
test3(-0.0000000001);
test3(-0.0000001);
test3(-0.0001);
test3(-0.1);
test3(-0.2);
test3(-0.3);
test3(-0.346573590279972); // Within domain.
test3(-0.346573590279973); // Outside of domain.
test3( 1.0);
return 0;
}
|
