Date: Sun, 16 Sep 2012 18:23:06 +1000 (EST) From: Bruce Evans <brde@optusnet.com.au> To: Stephen Montgomery-Smith <stephen@missouri.edu> Cc: freebsd-numerics@freebsd.org Subject: Re: Complex arg-trig functions Message-ID: <20120916174306.H1527@besplex.bde.org> In-Reply-To: <505558A8.6040600@missouri.edu> References: <5017111E.6060003@missouri.edu> <50297CA5.5010900@missouri.edu> <50297E43.7090309@missouri.edu> <20120814201105.T934@besplex.bde.org> <502A780B.2010106@missouri.edu> <20120815223631.N1751@besplex.bde.org> <502C0CF8.8040003@missouri.edu> <20120906221028.O1542@besplex.bde.org> <5048D00B.8010401@missouri.edu> <504D3CCD.2050006@missouri.edu> <504FF726.9060001@missouri.edu> <20120912191556.F1078@besplex.bde.org> <20120912225847.J1771@besplex.bde.org> <50511B40.3070009@missouri.edu> <20120913204808.T1964@besplex.bde.org> <5051F59C.6000603@missouri.edu> <20120914014208.I2862@besplex.bde.org> <50526050.2070303@missouri.edu> <20120914212403.H1983@besplex.bde.org> <50538E28.6050400@missouri.edu> <20120915231032.C2669@besplex.bde.org> <50548E15.3010405@missouri.edu> <5054C027.2040008@missouri.edu> <5054C200.7090307@missouri.edu> <20120916041132.D6344@besplex.bde.org> <50553424.2080902@missouri.edu> <505558A8.6040600@missouri.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
On Sat, 15 Sep 2012, Stephen Montgomery-Smith wrote:
> Hey guys - I have a piece of code like this:
>
> if (ax < DBL_EPSILON && ay < DBL_EPSILON)
> if ((int)ax==0 && (int)ay==0) { /* raise inexact */
> if (sy == 0)
> return (cpack(m_pi_2 - x, copysign(ay, -1)));
> return (cpack(m_pi_2 - x, ay));
> }
>
> Is there a good reason I didn't code it like this?
>
> if (ax < DBL_EPSILON && ay < DBL_EPSILON)
> if ((int)ax==0 && (int)ay==0) /* raise inexact */
> return (cpack(m_pi_2 - x, -y));
>
>
> I'm trying to remember if I coded it the second way, and one of you told me
> to code it the first way. Or maybe I came up with the first way myself -
> maybe I wasn't sure what would happen if y was 0 or -0.
I can only think of [fear of] -y not working right on +-0.
Combined with previous opttimizations and fixes, this gives:
if (ax < DBL_EPSILON && ay < DBL_EPSILON)
return (cpack(m_pi_2 + tiny, -y)); /* PI/2 with inexact...*/
cacos(0 + I*NaN) and several cases for catanh() should similarly add to
m_pi_2 to raise inexact when they return a part with an inexact PI/2.
Otherwise, catrig*.c is remarkably careful about raising inexact.
Refinement: be more careful with the rounding direction (as in fdlibm?):
(1) make sure that m_pi_2 is PI/2 rounded down for the above use (but
round to nearest for other uses). Or maybe, if rounding to nearest
happens to round up, use m_pi_2 - tiny instead of m_pi_2 + tiny so
that the runtime rounding goes in the right direction in hopefully
all rounding modes.
(2) add (or subtract) more than `tiny' to m_pi_2 if necessary to bump it
to the correct side of the infinite-precision PI/2, so that the runtime
rounding goes in the right direction. I'm not sure if this is necessary
or even possible.
Copying the values of PI/2 from the real functions should give both of these,
to the same extent that it gives them for the real functions. The spelling
of the variables should be copied too. The latter is pio2_hi + pio2_lo.
Using pio2_lo instead od `tiny' may be unnecessary and pessimal. pio2_lo
is declared volatile so that it is runtime here, but it is also used in
code where it doesn't need to be volatile. The real functions don't have
a `tiny' variable, and just re-use the general pio2_lo to get inexact here.
So it looks like (2) is unnecessary, with the real functions using pio2_lo
just because it is good enough.
Note that when you need to control the rounding direction or just have
a hi+lo decomposition, it is critical that the constant for the hi
part have a particular value in binary. When it is declared in decimal,
the decimal value should be rounded to match the desired binary value,
so its higher digits will be quite different from the ones of the
infinite- precision full value, even when the hi value is the best
approximation to the full value (and doesn't have bits in it killed
for technical reasons). I noticed this in the opposite direction when
I calculated the decimal and binary values to put in the constant
tables recently. Normally I round in binary and then print the rounded
value in decimal. This looked strange for m_e, so I switched to
printing a value with it rounded in decimal, with the binary rounding
only in a comment. The strangeness is largest when there are many
extra guard digits in the decimal value, like you had originally. It
is unclear whether these digits should match the infinite-precision
value or the expected rounded binary value.
Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20120916174306.H1527>