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>