Date: Sat, 04 Aug 2012 11:46:15 -0500 From: Stephen Montgomery-Smith <stephen@missouri.edu> To: Bruce Evans <brde@optusnet.com.au> Cc: freebsd-numerics@FreeBSD.org Subject: Re: Complex arg-trig functions Message-ID: <501D51D7.1020101@missouri.edu> In-Reply-To: <20120804165555.X1231@besplex.bde.org> References: <5017111E.6060003@missouri.edu> <501C361D.4010807@missouri.edu> <20120804165555.X1231@besplex.bde.org>
next in thread | previous in thread | raw e-mail | index | archive | help
On 08/04/2012 03:49 AM, Bruce Evans wrote: > On Fri, 3 Aug 2012, Stephen Montgomery-Smith wrote: > I made good progress with clog[fl](). The accuracy problems near 0 > are long solved and proved to be solved. I would be interested to see how you solved this. Would you be willing to post your code somewhere? > Several more special cases > were needed to avoid underflow. It now has about as many special cases > as hypot(), with underflow avoidance in similar places. I'm having > most difficulty with accuracy problems not so near 0 (x*x + y*y near > 1/4 or 2). I expected to get accuracies of 1-1.5 ulps but got 2-3 > ulps. For the inverse trig functions, my accuracy is about 4 ulps, and I am not looking to improve it at all. In particular, the algorithm for casin(h) and cacos(h) is so complicated that I don't see how this can be improved without a lot of work. However if you solve this problem for clog, I can see that it is very likely that your solution will also work for catan(h). If you get this problem licked, I will be interested to see how you do it. >> One of the hardest things is trying to avoid unwarranted underflows. >> For example, when |z| is large (bigger than 1e20) I use the approximation >> >> catanh(z) = 1/z + cpack(0, copysign(M_PI_2, cimag(z))) >> >> If the real part of 1/z underflows, then this is a legitimate >> underflow. But if the imaginary part of 1/z underflows, then this is >> not legitimate, because then this really small number is being added >> to plus or minus I*PI/2, and it should just disappear into rounding >> error. > > The same techniques as used in hypot() and now clog() internals should > work. > Here it is for clog(): > > % /* Avoid underflow when ax is not 1. Also handle zero args. */ > % if (ay < DBL_MAX / 0x1p53 && ay * 0x1p53 <= ax) > % return (cpack(log(ax), atan2(y, x))); > > [Underflow can occur when ax is 1, but in no other cases, because > x*x+y*y cannot be very close to 1 in other cases. The proof that > doubled precision is enough for x*x-1+y*y includes showing by a general > argument that the closeness is at most 2**-157.] > > hypot() does this using a more magic bit test (use difference of > exponents). The magic 53 here is not just DBL_MANT_DIG, but is a good > enough threshold for when log(ax*ax+ay*ay) = log(ax) + > log(1+(ay/ax)*(ay/ax)) rounds to almost the same as log(ax). hypot() > doesn't take logs so has a different magic number (60). 1/z sqares > the hypotenuse but has a factor of x or y, so it probably needs another > different magic number. Again, I would be interested to see the details when you are finished. I am guessing that you have many special cases to handle things when ax or ay are close to DBL_MAX or DBL_MIN. And then in the middle of the night (about when you sent your email), I realized that you can also have inputs like DBL_MIN/0x1p45. These denormals can really mess things up, and trigger false underflows. I see below that you have given thought to this as well. > >> I found this out by testing the underflow flag for lots of random >> inputs. I > > I didn't test the underflow flag, but noticed underflow from inaccuracies. > Gradual underflow (denormals) occurs more often, so it would be worth > tsting for it too. > >> ended up writing my own code to compute the real part of 1/z without >> computing the imaginary part. I have a working solution, but I'm sure >> it isn't optimal. (I see there is code in the IEEE standard PDF >> document which I might copy, or I might copy ideas from the hypot >> function.) > > I think the above will work -- just ignore x or y in x*x+y*y when it > is tiny relative to y or x, and ignore y in y/(x*x+y*y) when it is > tiny relative to x*x+y*y. > > Then there is the problem of x and y tiny or denormal but with their > ration not tiny enough to ignore one of them. 1/z is then large or > infinite, but it it should rarely be infinite, and underflow shouldn't > occur for just x*x+y*y. This is easily handled as in hypot() (not sure > about the method here) and now clog() by scaling: > > % /* Avoid underflow. */ > % if (ax < 0x1p-511) > % return (cpack( > % log(hypot(x * 0x1p511, y * 0x1p511)) - 511 * ln2_lo - > % 511 * ln2_hi, > % atan2(y, x))); Wow - this is really getting down and dirty! Why do you first subtract 511 * ln_2 lo and then 511 * ln2_hi? I am guessing that you are aiming at saving a few ulp by trying to represent M_LN_2 in some kind of 1.5 times double precision? This could quite likely help with my implementation of catan(h). And from a mathematical point of view, catan(h) is somehow more fundamental than casin(h) and cacos(h), so it might be worthwhile pushing for those extra few ulps of accuracy for catan(h). I'll wait until I can see your entire code, then I will peruse it for inspiration.
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?501D51D7.1020101>