## Re: [Mingw-users] [Mingw-w64-public] Math library discrepancies that surprised me.

 Re: [Mingw-users] [Mingw-w64-public] Math library discrepancies that surprised me. From: Peter Rockett - 2011-05-03 10:37:14 ```On 03/05/2011 10:11, Kai Tietz wrote: > 2011/5/3 JonY: >> On 5/3/2011 04:05, Kai Tietz wrote: >>> 2011/5/2 Peter Rockett: >>>> On 02/05/2011 17:00, Keith Marshall wrote: >>>>> On 02/05/11 16:30, Charles Wilson wrote: >>>>>>> So, we would like sqrt (x) and pow (x, 0.5) to agree. >>>> The point below about off-topic diversions into precision is valid but >>>> can I reiterate, there is no good reason why x^(0.5) calculated two >>>> different ways should agree since the two different routes will >>>> accumulate different rounding errors. Adding a special case to force >>>> agreement for x^n when n = 0.5 will introduce a discontinuity into the >>>> pow() function and is a very bad idea. The value of returned by pow(x, >>>> 0.5) might be less accurate compared to sqrt(x) but that is part of the >>>> fundamental issue with using floating-point numbers: they are only ever >>>> approximations. You have to take that fact on board. >>> Correct, the ways of calculation of pow (x, 0.5) and sqrt (x) are different. >>> Nevertheless is it an issue in gcc. As if you are writing pow (1.1, >>> 0.5) or sqrt (1.1) - which means gcc optimize it via gmp to a >>> precalculated result - it has identical values. But by using >>> math-library results in a - expectable - difference in result. >>> So IMHO pow should special case here to provide same result as pow >>> does for y = 0.5. >> I would agree with the others that its a horrible hack to treat "0.5" >> special, I don't remember reading anything about requiring pow() results >> to agree with sqrt(). >> >> Can we at least follow a standard spec like: >> ; >> >> We can remove the special case, right? > Well, as some people aren't able to use reply-all here, but have wise > suggestion about implementation details they don't satisfy themself, I > am a bit annoyed. > Anyway, I see that this special-casing is neither explicitly demanded, > but also it isn't explicitly forbidden. So I see here not much reason > for removing this. The reason is not introducing a discontinuity into pow() around 0.5. Forcing the two functions to agree may seem like a good idea but in reality, you are creating another much more insidious problem. You are embedding a potentially disastrous property into pow(). > It might be that I haven't seen the forbidden > thing on the spec. The sqrt special-case has at least one advantage > over the log2/exp2 variant. First it is faster, as just one FPU > command is used to calculate sqrt. Secondly, is provides better result > and compatible one to gcc's internal used soft-floating-point library. > As I don't see in spec that this special casing shouldn't be done, I > think we will keep it on mingw-w64's side. Sorry, but that's a really bad call IMO for the reason stated above. You are ignoring the consequences to try to (mistakenly) enforce mathematical correctness. Even basic FP ops like addition and subtraction deviate from pure mathematical notions of equality. For example, 1 + a, where a is less than the machine epsilon will return precisely unity! Floating point numbers ain't reals and you can't force them to behave exactly like reals... P. > Well, it might be a good idea on mingw.org's side to recheck the exp2, > exp2f, exp2l, and expl implementation, as here you can easily prove > that the calculated results have a rounding issue, which leads for > combined operations - like exp2 (y * log2 (fabs (x)) - to > mis-calculations over two eps. But well, this might be treated as > typical inaccuracy ... > > Regards, > Kai ```