From: Greg Chicares <chicares@co...>  20050118 12:28:25

log1p() would be much more useful if its companion expm1() [C99 7.12.6.3] were provided as well. Is there a reason not to have expm1(), other than "no one has contributed it yet"? For instance, Conklin, credited as the original author of /mingw/runtime/mingwex/math/log1pl.S?rev=1.1&view=markup said in an earlier version http://www.klid.dk/pub/gnu/gnuutils/src/libm/src/s_log1p.s /* * Since the fyl2xp1 instruction has such a limited range: * (1  (sqrt(2) / 2)) <= x <= sqrt(2)  1 * it's not worth trying to use it. */ but apparently later realized its worth. Did something like that happen for expm1()? Hmm...I find Conklin's log1p() in the netbsd src/lib/libm sources, but no corresponding expm1(). I realize expm1() will be an intrinsic in gcc4.0 http://www.gnu.org/software/gcc/gcc4.0/changes.html but apparently only with 'ffastmath', so there's still a use for an expm1() that conforms to the language and numeric standards. 
From: Danny Smith <dannysmith@cl...>  20050119 09:03:20

Greg Chicares wrote: > log1p() would be much more useful if its companion expm1() > [C99 7.12.6.3] were provided as well. Is there a reason not > to have expm1(), other than "no one has contributed it yet"? > > For instance, Conklin, credited as the original author of > /mingw/runtime/mingwex/math/log1pl.S?rev=1.1&view=markup > said in an earlier version > http://www.klid.dk/pub/gnu/gnuutils/src/libm/src/s_log1p.s > > /* > * Since the fyl2xp1 instruction has such a limited range: > * (1  (sqrt(2) / 2)) <= x <= sqrt(2)  1 > * it's not worth trying to use it. > */ > > but apparently later realized its worth. Did something like > that happen for expm1()? Hmm...I find Conklin's log1p() in the > netbsd src/lib/libm sources, but no corresponding expm1(). > > I realize expm1() will be an intrinsic in gcc4.0 > http://www.gnu.org/software/gcc/gcc4.0/changes.html > but apparently only with 'ffastmath', so there's still > a use for an expm1() that conforms to the language and > numeric standards. > 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. Danny > > >  > The SF.Net email is sponsored by: Beat the postholiday blues > Get a FREE limited edition SourceForge.net tshirt from ThinkGeek. > It's fun and FREE  well, almost....http://www.thinkgeek.com/sfshirt > _______________________________________________ > MinGWusers mailing list > MinGWusers@... > > You may change your MinGW Account Options or unsubscribe at: > https://lists.sourceforge.net/lists/listinfo/mingwusers 
From: Greg Chicares <chicares@co...>  20050119 13:27:41

On 20050119 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. I see plenty of expm1() implementations online, but some don't meet the FSF definition of 'free'; and those that do are GPL or LGPL, which I believe MinGW won't use. 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 eliminating the lines "fld1\n\t" /* 4 1.0 */ "faddp\n\t" /* 3 2^(fract(x * log2(e))) */ but let me look into that. 
From: Greg Chicares <chicares@co...>  20050120 19:24:17

On 20050119 08:27 AM, Greg Chicares wrote: > On 20050119 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 cleanroom 'notapatch' 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 inlineasm 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 copyrightfree expm1() derived from mingwex/math/expl.c . // Written by Gregory W. Chicares 20050120 and placed in the public domain. // Compiled and tested with 'ggdb W Wall std=gnu99 pedantic' // using MinGW gcc3.2.3 and 3.4.2 . // [define your NEW_FUNCTION here for testing] #include <float.h> #include <math.h> #include <stdio.h> // Dailyinterest 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; } 
From: Danny Smith <dannysmith@cl...>  20050126 07:50:23

Greg Chicares wrote: >> 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 cleanroom 'notapatch' > 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 inlineasm 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. > If we restrict domain of an our asm __expm1 to (.5*M_LN2, +.5*M_LN2), than much of the code in exp() is really unnecessary. We don't need to worry about the integral part of the separaton x log2(e) = i + f, f <= .5 Danny . 
From: Greg Chicares <chicares@co...>  20050126 08:59:40

On 20050126 02:48 AM, Danny Smith wrote: > Greg Chicares wrote: > >>>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 cleanroom 'notapatch' >>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 inlineasm 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. > > If we restrict domain of an our asm __expm1 to (.5*M_LN2, +.5*M_LN2), than much > of the code > in exp() is really unnecessary. We don't need to worry about the integral part > of the separaton > > x log2(e) = i + f, f <= .5 You're right. I was trying to change as little as possible in order to avoid FSF issues, but this evening I remembered that mingwex isn't part of the GNU project. I'll rewrite patch 1109674 to reflect your suggestion, which will save some cycles as well as making things clearer. I think the domain can be made twice as broad, too. 
From: Greg Chicares <chicares@co...>  20050126 19:49:51

On 20050126 02:48 AM, Danny Smith wrote: > > If we restrict domain of an our asm __expm1 to (.5*M_LN2, +.5*M_LN2), than much > of the code > in exp() is really unnecessary. We don't need to worry about the integral part > of the separaton > > x log2(e) = i + f, f <= .5 How about this? static long double expm1l (long double x) { if (M_LN2 < x && x < M_LN2) { x *= M_LOG2E; __asm__("f2xm1\n\t" : "=t" (x) : "0" (x)); return x; } else return expl(x)  1.0L; } It passes my unit tests with MinGW gcc3.4.2 . The domain is (M_LN2, +M_LN2) because F2XM1's input is constrained to (1, +1). Maybe intel means [1, +1], but their documentation is unclear as to whether the interval is open or not, so it's safer to assume it's open. 
From: Danny Smith <dannysmith@cl...>  20050127 08:31:17

Greg Chicares wrote: > How about this? > > static long double expm1l (long double x) > { > if (M_LN2 < x && x < M_LN2) > { > x *= M_LOG2E; > __asm__("f2xm1\n\t" : "=t" (x) : "0" (x)); > return x; > } > else > return expl(x)  1.0L; > } > > It passes my unit tests with MinGW gcc3.4.2 . > > The domain is (M_LN2, +M_LN2) because F2XM1's input is > constrained to (1, +1). Maybe intel means [1, +1], but > their documentation is unclear as to whether the interval > is open or not, so it's safer to assume it's open. > Simple That works in my tests too. I also tried out using "fldl2e\n\t" "fmul %st(1)\n\t" in the asm to replace x *= M_LOG2E; but at least with O2 fomitframepointer, I couldn't see any perfomance difference. (Using the M_LOG2E constant may be faster in loops if inline version is visible) BTW, do we need long double versions of the M_* defines to ensure enough precision? Will you update your patch on the Sourceforge patch tracker? Danny > > >  > This SF.Net email is sponsored by: IntelliVIEW  Interactive Reporting > Tool for open source databases. Create drag&drop reports. Save time > by over 75%! Publish reports on the web. Export to DOC, XLS, RTF, etc. > Download a FREE copy at http://www.intelliview.com/go/osdn_nl > _______________________________________________ > MinGWusers mailing list > MinGWusers@... > > You may change your MinGW Account Options or unsubscribe at: > https://lists.sourceforge.net/lists/listinfo/mingwusers 
From: Greg Chicares <chicares@co...>  20050128 20:16:55

On 20050127 03:29 AM, Danny Smith wrote: > Greg Chicares wrote: > >>How about this? >> >>static long double expm1l (long double x) >>{ >> if (M_LN2 < x && x < M_LN2) >> { >> x *= M_LOG2E; >> __asm__("f2xm1\n\t" : "=t" (x) : "0" (x)); >> return x; >> } >> else >> return expl(x)  1.0L; >>} >> >>It passes my unit tests with MinGW gcc3.4.2 . >> >>The domain is (M_LN2, +M_LN2) because F2XM1's input is >>constrained to (1, +1). Maybe intel means [1, +1], but >>their documentation is unclear as to whether the interval >>is open or not, so it's safer to assume it's open. >> > > Simple > That works in my tests too. Great. Thanks. > I also tried out using > > "fldl2e\n\t" > "fmul %st(1)\n\t" > > in the asm > to replace x *= M_LOG2E; > > but at least with O2 fomitframepointer, I couldn't see any perfomance > difference. (Using the M_LOG2E constant may be faster in loops if inline > version is visible) > > BTW, do we need long double versions of the M_* defines to ensure enough > precision? Gosh, you're rightthat's a real problem. We want DECIMAL_DIG digits in a long double context. That means 21 digits according to C99 5.2.4.2.2/13: "if the widest type were to use the minimalwidth IEC 60559 doubleextended format (64 bits of precision), then DECIMAL_DIG would be 21." In math.h the M_* defines have this comment: /* Traditional/XOPEN math constants (double precison) */ but actually they're given to twenty digits after the decimal point, so some have twentyone significant digits while others have only twenty: 12345678901234567890 #define M_LN2 0.69314718055994530942 /* not ideal */ #define M_LN10 2.30258509299404568402 /* ideal */ and they need to have 'L' at the end to preserve precision. Moshier's cephes/lldouble/constll.c gives long double LOGE2L = 0.6931471805599453094172321214581765680755001L; which, significantly for our purpose here, is less than the M_LN2 value in <math.h>. I assume cephes is acceptable because you've used it elsewhere in mingwex, but it gives values only for a few M_* constants: M_LN2 M_LOG2E M_PI M_PI_2 M_PI_4 Published values are easy to find, but not copyrightfree values. Any thoughts? The 'enquire' program doesn't give us these constants. I tried using mingw gcc2.95.x's ability to print long doubles to to generate them, e.g.: #include <iomanip> #include <iostream> int main() { volatile unsigned short int control_word = 0x037f; asm volatile("fldcw %0" : : "m" (control_word)); long double m_ln2; __asm__("fldln2\n\t" : "=t" (m_ln2) : "0" (m_ln2)); std::cout << std::setprecision(21); std::cout << m_ln2 << std::endl; } but got 000000000111111111122 123456789012345678901 0.693147180559945286227 /* program output */ 0.693147180559945309417... /* cephes */ ^ discrepancy after 16 digits so I don't trust that technique. Copyrighted sources like djgpp agree with Moshier here. And adding this to the program above long double x = (1000000000 * m_ln2)  (int)(1000000000 * m_ln2); std::cout << x << std::endl; gives the trailing digits as ...559945309418... which again agrees with Moshier's value. IIRC, the old 2.95 code limited setprecision() to something like eighteen. > Will you update your patch on the Sourceforge patch tracker? I will. First, though, let's try to solve the M_LN2 problem. I'd hate to hardcode a constant value in expm1() now, when what we really want is to use that macro, but with a value that's not overstated. Does mingwex have to be utterly copyright free? Paul Bristow proposed a boost library for math constants (rejected, but still in their archives); it uses the boost license here: http://www.boost.org/LICENSE_1_0.txt It seems more liberal than the BSD license in include/getopt.h; is that OK for mingwex? 
From: Danny Smith <dannysmith@cl...>  20050129 05:42:46

Greg Chicares wrote: > On 20050127 03:29 AM, Danny Smith wrote: > >> Greg Chicares wrote: >> > > I will. First, though, let's try to solve the M_LN2 problem. > I'd hate to hardcode a constant value in expm1() now, when what > we really want is to use that macro, but with a value that's > not overstated. On second thought, hardcoded values of the constants may be the way to go, particulary if we put inline versions of the functions in math.h. The M_* defines are !__STRICT_ANSI__, so we can't be sure that they will be visible (eg, if std=c99). Danny. 
From: Greg Chicares <chicares@co...>  20050129 16:46:31

On 20050129 00:41 AM, Danny Smith wrote: > Greg Chicares wrote: > >>On 20050127 03:29 AM, Danny Smith wrote: >> >>>Greg Chicares wrote: >> >>I will. First, though, let's try to solve the M_LN2 problem. >>I'd hate to hardcode a constant value in expm1() now, when what >>we really want is to use that macro, but with a value that's >>not overstated. > > On second thought, hardcoded values of the constants may be the way to go, > particulary if we put inline versions of the functions in math.h. The M_* > defines are !__STRICT_ANSI__, so we can't be sure that they will be visible (eg, > if std=c99). Good point. I've updated http://sourceforge.net/tracker/index.php?func=detail&aid=1109674&group_id=2435&atid=302435 and I believe it reflects all our discussions to date. Please take a look at the included test program. It shows differences of exactly one between expm1() and exp()1 for arguments extreme enough to return infinity or indefinite. C99 7.12.6.3/2 says "A range error occurs if x is too large." but 7.12.1/3 doesn't seem to say what an implementation should do when a range error occurs, so I guess anything goes for extremely high arguments. I'm not sure what should happen for an argument of INF. I get an indefinite for both expm1(INF) and exp(INF)1. 7.12.1/5 seems to allow that if we regard indefinites as representable "without extraordinary roundoff error". 
From: Earnie Boyd <earnie@us...>  20050129 17:59:54

<quote who="Greg Chicares"> > > Does mingwex have to be utterly copyright free? Paul Bristow > proposed a boost library for math constants (rejected, but > still in their archives); it uses the boost license here: > http://www.boost.org/LICENSE_1_0.txt > It seems more liberal than the BSD license in include/getopt.h; > is that OK for mingwex? > This license appears fine to me so if Danny's ok with it, go for it. You may wish to contact Paul Bristow and see if he'll release it to PD first though. Earnie  http://www.mingw.org http://sourceforge.net/projects/mingw https://sourceforge.net/donate/index.php?user_id=15438 
From: Danny Smith <dannysmith@cl...>  20050129 01:19:07

Greg Chicares wrote: > On 20050127 03:29 AM, Danny Smith wrote: >> BTW, do we need long double versions of the M_* defines to ensure enough >> precision? > > Gosh, you're rightthat's a real problem. > > We want DECIMAL_DIG digits in a long double context. That means > 21 digits according to C99 5.2.4.2.2/13: > <snip> > >> Will you update your patch on the Sourceforge patch tracker? > > I will. First, though, let's try to solve the M_LN2 problem. > I'd hate to hardcode a constant value in expm1() now, when what > we really want is to use that macro, but with a value that's > not overstated. > If we really want exact 64bit representation we could use float hex format, but IIRC 2.95 had a bug with hex format (S Moshier published a patch to fix the gcc bug, and I think I put it into some of the later 2.95.3x distros).. The hex float numbers could be generated using a union of shorts and long double. I could use the Cephes ldtoa code to get decimal values as well, and check against published values In any event, I think we should keep the OPENX defines and add new ones, like _M_EL, M_LOG2EL, etc. Danny 