Date: Sun, 8 Jul 2012 16:58:48 -0700 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: Stephen Montgomery-Smith <stephen@missouri.edu> Cc: freebsd-current@freebsd.org Subject: Re: Use of C99 extra long double math functions after r236148 Message-ID: <20120708235848.GB53462@troutmask.apl.washington.edu> In-Reply-To: <4FF9DA46.2010502@missouri.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>
next in thread | previous in thread | raw e-mail | index | archive | help
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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20120708235848.GB53462>