Re: [Mingw-users] expm1()

[Greg C. <chi...@co...>]

On 2005-01-27 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 gcc-3.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 -fomit-frame-pointer, 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 right--that'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 minimal-width IEC 60559
   double-extended 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 twenty-one 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 copyright-free values.
Any thoughts? The 'enquire' program doesn't give us these constants.

I tried using mingw gcc-2.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?