2011/4/30 Greg Chicares <[email protected]>:
> 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
------------------------------------------------------------------------------
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
[email protected]
https://lists.sourceforge.net/lists/listinfo/mingw-w64-public
