From owner-freebsd-numerics@FreeBSD.ORG Sun Sep 15 19:48:31 2013 Return-Path: Delivered-To: freebsd-numerics@FreeBSD.org Received: from mx1.freebsd.org (mx1.freebsd.org [8.8.178.115]) (using TLSv1 with cipher ADH-AES256-SHA (256/256 bits)) (No client certificate requested) by hub.freebsd.org (Postfix) with ESMTP id 7BADEA67 for ; Sun, 15 Sep 2013 19:48:31 +0000 (UTC) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.95.76.21]) (using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits)) (No client certificate requested) by mx1.freebsd.org (Postfix) with ESMTPS id 51EA9223B for ; Sun, 15 Sep 2013 19:48:31 +0000 (UTC) Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu [127.0.0.1]) by troutmask.apl.washington.edu (8.14.7/8.14.7) with ESMTP id r8FJmL3D025994; Sun, 15 Sep 2013 12:48:21 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.7/8.14.7/Submit) id r8FJmLH1025993; Sun, 15 Sep 2013 12:48:21 -0700 (PDT) (envelope-from sgk) Date: Sun, 15 Sep 2013 12:48:21 -0700 From: Steve Kargl To: Bruce Evans Subject: Re: (2nd time) tweaks to erff() threshold values Message-ID: <20130915194821.GA25940@troutmask.apl.washington.edu> References: <20130823202257.Q1593@besplex.bde.org> <20130823165744.GA47369@troutmask.apl.washington.edu> <20130824202102.GA54981@troutmask.apl.washington.edu> <20130825023029.GA56772@troutmask.apl.washington.edu> <20130825171910.GA60396@troutmask.apl.washington.edu> <20130827003147.R2328@besplex.bde.org> <20130826174307.GA67511@troutmask.apl.washington.edu> <20130827202751.X1787@besplex.bde.org> <20130827184701.GA76013@troutmask.apl.washington.edu> <20130828205705.S1447@besplex.bde.org> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20130828205705.S1447@besplex.bde.org> User-Agent: Mutt/1.5.21 (2010-09-15) Cc: freebsd-numerics@FreeBSD.org X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.14 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: Sun, 15 Sep 2013 19:48:31 -0000 On Wed, Aug 28, 2013 at 09:31:21PM +1000, Bruce Evans wrote: > On Tue, 27 Aug 2013, Steve Kargl wrote: > > > On Tue, Aug 27, 2013 at 09:24:39PM +1000, Bruce Evans wrote: > >> On Mon, 26 Aug 2013, Steve Kargl wrote: > > >>> ... > >>> /* > >>> - * Coefficients for approximation to erfc in [1.25,1/0.35] > >>> + * Domain [1.25,1/0.35], range [-7.043e-10,7.457e-10]: > >>> + * |log(x*erfc(x)) + x**2 + 0.5625 - r(x)/s(x)| < 2**-30 > >>> */ > >> > >> It is even less clear than for the others that the scaling of the error > >> is correct. The comment could make this clearer by giving the ratio of > >> the correct error scale with the one used in the approximation, at the > >> endpoints (max and min for this, usually at the endpoints). > > > > I'm not sure what you mean here. In looking at error curves, > > it turns out that one of the 2 limits in the range is taken > > from an endpoint. > > I mean that these are absolute errors for a function related to the > result that is especially obscure in this case -- the function is > log(x*result). The error that needs to be minimized is the relative > error of the result, where "relative" is not quite the ratio of the > error and the result -- it is the ratio of the error and [result > rounded to a nearby power of 2; I'm not sure if the rounding is > up or down]. For the absolute error to be close enough to the > relative error, the relative error must not vary much across the > interval, AND the scaling for the abslute error must be compatible. > This is far from clear here. I think taking logs changes the scale > for the absolute error to logarithmic, so this error is not really > absolute. Hopefully the scale for the result is similarly logarithmic, > so that the scales are compatible. There is also a factor of x before > taking logs, so the scaling is not logarithmic either. Before I go down the wrong road, I thought I would check if I'm interpreting the above correctly. To be concrete, in the interval [1/0.35, 11], we consider the function (1) f(x) = log(x * erfc(x)) + x**2 + 0.5625 and I've minimized the maximum error defined by (2) e(x) = f(x) - p(x) / (1 + q(x)) where p(x) = p0 + p1 * z + ..., q(x) = q1 * z + q2 * z**2 + ..., and z = 1 / x**2. In my current code, I get the following error curve: 4e-12 ++-----+-------+------+-------+------+-------+------+-------+-----++ + + + + + **** + + + * 3e-12 ++ ** ** ** ** | * ** *** ** ** *| | * ** ** * ** ** **| 2e-12 ++ * ** * ** * ** *++ | * ** * * ** ** * | 1e-12 ++ * ** * * * ** **++ | * ** ** ** ** * * | | ***** * * * ** ** | 0 ++ * * * * ** ** ** ++ | * * * ** * * * | -1e-12 ++ * * * * ** ** ** ++ | * * ** * * ** * | -2e-12 ++ * * * ** ** ** ** ++ | * *** * ** ** ** | | * ** ** * ** ** | -3e-12 ++ **** ** *** ++ + + + + + + + ****** + + -4e-12 ++-----+-------+------+-------+------+-------+------+-------+-----++ 2 3 4 5 6 7 8 9 10 11 If I rewrite Eq. (1) to isolate erfc(x), I arrive at (3) erfc(x) = F(x) = exp(p(x) / (1+q(x)) - x**2 - 0.5625) / x and if I understand your comment above, I should define a new error curve as (4) E(x) = (erfc(x) - F(x)) / erfc(x) where erfc(x) is the assumed exact result and E(x) is a relative error. The relative error curve, using coefficients determined from (2), then looks like 4e-12 ++-----+-------+------+-------+------+-------+------+-------+-----++ + + + + + **** + + + * 3e-12 ++ ** ** ** ** | * ** *** ** ** *| | * ** ** * ** ** **| 2e-12 ++ * ** * ** * ** *++ | * ** * * ** ** * | 1e-12 ++ * ** * * * ** **++ | * ** ** ** ** * * | | ***** * * * ** ** | 0 ++ * * * * ** ** ** ++ | * * * ** * * * | -1e-12 ++ * * * * ** ** ** ++ | * * ** * * ** * | -2e-12 ++ * * * ** ** ** ** ++ | * *** * ** ** ** | | * ** ** * ** ** | -3e-12 ++ **** ** *** ++ + + + + + + + ****** + + -4e-12 ++-----+-------+------+-------+------+-------+------+-------+-----++ 2 3 4 5 6 7 8 9 10 11 In fact, if I plot e(x) and E(x) on the same (high resolution) figure, the curves are indistingiushable. So, is (4) not the right relative error curve? -- Steve From owner-freebsd-numerics@FreeBSD.ORG Mon Sep 16 11:06:49 2013 Return-Path: Delivered-To: freebsd-numerics@FreeBSD.org Received: from mx1.freebsd.org (mx1.freebsd.org [8.8.178.115]) (using TLSv1 with cipher ADH-AES256-SHA (256/256 bits)) (No client certificate requested) by hub.freebsd.org (Postfix) with ESMTP id 7FB49BFF for ; Mon, 16 Sep 2013 11:06:49 +0000 (UTC) (envelope-from owner-bugmaster@FreeBSD.org) Received: from freefall.freebsd.org (freefall.freebsd.org [IPv6:2001:1900:2254:206c::16:87]) (using TLSv1.2 with cipher ECDHE-RSA-AES256-GCM-SHA384 (256/256 bits)) (No client certificate requested) by mx1.freebsd.org (Postfix) with ESMTPS id 6C705215D for ; Mon, 16 Sep 2013 11:06:49 +0000 (UTC) Received: from freefall.freebsd.org (localhost [127.0.0.1]) by freefall.freebsd.org (8.14.7/8.14.7) with ESMTP id r8GB6ntM089700 for ; Mon, 16 Sep 2013 11:06:49 GMT (envelope-from owner-bugmaster@FreeBSD.org) Received: (from gnats@localhost) by freefall.freebsd.org (8.14.7/8.14.7/Submit) id r8GB6mlt089698 for freebsd-numerics@FreeBSD.org; Mon, 16 Sep 2013 11:06:48 GMT (envelope-from owner-bugmaster@FreeBSD.org) Date: Mon, 16 Sep 2013 11:06:48 GMT Message-Id: <201309161106.r8GB6mlt089698@freefall.freebsd.org> X-Authentication-Warning: freefall.freebsd.org: gnats set sender to owner-bugmaster@FreeBSD.org using -f From: FreeBSD bugmaster To: freebsd-numerics@FreeBSD.org Subject: Current problem reports assigned to freebsd-numerics@FreeBSD.org X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.14 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, 16 Sep 2013 11:06:49 -0000 Note: to view an individual PR, use: http://www.freebsd.org/cgi/query-pr.cgi?pr=(number). The following is a listing of current problems submitted by FreeBSD users. These represent problem reports covering all versions including experimental development code and obsolete releases. S Tracker Resp. Description -------------------------------------------------------------------------------- o stand/175811 numerics libstdc++ needs complex support in order use C99 o bin/170206 numerics [msun] [patch] complex arcsinh, log, etc. o stand/82654 numerics C99 long double math functions are missing 3 problems total.