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? 