Date: Sat, 4 Aug 2012 18:49:48 +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: <20120804165555.X1231@besplex.bde.org> In-Reply-To: <501C361D.4010807@missouri.edu> References: <5017111E.6060003@missouri.edu> <501C361D.4010807@missouri.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
On Fri, 3 Aug 2012, Stephen Montgomery-Smith wrote: > On 07/30/12 17:56, Stephen Montgomery-Smith wrote: >> I will be posting my updates on this web page: >> >> http://people.freebsd.org/~stephen/ >> >> The file catrig.c contains implementations of cacos(h), casin(h), >> catan(h). I have been working on it everyday - I seem unable to put it >> down. > > I know you guys are busy and haven't had time to look this over. But in many > ways this is a good thing, as I am constantly finding problems with the code. :-). I made good progress with clog[fl](). The accuracy problems near 0 are long solved and proved to be solved. 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. > 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. > 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))); Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20120804165555.X1231>