From owner-freebsd-current@FreeBSD.ORG Sun Jul 8 23:58:49 2012 Return-Path: Delivered-To: freebsd-current@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [69.147.83.52]) by hub.freebsd.org (Postfix) with ESMTP id 32769106567B for ; Sun, 8 Jul 2012 23:58:49 +0000 (UTC) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.95.76.21]) by mx1.freebsd.org (Postfix) with ESMTP id 0AF0F8FC14 for ; Sun, 8 Jul 2012 23:58:49 +0000 (UTC) Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu [127.0.0.1]) by troutmask.apl.washington.edu (8.14.5/8.14.5) with ESMTP id q68NwmjE053660; Sun, 8 Jul 2012 16:58:48 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.5/8.14.5/Submit) id q68NwmKO053659; Sun, 8 Jul 2012 16:58:48 -0700 (PDT) (envelope-from sgk) Date: Sun, 8 Jul 2012 16:58:48 -0700 From: Steve Kargl To: Stephen Montgomery-Smith Message-ID: <20120708235848.GB53462@troutmask.apl.washington.edu> References: <20120528221731.GA76723@troutmask.apl.washington.edu> <4FC40449.3040602@missouri.edu> <20120528233035.GA77157@troutmask.apl.washington.edu> <4FC40DEA.8030703@missouri.edu> <20120529000756.GA77386@troutmask.apl.washington.edu> <4FC43C8F.5090509@missouri.edu> <20120529045612.GB4445@server.rulingia.com> <20120708124047.GA44061@zim.MIT.EDU> <210816F0-7ED7-4481-ABFF-C94A700A3EA0@bsdimp.com> <4FF9DA46.2010502@missouri.edu> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <4FF9DA46.2010502@missouri.edu> User-Agent: Mutt/1.4.2.3i Cc: freebsd-current@freebsd.org Subject: Re: Use of C99 extra long double math functions after r236148 X-BeenThere: freebsd-current@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Discussions about the use of FreeBSD-current List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Sun, 08 Jul 2012 23:58:49 -0000 On Sun, Jul 08, 2012 at 02:06:46PM -0500, Stephen Montgomery-Smith wrote: > Here is a technical question. I know that people always talk about > ulp's in the context of how good a function implementation is. I think > the ulp is the number of base 2 digits at the end of the mantissa that > we cannot rely on. There are a few definitions of ULP (units in the last place). Google 'Muller ULP' and check goldberg's paper (google 'goldberg floating point'. Here's how I define ULP for my testing and the MPFR code. /* From a das email: ulps = err * (2**(LDBL_MANT_DIG - e)), where e = ilogb(z). I use: mpfr_sub(err, acc_out, approx_out, GMP_RNDN); mpfr_abs(err, err, GMP_RNDN); mpfr_mul_2si(ulpx, err, -mpfr_get_exp(acc_out) + LDBL_MANT_DIG, GMP_RNDN); ulp = mpfr_get_d(ulpx, GMP_RNDN); if (ulp 0.506 && mpfr_cmpabs(err, ldbl_minx) 0) { print warning...; } */ #include "gmp.h" #include "mpfr.h" #include "sgk.h" static mpfr_t mp_ulp_err; /* * Set the precision of the ULP computation. */ void mp_ulp_init(int prec) { mpfr_init2(mp_ulp_err, (mpfr_prec_t)prec); } /* * Given the 'approximate' value of a function and * an 'accurate' value compute the ULP. */ double mp_ulp(mpfr_t approximate, mpfr_t accurate, int dig) { double ulp; mpfr_sub(mp_ulp_err, accurate, approximate, RD); mpfr_abs(mp_ulp_err, mp_ulp_err, RD); mpfr_mul_2si(mp_ulp_err, mp_ulp_err, - mpfr_get_exp(accurate) + dig, RD); ulp = mpfr_get_d(mp_ulp_err, RD); return (ulp); } > So if one were to write a naive implementation of lexp(x) that used > Taylor's series if x is positive, and 1/lexp(-x) is x is negative - well > one could fairly easily estimate an upper bound on the ulp, but it > wouldn't be low (like ulp=1 or 2), but probably rather higher (ulp of > the order of 10 or 20). I've already written an ld80 expl(). I have tested billions if not trillions of values. Don't recall the max ULP I saw, but I know it was less than 1. I don't have an ld128 version, so I won't submit a patch for libm. > So do people really work hard to get that last drop of ulp out of their > calculations? I know very few scientist who work hard to reduce the ULP. Most have little understanding of floating point. > Would a ulp=10 be considered unacceptable? Yes, it is unacceptable for the math library. Remember ULPs propagate and can quickly grow if the initial ulp for a result is large. > Also, looking through the source code for the FreeBSD implementation of > exp, I saw that they used some rather smart rational function. (I don't > know how they came up with it.) See Steven Moshier's book, which is the basis for CEPHES. Instead of using a minimax polynomial approximation for the function in some interval, one uses a minimax approximation with a rational function. -- Steve