Date: Wed, 15 Aug 2012 15:42:00 -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: <502C0998.7040004@missouri.edu> In-Reply-To: <20120814201105.T934@besplex.bde.org> 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> <20120814201105.T934@besplex.bde.org>
next in thread | previous in thread | raw e-mail | index | archive | help
On 08/14/2012 05:46 AM, Bruce Evans wrote:
> 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
>
>
All your solutions depend upon using (1-tiny) with the result being
used. But what if FE_DOWNWARD is set? Then 1-tiny becomes
1-DBL_EPSILON. And then if the result is used, everything is off by 1 ulp.
And
if ((int)(1 - tiny) == 1)
will fail.
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?502C0998.7040004>