Date: Fri, 21 Sep 2012 17:23:29 +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: <20120921161532.R945@besplex.bde.org> In-Reply-To: <505BD5DA.1070302@missouri.edu> References: <5017111E.6060003@missouri.edu> <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> <20120916134730.Y957@besplex.bde.org> <5055ECA8.2080008@missouri.edu> <20120917022614.R2943@besplex.bde.org> <50562213.9020400@missouri.edu> <20120917060116.G3825@besplex.bde.org> <50563C57.60806@missouri.edu> <20120918012459.V5094@besplex.bde.org> <5057A932.3000603@missouri.edu> <5057F24B.7020605@missouri.edu> <20120918162105.U991@besplex.bde.org> <20120918232850.N2144@besplex.bde! .org> <505BD5DA.1070302@missouri.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
On Thu, 20 Sep 2012, Stephen Montgomery-Smith wrote: > What I did was to make constants called SQRT_6_EPSILON, etc, and then make > your suggested optimizations to float also to double and long double. I'm just using SQRT_EPSILON here in all cases. Partial diffs between your new version of catrig.c and my unmerged version. 50,51c52 < SQRT_3_EPSILON = sqrt(3*DBL_EPSILON), < SQRT_6_EPSILON = sqrt(6*DBL_EPSILON), --- > SQRT_EPSILON = 0x1p-27, /* <= sqrt(DBL_EPSILON) */ Your version depends on sqrt(N*DBL_EPSILON) being a constant enough expression for it to be evaluated at compile time. A fairly large optimization. e_atanh.c uses 2**-28 here. This is significantly smaller than any of the above. It has a large saftely margin. But not after dividing the above by 8. e_atanhf.c uses the same 2**-28 here. This is nonsense. Properly translating the 2**-28 to float precision would have given about 2**-14. Exhaustive testing shows that 2**-13 gives the same results. sqrtf(FLT_EPSILON) is much larger (2**-11.5). That has a negative safety margin -- exhaustive testing shows that 2**-12 loses a little bit of accuracy compared with 2**-13. For catanh*(), we have to bound both x and y, and should have a larger safety margin for both. Non-exhaustive testing shows that 2**-12 works OK in float precision. My previous values had a negative safety marging. In double precision, sqrt(DBL_EPSILON) is not an integer power of 2, and the above gives an additional safety margin by rounding down to an integer power of 2. 304,309c313,315 < ... < if (ax < SQRT_6_EPSILON/8 && ay < SQRT_6_EPSILON/8) < return (z); --- > if (ax < SQRT_EPSILON && ay < SQRT_EPSILON) > if ((int)ax == 0 && (int)ay == 0) /* raise inexact */ > return (z); The divisions by 8 give a larger safety margin than my version. 384,389c390,407 < if (ax < DBL_EPSILON/8 && ay < SQRT_6_EPSILON/8) < return (cpack(m_pi_2, -y)); --- > if (ax < SQRT_EPSILON && ay < SQRT_EPSILON) > if ((int)ax == 0 && (int)ay == 0) > return (cpack(pio2_hi - (x - pio2_lo), -y)); I restored your z term in the approximation so that I could use the same threshold for x and y. This is more accurate and covers more cases. The approximation is now _better_ than the corresponding one in acos*() -- they should be using the extra term too. This has other subtlties involving rounding of Pi/2 -- see later mail. 580,581c596,598 < if (ax < SQRT_3_EPSILON/8 && ay < SQRT_3_EPSILON/8) < return (z); --- > if (ax < SQRT_EPSILON && ay < SQRT_EPSILON) > if ((int)ax == 0 && (int)ay == 0) /* raise inexact */ > return (z); > > I also wrote my own atanhl function so that your inexact optimizations could > be applied to long double as well as double and float. Hmm, I didn't notice that atanhl() was missing. I found that atanh[f] uses an inaccurate approximation for small |x|, so returning atanh*() early for y == 0 and |x| <= 1 breaks not only optimality of the above approximation for small |z|, but also its accuracy. I made a similar real function call to atan() for x == 0 (only implemented in float precision, and the equivalent for cacos and casinh() not tried). Now atanl() is not missing, and atan*(x) is not inaccurate for small x, so calling this early only breaks the optimality of the above. To preserve the optimality, I had to put most of the new special cases later in the function instead of earlier as planned. This makes them less good for avoiding special settings of inexact. Setting inexact early is also bad for optimality, so I no longer try to do it. See the next mail. Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20120921161532.R945>