From owner-freebsd-numerics@FreeBSD.ORG Sun Aug 12 23:03:58 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 7F6CC106564A for ; Sun, 12 Aug 2012 23:03:58 +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 DEF768FC08 for ; Sun, 12 Aug 2012 23:03:57 +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 q7CN3vMQ075647 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Mon, 13 Aug 2012 09:03:57 +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 q7CN3pmH021211 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Mon, 13 Aug 2012 09:03:51 +1000 (EST) (envelope-from peter@server.rulingia.com) Received: (from peter@localhost) by server.rulingia.com (8.14.5/8.14.5/Submit) id q7CN3puq021210 for freebsd-numerics@freebsd.org; Mon, 13 Aug 2012 09:03:51 +1000 (EST) (envelope-from peter) Resent-From: Peter Jeremy Resent-Date: Mon, 13 Aug 2012 09:03:51 +1000 Resent-Message-ID: <20120812230351.GK20453@server.rulingia.com> Resent-To: freebsd-numerics@freebsd.org From: Peter Jeremy Mail-Followup-To: freebsd-numerics@freebsd.org To: Bruce Evans Message-ID: <20120727210559.GF31169@server.rulingia.com> References: <20120717225328.GA86902@server.rulingia.com> <20120717232740.GA95026@troutmask.apl.washington.edu> <20120718001337.GA87817@server.rulingia.com> <20120718123627.D1575@besplex.bde.org> <20120722121219.GC73662@server.rulingia.com> <500DAD41.5030104@missouri.edu> <20120724113214.G934@besplex.bde.org> <501204AD.30605@missouri.edu> <20120727032611.GB25690@server.rulingia.com> <20120727155521.E4712@besplex.bde.org> MIME-Version: 1.0 Content-Type: multipart/signed; micalg=pgp-sha1; protocol="application/pgp-signature"; boundary="wq9mPyueHGvFACwf" Content-Disposition: inline In-Reply-To: <20120727155521.E4712@besplex.bde.org> X-PGP-Key: http://www.rulingia.com/keys/peter.pgp User-Agent: Mutt/1.5.21 (2010-09-15) Cc: Diane Bruce , John Baldwin , David Chisnall , Stephen Montgomery-Smith , Bruce Evans , Steve Kargl , David Schultz , 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:03:58 -0000 X-Original-Date: Sat, 28 Jul 2012 07:05:59 +1000 X-List-Received-Date: Sun, 12 Aug 2012 23:03:58 -0000 --wq9mPyueHGvFACwf Content-Type: text/plain; charset=us-ascii Content-Disposition: inline Content-Transfer-Encoding: quoted-printable On 2012-Jul-27 17:27:45 +1000, Bruce Evans wrote: >On Fri, 27 Jul 2012, Peter Jeremy wrote: > >> On 2012-Jul-26 22:02:05 -0500, Stephen Montgomery-Smith wrote: >>> I am not getting many responses to the programs I submitted. I >>> understand that you may be all very busy. > >I'm still working on testing and fixing clog. Haven't got near the more >complex functions. I'd suggest that clog() is the most difficult complex function. The others mostly build from other real and complex functions. >For clog, the worst case that I've found so far has x^2+y^2-1 ~=3D 1e-47: =2E.. >so it needs more than tripled double precision for a brute force >evaluation, and the general case is probably worse. I'm working >on a rearrangement so that doubled double precision works in the >general case. Whilst FreeBSD includes ld128 soft-float primitives, clogl() on sparc needs to provide a ld128 result - which implies ld240 in required for intermediate calculations (and, even on x86, ld128 won't provide sufficient additional precision for clogl()). How easy is it to determine which cases need higher precision? Can the majority of the cases be done in double precision or does everything need double double precision? Would the following approach work (all variables are double): x1 =3D (double)(float)x; /* x1 is top 24 bits of fraction */ x2 =3D x - x1; /* x2 is bottom 29 bits of fraction */ x3 =3D 2.0 * x1 * x2; /* 24b * 29b =3D 53b therefore this is exact = */ x1 *=3D x1; /* x1 now has 48 bits */ x2 *=3D x2; /* 58 bits rounded to 53 bits */ y1 =3D (double)(float)y; y2 =3D y - y1; y3 =3D 2.0 * y1 * y2; y1 *=3D y1; y2 *=3D y2; /* add low to high to minimise cancellation */ result =3D x2 + y2; result +=3D x3 + y3; result +=3D x1 + y1 - 1.0; The difficulty is that the final result summation needs to be rearranged to if the magnitudes of x and y are sufficiently different. For the case you presented, you need: result =3D y2; result +=3D x2 + y3; result +=3D x3 + y1; result +=3D x1 - 1.0; However this gives a result 1.094764425253763337e-47 without needing any more than double precision. >> just Appendix G.6 of WG14/N1256 turned into a C array, plus code to >> actually run the tests & interpret the results. So far, it's about >> 1100 lines of which about 1/3 is the test cases and is intended to run > >Yikes. My basic test program is getting too large and complex at 412 >lines. It basically just compares with a known good or different bad >function (with zillions of parameters to select the function and args). >It tests exceptional cases as a side effect. Portably supporting arm, sparc & x86 means copying structures from fpmath.h as well as various _fpmath.h and each family needs different code to display long double values. >Do you do tests for fenv (i.e, that the specified exceptions and no >others are raised)? This might be a problem with external libraries. I have included details of expected exceptions in the data array but haven't written code to actually check for them. >The can deliver values and NaNs for comparison, but nothing requires >them to have the same fenv behaviour as the C library. Currently, I just check that an expected NaN result is a NaN. I don't bother to verify that the NaN returned is tho "correct" NaN. On 2012-Jul-27 08:21:13 -0500, Stephen Montgomery-Smith wrote: >If I remember correctly, the C99 specification is very liberal in its=20 >requirements for cpow, so that >cpow(x,y) =3D cexp(y*clog(x)) >is compliant. The main reason I skipped cpow() is that it's dyadic, which means it needs a different test harness. And the lack of explicit requirements makes very difficult to determine when it should fail. One problem with cexp(y*clog(x)) is that you automatically lose about exponent-bits of precision. Ideally, you need to decompose it in much the same way as pow() - though it's messier. --=20 Peter Jeremy --wq9mPyueHGvFACwf Content-Type: application/pgp-signature -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.19 (FreeBSD) iEYEARECAAYFAlATArcACgkQ/opHv/APuIfMqQCguxTeQRKjIGXCoxC34wy6/RbU CswAoIrrosPW/oTy5+2tUrAKlVkJVM/j =4uf5 -----END PGP SIGNATURE----- --wq9mPyueHGvFACwf--