From owner-freebsd-numerics@FreeBSD.ORG Wed Aug 15 20:42: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 8F028106566C for ; Wed, 15 Aug 2012 20:42:08 +0000 (UTC) (envelope-from stephen@missouri.edu) Received: from wilberforce.math.missouri.edu (wilberforce.math.missouri.edu [128.206.184.213]) by mx1.freebsd.org (Postfix) with ESMTP id 5DD4A8FC12 for ; Wed, 15 Aug 2012 20:42:08 +0000 (UTC) Received: from [127.0.0.1] (wilberforce.math.missouri.edu [128.206.184.213]) by wilberforce.math.missouri.edu (8.14.5/8.14.5) with ESMTP id q7FKg0Ks061106; Wed, 15 Aug 2012 15:42:01 -0500 (CDT) (envelope-from stephen@missouri.edu) Message-ID: <502C0998.7040004@missouri.edu> Date: Wed, 15 Aug 2012 15:42:00 -0500 From: Stephen Montgomery-Smith User-Agent: Mozilla/5.0 (X11; Linux i686; rv:14.0) Gecko/20120714 Thunderbird/14.0 MIME-Version: 1.0 To: Bruce Evans 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> In-Reply-To: <20120814201105.T934@besplex.bde.org> Content-Type: text/plain; charset=ISO-8859-1; format=flowed Content-Transfer-Encoding: 7bit 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: Wed, 15 Aug 2012 20:42:08 -0000 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.