From owner-freebsd-numerics@FreeBSD.ORG Mon Aug 13 21:46:08 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 E60D7106564A for ; Mon, 13 Aug 2012 21:46:07 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail06.syd.optusnet.com.au (mail06.syd.optusnet.com.au [211.29.132.187]) by mx1.freebsd.org (Postfix) with ESMTP id 77C9A8FC12 for ; Mon, 13 Aug 2012 21:46:06 +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 mail06.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id q7DLjx9G025175 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Tue, 14 Aug 2012 07:45:59 +1000 Date: Tue, 14 Aug 2012 07:45:59 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Stephen Montgomery-Smith In-Reply-To: <50295F5C.6010800@missouri.edu> Message-ID: <20120814072946.S5260@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> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed Cc: freebsd-numerics@freebsd.org, Bruce Evans 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: Mon, 13 Aug 2012 21:46:08 -0000 On Mon, 13 Aug 2012, Stephen Montgomery-Smith wrote: > On 08/13/2012 11:57 AM, Bruce Evans wrote: > >> @ if (sy == 0) >> @ ! return (cpack(cimag(w), -creal(w))); >> @ ! return (cpack(-cimag(w), creal(w))); >> >> The sign of creal(cacos()) is always 1, but this makes it +- the sign >> of atan2(x, y). > > Yes, but the sign of atan2(y,x) will always be the same as the sign of y. So > the two negatives will cancel. 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. > But your code works just as well (and your code doesn't clobber the -0's in > the imaginary part). fabs() forces the sign even for NaNs. >> @ --- 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. ((int)x == 0)' is a safer method of raising inexact for certain x. Expressions using something like huge+x-huge caused me problems in Steve's s_expm1l.c recently. They not only didn't set inexact right, but they also didn't work if the expression is evaluated in extra precision. huge*huge also doesn't work in extra exponent range. i386 by default doesn't have extra precision for doubles, but it does have extra exponent range. So when huge = 1e300, 'return huge*huge;' returns the long double 1e600 in a register. Assignment of this to a double would overflow, possibly much later where it is harder to debug. But sometimes the value is just used from the register and never causes overflow. This is a smaller problem than for huge+x-huge, since the latter is an example where the result of the extra exponent range is used immediately to give another result that doesn't need the extra exponent range to represent. Bruce