Complex arg-trig functions

[Stephen Montgomery-Smith]

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.