From owner-freebsd-numerics@FreeBSD.ORG Tue Aug 14 10:46:13 2012 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 36F33106566B for ; Tue, 14 Aug 2012 10:46:13 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail08.syd.optusnet.com.au (mail08.syd.optusnet.com.au [211.29.132.189]) by mx1.freebsd.org (Postfix) with ESMTP id 4DD758FC0C for ; Tue, 14 Aug 2012 10:46:11 +0000 (UTC) Received: from c122-106-171-246.carlnfd1.nsw.optusnet.com.au (c122-106-171-246.carlnfd1.nsw.optusnet.com.au [122.106.171.246]) by mail08.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id q7EAk7Qw005407 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Tue, 14 Aug 2012 20:46:08 +1000 Date: Tue, 14 Aug 2012 20:46:07 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Stephen Montgomery-Smith In-Reply-To: <50297E43.7090309@missouri.edu> Message-ID: <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> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed Cc: freebsd-numerics@freebsd.org Subject: Re: Complex arg-trig functions X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Tue, 14 Aug 2012 10:46:13 -0000 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