From owner-freebsd-numerics@FreeBSD.ORG Sun Aug 12 23:10:56 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 12152106566B for ; Sun, 12 Aug 2012 23:10:56 +0000 (UTC) (envelope-from peter@rulingia.com) Received: from vps.rulingia.com (host-122-100-2-194.octopus.com.au [122.100.2.194]) by mx1.freebsd.org (Postfix) with ESMTP id 728A98FC08 for ; Sun, 12 Aug 2012 23:10:55 +0000 (UTC) Received: from server.rulingia.com (c220-239-249-137.belrs5.nsw.optusnet.com.au [220.239.249.137]) by vps.rulingia.com (8.14.5/8.14.5) with ESMTP id q7CNAtm4075844 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Mon, 13 Aug 2012 09:10:55 +1000 (EST) (envelope-from peter@rulingia.com) X-Bogosity: Ham, spamicity=0.000000 Received: from server.rulingia.com (localhost.rulingia.com [127.0.0.1]) by server.rulingia.com (8.14.5/8.14.5) with ESMTP id q7CNAm9B021665 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Mon, 13 Aug 2012 09:10:49 +1000 (EST) (envelope-from peter@server.rulingia.com) Received: (from peter@localhost) by server.rulingia.com (8.14.5/8.14.5/Submit) id q7CNAmM7021664 for freebsd-numerics@freebsd.org; Mon, 13 Aug 2012 09:10:48 +1000 (EST) (envelope-from peter) Resent-From: Peter Jeremy Resent-Date: Mon, 13 Aug 2012 09:10:48 +1000 Resent-Message-ID: <20120812231048.GC20453@server.rulingia.com> Resent-To: freebsd-numerics@freebsd.org Received: from vps.rulingia.com (host-122-100-2-194.octopus.com.au [122.100.2.194]) by server.rulingia.com (8.14.5/8.14.5) with ESMTP id q6N6ONBn012304 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Mon, 23 Jul 2012 16:24:23 +1000 (EST) (envelope-from brde@optusnet.com.au) Received: from mail34.syd.optusnet.com.au (mail34.syd.optusnet.com.au [211.29.133.218]) by vps.rulingia.com (8.14.5/8.14.5) with ESMTP id q6N6ONJk011191 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Mon, 23 Jul 2012 16:24:23 +1000 (EST) (envelope-from brde@optusnet.com.au) 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 mail34.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id q6N6O3Ce002338 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Mon, 23 Jul 2012 16:24:04 +1000 From: Bruce Evans Mail-Followup-To: freebsd-numerics@freebsd.org X-X-Sender: bde@besplex.bde.org To: Stephen Montgomery-Smith In-Reply-To: <500CD87D.9060804@missouri.edu> Message-ID: <20120723154108.O1647@besplex.bde.org> References: <20120714120432.GA70706@server.rulingia.com> <20120717084457.U3890@besplex.bde.org> <5004A5C7.1040405@missouri.edu> <5004DEA9.1050001@missouri.edu> <20120717200931.U6624@besplex.bde.org> <5006D13D.2080702@missouri.edu> <20120718205625.GA409@troutmask.apl.washington.edu> <500725F2.7060603@missouri.edu> <20120719025345.GA1376@troutmask.apl.washington.edu> <50077987.1080307@missouri.edu> <20120719032706.GA1558@troutmask.apl.washington.edu> <5007826D.7060806@missouri.edu> <5007AD41.9070000@missouri.edu> <20120719205347.T2601@besplex.bde.org> <50084322.7020401@missouri.edu> <20120720035001.W4053@besplex.bde.org> <50085441.4090305@missouri.edu> <20120720162953.N2162@besplex.bde.org> <20120720184114.B2790@besplex.bde.org> <50095CDE.4050507@missouri.edu> <20120723044308.X6145@besplex.bde.org> <500C5EE5.4090602@missouri.edu> <20120723131233.U1189@besplex.bde.org> <500CD87D.9060804@missouri.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed X-Mailman-Approved-At: Sun, 12 Aug 2012 23:56:00 +0000 Cc: Diane Bruce , Bruce Evans , John Baldwin , David Chisnall , Bruce Evans , Steve Kargl , David Schultz , Peter Jeremy , Warner Losh Subject: Re: Use of C99 extra long double math functions after r236148 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: , Date: Sun, 12 Aug 2012 23:10:56 -0000 X-Original-Date: Mon, 23 Jul 2012 16:24:03 +1000 (EST) X-List-Received-Date: Sun, 12 Aug 2012 23:10:56 -0000 On Sun, 22 Jul 2012, Stephen Montgomery-Smith wrote: > On 07/22/2012 11:12 PM, Bruce Evans wrote: >> On Sun, 22 Jul 2012, Stephen Montgomery-Smith wrote: >> >>> But I will say that your latest version of clog doesn't do as well as >>> mine with this input: >>> >>> x = unur_sample_cont(gen); >>> y = unur_sample_cont(gen); >>> h = hypot(x,y); >>> x = x/h; >>> y = y/h; >>> >>> I was able to get ULPs less than 2. Your program gets ULPs more like >>> up to 4000. >> ... >> What are the actual x and y? I'm not set up to use mpfr. > > The code segment didn't I showed didn't use mpfr. It used unuran. Basically > I am generating random numbers uniformly distributed on the disk |z|=1. You > could also do it using > x = cos(t) > y = sin(t) > where t is a random real number in the interval [0,2 pi]. I see. You are generating |z| near 1 as a side effect of the inaccuracies in sin^2 + cos^2 == 1 :-). I was thinking of generating such numbers using sqrt(). sqrt() is perfectly rounded, so this gives more control. Start with x in [sqrt(2), 1). (I would first try stepping (uniformly) through all of float space). Let y = sqrt(1 - x^2). Then x^2 + y^2 should be 1, but due to inaccuracies it won't be. Since sqrt(2) is irrational and sqrt() is perfectly rounded, in round-to-nearest mode y should always be too large or too small by less than half an ulp. x^2 is imprecise too, so the error for y may be larger. Check a few values on both sides to find the correctly rounded value y. Then this y or the `nextafter' value on one side if it makes x^2 + y^2 (in infinite precision) as close as possible to 1 as possible for this x (by monotonicity). I expect the differences to have a distribution that is somewhat unifom in the mantissa bits, so that the bad cases turn up fast with almost any distribution of x, but worst cases never turn up with blind nonexhaustive testing (testing twice as many cases tends to double the worst case error measured in ulps). There is no obvious reason why the worst cases for the size of |x^2 + y^2 - 1| should only be near x = 1. However, as x moves away from 1 towards sqrt(2), it becomes easier to calculate x^2 + y^2 accurately enough using only doubled precision. > I have been putting a lot of work into casinh, casin, cacosh and cacos, > getting the branches correct. That has exhausted me. I got exhausted with just clog :-). Bruce