Complex arg-trig functions

Stephen Montgomery-Smith stephen at missouri.edu
Mon Aug 13 22:16:06 UTC 2012


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) ...




More information about the freebsd-numerics mailing list