Date: Tue, 14 Aug 2012 20:46:07 +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: <20120814201105.T934@besplex.bde.org> In-Reply-To: <50297E43.7090309@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> <20120814072946.S5260@besplex.bde.org> <50297CA5.5010900@missouri.edu> <50297E43.7090309@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 05:16 PM, Stephen Montgomery-Smith wrote:
>> On 08/13/2012 04:45 PM, Bruce Evans wrote:
>>
>>> 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.
>>
>> OK.
>>
>>>>> @ --- 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.
>>
>> Yes, I didn't think of that!
>>
>>> ((int)x == 0)' is a safer method of raising inexact for certain x.
>>
>> But this only works if x is less than 1.
>>
>> OK, how about this:
>>
>> sqrt_huge = 1e150;
>> if (sqrt_huge+x>one || sqrt_huge+y>one) ...
>
> Oops
>
> if (sqrt_huge+x>one && sqrt_huge+y>one)
x and y can be DBL_MAX, giving overflow. I think raising overflow is
never correct, since clog() never overflows for large values, and
ccacos() apparently reduces to a rearrangement of clog() for large
values.
BTW, you can probably omit the ISFINITE() tests in:
>>> @ if (ISFINITE(bx) && ISFINITE(by) && (x > RECIP_SQRT_EPSILON_100
>>> || y > RECIP_SQRT_EPSILON_100)) {
since if bx or by is NaN, then it isn't > RECIP_SQRT_EPSILON_100, and
if it is Inf then I think handling it the same as DBL_MAX gives the
correct result. NaNs and Infs now fall through to do_hard_work().
Wouldn't it be easier to never pass them to do_hard_work()?
For just setting inexact, try an expression using `tiny'. There are
many examples to choose from. According to $(grep tiny.*inex *.c):
% e_sinh.c: if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */
% e_sinhf.c: if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */
Ones like you have.
% e_sqrt.c: z = one-tiny; /* trigger inexact flag */
% e_sqrtf.c: z = one-tiny; /* trigger inexact flag */
Works generally, modulo compiler bugs and extra precision, provided z is
used.
% s_erf.c: * erf(x) = sign(x) *(1 - tiny) (raise inexact)
% s_expm1.c: if(x+tiny<0.0) /* raise inexact */
% s_expm1f.c: if(x+tiny<(float)0.0) /* raise inexact */
% s_tanh.c: if(huge+x>one) return x; /* tanh(tiny) = tiny with inexact */
3 more that depend too much on x.
% s_tanh.c: z = one - tiny; /* raise inexact flag */
% s_tanhf.c: if(huge+x>one) return x; /* tanh(tiny) = tiny with inexact */
% s_tanhf.c: z = one - tiny; /* raise inexact flag */
To get z used, try `if ((int)(1 - tiny) == 1)'. To avoid compiler
bugs, it is necessary for `tiny' to be static const volatile (where
`tiny' is already static const). Only a few places in msun use a
volatile `tiny', so you could not worry about the compiler bugs equally
and wait for them to go away or for someone to notice that inexact is
not set properly. clang has similar bugs for huge*huge. gcc doesn't
evaluate huge*huge at compile time, but clang does. Both evaluate
tiny*tiny and 1-tiny at compile time. Spelling 1 as `one' has no
effect on the compiler bugs.
Note that the expressions that mix in x only do so to avoid setting
inexact when x = +-0, or maybe to preserve the sign of x, without using
a branch to classify this x. Here we already have branches to classify
x as large.
Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20120814201105.T934>