Date: Tue, 14 Aug 2012 07:45:59 +1000 (EST) From: Bruce Evans <brde@optusnet.com.au> To: Stephen Montgomery-Smith <stephen@missouri.edu> Cc: freebsd-numerics@freebsd.org, Bruce Evans <brde@optusnet.com.au> Subject: Re: Complex arg-trig functions Message-ID: <20120814072946.S5260@besplex.bde.org> In-Reply-To: <50295F5C.6010800@missouri.edu> References: <5017111E.6060003@missouri.edu> <501C361D.4010807@missouri.edu> <20120804165555.X1231@besplex.bde.org> <501D51D7.1020101@missouri.edu> <20120805030609.R3101@besplex.bde.org> <501D9C36.2040207@missouri.edu> <20120805175106.X3574@besplex.bde.org> <501EC015.3000808@missouri.edu> <20120805191954.GA50379@troutmask.apl.washington.edu> <20120807205725.GA10572@server.rulingia.com> <20120809025220.N4114@besplex.bde.org> <5027F07E.9060409@missouri.edu> <20120814003614.H3692@besplex.bde.org> <50295F5C.6010800@missouri.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
On Mon, 13 Aug 2012, Stephen Montgomery-Smith wrote:
> On 08/13/2012 11:57 AM, Bruce Evans wrote:
>
>> @ if (sy == 0)
>> @ ! return (cpack(cimag(w), -creal(w)));
>> @ ! return (cpack(-cimag(w), creal(w)));
>>
>> The sign of creal(cacos()) is always 1, but this makes it +- the sign
>> of atan2(x, y).
>
> Yes, but the sign of atan2(y,x) will always be the same as the sign of y. So
> the two negatives will cancel.
y can have any sign I think. But the problem only seemed to happen with
denormals and/or NaNs. There might be a problem with NaNs not giving one
of the canceling negatives.
> But your code works just as well (and your code doesn't clobber the -0's in
> the imaginary part).
fabs() forces the sign even for NaNs.
>> @ --- 408,420 ----
>> @ @ if (ISFINITE(bx) && ISFINITE(by) && (x >
>> RECIP_SQRT_EPSILON_100 || y > RECIP_SQRT_EPSILON_100)) {
>> @ ! /* XXX following can also raise overflow */
>
> I don't see how the code could raise an overflow. The output of clog should
> always be very much less than DBL_MAX. (Originally I had clog(2*z), and that
> could raise an unwarranted overflow.)
@ if (ISFINITE(bx) && ISFINITE(by) && (x > RECIP_SQRT_EPSILON_100 || y > RECIP_SQRT_EPSILON_100)) {
@ ! /* XXX following can also raise overflow */
@ ! if (huge+x+y>one) { /* raise inexact */
@ ! w = clog_for_large_values(z);
@ ! /* Can't add M_LN2 to w since it should clobber -0*I. */
@ ! rx = fabs(cimag(w));
@ ! ry = creal(w) + M_LN2;
@ if (sy == 0)
@ ! ry = -ry;
@ ! return (cpack(rx, ry));
@ }
@ }
clog() won't overflow spuriously, but huge+x+y might.
((int)x == 0)' is a safer method of raising inexact for certain x.
Expressions using something like huge+x-huge caused me problems in
Steve's s_expm1l.c recently. They not only didn't set inexact right,
but they also didn't work if the expression is evaluated in extra
precision.
huge*huge also doesn't work in extra exponent range. i386 by default
doesn't have extra precision for doubles, but it does have extra
exponent range. So when huge = 1e300, 'return huge*huge;' returns
the long double 1e600 in a register. Assignment of this to a double
would overflow, possibly much later where it is harder to debug. But
sometimes the value is just used from the register and never causes
overflow. This is a smaller problem than for huge+x-huge, since the
latter is an example where the result of the extra exponent range is
used immediately to give another result that doesn't need the extra
exponent range to represent.
Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20120814072946.S5260>