MinGW - Minimalist GNU for Windows

Hi TSalm!

On Sat, Apr 30, 2011 at 5:55 PM, TSalm <TSalm@...> wrote:
> Hi,
>
> Is there away to catch an exception and display his stacktrace ?

I believe that the answer is no, if I understand your question.

Ordinary, built-in c++ exceptions do not take a snapshot of the
stack or build a stacktrace as they unwind the stack.  This is
by design, as c++ doesn't want to impose on you the cost of
building and maintaining this stacktrace.

The stack gets unwound between where you throw the exception
and where you catch it.  This can be inconvenient when debugging,
because if your code throws an exception deep down in some nested
function calls, seeing the stacktrace (at the time the exception was
thrown) can help make the bug apparent.  (I think java exceptions do
give you the stacktrace, and this can be very helpful.)

When I'm doing development and debugging, I generally put in assets.
These cause the program to core dump (or send control back to the
debugger) right at the point where the error is detected (i.e., the assert
fails), and then you can look at the state of the stack when the error
occurred (unless your bug has corrupted the stack...).

> TIA,
> -TSalm

Good luck.


K. Frank

On 2011-04-29 19:55Z, K. Frank wrote:
> 
> By the way, could someone confirm that mingw does use msvcrt for
> sqrt and pow?

http://cygwin.com/cgi-bin/cvsweb.cgi/src/winsup/mingw/mingwex/math/?cvsroot=src

Looks to me like sqrt(double) calls into msvcrt. But pow() is a
different story. Danny implemented it in libmingwex because of
some quality-of-implementation problem with msvcrt; then people
complained that the libmingwex version was slower, and I can't
remember where it wound up, but it's all in the archives.

> sqrt (x) and pow (x, 0.5) ought to give the same result (even if not required
> to by IEEE-754).

If they should give exactly the same result always, then the
first thing that needs to be fixed is any standard that allows
them to differ. The language standards are so much looser than
the numerical standard, even now, that I'd rather see C and C++
catch up to IEEE-754-1985 before asking them to go beyond.

> Browsing some gcc lists, I did see some comments that suggest that
> gcc (on linux) does try to transform pow (x, 0.5) to sqrt (x).  This would
> make sqrt and pow consistent (whether or not technically correct).

Not so many years ago (perhaps when msvcrt was written), it was
thought that correctly-rounded transcendental functions weren't
feasible in practice, so library authors had a looser attitude
toward exactness. If you have an ideal implementation of sqrt()
but a pow() that's only approximate, should you special-case
the power of exactly 0.5? If you do, then pow() and sqrt() will
be consistent, but then the powers
  x^(0.5*(1-epsilon))
  x^(0.5            )
  x^(0.5*(1+epsilon))
might not be monotone. Some would call that a bad tradeoff.

Hi Greg!

On Fri, Apr 29, 2011 at 7:38 PM, Greg Chicares <gchicares@...> wrote:
> On 2011-04-29 19:55Z, K. Frank wrote:
>>
>> By the way, could someone confirm that mingw does use msvcrt for
>> sqrt and pow?
>
> http://cygwin.com/cgi-bin/cvsweb.cgi/src/winsup/mingw/mingwex/math/?cvsroot=src
>
> Looks to me like sqrt(double) calls into msvcrt. But pow() is a
> different story. Danny implemented it in libmingwex because of
> some quality-of-implementation problem with msvcrt; then people
> complained that the libmingwex version was slower, and I can't
> remember where it wound up, but it's all in the archives.

Thanks for that -- good to know.  From the source link you gave:

   http://cygwin.com/cgi-bin/cvsweb.cgi/src/winsup/mingw/mingwex/math/pow.c?annotate=1.1.10.1&cvsroot=src

the code for pow (the regular double pow (double, double)) has these
comments about accuracy.

      28:  * ACCURACY:
      29:  *
      30:  *                      Relative error:
      31:  * arithmetic   domain     # trials      peak         rms
      32:  *    IEEE     -26,26       30000      4.2e-16      7.7e-17
      33:  *    DEC      -26,26       60000      4.8e-17      9.1e-18
      34:  * 1/26 < x < 26, with log(x) uniformly distributed.
      35:  * -26 < y < 26, y uniformly distributed.
      36:  *    IEEE     0,8700       30000      1.5e-14      2.1e-15
      37:  * 0.99 < x < 1.01, 0 < y < 8700, uniformly distributed.

This commentary suggest that pow is not "exact."

This would reasonably explain Arthur seeing pow not match sqrt.
But it still leaves the puzzle of why the different mingw versions
don't behave the same.

One minor puzzle:  I don't see anything in the pow code that would
let the compiler / library take advantage of hardware support for pow,
which I think has been available in the x86 architecture for quite
some time now.

>
>> sqrt (x) and pow (x, 0.5) ought to give the same result (even if not required
>> to by IEEE-754).
>
> If they should give exactly the same result always, then the
> first thing that needs to be fixed is any standard that allows
> them to differ.

Well, it's a question of engineering compromise whether they should always
agree.  In an ideal world, it would be nice if pow were "exact" (in the sense
of closet floating-point number).  Then if pow were "exact," pow and sqrt
should agree.  But there is expense in getting pow to be "exact," so it's
arguable that the language / IEEE standards should give implementers some
flexibility.  (And, as you point out below, if pow isn't exact, then forcing pow
and sqrt to agree could actually cause other problems.)

> The language standards are so much looser than
> the numerical standard, even now, that I'd rather see C and C++
> catch up to IEEE-754-1985 before asking them to go beyond.
>
>> Browsing some gcc lists, I did see some comments that suggest that
>> gcc (on linux) does try to transform pow (x, 0.5) to sqrt (x).  This would
>> make sqrt and pow consistent (whether or not technically correct).
>
> Not so many years ago (perhaps when msvcrt was written), it was
> thought that correctly-rounded transcendental functions weren't
> feasible in practice, so library authors had a looser attitude
> toward exactness. If you have an ideal implementation of sqrt()
> but a pow() that's only approximate, should you special-case
> the power of exactly 0.5? If you do, then pow() and sqrt() will
> be consistent, but then the powers
>  x^(0.5*(1-epsilon))
>  x^(0.5            )
>  x^(0.5*(1+epsilon))
> might not be monotone. Some would call that a bad tradeoff.

I actually agree with your concern.  The transformation of pow (x, 0.5)
to sqrt (x) that it sounded like gcc / linux was doing sounds like a
bad idea if it's being done this simply and naively.

Any speculation on why Arthur gets different results with different
mingw's?

Thanks for your perspective.


K. Frank

2011/4/30 Greg Chicares <gchicares@...>:
> On 2011-04-29 19:55Z, K. Frank wrote:
>>
>> By the way, could someone confirm that mingw does use msvcrt for
>> sqrt and pow?
>
> http://cygwin.com/cgi-bin/cvsweb.cgi/src/winsup/mingw/mingwex/math/?cvsroot=src
>
> Looks to me like sqrt(double) calls into msvcrt. But pow() is a
> different story. Danny implemented it in libmingwex because of
> some quality-of-implementation problem with msvcrt; then people
> complained that the libmingwex version was slower, and I can't
> remember where it wound up, but it's all in the archives.
>
>> sqrt (x) and pow (x, 0.5) ought to give the same result (even if not required
>> to by IEEE-754).
>
> If they should give exactly the same result always, then the
> first thing that needs to be fixed is any standard that allows
> them to differ. The language standards are so much looser than
> the numerical standard, even now, that I'd rather see C and C++
> catch up to IEEE-754-1985 before asking them to go beyond.
>
>> Browsing some gcc lists, I did see some comments that suggest that
>> gcc (on linux) does try to transform pow (x, 0.5) to sqrt (x).  This would
>> make sqrt and pow consistent (whether or not technically correct).
>
> Not so many years ago (perhaps when msvcrt was written), it was
> thought that correctly-rounded transcendental functions weren't
> feasible in practice, so library authors had a looser attitude
> toward exactness. If you have an ideal implementation of sqrt()
> but a pow() that's only approximate, should you special-case
> the power of exactly 0.5? If you do, then pow() and sqrt() will
> be consistent, but then the powers
>  x^(0.5*(1-epsilon))
>  x^(0.5            )
>  x^(0.5*(1+epsilon))
> might not be monotone. Some would call that a bad tradeoff.

Well, I spent some efforts into log, exp, and pow implementation. Also
some other base-math functions I improved in mingw-w64's math library
in two terms. a) Make those functions behave as ISO-C99 specification
tells in corner-cases. and b) Improve accuracy in terms of gcc's
internal used gmp. As we noticed here result differences.

The following example illustrates the issues pretty well:

#include <stdio.h>
#include <math.h>

double abc[3][2] = { { 2.2, 3.1 }, { 2.0, 3.0 }, { 4.0, 0.5 } };
long double ab[3][2] = { { 2.2L, 3.1L }, { 2.0L, 3.0L }, { 4.0L, 0.5L } };
int main()
{
  double r[3];
  long double rl[3];
  int i;
  for (i = 0; i < 3; i++)
  __mingw_printf ("%.20F^%.20F = %.20F\n", abc[i][0], abc[i][1],
         (r[i] = pow (abc[i][0], abc[i][1])));
  __mingw_printf ("%.20F %.20F %.20F\n", pow (2.2, 3.1), pow (2.0,
3.0), pow (4.0, 0.5));
  r[0] -= pow (2.2, 3.1);
  r[1] -= pow (2.0, 3.0);
  r[2] -= pow (4.0, 0.5);
  __mingw_printf ("%.20F %.20F %.20F\n", r[0], r[1], r[2]);

  for (i = 0; i < 3; i++)
  __mingw_printf ("%.20LF^%.20LF = %.20LF\n", ab[i][0], ab[i][1],
         (rl[i] = powl (ab[i][0], ab[i][1])));
  __mingw_printf ("%.20LF %.20LF %.20LF\n",
                         powl (2.2L, 3.1L), powl (2.0L, 3.0L) , powl
(4.0L, 0.5L));
  rl[0] -= powl (2.2L, 3.1L);
  rl[1] -= powl (2.0L, 3.0L);
  rl[2] -= powl (4.0L, 0.5L);
  __mingw_printf ("%.20LF %.20LF %.20LF\n", rl[0], rl[1], rl[2]);
  return 1;
}

The interesting issue is here that gcc uses for constant
math-calculations gmp-library to get result. For more complex
variants, it uses crt's math.  By this it is important for a gcc based
runtime to be in the IEEE 754 floating point compatible to gmp's
results.

Btw new pow, log, and exp variant are slightly faster then variant in
msvcrt, but this is more or less just a side-effect.

I admit that this might be not that obvious to users, why results here
are different, but for a gcc based toolchain we need to play nice with
gcc's internal assumptions.

Regards,
Kai

2011/4/30 Kai Tietz <ktietz70@...>:
> 2011/4/30 Greg Chicares <gchicares@...>:
>> On 2011-04-29 19:55Z, K. Frank wrote:
>>>
>>> By the way, could someone confirm that mingw does use msvcrt for
>>> sqrt and pow?
>>
>> http://cygwin.com/cgi-bin/cvsweb.cgi/src/winsup/mingw/mingwex/math/?cvsroot=src
>>
>> Looks to me like sqrt(double) calls into msvcrt. But pow() is a
>> different story. Danny implemented it in libmingwex because of
>> some quality-of-implementation problem with msvcrt; then people
>> complained that the libmingwex version was slower, and I can't
>> remember where it wound up, but it's all in the archives.
>>
>>> sqrt (x) and pow (x, 0.5) ought to give the same result (even if not required
>>> to by IEEE-754).
>>
>> If they should give exactly the same result always, then the
>> first thing that needs to be fixed is any standard that allows
>> them to differ. The language standards are so much looser than
>> the numerical standard, even now, that I'd rather see C and C++
>> catch up to IEEE-754-1985 before asking them to go beyond.
>>
>>> Browsing some gcc lists, I did see some comments that suggest that
>>> gcc (on linux) does try to transform pow (x, 0.5) to sqrt (x).  This would
>>> make sqrt and pow consistent (whether or not technically correct).
>>
>> Not so many years ago (perhaps when msvcrt was written), it was
>> thought that correctly-rounded transcendental functions weren't
>> feasible in practice, so library authors had a looser attitude
>> toward exactness. If you have an ideal implementation of sqrt()
>> but a pow() that's only approximate, should you special-case
>> the power of exactly 0.5? If you do, then pow() and sqrt() will
>> be consistent, but then the powers
>>  x^(0.5*(1-epsilon))
>>  x^(0.5            )
>>  x^(0.5*(1+epsilon))
>> might not be monotone. Some would call that a bad tradeoff.
>
> Well, I spent some efforts into log, exp, and pow implementation. Also
> some other base-math functions I improved in mingw-w64's math library
> in two terms. a) Make those functions behave as ISO-C99 specification
> tells in corner-cases. and b) Improve accuracy in terms of gcc's
> internal used gmp. As we noticed here result differences.
>
> The following example illustrates the issues pretty well:
>
> #include <stdio.h>
> #include <math.h>
>
> double abc[3][2] = { { 2.2, 3.1 }, { 2.0, 3.0 }, { 4.0, 0.5 } };
> long double ab[3][2] = { { 2.2L, 3.1L }, { 2.0L, 3.0L }, { 4.0L, 0.5L } };
> int main()
> {
>  double r[3];
>  long double rl[3];
>  int i;
>  for (i = 0; i < 3; i++)
>  __mingw_printf ("%.20F^%.20F = %.20F\n", abc[i][0], abc[i][1],
>         (r[i] = pow (abc[i][0], abc[i][1])));
>  __mingw_printf ("%.20F %.20F %.20F\n", pow (2.2, 3.1), pow (2.0,
> 3.0), pow (4.0, 0.5));
>  r[0] -= pow (2.2, 3.1);
>  r[1] -= pow (2.0, 3.0);
>  r[2] -= pow (4.0, 0.5);
>  __mingw_printf ("%.20F %.20F %.20F\n", r[0], r[1], r[2]);
>
>  for (i = 0; i < 3; i++)
>  __mingw_printf ("%.20LF^%.20LF = %.20LF\n", ab[i][0], ab[i][1],
>         (rl[i] = powl (ab[i][0], ab[i][1])));
>  __mingw_printf ("%.20LF %.20LF %.20LF\n",
>                         powl (2.2L, 3.1L), powl (2.0L, 3.0L) , powl
> (4.0L, 0.5L));
>  rl[0] -= powl (2.2L, 3.1L);
>  rl[1] -= powl (2.0L, 3.0L);
>  rl[2] -= powl (4.0L, 0.5L);
>  __mingw_printf ("%.20LF %.20LF %.20LF\n", rl[0], rl[1], rl[2]);
>  return 1;
> }
>
> The interesting issue is here that gcc uses for constant
> math-calculations gmp-library to get result. For more complex
> variants, it uses crt's math.  By this it is important for a gcc based
> runtime to be in the IEEE 754 floating point compatible to gmp's
> results.
>
> Btw new pow, log, and exp variant are slightly faster then variant in
> msvcrt, but this is more or less just a side-effect.
>
> I admit that this might be not that obvious to users, why results here
> are different, but for a gcc based toolchain we need to play nice with
> gcc's internal assumptions.
>
> Regards,
> Kai
>


On closer look for long double, we can have rouding issues for long
double precission. I see that powl (0.01L, 0.5L) we have a rounding
difference to sqrtl (0.01L). Interestingly not for float or double
types.  So we we might need to special-case in pow the case for y ==
0.5 to solve this.

Kai

My comments don't seem to be making it into the thread this week, possibly
because of a switch between my reply-to email address and another.  So, if
it gets lost, perhaps nightstrike can intervene, or echo the information.
Perhaps not.

The sqrt() function can be implemented with excellent speed to within a bit
of full mantissa accuracy using iterations of Newton's method, as is
classically done in libraries.  Not so with ln() or exp(), which are
classically done by extracting the exponent and using the reduced-range
mantissa in a modified Chebychev approximation.  Arithmetic co-processors do
sqrt(), ln(), and exp() internally with 80-bit floating point, then store in
32, 64, or 80 bit floating point format.  So, the co-processors can get full
mantissa accuracy with 32-bit and 64-bit sqrt() and pow(x,0.5).  So, the
problem is non-existent for Intel/PPC and most other common workstation
processors.  Yet we are looking at software transcendental function
libraries.

James K Beard

-----Original Message-----
From: Greg Chicares [mailto:gchicares@...] 
Sent: Friday, April 29, 2011 7:38 PM
To: MinGW Users List
Cc: mingw64
Subject: Re: [Mingw-w64-public] [Mingw-users] Math library discrepancies
that surprised me.

On 2011-04-29 19:55Z, K. Frank wrote:
> 
> By the way, could someone confirm that mingw does use msvcrt for
> sqrt and pow?

http://cygwin.com/cgi-bin/cvsweb.cgi/src/winsup/mingw/mingwex/math/?cvsroot=
src

Looks to me like sqrt(double) calls into msvcrt. But pow() is a
different story. Danny implemented it in libmingwex because of
some quality-of-implementation problem with msvcrt; then people
complained that the libmingwex version was slower, and I can't
remember where it wound up, but it's all in the archives.

> sqrt (x) and pow (x, 0.5) ought to give the same result (even if not
required
> to by IEEE-754).

If they should give exactly the same result always, then the
first thing that needs to be fixed is any standard that allows
them to differ. The language standards are so much looser than
the numerical standard, even now, that I'd rather see C and C++
catch up to IEEE-754-1985 before asking them to go beyond.

> Browsing some gcc lists, I did see some comments that suggest that
> gcc (on linux) does try to transform pow (x, 0.5) to sqrt (x).  This would
> make sqrt and pow consistent (whether or not technically correct).

Not so many years ago (perhaps when msvcrt was written), it was
thought that correctly-rounded transcendental functions weren't
feasible in practice, so library authors had a looser attitude
toward exactness. If you have an ideal implementation of sqrt()
but a pow() that's only approximate, should you special-case
the power of exactly 0.5? If you do, then pow() and sqrt() will
be consistent, but then the powers
  x^(0.5*(1-epsilon))
  x^(0.5            )
  x^(0.5*(1+epsilon))
might not be monotone. Some would call that a bad tradeoff.

----------------------------------------------------------------------------
--
WhatsUp Gold - Download Free Network Management Software
The most intuitive, comprehensive, and cost-effective network 
management toolset available today.  Delivers lowest initial 
acquisition cost and overall TCO of any competing solution.
http://p.sf.net/sfu/whatsupgold-sd
_______________________________________________
Mingw-w64-public mailing list
Mingw-w64-public@...
https://lists.sourceforge.net/lists/listinfo/mingw-w64-public

Hello Kai!

On Sat, Apr 30, 2011 at 5:51 AM, Kai Tietz <ktietz70@...> wrote:
> 2011/4/30 Kai Tietz <ktietz70@...>:
>> 2011/4/30 Greg Chicares <gchicares@...>:
>>> On 2011-04-29 19:55Z, K. Frank wrote:
>>>>
>>>> By the way, could someone confirm that mingw does use msvcrt for
>>>> sqrt and pow?
>>>
>>> http://cygwin.com/cgi-bin/cvsweb.cgi/src/winsup/mingw/mingwex/math/?cvsroot=src
>>>
>>> Looks to me like sqrt(double) calls into msvcrt. But pow() is a
>>> different story. Danny implemented it in libmingwex because of
>>> some quality-of-implementation problem with msvcrt; then people
>>> complained that the libmingwex version was slower, and I can't
>>> remember where it wound up, but it's all in the archives.
>>>
>>>> sqrt (x) and pow (x, 0.5) ought to give the same result (even if not required
>>>> to by IEEE-754).
>>> ...
>>> Not so many years ago (perhaps when msvcrt was written), it was
>>> thought that correctly-rounded transcendental functions weren't
>>> feasible in practice, so library authors had a looser attitude
>>> toward exactness. If you have an ideal implementation of sqrt()
>>> but a pow() that's only approximate, should you special-case
>>> the power of exactly 0.5? If you do, then pow() and sqrt() will
>>> be consistent, but then the powers
>>>  x^(0.5*(1-epsilon))
>>>  x^(0.5            )
>>>  x^(0.5*(1+epsilon))
>>> might not be monotone. Some would call that a bad tradeoff.
>>
>> Well, I spent some efforts into log, exp, and pow implementation. Also
>> some other base-math functions I improved in mingw-w64's math library
>> in two terms. a) Make those functions behave as ISO-C99 specification
>> tells in corner-cases. and b) Improve accuracy in terms of gcc's
>> internal used gmp. As we noticed here result differences.
>>
>> The following example illustrates the issues pretty well:
>>
>> #include <stdio.h>
>> #include <math.h>
>>
>> double abc[3][2] = { { 2.2, 3.1 }, { 2.0, 3.0 }, { 4.0, 0.5 } };
>> long double ab[3][2] = { { 2.2L, 3.1L }, { 2.0L, 3.0L }, { 4.0L, 0.5L } };
>> int main()
>> {
>>  double r[3];
>>  long double rl[3];
>>  int i;
>>  for (i = 0; i < 3; i++)
>>  __mingw_printf ("%.20F^%.20F = %.20F\n", abc[i][0], abc[i][1],
>>         (r[i] = pow (abc[i][0], abc[i][1])));
>>  __mingw_printf ("%.20F %.20F %.20F\n", pow (2.2, 3.1), pow (2.0,
>> 3.0), pow (4.0, 0.5));
>>  r[0] -= pow (2.2, 3.1);
>>  r[1] -= pow (2.0, 3.0);
>>  r[2] -= pow (4.0, 0.5);
>>  __mingw_printf ("%.20F %.20F %.20F\n", r[0], r[1], r[2]);
>>
>>  for (i = 0; i < 3; i++)
>>  __mingw_printf ("%.20LF^%.20LF = %.20LF\n", ab[i][0], ab[i][1],
>>         (rl[i] = powl (ab[i][0], ab[i][1])));
>>  __mingw_printf ("%.20LF %.20LF %.20LF\n",
>>                         powl (2.2L, 3.1L), powl (2.0L, 3.0L) , powl
>> (4.0L, 0.5L));
>>  rl[0] -= powl (2.2L, 3.1L);
>>  rl[1] -= powl (2.0L, 3.0L);
>>  rl[2] -= powl (4.0L, 0.5L);
>>  __mingw_printf ("%.20LF %.20LF %.20LF\n", rl[0], rl[1], rl[2]);
>>  return 1;
>> }

Do I understand your point correctly that you are comparing gcc's
compile-time and run-time values for pow?

For the record, when I run this on a mingw64 build:

   gcc (GCC) 4.5.2 20101002 (prerelease) [svn/rev.164902 - mingw-w64/oz]

I get:

   C:\>gcc -o pow_test pow_test.c

   C:\>pow_test
   2.20000000000000017764^3.10000000000000008882 = 11.52153412678571875460
   2.00000000000000000000^3.00000000000000000000 = 8.00000000000000000000
   4.00000000000000000000^0.50000000000000000000 = 2.00000000000000000000
   11.52153412678571697825 8.00000000000000000000 2.00000000000000000000
   0.00000000000000177636 0.00000000000000000000 0.00000000000000000000
   2.20000000000000000004^3.09999999999999999991 = 11.52153412678571414718
   2.00000000000000000000^3.00000000000000000000 = 8.00000000000000000000
   4.00000000000000000000^0.50000000000000000000 = 2.00000000000000000000
   11.52153412678571414718 8.00000000000000000000 2.00000000000000000000
   0.00000000000000000000 0.00000000000000000000 0.00000000000000000000

As I understand it, this means I am seeing a small difference between
compile and run-time evaluation only for powl (2.2, 3.1), i.e., the
long-double case.

>> The interesting issue is here that gcc uses for constant
>> math-calculations gmp-library to get result. For more complex
>> variants, it uses crt's math.  By this it is important for a gcc based
>> runtime to be in the IEEE 754 floating point compatible to gmp's
>> results.

Yes, that's a good point.  I hadn't considered the issue of compile-time
evaluation, and, ideally, it should agree with run-time evaluation, adding
yet another constraint to our wish list.

Furthermore, if we take the idea of a cross-compiler seriously, we
would also like the compile-time evaluations of a cross-compiler and
native compiler to agree.

Hmm...

So, we would like sqrt (x) and pow (x, 0.5) to agree.  We would like
compile-time and run-time evaluations to agree.  We would like
cross-compilers and native compilers to agree.

Taken at face value, this would seem to preclude using the platform
library and platform hardware to compute things like sqrt and pow
(unless we know they comply with some agreed-upon standard,
such as IEEE-754).

Just to make the reasoning explicit, if I perform my run-time evaluation
on platform A, e.g., using platform A's FPU, then a B-->A cross-compiler
running on platform B won't have platform A's FPU available to perform
consistent compile-time evaluations.

So you either give up some of the above-wished-for consistency, you
forgo the potential performance benefits of using the platform's library
and / or hardware, or you get everybody to use a common standard.

>> ...
>> I admit that this might be not that obvious to users, why results here
>> are different, but for a gcc based toolchain we need to play nice with
>> gcc's internal assumptions.
>>
>> Regards,
>> Kai
>
> On closer look for long double, we can have rouding issues for long
> double precission. I see that powl (0.01L, 0.5L) we have a rounding
> difference to sqrtl (0.01L). Interestingly not for float or double
> types.  So we we might need to special-case in pow the case for y ==
> 0.5 to solve this.

Be, however, cognizant of Greg's concern.  Forcing pow to agree with
sqrt just for y = 0.5 has the potential to make pow, in a sense, internally
inconsistent,  It's a tricky trade-off.  (Ideally, both pow and sqrt would
be "exact" and therefore agree with one another automatically.)

By the way, I also see a problem with run-time powl for powl (2.2L, 0.5L).
It disagrees slightly from sqrtl (2.2L) and from the compile-time value
for powl (2.2L, 0.5L).

I modified your test program to compare sqrt with pow (x, 0.5) (and sqrtl
with powl).  For the record:

   C:\>gcc --version
   gcc (GCC) 4.5.2 20101002 (prerelease) [svn/rev.164902 - mingw-w64/oz]

   C:\>gcc -o sqrt_pow sqrt_pow.c

   C:\>sqrt_pow
   1.04880884817015163080, 1.04880884817015163080, 0.00000000000000000000
   1.48323969741913264109, 1.48323969741913264109, 0.00000000000000000000
   1.81659021245849494619, 1.81659021245849494619, 0.00000000000000000000
   1.04880884817015163080 1.04880884817015163080
   1.48323969741913264109 1.48323969741913264109
   1.81659021245849494619 1.81659021245849494619
   0.00000000000000000000 0.00000000000000000000 0.00000000000000000000
   0.00000000000000000000 0.00000000000000000000 0.00000000000000000000
   1.04880884817015154699, 1.04880884817015154699, 0.00000000000000000000
   1.48323969741913258980, 1.48323969741913258970, -0.00000000000000000011
   1.81659021245849499921, 1.81659021245849499921, 0.00000000000000000000
   1.04880884817015154699 1.04880884817015154699
   1.48323969741913258980 1.48323969741913258980
   1.81659021245849499921 1.81659021245849499921
   0.00000000000000000000 0.00000000000000000000 0.00000000000000000000
   0.00000000000000000000 -0.00000000000000000011 0.00000000000000000000

where sqrt_pow.c is:

   #include <stdio.h>
   #include <math.h>

   double abc[3] = { 1.1, 2.2, 3.3 };
   long double abcl[3] = { 1.1L, 2.2L, 3.3L };
   int main()
   {
    double r1[3], r2[3], r3[3];
    long double rl1[3], rl2[3], rl3[3];
    int i;
    for (i = 0; i < 3; i++)
    __mingw_printf ("%.20F, %.20F, %.20F\n",
    (r1[i] = sqrt (abc[i])),
    (r2[i] = pow (abc[i], 0.5)),
    (r3[i] = pow (abc[i], 0.5) - sqrt (abc[i])));
    __mingw_printf ("%.20F %.20F\n", sqrt (1.1), pow (1.1, 0.5));
    __mingw_printf ("%.20F %.20F\n", sqrt (2.2), pow (2.2, 0.5));
    __mingw_printf ("%.20F %.20F\n", sqrt (3.3), pow (3.3, 0.5));
    r1[0] -= sqrt (1.1);
    r1[1] -= sqrt (2.2);
    r1[2] -= sqrt (3.3);
    __mingw_printf ("%.20F %.20F %.20F\n", r1[0], r1[1], r1[2]);
    r2[0] -= pow (1.1, 0.5);
    r2[1] -= pow (2.2, 0.5);
    r2[2] -= pow (3.3, 0.5);
    __mingw_printf ("%.20F %.20F %.20F\n", r2[0], r2[1], r2[2]);

    for (i = 0; i < 3; i++)
    __mingw_printf ("%.20LF, %.20LF, %.20LF\n",
    (rl1[i] = sqrtl (abcl[i])),
    (rl2[i] = powl (abcl[i], 0.5L)),
    (rl3[i] = powl (abcl[i], 0.5L) - sqrtl (abcl[i])));
    __mingw_printf ("%.20LF %.20LF\n", sqrtl (1.1L), powl (1.1L, 0.5L));
    __mingw_printf ("%.20LF %.20LF\n", sqrtl (2.2L), powl (2.2L, 0.5L));
    __mingw_printf ("%.20LF %.20LF\n", sqrtl (3.3L), powl (3.3L, 0.5L));
    rl1[0] -= sqrtl (1.1L);
    rl1[1] -= sqrtl (2.2L);
    rl1[2] -= sqrtl (3.3L);
    __mingw_printf ("%.20LF %.20LF %.20LF\n", rl1[0], rl1[1], rl1[2]);
    rl2[0] -= powl (1.1L, 0.5L);
    rl2[1] -= powl (2.2L, 0.5L);
    rl2[2] -= powl (3.3L, 0.5L);
    __mingw_printf ("%.20LF %.20LF %.20LF\n", rl2[0], rl2[1], rl2[2]);

    return 1;
   }

>
> Kai

Thanks, Kai, for your insight and explanations.


K. Frank

2011/4/30 K. Frank <kfrank29.c@...>:
> Hello Kai!
>
> On Sat, Apr 30, 2011 at 5:51 AM, Kai Tietz <ktietz70@...> wrote:
>> 2011/4/30 Kai Tietz <ktietz70@...>:
>>> 2011/4/30 Greg Chicares <gchicares@...>:
>>>> On 2011-04-29 19:55Z, K. Frank wrote:
>>>>>
>>>>> By the way, could someone confirm that mingw does use msvcrt for
>>>>> sqrt and pow?
>>>>
>>>> http://cygwin.com/cgi-bin/cvsweb.cgi/src/winsup/mingw/mingwex/math/?cvsroot=src
>>>>
>>>> Looks to me like sqrt(double) calls into msvcrt. But pow() is a
>>>> different story. Danny implemented it in libmingwex because of
>>>> some quality-of-implementation problem with msvcrt; then people
>>>> complained that the libmingwex version was slower, and I can't
>>>> remember where it wound up, but it's all in the archives.
>>>>
>>>>> sqrt (x) and pow (x, 0.5) ought to give the same result (even if not required
>>>>> to by IEEE-754).
>>>> ...
>>>> Not so many years ago (perhaps when msvcrt was written), it was
>>>> thought that correctly-rounded transcendental functions weren't
>>>> feasible in practice, so library authors had a looser attitude
>>>> toward exactness. If you have an ideal implementation of sqrt()
>>>> but a pow() that's only approximate, should you special-case
>>>> the power of exactly 0.5? If you do, then pow() and sqrt() will
>>>> be consistent, but then the powers
>>>>  x^(0.5*(1-epsilon))
>>>>  x^(0.5            )
>>>>  x^(0.5*(1+epsilon))
>>>> might not be monotone. Some would call that a bad tradeoff.
>>>
>>> Well, I spent some efforts into log, exp, and pow implementation. Also
>>> some other base-math functions I improved in mingw-w64's math library
>>> in two terms. a) Make those functions behave as ISO-C99 specification
>>> tells in corner-cases. and b) Improve accuracy in terms of gcc's
>>> internal used gmp. As we noticed here result differences.
>>>
>>> The following example illustrates the issues pretty well:
>>>
>>> #include <stdio.h>
>>> #include <math.h>
>>>
>>> double abc[3][2] = { { 2.2, 3.1 }, { 2.0, 3.0 }, { 4.0, 0.5 } };
>>> long double ab[3][2] = { { 2.2L, 3.1L }, { 2.0L, 3.0L }, { 4.0L, 0.5L } };
>>> int main()
>>> {
>>>  double r[3];
>>>  long double rl[3];
>>>  int i;
>>>  for (i = 0; i < 3; i++)
>>>  __mingw_printf ("%.20F^%.20F = %.20F\n", abc[i][0], abc[i][1],
>>>         (r[i] = pow (abc[i][0], abc[i][1])));
>>>  __mingw_printf ("%.20F %.20F %.20F\n", pow (2.2, 3.1), pow (2.0,
>>> 3.0), pow (4.0, 0.5));
>>>  r[0] -= pow (2.2, 3.1);
>>>  r[1] -= pow (2.0, 3.0);
>>>  r[2] -= pow (4.0, 0.5);
>>>  __mingw_printf ("%.20F %.20F %.20F\n", r[0], r[1], r[2]);
>>>
>>>  for (i = 0; i < 3; i++)
>>>  __mingw_printf ("%.20LF^%.20LF = %.20LF\n", ab[i][0], ab[i][1],
>>>         (rl[i] = powl (ab[i][0], ab[i][1])));
>>>  __mingw_printf ("%.20LF %.20LF %.20LF\n",
>>>                         powl (2.2L, 3.1L), powl (2.0L, 3.0L) , powl
>>> (4.0L, 0.5L));
>>>  rl[0] -= powl (2.2L, 3.1L);
>>>  rl[1] -= powl (2.0L, 3.0L);
>>>  rl[2] -= powl (4.0L, 0.5L);
>>>  __mingw_printf ("%.20LF %.20LF %.20LF\n", rl[0], rl[1], rl[2]);
>>>  return 1;
>>> }
>
> Do I understand your point correctly that you are comparing gcc's
> compile-time and run-time values for pow?
>
> For the record, when I run this on a mingw64 build:
>
>   gcc (GCC) 4.5.2 20101002 (prerelease) [svn/rev.164902 - mingw-w64/oz]
>
> I get:
>
>   C:\>gcc -o pow_test pow_test.c
>
>   C:\>pow_test
>   2.20000000000000017764^3.10000000000000008882 = 11.52153412678571875460
>   2.00000000000000000000^3.00000000000000000000 = 8.00000000000000000000
>   4.00000000000000000000^0.50000000000000000000 = 2.00000000000000000000
>   11.52153412678571697825 8.00000000000000000000 2.00000000000000000000
>   0.00000000000000177636 0.00000000000000000000 0.00000000000000000000
>   2.20000000000000000004^3.09999999999999999991 = 11.52153412678571414718
>   2.00000000000000000000^3.00000000000000000000 = 8.00000000000000000000
>   4.00000000000000000000^0.50000000000000000000 = 2.00000000000000000000
>   11.52153412678571414718 8.00000000000000000000 2.00000000000000000000
>   0.00000000000000000000 0.00000000000000000000 0.00000000000000000000
>
> As I understand it, this means I am seeing a small difference between
> compile and run-time evaluation only for powl (2.2, 3.1), i.e., the
> long-double case.
>
>>> The interesting issue is here that gcc uses for constant
>>> math-calculations gmp-library to get result. For more complex
>>> variants, it uses crt's math.  By this it is important for a gcc based
>>> runtime to be in the IEEE 754 floating point compatible to gmp's
>>> results.
>
> Yes, that's a good point.  I hadn't considered the issue of compile-time
> evaluation, and, ideally, it should agree with run-time evaluation, adding
> yet another constraint to our wish list.
>
> Furthermore, if we take the idea of a cross-compiler seriously, we
> would also like the compile-time evaluations of a cross-compiler and
> native compiler to agree.
>
> Hmm...
>
> So, we would like sqrt (x) and pow (x, 0.5) to agree.  We would like
> compile-time and run-time evaluations to agree.  We would like
> cross-compilers and native compilers to agree.
>
> Taken at face value, this would seem to preclude using the platform
> library and platform hardware to compute things like sqrt and pow
> (unless we know they comply with some agreed-upon standard,
> such as IEEE-754).
>
> Just to make the reasoning explicit, if I perform my run-time evaluation
> on platform A, e.g., using platform A's FPU, then a B-->A cross-compiler
> running on platform B won't have platform A's FPU available to perform
> consistent compile-time evaluations.
>
> So you either give up some of the above-wished-for consistency, you
> forgo the potential performance benefits of using the platform's library
> and / or hardware, or you get everybody to use a common standard.
>
>>> ...
>>> I admit that this might be not that obvious to users, why results here
>>> are different, but for a gcc based toolchain we need to play nice with
>>> gcc's internal assumptions.
>>>
>>> Regards,
>>> Kai
>>
>> On closer look for long double, we can have rouding issues for long
>> double precission. I see that powl (0.01L, 0.5L) we have a rounding
>> difference to sqrtl (0.01L). Interestingly not for float or double
>> types.  So we we might need to special-case in pow the case for y ==
>> 0.5 to solve this.
>
> Be, however, cognizant of Greg's concern.  Forcing pow to agree with
> sqrt just for y = 0.5 has the potential to make pow, in a sense, internally
> inconsistent,  It's a tricky trade-off.  (Ideally, both pow and sqrt would
> be "exact" and therefore agree with one another automatically.)
>
> By the way, I also see a problem with run-time powl for powl (2.2L, 0.5L).
> It disagrees slightly from sqrtl (2.2L) and from the compile-time value
> for powl (2.2L, 0.5L).
>
> I modified your test program to compare sqrt with pow (x, 0.5) (and sqrtl
> with powl).  For the record:
>
>   C:\>gcc --version
>   gcc (GCC) 4.5.2 20101002 (prerelease) [svn/rev.164902 - mingw-w64/oz]
>
>   C:\>gcc -o sqrt_pow sqrt_pow.c
>
>   C:\>sqrt_pow
>   1.04880884817015163080, 1.04880884817015163080, 0.00000000000000000000
>   1.48323969741913264109, 1.48323969741913264109, 0.00000000000000000000
>   1.81659021245849494619, 1.81659021245849494619, 0.00000000000000000000
>   1.04880884817015163080 1.04880884817015163080
>   1.48323969741913264109 1.48323969741913264109
>   1.81659021245849494619 1.81659021245849494619
>   0.00000000000000000000 0.00000000000000000000 0.00000000000000000000
>   0.00000000000000000000 0.00000000000000000000 0.00000000000000000000
>   1.04880884817015154699, 1.04880884817015154699, 0.00000000000000000000
>   1.48323969741913258980, 1.48323969741913258970, -0.00000000000000000011
>   1.81659021245849499921, 1.81659021245849499921, 0.00000000000000000000
>   1.04880884817015154699 1.04880884817015154699
>   1.48323969741913258980 1.48323969741913258980
>   1.81659021245849499921 1.81659021245849499921
>   0.00000000000000000000 0.00000000000000000000 0.00000000000000000000
>   0.00000000000000000000 -0.00000000000000000011 0.00000000000000000000
>
> where sqrt_pow.c is:
>
>   #include <stdio.h>
>   #include <math.h>
>
>   double abc[3] = { 1.1, 2.2, 3.3 };
>   long double abcl[3] = { 1.1L, 2.2L, 3.3L };
>   int main()
>   {
>    double r1[3], r2[3], r3[3];
>    long double rl1[3], rl2[3], rl3[3];
>    int i;
>    for (i = 0; i < 3; i++)
>    __mingw_printf ("%.20F, %.20F, %.20F\n",
>    (r1[i] = sqrt (abc[i])),
>    (r2[i] = pow (abc[i], 0.5)),
>    (r3[i] = pow (abc[i], 0.5) - sqrt (abc[i])));
>    __mingw_printf ("%.20F %.20F\n", sqrt (1.1), pow (1.1, 0.5));
>    __mingw_printf ("%.20F %.20F\n", sqrt (2.2), pow (2.2, 0.5));
>    __mingw_printf ("%.20F %.20F\n", sqrt (3.3), pow (3.3, 0.5));
>    r1[0] -= sqrt (1.1);
>    r1[1] -= sqrt (2.2);
>    r1[2] -= sqrt (3.3);
>    __mingw_printf ("%.20F %.20F %.20F\n", r1[0], r1[1], r1[2]);
>    r2[0] -= pow (1.1, 0.5);
>    r2[1] -= pow (2.2, 0.5);
>    r2[2] -= pow (3.3, 0.5);
>    __mingw_printf ("%.20F %.20F %.20F\n", r2[0], r2[1], r2[2]);
>
>    for (i = 0; i < 3; i++)
>    __mingw_printf ("%.20LF, %.20LF, %.20LF\n",
>    (rl1[i] = sqrtl (abcl[i])),
>    (rl2[i] = powl (abcl[i], 0.5L)),
>    (rl3[i] = powl (abcl[i], 0.5L) - sqrtl (abcl[i])));
>    __mingw_printf ("%.20LF %.20LF\n", sqrtl (1.1L), powl (1.1L, 0.5L));
>    __mingw_printf ("%.20LF %.20LF\n", sqrtl (2.2L), powl (2.2L, 0.5L));
>    __mingw_printf ("%.20LF %.20LF\n", sqrtl (3.3L), powl (3.3L, 0.5L));
>    rl1[0] -= sqrtl (1.1L);
>    rl1[1] -= sqrtl (2.2L);
>    rl1[2] -= sqrtl (3.3L);
>    __mingw_printf ("%.20LF %.20LF %.20LF\n", rl1[0], rl1[1], rl1[2]);
>    rl2[0] -= powl (1.1L, 0.5L);
>    rl2[1] -= powl (2.2L, 0.5L);
>    rl2[2] -= powl (3.3L, 0.5L);
>    __mingw_printf ("%.20LF %.20LF %.20LF\n", rl2[0], rl2[1], rl2[2]);
>
>    return 1;
>   }
>
>>
>> Kai
>
> Thanks, Kai, for your insight and explanations.
>
>
> K. Frank

I add a special case to mingw-w64's pow(f|l) function about y=0.5 and
get now for your adjusted testcase the following results, which are
looking much better IMHO.

$ ./tst.exe
1.04880884817015163080, 1.04880884817015163080, 0.00000000000000004142
1.48323969741913264109, 1.48323969741913264109, -0.00000000000000000857
1.81659021245849494619, 1.81659021245849494619, -0.00000000000000000412
1.04880884817015163080 1.04880884817015163080
1.48323969741913264109 1.48323969741913264109
1.81659021245849494619 1.81659021245849494619
0.00000000000000000000 0.00000000000000000000 0.00000000000000000000
0.00000000000000000000 0.00000000000000000000 0.00000000000000000000
1.04880884817015154699, 1.04880884817015154699, 0.00000000000000000000
1.48323969741913258980, 1.48323969741913258980, 0.00000000000000000000
1.81659021245849499921, 1.81659021245849499921, 0.00000000000000000000
1.04880884817015154699 1.04880884817015154699
1.48323969741913258980 1.48323969741913258980
1.81659021245849499921 1.81659021245849499921
0.00000000000000000000 0.00000000000000000000 0.00000000000000000000
0.00000000000000000000 0.00000000000000000000 0.00000000000000000000

I used for this a 32-bit variant and 64-bit variant, which made no
differences in output.
The issue here is that 1.0 branch still uses the inaccurate variant
from Danny, but nowaday's trunk uses for pow, sqrt, log, and exp
functions new implementations, which are treating rounding and
short-path calculation in an ISO-C variant (as gcc expects internally
such a math). Also complex-math part is completely rewritten on
mingw-w64's trunk version.

IMHO it might be a good idea to replace step by step the classical
runtime functions, so that they are not using here msvcrt in first
place. Of course the import-symbols of msvcrt (for compatibility with
VC generated objects need to remain in .def exports, but the actual
functions should come directly from mingwex.

2011/4/30 Kai Tietz <ktietz70@...>:
> 2011/4/30 K. Frank <kfrank29.c@...>:
>> Hello Kai!
>>
>> On Sat, Apr 30, 2011 at 5:51 AM, Kai Tietz <ktietz70@...> wrote:
>>> 2011/4/30 Kai Tietz <ktietz70@...>:
>>>> 2011/4/30 Greg Chicares <gchicares@...>:
>>>>> On 2011-04-29 19:55Z, K. Frank wrote:
>>>>>>
>>>>>> By the way, could someone confirm that mingw does use msvcrt for
>>>>>> sqrt and pow?
>>>>>
>>>>> http://cygwin.com/cgi-bin/cvsweb.cgi/src/winsup/mingw/mingwex/math/?cvsroot=src
>>>>>
>>>>> Looks to me like sqrt(double) calls into msvcrt. But pow() is a
>>>>> different story. Danny implemented it in libmingwex because of
>>>>> some quality-of-implementation problem with msvcrt; then people
>>>>> complained that the libmingwex version was slower, and I can't
>>>>> remember where it wound up, but it's all in the archives.
>>>>>
>>>>>> sqrt (x) and pow (x, 0.5) ought to give the same result (even if not required
>>>>>> to by IEEE-754).
>>>>> ...
>>>>> Not so many years ago (perhaps when msvcrt was written), it was
>>>>> thought that correctly-rounded transcendental functions weren't
>>>>> feasible in practice, so library authors had a looser attitude
>>>>> toward exactness. If you have an ideal implementation of sqrt()
>>>>> but a pow() that's only approximate, should you special-case
>>>>> the power of exactly 0.5? If you do, then pow() and sqrt() will
>>>>> be consistent, but then the powers
>>>>>  x^(0.5*(1-epsilon))
>>>>>  x^(0.5            )
>>>>>  x^(0.5*(1+epsilon))
>>>>> might not be monotone. Some would call that a bad tradeoff.
>>>>
>>>> Well, I spent some efforts into log, exp, and pow implementation. Also
>>>> some other base-math functions I improved in mingw-w64's math library
>>>> in two terms. a) Make those functions behave as ISO-C99 specification
>>>> tells in corner-cases. and b) Improve accuracy in terms of gcc's
>>>> internal used gmp. As we noticed here result differences.
>>>>
>>>> The following example illustrates the issues pretty well:
>>>>
>>>> #include <stdio.h>
>>>> #include <math.h>
>>>>
>>>> double abc[3][2] = { { 2.2, 3.1 }, { 2.0, 3.0 }, { 4.0, 0.5 } };
>>>> long double ab[3][2] = { { 2.2L, 3.1L }, { 2.0L, 3.0L }, { 4.0L, 0.5L } };
>>>> int main()
>>>> {
>>>>  double r[3];
>>>>  long double rl[3];
>>>>  int i;
>>>>  for (i = 0; i < 3; i++)
>>>>  __mingw_printf ("%.20F^%.20F = %.20F\n", abc[i][0], abc[i][1],
>>>>         (r[i] = pow (abc[i][0], abc[i][1])));
>>>>  __mingw_printf ("%.20F %.20F %.20F\n", pow (2.2, 3.1), pow (2.0,
>>>> 3.0), pow (4.0, 0.5));
>>>>  r[0] -= pow (2.2, 3.1);
>>>>  r[1] -= pow (2.0, 3.0);
>>>>  r[2] -= pow (4.0, 0.5);
>>>>  __mingw_printf ("%.20F %.20F %.20F\n", r[0], r[1], r[2]);
>>>>
>>>>  for (i = 0; i < 3; i++)
>>>>  __mingw_printf ("%.20LF^%.20LF = %.20LF\n", ab[i][0], ab[i][1],
>>>>         (rl[i] = powl (ab[i][0], ab[i][1])));
>>>>  __mingw_printf ("%.20LF %.20LF %.20LF\n",
>>>>                         powl (2.2L, 3.1L), powl (2.0L, 3.0L) , powl
>>>> (4.0L, 0.5L));
>>>>  rl[0] -= powl (2.2L, 3.1L);
>>>>  rl[1] -= powl (2.0L, 3.0L);
>>>>  rl[2] -= powl (4.0L, 0.5L);
>>>>  __mingw_printf ("%.20LF %.20LF %.20LF\n", rl[0], rl[1], rl[2]);
>>>>  return 1;
>>>> }
>>
>> Do I understand your point correctly that you are comparing gcc's
>> compile-time and run-time values for pow?
>>
>> For the record, when I run this on a mingw64 build:
>>
>>   gcc (GCC) 4.5.2 20101002 (prerelease) [svn/rev.164902 - mingw-w64/oz]
>>
>> I get:
>>
>>   C:\>gcc -o pow_test pow_test.c
>>
>>   C:\>pow_test
>>   2.20000000000000017764^3.10000000000000008882 = 11.52153412678571875460
>>   2.00000000000000000000^3.00000000000000000000 = 8.00000000000000000000
>>   4.00000000000000000000^0.50000000000000000000 = 2.00000000000000000000
>>   11.52153412678571697825 8.00000000000000000000 2.00000000000000000000
>>   0.00000000000000177636 0.00000000000000000000 0.00000000000000000000
>>   2.20000000000000000004^3.09999999999999999991 = 11.52153412678571414718
>>   2.00000000000000000000^3.00000000000000000000 = 8.00000000000000000000
>>   4.00000000000000000000^0.50000000000000000000 = 2.00000000000000000000
>>   11.52153412678571414718 8.00000000000000000000 2.00000000000000000000
>>   0.00000000000000000000 0.00000000000000000000 0.00000000000000000000
>>
>> As I understand it, this means I am seeing a small difference between
>> compile and run-time evaluation only for powl (2.2, 3.1), i.e., the
>> long-double case.
>>
>>>> The interesting issue is here that gcc uses for constant
>>>> math-calculations gmp-library to get result. For more complex
>>>> variants, it uses crt's math.  By this it is important for a gcc based
>>>> runtime to be in the IEEE 754 floating point compatible to gmp's
>>>> results.
>>
>> Yes, that's a good point.  I hadn't considered the issue of compile-time
>> evaluation, and, ideally, it should agree with run-time evaluation, adding
>> yet another constraint to our wish list.
>>
>> Furthermore, if we take the idea of a cross-compiler seriously, we
>> would also like the compile-time evaluations of a cross-compiler and
>> native compiler to agree.
>>
>> Hmm...
>>
>> So, we would like sqrt (x) and pow (x, 0.5) to agree.  We would like
>> compile-time and run-time evaluations to agree.  We would like
>> cross-compilers and native compilers to agree.
>>
>> Taken at face value, this would seem to preclude using the platform
>> library and platform hardware to compute things like sqrt and pow
>> (unless we know they comply with some agreed-upon standard,
>> such as IEEE-754).
>>
>> Just to make the reasoning explicit, if I perform my run-time evaluation
>> on platform A, e.g., using platform A's FPU, then a B-->A cross-compiler
>> running on platform B won't have platform A's FPU available to perform
>> consistent compile-time evaluations.
>>
>> So you either give up some of the above-wished-for consistency, you
>> forgo the potential performance benefits of using the platform's library
>> and / or hardware, or you get everybody to use a common standard.
>>
>>>> ...
>>>> I admit that this might be not that obvious to users, why results here
>>>> are different, but for a gcc based toolchain we need to play nice with
>>>> gcc's internal assumptions.
>>>>
>>>> Regards,
>>>> Kai
>>>
>>> On closer look for long double, we can have rouding issues for long
>>> double precission. I see that powl (0.01L, 0.5L) we have a rounding
>>> difference to sqrtl (0.01L). Interestingly not for float or double
>>> types.  So we we might need to special-case in pow the case for y ==
>>> 0.5 to solve this.
>>
>> Be, however, cognizant of Greg's concern.  Forcing pow to agree with
>> sqrt just for y = 0.5 has the potential to make pow, in a sense, internally
>> inconsistent,  It's a tricky trade-off.  (Ideally, both pow and sqrt would
>> be "exact" and therefore agree with one another automatically.)
>>
>> By the way, I also see a problem with run-time powl for powl (2.2L, 0.5L).
>> It disagrees slightly from sqrtl (2.2L) and from the compile-time value
>> for powl (2.2L, 0.5L).
>>
>> I modified your test program to compare sqrt with pow (x, 0.5) (and sqrtl
>> with powl).  For the record:
>>
>>   C:\>gcc --version
>>   gcc (GCC) 4.5.2 20101002 (prerelease) [svn/rev.164902 - mingw-w64/oz]
>>
>>   C:\>gcc -o sqrt_pow sqrt_pow.c
>>
>>   C:\>sqrt_pow
>>   1.04880884817015163080, 1.04880884817015163080, 0.00000000000000000000
>>   1.48323969741913264109, 1.48323969741913264109, 0.00000000000000000000
>>   1.81659021245849494619, 1.81659021245849494619, 0.00000000000000000000
>>   1.04880884817015163080 1.04880884817015163080
>>   1.48323969741913264109 1.48323969741913264109
>>   1.81659021245849494619 1.81659021245849494619
>>   0.00000000000000000000 0.00000000000000000000 0.00000000000000000000
>>   0.00000000000000000000 0.00000000000000000000 0.00000000000000000000
>>   1.04880884817015154699, 1.04880884817015154699, 0.00000000000000000000
>>   1.48323969741913258980, 1.48323969741913258970, -0.00000000000000000011
>>   1.81659021245849499921, 1.81659021245849499921, 0.00000000000000000000
>>   1.04880884817015154699 1.04880884817015154699
>>   1.48323969741913258980 1.48323969741913258980
>>   1.81659021245849499921 1.81659021245849499921
>>   0.00000000000000000000 0.00000000000000000000 0.00000000000000000000
>>   0.00000000000000000000 -0.00000000000000000011 0.00000000000000000000
>>
>> where sqrt_pow.c is:
>>
>>   #include <stdio.h>
>>   #include <math.h>
>>
>>   double abc[3] = { 1.1, 2.2, 3.3 };
>>   long double abcl[3] = { 1.1L, 2.2L, 3.3L };
>>   int main()
>>   {
>>    double r1[3], r2[3], r3[3];
>>    long double rl1[3], rl2[3], rl3[3];
>>    int i;
>>    for (i = 0; i < 3; i++)
>>    __mingw_printf ("%.20F, %.20F, %.20F\n",
>>    (r1[i] = sqrt (abc[i])),
>>    (r2[i] = pow (abc[i], 0.5)),
>>    (r3[i] = pow (abc[i], 0.5) - sqrt (abc[i])));
>>    __mingw_printf ("%.20F %.20F\n", sqrt (1.1), pow (1.1, 0.5));
>>    __mingw_printf ("%.20F %.20F\n", sqrt (2.2), pow (2.2, 0.5));
>>    __mingw_printf ("%.20F %.20F\n", sqrt (3.3), pow (3.3, 0.5));
>>    r1[0] -= sqrt (1.1);
>>    r1[1] -= sqrt (2.2);
>>    r1[2] -= sqrt (3.3);
>>    __mingw_printf ("%.20F %.20F %.20F\n", r1[0], r1[1], r1[2]);
>>    r2[0] -= pow (1.1, 0.5);
>>    r2[1] -= pow (2.2, 0.5);
>>    r2[2] -= pow (3.3, 0.5);
>>    __mingw_printf ("%.20F %.20F %.20F\n", r2[0], r2[1], r2[2]);
>>
>>    for (i = 0; i < 3; i++)
>>    __mingw_printf ("%.20LF, %.20LF, %.20LF\n",
>>    (rl1[i] = sqrtl (abcl[i])),
>>    (rl2[i] = powl (abcl[i], 0.5L)),
>>    (rl3[i] = powl (abcl[i], 0.5L) - sqrtl (abcl[i])));
>>    __mingw_printf ("%.20LF %.20LF\n", sqrtl (1.1L), powl (1.1L, 0.5L));
>>    __mingw_printf ("%.20LF %.20LF\n", sqrtl (2.2L), powl (2.2L, 0.5L));
>>    __mingw_printf ("%.20LF %.20LF\n", sqrtl (3.3L), powl (3.3L, 0.5L));
>>    rl1[0] -= sqrtl (1.1L);
>>    rl1[1] -= sqrtl (2.2L);
>>    rl1[2] -= sqrtl (3.3L);
>>    __mingw_printf ("%.20LF %.20LF %.20LF\n", rl1[0], rl1[1], rl1[2]);
>>    rl2[0] -= powl (1.1L, 0.5L);
>>    rl2[1] -= powl (2.2L, 0.5L);
>>    rl2[2] -= powl (3.3L, 0.5L);
>>    __mingw_printf ("%.20LF %.20LF %.20LF\n", rl2[0], rl2[1], rl2[2]);
>>
>>    return 1;
>>   }
>>
>>>
>>> Kai
>>
>> Thanks, Kai, for your insight and explanations.
>>
>>
>> K. Frank
>
> I add a special case to mingw-w64's pow(f|l) function about y=0.5 and
> get now for your adjusted testcase the following results, which are
> looking much better IMHO.
>
> $ ./tst.exe
> 1.04880884817015163080, 1.04880884817015163080, 0.00000000000000004142
> 1.48323969741913264109, 1.48323969741913264109, -0.00000000000000000857
> 1.81659021245849494619, 1.81659021245849494619, -0.00000000000000000412
> 1.04880884817015163080 1.04880884817015163080
> 1.48323969741913264109 1.48323969741913264109
> 1.81659021245849494619 1.81659021245849494619
> 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000
> 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000
> 1.04880884817015154699, 1.04880884817015154699, 0.00000000000000000000
> 1.48323969741913258980, 1.48323969741913258980, 0.00000000000000000000
> 1.81659021245849499921, 1.81659021245849499921, 0.00000000000000000000
> 1.04880884817015154699 1.04880884817015154699
> 1.48323969741913258980 1.48323969741913258980
> 1.81659021245849499921 1.81659021245849499921
> 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000
> 0.00000000000000000000 0.00000000000000000000 0.00000000000000000000
>
> I used for this a 32-bit variant and 64-bit variant, which made no
> differences in output.
> The issue here is that 1.0 branch still uses the inaccurate variant
> from Danny, but nowaday's trunk uses for pow, sqrt, log, and exp
> functions new implementations, which are treating rounding and
> short-path calculation in an ISO-C variant (as gcc expects internally
> such a math). Also complex-math part is completely rewritten on
> mingw-w64's trunk version.
>
> IMHO it might be a good idea to replace step by step the classical
> runtime functions, so that they are not using here msvcrt in first
> place. Of course the import-symbols of msvcrt (for compatibility with
> VC generated objects need to remain in .def exports, but the actual
> functions should come directly from mingwex.
>

PS: The differences shown in results for 32-bit aren't caused by
inaccuarcy of used math directly. Instead it is caused by the fact
that 32-bit use FPU registers (which are 80-bit always). By this the
inline difference without using store, have for double small
differences.  By first storing values and then printing those stored
values you see also for 32-bit a delta of zero.

Interesting this behavior is on 64-bit x86_64 different, as it uses
always precise accuracy to store results in registers. By this you
have really just a 64-bit result for double type, and not indirecty
80-bit (as shown in first case). I altered your testcase a bit in
following lines

...
   for (i = 0; i < 3; i++)
    {
    r1[i] = sqrt (abc[i]);
    r2[i] = pow (abc[i], 0.5);
    r3[i] = pow (abc[i], 0.5) - sqrt (abc[i]);
   __mingw_printf ("%.20F, %.20F, %.20F\n",
   r1[i], r2[i], r3[i]);
    }
...

and get now
$ ./tst.exe
1.04880884817015163080, 1.04880884817015163080, 0.00000000000000000000
1.48323969741913264109, 1.48323969741913264109, 0.00000000000000000000
1.81659021245849494619, 1.81659021245849494619, 0.00000000000000000000
1.04880884817015163080 1.04880884817015163080
1.48323969741913264109 1.48323969741913264109
1.81659021245849494619 1.81659021245849494619
0.00000000000000000000 0.00000000000000000000 0.00000000000000000000
0.00000000000000000000 0.00000000000000000000 0.00000000000000000000
1.04880884817015154699, 1.04880884817015154699, 0.00000000000000000000
1.48323969741913258980, 1.48323969741913258980, 0.00000000000000000000
1.81659021245849499921, 1.81659021245849499921, 0.00000000000000000000
1.04880884817015154699 1.04880884817015154699
1.48323969741913258980 1.48323969741913258980
1.81659021245849499921 1.81659021245849499921
0.00000000000000000000 0.00000000000000000000 0.00000000000000000000
0.00000000000000000000 0.00000000000000000000 0.00000000000000000000

as result. it is an interesting fact that this happens.

Regards,
Kai

On Apr 30, 2011, at 10:15 AM, K. Frank wrote:

> I get:
> 
>   C:\>gcc -o pow_test pow_test.c
> 
>   C:\>pow_test
>   2.20000000000000017764^3.10000000000000008882 = 11.52153412678571875460


Victor Shoup's NTL and my own fp, to 22 digit accuracy, give:

2.20000000000000017764^3.10000000000000008882 = 11.52153412678571783832

so your result is accurate to only 16 digits. Your

2.20000000000000000004^3.09999999999999999991 = 11.52153412678571414718

should be:

2.20000000000000000004^3.09999999999999999991 = 11.52153412678571414732

so your result is accurate to about 20 digits.

Bob

2011/4/30 Bob Delaney <delaneyrm@...>:
>
> On Apr 30, 2011, at 10:15 AM, K. Frank wrote:
>
>> I get:
>>
>>   C:\>gcc -o pow_test pow_test.c
>>
>>   C:\>pow_test
>>   2.20000000000000017764^3.10000000000000008882 = 11.52153412678571875460
>
>
> Victor Shoup's NTL and my own fp, to 22 digit accuracy, give:
>
> 2.20000000000000017764^3.10000000000000008882 = 11.52153412678571783832
>
> so your result is accurate to only 16 digits. Your
>
> 2.20000000000000000004^3.09999999999999999991 = 11.52153412678571414718
>
> should be:
>
> 2.20000000000000000004^3.09999999999999999991 = 11.52153412678571414732
>
> so your result is accurate to about 20 digits.
>
> Bob

Well, what is to be expected as we used here in test application an
%.20F. By specifying here a %.40F its calculated result is
11.5215341267857141471786519559827866032720 and it is has binary
identity for what gmp calculates for 80-bit IEEE 754 floating point.
Which means for our purposes in gcc, it has right accuracy.

Kai

Hello Bob and Kai!

On Sat, Apr 30, 2011 at 5:41 PM, Kai Tietz <ktietz70@...> wrote:
> 2011/4/30 Bob Delaney <delaneyrm@...>:
>>
>> On Apr 30, 2011, at 10:15 AM, K. Frank wrote:
>>
>>> I get:
>>>
>>>   C:\>gcc -o pow_test pow_test.c
>>>
>>>   C:\>pow_test
>>>   2.20000000000000017764^3.10000000000000008882 = 11.52153412678571875460
>>
>>
>> Victor Shoup's NTL and my own fp, to 22 digit accuracy, give:
>>
>> 2.20000000000000017764^3.10000000000000008882 = 11.52153412678571783832
>>
>> so your result is accurate to only 16 digits. Your

This is exactly what I would expect (short of looking at the actual
binary representation).  An IEEE double is 64 bits, with a 52-bit
mantissa, which, with the "implicit" bit, gives 53 bits of precision.
This gives just under 16 decimal digits of precision, which is what
you see.

>>
>> 2.20000000000000000004^3.09999999999999999991 = 11.52153412678571414718
>>
>> should be:
>>
>> 2.20000000000000000004^3.09999999999999999991 = 11.52153412678571414732
>>
>> so your result is accurate to about 20 digits.

As I understand it, gcc defines long double as the x87 80-bit floating point
number.  This has a 64-bit mantissa, and I believe that it does not use an
implicit bit, so this gives slightly more than 19 decimal digits of precision,
again consistent with your result.

>>
>> Bob
>
> Well, what is to be expected as we used here in test application an
> %.20F. By specifying here a %.40F its calculated result is
> 11.5215341267857141471786519559827866032720 and it is has binary
> identity for what gmp calculates for 80-bit IEEE 754 floating point.
> Which means for our purposes in gcc, it has right accuracy.

(Kai's comment seems to confirm that gcc's long double is, indeed,
the 80-bit floating-point number.)

> Kai

Thanks to all for the insight.


K. Frank

2011/4/30 Kai Tietz <ktietz70@...>:
> 2011/4/30 Bob Delaney <delaneyrm@...>:
>>
>> On Apr 30, 2011, at 10:15 AM, K. Frank wrote:
>>
>>> I get:
>>>
>>>   C:\>gcc -o pow_test pow_test.c
>>>
>>>   C:\>pow_test
>>>   2.20000000000000017764^3.10000000000000008882 = 11.52153412678571875460
>>
>>
>> Victor Shoup's NTL and my own fp, to 22 digit accuracy, give:
>>
>> 2.20000000000000017764^3.10000000000000008882 = 11.52153412678571783832
>>
>> so your result is accurate to only 16 digits. Your
>>
>> 2.20000000000000000004^3.09999999999999999991 = 11.52153412678571414718
>>
>> should be:
>>
>> 2.20000000000000000004^3.09999999999999999991 = 11.52153412678571414732
>>
>> so your result is accurate to about 20 digits.
>>
>> Bob
>
> Well, what is to be expected as we used here in test application an
> %.20F. By specifying here a %.40F its calculated result is
> 11.5215341267857141471786519559827866032720 and it is has binary
> identity for what gmp calculates for 80-bit IEEE 754 floating point.
> Which means for our purposes in gcc, it has right accuracy.
>
> Kai
>


To be more precise here IEEE 754 floats for Intel/Amd architecture
x86/x86_64 has following definitions:

float (32-bit):
 digits for mantissa:24
 digit:6, min exp:-125 (-37 10th)
 max exp: 128 (38 10th)
 Epsilon:1.19209e-007
double (64-bit):
 digits for mantissa:53
 digit:15, min exp:-1021 (-307 10th)
 max exp: 1024 (308 10th)
 Epsilon:2.22045e-016
long double (80-bit):
 digits for mantissa:64
 digit:18, min exp:-16381 (-4931 10th)
 max exp: 16384 (4932 10th)
 Epsilon:1.0842e-019

Which means that indeed the accuracy of 20 digits after decimal point
are pretty precise here.

Regards,
Kai

On 2011-04-30 02:04Z, Ocean wrote:
[...using C++...]
> 1.  Does MinGW support RegEx "out of the box", or do I need to install a
> RegEx library seperately?

Yes: a regex facility was added to the C++ standard library in TR1.
To use it, #include <regex>.

> -----Original Message-----
> From: Greg Chicares [mailto:gchicares@...]
> Sent: Saturday, April 30, 2011 1:09 AM
> To: mingw-users@...
> Subject: Re: [Mingw-users] A beginner's questions about RegEx support
> 
> On 2011-04-30 02:04Z, Ocean wrote:
> [...using C++...]
> > 1.  Does MinGW support RegEx "out of the box", or do I need to install a
> > RegEx library seperately?
> 
> Yes: a regex facility was added to the C++ standard library in TR1.
> To use it, #include <regex>.
> 


	Thanks!  From what I can find out, the standard library in TR1 for
RegEx is somewhat limited.  How can I use a different RegEx libary (such as
PCRE)?

Try look at boost

/lars

On Sat, Apr 30, 2011 at 12:26 PM, Ocean <Ocean@...> wrote:
>
>> -----Original Message-----
>> From: Greg Chicares [mailto:gchicares@...]
>> Sent: Saturday, April 30, 2011 1:09 AM
>> To: mingw-users@...
>> Subject: Re: [Mingw-users] A beginner's questions about RegEx support
>>
>> On 2011-04-30 02:04Z, Ocean wrote:
>> [...using C++...]
>> > 1.  Does MinGW support RegEx "out of the box", or do I need to install a
>> > RegEx library seperately?
>>
>> Yes: a regex facility was added to the C++ standard library in TR1.
>> To use it, #include <regex>.
>>
>
>
>        Thanks!  From what I can find out, the standard library in TR1 for
> RegEx is somewhat limited.  How can I use a different RegEx libary (such as
> PCRE)?
>
>
>
> ------------------------------------------------------------------------------
> WhatsUp Gold - Download Free Network Management Software
> The most intuitive, comprehensive, and cost-effective network
> management toolset available today.  Delivers lowest initial
> acquisition cost and overall TCO of any competing solution.
> http://p.sf.net/sfu/whatsupgold-sd
> _______________________________________________
> MinGW-users mailing list
> MinGW-users@...
>
> This list observes the Etiquette found at
> http://www.mingw.org/Mailing_Lists.
> We ask that you be polite and do the same.  Disregard for the list etiquette may cause your account to be moderated.
>
> _______________________________________________
> You may change your MinGW Account Options or unsubscribe at:
> https://lists.sourceforge.net/lists/listinfo/mingw-users
> Also: mailto:mingw-users-request@...
>

On Sat, Apr 30, 2011 at 05:16, Jim Bell <Jim@...> wrote:

> So to purge msvcrt.dll from a non-trivial project, I need to rebuild
> libstdc++-6.dll and libgcc_s_dw2-1.dll with msvcr100 instead of msvcrt,
> right?
>
> If I were to be foolish enough to attempt this, where would I start?
> Download gcc and work through its docs & dependencies? Or are there
> MinGW instructions/hints? (I couldn't find any.)
>
>
I think the most simple way is:
1. cp mingw-w64-crt/lib32/msvcr100.def  mingw-w64-crt/lib32/msvcrt.def
2. cp mingw-w64-crt/lib64/msvcr100.def mingw-w64-crt/lib64/msvcrt.def

Then do a full mingw-w64 build, that's all.

Or you can modify gcc source, change the default library from '-lmsvcrt' to
'-lmsvcr100',
Then do a full mingw-w64 build, but I can not assume it can works for all
case.

--
Dongsheng

Le 29/04/2011 15:12, LRN a écrit :
> On 29.04.2011 16:11, LRN wrote:
>> On 29.04.2011 14:56, TSalm wrote:
>>> Le 29/04/2011 13:28, Greg Chicares a écrit :
>>>> On 2011-04-29 10:49Z, TSalm wrote:
>>>>> Le 29/04/2011 10:44, Gisle Vanem a écrit :
>>>> [...]
>>>>>> I think it's simpler (as the
>>>>>> comment says) to just define "-D_WIN32_WINNT=0x501" on the
>>>>>> gcc command-line.
>>>> [...]
>>>>> In fact, I've added this in a "Configure.h" file , and include it
>>>>> before
>>>>> every #include block :
>>>>>         #define WINVER 0x0501
>>>>>         #define _WIN32_WINNT 0x0501
>>>>> But G++ doesn't find  GetProcessImageFileNameA :
>>>>> D:\dev\C++\WinList\WinList\src\main.cpp:98: error:
>>>>> `GetProcessImageFileNameA' was not declared in this scope
>>>>> I don't understand where I'm wrong...
>>>> Try specifying "-D_WIN32_WINNT=0x501" on the gcc command line instead:
>>>> that's guaranteed to be parsed before any header is included. If that
>>>> works, then "Configure.h" wasn't working the way you expected.
>>>>
>>>> When the preprocessor doesn't behave as hoped, the gcc command-line
>>>> option '-E' is helpful for determining the cause.
>>> No, it's not resolv the problem too :-(
>>> [...]
>>>
>> Whoa-whoa, are you sure you're using MinGW-gcc and not Msys-gcc?
> I mean, mingw-w32api, not msys-w32api
>
Yes, it's in fact an old version of MinGW...(gcc-3.4.5).
I've just installing the latest version, and now, the problem disappear.

Thanks to everybody for your help.

____________
Analyse antivirus effectuée par G Data AntiVirus