There have been special cases in scientific libraries for quite a long time.
The special case for atan(1.0) is used by many programmers to extract the
system stored value for pi/4, knowing that that the math library (or
co-processor) wouldnt compute it instead of simply retrieving the known
result. The simple fact that a function like pow(*,*) will have an accuracy
requirement that allows a bit or two of numerical error while the nature of
the special case sqrt(*) allows full mantissa accuracy to within a possible
toggle of the last bit, is not the same thing. The arc tangent is computed
by libraries by expansion about the argument points 0 and 1, and identities
arctan(|x|) = pi/2 - arctan(1/|x|)
when |x|>1 and such. But this doesn't apply to pow(*,*) because, other than
extraction of the exponents of the two arguments to limit the range of the
logarithm and exponential in the identity
a^b = exp(b*ln(a))
there is no breaking up of things with a natural boundary at b=0.5, although
special cases for b=0 and 1 are often taken, and sometimes checks are made
to determine whether or not b is an integer. The usual process is breaking
up a into exponent and mantissa,
a = (2^ea)*ma
then computation of the logarithm of the result,
ln(a^b) = lr = b*ln(a) = b*(eb*ln(2)+ln(ma)) = ilr + mlr, |mlr|<1
and finally exponentiation of the log result,
a^b = (e^ilr)*(e^mlr).
There is no obvious advantage for checking for any particular value of b
other than 0 and 1 and possibly integers.
Thus it seems to me, the obvious thing to do is to accept the specified
accuracy of pow(*,*) and not do special cases other than 0, 1, and possibly
integers. One might stretch things a bit by evaluating what we would call
pow(a^2, b/2) but that would not likely survive an analysis of speed and
accuracy as compared to pow(a,b).
James K Beard
From: K. Frank [mailto:kfrank29.c@...]
Sent: Tuesday, May 03, 2011 9:32 AM
To: mingw64; MinGW Users List
Subject: Re: [Mingw-users] [Mingw-w64-public] Math library discrepancies
that surprised me.
Hello Kai and Everybody!
A couple of comments about whether to use special cases for
pow (x, y) for things like y = 0.5 or y = integer...
On Tue, May 3, 2011 at 6:58 AM, Kai Tietz wrote:
> 2011/5/3 Peter Rockett <p.rockett@...>:
>> On 03/05/2011 10:11, Kai Tietz wrote:
>> 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.
> I am aware of this what you are telling me, but there is in pow
> function already special-casing for x^y for y with integer values in
> range for y >= INT_MIN and y <= INT_MAX. As here more precise
> calculation via powi is used. This makes that results of pow having
> already discontinuity.
> So pow functions doesn't have the requirement of having continuity.
I'm not aware that either the C standard or the IEEE-754 standard specifies
one behavior or the other. (I think that they don't, but I could be wrong.)
I believe that is has been commonplace over the years to use these special
cases for pow. (This is just the sense I have from personal recollection
oral tradition.) I think speed optimization has been the reason, rather
improved accuracy or agreement with things like sqrt, but perhaps it's been
a mixture of those. Of course, just because there is precedent, doesn't
mean it's the right way to go.
If I had to choose between pow (x, 0.5) agreeing with sqrt (x) and pow (x,
being monotone in y, I think I would choose the latter (i.e., choose not to
use the special cases). But there are good arguments on both sides.
One solution that would be attractive to me would be to implement pow
"exactly" (i.e., bankers' rounding to the nearest representable value), and
then use all the special cases you want as a speed optimization. Now
the special cases won't change the value of pow (assuming sqrt is properly
implemented), so the special cases won't introduce any non-monotonicity
or "discontinuities" into pow.
What I don't know is how hard or expensive this is. I'm pretty sure it's
not too hard to calculate exp "exactly." But I think that calculating log
is harder. Perhaps James Beard knows whether it's practical to calculate
pow "exactly," and maybe even how to do it.
(Another approach, which I very much don't recommend, is to use pow
to calculate sqrt. You get reduced speed and reduced accuracy for
sqrt, but at least you do get sqrt and pow to agree without introducing
non-monotonicity into pow.)
Lastly, I have to believe that this is a well-understood and "solved"
in the numerical-analysis community (of which I am not a part). Perhaps
it would make sense to find out how the experts who have doing this for
fifty-plus years deal with this issue.
And thanks to the mingw and mingw64 teams who have invested the
time and effort to actually write the code we're so breezily discussing.
My comments are in no way meant to be a criticism of those efforts.
I'm just offering up one man's opinion from the user community in the
hope that it might be helpful.
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.
MinGW-users mailing list
This list observes the Etiquette found at
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: